In Vivo Characterization of the Toxicological Properties of DPhP, One of the Main Degradation Products of Aryl Phosphate Esters.

BACKGROUND: Aryl phosphate esters (APEs) are widely used and commonly present in the environment. Health hazards associated with these compounds remain largely unknown and the effects of diphenyl phosphate (DPhP), one of their most frequent derivatives, are poorly characterized.
OBJECTIVE: Our aim was to investigate whether DPhP per se may represent a more relevant marker of exposure to APEs than direct assessment of their concentration and determine its potential deleterious biological effects in chronically exposed mice.
METHODS: Conventional animals (FVB mice) were acutely or chronically exposed to relevant doses of DPhP or to triphenyl phosphate (TPhP), one of its main precursors. Both molecules were measured in blood and other tissues by liquid chromatography-mass spectrometry (LC-MS). Effects of chronic DPhP exposure were addressed through liver multi-omics analysis to determine the corresponding metabolic profile. Deep statistical exploration was performed to extract correlated information, guiding further physiological analyses.
RESULTS: Multi-omics analysis confirmed the existence of biological effects of DPhP, even at a very low dose of 0.1mg/mL in drinking water. Chemical structural homology and pathway mapping demonstrated a clear reduction of the fatty acid catabolic processes centered on acylcarnitine and mitochondrial β-oxidation in mice exposed to DPhP in comparison with those treated with vehicle. An interesting finding was that in mice exposed to DPhP, mRNA, expression of genes involved in lipid catabolic processes and regulated by peroxisome proliferator-activated receptor alpha (PPARα) was lower than that in vehicle-treated mice. Immunohistochemistry analysis showed a specific down-regulation of HMGCS2, a kernel target gene of PPARα. Overall, DPhP absorption disrupted body weight-gain processes.
CONCLUSIONS: Our results suggest that in mice, the effects of chronic exposure to DPhP, even at a low dose, are not negligible. Fatty acid metabolism in the liver is essential for controlling fast and feast periods, with adverse consequences on the overall physiology. Therefore, the impact of DPhP on circulating fat, cardiovascular pathologies and metabolic disease incidence deserves, in light of our results, further investigations. https://doi.org/10.1289/EHP6826.

Introduction

Diphenyl phosphate (DPhP) has been used as a main biomarker for assessing exposure to aryl phosphate esters (APEs), especially triphenyl phosphate (TPhP), a molecule suspected of presenting human health hazards. However, DPhP can be produced from degradation of several APEs, including ethylhexyl diphenyl phosphate (EHDPhP) or the resorcinol bis(diphenyl phosphate) (RDP) (Ballesteros-Gómez et al. 2015a, 2015b) and tert-butylphenyl diphenyl phosphate (BPDP) (Heitkamp et al. 1985). Moreover, DPhP itself is largely present in the environment worldwide (Björnsdotter et al. 2018; Fu et al. 2017; Li et al. 2019b; Wang and Kannan 2018; Wang et al. 2019), either owing to its spontaneous/microorganism production from known APEs (Fu et al. 2017; Su et al. 2016), or to its direct use in industry (Mitchell et al. 2019). Most APEs are used as flame retardants. They are added to consumer products and raw materials to delay combustion and meet flammability standards (e.g., Fire initiation and growth, Fire containment, and Non-combustibility test standards), such as the ISO/TC92 Fire Safety code, and the TC89 Fire Hazard Testing code existing in Europe. Moreover, unlike other flame retardants, TPhP and EHDPhP are also largely used as plasticizers and lubricants in hydraulic fluids, rubber, paints, textile coatings, food packaging, and PVC, drastically increasing their presence in the environment. These compounds are not usually covalently linked to plastic materials and can easily leach into the environment (Jamarani et al. 2018).

High vapor pressure of TPhP is also likely to facilitate its release in the air once it is freed from its original material (Meeker et al. 2013). It is not surprising that TPhP and EHDPhP have thus been shown to be ubiquitous components of the human indoor environment where their sources and exposure pathways are quite diverse and heterogeneous with regard to other flame retardants. Indeed, TPhP and EHDPhP quantification in food (Li et al. 2019b; Wang and Kannan 2018), house dust (Björnsdotter et al. 2018; Meeker and Stapleton 2010; Meeker et al. 2013), water (Li et al. 2019a), or air (Wei et al. 2015) has systematically demonstrated their presence in these different matrices, raising awareness on the safety of these compounds. A study characterizing the direct biological effects of DPhP and their relationship to TPhP exposure thus appeared to be of particular relevance to better define the effects and mechanisms of action associated with exposure to APEs in a more comprehensive way.

Historically, DPhP was believed to be produced from TPhP in the liver by oxidase and aryl esterase (Sasaki et al. 1984). However, recent analyses obtained from in vitro cultured hepatocytes revealed that the main metabolites derived from TPhP were hydroxylated and glucuronated forms of TPhP (Su et al. 2014). An important finding is that these metabolites and their equivalents for EHDPhP have recently been detected in human urine samples (Zhao et al. 2019), although DPhP remained the most abundant metabolite in these samples. This finding may indicate either that degradation of APEs does not primarily occur in the human liver or that APEs are rapidly degraded and absorbed as DPhP by the environment. Finally, the presence of DPhP in the environment could be directly due to its importance as a catalyst in polymerization processes.

In line with these hypotheses, a recent study showed that serum hydrolase significantly contributed to TPhP and EHDPhP clearance and production of DPhP (Van den Eede et al. 2016). Bacteria present in the environment and the microbiota may participate in this type of transformation because bacterial metabolism has been shown to fully degrade TPhP, with DPhP initially being the main metabolite released into the biofilm (Hou et al. 2019). Similarly, microsomal preparation of human skin demonstrated the ability of carboxylesterases to efficiently generate DPhP from TPhP (Abdallah et al. 2019), indicating that the likelihood of TPhP reaching subcutaneous fat and blood through this route of exposure is very low. Rapid detection of DPhP in urine samples of women exposed to TPhP through nail polish tends to confirm this hypothesis (Mendelsohn et al. 2016). Finally, DPhP concentrations in the environment were reported to be strongly correlated with TPhP levels present in the same environment (Björnsdotter et al. 2018), hence raising concerns about these potentially hazardous molecules for human health.

The complexity of the routes of exposure described above can cast doubts as to the relevance of in vitro and in vivo studies describing the toxicities associated with APEs such as TPhP. For instance, very high doses of TPhP administered via oral gavage () in adult mice (Chen et al. 2015) or through direct subcutaneous injection () in embryos/pups (Wang et al. 2018) may artificially expose the organism to an irrelevant dose of TPhP and its hydroxylated forms, thus misrepresenting the most common route of DPhP exposure, i.e., through the environment. These types of protocols have mainly led to the conclusion that TPhP has an obesogenic endocrine-disrupting activity. These findings were reinforced in vitro by studies showing that high doses () of TPhP disturb the activity of peroxisome proliferator–activated receptor gamma () (Belcher et al. 2014; Tung et al. 2017), or by the ability of TPhP to enhance the lipogenic activity of the thyroid hormone on isolated chicken embryo hepatocytes (Su et al. 2014). Of note, these doses were highly toxic for mammalian cells, casting doubts on the relevance of these results for human physiology. In addition, cell-based transactivation assays somehow failed to confirm agonistic or antagonistic activities of TPhP on either PPAR or TR nuclear receptors (Kojima et al. 2013). Moreover, a treatment combining four APEs administered at lower individual doses of decreased the body weight (BW) gain of these animals instead of increasing it. Similarly, recent reports indicated that exposure to DPhP and TPhP could disrupt the metabolism in opposite manners (Wang et al. 2018).

Based on these controversial findings, we herein present a large study conducted in mice, in which the routes and doses of exposure to DPhP could mimic those found in humans. Environmentally- and human-relevant doses of DPhP were established from an extensive review of the literature estimating the human mean exposure doses from the environment (Li et al. 2019b; Wei et al. 2015) and the known correlation between environmental measurements and biological fluid concentrations (Björnsdotter et al. 2018; Meeker et al. 2013; Pouech et al. 2015; Wei et al. 2015; Zhao et al. 2019). To validate our choice of using DPhP rather than TPhP or another APE in our toxicity study, we first analyzed DPhP concentrations in blood of mice treated with various doses of both molecules via different routes of exposure. We hypothesized that humans are more likely to be continuously/chronically exposed to TPhP and DPhP owing to the presence of TPhP in the air and dust, rather than temporarily/acutely exposed through nutrition. We have thus analyzed the effects of continuous absorption of these two molecules through drinking water via kinetic measurement of their transformation in mice (bioaccumulation and tissue distribution), in comparison with other acute modes of administration such as oral gavage or tail-vein injection. We then present the data reporting the bioaccumulation and distribution of these molecules in mice. Finally, because our aim was to analyze the biological consequences of a relevant DPhP exposure, we have conducted dual-omics analyses combining metabolomics and transcriptomic measurement on tissue extracts obtained from independent experiments, as well as a histological validation.

Methods

Reagents and Chemicals

Diphenyl phosphate, triphenyl phosphate, and were purchased from Merck with a purity higher than 98%. The confirmation standards of carnitine, acetylcarnitine, palmitoylcarnitine, lauroylcarnitine, decanoylcarnitine, hexanoylcarnitine, stearic acid, oleic acid and linoleic acid, 7-ketocholesterol, , Hypoxanthine, Ethyl oleate, cis-11-Eicosenoic acid, cis-11,14-Eicosadienoic acid, (R)-3-hydroxybutyric acid, Taurocholic acid, Taurodeoxycholate, Taurochenodeoxycholate, Tauroursodeoxycholic acid were obtained from Merck, all of them with a purity higher than 97%. Corn oil for DPhP and TPhP resuspension and animal exposure was also obtained from Merck.

Acetonitrile and heptane of LC-MS quality grade and ammonium formate were purchased from BioSolve. Water and formic acid quality grade optima LC-MS were purchased from Fisher Scientific.

Animal Experiments and Animal Care

Female FVB mice (5 wk old; ) from Charles River Laboratories were used for all experiments. Animals were housed (from 4 wk) in the ANICAN (Center de Recherche en Cancérologie de Lyon) animal facilities accredited by the French Ministry of Agriculture. Food and water were provided ad libitum (lights on: between 8 pm and 8 am; temperature: ; humidity: ). Food was obtained from Safe (150-SP-25) and contained low levels of phytoestrogens and isoflavones (no soybean and alfalfa).

Four procedures were used to treat animals with DPhP and TPhP. Minimal and maximal relevant dose were estimated as indicated in the introduction and from preliminary analyses measuring the quantity of TPhP and DPhP found in mouse biological fluids from a wide range of tested concentrations ( to ). Our aim here was to obtain, from environmentally relevant dose, blood DPhP concentrations close to those found in humans. Three procedures mimicked acute exposure either through vein-tail injection (, 1, 10, or and or ), oral gavage (, 1, 10, or , and or ), or drinking water [TPhP and (C1), (C2), and (C3)]. Experiments were stopped 1 h after treatment, except for that involving drinking water, where mice were euthanized after an overnight exposure. The amount of DPhP ingested through drinking water was higher than the amount ingested through the other routes, though it was spread over a longer period of time. The administered concentrations via drinking water were retained because they provided results that were similar to the acute oral gavage or vein-tail injection. Finally, one procedure mimicked chronic exposure, because mice were continuously treated from 5 to 17 wk through drinking water (, C2, C3). Here, drinking water containing DPhP or the vehicle was changed twice a week until the end of the treatment. The concentrations administered were not adapted to the weight gain nor to the amount of water ingested/d/mouse. During the course of the experiment, drinking water consumption was verified at the group level and fluctuated from to per mouse between the start and the end of the experiment. We thus calculated a mean consumption per day and per mouse of , , and of DPhP with the three different concentrations of DPhP used in drinking water along the duration of the experiments. As reported in the results section, BW gain was slightly different between the treated groups; therefore exposure by weight during the experiment increased by an average of 13% with C3 and was reduced by 4% and 7% in C2 and C1, respectively.

At least four animals were used in each group for acute exposure for each series of experiments described (five for oral gavage). For chronic exposure, 2 independent experiments were carried out with 10 animals per group. Metabolomics and transcriptomic analyses were performed on these separate and independent experiments, reinforcing the strength of the correlations observed. Between groups, animals were randomized according to their weight at the time of reception measured with an electronic balance and acclimated during 1 wk before the start of experiments.

For all exposures, DPhP and TPhP were first solubilized in an emulsion containing water and corn oil (Merck) (1:1 v/v) at . These molecules were then diluted in water to obtain the indicated quantity in for oral gavage and vein-tail injection and the indicated concentration in drinking water. Control animals were always given the same amount of diluted corn oil in water ( in of water for acute exposure and in of drinking water for chronic exposure).

For all conditions, blood was collected by cardiac puncture using a syringe prefilled with heparin (arterial blood collection syringe, BD Preset) after mice were anesthetized with a mix of lethal dose of ketamine () and xylazine (). All mice were finally euthanized by cervical dislocation. Blood was immediately frozen in liquid nitrogen and store at before metabolite extraction. Liver, mammary glands, and visceral fat were also harvested immediately after cervical dislocation. Gall bladder was first removed from the liver. All tissues were briefly rinsed twice in a bath of NaCl solution (0.9%) and once in distilled water to remove blood and bile. Organs were then immediately frozen in liquid nitrogen and store at before metabolite extraction.

Animal experiments were performed in compliance with French and European regulations on the protection of animals used for scientific purposes (EC Directive 2010/63/EU and French Decree 2013–118). They were approved by the Ethics Committee and authorized by the French Ministry of Research (APAFIS#3,680-2016010509529577v5).

Sample Preparation for Targeted Analysis and Metabolomics

Three different extraction methods were used for the four matrices, i.e., a solid–liquid extraction for the liver, a liquid–liquid microextraction for whole blood, and a derived-QuEChERS for visceral fat and mammary glands. For each extraction, an internal standard solution of DPP-d10 and TPP-d15 was used. The same liver extraction was performed for targeted APE measurement and metabolomics analyses.

A mass of of liver was extracted with of acetonitrile and of heptane, assisted by three zirconium balls. After centrifugation, the heptane was discarded, of acetonitrile were transferred to a vial, and a second extraction with of acetonitrile was conducted. Finally, all the extracts (total of of acetonitrile) were pooled, split into two portions of , evaporated at 35°C during approximately 90 min, and stored at . The samples were reconstituted with of water/acetonitrile 90:10 (v/v), prior to UHPLC-HRMS analysis reverse phase chromatography, with of acetonitrile/water 95:5 (v/v) for UHPLC-HRMS analysis HILIC chromatography or with MeOH for targeted UHPLC-MS/MS analysis. For metabolomics, the extraction was conducted from eight samples of each administered concentration to obtain suitable, reliable, and reproducible statistical results.

A volume of of whole blood was extracted with of ACN, assisted by six zirconium beads and microtube homogenization for cell lysis. After centrifugation, of supernatant was recovered. This extraction was repeated twice with recovery volumes of of supernatant. The combined extract was evaporated to dryness under a stream of nitrogen at 40°C and reconstituted with of distilled mono(2-ethyl-5-oxohexyl) phthalate (MeOH) for targeted UHPLC-MS/MS analysis.

A mass of of visceral fat or mammary glands was introduced into a centrifugation tube containing three zirconium balls. For the QuEChERS extraction, of water, of ACN, of heptane, of sodium acetate, and of magnesium sulfate were added into the tube. After a vigorous manual agitation and centrifugation, the heptane phase was removed, and of supernatant ACN phase was recovered, evaporated to dryness under a stream of nitrogen at ambient, and reconstituted with of distilled MeOH for targeted UHPLC-MS/MS analysis.

UHPLC-ESI-MS/MS Analysis

For targeted analysis, separation was carried out using a 1200 HPLC system from Agilent Technologies. The chosen column was a Kinetex XB C18 column ( ID, porosity, ) from Phenomenex. The mobile phase was composed of 0.01% formic acid (FA) and AmAc in distilled pure water (solution A) and AmAc in distilled MeOH (solution B) with a gradient from 5% to 100% (solution B) in 2 min. The columns were maintained at 60°C during the analysis. The chromatographic system was coupled to a triple-stage 5500 QTRAP® from ABSciex with electrospray ionization interface (ESI). The multiple reaction mode monitoring was used for the identification and quantification of DPhP and TPhP (limit of detection, , , , ). The Analyst software (version 1.6.2) was used for instrument control, data acquisition, and data processing.

For metabolomics, separation was carried out using an UltiMate 3000 UHPLC system (Thermo Scientific). The chosen column was a Luna Omega polar (, particle size) (Phenomenex), because it was reported to have a greater affinity for polar compounds, to detect and quantify a larger number of metabolites. A nucleodur HILIC column (, particle size) (Macherey Nagel) was also used to confirm the identified compounds or to detect other compounds that the column did not highlight due to their different polar affinities. The columns were maintained at 30°C during the analysis.

With the column the mobile phase was water/acetonitrile 90/10 (v/v) ammonium formate and 0.01% formic acid (solution A) and acetonitrile ammonium formate and 0.01% formic acid (solution B). The gradient elution had a flow rate of , and it started at 100% of solution A and was maintained for 1 min. The percentage of solution A then decreased until reaching 0% within 10 min. The gradient was held at this percentage for 4 min prior to finally returning to 100% of A and held for 3 min to condition the column for the next injection. The total running time was 18 min. The injection volume was .

For the HILIC column, the mobile phase was water ammonium formate and 0.01% formic acid (solution A) and acetonitrile:water 95/5 (v/v) ammonium formate and 0.01% formic acid (solution B). The elution gradient had a flow rate of and started with 95% of B and was held for 2 min. It was then decreased to 70% in 7 min and 50% in another 2 min. The percentage was held for 4 min to return to the initial percentage in 0.1 min and equilibrated during 10 min. The total running time was 25 min. The injection volume was .

The chromatographic system was coupled to a QTOF mass spectrometer (Maxis Plus, Bruker Daltonics®) with electrospray ionization interface (ESI) operating in positive and negative mode. The following settings were used: capillary voltage of , end plate offset of , nebulizer pressure of 3 bar (N2), drying gas of (N2), and drying temperature of 200°C. A solution of sodium formate and acetate () clusters was used for external calibration at the beginning of each run. The analysis was performed in a full scan over the mass range of 50–1,000 daltons (Da) with a scan rate of . Moreover, the analysis was carried out in profile mode with the following transfer parameters: funnel 1 RF of 200 Vpp, multipole RF of 50 Vpp, quadrupole energy of , collision energy of , stepping basic, and a pre-pulse storage of . The instrument resolution was estimated at 21,244 (full width at half maximum, FWHM) at .

All samples and quality controls (QCs) were injected in the MS mode. QC was prepared by mixing of each sample and was injected every 13 samples, so as to represent a ratio sample/QC of 10%. The data-dependent acquisition (DDA) mode was conducted twice on the QC to obtain MS/MS spectra. It was performed with a cycle time of 3 s and a spectral rate between and in order to record low and high intensity precursors. The samples were analyzed randomly to ensure the representativeness of the results.

The software used for instrument control and to acquire data in positive and negative mode for both chromatographic platforms were OTOF control 4.1, Hystar 4.1 (Bruker Daltonics®), and Data Analysis® 4.4. Preprocessing data (including bucketing, filtering, scaling, and data normalization) and alignment were performed with the algorithm Time aligned Region complete eXtraction 3D (T-Rex 3D) of Metaboscape 4.0 (Bruker Daltonics®). The principal parameters used to create the bucket tables were: ; ; perform MS/MS import, group by collision energy; ; ; electron ion collider ; , , , , in positive mode and , , , , . The bucket tables contained information regarding retention time, ion m/z ratio, neutral mass, detected ions, MS/MS spectrum, and relative peak intensity of features in each sample.

Untargeted Analysis and Annotation Workflow

Multivariate procedure used data of four analyses for all samples and the QCs. Features present simultaneously in QCs and in at least two samples of each group of animals were conserved (7,744 out of 11,168). Data retrieved from each type of metabolomics analysis (HILIC and C18 in both negative and positive mode) were combined and used to perform an unsupervised principal component analysis using the R package ade4 (R Development Core Team). The first four components were tested for sample discrimination and measuring the quality of our analytical procedure.

To determine relevant features for the analysis and sample discrimination, we calculated adjusted -value (false discovery rate, FDR) with the Benjamini-Hochberg method. Fold change and -value between groups were used to build volcano plot and represented which features were significantly different between group (FDR at least ) using the R package EnhanceVolcano (Blighe 2018).

For compound annotations, we selected putative compounds based on a mass deviation inferior to . The discovered formulas were from different databases (analyte DB, ChEBI, ChemSpider, PubChem, HMDB, LIPID MAPS®). This process was done using the compound crawler tool, included in MetaboScape 4.0 and the search tools associated to HMDB and LipidMaps website. Priority was given to the compounds arising from biological databases and previously found and quantified in human blood, urine, and feces (Data obtained from HMDB, ChEBI, and LIPID MAPS®). All features with putative annotations and MS/MS spectra were then tested for confirmation (see next section). The selection of conserved features was initially carried out at this step. These features and all those putatively annotated where no MS/MS spectra were available were then manually curated through physical and chemical properties of annotated metabolites. XLogP3, LogD at , Charge at , Topological surface area and molar mass were taken into consideration. From these data, features with a putative metabolite annotation but a discordant retention time were excluded. When possible, this concordant retention time was verified and cross-validated with both ionization modes ( for a specific chromatography) and eventually one of the feature was excluded from the analysis. Also, we performed the same type of comparison between C18 and HILIC systems (same m/z with concordant RT specific to each chromatographic system) for putative metabolites annotated in both analyses. Moreover, metabolites from a same chemical series (e.g., acylcarnitine, fatty acid) were re-analyzed together to verify the right order of elution according to these properties by using those confirmed by MS/MS as beacons (Figure S1). Finally, when doubts remained between several potential annotations, chemical functions present on the molecules were also used to predict the preferential ionisation mode and exclude unlikely metabolites. Compounds with a carboxylic acid function were for example considered to be preferentially ionized in negative mode, whereas amine, alcohol, and ketone were more likely to be ionized in positive mode.

Metabolite Validation by MS/MS Spectra and Standard

Some of the annotated masses had MS/MS spectra. These experimental spectra were compared to theoretical compounds with the same masses () in silico fragmented using MetFrag (MetaboScape) but also to a bank of experimental spectra available in MetFrag. Tested compounds were retrieved preferentially from HMDB, CHEBI, and LIPID MAPS® databases, or from PubChem when the number of potential candidates was inferior to three with the other databases. All potential compounds tested were scored by comparing their in silico/bank spectrum to our experimental data and by measuring their internal coverage. An example of in silico fragmentation is presented in Figure S2. Annotation was considered validated when its corresponding compounds had the highest score or was in the set of compounds presenting one the highest scores not significantly different [ with score computed from spectral similarity and percentage of internal coverage]. In this case, other parameters were used to finish the validation (presence in biological fluid, RT, preferred ionization mode). Annotations verified by MS/MS validation that were however not sufficient to firmly invalidate other possibilities have been indicated. Finally, to confirm these results, the analytical standards of some compounds with a logical retention time in both columns were purchased. Samples were spiked with these analytical standards and were injected in the same sequence as the unspiked samples in DDA mode with selection of the accurate mass (precursor ion) of the suspected compounds. Then, the extracted-ion chromatogram (EIC) and the MS/MS spectra of the suspected compounds for spiked and unspiked samples were compared. All standard validations (MS/MS and RT comparison) are reported in the Table S1.

ChemRICH, MetaMapp and Hierarchical Clustering Analysis

KEGG ID, PubChem ID, SMILE, and InchiKeys of annotated metabolites were retrieved from public databases (PubChem and KEGG) and used to run ChemRICH and MetaMapp algorithms as indicated in the former publication of these statistical tools (Barretina et al. 2012; Barupal and Fiehn 2017). Adjusted -value and fold enrichment were obtained as indicated previously from the signal intensities associated with the multiple injections of each group of samples. When several features were associated with a same metabolite, we conserved the one with the lowest adjusted -value. ChemRICH analyses were then obtained directly from the ChemRICH interface website, whereas MetaMapp results were loaded in the cytoscape software to represent the network with an organic layout. Presented enrichment plot used Chemical enrichment statistics calculated by the Kolmogorov–Smirnov test. Cluster names were retrieved in an unbiased way from ChemRICH results tables. A couple of clusters could not be directly named by ChemRICH and were in this case curated manually (named through chemical similarity and pathway mapping). All tables associated with these analyses are available in the Table S9-S11.

For hierarchical clustering, Spearman rank correlation and city-block distance were used, respectively, for metabolites and samples sorting. Cluster3 and Java Tree View software were used to generate the clustering and the heat map/tree associated to the result of this algorithm. A table listing the results is available in Table S12.

Sample Preparation for Transcriptomic Analyses and Next-Generation Sequencing (NGS)

RNA extraction for tissue.

Total RNA was extracted and purified using the RNeasy Mini Kit (Qiagen # 74106) from mouse liver. Furthermore, the RNeasy procedure enriches and excludes 5S rRNA, tRNAs, or other low molecular weight RNAs. RNA was isolated on the silica membrane in trusted RNeasy spin columns, with binding capacities of of RNA, according to the supplier’s recommendations (Qiagen).

RNAseq library sequencing and analysis.

For the preparation of the NGS RNA library, RNA concentration was measured using the GE NanoView Spectrophotometer (Biochrom US). The quality of RNA samples was analyzed using the RNA 6000 Pico Kit running on the 2100 BioAnalyzer (Agilent). Total RNA was diluted in a final volume of for a total input of . Only the RNA pools with a RIN score higher than 7 were used in the NGS library preparation prior to sequencing.

mRNAs arising from the pooled livers of 10 animals tested in each group were isolated using the NEBNext Poly(A) mRNA Magnetic Isolation Module from of total RNA. Three analytical preparations were performed at this step. The isolation procedure was based on the selection of mRNA using oligo dT beads directed against polyA tails of intact mRNA. Second, the NGS libraries were created from mRNA isolated using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England BioLabs). Final libraries were sequenced on an Illumina NextSeq 500 on a high output flow cell with paired-end read lengths.

The sequencing reads were obtained after demultiplexing the raw sequencing data using bcl2fastq v2.19.1.403 (version v2.15.0 for NextSeq™ 500 and HiSeq® X Systems; Illumina), after having validated the quality controls of each sample using the FastQC version 0.11.5 software (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The alignment files were generated with STAR version 2.5.2b (University of Birmingham) in the 2-pass mode. We used the GRCm38, version M16 (Ensembl 91) as reference. This mode is known to improve the detection of more reads mapping novel splice junctions.

Gene Ontology and Statistical Analyses of NGS

Genes were filtered through their obtained RPKM. Only genes sufficiently expressed in control conditions were conserved for further analyses (). Principal component analysis (PCA) and Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) were performed on the R platform, using ROPLS package and 3D plots for visualization of the results. Gene ontology analyses were performed by using the Eigen values associated with the most significant axis used to construct PCA and OPLS-DA. With these methods, all conditions, namely the three different DPhP concentrations and the vehicle, were thus taken into consideration and simultaneously compared. Threshold for retained Eigen values used to select genes in these analyses are indicated directly in the figures. Functional enrichment analysis (STRING© Consortium 2019) was then performed by sorting the selected genes according to their Eigen values. A volcano plot (R package EnhanceVolcano Plot) was used to represent the discriminating functions for the studied axis through the enrichment score and the false discovery rate generated by the STRINGv11 algorithm. Only significant functions were retrieved.

Alternatively, paired analyses were performed using each concentration of DPhP against the vehicle. In this case, retained genes were selected through the fold change between both conditions. Gene ontology analysis was then performed identically with Functional Enrichment Analysis (STRING© Consortium 2019).

Finally, in a complementary approach, density of the protein networks associated with discriminating functions and encoded by the genes selected with these different approaches were measured through protein-protein association networks (STRING© Consortium 2019). High number of nodes and edges indicate that large functional complexes are disturbed through multiple genetic regulations induced by DPhP exposure, reinforcing the probability of that function being disturbed by this condition. We used high-confidence settings to retain the experimentally validated interactions. A protein–protein interaction enrichment -value was calculated against an identical number of random proteins [PPI enrichment -value (Szklarczyk et al. 2019)].

Bipartite View of Enrichment Network

NetworkAnalysis server (Mc Gill University) was used through the List Enrichment function. Functions belonging to the KEGG database were used and only those with a were conserved (Zhou et al. 2019). Metabolic pathways and carbon metabolism functions were removed due to their very high coverage of the genome. The bipartite view was then selected, and an auto-layout was applied. Clustering between functions was then enhanced through the Force Atlas tool. Clusters were then manually highlighted and default encoding of the functions through their -value was converted into a continuous green-red scale.

Immunohistochemistry Analyses

For histological examination, liver samples were fixed in 10% buffered formalin and embedded in paraffin. tissue sections of formalin-fixed, paraffin-embedded tissue were prepared according to conventional procedures. Sections were then stained with hematoxylin and eosin.

Immunohistochemistry (IHC) was performed on an automated immunostainer (Ventana Discovery XT, Roche) using Omnimap DAB Kit according to the manufacturer’s instructions. Sections were incubated for 60 min at 37°C with the following antibodies: Anti-Perilipin2 (AtlasAntibodies-HPA016607), Anti-Hmgcs2 (Santa Cruz-sc-376,092), Anti-PPARalpha (Abnova-MAB12349), Anti-Pck1 (AtlasAntibodies-HPA006507), Anti-FBP1 (Abcam-ab109020) (all diluted at 1:100). An antirabbit/mouse – HRP was applied on sections during 16 min at 37°C. Staining was visualized with DAB solution with 3,3-diaminobenzidine as a chromogenic substrate. Sections were counterstained with Gill’s hematoxylin and finally scanned with panoramic scan II (3D Histech) at 20×. Scoring was performed by three independent investigators using a staining scale ranging from 0 to 10. The means were then calculated and encoded as follows. , , , .

Statistical Test and Accession Numbers

Statistical tests were performed as indicated in the caption of the figure. Briefly, for simple comparison between two experimental groups, Student’s -test was used. For more complex statistical procedures, test already integrated in the previously described algorithms were conserved without any modification.

Metabolomics data have been deposited to the EMBL-EBI MetaboLights database (doi: 10.1093/nar/gkz1019; PMID: 31691833) with the identifier MTBLS2038. The complete data set can be accessed at https://www.ebi.ac.uk/metabolights/MTBLS2038.

Microarray analyses have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) and are accessible through GEO accession number GSE158265.

Results

Correlation Analysis between TPhP/DPhP Exposure and Their Level in Blood and Liver

When 0.1 or DPhP was injected directly into vein-tail or administered by oral gavage, DPhP was readily quantified after 1 h in a dose-dependent manner in whole blood (Figure 1A and 1B). Conversely, in animals exposed to these same concentrations of TPhP, this molecule could not be detected in their blood (). After administration of or TPhP, TPhP was only quantified above the LOD in the blood of two animals at and , exposed to following intravenous injection and oral gavage, respectively. In all other animals (18 out of 20 animals), TPhP remained undetected. To determine whether TPhP transformation into DPhP was the reason for the lack of detection of TPhP in the bloodstream, we also quantified DPhP in these same experiments. DPhP could not be detected in the blood of animals treated with 0.1 or of TPhP (, ). At the highest doses of TPhP (10 and ), we were able to quantify DPhP in the blood of animals, but in a level comparable to that obtained after exposure to DPhP a hundred times lower (Figure 1A, 1B). Therefore, the large majority of injected or force-fed TPhP seemed not to have been metabolized into DPhP in the animals.

Figure 1.

Acute exposure to DPhP and TPhP. (A–D) Four or 5 animals were exposed to the indicated doses of DPhP or TPhP via the indicated routes. After 1 h, DPhP was quantified from the whole blood of each animal. Concentrations obtained have been plotted on the box-and-whisker plot, indicating significant -values between tested conditions and the vehicle (Student’s -test). (E–F) Identical experiment with DPhP quantification performed on liver extracts from each animal. (, and 3rd quartile, , and smallest value within 1.5 × the interquartile range. ummary data in Table S25.) Note: LOD, limit of detection; LOQ, limit of quantification; ns, nonsignificant.

We further compared these results to a more continuous exposure to both molecules present in the drinking water of animals to mimic their chronic ingestion from swallowed dust. We used 3 concentrations of (C1), (C2) and (C3), equivalent to , , and of each molecule ingested overnight (active period for mice). These quantities were comparable to those used with other routes of exposure previously given through oral gavage or vein-tail injection. At the two highest concentrations, DPhP was still dose-dependently quantified in whole blood with lower concentrations measured than those obtained previously (Figure 1C). This finding was expected because the dose was now spread over a much longer period. Conversely, we did not record any measurable amount of TPhP in animals exposed to this same molecule through drinking water (). Most important, the amount of DPhP retrieved from the whole blood of animals exposed by this route to TPhP was not significantly different from those of animals treated with the vehicle (Figure 1D).

Because the liver is the first organ to process exogenous molecules absorbed from the digestive tract, we further determined whether these molecules were present in this organ using the same exposure doses via drinking water. DPhP was detected and quantified in a dose-dependent manner at all concentrations tested (Figure 1E). Conversely, neither TPhP nor DPhP were measured at a significant different amount in the liver of mice exposed to TPhP or the vehicle through drinking water [ for all animals and Figure 1F].

Bioaccumulation and Distribution of DPhP in Chronically Exposed Mice

Based on these analyses, we next focused on the consequences of chronic exposure to DPhP. Mice were exposed daily to similar doses (C1, C2, or C3) of DPhP in their drinking water over a 12-wk period. We analyzed four tissues, namely whole blood, the liver, visceral fat, and the mammary glands, potentially presenting a tropism for the molecule, either due to their chemical characteristics (presence of hydrophobic constituent for mammary gland and visceral fat) or to the route of exposure used for these experiments (whole blood and liver). In the blood, DPhP concentrations were significantly higher in all animals chronically exposed to DPhP in drinking water at and in comparison with control animals (Figure 2A). DPhP was also measured in the liver, visceral fat, and mammary glands at all tested concentrations in a dose-dependent manner (Figure 2B–D). These results confirmed the ability of chronically ingested DPhP to reach various tissues in mice.

Figure 2.

Chronic exposure to DPhP. (A–D) Ten animals were chronically exposed to the indicated concentrations of DPhP for 12 wk through drinking water. DPhP was then quantified in the indicated biological tissue. Blood of each animal was analyzed separately, whereas analytical replicates were performed from pooled tissue extracts. Concentrations obtained have been plotted on the box-and-whisker plot indicating significant -values between tested conditions and the vehicle (Student’s -test). (, and 3rd quartile, , and smallest value within 1.5 × the interquartile range. Summary data in Table S26). Note: ns, nonsignificant.

Biological Effects of Exposure to DPhP through Metabolomics Analyses of the Liver

Because DPhP was abundant in the liver of all treated animals, we focused our subsequent experiments on this organ to measure the possible biological consequences of the presence of DPhP in this tissue. The metabolite profile of each group [untreated control (CTRL) or treated with three doses (C1, C2, or C3) of DPhP] was obtained with two chromatographic methods (C18 and HILIC) using two ionization modes for the mass spectrometry analyses of eluted metabolites.

Unsupervised PCA was first performed with features present simultaneously in QCs and samples (7,744 out of 11,168 by combining all acquisitions). The consistency of our analytical process was verified by confirming the spatial clustering of QCs on the first four principal components (PCs) (Figure 3A) and their intermediate position in comparison with real samples as expected from their preparation (see “Methods”). Acquired metabolic profiles discriminated the four groups of samples even for the lowest concentration C1. An interesting finding was that samples corresponding to C2 and C3 were clearly not separated by the same PCs. PC1–2 strongly discriminated C2 from CTRL, whereas C3 was only separated from the CTRL using PC3–4. Therefore, two sets of metabolites were differently altered by DPhP according to concentration at which it was administered.

Figure 3.

Hepatic metabolomics analysis of animals exposed to DPhP. (A) Principal Component Analysis (PCA) was performed to discriminate groups of samples based on the amount of retained ion m/z (see “Methods”). Explained variance with the four first axes has been indicated. , , , , . (B) An UpSet plot was calculated reporting the number of metabolites annotated in each analytical procedure and the number of each type of co-occurrence between these procedures. (C) An Euler diagram was created reporting the number of annotated metabolites and their level of validation (expected RT and m/z, MS/MS, spiked chemical standard). (D–G) Volcano plots were created through observed median fold change and adjusted -value algorithm), comparing the metabolite level for the C2 group used against control values for each type of indicated analytical platform. Most significant metabolites were annotated when possible and indicated directly on the plot. Cutoff for for (D), 3/0.002 for (E), 2/0.05 for (F), and 2/0.05 for (G).

All features were further annotated at this stage through three procedures of validation (see methods). A total of 633 metabolites were conserved from this process (Tables S2–S8) with their level of validation and their co-occurrence in the different types of acquisitions indicated in Figure 3B,C and the Tables S2–S5. Volcano plot representations were then used to compare each treated group to control animals for each type of chromatography and ionization to capture potential features and metabolites present in several analyses and thus highly relevant (Figure 3D–G and S3A–SD). From this first analysis, we noted that a subset of acylcarnitines, key metabolites of fatty acid degradation, were significantly depleted in exposed animals at all concentrations when compared to those in vehicle-exposed animals. An important observation was that this effect could be observed with both types of chromatography confirming the robustness of the identified modification. From these first hits, we decided to perform an in-depth analysis combining our 633 annotated metabolites. Due to the heterogeneity of these metabolites and, for some of them, their nonappurtenance to the endogenous metabolism, we performed an initial comparison based on their chemical ontology that can be represented through a circular tree plot of annotated compounds (Figure 4A and S4A and SB). This type of analysis was shown to eliminate bias due to incomplete mapping and network size heterogeneity (Barupal and Fiehn 2017). Moreover, it is weakly sensitive to isomer confusion belonging to a similar chemical series, providing some flexibility for the analysis of highly similar metabolites (e.g., phospholipids) with incomplete validation. From this statistical analysis, an enrichment plot could be built, displaying the significance of the observed change for structural homologies, the size of the cluster, and the homogeneity of their variation (Figure 4B, Figure S4C and S4D, and Table S9–S11). In line with the multivariate analysis, the intermediate concentration C2 induced the largest number of altered chemical clusters. Among the different disturbed clusters, acylcarnitines were the most affected in liver extracts at all DPhP doses. Lower levels of ethyl and choline-esterified fatty acids were also highlighted using this analysis. An important finding was that clusters associated with ketone bodies usually produced from fatty acid degradation also showed a uniform lower level of these ketone bodies after exposure to DPhP. Finally, multiple metabolic pools without clear functional connection with acylcarnitines were also altered.

Figure 4.

Metabolic network analysis through structural identity and pathway mapping. (A) Data associated with 635 annotated metabolites were used to perform a Chemical Similarity Enrichment Analysis (ChemRICH). A Tanimoto chemical similarity mapping form a clustered circular similarity tree. Dark black lines indicate boundaries of clusters that were significantly different in exposed animals at the concentration C2 vs. control mice (). Higher metabolite levels in DPhP-exposed mice (compared with those exposed to vehicle) have been labeled as red nodes; lower levels have been marked in blue. Acylcarnitine Cluster label has also been indicated. (B) ChemRICH set enrichment statistics plot for the same metabolites as (A), extracted from exposed animals with the C2 concentration of DPhP vs. control. Each node reflects a significantly altered chemical cluster of metabolites. Enrichment -values have been obtained using the Kolmogorov-Smirnov test, and median XlogP3 is calculated from individual value associated to each metabolite present in the cluster. Node sizes represent the total number of metabolites in each set of clusters. Continuous scale color for the nodes shows the proportion of fully increased (red) or fully decreased (blue) family of compounds in exposed mice compared to control mice. Nodes with an intermediate color (pink to purple) have both increased and decreased metabolites. (C) MetaMapp visualization of the same metabolomics data highlighting the differential metabolic regulation and the organization of metabolic clusters based on KEGG reactant pair information and Tanimoto chemical similarity matrix. Only metabolites significantly different (adjusted ) are used to build the network. Metabolite levels that were higher in exposed mice have been labeled as red nodes; those with lower levels have been marked in blue (continuous scale with white intermediate). -Values are encoded in node size (lowest ). Cluster label has been indicated. Metabolites previously clustered together based on their structural similarity could now be separated according to their different pathway mapping. (D–F) Identical data obtained for the same 633 metabolites and the 4 groups of exposed and control animals were clustered hierarchically through their relative level (centroid linkage with Spearman correlation (metabolites) and city-block distance (samples). Features corresponding to a same metabolite present in different analytical procedure have been left separated, explaining why a metabolite can appear several times in the heat map. Yellow-blue encoding has been used to represent these metabolites according to their absolute amounts. Distance between levels of samples and metabolites have been shown as two tree plots. The complete heat map is presented in the Figure S5. Subclusters containing metabolites globally less present in exposed animal such as the acylcarnitines, (R)-3-hydroxybutyric acid have been highlighted and enlarged in the indicated left-hand panels. Similarly, a subcluster containing metabolites more present in the C1–3 groups such as N-acylethanolamide is presented in the right-hand panel.

Stress markers, such as N-acylethanolamides were homogeneously more concentrated in extracts from exposed animals. Conversely, sulphated phenol group, 3-oxosteroid and bile acids were more specifically disturbed in the C2 group of samples, whereas groups of metabolites such as those containing a naphthalene or an indole ring showed perturbations with opposite directions between C2 and C3. Next, to determine how these different clusters were organized around endogenous metabolic networks, we combined chemical and biochemical mapping in a joint analysis using the MetaMapp algorithm (Barupal et al. 2012) (Figure 4C and S4E and S4F). Strong connections were observed among the reduced concentration of acylcarnitine, ketone bodies, and the esterified fatty acids in treated animals. An important finding was that metabolites highlighted by the previous analyses and related to other fatty acid derivatives were globally weakly altered (most phospholipids, unsaturated fatty acids), even if their concentration was statistically significantly different. However, some exceptions could be noted for instance for the main lysophosphatidic acid (16:0/0:0). Other clusters of metabolites were more scattered in the network and less well connected. Subnetworks were however present, such as the steroids, clearly less concentrated in exposed animals, especially in the C2 group. Similarly, disturbed pools of homocyclic aromatic compounds usually corresponding to xenometabolites were also spatially clustered. Moreover, perturbation of indole-related molecules normally connected to tryptophan and nicotinamide metabolism was confirmed through this analysis.

We performed a hierarchical clustering of our annotated metabolites associated with their measured level in each condition, the aim being to combine on the same graph the results obtained with the different concentrations of DPhP used. It is interesting to note that this analysis performed with only 633 features clearly discriminated the three exposure conditions from the control, C2 group being the most different as observed with the PCA that was based on the entire data set of detected mass (Figure S5 and Table S12). Careful inspection of the obtained clusters confirmed lower concentrations of the acylcarnitine pool and the ketone bodies in exposed animals, even with the lowest dose of DPhP and without a clear dose-dependent response (Figure 4D and S5). The ethyl-esterified form of the most common fatty acids (palmitate, linolenate, oleate) were also present in this cluster. Furthermore, we noted the same variation for trans-2-hexenedioic acid, a metabolite associated to the and the degradation of fatty acid. Conversely, a cluster highlighted a group of more concentrated metabolites in all treated animals (Figure 4E). As expected from previous analyses, N-acylethanolamides were present in this cluster as well as several polyunsaturated fatty acids (linolenic acid, linoleic acid, eicosadienoic acid). Also, several monoacylglycerols were more concentrated in liver extracts of exposed animals, most of them being polyunsaturated.

Finally, from the previous analysis, we looked for position of homocyclic aromatic compounds and indole-derivative. These compounds were more scattered in the hierarchical tree (Figure S5 and Table S12). However, sulphated derivatives were globally less concentrated in exposed animals (phenol sulfate, -cresol sulfate, indoxyl sulfate). At the opposite, the picture for nonprocessed forms of homocyclic aromatic compounds and indole-related compounds was more complicated. Phenol and 1,4,5 naphthalene-triol displayed for example a classic dose–response relationship (Figure 4D). Several indole-related compounds behaved differently between low C1/C2 dose and high C3 dose where they were markedly less concentrated (indole-3-carbinol, indole-3-acetamide, indole-3-carboxaldehyde, quinolin-2-ol). At the opposite, 5-hydroxyindole and indole-3-ylacetaldoxime were homogeneously less concentrated in exposed animals in comparison of the control.

In a final experiment, the robustness of observed differences was assessed by comparing observed fold change between the different analytical chromatography techniques. This analysis focused on acylcarnitines because they were the most altered metabolites. A strong correlation was evident between both chromatographic systems, globally for all annotated metabolites and more specifically for acylcarnitine (Figure S6). An important finding is that the same results were observed for metabolites connected to acylcarnitines such as fatty acids and ketone bodies (Figure S7).

Biological Effects of DPhP Exposure through Transcriptomic Analyses of the Liver

The second part of the dual-omics analyses was performed on four distinct batches of 10 animals untreated or treated chronically with one of the three DPhP concentration C1, C2, or C3 during 12 wk. mRNAs were extracted from pooled liver of each group and analyzed by NGS using four analytical replicates. Genes were then filtered, and those displaying a mean in the control condition were conserved (Table S13). PCA was performed using the mean expression associated with the four experimental conditions. When the first axis was considered, the three treated conditions were significantly different from the control, the most discriminating treatment being associated with the intermediate C2 concentration of DPhP (Figure 5A, left panel). Interestingly, this pattern was correlated with the metabolic difference observed, because this concentration had the greatest effects. We noticed that 44% of the total variance was attributable to the first axis (Figure 5A, right panel), indicating that the genes with the highest Eigen values on this axis were the most relevant for the biological effects of DPhP. We consequently selected these most discriminating genes (Table S14, ) and performed a gene ontology analysis using the STRING software (version 11). For the analysis, genes were ranked according to their Eigen value for generating a functional enrichment score/false discovery rate, then used to construct a volcano plot (Figure 5B and Table S15). Among the highly depleted processes in the exposed vs. control conditions, those related to lipid metabolism and more specifically to fatty acid oxidation represented the largest part of altered functions. However, we also noticed a significant inhibition of genetic responses associated with xenobiotic metabolism. Moreover, when we retrieved the gene list associated with the fatty acid catabolic process, and selected the most discriminating genes in our dataset (), we confirmed that these genes encoded functional protein networks related to mitochondrial and peroxisomal fatty acid oxidation with a very high confidence rate (Figure 5C). Since fatty acid oxidation is strongly regulated by the PPAR transcription factor, we verified that PPAR signaling was also among the significant terms associated with our ranked gene list (Figure 5D). The network reconstructed from this last term demonstrated that target genes were at the heart of the dysregulation process.

Figure 5.

Transcriptomic analysis of livers belonging to animals exposed to DPhP. (A) Principal component analysis (PCA) was performed using pareto algorithm as an observatory method to discriminate group of sample replicates based on the amount of mRNA expression obtained from indicated animals, through reverse transcription and next generation sequencing (NGS). Considered genes and explained variance with the two first axes of the PCA have been indicated. (B) Gene ontology analysis of the most discriminant genes () used to build the first axis of the previous PCA performed through functional enrichment analysis (STRING Consortium 2019). Results have been presented as a volcano plot of the significant discriminating functions associated with this principal component. Term related to lipid oxidation (red square), lipid metabolic processes (blue triangle), and xenobiotics metabolism (yellow diamond) have been highlighted as indicated. The highest significant functions belonging to these terms have been listed in the right-hand panel. (Cutoff for ). (C) Identical genes with a significant Eigen value () and overlapping the indicated GO term were used to build a protein–protein interaction network (STRING© Consortium 2019) with a high level of confidence setting. Genes belonging to a particular organelle network have been highlighted in blue (mitochondria) and red (peroxisome). Numbers of edges and nodes have been indicated, as well as Protein-Protein-Interaction (PPI) enrichment (see “Methods”). (D) Identical genes with a significant Eigen value () and overlapping the indicated GO term were used to build a protein–protein interaction network (STRING© Consortium 2019) with a high level of confidence. Genes belonging to a particular network of organelles have been highlighted in blue ( specific target genes), green ( specific target genes) and red (any PPAR target gene). Number of edges and node numbers have been indicated, as well as PPI enrichment (see “Methods”).

Next, we verified for each treatment dose that the same trends were observed by comparing their RPKM values to those of the control. In this case, we directly used the calculated fold change of genes with an . Based on an identical approach combining a volcano plot and a protein network reconstruction for each analysis, we observed that genes related to fatty acid oxidation and lipid catabolism were expressed at lower levels in treated animals, independently of the doses used (Figure S8A–C and S9A–C, Tables S16-S18). The enrichment scores of the observed differences increased according to DPhP doses, whereas the density of the protein network encoded by genes altered and related to lipid catabolism was the highest with the intermediate concentration (edge and node numbers, see methods).

Finally, to exclude a possible artifact arising from poor analytical replicates, we performed a PCA and an OPLS-DA analysis with the individual RPKM values of the 16 analytical samples (instead of the mean expression for each group) and compared each treated group with control animals. Adequate separation between replicates was obtained by using the second, the fourth, and the fifth axis of the PCA (Figure 6A). For each principal component, we retrieved the Eigen values and performed a gene ontology analysis with these scores (Tables S19-S21). Using the second axis in which C2 is the most distant from the control, this analysis confirmed the lower expression of genes controlling fatty acid catabolism in samples arising from animals exposed to DPhP compared with those exposed to vehicle (Figure 6B). An interesting result is that, by analyzing the two other components, we noticed that the lowest DPhP concentration also led to the higher expression of genes related to tryptophan and xenobiotic metabolism (Figure 6C).

Figure 6.

Transcriptomic functions significantly associated with each group of exposed animals. (A) principal component analysis (PCA) was performed using pareto algorithm as an observatory method to discriminate individual sample replicates based on the amount of mRNA expression obtained from indicated animals, through reverse transcription and Next-Generation-Sequencing. Individual replicates could be grouped efficiently through the indicated PCA axis. (B–D) Gene ontology analysis of the most discriminating genes () used to build the indicated axis of the previous PCA performed with STRING. Results have been presented as a volcano plot of the significant discriminating function for the indicated axis. Note that axis 2, axis 4, and axis 5, respectively, correspond to this order of gene expression, , , and . Cutoff for for (B), 2/0.05 for (C), and 2/0.05 for (D). (E) Orthogonal Projections to Latent Structures-Discriminant Analysis (OPLS-DA) analysis was performed to discriminate individual sample replicates based on the amount of mRNA expression obtained from indicated animals, through reverse transcription and Next-Generation-Sequencing, and in a supervised manner. A 3D dot plot has been presented showing the efficient replicate clustering and the discriminating power attributed to the first three predictors. (F) Using the predictor (x axis) that classified the samples in a classical dose–response, a bipartite view of enrichment network based on the 2,000 most discriminating genes is represented. Small dots and large dots represent individual genes and enriched functions (KEGG based), respectively. Genes or functions have been encoded by a green-red continuous scale according to their fold change or their -value, respectively, as well as the direction (repressed–activated) of the regulation. Dot size of the functions represents the percentage of genes used and matching the full list associated to this function in the KEGG database. Three clusters of functions are highlighted: , , and acute phase response. Note: CTL, control.

Finally, we observed that the highest concentrations of DPhP also induced a stronger expression of genes associated with the synthesis of sterol-derived metabolites, indicating a more complex alteration of lipid metabolism (Figure 6D). To improve sample clustering, we then tested by OPLS-DA all possible combinations able to discriminate our four groups of replicates, using the control replicates as negative controls. Four predictive components were determined (Table S22), three of which were plotted as indicated in the Figure 6E and confirmed the reliable separation of the different groups of replicates. The predictive component 3, classifying the samples in a dose–response manner was then used to build a bipartite view of the enrichment network associating the 2,000 most relevant genes to this component and the KEGG function related to these genes (Figure 6F and Tables S23, S24). The network highlighted two interconnected clusters of genetic programmes repressed in a dose-dependent manner by DPhP, namely fatty acid metabolism and the xenobiotic response. Conversely, one cluster was related to weakly activated functions belonging to the control of mRNA processes and the acute phase response. The fatty acid cluster encompassed highly significant programmes related to fatty acid oxidation, peroxisome, and PPAR transcriptional response (Figure 6F). Moreover, fatty acid biosynthetic processes such as fatty acid elongation also contributed to this large cluster, these latter functions being inhibited by DPhP treatment. Similarly, the xenobiotic cluster included several types of xenobiotic responses, such as those related to the cytochrome P450, the glucuronidation through the aldarate metabolism, or glutathione metabolism. Last, although not present directly in these clusters, we confirmed that the genetic program controlling the tryptophan metabolism was repressed by increasing doses of DPhP and established the existence of a genetic connection (Maob, Hadh, Aldh3a2, etc) between this response and the two previously mentioned repressed clusters.

Liver Protein Expression and Physiological Analysis in Mice Treated Chronically with DPhP

To further investigate the observed abnormalities in the processes related to lipid catabolism and signaling, we performed a series of IHC analyses revolving around key enzymes of liver physiology, involved in fatty acid catabolism and ketone body formation. In addition, we stained the lipid droplets using Perilipin 2 staining, a protein surrounding these vesicles and correlated with their abundance (Straub et al. 2013). Furthermore, we analyzed enzymes controlling the gluconeogenesis to verify that an overall nonspecific dysregulation of liver metabolism did not appear following DPhP treatment (Figure 7A). Expression of phosphoenolpyruvate carboxykinase1 (PCK1), an important enzyme for gluconeogenesis was thus not lower and was even higher in some animals after DPhP exposure. Conversely, the expression of Hmgcs2, a kernel target gene, controlling the ketone bodies synthesis from fatty acid oxidation was strongly inhibited compared with control, even at the lowest dose of DPhP (Figure 7A). Therefore, only specific enzymes were altered after DPhP exposure. This inhibition occurred more specifically in the peri-portal area, where Hmgcs is normally more strongly expressed and active (Figure 7B and S10A). At the highest dose, 100% of animals presented this physiological alteration. Interestingly, this weaker expression of HMGCS2 was highly correlated with a lower amount of lipid droplets in animals exposed to of DPhP (Figure 7A–C and S10B). However, this was not the case for the lowest doses where the amount of Perilipin 2 staining was more intense upon exposure to DPhP, demonstrating here that lower levels of acylcarnitine in these conditions was not the consequence of a general depletion in lipid stores. Counterintuitively, staining was also uncorrelated with Hmgcs2 expression because exposed animals displayed a stronger nuclear and cytosolic expression (Figure 7A).

Figure 7.

Histological and physiological alterations induced by exposure to DPhP. (A) Liver sections at the indicated scale of mice exposed to the indicated concentration of DPhP or a vehicle were immunostained with the indicated antibodies. Star and hash symbols denote the centrolobular and the portal area of the liver, respectively. Inset represents a magnification of the original image. (B) Histogram quantifying immunohistochemistry (IHC) scores associated with Hmgcs staining of five animals in each indicated group (representative of two independent experiments, ). Score associated with the centrilobular and the portal area were dissociated as indicated. , , , (see “Methods”). (C) Histogram quantifying IHC scores associated with Perilipin 2 staining of five animals in each indicated group (representative of two independent experiments, ). Score associated with the intermediate zone and the portal area were dissociated as indicated (no staining was present in the zone contiguous to the centrilobular vein). , , , (see methods). (D) Box-and-whisker plot quantifying the body weight of 10 animals representative of two independent experiments (), after 8 wk of DPhP exposure at the indicated concentrations or with a vehicle. (, and 3rd quartile, and smallest value within 1.5 × the interquartile range; summary data in Table S27). Outliers and significant -values have been indicated (Student’s -test, compared to vehicle).

Finally, because all of our results highlighted a disturbance in lipid metabolism, we collected the data obtained from animals used for the metabolomics and the transcriptomic analyses to construct the growth curve of these animals according to their overall weight gain. No significant difference was observed at the lowest doses, but a trend for weight loss was clearly evident in a dose-dependent manner (Figure 7D). Moreover, weight was significantly lowered in animals treated with the highest dose of DPhP.

Discussion

Only few studies have been conducted to directly test the effect of DPhP, the most common APE derivative in human samples. Our results showed that at least in mice, DPhP levels in biological fluids are unlikely to represent a surrogate of direct APE ingestion. We demonstrated that even through direct and acute exposure by IV or oral gavage, only a minor fraction of a parent APE such as TPhP were converted in vivo into DPhP. Moreover, during an overnight exposure to low TPhP concentrations in drinking water, conversion of TPhP into DPhP could no longer be measured. In vivo DPhP may thus be a surrogate of spontaneous degradation of APEs in the environment, suggesting that experimental procedures revolving around DPhP are likely more relevant for assessing APE toxicity.

Using a dual-omics analyses approach based on this strategy, our results demonstrated that in mice, DPhP exposure, even at doses close to the lowest doses encountered in the environment, disturbed important parameters associated with the control of liver metabolism. A substantial body of evidence consistently pointed toward a dysregulation of lipid catabolism. As a key result, we observed an important perturbation, in terms of intensity and distribution, of the acylcarnitine pools, the main metabolites associated with the degradation of fatty acids by . Moreover, additional modifications could be connected to this perturbation, reinforcing the likelihood that DPhP induced alteration of lipid metabolism. The most important of them was the simultaneous reduction of ketone bodies, especially (R)-3-hydroxybutyric acid, the final product of this synthetic pathway. Indeed, ketone bodies are typically produced from fatty acid oxidation during fasting (Grabacka et al. 2016), their diminution supporting a lower activity of this pathway. In addition, small dicarboxylic acids such as the trans-2-hexenedioic acid, known for a long time to be excreted in correlation with fatty acid oxidation activity (Tserng and Jin 1991), were also less concentrated in liver extracts of exposed animals. Finally, even if we could not confirm a significant modification in the main endogenously produced free fatty acids, we noted that acylcholine and ethyl-fatty acids, two of their esterified forms produced nonenzymatically from available fatty acids, were significantly depleted in liver extracts of exposed animals. Similarly, in phosphatidylcholine, the occurrence of even fatty acid chains produced endogenously was reduced in favor of odd and polyunsaturated carbon chain. These observations suggested that fatty acid production was also likely to be altered.

Aside from metabolic evidence, our results demonstrated that the genetic programs associated with the oxidation of fatty acids were also significantly repressed, even at the lowest DPhP concentration. target genes, such as Hmgcs2 or Cpt1a were specifically less expressed in exposed animals. It is important to note that these genes are directly connected to the regulation of perturbed pool of metabolites, Hmgcs2 controlling the synthesis of (R)-3-Hydroxybutyric acid (Grabacka et al. 2016), and Cpt1a, catalyzing acylcarnitine production and rate limiting for fatty acid oxidation (Foster 2012). Finally, reduced levels of fat in the liver and smaller BW could be observed in exposed animals, especially at the highest dose of DPhP. In line with this result, the transcriptional control of fatty acids synthesis/uptake was apparently also disturbed, especially at the highest dose of DPhP. Consistently, Cd36, Scd1, Fads2, all canonical target genes of involved in fatty acid biosynthesis, were found to be less expressed in exposed animals. Taken together, these results strongly pointed toward a significant dysregulation of lipid homeostasis upon exposure to DPhP and suggested an important role for activity in this perturbation.

is considered to be a key coordinator of fast-fed transition at the hepatic level with paradoxical effects (Dubois et al. 2017). is thus the main activator of fatty acid oxidation and ketogenesis during adaptation to long-term fasting (Janssen et al. 2015). However, it is now recognized that also acts in normal conditions to adapt metabolism to circadian feeding (Gachon et al. 2011). Therefore, controls on the one hand de novo lipid synthesis and fatty acid uptake during feeding periods, supplying store droplets, and on the other hand activates triglyceride and cholesteryl-ester lysis and oxidation of released lipids between meals (Chen and Yang 2014; Gachon et al. 2011). Accordingly, our results showing a double perturbation of both mechanisms, fatty acid synthesis and degradation, by DPhP exposure were consistent with this known role of . In absence of any observed reduction of expression, our results match a model of inhibition of transcriptional activity upon exposure to DPhP. Conversely, our findings that lipid stores in the liver of animals exposed to low concentrations of DPhP are not depleted invalidated a model of low fatty acid oxidation by an overall lack of lipids and was consistent with a genetic perturbation of this metabolic pathway. Finally, the absence of a strong rise in lipid droplets in these same conditions was also coherent with an inhibition of activity, because an equilibrium may be found between the lack of fatty acid storage and a lower use of these fatty acids.

Interestingly, our study highlighted the perturbation of several pools of metabolites and genetic controls susceptible to alter this equilibrium, directly in the liver or involving a more complex interaction with other organs. Among them, N-acylethanolamide, significantly more represented in extracts of animals exposed to DPhP, are known lipid stress markers enzymatically produced and controlling activity (Piomelli 2013). Similarly, depletion of bile acids may participate in the control of lipid metabolism, either directly through a decrease in lipid absorption, or via the regulation of nuclear receptors such as FXR protein. Finally, we noticed the alteration of several metabolites, endogenous and exogenous, containing homocyclic aromatic rings and/or indole chemical functions associated with tryptophan metabolism. Several perturbed indole-containing compounds (indoxyl-sulfate, quinolin-2-ol) (Wikoff et al. 2009) and tryptophan derivatives (Venkateswaran et al. 2019; Yamamoto et al. 2019) were shown to be directly implicated in the regulation of AHR protein, an important xenobiotic receptor (Bock 2019b) involved in the elimination of aromatic xenobiotics. Moreover, even though DPhP does not contain the classical polycyclic structure associated with known ligands of these receptors, it possesses two aromatic rings susceptible of stacking with these structures and disturbing their interaction with AHR (Soshilov and Denison 2014). Importantly, we also noticed at the genetic level that exposure to DPhP disturbed xenobiotic response and tryptophan metabolism reinforcing the possibility that DPhP disturbed AHR activity in our experiments. Recent studies have suggested that AHR may be involved in fatty acid metabolism (Bock 2019a), independently of its role as a xenobiotic regulator but through the control of activity (Wang et al. 2011).

Overall, these regulations linked to exposure to DPhP may be highly relevant to human health with regards to long-term exposure. Of note, our study used relatively young animals and did not analyze long-term consequences of chronic exposure to DPhP. Therefore, it is difficult to predict the exact consequences of inhibition by DPhP on older animals. Moreover, it is likely that this type of perturbation intersects with other environmental constraints such as the type of diet susceptible to tilt the equilibrium controlled by in one direction or the other. Finally, the dose of DPhP exposure may strongly impact the final outcome of perturbation in an inverted U-shape type of dose–response, because we have seen in mice that the intermediate concentration had a stronger impact on liver metabolism than the highest. However, several data could support a role for inhibition by DPhP in the development of lipid disorders. The lack of activity leads for instance to a higher amount of plasma triglyceride and LDL (Schoonjans et al. 1996; Watts et al. 2003) and is the target of fibrates, a therapeutic class of molecules extensively used to treat human metabolic syndromes. Moreover, transgenic mice with a liver-specific knockout were more sensitive to high-fat diet, displaying specific inflammation and hyperlipidemia at the hepatic level (Stec et al. 2019). Therefore, inhibition in individuals exposed to DPhP may contribute to the development of metabolic syndromes and their consequences, such as cardiovascular disease and even cancer.

Strengths and Limitations

As described in the “Methods” section, we used pooled samples of livers for the omics analyses. This approach was chosen due to technical reasons with the aim of obtaining sufficient material for each group/concentration to perform the multiple analytical experiments used for the metabolomics analyses. The robustness of our analytical procedure and our metabolic coverage were strongly improved by this decision. We assumed that a divergent metabolite profile of one animal could overlap those of others pooled in the same extraction procedure. However, to compensate this limitation, we performed our transcriptomic analyses on a distinct batch of animals. A strong correlation and integration could be observed and built between both types of omics analysis. Therefore, it was extremely unlikely that observed differences in both omics analyses were the consequences of some animals presenting highly divergent biological responses. Moreover, immunohistochemistry (IHC) analyses were performed on separate animals and confirmed the homogeneous response to the DPhP treatment.

Conclusion

Our results raised many questions on the use, the safety and the presence in the environment of APEs, most of which may expose humans to DPhP. We did not fully characterize the molecular mechanisms underlying the apparent alterations in lipid homeostasis and activities by these compounds because this was beyond the scope of our study. However, the known function of these factors and their association with the metabolic syndromes should constitute a sufficiently strong risk factor to measure more precisely the health hazards associated with their presence in the environment, by taking into consideration diet in future epidemiological studies. Moreover, the potential interaction of these compounds with known activators of AHR should be investigated with the aim of determining the possible existence of synergistic effects and to characterize the mechanisms, directly or indirectly, inhibiting activities.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.