Transcriptome-Wide Analysis Reveals the Role of PPARγ Controlling the Lipid Metabolism in Goat Mammary Epithelial Cells
To explore the large-scale effect of peroxisome proliferator-activated receptor (PPARG) in goat mammary epithelial cells (GMEC), an oligonucleotide microarray platform was used for transcriptome profiling in cells overexpressing PPARG and incubated with or without rosiglitazone (ROSI, a PPARγ agonist). A total of 1143 differentially expressed genes (DEG) due to treatment were detected. The Dynamic Impact Approach (DIA) analysis uncovered the most impacted and induced pathways “fatty acid elongation in mitochondria,” “glycosaminoglycan biosynthesis-keratan sulfate,” and “pentose phosphate pathway.” The data highlights the central role of PPARG in milk fatty acid metabolism via controlling fatty acid elongation, biosynthesis of unsaturated fatty acid, lipid formation, and lipid secretion; furthermore, its role related to carbohydrate metabolism promotes the production of intermediates required for milk fat synthesis. Analysis of upstream regulators indicated that PPARG participates in multiple physiological processes via controlling or cross talking with other key transcription factors such as PPARD and NR1H3 (also known as liver-X-receptor-α). This transcriptome-wide analysis represents the first attempt to better understand the biological relevance of PPARG expression in ruminant mammary cells. Overall, the data underscored the importance of PPARG in mammary lipid metabolism and transcription factor control.
Ruminant milk products are now common and popular throughout the world. Milk fat is an important component of dairy products and is a major contributor to dietary energy density. The higher concentrations of unsaturated and medium-chain fatty acids are responsible for the characteristic “goaty” odour of goat milk and also confer unique organoleptic properties . Therefore, understanding the mechanisms for altering the milk fatty acid composition of goat milk may lead to further improvements in nutritional value. Recent evidence indicates that milk fat biosynthesis is regulated by key transcription factors including peroxisome proliferator-activated receptor (PPARG) [2, 3].
It is well established that PPARG is a critical transcription factor controlling adipogenesis and glucose metabolism in various cells in nonruminants [4–6]. After binding of ligands (e.g., rosiglitazone (ROSI) or pioglitazone), PPARG causes conformational changes in the receptor [7, 8] and then forms a heterodimeric complex with RXR proteins and binds to PPAR response element (PPRE) upstream of target genes . Through controlling the downstream genes, PPARG regulates adipocyte differentiation and promotes insulin sensitivity in human and rodents . The activation of PPARG also enhances macrophage lipid uptake as well as lipid export and has anti-inflammatory effects .
In bovine cells, the activation of PPARG with rosiglitazone provided a demonstration that PPARG could control expression of genes involved in milk fat synthesis . The current data from goats indicates that PPARG regulates genes involved in triacylglycerol synthesis and secretion in mammary gland epithelial cells . It was also demonstrated that PPARG stimulates the synthesis of monounsaturated fatty acids in dairy goat mammary epithelial cells (GMEC) via the control of stearoyl-coenzyme A desaturase (SCD) . Furthermore, our recent data revealed that PPARG could modulate lipid accumulation via regulation of Perilipin 2 (PLIN2) gene expression in GMEC . Although some work [2, 3] has been performed to study the function of PPARG in ruminant mammary cells, a comprehensive dataset on gene profiles altered by PPARG is not available.
Microarray analysis provides an efficient tool to simultaneously study the expression of multiple genes in tissues or cells in response to a given treatment or physiological condition. It has been widely used in the bovine to study the differential gene expression among different treatments or physiological conditions [14–16]. Structural genomic studies of domestic animals have indicated that goats are closely related to bovine species . Previous evidences were highly suggestive that cross-species hybridization is possible using a bovine cDNA microarray to study goat gene expression [18–20].
The primary aim of this study was to assess the potential role of PPARG in GMEC at global scale. To that aim, a microarray analysis was used to detect the transcriptome alterations of GMEC after overexpression of PPARG. The results indicated that PPARG gain of function induced more than 1,000 differentially expressed genes (DEG), most of which are related to metabolism pathways.
2. Experimental Section
2.1. Cell Culture and Treatments
The mammary epithelial cells were isolated from peak lactation Xinong Saanen goats as described previously . Details of cell culture were described recently [3, 12]. Cultures of GMEC at approximately 80% confluence were transfected with one of the adenovirus supernatants (Ad-PPARG or Ad-GFP). Transfected GMEC were cultured with the PPARG-specific ligand ROSI (BioVision, USA) (PPARG+ROSI) or control [dimethyl sulfoxide (DMSO)] (Sigma, St. Louis, MO, USA) (PPARG+DMSO and Ad-GFP+DMSO) at 50 μM after 24 h of the initial culture and then harvested at 48 h (24 h later) for RNA extraction. The generation and application of the adenovirus expression PPARG (Ad-PPARG) were described elsewhere . Each treatment was performed in triplicate.
2.2. Total RNA Extraction
The procedures for total RNA extraction, purification, and qPCR were recently described . Total RNA from GMEC was extracted using the RNA Prep pure cell kit (Tiangen Biotech Co. Ltd., Beijing, China) according to the manufacturer’s protocol. The RNA used in the qPCR was treated with DNAase (Tiangen Biotech Co. Ltd., Beijing, China) to remove genomic DNA contamination. Synthesis of cDNA was conducted using the Prime Script™ RT kit (Takara Bio Inc., Otsu, Japan) according to the manufacturer’s instructions.
An Agilent platform was chosen to conduct the microarray experiment (44K Bovine (V2) gene expression microarray chip, Agilent Technologies Inc.) following the manufacturer’s protocols. Briefly, a total of 200 ng of RNA per sample were used to generate first-strand cDNA, which was reverse transcribed to cRNA using the low-input quick amp labeling kit (Agilent Technologies Inc.). The resulting cRNA was labeled with either Cy3 or Cy5 fluorescent dye, purified using RNeasy minispin columns (Qiagen), and subsequently eluted in 30 μL of DNase-RNase-free water. The NanoDrop ND-1000 (Thermo Fisher Scientific Inc., Waltham, MA) and a Bioanalyzer 2100 (Agilent Technologies) were used to confirm the manufacturer’s recommended criteria for yield of at least 0.825 μg/μL and RNA integrity , respectively.
2.4. Quantitative Real-Time PCR (qPCR)
The results from microarray were validated via qPCR for a selected panel of 15 genes considered important for fatty acid metabolism. The gene names and primers used in this study are reported in Supporting File 1 (in Supplementary Material available online at http://dx.doi.org/10.1155/2016/9195680). Methods for primer pair design and validation and qPCR were as previously described . Data of qPCR were normalized to three internal control genes, Ubiquitously Expressed, Prefoldin-Like Chaperone (UXT), Mitochondrial Ribosomal Protein L39 (MRPL39), and Ribosomal Protein S9 (RPS9).
2.5. Data Analysis
Data from microarrays were normalized using Lowess prior to statistical analysis using ANOVA in GeneSpring (Agilent Technologies). Differences in relative expression between PPARG versus CON, PPARG+ROSI versus CON, and PPARG+ROSI versus PPARG were considered significant at an unadjusted and a fold change greater or lower than 2 . The qPCR data were transformed prior to statistical analysis. The data were analyzed using a Generalized Linear Model (GLM) using SAS with treatments (CON, PPARG, and PPARG+ROSI) as the main effect. Significance was declared at .
2.6. Data Mining
Data were mined by an integrative systems biology approach applying the newly developed Dynamic Impact Approach (DIA)  and an upstream gene network analysis using Ingenuity Pathway Analysis (IPA) . The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) biological process category database of bovine were used for functional analysis with the DIA. The detailed methodology for data analysis using DIA was described previously . The IPA Knowledgebase is used to predict the expected causal effects between upstream regulators and targets (i.e., DEG).
3.1. Number of Differentially Expressed Genes (DEG) in the Microarray Data
Overall, there were more than 1,398 DEG detected by microarray. Among these, only the genes (1143) annotated with a bovine Entrez gene ID with a significant difference () and 2-fold change ratio were used for the analysis. The number of DEG indicated a marked difference in expression in the cells overexpressing PPARG with ROSI compared with cells without ROSI (Figure 1). Compared with control, there were 464 DEG upregulated and 536 DEG downregulated in PPARG+ROSI versus CON. The overexpression of PPARG alone did not markedly alter the transcriptome, but there were 72 upregulated and 22 downregulated genes. When compared with cells expressing PPARG with and without ROSI, the analysis indicated that the number of upregulated and downregulated DEG was 221 and 483, respectively.
3.2. Overall Summary of KEGG Categories
Using the DIA, the estimate of the perturbation in a biological pathway is represented by the “impact” while the overall direction of the perturbation is represented by the “flux” (or Direction of the Impact) . The DIA provides a summary of the KEGG pathways in the form of categories and subcategories (Figure 3) which are altered by treatments. The details of each pathway are reported in Supporting File S3.
In accordance with the number of DEG in Figure 1, KEGG pathway categories were more impacted in the two comparisons related to cells treated with ROSI. Among these pathways, the category “metabolism” was the most impacted (Figure 2). With the exception of the subcategories of pathways within “biosynthesis of other second metabolites,” “nucleotide metabolism,” and “amino acid metabolism,” all the other subcategories within metabolism had an impact value >25 in the comparison of PPARG+ROSI with CON.
A similar induction effect was uncovered in the comparison of PPARG+ROSI with PPARG. Except for the minor inhibition of “glycan biosynthesis and metabolism,” in the comparison of PPARG+ROSI with CON, most of the metabolic pathways were markedly activated including “carbohydrate metabolism,” “energy metabolism,” “lipid metabolism,” “amino acid metabolism,” “metabolism of other amino acids,” “glycan biosynthesis and metabolism,” “metabolism of cofactors and vitamins,” “metabolism of terpenoids and polyketides,” and “xenobiotics biodegradation and metabolism.” Compared with the control group, only the overexpression of PPARG had a weaker impact on pathway categories except “metabolism.”
According to the impact value, the categories “genetic information processing,” “environment information processing,” “cellular process,” and “organismal system” also were altered in the comparison of PPARG+ROSI with CON. However, most of their flux values were slightly activated or did not change. In PPARG+ROSI versus PPARG, the fluxes in the four categories were inhibited or exhibited no change.
3.3. Most Impacted KEGG Pathways
The DIA analysis revealed that the most impacted pathway was “fatty acid elongation in mitochondria” with flux >60, followed by “glycosaminoglycan biosynthesis” (Figure 3). The categories containing “fatty acid elongation in mitochondria,” “pentose phosphate pathway,” “glyoxylate and dicarboxylate metabolism,” “riboflavin metabolism,” “nicotinate and nicotinamide metabolism,” “PPAR signaling pathway,” and “pantothenate and CoA biosynthesis” were highly activated. In contrast, the pathways “glycosphingolipid biosynthesis-globoseries” and “folate biosynthesis” were inhibited.
Even though “glycosaminoglycan biosynthesis” was the second most impacted pathway, it was slightly inhibited by the activation of PPARG. “Glycosphingolipid biosynthesis” was highly inhibited with the activation of PPARG. Among the top ten overall most impacted terms, only PPARG belonged to “endocrine system”; the rest of them belonged to “metabolism” (Figure 3).
3.4. Expression of Selected Genes by qPCR
Fifteen genes considered important for fatty acid metabolism were selected to assess the reliability of the microarray data. Overall, >80% of genes measured by qPCR had a result deemed similar to microarray data. Compared with the control group, the cells overexpressing PPARG plus ROSI altered more genes compared with PPARG without ROSI. Among the genes involved in the upstream transcription factor regulation network, the expression level of NR1H3, PPARG, SREBF2, and PPARD by qPCR was similar to microarray, whereas data of SREBF1 and PPARGC1A were less sensitive by microarray compared with qPCR (Figure 4). A contrasting response between microarray and qPCR was also observed for FASN.
3.5. Upstream Regulators
Consistent with the number of DEG, there were a high number of upstream transcription regulators in the comparisons of PPARG+ROSI versus CON and PPARG+ROSI versus PPARG. All the upstream upregulated transcription regulators and their potential targets are depicted in Figures 5, 6, and 7. Among the upsteam transcription regulators with PPARG versus CON, there were two related to lipid metabolism including PPARG and CCAAT/enhancer binding protein (C/EBP), alpha (CEBPA). Comparing PPARG+ROSI with CON, eight upstream transcription regulators were upregulated: activating transcription factor 3 (ATF3), CEBPA, Jun protooncogene (JUN), homeobox A9 (HOXA9), hypoxia inducible factor 1, alpha (HIF1A), NR1H3, peroxisome proliferator-activated receptor delta (PPARD), and PPARG (Figure 6). Among them, CEBPA, NR1H3, PPARD, and PPARG are classical transcription factors related to lipid metabolism. Compared with PPARG+ROSI versus CON, the comparison of PPARG+ROSI with PPARG had a lower number of upregulated upstream transcription regulators including HIFA, nuclear factor, erythroid 2-like 3 (NFE2L3), NR1H3, and PPARG (Figure 7).
A few upstream transcription regulators were inhibited in the comparison of PPARG+ROSI versus CON and PPARG+ROSI versus PPARG (Figures S1 and S2). In comparison of PPARG+ROSI with CON, the transcription regulators early growth response 1 (EGR1), CXXC finger protein 1 (CXXC1), neurogenin 1 (NEUROG1), Protein Inhibitor of Activated STAT, 1 (PIAS1), pleomorphic adenoma gene-like 1 (PLAGL1), Kruppel-Like factor 4 (KLF4), RXRA, GFI1B, KLF5, KLF6, RARG, MYOD1, and SOX2 were inhibited. Similar to PPARG+ROSI versus CON, the genes expression of NEUROD1, KLF4, KLF6, CXXC1, KLF5, RXRA, myogenic differentiation 1 (MYOD1), PIAS1, sex determining region Y box 2 (SOX2), and growth factor independent 1B transcription repressor (GFI1B) was also inhibited in the comparison of PPARG+ROSI with PPARG.
Due to the unavailability of goat microarrays and the fact that structural genome of goats is closely related to that of bovine species, bovine arrays have been successfully adapted and applied in studies with goat mammary tissue [20, 25], goat ovary , and goat milk leukocytes . To further explore the transcriptome alteration by PPARG gain of function, a commercial whole-transcriptome bovine microarray was used in the present study. The data revealed close to 1,000 DEG altered by overexpression of PPARG plus the chemical agonist ROSI. The most impacted category by PPARG was related to metabolism, which agrees with the previous findings demonstrating that PPARG plays a central role in adipogenesis . Furthermore, analysis of a subset of genes by qPCR revealed a high degree of agreement with microarray data.
The DIA is efficient for the analysis of data from multiple treatment comparisons [24, 29]. Among the overall most impacted pathways in present study, the “fatty acid elongation in mitochondria,” “glycosaminoglycan biosynthesis-keratan sulfate,” and “glycosphingolipid biosynthesis-globo series” are novel and of biological interest. In adipose cells, PPARG promotes the uptake of fatty acids and storage as energy . Our previous data also revealed that PPARG stimulated the expression of genes related to triacylglycerol (TAG) synthesis in GMEC [3, 12]. Thus, we expected to find that TAG synthesis would be the most impacted pathway in the present study. The finding that “fatty acid elongation in mitochondria” was the most impacted is supported by the high expression of hydroxyacyl-CoA dehydrogenase, alpha subunit (HADHA), and hydroxyacyl-CoA dehydrogenase, beta subunit (HADHB), both of which are the rate-limiting enzymes for fatty acid elongation. Further, these genes appear to be potential PPARG target genes in GMEC. Consistent with promoting fatty acid elongation, the uptake of long-chain fatty acid was also induced because CD36  and solute carrier family 27 (fatty acid transporter), member 6 (SLC27A6), were upregulated (Figure 2).
Both qPCR and microarray revealed that the expression of long-chain acyl-CoA synthetase 1 (ACSL1) was enhanced by overexpression of PPARG with ROSI (Figure 2). ACSL1 catalyzes the conversion of free fatty acids (FFAs) into their activated acyl-CoA derivatives, which are in turn used in the cell for β-oxidation, synthesis, or reacylation of many different cellular lipids or other cellular processes. Previous data suggested that FA activation in bovine mammary tissue occurs primarily via ACSL1 due to the fact that its mRNA is the most predominant among ACSL isoforms [30, 31].
Synthesis of very-long-chain FA is carried out by fatty acid desaturases 1 (FADS1) and 2 (FADS2), which add double bonds at the and position of PUFA and synthesize eicosapentaenoic acid (20:5n-3) and docosahexaenoic acid (22:6n-3). In this study, the fact that the expression of FADS1 was significantly upregulated in PPARG-overexpressing cells indicated that this nuclear receptor may enhance the biosynthesis of polyunsaturated fatty acids. These data suggested that FADS1 may be a target of PPARG. Hence, we hypothesize that the increase of omega-3/omega-6 ratio in milk fat could be achieved through the activation of PPARG in mammary cells. In fact, the hypothesis is supported by the fact that the subcategory “biosynthesis of unsaturated fatty acid” was among the top 30 categories in this study (File S3).
The perilipin (PAT) family [32–34] and cell death-inducing DFF45 like effector (CIDE) family [35, 36] play a pivotal role in lipid formation. In the present study, the marked upregulation of PLIN2 in PPARG-overexpressing GMEC is consistent with recent data indicating that PPARG could directly bind to the promoter of PLIN2 and modulate the lipid formation in GMEC . Less is known about the role of CIDEA in lipid droplet formation in ruminant mammary cells; however, it was the only CIDE isoform which was upregulated significantly after overexpression of PPARG plus ROSI. This indicates that CIDEA is a target of PPARG.
In addition to the PPAR family, the key transcription factors SREBF1, NR1H3, CEBPA, H1F1A, JUN, and HOXA9 also had a significant change in response to the PPARG gain of function with or without ROSI. The cross talk between PPARG and SREBF1 and NR1H3 was described in our recent papers [2, 3, 12] and completely agrees with the present data that expression of SREBF1 and NR1H3 was enhanced by the overexpression of PPARG plus ROSI. The data from the IPA analysis indicating that overexpression of PPARG down- or upregulated these upstream transcription factors further supports our previous hypothesis that PPARG regulates the gene network related to fatty acid metabolism in a direct or indirect manner [3, 12]. Overall, the results indicated that goat mammary tissue relies heavily on PPARG regulation of genes to induce copious milk fat synthesis and secretion.
The pathway “glycosaminoglycan biosynthesis-keratan sulfate” is involved in the synthesis of keratan sulfate (KS); thus, its marked activation indicated that PPARG could control inflammatory response via regulating the synthesis of keratan sulfate. The hypothesis is consistent with the role of PPARG in inflammation in nonruminants [37–40]. Thus, this finding is novel and more research on KS synthesis seems warranted to better understand its role in the process inflammation, for example, during onset of mastitis.
The high activation of “pentose phosphate pathway” after overexpression of PPARG suggests that it promoted the efficient utilization of glucose in GMEC to generate substrates supporting other cellular processes. The high activation of glucose oxidation or other carbohydrate metabolism pathways in PPARG-overexpressing GMEC supports the view of a mechanism whereby PPARG alters metabolic pathways in lactating mammary gland; that is, PPARG promotes carbohydrate metabolism to produce intermediates to serve other aspects of milk fatty acid metabolism and lactose synthesis .
Due to the limitation of the microarray platform used , the interpretation of the findings from the present study has some limitations. For instance, the microarray platform used could not completely cover the genes with functional annotation in the goat genome. In addition, the difference between the goat and bovine genome will unavoidably miss some genes. In the future, goat specific oligo microarrays or next-generation sequencing should be used to confirm the transcriptome alterations caused by the PPARG gain of function. In that context, however, the present transcriptome analysis provides an initial global insight into the biological processes altered by PPARG in ruminant mammary cells.
Using cross-species hybridization microarray data, the present data support a role for PPARG activation on biological processes including and going beyond milk fat synthesis. The data indicated an overall increase in metabolism with large increase in anabolism, particularly involving fatty acid synthesis and glucose utilization. Most impacted terms underscored the regulatory role of PPARG in fatty acid elongation. The fact that pentose phosphate pathway was highly activated by PPARG suggests an important role in carbohydrate metabolism to produce intermediates for milk fatty metabolism.
The upstream regulator analysis indicated that PPARG controls molecular processes through an extensive level of cross talk with other signaling pathways, for example, JUN and CEBPA. All these data support our previous hypothesis that PPARG plays a central role in milk fatty metabolism in GMEC. In addition, the data also uncovered a likely role of PPARG in the GMEC response to inflammation via the “glycosaminoglycan biosynthesis-keratan sulfate” pathway. In conclusion, the data highlighted a strong transcriptional regulation of PPARG in the metabolism in GMEC.
The authors declare that there is no conflict of interests regarding the publication of this paper.
Hengbo Shi and Wangsheng Zhao contributed equally to this paper.
This work is jointly supported by the Transgenic New Species Breeding Program of China (2014ZX08009-051B).
Selected gene names and primers used in this study.
H. B. Shi, J. Luo, D. W. Yao et al., “Peroxisome proliferator-activated receptor-γ stimulates the synthesis of monounsaturated fatty acids in dairy goat mammary epithelial cells via the control of stearoyl-coenzyme A desaturase,” Journal of Dairy Science, vol. 96, no. 12, pp. 7844–7853, 2013.View at: Publisher Site | Google Scholar
A. K. G. Kadegowda, M. Bionaz, L. S. Piperova, R. A. Erdman, and J. J. Loor, “Peroxisome proliferator-activated receptor-γ activation and long-chain fatty acids alter lipogenic gene networks in bovine mammary epithelial cells to various extents,” Journal of Dairy Science, vol. 92, no. 9, pp. 4276–4289, 2009.View at: Publisher Site | Google Scholar
K. Shahzad, M. Bionaz, E. Trevisi, G. Bertoni, S. L. Rodriguez-Zas, and J. J. Loor, “Integrative analyses of hepatic differentially expressed genes and blood biomarkers during the peripartal period between dairy cows overfed or restricted-fed energy prepartum,” PLoS ONE, vol. 9, no. 6, Article ID e99757, 2014.View at: Publisher Site | Google Scholar
L. Schibler, D. Vaiman, A. Oustry, C. Giraud-Delville, and E. P. Cribiu, “Comparative gene mapping: a fine-scale survey of chromosome rearrangements between ruminants and humans,” Genome Research, vol. 8, no. 9, pp. 901–915, 1998.View at: Google Scholar
S. Ollier, C. Robert-Granié, L. Bernard, Y. Chilliard, and C. Leroux, “Mammary transcriptome analysis of food-deprived lactating goats highlights genes involved in milk secretion and programmed cell death,” The Journal of Nutrition, vol. 137, no. 3, pp. 560–567, 2007.View at: Google Scholar
S. Ollier, C. Leroux, A. de la Foye, L. Bernard, J. Rouel, and Y. Chilliard, “Whole intact rapeseeds or sunflower oil in high-forage or high-concentrate diets affects milk yield, milk composition, and mammary gene expression profile in goats,” Journal of Dairy Science, vol. 92, no. 11, pp. 5544–5560, 2009.View at: Publisher Site | Google Scholar
Z. Wang, J. Luo, W. Wang, W. Zhao, and X. Lin, “Characterization and culture of isolated primary dairy goat mammary gland epithelial cells,” Chinese Journal of Biotechnology, vol. 26, no. 8, pp. 1123–1127, 2010.View at: Google Scholar
X.-Z. Lin, J. Luo, L.-P. Zhang, W. Wang, H.-B. Shi, and J.-J. Zhu, “MiR-27a suppresses triglyceride accumulation and affects gene mRNA expression associated with fat metabolism in dairy goat mammary gland epithelial cells,” Gene, vol. 521, no. 1, pp. 15–23, 2013.View at: Publisher Site | Google Scholar
M. Bionaz, K. Periasamy, S. L. Rodriguez-Zas, W. L. Hurley, and J. J. Loor, “A novel dynamic impact approach (DIA) for functional analysis of time-course omics studies: validation using the bovine mammary transcriptome,” PLoS ONE, vol. 7, no. 3, Article ID e32455, 2012.View at: Publisher Site | Google Scholar
F. Wilfling, J. T. Haas, T. C. Walther, and R. V. Farese Jr., “Lipid droplet biogenesis,” Current Opinion in Cell Biology, vol. 29, pp. 39–45, 2014.View at: Google Scholar