- Research article
- Open Access
The direction of cross affects obesity after puberty in male but not female offspring
BMC Genomics volume 16, Article number: 904 (2015)
The Erratum to this article has been published in BMC Genomics 2015 16:1018
We investigated parent-of-origin and allele-specific expression effects on obesity and hepatic gene expression in reciprocal crosses between the Berlin Fat Mouse Inbred line (BFMI) and C57Bl/6NCrl (B6N).
We found that F1-males with a BFMI mother developed 1.8 times more fat mass on a high fat diet at 10 weeks than F1-males of a BFMI father. The phenotype was detectable from six weeks on and was preserved after cross-fostering. RNA-seq data of liver provided evidence for higher biosynthesis and elongation of fatty acids (p = 0.00635) in obese male offspring of a BFMI mother versus lean offspring of a BFMI father. Furthermore, fatty acid degradation (p = 0.00198) and the peroxisome pathway were impaired (p = 0.00094). The circadian rhythm was affected as well (p = 0.00087). Among the highest up-regulated protein coding genes in obese males were Acot4 (1.82 fold, p = 0.022), Cyp4a10 (1.35 fold, p = 0.026) and Cyp4a14 (1.32 fold, p = 0.012), which hydroxylize fatty acids and which are known to be increased in liver steatosis. Obese males showed lower expression of the genetically imprinted and paternally expressed 3 (Peg3) gene (0.31 fold, p = 0.046) and higher expression of the androgen receptor (Ar) gene (2.38 fold, p = 0.068). Allelic imbalance was found for expression of ATP-binding cassette transporter gene Abca8b. Several of the differentially expressed genes contain estrogen response elements.
Parent-of-origin effects during gametogenesis and/or fetal development in an obese mother epigenetically modify the transcription of genes that lead to enhanced fatty acid synthesis and impair β-oxidation in the liver of male, but not female F1 offspring. Down-regulation of Peg3 could contribute to trigger this metabolic setting. At puberty, higher amounts of the androgen receptor and altered access to estrogen response elements in affected genes are likely responsible for male specific expression of genes that were epigenetically triggered. A suggestive lack of estrogen binding motifs was found for highly down-regulated genes in adult hepatocytes of obese F1 males (p = 0.074).
Besides nutrition and sedentary lifestyle, genetic predisposition contributes to obesity. To assess the genetic determinants of obesity, many large-scale genome-wide association studies were performed in different human populations. These studies led to the identification of about 350 genomic loci contributing to obesity . However, their effects are small and the overall contribution of all detected loci to the variance in a population was small (0.01 % to 0.34 %) [2–6]. Therefore, additional factors could contribute to obesity in an individual. Missing heritability can be hidden in genetic effects of rare alleles, epistatic interaction and gene-by-environment interaction, but can be also contributed to by parent-of-origin specific effects. Parent-of-origin effects include the Y chromosome, mitochondria, parental-specific genetic imprinting, and maternal environmental effects during pregnancy and the suckling period.
In a mouse model of juvenile obesity, we tested, if parental obesity influences the transmission of predisposition for adiposity to F1 offspring. We used the Berlin Fat Mouse Inbred line (BFMI) which carries a recessive mutation responsible for high fat deposition in particular during early development between four and ten weeks . At ten weeks BFMI males are 5.5 times as obese as C57BL/6Nrl mice (B6N).
In reciprocal crosses between BFMI and B6N mice we tested the development of obesity in F1 males and females and their response to the obese parent with respect to gene expression in the liver. The results provide evidence for regulation of genes by the influence of the direction-of-cross, which leads to obesity and fatty liver exclusively in male offspring of obese BFMI mothers.
Reciprocal crosses were generated from the Berlin Fat Mouse Inbred line BFMI860-12/Hber (BFMI) and C57BL/6NCrl (B6N). Seven BFMI males were mated with seven B6N females and six B6N males with six BFMI860 females to generate 48 F1 animals (deemed patBFMI and matBFMI, respectively). Each reciprocal F1 offspring groups consisted of 12 males and 12 females. From these animals three F1 males (matBFMI) from three BFMI females and three F1 males (patBFMI) from three B6N females were used for RNA-seq transcriptome measurements. One high fat diet raised male and female from each inbred parental strain (BFMI and B6N) were used as parental strain control animals and measured via RNA-seq as well. For the cross-fostering experiment, six BFMIxB6N and six B6NxBFMI crosses generated two additional groups of reciprocal F1 animals with a total of 88 animals, 46 with a BFMI father (23 males, 23 females) and 42 with a BFMI mother (24 males, 18 females). For cross-fostering, litters were exchanged between B6N and BFMI mothers within 12 hours after birth. The litter size ranged between four and eleven offspring. Two animals per cross-fostering group were excluded due to outlier testing (R-package “outliers”).
Parents were fed a standard breeding diet (SBD). F1-animals were fed a high fat diet (HFD) from weaning at 21 days until 70 days when offspring was euthanized and dissected. The SBD contained 16.7 MJ/kg of metabolisable energy, 11 % from fat, 36 % from proteins and 53 % from carbohydrates (V1534-000, ssniff EF R/M, Ssniff GmbH, Soest/Germany). The HFD contained 19.5 MJ/kg of metabolisable energy, 45 % from fat, 24 % from proteins and 31 % from carbohydrates (E15103–34, ssniff EF R/M, Ssniff GmbH, Soest/Germany). All experimental treatments of animals were approved by the German Animal Welfare Authorities (approval no. G016/11).
Body weight and body composition (total lean mass, total fat mass) were recorded weekly from week three and four. Body composition was measured using EchoMRI whole-body composition analyzer (Echo Medical Systems, Houston, TX, USA). At ten weeks blood glucose levels were measured prior to dissection after 2 hours of fasting using a glucometer (Contour, Bayer, Leverkusen, Germany). Liver triglycerides and proteins were measured in homogenized liver samples by commercially available kits (DC Protein Assay, Bio-Rad, München, Germany; TRIGS Triglycerides Assay, Randox, Crumlin, United Kingdom). For test statistics between the matBFMI and patBFMI, we applied a linear model using anova in R. Since litter size had a significant effect on all traits, phenotypic values were adjusted for litter size before statistic tests were performed, except for blood glucose and liver triglycerides.
RNA preparation and sequencing
RNA sequencing was performed on 10 liver samples: 1 male and 1 female B6N mouse, 1 male and 1 female BFMI mouse, 3 F1 (BFMI x B6N) males (patBFMI), and 3 F1 (B6N x BFMI) males (matBFMI). The F1 males were of non-cross-fostered mice. Tissue samples were frozen in liquid nitrogen and stored at –80 °C until RNA preparation. RNA was prepared using the RiboZero RNA library protocol. Paired end RNA sequencing was run on an Illumina HiSeq2000 machine.
Analysis of RNA-Seq data
After sequencing, bcl2fastq v1.8.4  was used for base calling and demultiplexing. Paired end reads were trimmed using Trimmomatic  using standard settings. TruSeq3-PE adapter sequences provided by Illumina were used for adapter trimming. Reads were aligned against Mus musculus GRCm38 reference genome (Ensembl) using BWA , after which GATK  was applied to do base quality score recalibration, indel realignment, duplicate removal, SNP/INDEL discovery and genotyping across all 10 samples using standard hard filtering parameters or variant quality score recalibration according to GATK Best Practices recommendations .
Raw reads from RNA-Seq data were extracted from BAM files using the GenomicFeatures package using the Mus_musculus.GRCm38.81.gtf file obtained from ensembl, overlap between genes and raw reads was calculated in ‘union’ mode using the Rsamtools package . Reads were quantile normalized between samples using the preprocessCore package. After sample quantile normalization, RPKM values for each gene were calculated, which were then log transformed ( ln(n + 1) ) to reduce the effect of extreme RPKM values in subsequent statistical tests. The biomaRt R package was used to provide additional gene annotation such as MGI descriptions (using the M. musculus ensembl database GRCm38.p4).
We analyzed the gene expression of the reciprocal F1 populations with respect to differentially expressed genes, pathways, and allele specific expression. Testing for differences between the two groups of reciprocal crosses was done using two sided t-tests. Due to limited sample sizes of three animals per F1 group, adjustment for multiple testing was not applied for differential gene expression analysis. The original transcriptome sequencing data for all ten animals tested are available at the GEO database (http://www.ncbi.nlm.nih.gov/geo/).
Genes that were differentially expressed (p < 0.05) between the two groups of F1 males with RPKM values > 0.5 in at least one of the F1 groups and gene expression differences of > = 20 % were subjected to overrepresentation analyses in KEGG pathways using innateDB  to identify pathways associated with these genes . We used the default hypergeometric algorithm and the Benjamini Hochberg algorithm for multiple testing corrections. Heatmaps were created for selected pathways using the mean gene expression per group.
Analysis of allele specific expression (ASE)
To assess ASE in the reciprocal F1 crosses, bcftools was used to perform SNPs calling on the entire population (all 10 animals: 2 BFMI, 2 B6N, 3 matBFMI, 3 patBFMI). Detected SNPs underwent a rigorous data quality checks to make sure detected SNPs showed sufficient coverage in the entire population. The individuals from the BFMI and B6N pure lines were analyzed to detect homozygous SNPs that showed different alleles between B6N and BFMI. We demanded that all SNPs found in the four parental strain animals had non-missing data (neither low confidence nor low coverage). SNPs were called using the samtools mpileup tool. Subsequently, the generated bcf file was processed using BCFtools  and filtered using the provided vcfutils.pl script for SNPs with a 100 reads minimum for the entire population (over all 10 animals). After using the entire population for SNP calling, we reanalyzed all SNPs per individual using the samtools mpileup tool to obtain DP4 measurements for each F1 individual. SNPs with low confidence (Phred probability threshold < 50) or less than 5 high quality reads in an individual were excluded from further analyses. Allele ratios alt/(ref + alt) were calculated for all remaining SNPs. Additionally, we filtered for consistency between the different individuals from the same F1 cross direction. If the standard deviation of the alt/(ref + alt) ratios across the 3 individuals per F1 group was > 20 %, the SNP was removed due to the fact that the data was not consistent across the individuals.
After filtering for high quality SNPs which are homozygous (but different) between the parental strains, \( \chi \) 2 scores were calculated for each individual. \( \chi \) 2 scores were adjusted to account for direction of the effect (< 50 % vs. > 50 %) and scores were square root transformed to minimize the impact of a single F1 individual (in the F1 groups) showing a highly significant deviation from the expected 50:50 ratio. Summed scores across all individuals in a F1 cross direction (matBFMI and patBFMI) were calculated to obtain a group based \( \chi \) 2 score, which were then subtracted from each other to get a group difference score per SNP. To account for over dispersion in the read counts obtained from the DP4 values, we permuted our group based \( \chi \) 2 scores to obtain a 5 % FDR threshold for SNPs that show significant ASE in either the patBFMI or the matBFMI group.
Analysis of transcription factor binding sites
Differentially expressed genes (p < 0.05) were submitted to transcription factor binding site (TFBS) analysis. For each gene we extracted the upstream (2000 bp) and downstream (1000 bp) region around the transcription start using the TxDb.Mmusculus.UCSC.mm10 package from bioconductor . All known mouse TFBS motifs were downloaded from the Jaspar database into R using the MotifDb package [18, 19].
Each gene region of 3000 basepairs was scanned for all estrogen receptor related motifs in the Jasper database using the Biostrings package , for each motif we recorded the number of matches within the 3000 bp region (Mmusculus-jolma2013-Esrra-2 -UniPROBE-Esrra.UP00079, -JASPAR_CORE-Esrrb-MA0141.1, -JASPAR_2014-Esrrb-MA0141.2). Additionally, position weight matrices for estrogen receptor alpha (ER), estrogen response element (ERE) and the binding site for glucocorticoid receptor (GR)[21, 22] were applied as well.
Since only one weight matrices (jolma2013-Ar) was available in the Jasper database for androgen response elements (AREs) in mouse, we created position weight matrices for androgen response elements using the motifs from Shaffer et al. for ARE, ARE2, ADR3, IR3 . We used a conservative threshold of 90 % for sensitivity and specificity. A list of the applied position weight matrices is given in Additional file 1: Table S3.
Additionally, we compared the occurrence of TFBS in both differentially and allele specific expressed genes with 1,000 permutations of random sets of genes of the respective same group size.
RNA from ten male matBFMI and ten male patBFMI was extracted from liver tissue. Total RNA was isolated from liver tissue using the NucleoSpin RNA II kit (Macherey-Nagel, Düren, Germany). One μg of total RNA was reverse-transcribed into cDNA using AccuScript High Fidelity Transcriptase (Stratagene Europe, Amsterdam, Netherlands) and Oligo (dT) primers (Ark Scientific, Sigma-Aldrich, München, Germany). The cDNA integrity was checked in a polymerase chain reaction (PCR). For semi-quantitative real-time PCR 20 ng cDNA was mixed with the qPCR mastermix of the Plus Kit for SYBR Green (Eurogentec, Cologne, Germany) and 300 nM of the primers at a final volume of 10 μl per reaction. The amplification protocol included a 10 min DNA denaturation step at 95 °C, followed by 40 cycles consisting of 20 s at 60 °C and 40 s at 72 °C. Quantification of mRNA transcript amounts was performed using the Applied Biosystems® ViiA™ 7 Real-Time PCR System and ViiA™ Software v.1.2.1. Primers for target genes Ar, Peg3 and Srebf1 were designed using Primer3web . Sequences for reference gene primers (Actb, Mrpl46, Srp72) were taken from Hruz et al. . The geometric mean of the three reference genes was used to normalize target gene expression. Relative quantification of gene expression was performed according to ABI's user guide . Further details are provided in Additional file 2: Table S4.
Results and discussion
Since BFMI animals carry a major recessive mutation leading to juvenile obesity in homozygous carriers , we expected lean F1-offspring in either cross direction. However, data provided evidence for the development of obesity in F1 males, if the mother originates from the obese BFMI line (matBFMI), while females of the same cross direction and both sexes of the reciprocal cross had normal body weight and fat mass (Table 1). Body mass did not differ significantly between offspring of the reciprocal crosses, but F1 males from BFMI mothers (matBFMI) had 1.83 times more adipose tissue mass at 10 weeks (p = 0.0037) than F1 males from BFMI fathers (patBFMI). Since the lean mass was only marginally lower in the obese matBFMI versus the lean patBFMI F1 males (0.98 fold, p = 0.054), the fat to lean mass ratio was also 1.86 times higher in matBFMI males (p = 0.0009). Moreover, matBFMI also showed ectopic fat storage in the liver with fat content being 1.30 times higher than in patBFMI males from the reciprocal cross (p = 0.084). The higher fat mass in matBFMI males began at puberty and was significant from seven weeks on (Fig. 1).
A cross-fostering experiment was performed to assess the influence of postnatal maternal effects such as milk energy content and behavior during the suckling period on later development of fat deposition in offspring. This experiment validated the results obtained before, but fat deposition was a little lower (Fig. 1). Cross-fostering led to a reduction of adipose tissue mass in male and female F1 offspring of BFMI mothers by approximately 22 % if they were raised by lean B6N foster mothers. Cross-fostered F1 males of BFMI mothers still developed 1.36 times higher fat mass (p = 0.0142), had a 1.36 times higher fat to lean ratio (p = 0.0067), 1.40 times higher hepatic triglyceride levels (p = 0.03070) and higher blood glucose levels (1.10 fold, p = 0.03073) than cross-fostered F1 males of the reciprocal cross (Table 2). These data support an effect of the suckling period on later obesity and diabetes development that has been shown before [28–30].
Differential gene expression between F1 males of BFMI mothers and BFMI fathers
Since lipid accumulation of F1 males differed between the reciprocal crosses, we investigated the direction-of-cross effect on hepatic gene expression profiles in males. Overrepresentation analyses of genes that were differentially expressed between the two groups of reciprocal F1 males at p < 0.05 and a fold-change difference of at least 20 % provided evidence for alterations of pathways for lipid metabolism, peroxisome, circadian rhythm, cellular transport processes and the endocrine system as a results of the direction of cross (Table 2, Table 3). Most genes with higher expression in the obese F1 males control the cellular utilization of fatty acids for energy production or fat synthesis including the hydrolysis of triglycerides into free fatty acids and glycerol, the conversion of free fatty acids into fatty acyl-CoA esters, the biosynthesis of saturated fatty acids or the degradation and conversion into acetyl-CoA by β-oxidation  (Fig. 2). But also genes for the transport between peroxisome and mitochondria via carnitine [32, 33] and the electron carrier generation by the citrate cycle and the respiratory chain in mitochondria were differentially expressed.
The metabolic changes in obese males were associated with a shift of the circadian rhythm, in which the gene Arntl (0.54 fold, p = 0.01), was down- and Per2 (1.32 fold, p =0.02) and Per3 (1.71 fold, p = 0.03) were up-regulated (Table 2). Disruptions of the normal circadian rhythm have been linked to obesity in Clock  and Arntl (alias Bmal1) knock-out mice and feeding experiments . A shift of the circadian rhythm in offspring as a result of in utero exposure to maternal obesity has been shown in rats .
Indications for impaired β-oxidation in the obese males were found from up-regulated expression of cytoplasmic acetyl-CoA carboxylase (Acaca, 1.20 fold, p < 0.06) and citrate lyase (Acly, fold 1.11, p < 0.08), which usually inhibits the cytosol-mitochondria transport of the fatty acid derived substrates via allosteric inhibition of Crat (1.27 fold, p = 0.01) by malonyl-CoA [37–39]. Furthermore, acyl-CoA thioesterases (Acot) are transporters that regulate the amounts of intermediates from fatty acid oxidation inside mitochondria. The mitochondrial matrix bound Acot2 (1.62 fold, p = 0.05), for example, enhances hepatic fatty acid oxidation by its transporter activity to prevent substrate overload , which otherwise could impair β-oxidation . In obese F1 males also Acot4 (1.82 fold, p < 0.02) and Acot3 (6.49 fold, p < 0.08) were up-regulated, which could indicate high amounts of such intermediates inside the mitochondrial matrix, leading to impaired fatty acid oxidation and increased levels of free fatty acids.
Genes encoding peroxisome proliferator-activated receptors (Ppara, Pparg, Ppargc1a, Ppargc1b) that promote β-oxidation and inhibit triglyceride synthesis were not differentially expressed between F1 males of reciprocal crosses. In contrast, the opponent sterol regulatory element–binding factor 1 (Srebf1) which increases fatty acid and triglyceride synthesis  was 1.17 fold higher expressed (p = 0.03) in the obese F1 males. Many Srebf1 target genes were increased as well (Ar: 2.38 fold, p = 0.07; Elovl6: 1.32 fold, p = 0.02; Fads1: 1.07 fold, p = 0.02; Fasn: 1.21 fold, p = 0.01; Gpam: 1.13 fold, p = 0.02; Hmgcr: 1.32 fold, p = 0.003; Hsd17b7: 1.14 fold, p = 0.05; Ldlr: 1.10 fold, p = 0.005; Me1: 1.16 fold, p = 0.01; Pgd: 1.15 fold, p = 0.04). Certain alleles of Srebf1 were found in humans to promote hepatic fibrogenesis , which can develop from nonalcoholic fatty liver disease. Consistently with these differences in expression, the hepatic lipase (Lipc) that degrades triglycerides was by trend down-regulated in obese F1 males (0.96 fold, p = 0.01).
In a situation, in which fatty acids are preferentially synthesized and accumulated, while their β-oxidation is impaired, the cell has to increase energy generation from alternative fuels. Previous measurements of the respiratory quotient of the BFMI inbred mice provided evidence for a higher use of carbohydrates as the main energy fuel . Succinyl-CoA, for example, is a citrate cycle intermediate that can be produced from degradation of certain amino acids but also from the propanoate pathway. The moderate up-regulation of Mgll (1.17 fold, p = 0.01), Aldh3a2 (1.11 fold, p = 0.01) and Dak (1.09 fold, p = 0.02) from the glycerolipid pathway aims at the conversion of triglycerides to precursor substrates for subsequent glycolysis. Increased carbohydrate metabolism is further supported by up-regulation of the carbohydrate transporters Slc16a7 (1.20 fold, p = 0.0007) and Slc16a12 (1.13 fold, p = 0.02) in matBFMI (Table 3, Additional file 3: Table S1).
In humans, Slc16a7 catalyzes the rapid transport of lactate, pyruvate, branched-chain oxo acids and other substrates across the plasma membrane of mitochondria and is regulated by estradiol [45, 46]. This regulation could contribute to the fact that female F1 offspring of the BFMI mother are protected against elevated hepatic triglyceride storage.
Among the highest up-regulated genes in matBFMI were the cytochrome P450, family 4, subfamily A members Cyp4a10 (1.35 fold, p = 0.03) and Cyp4a14 (1.32 fold, p = 0.01), which together with an increase of their rate limiting regulator Por (1.16 fold, p = 0.06) is an indicator for increased ω-oxidation , a prerequisite for subsequent peroxisomal β-oxidation. These genes were found to play a part in the development of liver steatosis [48, 49]. The cytochrome P450 proteins are involved in the synthesis of cholesterol, steroids, and other lipids. The proteins encoded by Cyp4a10 and Cyp4a14 hydroxylize fatty acids [50, 51]. Recently, fatty acid-hydroxy fatty acids (FAHFAs) have been discovered as a class of lipids with anti-diabetic and anti-inflammatory effects . Higher concentrations of FAHFAs in the obese F1 males due to increased expression of the two Cyp4A genes would be consistent with similar blood glucose levels in male F1 offspring of both cross directions (Table 1).
The hepatic expression of the androgen receptor (Ar) gene was higher in the obese matBFMI F1 males compared to the lean counterparts (2.38 fold, p = 0.07). The androgen receptor functions as a steroid-hormone activated transcription factor . Therefore, Ar could contribute to the sex specific obesity effect. The androgen receptor affects the liver metabolism differently in males and females . Besides the liver effects, the androgen receptor protects females against high fat diet-induced obesity and dyslipidemia . Additional studies showed that mice lacking the androgen receptor develop hepatic steatosis and insulin resistance in males but not females under high fat diet .
Among the genes that were highly differentially expressed between the two groups of F1 males (p < 0.05, fold-change difference min. 20 %, Table 3), six and two genes contain binding sites for estrogen (Acnat2, Cyp2c38, Cyp2c39, Cyp4a10, Fasn, Nlrp12) and androgen (Acnat2, Cyp4a10), respectively. Considering all differentially expressed genes under p < 0.05 estrogen (ER), androgen (AR) or glucocorticoid response (GR) elements occur in 220 up- (55 ER, 14 AR, 3 GR) and 189 down-regulated (30 ER, 11 AR, 6 GR) genes in the obese versus lean F1 males (Additional file 3: Table S1).
We compared the occurrence of these transcription factor binding sites in genes being differentially expressed (p < 0.05, >20 % fold-change difference) between the obese and the non-obese F1 males with randomly selected genes. We found that binding motifs for estrogen seemed to be under-represented in down-regulated genes (p = 0.074). However, when relaxing our restraints and define our down-regulated genes as p < 0.05 and do not filter for only the highest fold-changes (> 20 %), this under-representation of estrogen binding sites becomes highly significant (p = 0.004) while androgen binding sites tended to be over-represented in up-regulated genes (p = 0.056). Therefore, modifications of these genes as a result of epigenetic reprogramming during genetic imprinting or during pregnancy could modify their expression in response to sex-hormone production beginning at puberty.
Gene expression of hemizygous genetic elements
Differential gene expression on the X and Y chromosome and mitochondria was examined in particular because these homozygous genomic units can be polymorphic between the reciprocal crosses. Among them are the androgen receptor (Ar: 2.38 fold, p = 0.07), zinc finger protein X-linked (Zfx: 0.84 fold, p = 0.09) and polyglutamine binding protein 1 (Pqbp1: 1.07 fold, p < 0.1), which act as transcription factors. Differentially expressed X chromosomal genes that are involved in the lipid metabolism are ATP-binding cassette, sub-family D (ALD), member 1 (Abcd1: 1.22 fold, p = 0.04), HAUS augmin-like complex, subunit 7 (Haus7: 1.21 fold, p = 0.05), NAD(P) dependent steroid dehydrogenase-like (Nsdhl: 1.10 fold, p = 0.03), monocarboxylic acid transporter (Slc16a2) and glycerol kinase (Gyk) (Table 4, Additional file 3: Table S1). Zfx is well studies in cancer as a proliferation factor; it is also involved in animal size and fertility and lower gene expression was found to inhibit cell growth [57, 58]. Dysfunctional Pqbp1 was associated to lipid metabolism in C. elegans. Repression of this gene led to reduction of lipid content in mammalian primary white adipocytes and also triglyceride synthesis from fatty acids was impaired . Since no differences in expression were found for genes on the Y-chromosome, we provide further evidence that the Y chromosome can most likely ruled out as cause for differences in obesity of males. In the class of mitochondrial genes we found mitochondrial NADH dehydrogenase 1 (mt-Nd1) expressed at a lower level (0.95 fold, p = 0.02) in obese matBFMI F1 males (Table 4).
Quantitative PCR validation
We verified our gene expression fold-change results from the RNA-seq analyses performing qPCR for Ar, Peg3 and Srebf1. The qPCR results confirmed the direction of fold-changes in the obese matBFMI animals for the most genes with Ar being 3.27 fold (p = 0.004) up-regulated and Peg3 being severely down-regulated (0.36 fold, p = 0.002; Additional file 2: Table S4). Srebf1 was found to be 1.66 fold upregulated, although with lower significance (p = 0.065) as compared to the RNA-seq results (p = 0.027). It is possible that the upregulation of the androgen receptor is a reaction to low androgen levels that are known to be caused by obesity, though a study in cell culture found androgens to inhibit Peg3 expression .
Allele specific expression
Since allele-specific expression could result from parental imprinting, we exploited the sequence information of expressed genes. After rigorous quality check, we found 47,566 SNPs on autosomes (7,470 not in genes, and 40,096 SNPs in 6,260 genes). After selecting SNPs for being different between parental strains BFMI and B6N 26,101 SNPs were left (1,846 not in genes, and 24,255 SNPs in 3,961 genes). After selecting for ASE there were 1,362 SNPs left (139 not in genes, and 1,223 SNPs in 415 genes). Filtering those for concordance yielded 1,315 SNPs (137 not in genes, and 1,178 SNPs in 402 genes, see Additional file 4: Table S2). These SNPs were used to search for parent-of origin-specific expression in the liver of 10 week old males. Deviations from balanced biallelic expression can cause different phenotypes by causing total gene expression differences but also independently since allele functions can be altered due to sequence variants, altering the function of gene products. We found evidence for allelic imbalances for protein coding genes Abca8b, Commd1, H13, Harbi1, Mcts2, Peg3, Plg, Prkd3 and Zrsr1 (Table 5).
Paternal allele expression was identified for paternally expressed 3 (Peg3) and zinc finger (CCCH type), RNA binding motif and serine/arginine rich 1 (Zrsr1). Maternal allele expression was found for the gene H13. Previously, all three genes have been described to be genetically imprinted during gametogenesis and differentially expressed in offspring [61–64]. However, even if parent-of-origin specific expression is proven, the expression level of the paternal allele could differ between offspring due to cis or trans acting regulation that cannot be derived from RNA-seq data . This was the case for Peg3 in our study. Interestingly, the obese matBFMI males showed only 0.31 fold expression of Peg3 (p = 0.046), which was one of the genes with the most extreme decrease in expression. Peg3 is a Krüppel C2H2-type zinc finger protein, which plays a role in the olfactory system related sexual learning, cell proliferation and p53-mediated apoptosis [65, 66]. In zebrafish, Peg3 inhibits the Wnt pathway . In our data, lower expression of Peg3 in in obese males was accompanied by lower expressions of the Wnt target genes cadherin 1 (Cdh1: 0.75 fold, p = 0.03), epidermal growth factor receptor (Egfr: 0.85 fold, p = 0.05) and neuropilin 1 (Nrp1: 0.90 fold, p = 0.02). High Peg3 expression was recently suggested to protect from high-fat diet induced obesity in F1 offspring of reciprocal crosses between PWK and B6 . The ‘paternal transmission’ of the Peg3 effect was deduced from the observation that B6 develop obesity under high fat diet, whereas PWK mice do not. Our data suggest that different paternal alleles and expression levels of Peg3 act differently in heterozygous BFMI/B6N mice.
Very little is known about the function of the imprinted gene Zrsr1 although it was associated with the pluripotency state in induced pluripotent stem cells . Furthermore, it could be shown that trans-acting DNA binding proteins, such as methyl-CpG binding domain proteins 2 (Mecp2), methyl-CpG binding domain protein Mbd1 and Mbd3 bind at maternal and paternal Zrsr1 alleles differently . In our F1 male mice, Zrsr1 was 0.80 fold expressed in the obese F1 males (p = 0.08). The predominant expression of the maternal H13 allele was confirmed by our data and the gene was only 0.94 fold expressed (p < 0.05) in the matBFMI males.
Differences in expression of imprinted genes could be causal for regulation of down-stream genes contributing to obesity. Recently, paternal intergenerational metabolic reprogramming was found in drosophila, where paternal fat-diet led to changes in the germ cell chromatin state, transcription patterns and modified offspring metabolism. The genetic factors identified in drosophila have been verified in mouse and human studies .
We also found a lower proportion of the BFMI Abca8b allele expressed by the obese matBFMI males. While the patBFMI expressed between 58–68 %, the matBFMI expressed only 31–46 % of the BFMI allele. This corresponds to a mean difference for the alternative allele proportion of 24 % under a fold-change of 0.87 (p = 0.06) in total gene expression. ATP-binding cassette (ABC) transporter genes represent the largest family of transmembrane proteins, which arrange the transport of various molecules across all cell membranes utilizing energy from bound ATP. Most ABC genes move compounds from the cytoplasm to the outside of the cell or into an extracellular compartment (endoplasmic reticulum, mitochondria, peroxisome). ABC transporters shuttle hydrophobic compounds within the cell as part of a metabolic process or outside the cell for transport to other organs or secretion from the body. The Abca subfamily includes the largest Abc genes (several >2000 amino acids) and its members are involved in transport of vitamin A (retinol, etc.) derivatives as well as in disorders of cholesterol transport and high-density lipoproteins (HDL) biosynthesis . Our data from pathway analyzes and ASE indicate that altered functions of ABC transporters (mainly subfamily Abca and Abcd) contribute to the accumulation of body fat in the obese matBFMI males.
In reciprocal crosses between the obese line BFMI and the lean line B6N, we observed high fat deposition with and without cross-fostering beginning at puberty only in male F1 offspring of obese BFMI mothers, but not in females of the same cross and males and females of the reciprocal cross. Gene expression analysis of liver RNA provided evidence for higher expression of genes controlling fat deposition and impaired β-oxidation of lipids as well as effects on the circadian rhythm in the obese F1 males. The obese F1 males likely use more carbohydrates and amino acids for energy production than lean males. Among the extremely up-regulated genes were Cytochrome A4 genes leading to hydroxylation of fatty acids. Furthermore, typical imprinted genes were differentially expressed between the two groups of males in reciprocal crosses.
The observations led us to the conclusion that parent-of-origin effects during gametogenesis and pregnancy play a significant role for the development of obesity in later life. Postnatal maternal effects significantly influence fat deposition, but could not entirely overwrite perinatal parent-of-origin-effects. Reprogramming of gametes and metabolic imprints during pregnancy set epigenetic marks that lead to differential expression of genes in response to sex-hormone activation beginning at puberty. We suggest that the obese F1 males store more fat due to impaired β-oxidation and cellular compound transport. They compensate the lower energy production from fatty acids by an increased use of alternative fuels for the mitochondrial and non-mitochondrial energy production. The hydroxylases Cyt4a10 and Cyp4a14 are suggested to contribute to increased fatty liver, but the hydroxylated fatty acids are likely nonhazardous. The imprinted genes Peg3, Zrsr1, und H13 and the androgen receptor might play important sex-specific roles for affected down-stream signaling and metabolic pathways. We further suggest that the Y and X chromosomes as well as mitochondria play a minor role in the identified phenotypic differences.
Since the higher adipose tissue mass in F1 males of the BFMI mother was maintained during cross-fostering, parent-of-origin effects of genomic reprogramming during gametogenesis and maternal effects during pregnancy could be causal for the elevated fat accumulation.
Although, we found indicators for substrate inhibition circuits in mitochondria, the hypothesis of impaired substrate selection or mitochondrial substrate overload by fatty acids in BFMI mitochondria needs further research. Furthermore, we like to point out that allele-specific-expression analysis with RNA-seq data requires high control standards. Allele-specific expression could only be tested, if sequence variation was identified between the examined mouse strains. Since the analysis performed in this study was very conservative, the data we present here are highly reliable. Nonetheless, we miss information on genes without variation in the coding sequence and genes with very low transcript levels in liver. Moreover, the molecular mechanisms necessary for allele specific expression in adult tissues need further evaluation. A controlling function of sex hormones in metabolism was shown before and is supported by reduced counts of binding sites for estrogens in down-regulated genes in our data.
allele specific expression
androgen receptor binding site
androgen direct repeats of the hexameric half-site
androgen response element
Berlin Fat Mouse Inbred strain 860-12/Hber
estrogen related receptor a/b
estrogen receptor alpha
estrogen response element
androgen inverted repeats of hexameric half-site separated by 3 bp of spacer
maternal BFMI F1 animals (male B6N x female BFMI), obese
overrepresentation analysis (via innateDB)
paternal BFMI F1 animals (male BFMI x female B6N), non-obese
Reads Per Kilobase per Million mapped reads
transcription factor bindings site
National Genome Research Institute [http://www.genome.gov/gwastudies/]
Speliotes EK, Willer CJ, Berndt SI, Monda KL, Thorleifsson G, Jackson AU, et al. Association analyses of 249,796 individuals reveal 18 new loci associated with body mass index. Nat Genet. 2010;42:937–48.
Willer CJ, Speliotes EK, Loos RJF, Li S, Lindgren CM, Heid IM, et al. Six new loci associated with body mass index highlight a neuronal influence on body weight regulation. Nat Genet. 2009;41:25–34.
Zeggini E, Scott LJ, Saxena R, Voight BF, Marchini JL, Hu T, et al. Meta-analysis of genome-wide association data and large-scale replication identifies additional susceptibility loci for type 2 diabetes. Nat Genet. 2008;40:638–45.
Locke AE, Kahali B, Berndt SI, Justice AE, Pers TH, Day FR, et al. Genetic studies of body mass index yield new insights for obesity biology. Nature. 2015;518:197–206.
Shungin D, Winkler TW, Croteau-Chonka DC, Ferreira T, Locke AE, Mägi R, et al. New genetic loci link adipose and insulin biology to body fat distribution. Nature. 2015;518:187–96.
Wagener A, Schmitt AO, Aksu S, Schlote W, Neuschl C, Brockmann GA. Genetic, sex, and diet effects on body weight and obesity in the Berlin Fat Mouse Inbred lines. Physiol Genomics. 2006;27:264–70.
bcl2fastq Conversion. 2013.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754–60.
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303.
DePristo M a, Banks E, Poplin R, Garimella K V, Maguire JR, Hartl C, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet 2011, 43:491–8.
Morgan, Martin, Pages, Herve, Obenchain, Valerie HN: Rsamtools: Binary alignment (BAM), FASTA, variant call (BCF), and tabix file import. http://bioconductor.org/packages/release/bioc/html/Rsamtools.html 2014:1.
Breuer K, Foroushani AK, Laird MR, Chen C, Sribnaia A, Lo R, et al. InnateDB: systems biology of innate immunity and beyond--recent updates and continuing curation. Nucleic Acids Res. 2013;41(Database issue):D1228–33.
Ogata H, Goto S, Fujibuchi W, Kanehisa M: Computation with the KEGG pathway database. Biosystems, 47:119–28.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Carlson M. Annotation package for TxDb object(s). 2015.
Bryne JC, Valen E, Tang M-HE, Marstrand T, Winther O, da Piedade I, et al. JASPAR, the open access database of transcription factor-binding profiles: new content and tools in the 2008 update. Nucleic Acids Res. 2008;36(Database issue):D102–6.
Shannon P. MotifDb: An Annotated Collection of Protein-DNA Binding Sequence Motifs. 2014.
DebRoy HP and PA and RG and S: Biostrings: String objects representing biological sequences, and matching algorithms.
Matys V, Kel-Margoulis OV, Fricke E, Liebich I, Land S, Barre-Dirrie A, et al. TRANSFAC and its module TRANSCompel: transcriptional gene regulation in eukaryotes. Nucleic Acids Res. 2006;34(Database issue):D108–10.
Lin C-Y, Vega VB, Thomsen JS, Zhang T, Kong SL, Xie M, et al. Whole-genome cartography of estrogen receptor alpha binding sites. PLoS Genet. 2007;3, e87.
Shaffer PL, Jivan A, Dollins DE, Claessens F, Gewirth DT. Structural basis of androgen receptor binding to selective androgen response elements. Proc Natl Acad Sci U S A. 2004;101:4758–63.
Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, et al. Primer3--new capabilities and interfaces. Nucleic Acids Res. 2012;40, e115.
Hruz T, Wyss M, Docquier M, Pfaffl MW, Masanetz S, Borghi L, et al. RefGenes: identification of reliable and condition specific reference genes for RT-qPCR data normalization. BMC Genomics. 2011;12:156.
Guide to Performing Relative Quantitation of Gene Expression Using Real-Time Quantitative PCR
Schäfer N, Yu Z, Wagener A, Millrose MK, Reissmann M, Bortfeldt R, et al. Changes in metabolite profiles caused by genetically determined obesity in mice. Metabolomics. 2014;10:461–72.
Hager R, Cheverud JM, Wolf JB. Maternal effects as the cause of parent-of-origin effects that mimic genomic imprinting. Genetics. 2008;178:1755–62.
Ozanne SE, Hales CN. Lifespan: catch-up growth and obesity in male mice. Nature. 2004;427:411–2.
Reifsnyder PC, Churchill G, Leiter EH. Maternal environment and genotype interact to establish diabesity in mice. Genome Res. 2000;10:1568–78.
Dirkx R, Meyhi E, Asselberghs S, Reddy J, Baes M, Van Veldhoven PP. Beta-oxidation in hepatocyte cultures from mice with peroxisomal gene knockouts. Biochem Biophys Res Commun. 2007;357:718–23.
Van Roermund CWT, Ijlst L, Wagemans T, Wanders RJA, Waterham HR. A role for the human peroxisomal half-transporter ABCD3 in the oxidation of dicarboxylic acids. Biochim Biophys Acta. 1841;2014:563–8.
Westin MAK, Hunt MC, Alexson SEH. Short- and medium-chain carnitine acyltransferases and acyl-CoA thioesterases in mouse provide complementary systems for transport of beta-oxidation products out of peroxisomes. Cell Mol Life Sci. 2008;65:982–90.
Turek FW, Joshu C, Kohsaka A, Lin E, Ivanova G, McDearmon E, et al. Obesity and metabolic syndrome in circadian Clock mutant mice. Science. 2005;308:1043–5.
Arble DM, Bass J, Laposky AD, Vitaterna MH, Turek FW. Circadian timing of food intake contributes to weight gain. Obesity (Silver Spring). 2009;17:2100–2.
Borengasser SJ, Kang P, Faske J, Gomez-Acevedo H, Blackburn ML, Badger TM, et al. High fat diet and in utero exposure to maternal obesity disrupts circadian rhythm and leads to metabolic programming of liver in rat offspring. PLoS One. 2014;9.
McGarry JD, Brown NF. The mitochondrial carnitine palmitoyltransferase system. From concept to molecular analysis. Eur J Biochem. 1997;244:1–14.
Muoio DM. Metabolic Inflexibility: When Mitochondrial Indecision Leads to Metabolic Gridlock. Cell. 2014;159:1253–62.
Kerner J, Hoppel C. Fatty acid import into mitochondria. Biochim Biophys Acta. 2000;1486:1–17.
Moffat C, Bhatia L, Nguyen T, Lynch P, Wang M, Wang D, et al. Acyl-CoA thioesterase-2 facilitates mitochondrial fatty acid oxidation in the liver. J Lipid Res. 2014;55:2458–70.
Van Eunen K, Simons SMJ, Gerding A, Bleeker A, den Besten G, Touw CML, et al. Biochemical competition makes fatty-acid β-oxidation vulnerable to substrate overload. PLoS Comput Biol. 2013;9, e1003186.
Teodoro JS, Rolo AP, Palmeira CM. Hepatic FXR: key regulator of whole-body energy metabolism. Trends Endocrinol Metab. 2011;22:458–66.
Krawczyk M, Grünhage F, Lammert F. Identification of combined genetic determinants of liver stiffness within the SREBP1c-PNPLA3 pathway. Int J Mol Sci. 2013;14:21153–66.
Meyer CW, Wagener A, Rink N, Hantschel C, Heldmaier G, Klingenspor M, et al. High energy digestion efficiency and altered lipid metabolism contribute to obesity in BFMI mice. Obesity (Silver Spring). 2009;17:1988–93.
Bektic J, Wrulich OA, Dobler G, Kofler K, Ueberall F, Culig Z, et al. Identification of genes involved in estrogenic action in the human prostate using microarray analysis. Genomics. 2004;83:34–44.
Dahm AEA, Eilertsen AL, Goeman J, Olstad OK, Ovstebø R, Kierulf P, et al. A microarray study on the effect of four hormone therapy regimens on gene transcription in whole blood from healthy postmenopausal women. Thromb Res. 2012;130:45–51.
den Hil EF H-v, Keijer J, Bunschoten A, Vervoort JJM, Stankova B, Bekkenkamp M, et al. Quercetin induces hepatic lipid omega-oxidation and lowers serum lipid levels in mice. PLoS One. 2013;8:e51588.
Hardwick JP. Cytochrome P450 omega hydroxylase (CYP4) function in fatty acid metabolism and metabolic diseases. Biochem Pharmacol. 2008;75:2263–75.
Hardwick JP, Osei-Hyiaman D, Wiland H, Abdelmegeed MA, Song B-J. PPAR/RXR Regulation of Fatty Acid Metabolism and Fatty Acid omega-Hydroxylase (CYP4) Isozymes: Implications for Prevention of Lipotoxicity in Fatty Liver Disease. PPAR Res. 2009;2009:952734.
Holla VR, Adas F, Imig JD, Zhao X, Price E, Olsen N, et al. Alterations in the regulation of androgen-sensitive Cyp 4a monooxygenases cause hypertension. Proc Natl Acad Sci U S A. 2001;98:5211–6.
Enriquez A, Leclercq I, Farrell GC, Robertson G. Altered expression of hepatic CYP2E1 and CYP4A in obese, diabetic ob/ob mice, and fa/fa Zucker rats. Biochem Biophys Res Commun. 1999;255:300–6.
Yore MM, Syed I, Moraes-Vieira PM, Zhang T, Herman MA, Homan EA, et al. Discovery of a Class of Endogenous Mammalian Lipids with Anti-Diabetic and Anti-inflammatory Effects. Cell. 2014;159:318–32.
Askew EB, Gampe RT, Stanley TB, Faggart JL, Wilson EM. Modulation of androgen receptor activation function 2 by testosterone and dihydrotestosterone. J Biol Chem. 2007;282:25801–16.
Rando G, Wahli W. Sex differences in nuclear receptor-regulated liver metabolic pathways. Biochim Biophys Acta. 1812;2011:964–73.
Fagman JB, Wilhelmson AS, Motta BM, Pirazzi C, Alexanderson C, De Gendt K, et al. The androgen receptor confers protection against diet-induced atherosclerosis, obesity, and dyslipidemia in female mice. FASEB J. 2014.
Lin H-Y, Yu I-C, Wang R-S, Chen Y-T, Liu N-C, Altuwaijri S, et al. Increased hepatic steatosis and insulin resistance in mice lacking hepatic androgen receptor. Hepatology. 2008;47:1924–35.
Fang Q, Fu W-H, Yang J, Li X, Zhou Z-S, Chen Z-W, et al. Knockdown of ZFX suppresses renal carcinoma cell growth and induces apoptosis. Cancer Genet. 2014.
Luoh SW, Bain PA, Polakiewicz RD, Goodheart ML, Gardner H, Jaenisch R, et al. Zfx mutation results in small animal size and reduced germ cell number in male and female mice. Development. 1997;124:2275–84.
Takahashi K, Yoshina S, Masashi M, Ito W, Inoue T, Shiwaku H, et al. Nematode homologue of PQBP1, a mental retardation causative gene, is involved in lipid metabolism. PLoS One. 2009;4, e4104.
Ulrix W, Swinnen JV, Heyns W, Verhoeven G. Androgens down-regulate the expression of the human homologue of paternally expressed gene-3 in the prostatic adenocarcinoma cell line LNCaP. Mol Cell Endocrinol. 1999;155:69–76.
Blake A, Pickford K, Greenaway S, Thomas S, Pickard A, Williamson CM, et al. MouseBook: an integrated portal of mouse resources. Nucleic Acids Res. 2010;38(Database issue):D593–9.
Morison IM, Paton CJ, Cleverley SD. The imprinted gene and parent-of-origin effect database. Nucleic Acids Res. 2001;29:275–6.
Crowley JJ, Zhabotynsky V, Sun W, Huang S, Pakatci IK, Kim Y, et al. Analyses of allele-specific gene expression in highly divergent mouse crosses identifies pervasive allelic imbalance. 2015(January).
Swaney WT, Curley JP, Champagne FA, Keverne EB. Genomic imprinting mediates sexual experience-dependent olfactory learning in male mice. Proc Natl Acad Sci U S A. 2007;104:6084–9.
Broad KD, Curley JP, Keverne EB. Increased apoptosis during neonatal brain development underlies the adult behavioral deficits seen in mice lacking a functional paternally expressed gene 3 (Peg3). Dev Neurobiol. 2009;69:314–25.
Jiang X, Yu Y, Yang HW, Agar NYR, Frado L, Johnson MD. The imprinted gene PEG3 inhibits Wnt signaling and regulates glioma growth. J Biol Chem. 2010;285:8472–80.
Morita S, Horii T, Kimura M, Arai Y, Kamei Y, Ogawa Y, et al. Paternal allele influences high fat diet-induced obesity. PLoS One. 2014;9, e85477.
Gao S, Chang G, Tian J, Gao S, Cai T. Identification of the new gene Zrsr1 to associate with the pluripotency state in induced pluripotent stem cells (iPSCs) using high throughput sequencing technology. Genomics Data. 2014;2:73–7.
Fournier C, Goto Y, Ballestar E, Delaval K, Hever AM, Esteller M, et al. Allele-specific histone lysine methylation marks regulatory regions at imprinted mouse genes. EMBO J. 2002;21:6560–70.
Öst A, Lempradl A, Casas E, Weigert M, Tiko T, Deniz M, et al. Paternal Diet Defines Offspring Chromatin State and Intergenerational Obesity. Cell. 2014;159:1352–64.
Dean M, Rzhetsky A, Allikmets R. The human ATP-binding cassette (ABC) transporter superfamily. Genome Res. 2001;11:1156–66.
The project was funded by the DFG (BR 1285/12-1). We also acknowledge the support of the Max Planck Society.
We declare no conflicts of interest or competing interests.
SK designed the experiments, analyzed the phenotypes, performed statistical analyses, and drafted the manuscript. DA carried out NGS-bioinformatics, gene expression and ASE analysis and performed statistical analysis. SH and JT participated in analyses of phenotypes, data analyses and qPCR. MLY, VA and TR performed the NGS and participated in sequencing analyses and gene expression statistics. HL participated in the design of the study. GAB conceived of the study, and participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.
Stefan Kärst, Danny Arends and Sebastian Heise contributed equally to this work.
An erratum to this article is available at http://dx.doi.org/10.1186/s12864-015-2204-y.
Table S3. Position weight matrices for scanned TFBS. (XLSX 12 kb)
Table S4. Validation of RNA-seq derived fold-changes via qPCR. (XLSX 10 kb)
Table S1. All genes being differentially expressed between both reciprocal crosses (p < 0.1). (XLSX 96 kb)
Table S2. Sequence variants measured as monoallelic in one or both reciprocal crosses. (XLSX 814 kb)