Proteomics and metabolic phenotyping define principal roles for the aryl hydrocarbon receptor in mouse liver.

Dioxin-like molecules have been associated with endocrine disruption and liver disease. To better understand aryl hydrocarbon receptor (AHR) biology, metabolic phenotyping and liver proteomics were performed in mice following ligand-activation or whole-body genetic ablation of this receptor. Male wild type (WT) and Ahr -/- mice (Taconic) were fed a control diet and exposed to 3,3',4,4',5-pentachlorobiphenyl (PCB126) (61 nmol/kg by gavage) or vehicle for two weeks. PCB126 increased expression of canonical AHR targets (Cyp1a1 and Cyp1a2) in WT but not Ahr -/-. Knockouts had increased adiposity with decreased glucose tolerance; smaller livers with increased steatosis and perilipin-2; and paradoxically decreased blood lipids. PCB126 was associated with increased hepatic triglycerides in Ahr -/-. The liver proteome was impacted more so by Ahr -/- genotype than ligand-activation, but top gene ontology (GO) processes were similar. The PCB126-associated liver proteome was Ahr-dependent. Ahr principally regulated liver metabolism (e.g., lipids, xenobiotics, organic acids) and bioenergetics, but it also impacted liver endocrine response (e.g., the insulin receptor) and function, including the production of steroids, hepatokines, and pheromone binding proteins. These effects could have been indirectly mediated by interacting transcription factors or microRNAs. The biologic roles of the AHR and its ligands warrant more research in liver metabolic health and disease.

Introduction

The liver is the largest internal organ. It is responsible for myriad detoxification and synthetic processes which protect and nourish the body. Likely due to its prominent roles in xenobiotic and intermediary metabolism, the liver is the principal target of toxicities from alcohol, pharmaceuticals, environmental pollutants and obesity. The liver-related death rates from cirrhosis and liver cancer are rapidly increasing. Nonalcoholic fatty liver disease (NAFLD) is the most prevalent liver disease worldwide, and its more severe form is called non-alcoholic steatohepatitis. Although typically considered to be the hepatic manifestation of obesity and metabolic syndrome, NAFLD may also be caused or modulated by pharmaceuticals and environmental pollutants including some endocrine and metabolism disrupting chemicals2, 3, 4, 5, 6. The terms toxicant-associated fatty liver disease (TAFLD) and toxicant-associated steatohepatitis (TASH) have been proposed to describe the latter situation,.

Many liver physiologic and pathophysiologic processes are regulated by ligand-activated transcription factors, including the aryl hydrocarbon receptor (AHR). The AHR is a well-established master regulator of xenobiotic metabolism. Its canonical target genes include the cytochrome P4501A (Cyp1a) family, also implicated in carcinogenesis. AHR's high-affinity binding ligands include dioxins (e.g., 2,3,7,8-tetrachlorodibenzo-p-dioxin) and other dioxin-like molecules (e.g., 3,3′,4,4′,5-pentachlorobiphenyl also called polychlorinated biphenyl 126 or PCB126). The high thermodynamic stability of these persistent organic pollutants makes them resistant to degradation leading to sustained AHR activation. AHR's low-affinity ligands include rapidly metabolized dietary/endogenous molecules such as: flavonoids, bilirubin, and gut microbiome-derived tryptophan metabolites. Apart from its regulatory role in hepatic xenobiotic metabolism and carcinogenesis, the AHR has also been implicated in intermediary metabolism, metal homeostasis,, fibrosis and the production of hepatokines and liver-derived circulating pro-atherogenic molecules,. Given AHR's role in the regulation of intermediary metabolism, it is not surprising that dioxins and dioxin-like molecules have been associated with endocrine and metabolic disruption. Our group and others have been studying these foreign compounds in TAFLD/TASH,,,,16, 17, 18, 19.

The role of the AHR in hepatic lipid metabolism and obesity-related diseases like NAFLD is complex. Complicating matters, this role may be ligand-dependent and species/strain-dependent. While some dioxins and dioxin-like molecules may cause TAFLD, some endogenously-produced AHR ligands may be protective against at least some aspects liver disease. Moreover, while some strains of Ahr–/– mice spontaneously developed steatosis, others were protected against diet-induced fatty liver disease20, 21, 22, 23, 24. While a proteomics approach has previously been utilized to investigate PCBs in TAFLD, we could find no liver proteomics data for Ahr knockout mice in the published literature. Therefore, in this manuscript we performed paired liver proteomics and metabolic phenotyping in wild type and whole-body Ahr knockout mice with or without PCB126-induced ligand-activation of the AHR. The top liver pathways associated with these interventions were elucidated. The results will enhance the current understanding of AHR's principal roles in liver health and metabolic diseases.

Materials and methods

Animal studies

The animal protocol was approved by the University of Louisville Institutional Animal Care and Use Committee (Louisville, KY, USA). Adult male C57BL/6 mice (wild type, WT, catalog number: B6-M, 9 weeks old) and Ahr–/– mice (catalog number: 9166-M, 8–9 weeks old) were purchased from Taconic Biosciences Laboratory (Hudson, NY, USA). While additional strains of Ahr–/– mice have been generated by other laboratories, this strain was chosen because it is a whole-body knockout model generated on a C57BL/6 background, and it was commercially available. Importantly, C57BL/6 mice have a high-affinity AHR. Taconic specifically recommend the B6-M WT strain as their most appropriate control for the 9166-M knockout model, although B6-M is not a litter-mate control. Notably, our co-authoring toxicologist, Dr. Wahlang, previously published on the hepatotoxicity of PCB126 in Taconic's B6-M WT strain, albeit at a higher PCB126 dose and different diet. That manuscript documented induction of the canonical AHR target gene, Cyp1a1, and metabolic disruption by PCB126.

All mice were fed a control synthetic diet (20.0%, 69.8%, and 10.2% of total calories from protein, carbohydrate, and fat; TekLad TD 06416). At 9–10 weeks of age, the mice were administered either vehicle control (corn oil) or 20 μg/kg PCB126 (61 nmol/kg) via a one-time gavage for two weeks. The PCB126 dose and duration of exposure is justified in our prior publications,,. Briefly, this dose is lower than that typically used by other investigators, and we believe it to be relevant to human exposures while still activating key canonical AHR target genes in mice (e.g., Cyp1a1 and Cyp1a2). These procedures generated four groups of mice (n = 10 mice) which were designated as: the WT Vehicle group; WT PCB126 group; the Ahr–/– Vehicle group and the Ahr–/– PCB126 group. The study design is summarized graphically in Supporting Information Fig. S1. After two weeks, the mice were fasted overnight and dual energy X-ray absorptiometry scanning (Lunar PIXImus densitometer, WI, USA) was performed to analyze body composition. Tissues (e.g., liver, adipose, plasma, etc.) were then collected following administration of ketamine/xylazine (120/16 mg/kg body weight) given by intraperitoneal injection. Plasma was collected with EDTA as the anticoagulant.

Histological staining

Liver and adipose tissues were fixed in 10% neutral buffered formalin for 72 h and embedded in paraffin for routine histological examination. Hematoxylin-eosin (H&E) staining was performed to identify histopathological changes. In order to better evaluate histologic steatosis, Oil Red O stain was done in tissue that was placed in optimal cutting temperature reagent at the time of necropsy and snapped frozen in liquid nitrogen. Micrographic images were acquired by a high-resolution digital scanner (Olympus) with a digital camera (Olympus BX41).

Real-time PCR

Mouse liver tissues were homogenized and total RNA was extracted using RNA-STAT 60 (AMSBIO, Cambridge, MA, USA) according to the manufacturer's protocol. The purity and quantity of total RNA were assessed with a Nanodrop spectrometer (ND-1000, ThermoFisher Scientific, Waltham, MA, USA) using ND-1000 V3.8.1 software. cDNA was reverse transcribed from 1 μg RNA with a one-step cDNA synthesis reagent (QScript cDNA Supermix, QuantaBio, Beverly, MA, USA). Then RT-PCR was performed on the CFX384TM Real-Time System (Biorad, Hercules, CA, USA) using iTaq Universal probe Supermix and Taqman probes as described previously. All reactions were performed in triplicate. The relative mRNA expression was calculated using the comparative 2–ΔΔCt method and normalized against GAPDH mRNA. MicroRNAs (miRs) were isolated from mouse livers using the MagMAX™ mirVana™ Total RNA Isolation Kit (ThermoFisher Scientific) according to the manufacturer's protocol. Hepatic expression for selected miRs was measured using RT-PCR, like the gene expression (mRNA) method but using Taqman miR probes (ThermoFisher Scientific).

Measurement of hepatic lipids, plasma lipids and cytokines

The liver tissues were rinsed in 1× phosphate buffered saline and homogenized in 50 mmol/L NaCl solution. Hepatic lipids were extracted by a mixed solution of chloroform and methanol (2:1) according to a published protocol. Hepatic triglycerides and free fatty acids were assessed using commercial kits with final values normalized to liver weight. Plasma alanine transaminase (ALT), aspartate transaminase (AST), cholesterol, triglyceride, high-density lipoprotein (HDL), low-density lipoprotein (LDL), very low-density lipoprotein (VLDL), and non-HDL cholesterol (nHDLc) levels were determined with lipid panel plus kits on a Piccolo Xpress Chemistry Analyzer (Abbott Laboratories, Chicago, IL, USA). Plasma cytokine and adipokine levels were evaluated using a customized Milliplex MAP mouse adipokine panel (Millipore Sigma, Billerica, MA, USA) on a Luminex 100 system (Luminex Corp, Austin, TX, USA).

Proteomics analysis

Proteins were extracted from liver tissues in 1% SDS modified RIPA buffer with protease and phosphatase inhibitors, using a bead homogenizer, and protein amounts were measured by BCA protein assay. Protein lysates (200 μg) were trypsinized using the modified filter-aided sample preparation method. Protein samples were first reduced by dithiothreitol, denatured by 8 mol/L urea and alkylated by iodoacetamide, followed by centrifugation through a high molecular weight cutoff centrifugal filter (Millipore Sigma, 10k MWCO). Next, after overnight digestion with sequencing grade trypsin (Promega, Madison, WI, USA) at 37 °C, the digested peptides were collected and cleaned with a C18 Proto™ 300 Å ultra microspin column. Digested peptide samples (50 μg) were labeled with tandem mass tag (TMT) TMT10plex™ isobaric label reagent set (ThermoFisher Scientific). Next, they were concentrated and desalted with Oasis HLB extraction cartridges (Waters Corporation, Milford, MA, USA) using a modified protocol for extraction of the digested peptides. Samples were separated by high pH reversed phase separation with fraction concatenation on a Beckman System Gold LC system supplemented with 126 solvent module and 166 UV–Vis detector in tandem with a BioRad Model 2110 Fraction Collector.

Liquid chromatography/mass spectrometry was used to measure TMT-labeled peptides. Briefly, every high pH reversed phase fraction was dissolved in 50 μL solution of the combination of 2% v/v acetonitrile with 0.1% v/v formic acid. 1 μL of each fraction was analyzed on EASY-nLC 1000 UHPLC system (ThermoFisher Scientific) and an Orbitrap Elite–ETD mass spectrometer (Thermo Fisher Scientific). The raw data from the mass spectrometer were analyzed by Proteome Discoverer v2.2.0.388.

Statistical analysis and data sharing

Statistical evaluation was performed by two-way analysis of variance (ANOVA) using GraphPad Prism version 7.02 for Windows (GraphPad Software Inc., La Jolla, CA, USA), and Tukey's test was used as a post hoc test. The two factors analyzed were mouse genotype (genotype effect) and PCB126 exposure (PCB effect). An interaction effect was also determined. Results are reported as mean ± standard deviation (SD). P < 0.05 was considered statistically significant.

Statistical analysis for the proteomic data was carried out with R packages using a modified version of our biostatistician's previously published protocol. First, the raw data were transformed by taking logarithmic base 2 followed by quantile normalization. Then, the missing values were imputed using singular value decomposition method. Proteins with missing values > 40% were excluded from subsequent analysis. Finally, differentially abundant proteins (P < 0.05) were further filtered by fold-change (FC) criteria (−1 < log2FC < 1) and multiple comparisons testing with a false discovery rate (FDR) <0.05 yielding the final significant results for the differentially abundant proteins. Next, these proteins were imported into MetaCore software (Clarivate Analytics, Philadelphia, PA, USA) and the following analyses were performed: gene ontology (GO) processes, enrichment by protein function (EPF), and interaction by protein function (IPF). Proteomics data files were deposited with MassIVE (http://massive.ucsd.edu/) data repository, Center for Computational Mass Spectrometry at the University of California, San Diego and shared with the ProteomeXchange (www.proteomexchange.org).

Results

Impact of Ahr genotype and PCB126 exposure on the hepatic expression of key xenobiotic receptors and their target genes

The mRNA expression levels of Ahr, Pxr and Car, as well as their target genes were assessed by RT-PCR. These studies confirmed the absence of Ahr expression in the knockout mice (Fig. 1A). Consistent with this finding, expression of the canonical AHR target genes, Cyp1a1 and Cyp1a2, was significantly decreased in Ahr–/– (genotype effect, Fig. 1B and C). In WT mice, PCB126 exposure increased Cyp1a1 (∼2610-fold) and Cyp1a2 (∼10-fold) expression vs. vehicle consistent with AHR activation. Ahr–/– mice had significantly increased mRNA expression of Pxr (genotype effect, Fig. 1D) and its target gene, Cyp3a11 (genotype effect, Fig. 1E). Car expression was down-regulated in Ahr–/– mice and up-regulated by PCB126 exposure (genotype effect and PCB effect, Fig. 1F). However, expression of the CAR target gene, Cyp2b10, was paradoxically increased in Ahr–/– mice (genotype effect, Fig. 1G). Thus, Ahr–/– mice had deleted Ahr; increased Pxr expression and activity; and increased Car activity despite reduced Car expression. PCB126 activated the AHR, and to a lesser degree, CAR.

Figure 1

Hepatic expression of key receptors regulating xenobiotic metabolism and selected target genes. Expression of mRNA was assessed by RT-PCR with results (mean ± SD) normalized to the WT Vehicle group. (A) Ahr; (B) Cyp1a1 (AHR target gene); (C) Cyp1a2 (AHR target gene); (D) Pxr; (E) Cyp3a11 (PXR target gene); (F) Car; and (G) Cyp2b10 (CAR target gene). A complete list of P-values (determined by two-way ANOVA with Tukey's post-test) is provided in the accompanying table. In the figure panels, statistical significance is denoted by: a = genotype effect; b = PCB effect; c = interaction effect; d = WT Vehicle vs. WT PCB126; e = Ahr–/– Vehicle vs. Ahr–/– PCB126; f = WT Vehicle vs. Ahr–/– Vehicle; g = WT PCB126 vs. Ahr–/– PCB126. AHR, aryl hydrocarbon receptor; CAR, constitutive androstane receptor; CYP, cytochrome P450; PCB, polychlorinated biphenyl; PXR, pregnane xenobiotic receptor; WT, wild type.

Impact of Ahr genotype and PCB126 exposure on body composition, blood glucose and lipids

Body weight and composition varied by experimental group. Body weights over time are given in Supporting Information Fig. S2A. The Ahr–/– mice had significantly lower end of study body weights (genotype effect, Fig. S2B). However, because these mice also weighed less at the start of the study (Fig. S2A), the % change in body weight did not differ between the Ahr–/– and WT mice (Fig. 2A). However, PCB126 exposure was associated with reduced % change in body weight (PCB effect, Fig. 2A). Body composition varied according to genotype. The total % body fat (genotype effect, Fig. 2B) and the white adipose tissue to body weight ratio (genotype effect, Fig. S2C) were significantly increased in Ahr–/– mice, while % lean body mass was reduced (genotype effect, Fig. 2C).

Figure 2

Metabolic phenotyping. Body composition was determined by (A) % change in final total body weight relative to the initial total body weight; and using dual energy X-ray absorptiometry scan to obtain measurements for (B) % body fat; and (C) % lean body mass. (D) Fasting plasma glucose levels was determined by Piccolo Xpress chemical analyzer; and (E) Plasma insulin level was measured by Luminex® 100 system. A glucose tolerance test (GTT) was performed and (F) curves for GTT were plotted; in addition, (G) the area under curve (AUC) was determined. Circulating lipids, namely, fasting plasma (H) total cholesterol; (I) HDL cholesterol; and (J) triglycerides were measured using Piccolo Xpress chemical analyzer. Values are mean ± SD. A complete list of P-values (determined by two-way ANOVA with Tukey's post-test) is provided in the accompanying table. In the figure panels, statistical significance is denoted by: a = genotype effect; b = PCB effect; c = interaction effect; d = WT Vehicle vs. WT PCB126; e = Ahr–/– Vehicle vs. Ahr–/– PCB126; f = WT Vehicle vs. Ahr–/– Vehicle; g = WT PCB126 vs. Ahr–/– PCB126. Ahr, aryl hydrocarbon receptor; PCB, polychlorinated biphenyl; WT, wild type.

Regarding biomarkers of intermediary metabolism, Ahr–/– genotype was associated with significantly decreased fasting blood glucose (genotype effect, Fig. 2D) and insulin (genotype effect, Fig. 2E). Curves generated from the glucose tolerance test and area under the curve (AUC) data are provided in Fig. 2F and G. PCB126 exposure was associated with improved glucose tolerance (PCB effect), with a trend towards worsened glucose tolerance in Ahr–/– mice (P = 0.16, genotype effect). However, AUC was significantly increased in Ahr–/– Vehicle vs. WT Vehicle and in Ahr–/– PCB126 vs. WT PCB126. Ahr–/– was associated with significantly reduced total cholesterol (genotype effect, Fig. 2H), HDL-c (genotype effect, Fig. 2I) and triglycerides (genotype effect, Fig. 2J) with a trend towards reduced VLDL (P = 0.06, genotype effect, Fig. S2D). These lipid parameters were not affected by PCB126. No significant differences were seen for LDL-c (Fig. S2E). Additionally, histological analysis on white adipose sections was performed (Fig. S2F), and no differences were observed in adipocyte morphology between groups. In summary, Ahr–/– genotype was associated with increased adiposity and variably worsened glucose tolerance with paradoxically decreased fasting plasma glucose, insulin and lipid levels. PCB126 exposure was associated with decreased body weight gain and improved glucose tolerance.

Impact of Ahr genotype and PCB126 exposure on liver

Representative liver histology is provided in Fig. 3A (H&E stain) and Fig. 3B (Oil Red O stain). These images show qualitatively increased hepatic steatosis in the Ahr–/– mice. Consistent with the histology, hepatic triglycerides (genotype effect, Fig. 3C) and free fatty acids (genotype effect, Fig. 3D) were increased with Ahr ablation, although hepatic cholesterol was unchanged (Fig. 3E). PCB126 exposure was associated with increased hepatic triglycerides (PCB effect, Fig. 3C), but this effect was driven solely by the comparison of Ahr–/– PCB126 vs. Ahr–/– Vehicle. Despite increased steatosis, the liver to body weight ratio was reduced in Ahr–/– mice (genotype effect, Fig. 3F). Plasma ALT (genotype effect, Fig. 3G) and AST (genotype effect, Fig. 3H) activities were increased in the knockout mice. Hepatic mRNA expression of the scavenger receptor, cluster of differentiation 36 (Cd36), and the lipase, patatin like phospholipase domain containing protein 3 (Pnpla3), were measured. Cd36 (Fig. 3I) was increased by PCB126 exposure (PCB effect) and by Ahr–/– genotype (genotype effect). Pnpla3 was reduced in Ahr–/– (genotype effect, Fig. 3J). Overall, these data show that the Ahr–/– mice had increased liver steatosis and injury which could be related, in part, to increased blood lipid uptake by hepatic CD36 receptors.

Figure 3

Characterization of fatty liver disease. Histological analysis of liver sections was performed using (A) hematoxylin-eosin (H&E) stain (10× magnification) and (B) Oil red O stain (10× magnification). The inset is 40× magnification. Hepatic lipids were extracted and measured including (C) liver triglycerides; (D) liver free fatty acids (FFA); and (E) liver cholesterol. Whole liver was isolated at euthanasia and (F) liver weight to total body weight ratio was calculated. Activity of circulating liver enzymes, namely (G) alanine aminotransferase (ALT); and (H) aspartate aminotransferase (AST) were calculated using Piccolo Xpress chemical analyzer. RT-PCR was performed to measure hepatic mRNA expression levels of (I) Cd36; and (J) Pnpla3. Expression of mRNA was normalized to the WT Vehicle group. Values are mean ± SD. A complete list of P-values (determined by two-way ANOVA with Tukey's post-test) is provided in the accompanying table. In the figure panels, statistical significance is denoted by: a = genotype effect; b = PCB effect; c = interaction effect; d = WT Vehicle vs. WT PCB126; e = Ahr–/– Vehicle vs. Ahr–/– PCB126; f = WT Vehicle vs. Ahr–/– Vehicle; g = WT PCB126 vs. Ahr–/– PCB126. AHR, aryl hydrocarbon receptor; CD36, cluster of differentiation 36; PCB, polychlorinated biphenyl; PNPLA3, patatin-like phospholipase domain-containing protein 3; WT, wild type.

Impact of Ahr genotype and PCB126 exposure on plasma adipocytokines and hepatokine mRNA expression

Despite having increased obesity and NAFLD, the Ahr–/– mice had decreased circulating levels of pro-inflammatory cytokines including: resistin (genotype effect, Supporting Infromation Fig. S3A) and interleukin 6 (IL-6, genotype effect, Fig. S3B). Plasma monocyte chemoattractant protein-1 (MCP-1, Fig. S3C) and plasminogen activator inhibitor-1 (PAI-1, Fig. S3D) were unchanged. The adipokine, leptin, was increased in Ahr–/– mice (genotype effect, Fig. S3E) consistent with the increased % body fat in these mice. These adipocytokines were not changed by PCB126.

The mRNA expression of several hepatokines previously implicated in obesity, diabetes and NAFLD was determined. The AHR suppressed fibroblast growth factor (Fgf21) transcription via binding to a xenobiotic response element within the Fgf21 promoter. Therefore, it was not surprising that Fgf21 expression was dramatically increased in Ahr–/– mice (genotype effect, Fig. S3F). Insulin-like growth factor 1 (Igf1) was slightly increased by PCB126 (PCB effect, Fig. S3G), while betatrophin was reduced by Ahr–/– (genotype effect, Fig. S3H) or PCB126 exposure (PCB effect, Fig. S3H). The observed increases in Fgf21 and Igf1 expression should theoretically mitigate obesity-associated diseases, while the reduced betatrophin could lead to decreased pancreatic beta cell mass and eventually worsened diabetes over time. In summary, PCB126 was associated with increased Igf1 and decreased betatrophin. Ahr–/– increased leptin and Fgf21 while decreasing IL-6, resistin and betatrophin.

Impact of Ahr genotype and PCB126 exposure on the hepatic proteome

Primary proteomics results

5075 unique proteins and their isoforms were detected. Ahr–/– and exposure to the AHR agonist, PCB 126, produced distinct hepatic proteomes (Supporting Information Table S1). To facilitate interpretation of intergroup comparisons, volcano plots and a Venn diagram were constructed demonstrating differential abundance of hepatic proteins (Supporting Information Fig. S4).

In WT mice, PCB126 exposure (vs. vehicle control) was associated with 8 differentially abundant proteins (5 increased and 3 decreased, Fig. 4A). PCB126 increased hepatic abundance of CYP1A1 (Log2FC = 4.82) and CYP1A2 (Log2FC = 4.36) consistent with the mRNA expression data (Fig. 1B and C). These P450s were the proteins with the largest fold-changes and greatest statistical significances in this intergroup comparison.

Figure 4

Changes in the hepatic proteome. Alterations in hepatic proteins were demonstrated by volcano plots showing log2-transformed changes in protein abundance (vs. control) on the x-axis with log10-transformed P values on the y-axis with comparisons as follows: (A) WT Vehicle vs. WT PCB126; (B) WT Vehicle vs. Ahr–/– Vehicle; (C) WT Vehicle vs. Ahr–/– PCB126; (D) WT PCB126 vs. Ahr–/– PCB126; and (E) Ahr–/– Vehicle vs. Ahr–/– PCB126. Red denotes significantly up-regulated proteins and blue denotes significantly down-regulated proteins meeting the pre-defined fold-change (≥2-fold increase or decrease in abundance) and FDR (≤0.2) thresholds. The proteins in green met the threshold for statistical significance, but they did not meet the fold-change criterion. (F) A summarized table of the differentially expressed proteins across experimental groups. (G) Venn diagram of the differentially expressed proteins in each experimental group relative to the WT Vehicle group.

Ahr–/–vs. WT genotype was associated with 340 differentially abundant proteins (279 increased and 61 decreased) in vehicle-administered animals (Fig. 4B). CYP3A11 protein (Log2FC = 3.61) was increased in Ahr–/– consistent with its mRNA expression (Fig. 1E). Cyp1A2 was increased at the protein (Log2FC = 1.91) but not the mRNA level (Fig. 1C). The up-regulated proteins with either the greatest fold-change or statistical significance included: glutathione S-transferase A1 (Log2FC = 4.97) and aldehyde dehydrogenase X, mitochondrial (Log2FC = 3.14). Down-regulated proteins with either the greatest fold-changes or statistical significance included: major urinary proteins 1 (Log2FC = −5.58) and 17 (Log2FC = −5.00); thyroid hormone-inducible hepatic protein (Log2FC = −4.45); and 3 beta-hydroxysteroid dehydrogenase type 5 (Log2FC = −4.28). The apoptosis maker, annexin A5 (Log2FC = 1.93) and the lipid droplet protein, perilipin-2 (PLIN2, Log2FC = 2.24), were increased in Ahr–/–.

Ahr–/– PCB126 vs. WT Vehicle was associated with 301 differentially abundant proteins (209 increased and 92 decreased) (Fig. 4C). Ahr–/– vs. WT genotype was associated with 262 differentially abundant proteins (185 increased and 77 decreased) in PCB126-exposed animals (Fig. 4D). CYP1A1 (Log2FC = −4.90) and CYP1A2 (Log2FC = −6.43) protein levels were reduced and CYP3A11 was increased (Log2FC = 2.38) in Ahr–/– PCB126 vs. WT PCB126 consistent with the gene expression data (Fig. 1). Glutathione S-transferase A1 was increased to the greatest degree (Log2FC = 5.13), while major urinary protein 1 was decreased to the greatest degree (Log2FC = −6.46). Aldehyde dehydrogenase X, mitochondrial (Log2FC = 3.06) was the most significantly increased protein and 3 beta-hydroxysteroid dehydrogenase type 5 (Log2FC = −4.59) was the most significantly decreased protein. Annexin A5 (Log2FC = 2.10) and perilipin-2 (Log2FC = 1.97) were increased while thyroid hormone-inducible hepatic protein (Log2FC = −2.67) was decreased.

In Ahr–/– mice, PCB126 exposure (vs. vehicle control) was not associated with any differentially abundant proteins (Fig. 4E). Primary proteomics results are summarized in Fig. 4F. A Venn diagram for the three comparisons relative to WT Vehicle is provided as Fig. 4G. Only four proteins were shared between all experimental groups. These included: CYP1A2, CYP2C50 isoform 2, fatty acid-binding protein 5, and acyl-coenzyme A thioesterase 1 (Table S1).

Secondary proteomics analyses

Secondary proteomics data analyses were performed in MetaCore. These analyses included: enrichment by protein function, pathway enrichment by gene ontology processes and interactions by protein function.

Enrichment by protein function

EPF analysis was performed to evaluate the protein classes most impacted by genotype or PCB126 exposure (Table 1). ‘Enzymes’ were the most enriched protein class (z-score range 6.42–17.84) while the ‘other’ protein class was enriched less than expected (z-score range −3.31 to −9.24) for the three main comparisons [e.g., WT (Vehicle vs. PCB126); Vehicle (WT vs. Ahr–/–) and PCB126 (WT vs. Ahr–/–]. EPF could not be performed for Ahr–/– (Vehicle vs. PCB126) because no proteins were differentially abundant for that comparison.

Table 1

Enrichment by protein function analysis (z-score).

Protein class	WT (Vehicle vs. PCB126)	Vehicle (WT vs. Ahr–/–)	PCB 126 (WT vs. Ahr–/–)
Ligands		1.87	2.60
Phosphatases		0.90	−0.05
Proteases		−0.67	−1.05
Kinases		−0.90	−0.40
Transcription factors
Receptors		−2.91	−2.57
Enzymes	6.42	17.84	15.52
Other	−3.31	−9.24	−8.07

For a given protein class, a positive z-score indicates that more proteins in that class were altered more than expected. Likewise, a negative z-score means that fewer proteins in the class were altered than expected.

Enrichment by gene ontology processes

Enrichment by GO processes was performed. For the WT (Vehicle vs. PCB126) comparison, 240 significant processes were identified (not shown). For Vehicle (WT vs. Ahr–/–) and PCB126 (WT vs. Ahr–/–), 1663 and 1947 significant processes were identified, respectively (not shown). Because there were no differentially abundant proteins detected for Ahr–/– (Vehicle vs. PCB126), there were no GO processes enriched by this comparison. The top 20 enriched GO processes by P-value for each comparison are given in Fig. 5. Due to overlap, a total of 35 pathways are provided in the figure. Because this analysis does not determine directionality, it is not surprising that the same top processes were affected by Ahr–/– or the AHR agonist, PCB126. Indeed, for the three comparisons with differentially abundant proteins, 31 of the top 35 pathways were enriched in all three. This consistency yields insight into the functionality of the aryl hydrocarbon receptor in liver. Moreover, these top GO processes were consistent with the EPF results which demonstrated enrichment in the enzyme protein class. Broadly, the top GO processes were involved in oxidation-reduction and other reactions impacting the metabolism of lipids, xenobiotics and organic acids; as well as the generation of energy. These results are not surprising, given the large number of differentially abundant proteins implicated in lipid metabolism that were associated with Ahr ablation and/or activation [e.g., fatty acid synthase, desaturases (n = 2), and binding proteins (n = 2); acyl-coenzyme A thioesterases (n = 5); and lipid droplet proteins like perilipin-2 and apolipoprotein A4; etc.] (Table S1). Carbohydrate metabolism was also impacted. For example, glycogen phosphorylase liver form, the rate limiting step in glycolysis, was down-regulated in the Ahr–/– (Table S1). This could help explain the decreased fasting blood glucose observed in these mice. Other top GO processes involved steroid metabolism, response to hormones and antibiotic biosynthesis. Ahr–/– was associated with significant alterations proteins impacting steroid synthesis, including: the upregulated 17-beta-hydroxysteroid dehydrogenases type 6 and type 13 (isoforms 1 and 2); and the down-regulated 3 beta-hydroxysteroid dehydrogenase type 5. The “dibenzo-p-dioxin metabolic process” was significantly associated with the WT (Vehicle vs. PCB126) and the PCB126 (WT vs. Ahr–/–) comparisons but not Vehicle (WT vs. Ahr–/–) (Table S1).

Figure 5

Pathway enrichment by Gene Ontology (GO) processes. Heatmap showing selected top enriched GO processes by P-value.

Enrichment by interactions by protein function

IPF analysis was performed, and for WT (Vehicle vs. PCB126), 140 significant interactions were identified (not shown). For Vehicle (WT vs. Ahr–/–) and PCB126 (WT vs. Ahr–/–), 279 and 335 significant interactions were identified, respectively (not shown). The top twenty significant over-connected objects by z-score are provided in Fig. 6 and Supporting Information Fig. S4 for each intergroup comparison. Although they did not rank in the top 20 IPFs, AHR and PXR were included in Fig. 6 for internal validation. There were no significant interactions by protein function for Ahr–/– (Vehicle vs. PCB126) because there were no differentially abundant proteins for this comparison. Fig. 6 and Fig. S4 includes a total of 54 unique objects. These include transcription factors (n = 13); receptors (n = 2); kinases (n = 1); proteases (n = 4); ligands (n = 1); phosphatases (n = 1); enzymes (n = 10) and others (n = 22). In general, the objects identified by intergroup comparisons involving the Ahr–/– mice (e.g., Vehicle (WT vs. Ahr–/–) and PCB126 (WT vs. Ahr–/–) were more similar to each other than to WT (Vehicle vs. PCB126).

Figure 6

Interactions by protein function. Heatmap of top interactions by protein function by z-score.

Seven objects were common to all three comparisons including: the transcription factors, AHR, zinc finger protein 125 (ZFP125), liver X receptor α, peroxisome proliferator-activated receptor α (PPARα) and hepatocyte nuclear factor α; the insulin receptor; and miR-132-5p. The shared transcription factors possibly contributed to the observed alterations in lipid and xenobiotic metabolism. While not the top object by z-score, it was nonetheless reassuring that AHR was detected by IPF analysis. The objects, PXR and the glucocorticoid receptor (GCR), were over-connected only for the comparisons involving Ahr–/– mice [e.g., Vehicle (WT vs. Ahr–/–) and PCB126 (WT vs. Ahr–/–)]. The PXR IPF data are consistent with the Cyp3a11 gene expression and protein abundance data presented in Fig. 1 and Table S1. Thus, the AHR and PXR IPF data were validated by internal controls. MiR-132-5p was the top object (by z-score) enriched by PCB126 in wild type mice (vs. vehicle control).

“Steroid metabolic processes” and “response to hormone” were among the top enriched GO processes (Fig. 5). Indeed, GCR and insulin receptor were among the top objects identified by IPF analyses involving the Ahr–/– mice. Thus, GCR expression and activity were determined by measuring the mRNA expression levels of Gcr and its target genes, tyrosine aminotransferase (Tat), serum/glucocorticoid regulated kinase (Sgk1) and nuclear factor kappa-inhibitor alpha (Nfkbia) by RT-PCR (Fig. 7). Gcr mRNA expression was upregulated by PCB126 exposure but not Ahr–/– genotype. However, GCR activity was increased in Ahr–/– mice, because the mRNA expression of Tat, Sgk1, and Nfkbia were increased; while PCB126 upregulated only Nfkbia (Fig. 7A–D). These data demonstrate varying degrees of GCR activation in Ahr–/– mice and with PCB126-induced ligand-activation of AHR.

Figure 7

Validation of additional proteomics findings. Hepatic mRNA levels were determined by RT-PCR for glucocorticoid receptor and its targets including: (A) Gcr, (B)Tat; (C) Sgk1; and (D) Nfkbia. Likewise, (E) perilipin-2 expression was determined. Expression of mRNA was normalized to the WT Vehicle group. Values are mean ± SD. A complete list of P-values (determined by two-way ANOVA with Tukey's post-test) is provided in the accompanying table. In the figure panels, statistical significance is denoted by: a = genotype effect; b = PCB effect; c = interaction effect; d = WT Vehicle vs. WT PCB126; e = Ahr–/– Vehicle vs. Ahr–/– PCB126; f = WT Vehicle vs. Ahr–/– Vehicle; g = WT PCB126 vs. Ahr–/– PCB126. AHR, aryl hydrocarbon receptor; NFKBIA, nuclear factor kappa-inhibitor alpha; PCB, polychlorinated biphenyl; SGK1, serum/glucocorticoid regulated kinase; TAT, tyrosine aminotransferase; WT, wild type.

PPARα and hypoxia inducible factor-1α (HIF1α) which heterodimerizes with the aryl hydrocarbon receptor nuclear translocator have been shown to regulate the expression of perilipin-2,. IPF showed over-connected interactions with PPARα and AHR, and perilipin-2 (PLIN2) protein was more abundant in Ahr–/– mice relative to WT controls. Hepatic Plin2 mRNA expression was measured by RT-PCR (Fig. 7E). Plin2 expression was altered according to genotype and was significantly higher in PCB126-exposed Ahr–/– mice relative to PCB126-exposed WT mice. There was a trend towards higher Plin2 expression in vehicle-administered Ahr–/– mice relative to vehicle-administered WT mice (P = 0.08, Fig. 7E). Thus, the perilipin-2 data are consistent at the mRNA and protein levels.

MicroRNAs are small non-coding RNAs that maintain cellular homeostasis and potentially modulate responses to environmental exposures. IPF determined enrichment in a total of nine miRs associated with PCB126 (n = 7), Ahr–/– (n = 4), or the combination of PCB126 and Ahr–/– (n = 3) (Supporting Information Table S2). Therefore, RT-PCR was performed to assess the hepatic expression of these nine miRs across the four groups (Supporting Information Fig. S5). The expression of miR-132-5p, miR-222-3p, and miR-544-3p could not reproducibly be detected across most samples. While differential miR expression was observed for the remaining miRs, the pattern was different than predicted by IPF. Ahr–/– was associated with significantly up-regulated miR-142-3p, miR-221-3p, miR-223-3p, and miR-150-5p and significantly down-regulated miR-192-3p. PCB126 exposure was not associated with these miRs. However, an interaction between exposure and genotype significantly decreased miR-122-5p expression in Ahr–/– but not WT mice (Fig. S5).

Discussion

PCB126 activated canonical AHR signaling in WT, but not Ahr–/–, as demonstrated by increased Cyp1a1 and Cyp1a2 mRNA expression and protein abundance only in WT mice. The WT mice results confirm our previously reported findings,. The 2-way ANOVA analysis, which included all four experimental groups, demonstrated several significant PCB126 effects. These included: decreased body weight gain (%); improved glucose tolerance (AUC); and increased hepatic triglycerides with altered hepatokine expression. However, on post-hoc intergroup comparison testing, PCB126 was not associated with any significant differences in the metabolic or liver phenotype of the WT mice (vs. vehicle control). However, eight hepatic proteins were differentially abundant in WT PCB126 vs. WT Vehicle groups. All eight of these were AHR-dependent, as PCB126 exposure was associated with zero differentially abundant proteins in the knockouts (Ahr–/– PCB126 vs. Ahr–/– Vehicle). In a previous publication, we reported that PCB126 exposures were associated with 396 differentially abundant hepatic proteins. The variable proteomic results across studies may potentially be explained by different: (i) statistical methods; (ii) strains of mice (Taconic vs. Jackson); (iii) diets (control vs. high fat); and (iv) exposure durations (2 vs. 12 weeks). Most notably, the present study utilized much more stringent criteria to identify significant differentially abundant proteins (e.g., FDR < 0.05 vs. FDR < 0.2; and −1 < log2FC < 1 vs. −0.5 < log2FC < 0.5). Filtering the prior study's results by the new more stringent FC criteria alone reduces the number of differentially abundant proteins to fewer than fifty. Importantly, several key consistencies were noted across studies. These include similar degrees of CYP1A1 and CYP1A2 protein up-regulation as well as the enrichment of similar metabolic pathways by PCB126. Regardless, the potential dose-responsiveness of the liver proteome to dioxin-like molecules, including PCB126, warrants future investigation.

The impact of whole-body Ahr ablation was much more striking than the effects of the ligand-activation of this receptor by PCB126. Ahr–/– mice exhibited severe metabolic disruption. Adiposity (%) and circulating leptin levels were increased, while glucose tolerance was decreased (despite decreased fasting glucose and insulin). Multiple blood and liver lipid species were increased accordingly. The observed hepatic steatosis was associated with liver injury and cell death (e.g., elevated AST, ALT and annexin A5) and markedly increased mRNA expression of the hepatokine, Fgf21. Despite increased steatosis, liver weight/body weight ratio was decreased in the knockouts. Several circulating pro-inflammatory cytokines, most notably IL-6, were paradoxically decreased. Not surprisingly, the Ahr–/– mice had 340 differentially abundant hepatic proteins (Ahr–/– Vehicle vs. WT Vehicle).

Unexpectedly, the AHR appeared to possibly be protective against hepatic steatosis (both spontaneous steatosis and PCB126-induced steatosis). Specifically, hepatic triglycerides were higher in PCB-exposed (vs. Vehicle–administered) Ahr–/– mice (but not in WT mice) with a trend towards a significant genotype-PCB126 interaction (P = 0.09). It was also somewhat surprising that PCB126 did not increase hepatic steatosis in WT mice vs. vehicle control. We recently reported on PCB126-induced NAFLD using a similar exposure protocol, albeit in Jackson instead of Taconic mice. While several other studies also reported PCB126-induced hepatic steatosis, these used higher doses and sometimes different routes of administration, sexes, strains, species or durations,,,.

Secondary MetaCore analyses of the primary proteomics data were performed to characterize the principal actions of Ahr genotype or AHR ligand-activation by PCB126 in liver. Given the phenotyping data, it was not surprising that EPF analysis determined enzymes to be the major protein class effected by either activation or ablation of the Ahr, with the strongest effects occurring in the knockouts. By GO process analysis, these enzymes were chiefly involved in metabolism of lipids/steroids, xenobiotics and organic acids as well as the generation of energy. Amongst these processes, the largest number were involved with lipid metabolism. Other enriched processes included antibiotic synthesis and response to hormone.

The IPF data suggested that some of the observed effects of aryl hydrocarbon receptor ligand-activation and/or deletion may have been mediated by over-connected interactions with transcription factors previously implicated in metabolism and NAFLD. These included nuclear receptors (e.g., GCR, liver X receptor α, PPARα, thyroid receptor β, hepatocyte nuclear factor α and PXR) and other transcription factors (ZFP125, sterol regulatory element-binding transcription factor 1, and cyclic-AMP-responsive-element-binding protein H, etc.)36, 37, 38. Indeed, altered expression of selected PXR, PPARα, or GCR targets was documented at the mRNA and/or protein levels in this study. The PPARα/aryl hydrocarbon receptor nuclear translocator target, perilipin-2, was up-regulated in Ahr–/– mice. Previous studies have shown that PCB126 decreased expression of PPARα and/or its target genes,,, and PLIN2 has previously been positively associated with NAFLD severity. Ahr–/– increased PLIN2 expression, possibly contributing to the observed NAFLD phenotype. This effect may have been indirectly mediated, occurring via over-connected interactions with the transcription factors regulating PLIN2 expression, including PPARα. While AHR-PPARα interactions are well-documented in constitutively active and Ahr null mice,, the putative interaction between the AHR and ZFP125 appears novel. ZFP125 is a FOXO1-inducible hepatic transcriptional repressor that caused NAFLD by decreasing hepatocyte lipid secretion to cause liver steatosis. It regulates several differentially abundant proteins which were associated with ligand-activation of AHR by PCB126 and/or Ahr–/– including fatty acid-binding protein 5 and apolipoprotein A4.

It was previously demonstrated that the GCR potentiates rat CYP1A1 induction by polycyclic aromatic hydrocarbons via direct interaction of the GCR with glucocorticoid response elements located in intron 1 of Cyp1a1. Here, we show that either ligand-activation or ablation of the aryl hydrocarbon receptor variably increased GCR target gene expression, suggesting that the AHR may also modulate GCR activity. This finding may have implications for human health. For example, some AHR ligands have demonstrated anti-inflammatory actions in fatty liver disease models,,. Perhaps these AHR-ligands could increase glucocorticoid receptor signaling to decrease liver inflammation. However, in the present study, the Ahr–/– genotype was associated with reduced circulating pro-inflammatory cytokines, while PCB126 exposure did not alter these mediators. More data are clearly required to better understand interactions between the AHR and the transcription factors elucidated by the enrichment by IPF analysis. However, cross-regulation of Ahr and Gcr mRNA expression has been reported.

To our knowledge, miRs have not previously been associated with Ahr ablation, but they have been associated with PCB exposures. Here IPF elucidated enriched interactions between nine miRs. The majority of these miRs were confirmed to be differentially expressed by RT-PCR, although not usually in the same pattern predicted by the IPF analysis. While the discrepancy between IPF and RT-PCR results could be explained if one set of results was incorrect, alternative explanations may exist (e.g., small miR-quenching RNAs). Interestingly three of the miRs with significantly different liver expression by RT-PCR (e.g., miR-122-5p, miR-221-3p, and miR-192-3p) were also associated with liver disease in an abstract evaluating circulating miRs in a human cohort with high residential PCB exposures. The initial miR findings from the present manuscript require confirmation and additional investigation. Nonetheless, the potential role of miRs in AHR-regulated liver homeostasis is an unanticipated, yet interesting, future research direction.

Major urinary proteins (MUPs) 1 and/or 17 were among the most highly down-regulated proteins identified in the Ahr–/– mice. As such, they warrant discussion. MUPs are unique members of the lipocalin superfamily that mediate both chemical and metabolic signaling. MUPs bind and stabilize pheromones to regulate their transport and release into the air from urine. MUPs also regulate carbohydrate and lipid metabolism. Recombinant MUP1 has anti-diabetic actions. It suppressed hepatic gluconeogenesis while stimulating mitochondrial biogenesis in skeletal muscle to increase energy expenditure. Little is known about the AHR and MUPs, but 2,3,7,8-tetrachlorodibenzo-p-dioxin repressed the signal transducer and activator of transcription 5b (STAT5b) target gene, Mup2, in an AHR-dependent manner. These data suggest the possibility that pheromones could be disrupted by environmental chemicals interacting with the AHR. Perhaps the effects of pheromones on chemical signaling and metabolism warrant additional investigation in the environmental health sciences (pheromone disruption).

While genetic knockout approaches are generally strong, the approach taken in the manuscript may be limited by the differential performance of the various strains of Ahr–/– mice on liver and metabolic endpoints in the published literature. While these differences could be related to genetic differences caused by the specific approaches used to ablate the Ahr, they could also be due to other factors including genetic background, diet, age, sex, microbiome, etc. To clarify, the Ahr–/– mice (C57Bl/6 background) utilized in this study were obtained from Taconic. This strain was derived by CXR Biosciences from a humanized AHR Mouse line through a Cre-mediated deletion of the human AHR sequence, and it is a whole-body knockout. Mouse exon 3 is deleted in the knockout, resulting in an out of frame splicing of exons 2 to exon 4. Strain specific-differences have been reported between the whole-body Ahr–/– generated by the Bradfield and Gonzalez groups (both of which are on similar genetic backgrounds–mixed C57Bl/6 and 129SV),. While these strains have significantly smaller livers, only the Bradfield strain was reported to have increased early-life hepatic steatosis which improved over time. In a 14 week feeding study, male Bradfield mice were protected from high fat diet (HFD)-induced obesity and displayed improved insulin signaling. Liver-specific knockouts generated using Bradfield's conditional Ahr–/– model (on a 129SV background) had increased HFD-induced hepatic steatosis (vs. WT controls),. Liver cell-specific knockouts have likewise been generated from this conditional model, demonstrating that Ahr ablation in stellate cells was associated with increased CCl4-induced fibrosis. However, these models been criticized due to the lower-affinity AHR present in 129SV mice, prompting the recent development of a series of new conditional knockouts using CRISPR-Cas9. Likewise, a novel tamoxifen-inducible and hepatocyte-specific Ahr–/– mouse model, had significantly reduced adiposity attributed to increased production of the thermogenic hepatokine, FGF21. Based on this literature, it is possible that our results might not be generalizable to other Ahr–/– strains. Future studies should address several related limitations of the present manuscript. For example, tissue-specific (e.g., liver, intestine, etc.) Ahr knockouts are required to dissect the direct effects of hepatic AHR from the systemic effects of extrahepatic AHR indirectly impacting the liver proteome. These studies could also include integrated multi-‘omics analyses (e.g., transcriptomics, metagenomics, phosphoproteomics, etc.) to better clarify AHR's role in liver health and disease.

Overall, the metabolic phenotyping and proteomic data converged on Ahr's critical role in the coordination and maintenance of systemic energy and lipid homeostasis. Perturbing this system influenced the development of obesity-associated diseases. For example, Ahr–/– mice had increased adiposity, hepatic steatosis and glucose intolerance with paradoxically decreased blood lipids. Mechanistically, Ahr mainly targeted liver metabolism (e.g., lipids, xenobiotics, organic acids, etc.), bioenergetics and endocrine function (e.g., production of hepatokines, steroids and pheromone binding proteins). Some of these effects may have been indirectly mediated by the interacting miRs or the transcription factors identified by MetaCore. While the proteome was impacted to a greater degree by Ahr–/– than by PCB126-induced AHR activation, the top GO processes identified were similar. Interestingly, the PCB126-associated liver proteome was entirely Ahr-dependent. Based on the literature, it is possible that some results could be dependent on the specific strain of Ahr–/– mouse utilized. Regardless, the Ahr and its ligands warrant more research in metabolic health and disease.