Multi-omics insights into functional alterations of the liver in insulin-deficient diabetes mellitus.

OBJECTIVE: The liver regulates the availability of insulin to other tissues and is the first line insulin response organ physiologically exposed to higher insulin concentrations than the periphery. Basal insulin during fasting inhibits hepatic gluconeogenesis and glycogenolysis, whereas postprandial insulin peaks stimulate glycogen synthesis. The molecular consequences of chronic insulin deficiency for the liver have not been studied systematically.
METHODS: We analyzed liver samples of a genetically diabetic pig model (MIDY) and of wild-type (WT) littermate controls by RNA sequencing, proteomics, and targeted metabolomics/lipidomics.
RESULTS: Cross-omics analyses revealed increased activities in amino acid metabolism, oxidation of fatty acids, ketogenesis, and gluconeogenesis in the MIDY samples. In particular, the concentrations of the ketogenic enzyme 3-hydroxy-3-methylglutaryl-CoA synthase 2 (HMGCS2) and of retinol dehydrogenase 16 (RDH16), which catalyzes the first step in retinoic acid biogenesis, were highly increased. Accordingly, elevated levels of retinoic acid, which stimulates the expression of the gluconeogenic enzyme phosphoenolpyruvate carboxykinase (PCK1), were measured in the MIDY samples. In contrast, pathways related to extracellular matrix and inflammation/pathogen defense response were less active than in the WT samples.
CONCLUSIONS: The first multi-omics study of a clinically relevant diabetic large animal model revealed molecular signatures and key drivers of functional alterations of the liver in insulin-deficient diabetes mellitus. The multi-omics data set provides a valuable resource for comparative analyses with other experimental or clinical data sets.

Introduction

The liver is the central glucoregulatory organ since all insulin secreted by the beta cells enters the liver while only part of it reaches the peripheral circulation to ensure appropriate glucose uptake by the main insulin target tissues (muscle and adipose tissue) and to maintain physiological blood glucose levels (reviewed in [1], [2]). The liver is thus exposed to two- to four-fold higher levels of insulin than peripheral insulin target tissues (reviewed in [3]). During fasting, basal insulin inhibits gluconeogenesis and glycogenolysis, whereas postprandial insulin peaks stimulate glucose storage by the liver as glycogen (reviewed in [1]).

To determine consequences of missing insulin action in the liver, mice with a liver-specific insulin receptor gene (Insr) knockout (LIRKO) were generated [4]. Since LIRKO mice develop progressive hepatic dysfunction (reviewed in [5]), they are not suitable as a model for studying long-term effects of insulin deficiency in the liver.

We therefore used a physiologically more relevant pig model of insulin-deficient diabetes mellitus. Transgenic pigs expressing mutant insulin C94Y, a model for mutant INS gene-induced diabetes of youth (MIDY), reveal impaired insulin secretion, endoplasmic reticulum stress, and apoptosis of the beta cells [6]. MIDY pigs show diabetic complications, such as cataract development [6], reduced capillarization and pericyte investment in the myocardium [7], and diabetes-associated retinal changes [8]. MIDY pigs were maintained for two years with limited insulin treatment to represent poorly controlled diabetes in humans. Fasting plasma glucose and fructosamine concentrations of MIDY pigs were permanently elevated, C-peptide levels decreased with age and were undetectable at 2 years. Plasma glucagon and beta hydroxybutyrate levels were chronically elevated [9]. A comprehensive biobank was established from 4 female MIDY pigs and 5 female wild-type (WT) littermates [9]. To systematically address hepatic changes in response to chronic insulin deficiency and hyperglycemia, an integrative multi-omics analysis [10] covering transcripts, proteins and different metabolite/lipid classes was performed. The design of the study is shown in Figure 1.

Figure 1

Outline of the multi-omics study of consequences of insulin-deficient diabetes for the liver. MIDY pigs (n = 4) and WT littermate controls (n = 5) were maintained for two years under standardized conditions. MIDY pigs had significantly (p < 0.001) elevated fasting blood glucose (FBG) levels (310 ± 39 mg/dL vs. 120 ± 26 mg/dL in WT) and plasma beta hydroxybutyrate (BHB) concentrations (48 ± 20 μmol/L vs. 11 ± 6 μmol/L in WT). C-peptide was undetectable in plasma from MIDY pigs. A complex biobank was established, including liver samples taken by systematic random sampling [9].

Material and methods

Samples

This study used liver tissue samples of two-year-old female MIDY pigs (n = 4) and female WT littermates (n = 5) harvested by systematic random sampling [11] for different omics analyses. The pigs were fasted overnight before necropsy. The samples were shock-frozen on dry ice and stored at −80 °C in the Munich MIDY Pig Biobank [9] until analysis. All samples were processed in parallel to avoid variation related to different storage times and batch effects.

Transcriptomics

Liver samples were homogenized in Trizol, and total RNA was isolated with chloroform following manufacturer's protocol. Isolated total RNA was quantified (Nanodrop, ND1000) and quality controlled (Agilent, Bioanalyzer 2100). Good quality RNA (RIN >7.0) was used to construct sequencing libraries (Nugen, Encore Complete RNA-Seq library system). This kit enables the analysis of transcriptome profiles with reduced representation of ribosomal RNA because of not so random priming during cDNA synthesis. All libraries were sequenced on a HiSeq 1500 (Illumina) as 100 b single reads. Demultiplexing and quality control were performed on the obtained FastQ files followed by mapping to the S.scrofa 11.1 reference genome using the gapped-mapper STAR. HTSeq [12] using strict intersection mode and a minimum alignment quality of 10 was used to quantify the number of hits to each gene. DESeq2 [13] with outlier replacement and independent filtering was used to detect differentially abundant transcripts between MIDY and WT samples. Pre-ranked unweighted gene set enrichment analyses (GSEA) [14] were performed on the signed log transformed p-values as ranking metric using MSigDB [15] and the specific Sus scrofa KEGG pathways [16]. For network visualization, Cytoscape [17] was used with the ClueGO [18] and CluePedia apps to analyze significant genes.

Proteomics

For quantitative proteome analysis, liver samples were taken from the same localizations as for mRNA measurements and homogenized as previously described [19]. Protein concentrations were determined using the Pierce 660 nm Protein Assay (Thermo Scientific) [20]. 100 μg of protein were digested with Lys-C (Wako) for 4 h and trypsin (Promega) overnight at 37 °C [9]. For nano-LC-MS/MS analysis, a Q Exactive HF-X mass spectrometer equipped with an UltiMate 3000 nano LC system (Thermo Scientific) was used. Briefly, 2.5 μg of peptides were separated at 200 nL/min using consecutive linear gradients from 1% to 5% solvent B (0.1% formic acid in acetonitrile) in 10 min, from 5% to 25% B in 115 min and from 25% to 50% B in 20 min. Spectra were acquired using one survey scan at a resolution of 120,000 from 380 to 2000 m/z followed by MS/MS scans of the 24 most intense peaks at a resolution of 15,000. For protein identification (FDR < 1%) and label-free quantification, MaxQuant (v. 1.6.1.0) [21] and the NCBI RefSeq Sus scrofa database (v. 3-13-2018) was used. Identifications were filtered for at least three valid values in one group and missing values were replaced from normal distribution using the data imputation feature implemented in Perseus [22]. Functional annotation enrichment analyses were performed using STRING [23] and Proteomaps [24].

Western blot analysis of insulin receptor signaling

Concentrations and phosphorylation levels of insulin receptor (INSR)-related signaling molecules in the liver were evaluated by western blot analyses as described previously [25], [26]. Briefly, liver tissue samples were homogenized in Laemmli extraction buffer, and the protein content was determined by the bicinchoninic acid protein assay. Twenty micrograms of total protein were separated by SDS-PAGE and transferred to PDVF membranes (Millipore) by electro-blotting. Membranes were washed in TBS with 0.1% Tween-20 and blocked in 5% w/v fat-free milk powder (Roth) for 1 h. The membranes were then washed again and incubated in 5% w/v BSA (Roth) solution with specific primary antibodies overnight at 4 °C. After washing, the membranes were incubated in 5% w/v fat-free milk powder solution with the appropriate secondary antibodies for 1 h. The antibodies and concentrations used are listed in Supplementary Table 1. Bound antibodies were detected using the ECL Advance Western Blotting Detection Kit (GE Healthcare). Band intensities were quantified using the ImageQuant software package (GE Healthcare).

Targeted metabolomics

The targeted metabolomics approach was based on liquid chromatography-electrospray ionization-tandem mass spectrometry (LC-ESI-MS/MS) and flow injection analysis-electrospray ionization tandem mass spectrometry (FIA-ESI-MS/MS) measurements using the AbsoluteIDQTM p180 Kit (Biocrates Life Sciences AG). Liver tissue samples were processed, extracted, and quantified as described in full detail previously [27], [28]. To each mg of frozen wet liver tissue, 3 μL of a dry ice cooled mixture of ethanol/phosphate buffer (85/15 v/v) were added. Out of 10 μL liver tissue homogenate, 188 metabolites were quantified. Sample handling was performed by a Hamilton Microlab STAR™ robot (Hamilton Bonaduz AG) and an Ultravap nitrogen evaporator (Porvair Sciences), beside standard laboratory equipment. Mass spectrometric analyses were done on an API 4000 triple quadrupole system (Sciex Deutschland GmbH) equipped with a 1200 Series HPLC (Agilent Technologies Deutschland GmbH) and a HTC PAL auto sampler (CTC Analytics) controlled by the software Analyst 1.6.2. Data evaluation for quantification of metabolite concentrations and quality assessment was performed with the software MultiQuant 3.0.1 (Sciex) and the MetIDQ™ software package. Internal standards were used as reference for the calculation of metabolite concentrations. The concentrations in the tissue samples are given in pmol/mg wet tissue and the concentrations in tissue homogenate in μM. Limit of detection (LOD) for each metabolite in tissue homogenate was calculated by multiplying the median concentration of the three zero-samples (ethanol/phosphate buffer) times 3. Metabolite concentrations were log transformed and Pareto scaled to model them for two-tailed Student's t-test.

Lipidomics

Liver samples were homogenized on ice in ammonium-bicarbonate buffer (150 mM ammonium bicarbonate, pH 7) with an ultra-turrax homogenizer. Protein content was assessed using BCA Protein Assay Kit (Thermo Fisher). Equivalents of 20 μg of protein were taken for mass spectrometry analysis. Mass spectrometry-based lipid analysis was performed by Lipotype GmbH as described [29]. Briefly, lipids were extracted using a two-step chloroform/methanol procedure [30] and spiked with an internal lipid standard mixture. After extraction, the organic phase was transferred to an infusion plate and dried in a speed vacuum concentrator. 1st step dry extract was re-suspended in 7.5 mM ammonium acetate (Sigma) in chloroform/methanol/propanol (1:2:4, v:v:v) and 2nd step dry extract in 33% ethanol solution of methylamine in chloroform/methanol (0.003:5:1; v:v:v). All liquid handling steps were performed using a Hamilton Robotics STARlet robotic platform with the Anti Droplet Control feature for organic solvent pipetting. Samples were analyzed by direct infusion on a QExactive mass spectrometer (Thermo Scientific) equipped with a TriVersa NanoMate ion source (Advion Biosciences). Samples were analyzed in both positive and negative ion modes with a resolution of Rm/z=200 = 280,000 for MS and Rm/z=200 = 17,500 for MS/MS experiments, in a single acquisition. MS/MS was triggered by an inclusion list encompassing corresponding MS mass ranges scanned in 1 Da increments. Data were analyzed with a lipid identification software based on LipidXplorer [31]. Only lipid identifications with a signal-to-noise ratio >5, and a signal intensity 5-fold higher than in corresponding blank samples were considered for further data analysis. Individual lipid measurements were combined into family groups and normalized for further analysis.

Quantification of retinol, retinal and retinoic acid

Prior to lipid extraction 100 mg of liver tissue was homogenized in 200 μL citric acid (0.4 mol/L) using a pestle. Lipids were extracted following Folch's protocol. Briefly, 300 μL methanol and subsequently 2 × 600 μL chloroform were added to the liver homogenate. The mixture was shaken vigorously for 2 × 10 min and centrifuged. The lower organic layer was transferred into a brown glass vial and evaporated to dryness under nitrogen at 30 °C. The residue was resolved in 1 mL ethanol and analyzed using HPLC-MS/MS.

The analysis was performed using an Agilent 1290 HPLC coupled with a triple quadrupole mass spectrometer Agilent 6470 equipped with an electrospray jet stream ion source. A Zorbax SB-C18 (50 × 2.1 mm, 1.8 μm) was used as stationary phase. A gradient of ammonium formate (5 mM)/formic acid (0.05%) in water and acetonitrile (30–95% acetonitrile in 10 min) was the mobile phase. The mass spectrometer was operated in positive/negative switching mode because retinoic acid reveals a more selective signal in the negative mode. The list of transitions for each compound is shown in Supplementary Table 2.

Glutathione assay

Liver glutathione (GSH) and oxidized glutathione (GSSG) concentrations were determined using a commercial Glutathione Colorimetric Detection Kit (EIAGSHC, Invitrogen) following the manufacturer's instructions. Briefly, 40–70 mg of frozen liver tissue was homogenized with a rod homogenizer (Polytron® PT 2500 E) in ice-cold 1 x PBS solution and immediately centrifuged (14,000 rpm, 10 min, 4 °C). An aliquot of the supernatant was removed for protein quantification using the Pierce™ 660 nm Protein Assay (Thermo Scientific). The remaining supernatant was deproteinized with 5% 5-sulfo-salicylic acid dihydrate solution (SSA). For GSSG determination samples were pre-treated with 2-vinylpyridine (2VP) and incubated for one hour at room temperature. Colorimetric reaction was detected at a wavelength of 405 nm using a Tecan M200 Infinite Pro plate reader. Free glutathione concentration was calculated by subtracting GSSG from GSH. Glutathione concentration of the homogenized liver tissue was normalized for tissue protein content and expressed as μmol/g protein.

Quantification of IBA1-positive macrophages (Kupffer cells) in liver samples

Hepatic macrophages were detected by immunohistochemistry in sections of three systematically randomly sampled, formalin-fixed paraffin-embedded (FFPE) liver tissue samples per case (WT: n = 5; MIDY: n = 4), using a goat polyclonal anti-IBA1 (ionized calcium binding adaptor molecule 1) antibody (ab5076, abcam) and a biotinylated rabbit anti-goat Ig secondary antibody (BA-5000, Vector). Diaminobenzidine was used as chromogen and hemalum as nuclear counterstain. The volume density of IBA1-positive macrophages in the liver (VV(macrophages/liver)) was determined following the principle of Delesse and calculated as the sum of cross-sectional areas of IBA1-positive cells, divided by the sum of cross-sectional areas of liver tissue in 31 ± 2 systematically randomly sampled section areas per case. Area densities were determined by differential point counting [11]. In each case, >12,000 points were counted.

Quantification of non-esterified fatty acids (NEFA) in plasma

Frozen EDTA-plasma aliquots were thawed in a fridge at 2–4 °C for 1–2 h. Subsequently, samples were shortly vortexed, centrifuged (5,000× g, 10 min, 8 °C, Biofuge Fresco, Heraeus) and analyzed within one hour. NEFA measurements were performed using an AU480 clinical chemistry analyzer (Beckman–Coulter) with the NEFA HR reagent kit (Wako Diagnostics) with corresponding calibrator and controls.

Statistical analysis

Statistics and visualizations were performed in Perseus [22] and R [32] using the gplots [33] ggplot2 [34] packages. Differences between MIDY and WT were evaluated using two-tailed Student's t-tests, where appropriate. DESeq2 was used to detect differentially abundant transcripts. The Benjamini-Hochberg procedure was used for FDR calculation in the case of transcriptomics, metabolomics and lipidomics. A permutation-based FDR estimation was used for statistical evaluation of differentially abundant proteins. Values were considered significant at FDR < 0.05.

For the comparative analysis of proteomics and transcriptomics data, datasets comprising protein abundance ratios and DESeq2 normalized mRNA abundance ratios, respectively, were combined and condensed on common identifications. A scatter plot of matched abundance ratios was color-coded according to the significance of regulation on the transcriptome (DESeq2 corrected p-values) and on the proteome dimension (permutation-based FDR corrected p-values). 2D annotation enrichment was performed on the merged quantitative 2D data using the algorithms implemented in Perseus [35]. The significance cutoff restricting the correlating, non-correlating, and anti-correlating regions was calculated by a nonparametric two-sample test. For significant functional categories at p < 0.01, abundance ratios of the corresponding proteins were separately replaced by ranks in both transcriptomics and proteomics dimension and the average rank per category was rescaled to a 2D score (MIDY/WT transcriptome and proteome score) between −1 and 1.

Results

Overview of transcriptome differences

Comprehensive mRNA profiles of the MIDY (n = 4) and WT (n = 5) liver samples were generated by random-primed cDNA sequencing (RNA-Seq) with reduced representation of ribosomal RNA. The average depth of mapped reads was ∼32 million reads. In total, transcripts of 14,818 different genes were identified, 320 with significantly (DESeq2; adj. p-value < 0.05) higher (Supplementary Table 3A) and 213 with lower (Supplementary Table 3B) abundance in MIDY vs. WT pigs. Figure 2A shows an MA-plot of the differentially expressed genes in MIDY pig liver, Figure 2B a heatmap of the 20 most upregulated and the 20 most downregulated transcripts. Differences are expressed as log2 fold change (l2fc).

Figure 2

Transcriptome differences between liver tissue samples from MIDY and WT pigs. A: Plot showing the log fold change between MIDY/WT and the mean count abundance per gene. The red and blue colored dots indicate transcripts with significantly increased or decreased abundance (FDR < 0.05). B: Heat map showing the 20 transcripts with the most significant (FDR < 0.05) increase and the 20 transcripts with the most marked decrease in abundance in MIDY vs. WT pigs. Red indicates higher expression in MIDY compared to WT and blue the reverse. Genes with the “LOC” name could not be matched to a recognized gene name. C: A network showing the ClueGO functional enrichment analysis for differentially abundant transcripts (red squares: increased in MIDY; blue diamonds: decreased in MIDY). The circles represent enriched KEGG pathways. The size of the circles indicates the significance of enrichment, the color code from red to blue indicates the proportion of transcripts with increased (red) and decreased (blue) abundance in MIDY pigs.

GSEA using the KEGG database identified gene sets related to amino acid metabolism, gluconeogenesis/glycolysis, glucagon signaling, retinol metabolism, peroxisome proliferator activated receptor (PPAR) signaling, and peroxisome enriched in the MIDY samples, whereas gene sets associated with immune functions and extracellular matrix interactions were enriched in the WT samples (Supplementary Table 4A,B). Related to amino acid metabolism, the levels of transcripts for enzymes involved in the degradation of specific amino acids were significantly increased: glutamic-pyruvic transaminase 2 (GPT2; alanine), glutamic-oxaloacetic transaminase 1 (GOT1; aspartate), glutaminase (GLS2; glutamine), arylformamidase (AFMID; tryptophan), homogentisate 1,2-dioxygenase (HGD; tyrosine, phenylalanine), serine dehydratase (SDS; serine), histidine ammonia-lyase (HAL; histidine), aminoadipate-semialdehyde synthase (AASS; lysine), aldehyde dehydrogenase 7 family member A1 (ALDH7A1; lysine), and kynurenine aminotransferase 1 (KYAT1; tryptophan, cysteine conjugates).

A ClueGO functional annotation network analysis using the significant genes (adj. p-value < 0.05) revealed similar pathways as GSEA (Figure 2C; Supplementary Table 5).

Overview of proteome differences

Quantitative LC-MS/MS-based proteomics identified a total of 2,535 proteins with high confidence (FDR < 0.01, [36]) (Supplementary Table 6). MIDY and WT samples were clearly separated by hierarchical clustering (Figure 3A) and PCA (Figure 3B).

Figure 3

Quantitative proteome analysis of liver tissue from MIDY and WT pigs. A: Unsupervised hierarchical clustering of normalized LFQ intensity values. Liver proteomes of MIDY animals are separated from WT replicates. The color code indicates z-score normalized expression values. B: Principal component analysis (PCA) clearly separates proteomes from MIDY and WT pigs. Spots represent individual animals. C: Volcano plot of log2 fold changes (MIDY/WT). Red and blue dots indicate differentially abundant proteins. Black curves represent the permutation-based FDR significance cutoff. Prominent differentially abundant proteins are highlighted. D and E: Proteomaps illustrating functional changes in the MIDY pig liver proteome. Treemaps for proteins significantly increased (D), and decreased (E), in abundance in MIDY compared to WT pigs (FDR < 0.05) are shown. Functionally related proteins are annotated based on KEGG orthology. Related functional categories are arranged in adjoining locations and share similar colors. Polygon areas represent the mass fraction of the corresponding proteins, i.e. protein abundances weighted by protein size.

A two-sided t-test with a permutation-based FDR approach revealed 60 significantly (FDR < 0.05) more abundant proteins and 84 less abundant proteins in MIDY vs. WT samples (Figure 3C; Supplementary Table 7). Among the proteins with the highest abundance increase in MIDY liver were retinol dehydrogenase 16 (RDH16; l2fc 4.7, p = 0.0226) and 3-hydroxy-3-methylglutaryl-CoA synthase 2 (HMGCS2; l2fc 2.7, p = 0.0064) (Figure 3C). In addition, the abundances of phosphoenolpyruvate carboxykinase (PCK1; l2fc 1.0, p = 0.0015) and of several other enzymes involved in gluconeogenesis, i.e., glutamic-pyruvic transaminase 2 (GPT2; l2fc 1.2, p = 0.0316), l-lactate dehydrogenase B chain (LDHB; l2fc 1.1, p = 0.0063), and alanine-glyoxylate aminotransferase (AGXT; l2fc 1.2, p = 0.0022) were significantly increased in the MIDY samples. The set of proteins most decreased in abundance contained, among others, collagenous members of the extracellular matrix, e.g. collagen alpha-1(I) chain (COL1A1; l2fc −1.7, p = 0.0259) and collagen alpha-1(XIV) chain (COL14A1; l2fc −1.4, p = 0.0169) (Figure 3C).

STRING analysis targeting GO annotations and KEGG pathways was performed for the differentially abundant proteins (Supplementary Table 8). Results of the Proteomaps analysis are shown in Figure 3D,E. Proteins more abundant in MIDY samples are involved in amino acid metabolism, gluconeogenesis/glycolysis, and tricarboxylic acid (TCA) cycle; the less abundant proteins have functions in pathogen defense response, response to cellular stress, or in cell signaling and genetic information processing.

Cross-omics comparisons

An integrated analysis of transcriptome and proteome changes was performed to investigate extent and levels of transcriptional or post-transcriptional regulation. A total of 1,572 transcripts/proteins could be matched as intersection between both data sets and correlations were calculated for products of individual genes (Figure 4A) and on the basis of functional categories (Figure 4B).

Figure 4

Correlation between proteomics and transcriptomics data in MIDY vs. WT liver tissue. A: Scatter plot of protein abundance ratios against corresponding mRNA ratios. Proteins with significant abundance alterations (FDR < 0.05) but no significant change in the corresponding transcript level are marked in blue. Significant regulation (FDR < 0.05) only in the transcriptome but not in the proteome is indicated in red. Common regulation at both levels is depicted in purple. Selected hits are highlighted. R: Pearson R. B: 2D annotation enrichment analysis between proteome and transcriptome data. Abundance ratios between MIDY and WT pigs were rescaled and depicted as transcriptome and proteome score. Significant pathways and gene ontology categories with p < 0.01 are shown. Terms located close to the ascending diagonal indicate common regulation at the transcriptome and proteome level. Annotation categories are color coded. Selected functional groups are highlighted.

Overall, the correlation between mRNA and protein log2 fold changes was moderate (R = 0.27), in line with the general observation that transcript levels are not sufficient to predict protein levels (reviewed in [37]). However, changes in transcript and protein levels were strikingly concordant for the most significantly affected gene products (e.g. HMGCS2, RDH16, SLC22A7, COL1A1; Figure 4A).

In addition, transcriptome and proteome data were evaluated in a 2D annotation enrichment analysis [35]. Functional processes and pathways, such as “urea cycle”, “arginine biosynthesis,” “gluconeogenesis,” “glucagon signaling pathway,” and “biosynthesis of amino acids” were enriched in MIDY liver tissue (p < 0.01), while functional categories related to extracellular matrix organization and defense response were enriched in the WT samples, both at the transcriptome and proteome level (Figure 4B, Supplementary Table 9).

Insulin receptor activation and downstream signaling

While INSR transcript and protein levels were significantly increased, INSR phosphorylation was significantly reduced in MIDY samples (Figure 5). Phosphoinositide 3-kinase (PI3K) was as a trend reduced in abundance, but phosphorylated PI3K levels were not different between the two groups. The phosphorylation levels of 3-phosphoinositide-dependent protein kinase-1 (PDPK1), protein kinase B (PKB, AKT), and glycogen synthase 3 beta (GSK3B) were significantly reduced in MIDY liver samples. Phosphorylated forkhead box protein O1 (pFOXO1) was markedly reduced, but also total FOXO1 levels were lower in MIDY than in WT samples. Total concentrations and phosphorylation levels of mechanistic target of rapamycin (mTOR), AMP-activated protein kinase (AMPK), and ribosomal protein S6 were not different between MIDY and WT samples (Supplementary Fig. 1).

Figure 5

Western blot analysis of insulin receptor activation and downstream signaling molecules. Densitometric data were square root transformed to approximate normal distribution and evaluated using Student's t-tests. The bar diagrams show means and standard deviations. Significant differences between MIDY and WT pigs are indicated by asterisks: *p < 0.05; **p < 0.01; ***p < 0.001; °borderline significance (p < 0.07).

Overview of metabolome and lipidome differences

The results of the targeted metabolomics analysis are shown in Supplementary Table 10. The concentrations of lysine and methionine (∼170% of WT) and of the branched chain amino acids (BCAA) leucine, isoleucine, and valine (∼125% of WT) were increased in MIDY samples. In contrast, the concentration of serine was reduced (∼60% of WT) (Figure 6A). Arginine was not detected in MIDY samples, but detectable at low concentrations in 4/5 WT samples (0.45 ± 0.20 pmol/mg). Among the biogenic amines, the most prominent changes were decreased creatinine and serotonin concentrations and increased spermine and histamine levels in MIDY samples. Other compounds such as kynurenine, methionine-sulfoxide, and dimethylated arginine (DMA) were increased in abundance but with large variances (Figure 6B). No significant differences in liver glutathione (GSH), oxidized glutathione (GSSG), and free GSH concentrations were observed between MIDY pigs and WT controls (Supplementary Fig. 2).

Figure 6

Relative metabolite abundance changes between MIDY and WT liver tissue. The graphs show the relative abundances of amino acids (A), biogenic amines (B), and selected metabolic indicators (C), determined by targeted metabolomics. Graph (D) shows the relative abundance of lipid groups determined by lipidomics. The dashed line indicates the mean abundance of the WT samples, bars show the mean relative abundance in MIDY samples. Error bars represent the standard error of the mean (SEM). Abbreviations: Ac-Orn, N-acetylornithine; ADMA, asymmetric dimethylarginine; AAA, aromatic amino acid; Met-SO, methionine sulfoxide; PEA, phenethylamine; t4-OH-Pro, hydroxyproline; DMA, dimethylarginine; H1, hexose; AC, acyl carnitine; C(n), acyl carnitine chain length; DC, dicarboxylated, OH, hydroxylated; PC, glycerophosphocholines; SFA, saturated fatty acid; MUFA, mono-unsaturated FA; PUFA, polyunsaturated FA; SM, sphingomyelins; CE, cholesteryl ester; CL, cardiolipin; Cer, ceramide; Chol, cholesterol; DAG, diacylglycerol; HexCer, hexosylceramide; LP(C,E,S), lysophosphatidyl (choline, ethanolamine, serine); PA, phosphatidate; PE, phosphatidylethanolamine; O-, ether linked; P(C,E,G,I,S), phosphatidyl (choline, ethanolamine, glycerol, inositol, serine); TAG, triacylglycerol.

The concentrations of long-chain acylcarnitines (C16, C18) and the ratio of (C16 + C18) to free carnitine (C0) were significantly increased (Figure 6C). In contrast, the levels of short-chain acylcarnitines (C2, C3, C4, C5) and the ratios of acetylcarnitine (C2) to C0 and of short-chain acylcarnitines (C2 + C3) to C0 were significantly decreased in MIDY samples. In addition, the ratio of total acylcarnitines to C0 was significantly decreased, while the ratios of dicarboxy-acylcarnitines to total acylcarnitines and of hydroxy-acylcarnitines to total acylcarnitines were increased in MIDY vs. WT samples. Total sphingomyelin (SM) and hydroxy-sphingomyelin (SM-OH) levels as well as the ratio of SM to phosphatidylcholines (PC) were significantly decreased in MIDY samples (Figure 6C).

The lipidomics analysis confirmed the metabolomics data in that PC were unchanged, while SM were reduced in MIDY samples. Concentrations of cholesterol (Chol) and phosphatidylserine (PS) were slightly reduced and lyso-phosphatidylserine (LPS) and phosphatidic acid (PA) more markedly reduced. In contrast, diacylglyceride (DAG) and triacylglyceride (TAG) levels were increased in MIDY samples (Figure 6D). The complete lipidomics data set is provided in Supplementary Table 11.

Since studies in rats demonstrated that the rates of fatty acid esterification into hepatic triglyceride were dependent on the concentration of free fatty acids in plasma, but independent of plasma insulin concentrations and hepatocellular insulin signaling [38], we measured the concentrations of non-esterified fatty acids (NEFA) in plasma samples from MIDY and WT pigs. NEFA concentrations were as a trend, but not significantly, elevated in the MIDY samples (Supplementary Fig. 3).

To clarify if the markedly increased concentration of RDH16 in the MIDY samples affects the levels of retinoids (Figure 7A,B), we quantified retinol, retinal and retinoic acid by mass spectrometry. While retinol levels were not significantly different between the two groups, the concentrations of retinal and retinoic acid were significantly increased in MIDY vs. WT samples (Figure 7C).

Figure 7

Activity of RDH16 in the retinoid metabolism in MIDY liver. A: Schematic representation of the two-step reaction of retinol to retinoic acid. The first and rate-limiting step is catalyzed by retinol dehydrogenases. B: Increased abundance of retinol dehydrogenase 16 (RDH16) in MIDY compared to WT pigs. C: Quantification of retinol, retinal, and retinoic acid in extracts from MIDY and WT liver tissue. Concentrations are given in ng per g tissue. Differences were tested for statistical significance using Student's t-test.

Discussion

To systematically assess consequences of insulin deficiency for the liver, we analyzed liver samples from a genetically engineered pig model for mutant INS gene-induced diabetes of youth (MIDY) and WT littermate controls. In contrast to diabetes induction by pancreatectomy or treatment with streptozotocin, the primary cause of insulin deficiency in MIDY pigs, i.e., formation of misfolded insulin resulting in impaired insulin secretion and beta-cell apoptosis, is limited to the beta cells, thus excluding confounding effects by an invasive surgery or toxicity to other cell types (reviewed in [39]). Circulating glucagon levels were consistently elevated in MIDY pigs [9], in line with the progressive loss of beta cells and the lacking paracrine control of glucagon secretion from neighboring alpha cells by insulin (reviewed in [40]). All animals were maintained under standardized conditions, the MIDY pigs with limited subcutaneous insulin treatment. Studies in rats [41] and dogs (reviewed in [42]) demonstrated that subcutaneous insulin at therapeutic doses distributes to muscle and adipose tissue, but barely to the liver. This was supported by significantly reduced phosphorylation levels of INSR and downstream signaling molecules PDPK1, AKT and GSK3B in MIDY liver samples. The limited insulin treatment of MIDY pigs was thus no confounding factor in our study of hepatic consequences of insulin-deficient diabetes mellitus. A standardized biobank of two-year-old WT and MIDY pigs was established applying the principles of systematic random sampling [9]. Variation in sample quality due to sample collection and storage could thus be minimized.

A multi-omics analysis combining the high coverage of transcriptome profiling, the immediate functional relevance of the protein layer and quantitative readouts of relevant metabolite classes was performed, to reveal biological processes and pathways altered by insulin-deficient diabetes in the liver and to identify molecular key drivers of these alterations.

Increased abundance of gluconeogenic enzymes suggests stimulated hepatic gluconeogenesis

The abundances of PCK1, the rate limiting enzyme of gluconeogenesis (reviewed in [43]), and of several other enzymes involved in gluconeogenesis were significantly increased in the MIDY samples. In contrast, the transcript level of PFKFB3 encoding 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3, a key stimulator of glycolysis (reviewed in [44]), was significantly decreased (l2fc −1.2, p = 3.1e-6).

The expression of PCK1 is stimulated by retinoic acid, which is generated in a two-step reaction from retinol (Figure 7A). The first step is catalyzed by retinol dehydrogenases. In human hepatoma cells, serum starvation stimulated the expression of retinol dehydrogenase genes RDH10 and RDH16 while insulin in serum-free medium decreased the expression of both genes. The latter effect was dependent on PI3K and AKT, leading to phosphorylation and degradation of forkhead box O1 (FOXO1) [45], an essential transcription factor for the expression of RDH genes [43]. The marked abundance increase of RDH16 in MIDY samples emphasizes the role of insulin as a negative regulator of RDH16 expression. Significantly reduced phosphorylation levels of INSR, AKT, and FOXO1 in MIDY samples point to a similar regulation of RDH16 by insulin in vivo as previously observed in vitro [43]. Significantly increased concentrations of retinal and retinoic acid in MIDY samples (Figure 7C) demonstrate the biological relevance of increased RDH16 levels, since the abundance of RDH10 mRNA was not altered in MIDY vs. WT samples. In addition, the level of CYP26A1 mRNA encoding cytochrome 26 A1, the main retinoic acid hydroxylase [46], was significantly decreased in MIDY samples (l2fc −1.0, p = 0.0011).

Glycogen synthase 2 (GYS2), which catalyzes the rate-limiting step of hepatic glycogen synthesis [47], was more abundant in MIDY than in WT samples (l2fc mRNA 0.9; protein 0.7). However, since the phosphorylation level of GSK3B was reduced (which increases its inhibitory activity on GYS2), no increased glycogen synthesis in the liver of MIDY pigs was expected nor observed (data not shown). Poorly controlled type 1 diabetic patients showed reductions in both synthesis and breakdown of hepatic glycogen ([48]; reviewed in [49]).

Altered acylcarnitine homeostasis associated with increased beta-oxidation of fatty acids and ketogenesis

Stimulated ketogenesis in MIDY pigs was evident by significantly increased plasma concentrations of beta hydroxybutyrate [9]. In hepatic ketogenesis, fatty acids are metabolized to acetyl-CoA via mitochondrial or peroxisomal beta-oxidation (reviewed in [50]). Since the mitochondrial membrane is impermeable to acyl-CoAs, they are converted into acylcarnitines by carnitine palmitoyltransferase 1A (CPT1A). An increased ratio of long-chain acylcarnitines to free carnitine [(C16 + C18)/C0] indicates increased activity of this enzyme. Hydroxyacyl-CoA dehydrogenase trifunctional multienzyme complex subunit beta (HADHB), part of a complex that catalyzes the last three steps of mitochondrial beta-oxidation of long-chain fatty acids, was significantly increased in abundance in the MIDY samples (l2fc 0.5, p = 0.0046), as was acyl-CoA synthetase medium chain family member 2B (ACSM2B), which has medium-chain fatty acid-CoA ligase activity (l2fc 1.0; p = 0.0089). The increased concentration of dicarboxylated acylcarnitines in MIDY samples suggests a higher level of omega-oxidation of fatty acids compared to WT samples (reviewed in [51]).

In addition to their role as a shuttle for long-chain fatty acids into mitochondria, acylcarnitines are important for regulating the availability of free CoA. If mitochondrial disturbances lead to formation of excess acyl-CoA esters, they can be transesterified with l-carnitine, forming acylcarnitines and free CoA. The intramitochondrial ratio of acyl-CoA to free CoA is reflected by the extramitochondrial ratio of acylcarnitines to free carnitine (C0; reviewed in [52]). The latter ratio is an interesting diagnostic parameter for mitochondrial disturbances. Studies of human liver revealed ratios of acylcarnitines to C0 of 35%/65% [53] and 44%/56% [54]. In liver samples from WT pigs, the ratio was in a similar range (37%/63%). In MIDY liver, the ratios of total acylcarnitines to C0, of short-chain acylcarnitines to C0 and of acetylcarnitine to C0 were reduced, suggesting that acyl-CoAs were efficiently processed and transesterification with C0 was not necessary to maintain the pool of free CoA. The possibility that acetylcarnitine produced to economize CoA was exported from the hepatocytes is unlikely, since plasma acetylcarnitine concentrations were not different between MIDY pigs (1.07 ± 0.56 μM) and WT controls (1.05 ± 0.34 μM) [9].

Among the ketogenic enzymes, mitochondrial 3-hydroxy-3-methylglutaryl-CoA synthase 2 (HMGCS2), which catalyzes the first reaction of ketogenesis (reviewed in [50]), was strongly increased in abundance (l2fc 2.7, p = 0.0064) in MIDY samples, associated with a corresponding increase in HMGCS2 transcript levels. Activation of HMGCS2 transcription involves the forkhead box transcription factor FOXA2 (reviewed in [50]). Insulin signaling via the PI3K-AKT pathway leads to inactivation of FOXA2 via phosphorylation and nuclear export [55]. Reduced activation of AKT in MIDY liver is therefore a likely mechanism leading to markedly increased HMGCS2 levels. Interestingly, enzymes catalyzing the subsequent steps in ketogenesis, i.e. 3-hydroxy-3-methylglutaryl-CoA lyase (HMGCL), which liberates acetoacetate from HMG-CoA, and 3-hydroxybutyrate dehydrogenase 1 (BDH1), which metabolizes acetoacetate to beta hydroxybutyrate (reviewed in [50]), were not increased in abundance in MIDY vs. WT liver samples, supporting the notion that upregulation of HMGCS2 was sufficient for stimulated ketogenesis in the liver of MIDY pigs.

Increased amino acid metabolism providing fuels for ketogenesis and gluconeogenesis

In the transcriptome of MIDY pig liver samples, the enrichment of gene sets related to amino acid metabolism was most prominent. Specifically, the abundance of transcripts for enzymes involved in the degradation of specific amino acids were significantly increased, with a similar tendency on the proteome level. In addition, the abundance of branched chain keto acid dehydrogenase E1 subunit beta (BCKDHB) was significantly increased (l2fc 0.6, p = 0.0049), and BCKDHA was slightly more abundant in MIDY samples. The BCKDH complex catalyzes the second major step in the catabolism of the branched-chain amino acids leucine, isoleucine, and valine (reviewed in [56]). Nevertheless, the concentrations of these amino acids were slightly increased in the circulation [9] and in liver samples of MIDY pigs, suggesting increased rates of protein degradation in muscle and liver (reviewed in [56]).

Amino acids are deaminated before their carbon skeletons are used as substrates for gluconeogenesis or ketogenesis. The final acceptor of the α-amino group is α-ketoglutarate. The resulting glutamate undergoes oxidative deamination, releasing ammonia that is detoxified via the urea cycle. The mRNA levels for urea cycle enzymes, i.e. carbamoyl-phosphate synthase 1 (CPS1), ornithine carbamoyltransferase (OTC), argininosuccinate synthase 1 (ASS1), and arginase 1 (ARG1) were significantly increased in MIDY samples, and considerably increased protein concentrations of OTC (l2fc 0.4, p = 0.0061) and ARG1 (l2fc 0.8, p = 0.0003) were revealed. OTC enters ammonium into the urea cycle by catalyzing the reaction between carbamoyl phosphate and ornithine to form citrulline. ARG1 hydrolyzes the final intermediate arginine to urea and ornithine [57]. The increased abundance of ARG1 may explain why arginine was not detectable in MIDY samples, while it was detected in 4/5 WT samples.

Serine/methionine imbalance and consequences for sphingomyelin synthesis and glutathione homeostasis

Notably, we found decreased serine and increased methionine concentrations in MIDY liver. The decrease of serine can be explained by its use for gluconeogenesis after metabolization to pyruvate by serine dehydratase (SDS) [58]. Indeed, SDS mRNA levels were significantly increased (l2fc 1.3, p = 0.0001) in MIDY samples. A potential consequence of lower intercellular serine in MIDY samples could be the marked decrease in sphingomyelins and its precursor ceramide, as the first and rate-limiting step of sphingolipid synthesis is affected by serine concentration [59].

Methionine catabolism is more complicated, involving the methionine cycle and the transsulfuration pathway. Within the methionine cycle (reviewed in [60]), methionine is – via the intermediates S-adenosyl methionine (SAM) and S-adenosyl homocysteine (SAH) – metabolized to homocysteine. Glycine N-methyltransferase (GNMT) that accelerates the conversion of SAM into SAH was significantly upregulated in MIDY pig liver (l2fc 1.0, p = 0.0075). In addition, increased mRNA levels for adenosylhomocysteinase (AHCY; l2fc 0.6, p = 0.0001), which converts SAH into homocysteine, were observed. Homocysteine can be used to regenerate methionine by addition of a methyl group from the folate cycle. This reaction is catalyzed by methionine synthetase (MTR) that was found to be slightly upregulated (l2fc 0.3) in MIDY liver on the mRNA level, which may – at least in part – explain the increased methionine concentrations.

In the transsulfuration pathway, homocysteine can be used as donor of its sulfur group that is combined with the carbon skeleton of serine to produce cystathionine, which is then converted to cysteine (reviewed in [61]), a component of the antioxidative tripeptide glutathione that is mainly produced in the liver (reviewed in [62]). In a human hepatoma cell line, it was shown that about 50% of the cysteine in glutathione is derived by the transsulfuration pathway [63]. The transcript levels for enzymes catalyzing the two steps of the transsulfuration pathway, i.e. cystathionine beta synthase (CBS; l2fc 0.5, p = 0.0314) and cystathionine gamma lyase (CTH; l2fc 1.3, p = 2.1e-10) were significantly upregulated in MIDY samples, and the protein levels showed a trend of increasing (CBS: l2fc 0.3, p = 0.0507; CTH: l2fc 0.4, p = 0.0505). These changes may contribute to the fact that - in spite of reduced serine concentrations - the synthesis of glutathione in MIDY liver was apparently not impaired.

Upregulation of antioxidative mechanisms

Hyperglycemia is known to increase oxidative stress in the liver as well as in many other tissues (reviewed in [64]). The abundance of glutathione S-transferase mu 2 (GSTM2), which protects against endogenous oxidative stress (reviewed in [65]), was significantly increased in MIDY samples (l2fc 2.8, p = 0.0117). Glutathione peroxidase 1 (GPX1) was also more abundant in MIDY than in WT samples (l2fc 0.4, p = 0.0078) and GPX4 levels were as a tendency increased (l2fc 0.5, p = 0.0689). In addition, the transcript abundance of SLC25A47, which encodes a liver mitochondria specific uncoupling protein that facilitates proton leak in the mitochondrial matrix and thus leads to less ATP and H2O2 production [66], was significantly increased (l2fc 1.8, p = 2.4e-19).

The abundance of isocitrate dehydrogenase 1 (IDH1), which has antioxidant properties by producing NADPH for the regeneration of glutathione [67], was also increased in MIDY samples (l2fc 0.6, p = 0.0002). IDH1 requires cytosolic isocitrate, and accordingly the transcript level for the mitochondrial citrate/isocitrate transporter, SLC25A1, was significantly increased in MIDY samples (l2fc 0.6, p = 0.0191). In addition, the mRNA level for the plasma membrane transporter of citrate/isocitrate, SLC13A5, was significantly increased (l2fc 1.0, p = 0.0001).

Increased hepatic biosynthesis and/or storage of triacylglycerides

Triacylglyceride (TAG) accumulation in MIDY liver could result from increased synthesis and storage or from altered lipoprotein metabolism. Increased levels of diacylglycerides (DAG), the TAG precursors, suggest increased production of TAG. A study in rats varying plasma fatty acid and insulin concentrations independently demonstrated that the esterification of fatty acids into hepatic triglycerides was dependent on plasma free fatty acid concentrations and largely independent of hepatic insulin action [38]. Although not significant, the trend of higher plasma concentrations of non-esterified fatty acids in MIDY vs. WT pigs may contribute to their increased hepatic TAG levels. An additional factor may be the altered expression of several apolipoproteins. The abundance of APOA1 was significantly increased in MIDY liver (l2fc 1.0, p = 0.0145). In addition, increased transcript levels for APOA4 (l2fc 1.0, 0.0056) and APOA5 (l2fc 0.8, 0.0033) were revealed. APOA1 is the main protein of high-density lipoprotein (HDL) and thus central to cholesterol metabolism. As part of several lipoproteins, APOA5 is related to TAG transport and facilitates cytosolic storage of TAG in hepatocytes [68]. Nevertheless, there was no histological evidence for fatty liver disease in MIDY pigs (data not shown), although up to 40% of adult patients with type 1 diabetes were reported to have nonalcoholic fatty liver disease (NAFLD) (reviewed in [69]). This discrepancy is most likely due to the natural resistance of pigs against fatty liver disease, even in morbid obesity [70].

ADAMTS17 overexpression and reduced expression of collagens suggesting extracellular matrix alterations

The transcript with the highest abundance increase in MIDY liver was ADAMTS17 encoding ADAM metallopeptidase with thrombospondin type 1 motif 17. Members of the ADAMTS family of secreted zinc metalloproteases execute a plethora of functions in extracellular matrix (ECM) biology (reviewed in [71]). Although ADAMTS17 is known to be expressed in liver, its upregulation in diabetes mellitus has not been described yet.

The expression of collagen type I alpha 1 chain (COL1A1) was significantly reduced in MIDY samples, both on the transcript (l2fc −0.8, p = 0.0035) and the protein level (l2fc −1.7, p = 0.0259). Increased COL1A1 expression in the liver has been observed in the context of nonalcoholic steatohepatitis (NASH) (e.g. [72]), but reduced hepatic COL1A1 in insulin-deficient diabetes mellitus has not been described. A potential mechanism is the reduced activation of AKT, which is involved in the stimulation of collagen synthesis [73]. Accordingly, the levels of several other collagen chains were reduced in MIDY samples. Liver inflammation results in fibrosis due to activation of hepatic stellate cells (Ito cells) and consequently increased collagen expression (reviewed in [74]). Compared with WT, MIDY liver revealed signs of reduced inflammatory activation, and the volume density of macrophages in the liver was not different between both groups (Supplementary Fig. 4).

Reduced inflammatory and immune-related functions

Several pathways related to immune functions were found to be less active in MIDY vs. WT liver samples. The gut-liver axis permits – via the portal circulation – interactions between gut-derived substances and hepatocytes, other liver parenchymal cells, and liver immune cells (reviewed in [75]). Trace amounts of microbial products reach the liver and are scavenged by hepatocytes and Kupffer cells (resident hepatic macrophages). Innate immune activation and regulation is guided by different classes of pattern recognition receptors (PRRs), which respond to pathogen-associated molecular patterns (PAMPs) or damage-associated molecular patterns (DAMPs) (reviewed in [76], [77]). Interestingly, our analyses revealed significantly decreased transcript and/or protein abundances of key players in liver immune activation and the inflammatory signaling network in MIDY samples. For instance, the transcript abundance for C-reactive protein (CRP), an indicator of inflammation [76], [77], was significantly decreased (l2fc −1.2, p = 2.9e-7). While increased circulating concentrations of CRP have been reported in association with insulin resistance and type 2 diabetes [78], [79], primary insulin deficiency apparently leads to downregulation of hepatic CRP expression. Furthermore, the concentration of high mobility group protein B1 (HMGB1), an early inflammatory mediator and a well-established DAMP that activates the PRR Toll-like receptor 4 (TLR4), was decreased in MIDY samples (l2fc −0.7, p = 0.0150). A number of proteins known to be up-regulated upon TLR4 stimulation [80] were found at lower abundance in MIDY vs. WT liver samples. These include proteasome activator complex subunit 2 (PSME2; l2fc −1.1, p = 0.0006), GMP reductase 1 (GMPR; l2fc −1.0, p = 0.0018), protein transport protein Sec61 subunit beta (SEC61B; l2fc −1.9, p = 0.0391), and 2′-5′-oligoadenylate synthetase 2 (OAS2; l2fc −1.3, p = 0.0498). The differences on the transcript level were less pronounced, suggesting posttranscriptional regulation of these proteins. In addition, several other proteins known to be involved in or regulated by TLR signaling were decreased in abundance in MIDY samples: Rac family small GTPase 1 (RAC1; l2fc −0.8, p = 0.0083), protein phosphatase 2 scaffold subunit A alpha (PPP2R1A; l2fc −1.0, p = 0.0034), ubiquitin conjugating enzyme E2 D2 (UBE2D2; l2fc −1.8, p = 0.0113), S100 calcium binding protein A1 (S100A1; log2 fold change −0.3, p = 0.0488), legumain (LGMN; l2fc −0.4, p = 0.0347), and mitogen-activated protein kinase 3 (MAPK3; l2fc −0.9, p = 0.0373).

TLR engagement enhances PI3K-AKT-mTOR pathway activity, which is an important component in the regulation of the inflammatory immune response (reviewed in [81], [82], [83]). A reduced activation of AKT in the liver of MIDY pigs may therefore restrain TLR-mediated PI3K-AKT-mTOR pathway signaling and consequently influence innate immune homeostasis.

Transcriptome profiling revealed signs that macrophages and lymphocytes in the liver are altered in insulin-deficient diabetes, either in quantity or in the level of activation. The transcript levels of several major histocompatibility complex class 2 (MHCII) genes (SLA-DQA1, SLA-DQB1 and SLA-DRA) were significantly reduced in MIDY samples (l2fc between −0.4 and −0.7; p-values between 0.0350 and 0.0006). In addition, the mRNA concentration of RFX5 encoding regulatory factor X5 that regulates MHCII genes [84] was significantly reduced (l2fc −0.5, p = 0.0255), as was the mRNA concentration for MHCII stabilizing CD74 (l2fc −0.7, p = 0.0025). Transcript levels for SLA-1, SLA-2 and SLA-3 of the MHC class 1 (MHCI) were also reduced in MIDY samples. Statistical significance was found for SLA-2 on the transcript (l2fc −0.7, p = 0.0144) and protein level (l2fc −1.1, p = 0.0149).

Collectively, these findings indicate that inflammatory and immune-related functions were downregulated in MIDY liver. Future studies including analyses of the gut microbiome need to uncover the causes of this unexpected observation.

In conclusion, our study provides the first multi-omics analysis of liver in insulin-deficient diabetes mellitus and identified key drivers of known functional consequences of insulin deficiency. In addition, previously unknown consequences especially for inflammatory and immune functions of the liver were revealed. The multi-omics data set generated in this study provides a valuable resource for comparative studies with other experimental or clinical data sets.

Funding

This study was supported in part by a grant from the German Federal Ministry of Education and Research (BMBF) to the German Center for Diabetes Research (DZD e.V.) and by the German Research Council (Graduate School QBM; TRR127).

Duality of interest

M.R. is owner of Lipidomix GmbH. This does not alter the author's adherence to all policies on sharing data and materials. All other authors report no potential conflicts of interest relevant to this article.

Author contributions

E.W., A.B., G.J.A., and H.B. designed the study; A.B., S.R., R.W., and E.W. established the Munich MIDY pig biobank; M.B., J.P.-M., S.K., and H.B. performed transcriptome studies; F.F., E.L., T.F., and G.J.A. proteome studies; C.P. and J.A. targeted metabolomics; M.G. and Ü.C. lipidomics; M.R. retinoid measurements; B.R. and M.H.d.A. clinical-chemical studies; M.D. signaling studies; A.B. and R.W. pathological investigations. M.B., F.F., and E.W. wrote the manuscript with contributions from all authors. All authors read and approved the final manuscript. M.B, F.F., and E.W. are the guarantors of this work and, as such, had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.

Data and resource availability

The RNA-seq data generated and analyzed during the current study are available in the GEO repository, GSE122029. The mass spectrometry proteomics data generated and analyzed during the current study have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository [85], http://proteomecentral.proteomexchange.org; PXD011536.