Maternal nutrition induces gene expression changes in fetal muscle and adipose tissues in sheep

Ethics statement

The College of Agriculture and Life Sciences Animal Care and Use Committee of the University of Wisconsin-Madison approved the procedures used in this study.

Animals, experimental design, and maternal diets

Tissue collection

Ewes were necropsied on day 130 ± 1 of gestation. A total of 26 fetuses were removed from 15 dams. From the fetuses three different tissues were collected including one muscle, longissimus dorsi muscle from the left side of the fetus, and two adipose tissues, perirenal adipose depot and subcutaneous adipose depot near the shoulder blade. These tissues were immediately frozen at -80°C until RNA extraction was performed.

RNA extraction and preparation of pooled samples

Total RNA was extracted from tissue samples using the RNeasy Mini kit (QIAGEN, Valencia, CA, USA) and then treated with DNase-free DNase Set (QIAGEN) to avoid genomic DNA amplification. Concentrations and OD_260/280 ratios of RNA samples were measured with the Nanodrop ND-1000 spectrophotometer (Nanodrop Technologies, Montchanin, DE).

The RNA samples from the 26 fetuses were pooled to generate four biological replicates per maternal diet and tissue. In particular, per each diet and tissue, two RNA pools were created from male fetuses and two RNA pools were created from female fetuses. The description of the pools including number of fetuses and dams can be found in Additional file 1.

Library generation and RNA sequencing

A total of 50 ng of RNA from each pool was used to prepare sequencing libraries following Illumina’s mRNA-Seq protocol. Libraries were sequenced with Illumina’s HiSeq 2000 at the Biotechnology Center in the University of Wisconsin-Madison. The 36 libraries (i.e., 12 libraries per tissue, three tissues in total) were barcoded, multiplexed, and then sequenced. A read was defined as a 100 bp cDNA fragment sequenced from a single end. Approximately 30 million reads were sequenced from each library.

Mapping reads to the reference genome

Raw sequence reads were mapped directly to the reference sheep genome (Oar_v3.1) using the software package Tophat (v2.0.4) [33]. A two-step approach was used in order to maximize sensitivity to splice junction discovery: the first alignment was performed in each of the samples independently, and then novel splice junctions plus known splice junctions from Ensembl annotation were combined and supplied to Tophat for a second alignment. This approach allows for full utilization of the novel junctions identified in the samples. A maximum of two mismatches were allowed and reads that mapped equally well to more than 40 genomic locations were discarded.

Assembly of transcripts and estimation of abundance

The resulting alignments were used to reconstruct transcript models using the software package Cufflinks (v2.1.1) [34]. This program constructs a parsimonious set of transcript models that best explain the read alignments observed in the samples. In addition, the tool cuffmerge was used for merging together each of the assemblies with the reference sheep annotation file in order to combine novel transcripts with known annotated transcripts. This strategy maximizes the overall quality of the final assembly. Finally, abundances of transcripts were scaled via the median of the geometric means of fragment counts across all libraries (as described in [35]) and corrected for sequence bias [36] in order to improve expression estimates.

Overall gene expression analysis

Differentially expressed genes were detected using Cuffdiff, a companion tool of Cufflinks. Cuffdiff performs differential analysis at both gene and transcript-level resolution and controls for both variability across replicates and uncertainty in abundance expression estimates caused by ambiguously mapped reads [37]. In addition, Cuffdiff controls for cross-replicate variability and also read-mapping ambiguity by using a model for fragment counts based on the beta negative binomial distribution [37]. Here, the analyses were performed using the default settings for all parameters. Finally, to account for multiple hypothesis testing and control the false discovery rate, P-values reported by Cuffdiff were corrected using the Benjamini-Hochberg procedure [38].

Validation of differentially expressed genes

Gene set enrichment analysis

The significant enrichment of Gene Ontology (GO) terms and KEGG pathways with genes differentially expressed between maternal diets was analyzed using Fisher’s exact test, a test of proportions based on the cumulative hypergeometric distribution [40]. For each comparison of interest, genes that showed a FDR ≤0.20 and had ENSEMBL annotations were tested against the background set of all genes with ENSEMBL annotations. The procedure proposed by Benjamini and Hochberg [38] was applied in order to account for multiple testing. Functional categories with FDR ≤0.05 were considered significant. These analyses were performed using the procedure FatiGO [41], implemented on the platform Babelomics [42], and using the Bioconductor packages org.Bt.eg.db, GO.db, and KEGG.db that are available in the R language/environment.

Sequencing of the transcriptome of different fetal sheep tissues

To assess the effects of three different maternal diets on the transcriptome of three different fetal tissues in sheep, a total of 36 pooled samples (12 pooled samples per tissue with 4 biological replicates per diet) were analyzed using RNA-sequencing. Additional file 2 shows the results of the alignment of the sequencing reads to the sheep reference genome. Approximately, 30 million reads were sequenced for each sample. Sequencing reads were aligned against the recent ovine reference genome (Oar_v3.1) using the software package Tophat. On average, 83% of the total reads were successfully mapped by allowing no more than two mismatches and restricting the alignments to at most 40 genomic locations (Additional file 2). Among the aligned reads, approximately 86% and 65% were mapped to unique genomic regions in muscle and adipose tissues, respectively. Sequencing data can be accessed by GEO with the accession number GSE62938.

Differential gene expression in longissimus dorsi muscle tissue

A total of 18,393 genes were tested for differential expression between maternal diets in fetal longissimus dorsi muscle tissue. Marked differences in gene expression were found between CN diet and the other two maternal diets. Controlling FDR at 0.05, a total of 224, 823 and 29 genes showed differential expression between CN vs. DG, CN vs. HY, and DG vs. HY, respectively (Figure 1A). Notably, the expression of most of these genes was decreased by maternal CN diet. In particular, 168 out of 224 genes, and 600 out of 823 genes showed higher expression in fetuses from DG and HY compared with those derived from CN diet, respectively (Figure 1A). Maternal diets DG and HY induced a very similar pattern of gene expression in fetal muscle; indeed, only 29 genes showed differential expression at a FDR ≤0.05 between these two diets and none of them were significant at a FDR ≤0.01.

Figure 1B shows the overlap between genes that were differentially expressed in each of the three pairwise comparisons. Interestingly, a total of 166 genes differed simultaneously between CN diet and the other two maternal diets. Many of these genes are directly involved in embryonic and fetal development [e.g., ankyrin repeat domain 11 (ANKRD11), axin 1 (AXIN1), epsin 1 (EPN1), epsin 2 (EPN2), guanine nucleotide binding protein (G protein) alpha 12 (GNA12), methyl-CpG binding domain protein 3 (MBD3), WD and tetratricopeptide repeats 1 (WDTC1)], skeletal muscle cell and tissue differentiation [e.g., ankyrin repeat domain 1 (ANKRD1), B cell CLL/lymphoma 9-like (BCL9L), histone cell cycle regulator (HIRA), myogenic differentiation 1 (MYOD1), sarcoglycan alpha (SGCA), tripartite motif containing 72, E3 ubiquitin protein ligase (TRIM72)], and muscle myosin complex and sarcomere organization [e.g., myosin, heavy chain 13, skeletal muscle (MYH13), nebulin-related anchoring protein (NRAP), parvin beta (PARVB), synaptopodin 2-like (SYNPO2L), titin-cap (TCAP)]. Importantly, the expression of all these genes was decreased in fetuses derived from dams fed CN diet (FDR ≤0.05).

Differential gene expression in subcutaneous adipose tissue

Noticeable differences in gene expression in fetal subcutaneous adipose tissue were found between DG maternal dietary feed source and the other two maternal diets. Of the 17,449 genes evaluated for differential expression, a total of 745, 705 and 43 genes showed differential expression (FDR ≤0.05) between DG vs. CN, DG vs. HY, and CN vs. HY, respectively (Figure 1C). Most of the significant genes showed higher expression in fetuses derived from DG-fed dams compared with CN-fed dams, while similar proportion of highly and lowly expressed genes was found in the comparison between DG and HY diets (Figure 1C). Interestingly, only 43 genes showed differential expression at a FDR ≤0.05 between CN and HY dietary treatments indicating very similar pattern of overall gene expression (Figure 1C). None of these 43 genes was significant at a FDR ≤0.01.

Among the differentially-expressed genes (FDR ≤0.05) in DG vs. CN and DG vs. HY, a total of 304 genes were detected in both pairwise comparisons (Figure 1D). Many of these genes are closely related to embryonic and fetal development [e.g., angiomotin (AMOT), integrin, beta 6 (ITGB6), shisa family member 2 (SHISA2), SRY (sex determining region Y)-box 6 (SOX6), transforming growth factor, beta 2 (TGFB2), ubiquitin protein ligase E3 component n-recognin 3 (UBR3), WD and tetratricopeptide repeats 1 (WDTC1)], adipose tissue development [e.g., acetoacetyl-CoA synthetase (AACS), AT rich interactive domain 5B (MRF1-like) (ARID5B), patatin-like phospholipase domain containing 3 (PNPLA3)], cholesterol and fatty acid biosynthetic process [e.g., MLX interacting protein-like (MLXIPL), protein kinase, AMP-activated, alpha 2 catalytic subunit (PRKAA2), protein kinase, AMP-activated, beta 2 non-catalytic subunit (PRKAB2)], and metabolism of lipids and lipoproteins [e.g., angiopoietin-like 4 (ANGPTL4), cell death-inducing DFFA-like effector A (CIDEA), cell death-inducing DFFA-like effector C (CIDEC), carnitine palmitoyltransferase 1A (CPT1A), nuclear receptor subfamily 1, group D, member 2 (NR1D2), phospholipase A2, group VII (platelet-activating factor acetylhydrolase, plasma) (PLA2G7)]. The majority of these genes showed higher expression in the fetuses derived from dams fed DG diet.

Gene set enrichment analysis was performed to further characterize the processes that could be regulated differentially in fetal subcutaneous fat due to maternal nutrition. Interestingly, we found three functional categories, GO terms dephosphorylation, phosphatase activity, and phosphoric ester hydrolase activity that were enriched (FDR ≤0.05) with significant genes detected in both pairwise comparisons involving DG maternal diet, i.e., DG vs. CN and DG vs. HY comparisons (Table 4). The majority of the genes in these pathways showed higher expression in fetuses from dams fed DG diet; this pattern is more noticeable in the GO term dephosphorylation (Table 4). Moreover, 10 different GO terms that belong to the biological process domain showed significant enrichment (FDR ≤0.05) of differentially expressed genes detected between DG and CN maternal diets (Figure 2A). These functional categories are mainly involved in embryo and fetal development (e.g., tissue development, organ development, multicellular organismal development) and ion transport (e.g., cation transport, calcium ion transport, metal ion transport). Remarkably, the expression of most of the genes in these pathways was increased by maternal DG dietary treatment (Figure 2A).

GO term	DG vs CN			DG vs HY
	No. genes	No. DEG (↑ DG)	q-value	No. genes	No. DEG (↑ DG)	q-value
Dephosphorylation (0016311)	106	24 (21)	0.025	108	29 (21)	0.035
Phosphatase activity (0016791)	156	29 (18)	0.050	160	36 (21)	0.030
Phosphoric Ester Hydrolase Activity (0042578)	199	36 (25)	0.043	205	46 (27)	0.008

Diets fed to dams: HY = ad libitum fed alfalfa haylage; CN = limit-fed whole shell corn; DG = limit-fed corn dried distiller’s grains.

Differential gene expression in perirenal adipose tissue

A total of 18,310 genes were evaluated for differential expression between maternal diets in fetal perirenal adipose tissue. Among the genes analyzed, 208, 109 and 27 genes showed differential expression (FDR ≤0.05) between DG vs. CN, DG vs. HY, and CN vs. HY, respectively (Figure 1E). Marked differences were found only between DG and CN maternal diets; none of the other two pairwise comparisons detected significant genes at FDR ≤0.01. Most of the significant genes showed higher expression in fetuses derived from dams fed DG diet.

We further characterized the 208 genes that showed differential expression (FDR ≤0.05) between DG and CN maternal dietary feed sources. Many of these genes are associated with embryo and fetal development [e.g., integrin, beta 4 (ITGB4), latent transforming growth factor beta binding protein 4 (LTBP4), methyl-CpG binding domain protein 3 (MBD3), retinol binding protein 4 (RBP4), secreted frizzled-related protein 5 (SFRP5)], adipogenesis [e.g., CCAAT/enhancer binding protein (C/EBP), beta (CEBPB), delta-like 1 homolog (Drosophila) (DLK1), lamin A/C (LMNA)], and lipid metabolic process [e.g., acetyl-CoA acyltransferase 1 (ACAA1), apolipoprotein D (APOD), apolipoprotein E (APOE), low density lipoprotein receptor-related protein 1 (LRP1), prostaglandin I2 (prostacyclin) synthase (PTGIS)]. All these genes were up-regulated in the perirenal adipose tissue of fetuses derived from dams fed the DG diet (FDR ≤0.05).

Gene set enrichment analysis revealed 10 GO terms from the biological process domain significantly enriched with genes differentially expressed between DG and CN maternal diets (Figure 2B). These GO terms are closely related to fetal development (e.g., tissue development, organ development, anatomical structure development), cellular adhesion (e.g., cell adhesion, cell-substrate adhesion) and Wnt signaling pathway. The expression of most genes in these functional categories was increased by maternal DG diet (Figure 2B). Interestingly, many of these significant GO terms, especially those related to tissue and organ development, were also enriched with significant genes detected between DG and CN maternal diets in fetal subcutaneous fat tissue.

Validation of overall gene expression

To validate genes found to be significant in the RNA-Seq analysis, four differentially expressed genes detected in fetal subcutaneous adipose tissue (ANKRD1, KREMEN1, MURC, and SYNPO2L) were selected and their expression was assessed using qRT-PCR. RNA-Seq analysis revealed higher expression of these four genes in subcutaneous fat in fetuses derived from DG-fed dams compared with fetuses from CN-fed dams. Figure 3 displays the fold differences in gene expression measured by both RNA-Seq and qRT-PCR. The four genes showed similar patterns of mRNA abundance with both methods (Figure 3).

Maternal nutrition impacts fetal tissue development including myogenesis and adipogenesis

Remarkable gene expression differences in longissimus dorsi muscle were found between fetuses from CN-fed dams and fetuses from either DG-fed or HY-fed dams (Figure 1A-B). Many of the most significantly altered genes are associated with skeletal muscle development. For instance, MYOD1 encodes a nuclear protein that belongs to the myogenic basic helix-loop-helix transcription factor family, and plays key roles in skeletal muscle cell specification and differentiation [45]. ANKRD1, a member of the muscle ankyrin repeat protein family, is a transcription factor involved in the signaling cascade during myogenic differentiation, myofibrillar assembly, and muscle remodeling [46]. HIRA is a replication-independent histone chaperone that is linked to transcription and various developmental processes, including the activation of relevant genes during skeletal myogenesis [47]. BCL9 is an essential component of canonical Wnt signaling cascade that mediates myogenic differentiation during muscle development and regeneration [48]. Importantly, the expression of all these genes was decreased in fetuses derived from dams fed the CN diet. In addition, relevant genes associated with muscle myosin complex and sarcomere organization also showed differential gene expression due to CN maternal diet. For instance, NRAP is a muscle-specific scaffolding protein that is involved in myofibril assembly and sarcomere organization [49]. PARVB is localized at regions of the sarcolemma in skeletal muscle and plays both physical and signaling roles in muscle cells [50]. TCAP is a key component of the Z-disk in striated and cardiac muscle and has a relevant role in sarcomere assembly and integrity [51]. SYNPO2L is another component of the Z-disk with important roles in heart and skeletal muscle function and development [52]. These genes were also down-regulated in the muscle of fetuses derived from CN-fed dams. Overall, these results suggest that a prepartum maternal CN diet, i.e., starch-based diet, altered myogenesis, myofibrillogenesis and sarcomere organization in fetal longissimus dorsi muscle. Crude protein intake of ewes fed CN was less than that for ewes fed DG and HY. Interestingly, previous studies have shown that low protein maternal diets, compared against high protein diets fed at isocaloric intakes, adversely affected myogenesis, muscle differentiation, and muscle fiber number in newborn and fattening piglets [53, 54]. Although it seems that a high-starch maternal diet during pregnancy does not affect muscle weights in sheep progeny at birth [31], the possible effect of this diet on morphological and histological parameters of the muscle remains unknown and warrants future studies.

Fetal adipose tissues also showed notable gene expression changes due to different maternal dietary treatments (Figure 1C-F). Among the most significant genes, many are closely related to adipose tissue development. For instance, ARID5B, a member of the AT-rich interaction domain family of DNA binding proteins, regulates the transcription of key genes involved in adipogenesis, including adipogenic transcription factors C/EBPα, and PPARγ [55]. CEBPB, a member of CCAAT/enhancer binding protein family, encodes a transcription factor whose expression is induced directly by adipogenic hormones, and which plays essential early catalytic roles in the adipocyte differentiation program [56]. LMNA encodes an intermediate filament protein localized in the nuclear lamina and has crucial roles in early post-natal differentiation and maturation of adipose tissue [57]. All these genes showed higher expression in fat tissues from fetuses derived from DG-fed dams. In addition, gene set enrichment analysis revealed that GO functional categories related to tissue development and organogenesis, such as tissue development, organ development, and multicellular organismal development, were significantly enriched with differentially expressed genes (Figure 2). These gene expression changes were more pronounced between DG and CN diets, and in general, genes in these GO terms showed higher expression in fetuses derived from DG-fed dams. It is important to note that dams consuming the DG prepartum diet had greater protein and fat intake than dams consuming the CN diet, and hence our findings suggest that prepartum diets with greater fat and protein content may alter fetal adipogenesis and tissue development. Previous studies have investigated the impact that different prepartum diets have on fetal adipose tissues. Larson et al. [17] have shown that maternal protein supplementation during late gestation increased overall percentage of body fat in steer progeny. Similarly, Underwood et al. [58] have shown that a low protein maternal diet negatively impact subcutaneous adipose tissue deposition in beef. In contrast, Rehfeldt et al. [53] reported that either high or low protein maternal diets, fed at isocaloric intakes during pregnancy, reduced both subcutaneous and perirenal fat content in newborn piglets. Furthermore, there is evidence that ingestion of a high-fat diet during pregnancy induced hypertrophy and hyperplasia of fetal adipose tissue [59]. Interestingly, Radunz et al. [15] have shown that a prepartum diet based on dried distillers grains, similar to our DG maternal diet, tended to increase overall internal fat, including perirenal fat, in sheep progeny.

Of particular interest, the term Wnt receptor signaling pathway showed an overrepresentation of significant genes in fetal perirenal fat tissue (Figure 2B). The seven genes that revealed differential expression were up-regulated in fetuses from DG-fed dams. Interestingly, these seven genes, AXIN1, DKK3, LRP1, LRP4, SFRP1, SFRP2, and SFRP4, are all negative regulators of the Wnt receptor signaling pathway [60, 61]. It is well established that the Wnt receptor signaling pathway inhibits adipogenesis and hence the repression of this pathway is crucial for adipocyte differentiation and fat tissue development [61, 62]. Therefore, these findings provide further evidence of the potential impact of a DG maternal diet on fetal adipose tissue development.

Maternal nutrition alters metabolic processes in fetal tissues

Prepartum maternal diets impacted the energetic metabolism of fetal muscle. These changes were more marked between fetuses from CN-fed dams vs. HY-fed dams, i.e., starch-based vs. fiber-based maternal diets. Several genes involved in energetic metabolism, including those implicated in glycolysis and gluconeogenesis, citrate cycle, and valine, leucine and isoleucine degradation, were down-regulated by maternal CN diet (Table 3). Similarly, a recent study has shown that a low protein prepartum diet compared against a high protein control diet when fed at isoenergetic intake in mice, caused downregulation of several pathways related to energetic metabolism processes including glycolysis, citrate cycle, and amino acid degradation, in skeletal muscle in newborn offspring [63]. Additionally, fetuses derived from HY-fed dams showed increased expression of genes involved in propanoate metabolism and insulin signaling pathway (Table 3). Overall, these findings show that prepartum maternal isoenergetic diets, differing in the primary feed source, alter the energetic metabolism of fetal muscle.

Maternal dietary treatments also perturbed different biological and molecular processes in fetal adipose tissues. For instance, in subcutaneous adipose tissue, three functional categories closely related to metabolism of phosphate anhydrides and esters showed overrepresentation of differentially expressed genes (Table 4). GO terms dephosphorylation, phosphatase activity, and phosphoric ester hydrolase activity were mainly enriched with genes that showed significant higher expression in fetuses derived from DG-fed dams. Many relevant biochemical reactions in adipose tissue, including the formation of acyl CoA during the synthesis of fatty acids, involved phosphate anhydrides and esters. Other important processes, such as ion transport in subcutaneous fat depot and cell adhesion in perirenal fat tissue, were also altered by prepartum maternal diets, showing an up-regulation of gene expression induced by a DG diet (Figure 2). It is well documented that each fat depot has its own particular features, including differences in growth, cellularity and metabolism [64]. This may explain why maternal DG diet altered mainly metabolism processes in subcutaneous fat while it affected mainly cellularity and tissue integrity in perirenal fat depot.

Potential mechanisms underlying the observed transcriptomic changes

There is growing evidence that maternal nutrition can induce epigenetic alterations to change gene expression in the fetal genome, such as DNA methylation or histone modifications including methylation, acetylation, and phosphorylation [65]. These modifications depend on the availability of key compounds, such as methyl donors, supplied by different amino acids and vitamins present in the maternal diet. Indeed, this link between maternal nutrition and subsequent modification of gene expression in the fetal genome is one of the important molecular mechanisms proposed to explain the phenomenon of fetal programming [65–68]. We have recently reported that maternal DG and HY diets were associated with altered DNA methylation and expression of the imprinted genes IGF2R and H19 in fetal sheep tissues compared to maternal CN diet [32]. Interestingly, while DNA methylation at the gene promoter suppressed expression, gene body methylation was positively associated with gene expression [32]. Here, two functional terms related to chromatin biology and epigenetics, chromatin and nucleosome assembly, were altered in fetal muscle showing significantly lower expression in fetuses from CN-fed dams. Many genes from these pathways, such as H1FX, HIST2H2BE, MBD3, and SATB1, were also differentially expressed in fetal adipose tissues, showing lower expression in fetuses from maternal CN diet. Of special interest is MBD3, an essential member of the nucleosome remodeling and deacetylation (NuRD) complex, a multi-functional complex which has key roles in chromatin remodeling during development [69]. MBD3 is associated with sites in the genome with activated genes, i.e., sites devoid of methylated CpG islands and with chromatin modifications that are characteristic of an active chromatin state [70]. Although DNA methylation or histone modifications were not directly evaluated in our experiment, it seems likely that maternal diets altered chromatin marks of key regions of DNA and these changes probably underlie the transcriptomic changes observed in fetal tissues.