Email updates

Keep up to date with the latest news and content from BMC Genomics and BioMed Central.

Open Access Research article

Functional and gene network analyses of transcriptional signatures characterizing pre-weaned bovine mammary parenchyma or fat pad uncovered novel inter-tissue signaling networks during development

Paola Piantoni1, Massimo Bionaz1, Daniel E Graugnard1, Kristy M Daniels2, Robin E Everts14, Sandra L Rodriguez-Zas1, Harris A Lewin13, Hurley L Hurley1, Michael Akers5 and Juan J Loor1*

Author Affiliations

1 Department of Animal Sciences, University of Illinois, 1207 West Gregory Drive, Urbana, IL, 61801, USA

2 Department of Animal Science, The Ohio State University, Wooster, OH, 44691, USA

3 Institute for Genomic Biology, University of Illinois, 1206 West Gregory Drive, Urbana, 61801, USA

4 Sequenom, Inc., 3595 John Hopkins Court, San Diego, CA 92121, USA

5 Dairy Science Department, Virginia Tech, Blacksburg, VA, 24061, USA

For all author emails, please log on.

BMC Genomics 2010, 11:331  doi:10.1186/1471-2164-11-331


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2164/11/331


Received:6 November 2009
Accepted:26 May 2010
Published:26 May 2010

© 2010 Piantoni et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

The neonatal bovine mammary fat pad (MFP) surrounding the mammary parenchyma (PAR) is thought to exert proliferative effects on the PAR through secretion of local modulators of growth induced by systemic hormones. We used bioinformatics to characterize transcriptomics differences between PAR and MFP from ~65 d old Holstein heifers. Data were mined to uncover potential crosstalk through the analyses of signaling molecules preferentially expressed in one tissue relative to the other.

Results

Over 9,000 differentially expressed genes (DEG; False discovery rate ≤ 0.05) were found of which 1,478 had a ≥1.5-fold difference between PAR and MFP. Within the DEG highly-expressed in PAR vs. MFP (n = 736) we noted significant enrichment of functions related to cell cycle, structural organization, signaling, and DNA/RNA metabolism. Only actin cytoskeletal signaling was significant among canonical pathways. DEG more highly-expressed in MFP vs. PAR (n = 742) belong to lipid metabolism, signaling, cell movement, and immune-related functions. Canonical pathways associated with metabolism and signaling, particularly immune- and metabolism-related were significantly-enriched. Network analysis uncovered a central role of MYC, TP53, and CTNNB1 in controlling expression of DEG highly-expressed in PAR vs. MFP. Similar analysis suggested a central role for PPARG, KLF2, EGR2, and EPAS1 in regulating expression of more highly-expressed DEG in MFP vs. PAR. Gene network analyses revealed putative inter-tissue crosstalk between cytokines and growth factors preferentially expressed in one tissue (e.g., ANGPTL1, SPP1, IL1B in PAR vs. MFP; ADIPOQ, IL13, FGF2, LEP in MFP vs. PAR) with DEG preferentially expressed in the other tissue, particularly transcription factors or pathways (e.g., MYC, TP53, and actin cytoskeletal signaling in PAR vs. MFP; PPARG and LXR/RXR Signaling in MFP vs. PAR).

Conclusions

Functional analyses underscored a reciprocal influence in determining the biological features of MFP and PAR during neonatal development. This was exemplified by the potential effect that the signaling molecules (cytokines, growth factors) released preferentially (i.e., more highly-expressed) by PAR or MFP could have on molecular functions or signaling pathways enriched in the MFP or PAR. These bidirectional interactions might be required to coordinate mammary tissue development under normal circumstances or in response to nutrition.

Background

As reported by Connor and colleagues [1]: "The mammary gland is a complex organ of various tissue and cell types that will undergo multiple stages of growth, differentiation, secretory activity, and involution during the lifetime of a female mammal". Among the "various tissues" the parenchyma (PAR), which is, in lactating mammary gland, the tissue that synthesizes and secretes milk, and the fat pad (MFP), which is a matrix of connective and adipose tissue surrounding the PAR [2], are considered the most crucial during post-natal development.

Interactions between PAR and MFP during bovine mammary development are still not fully understood. It has been postulated that during mammary development the MFP surrounding PAR exerts proliferative effects on the PAR through secretion of local modulators of growth induced by the impacts of selected systemic hormones (e.g., growth hormone, estrogen) [1,3-5] or growth factors (e.g., IGF-1) [6,7]. It is believed that such an effect occurs because the epithelial tissue that is in direct contact with the MFP has a greater degree of proliferation compared with the more central epithelial tissue [8-10].

Local interaction between PAR and MFP could occur in both directions, i.e., MFP acts on PAR and PAR acts on surrounding MFP [11]. How these tissues could communicate through locally-produced modulators has not yet been studied in the pre-weaning prepubertal bovine mammary gland. Hovey and colleagues [12], using prepubertal ewes, showed that IGF1 mRNA expression was greater in MFP cells adjacent to PAR than in MFP cells with no PAR contact, which indicated the existence of a local "diffusible factor" secreted by PAR that could increase the expression of IGF isoform in MFP. Based on those findings, a potential crosstalk between the two tissues was suggested. It was proposed that MFP stimulates PAR and PAR then exerts a positive feedback on the MFP during development [2].

Mammary gland development and tissue interactions have been previously studied using gene expression analysis. For example, in a serial slaughter study [13] it was observed that peak expression of IGF1 in MFP and estrogen receptor-α (ESR1) in PAR from 100 kg body weight Holstein heifers coincided with peak mammary epithelial cell proliferation [14]. Expression of both genes decreased in mammary tissue in older animals. Regarding tissue interaction, Thorn and colleagues [15] hypothesized that the MFP could impact PAR through inflammation-related proteins, such as TNFα, IL-6, and IL-1β. They showed in vitro the inhibitory effect on proliferation of TNFα, but not IL-6 or IL-1β, on epithelial cell proliferation.

Li and colleagues [16] conducted a microarray study exploring the interaction between MFP and PAR in response to estrogen treatment in prepubertal heifers. Results indicated that MFP might affect PAR cell proliferation via the secretion of paracrine stimulators such as the stem cell growth factor precursor C-type lectin domain family 11 member A (CLEC11A) and IGF-1. Despite work conducted to date, there is still uncertainty regarding how PAR and MFP tissues interact during mammary development in prepubertal heifers prior to weaning.

In the present study, mammary glands of pre-weaned Holstein heifer calves were harvested at 65 d of age to extract total RNA for microarray analysis. Extensive bioinformatics analysis of microarray data was performed to (1) characterize differences in transcript profiles between mammary PAR and MFP, with the specific aim of uncovering predominant transcriptomic signatures, and (2) uncover predominant signaling molecules (e.g., cytokines and growth factors) in one tissue relative to the other. The latter would allow for the identification of potential targets among genes that are more highly expressed in the other tissue in order to discover novel inter-tissue signaling networks.

Results

Coverage of microarray elements in the IPA knowledge base

Over 10,000 oligonucleotides (ca. 76% of total) from the microarray (see details in Additional file 1) were mapped by Ingenuity Pathway Analysis® (IPA). Of these, > 7,500 genes were eligible for generating networks and >6,400 genes were associated with a function or pathway. Almost 90% of all annotated genes in our microarray were differentially expressed between PAR and MFP (Table 1 and Additional file 2). Of these, 16.3% had a difference between tissues of > 1.5-fold, with 8.1% being more highly expressed in PAR and 8.2% being more highly expressed in MFP (Table 1 and Additional file 2). Among DEG, 0.6% and 0.8% were 3-fold greater in PAR vs. MFP and MFP vs. PAR, respectively. The data mining analysis was performed on DEG exhibiting ≥1.5-fold difference between tissues. Among DEG, ca. 500 were eligible for generating networks in either PAR or MFP in IPA (Table 1).

Additional file 1. Methods and Figure S1. The file contains additional materials and methods as well as additional results. Extraction and purification of RNA, cDNA synthesis, microarrays, qPCR, Primer design, and GO data analysis accompanied by 4 tables: Table S1 titled 'Accession number, sequence, and amplicon size of genes used for qPCR', Table S2 titled 'Sequencing results of PCR products of genes designed for this study', Table S3 'Genes chosen on fold differences or highly enriched function by DEG between tissues to verify microarray data', and Table S4 titled 'Genes chosen on fold differences or highly enriched function by DEG between tissues to verify microarray data'.

Format: DOC Size: 204KB Download file

This file can be viewed with: Microsoft Word ViewerOpen Data

Additional file 2. Statistical analysis of genes on the microarray. The file contains statistical P-values for each gene on the microarray including Benjamini-Hochberg FDR P-values. In addition, separate sheets are included that contain the DEG and DEG sorted by fold-change both in PAR and MFP: 1,478 DEG with FDR ≤ 0.05 and fold difference (FD) >1.5; 736 DEG higher in parenchyma vs. fat pad with FDR ≤ 0.05 and FD >1.5; 742 DEG higher in fat pad vs. parenchyma with FDR ≤ 0.05 and FD >1.5; 59 DEG higher in parenchyma vs. fat pad with FDR ≤ 0.05 and FD >3; 75 DEG higher in fat pad vs. parenchyma with FDR ≤ 0.05 and FD > 3.0.

Format: XLS Size: 5.7MB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Table 1. Number of overall differentially expressed genes (DEG, Benjamini-Hochberg FDR ≤ 0.05) and DEG with cut-off of 1.5- and 3.0-fold difference in mRNA abundance between mammary parenchyma and fat pad from pre-weaned Holstein calves (ca. 65 d of age).

Functional analysis of DEG between PAR and MFP

The IPA analysis results using all DEG with ≥1.5-fold between PAR and MFP are reported in detail in Additional file 3. Briefly, cell movement, cell death, cell growth and proliferation, cell-to-cell signaling and interaction, and tissue development were the top 5 functions among DEG with ≥1.5-fold expression difference. Among canonical pathways, the top 9 were Aryl hydrocarbon receptor signaling, metabolism of xenobiotics by cytochrome P450, propanoate metabolism, pyruvate metabolism, LPS/IL-1 inhibition of RXR function, xenobiotic metabolism signaling, α-adrenergic signaling, p53 signaling, and acute phase response signaling. Interestingly, all of the mentioned pathways were primarily, if not completely, enhanced by genes that were more highly-expressed in MFP vs. PAR. The exception was p53 signaling, which was mostly enriched by genes more highly-expressed in PAR vs. MFP (Additional file 3)

Additional file 3. Functions and canonical pathways within DEG. Functions, canonical pathways, effect on function, and transcription factors within DEG (n = ca. 9,000) in both PAR and MFP using IPA (Ingenuity Systems, Inc.). IPA analysis DEG 1.5-fold, without separation between the two tissues and report images of overall functions and pathways affected by DEG with relative -logP-value and list of genes for each significant affected function/pathway: effect on function DEG 1.5-fold, without separation between the two tissues; pathways in all DEG 1.5-fold, reports images of significant affected pathways IPA analysis in DEG PAR, genes highly expressed in PAR vs. MFP and report images of overall functions and pathways affected by DEG with relative -logP-value and list of genes for each significant affected function/pathway; effect on function PAR, genes highly expressed in PAR vs. MFP; pathways in DEG PAR, reports images of significantly affected pathways; IPA analysis in DEG MFP, genes highly expressed in MFP vs. PAR and images of overall functions and pathways affected by DEG with relative -logP-value and list of genes for each significant affected function/pathway; effect on function MFP, genes highly expressed in MFP vs. PAR; pathways in DEG MFP, reports images of significant affected pathways; transcription factor PAR genes highly expressed in PAR vs. MFP; transcription factor MFP, genes highly expressed in MFP vs. PAR.

Format: XLS Size: 12.6MB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Results of the most-enriched biological processes from the Gene Ontology (GO) analysis that considered all DEG with ≥ 1.5-fold between the two tissues are reported in Figure 1 and in more detail in Additiona file 1. Cell signaling and development were among the most enriched biological processes. Among signaling-associated GO biological process categories (Figure 1) more highly expressed among DEG in PAR relative to MFP, those associated with the protein kinase cascade and cell activation were most predominant. The same types of molecules were enriched in most of development-associated categories. Other functions enriched among DEG with ≥ 1.5-fold in PAR vs. MFP were associated with cell death, cell organization and biogenesis, metabolism and homeostasis, and localization and transport (prevalently protein transport). Biological process categories enriched among DEG ≥ 1.5-fold in MFP vs. PAR were related to wound healing, catabolic processes, and regulation of localization and transport.

thumbnailFigure 1. GO analysis of DEG. Biological processes significantly-enriched with a Benjamini-Hochberg corrected-P-value ≤ 0.05 among all differentially expressed genes (DEG) with ≥ 1.5-fold expression between parenchyma (PAR) and mammary fat pad (MFP). Reported are the Biological processes clustered in pre-selected categories (left column) by the authors to simplify interpretation of the data, the number of DEG for each process (middle column), and main tissue affected by the function (right column): dark blue denotes high enrichment in PAR, orange denotes high enrichment in MFP, and light green denotes same magnitude of enrichment between tissues (see Additional Materials and Methods for explanation). Additional Biological processes and other GO categories (including also Cellular component and Molecular function) are available in Additional file 4.

Functions overrepresented in DEG ≥ 1.5-fold in PAR vs. MFP

The main results from the functional analysis with IPA are reported in Table 2. Complete details of the analysis and associated genes are reported in Additional file 3. Among 21 significantly-enriched functions (Benjamini-Hochberg FDR-corrected P ≤ 0.01), most pertained to functions related to cell development and structure, which became more evident when compared to the functional analysis of DEG more highly expressed in MFP vs. PAR (Table 3). For example, the most-enriched functions among DEG more highly expressed in PAR vs. MFP (with > 130 genes) were cell death, cell growth and proliferation, cellular development, cellular movement, and cell morphology. Enrichment of these opposing processes likely reflects the different cell types within PAR to accommodate the remodeling. Detailed functional analysis of DEG that were more highly-expressed in PAR vs. MFP suggested a greater degree of apoptosis, proliferation/growth/development, movement and adhesion of cells, and morphogenesis/shaping of cells in PAR vs. MFP (Table 2). Overall, angiogenesis, DNA metabolism, and survival of mammals functions were enriched (Table 2) in PAR when compared to MFP. However, detailed analysis of DEG did not indicate induction of gene expression associated with these specific functions.

Table 2. Significantly (Benjamini-Hochberg FDR ≤ 0.01) enriched functions among differentially expressed genes (DEG) highly-expressed in parenchyma relative to fat pad using Ingenuity Pathways Analysis® (IPA).

Table 3. Significantly (Benjamini-Hochberg FDR ≤ 0.01) enriched functions among differentially expressed genes (DEG) highly-expressed in fat pad relative to parenchyma tissue using Ingenuity Pathways Analysis® (IPA).

The most-enriched biological processes from GO analysis within DEG with ≥ 1.5-fold higher expression in PAR vs. MFP were associated with cell organization (chiefly chromosome/chromatin remodeling), development, differentiation, cell cycle, negative regulation of nucleotide metabolism, negative regulation of transcription, response to DNA damage, and antigen presentation and inhibition of immune system process (Figure 2). Among molecular functions, the most enriched related to binding, particularly protein (i.e., cytoskeletal and enzyme) and nucleotide binding (Additional file 4). Also, functions related to transcription, and particularly repression of transcription, zinc binding, kinases, and cytokine activity were significantly enriched (Additional file 4). The most-enriched cellular components in PAR were of cytosolic origin, particularly components associated with cytoskeleton and intracellular non-membrane-bound organelles. Components related to the nucleus, particularly chromosome allocation were also noted (Additional file 4).

thumbnailFigure 2. Comparison of GO enrichment between tissues. Biological processes significantly-enriched with a Benjamini-Hochberg corrected-P-value ≤ 0.05 and uniquely present in differentially expressed genes (DEG) with ≥ 1.5-fold expression in parenchyma (PAR) or mammary fat pad (MFP) relative to the other tissue. Reported are the Biological processes (center column) clustered in pre-selected categories by the authors to simplify interpretation of the data, and the number of DEG for each process (right column). Additional Biological processes and other GO categories (including also Cellular component and Molecular function) are available in Additional file 4 for the same analysis.

Additional file 4. GO analysis. This file contains results from GO analysis of analysis of DEG (overall genes with significant ≥ 1.5-fold difference in expression between PAR and MFP), DEG higher in PAR vs. MFP, and DEG higher in MFP vs. PAR. DEG are classified according to significant biological process, molecular function, and cellular components. Analysis of GO was conducted using DAVID http://david.abcc.ncifcrf.gov:8080/ webcite.

Format: XLS Size: 155KB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Functions overrepresented in DEG with ≥ 1.5-fold expression in MFP vs. PAR

The most-abundant functions among DEG that were more highly-expressed in MFP vs. PAR are reported in Table 3. Complete details of the functional analysis and associated genes are reported in the Additional file 3. Metabolic functions related to transport and lipid synthesis and oxidation were among the most significantly-enriched. In addition, IPA analysis of genes more highly expressed in MFP vs. PAR suggested that MFP compared with PAR had a greater degree of cell migration, carbohydrate metabolism, activation of cells, and a lower cellular growth and proliferation (Table 3).

The GO analysis indicated that DEG more highly-expressed in MFP vs. PAR enriched significantly metabolism, particularly lipid biosynthesis, catabolism, and oxidation (Figure 2). A more important role of signaling in MFP compared with PAR was suggested by the significant enrichment of signaling response-related genes, particularly for response to insulin and cytokines (Figure 2). Among cellular functions, the overall analysis indicated a large enrichment of enzyme activity and transport, with the former related to phosphate metabolism and oxidoreductase activity and the latter related to ion (particularly Mg, Fe, and K) and carbohydrate transport with a large enrichment of passive transport-related molecules (Additional file 4). The only significantly-enriched cellular component in MFP was the mitochondria and its related membranes (Additional file 4).

Main canonical pathways overrepresented among DEG with ≥ 1.5-fold between PAR and MFP

In PAR, the main canonical pathway overrepresented was actin cytoskeleton signaling. This pathway was significantly enriched with an FDR ≤ 0.05 and appeared that it was enhanced among DEG with > 1.5-fold expression in PAR vs. MFP (Table 4). Other pathways appeared to be enriched (Table 4) at a lower level of significance (i.e., FDR ≤ 0.12). Most of those are related to signaling (actin, p53, Wnt/β-catenin, PI3K/AKT) and cell cycle (Table 4). Detailed visualization of the pathways (Additional file 3) suggested an autocrine effect on actin cytoskeleton signaling elicited by FGF (fibroblast growth factor), PDGF (platelet-derived growth factor), and FN1 (fibronectin 1), all molecules that were more highly-expressed in PAR vs. MFP (Additional file 2). Among the canonical pathways enriched with a non-corrected P-value ≤ 0.01, a detailed analysis of Wnt/β-catenin and PI3K/AKT signaling revealed a potentially crucial role for these in increasing cell cycle activity (e.g., mitosis, Table 2) and apoptosis as well as increasing protein synthesis through mTOR (Additional file 3).

Table 4. Top canonical metabolic and signaling (in bold font) pathways uncovered by Ingenuity Pathway Analysis® (Benjamini-Hochberg FDR ≤ 0.05) within differentially expressed genes (DEG) with > 1.5-fold mRNA abundance between parenchyma and fat pad in mammary gland from pre-weaned Holstein heifers.

The DEG more highly-expressed in MFP vs. PAR had a significant (FDR ≤ 0.05) enrichment of 23 canonical pathways with a large presence of metabolic pathways (11 out of 23 enriched pathways; Table 4). Most of the significantly-enriched canonical pathways were involved in energy utilization, especially utilization of glucose, fatty acids, and several amino acids as sources of energy. Among signaling pathways, most pertained to immune response (e.g., LPS/IL-1 mediated inhibition of RXR function, acute phase response signaling, Complement system) and stress/catabolic response (e.g., mitochondrial dysfunction, NRF2-mediated oxidative stress response, β-adrenergic signaling, cAMP-mediated signaling, aryl hydrocarbon receptor signaling [AHR], xenobiotic metabolism signaling). The overall analysis of canonical pathways suggested an increase in oxidation of organic compounds as also suggested by functional IPA and GO analyses (Table 3 and Figures 1 and 2). The canonical pathway analysis also indicated an inherently-greater predisposition of MFP to mount an immune response compared with PAR.

One noteworthy canonical pathway highlighted by the analysis was IGF-1 signaling, which was enriched at an uncorrected P-value ≤ 0.01 (Table 4). The lack of significant enrichment potentially suggests a lower importance of this pathway, but due to the known importance of IGF-1 in mammary development and the likely existence of crosstalk between PAR and MFP [6,13,17,18] we have reported details related to this pathway. Detailed analysis revealed an enrichment of several genes that are downstream (insulin receptor substrate 2 [IRS2], v-akt murine thymoma viral oncogene homolog 1 and 2 [AKT1 and 2], protein kinases, Ras) effectors of IGF-1 and its signaling, which overall seems to suggest an increase in cell growth and survival (Additional file 3). However, a number of upstream effectors (e.g., IGFBP1, 3, 5, 6) also were enriched in MFP. The above agree with results from Daniels et al. [10] using the same tissues as in the present study. Due to the indication by the functional analysis (Table 3) of an apparent decrease in overall cell proliferation, it could be possible that the IGFBPs exerted some level of control on IGF-1 availability to MFP in these animals

Transcription factors potentially controlling DEG with ≥ 1.5-fold expression between PAR and MFP

Long-term transcriptomic adaptations of tissues are driven by transcription factors (TF), which sense external stimuli and allow cellular functions to adapt to the specific stimulus. Among DEG more highly-expressed in PAR vs. MFP, IPA identified 76 transcriptional regulators and 2 ligand-dependent nuclear receptors (ESR1 and nuclear receptor subfamily 2 group F member 2 [NR2F2], Additional file 3). Twenty-nine TF are potentially able to affect the expression of the other 126 DEG more highly expressed in PAR vs. MFP, based on the IPA knowledge base (Figure 3). Most of the genes in the transcriptional networks (Figure 3) are involved in cell cycle and proliferation (e.g., DEG affected by v-myc myelocytomatosis viral oncogene homolog (avian) [MYC] and TP53). Interestingly, several genes encoding cytokines and growth factors (e.g., IL7, SPP1, CXCL10; see description in Table 5) were present in the transcriptional networks, particularly as potential targets of MYC, TP53, and catenin beta 1 (CTNNB1) (Figure 3).

thumbnailFigure 3. Regulatory network expression in mammary parenchyma. Network analysis using differentially expressed genes (DEG) largely expressed (≥ 1.5-fold) in parenchyma (PAR) vs. mammary fat pad (MFP). The interactions shown involve the effect on expression (E), transcription (T), and protein-DNA interactions (PD) between transcription factors (blue background, white font) and other DEG largely expressed (≥ 1.5-fold) in PAR vs. MFP. Molecule types have their background highlighted as follows: enzymes, yellow; cytokines, red; growth factors, dark-green; membrane receptors, black (white font); kinases and phosphatases, gray; transporters, light-blue. The legend for the shape of the objects is reported in the figure. The intensity of the color in the object is proportional to the fold-difference in PAR vs. MFP.

Table 5. Parenchyma and fat pad cytokines and growth factors among DEG between tissues (1.5-fold difference in mRNA abundance).

Among DEG highly-expressed in MFP vs. PAR, we uncovered 40 TF and 5 ligand-dependent nuclear receptors (peroxisome proliferator-activated receptor γ [PPARG], nuclear receptor subfamily 1 group D member 1 [NR1D1], progesterone receptor membrane component 2 [PGRMC2], and retinoid × receptor beta [RXRB] and gamma [RXRG]; Additional file 3). RXRG (7.6-fold) and PPARG (2.8-fold) had the largest differences in expression in MFP vs. PAR (Additional file 3). IPA transcriptional networks indicated that 14 of those TF could potentially affect the expression of 77 more highly expressed genes in MFP vs. PAR (Figure 4). PPARG, hepatic nuclear factor 4 alpha [HNF4A], Kruppel-like factor 2 (lung) [KLF2], early growth response 2 [EGR2], and v-fos FBJ murine osteosarcoma viral oncogene homolog [FOS] had potentially larger effects on transcription of genes with greater expression in MFP vs. PAR (Figure 4). PPARG controls the expression of several enzymes involved in triacylglycerol synthesis and also expression of adipokines, such as adiponectin [ADIPOQ] and leptin [LEP]. The major functions of the transcriptional networks among more highly expressed DEG in MFP vs. PAR were lipid metabolism (synthesis and transport, mainly) and cellular movement.

thumbnailFigure 4. Regulatory network expression in mammary fat pad. Network analysis using differentially expressed genes (DEG) largely expressed (≥ 1.5-fold) in mammary fat pad (MFP) vs. parenchyma (PAR). The interactions shown involve the effect on expression (E), transcription (T), and protein-DNA interactions (PD) between transcription factors (blue background, white font) and other DEG largely expressed (≥ 1.5-fold) in MFP vs. PAR. Molecule types have their background highlighted as follows: enzymes, yellow; cytokines, red; growth factors, dark-green; membrane receptor, orange. Highlighted in dark-blue is THRSP, which is considered a key lipogenic transcription factor but IPA did not recognize it as such. The legend for the shape of the objects is reported in the figure. The intensity of the color in the object is proportional to the fold-difference in MFP vs. PAR.

Cytokines and growth factors among DEG with ≥ 1.5-fold expression between PAR and MFP

A total of 9 genes classified as growth factors and 10 genes classified as cytokines had ≥ 1.5-fold greater expression in PAR vs. MFP (Table 5). Among these were osteopontin [SPP1] and chemokine (C-X-C motif) ligand 10 [CXCL10], which were 22- and 3.5-fold greater in PAR vs. MFP (Table 5). There were 8 genes classified as growth factors and 7 classified as cytokines with greater expression in MFP vs. PAR (Table 5). Among these were ADIPOQ and FGF2, which were 42- and 3.5-fold greater in MFP vs. PAR.

Gene network analysis revealed that most of the signaling molecules identified can potentially elicit effects on 1.5-fold DEG between PAR and MFP (Table 5; Figure 5 and 6). Among these is interleukin 1 beta [IL1B], which had greater expression in PAR vs. MFP in qPCR analysis (Table S4 in Additional file 1), and had the largest number of connections with DEG that were highly-expressed in MFP vs. PAR (Figure 5). Among the various actions of IL1B, it can affect the expression of several TF, some of which were highly expressed in MFP vs. PAR: PPARG and THRSP expression, crucial for lipid synthesis, is inhibited by IL1B, whereas endothelial PAS domain protein 1 [EPAS1] and NR1D1 expression is increased by the same cytokine (Figure 5). IL1B also appears to control the expression of several cytokines and growth factors potentially more actively secreted by MFP compared with PAR (i.e., ≥ 1.5-fold more expressed in MFP vs. PAR), including inhibition of LEP and FGF2, and activation of IL1A (Figure 5).

thumbnailFigure 5. Cytokine and growth factor enrichment in mammary parenchyma. Interaction of cytokines and growth factors highly expressed in parenchyma (PAR) relative to mammary fat pad (MFP) with DEG highly expressed in MFP vs. PAR. Objects for cytokines and growth factors largely expressed in PAR vs. MFP are denoted by blue-filled objects. Objects for DEG largely expressed in MFP vs. PAR are depicted in shades of orange. The intensity of the color relates to the fold difference in PAR vs. MFP (blues) of MFP vs. PAR (orange). Molecule types have their background highlighted as follows: enzymes, yellow; cytokines and growth factors potentially secreted by PAR, red; cytokines and growth factors highly-expressed in MFP vs. PAR, green; membrane and G-coupled receptors, black (white font); phosphatase and kinases, gray; transcription factors and nuclear-dependent transcription regulators, dark-pink. All other molecules have white background. The intensity of the shaded color in the object relates to the fold-difference in MFP vs. PAR (orange) or PAR vs. MFP (blue). Red arrows denote effects on gene expression, purple arrows denote activation, dark-violet arrows denote phosphorylation, and green arrows denote protein-protein interactions. Arrow edges are expression (E), activation (A), modification (m), protein-protein interaction (PP), protein-RNA interaction (PR), phosphorylation (P), and translocation (TR). Legend for the shape of the objects is reported in the figure.

thumbnailFigure 6. Cytokine and growth factor enrichment in mammary fat pad. Interaction of cytokines and growth factors highly expressed in mammary fat pad (MFP) relative to parenchyma (PAR) with differentially expressed genes (DEG) highly expressed in PAR vs. MFP. Objects for cytokines and growth factors largely expressed in MFP vs. PAR have an orange shade. DEG largely expressed in PAR vs. MFP are depicted in blue shades. The intensity of the shaded color in the object relates to the fold-difference in MFP vs. PAR (orange) or PAR vs. MFP (blue). Molecule types have their background highlighted as follows: enzymes, yellow; cytokines and growth factors potentially secreted by MFP, red; cytokines and growth factors highly expressed in PAR vs. MFP light-green; membrane and G-coupled receptors, black (white font); phosphatase and kinases, gray; transcription factors and nuclear-dependent transcription regulators, dark-pink. All other molecules have white background. Red arrows denote effects on gene expression, purple arrows denote activation, dark-violet arrows denote phosphorylation, and green arrows denote protein-protein interaction. Arrow edges are expression (E), activation (A), modification (m), protein-protein interaction (PP), protein-RNA interaction (PR), phosphorylation (P), and translocation (TR). Legend for the shape of the objects is reported in the figure.

Other signaling molecules that are likely secreted in greater amounts by PAR than MFP include FGF7 and neuregulin 1 [NRG1]. Based on IPA annotations, FGF7 decreases the expression of stearoyl-CoA desaturase [SCD] while it increases the expression of fatty acid synthase [FASN]; whereas, NRG1 appears to regulate the expression of several TF such as FOS and EGR2 (Figure 5). These two molecules could potentially control the expression of several genes that were more highly-expressed in MFP vs. PAR (Figure 4). In addition, IL7 and chemokine (C-C motif) ligand 2 [CCL2], potentially released by PAR, may have determined the greater expression of IL13 in MFP vs. PAR (Additional file 2).

Among cytokines and growth factors potentially released to a greater extent from MFP compared with PAR, IL13, FGF2, IL1A, and LEP had the largest number of potential interactions with DEG that had greater expression in PAR vs. MFP (Figure 6). FGF2 may affect a number of biological events through increasing the expression of IL1B and PDGFA, decreasing the expression of CCL2, and decreasing the activation of TP53 in PAR. The transcription factor TP53 might control the expression of many genes in PAR (Figure 3). IL13 also might control the expression of several of the same cytokines and growth factors controlled by FGF2 (Figure 6). LEP might increase expression of several enzymes involved in metabolism (e.g., aconitase 1 [ACO1], carnitine O-octanoyltransferase [CROT]) as well as decrease the expression of ESR1, which seems to have a central role in regulating expression of several other transcription factors (e.g., CTNN1B, FOXA1, MYC) as well as the prolactin receptor (PRLR) in PAR (Figure 3). Lastly, it is noteworthy to highlight the cytokine IL1A for its effects on increasing expression of MYC, one of the most studied cytokines with a demonstrated role in many functions, chiefly growth and development [19](Figure 3).

Integrative model of potential interactions between MFP and PAR

Results from the functional and gene networks analyses of microarray data were used to develop an integrative model of putative interactions between MFP and PAR (Figure 7). The model highlights growth factors and cytokines that seem to be preferentially expressed in one tissue versus the other. Based on our analysis, which relied on data within the IPA knowledge base, it appears that many of these molecules could interact with genes preferentially expressed in the other tissue and affect a wide range of molecular and cellular functions (Figure 7).

thumbnailFigure 7. Model of inter-tissue crosstalk. Summary of potential effects of cytokines and growth factors predominantly secreted by parenchyma (PAR) or mammary fat pad (MFP) tissue on molecular and cellular functions and signaling pathways suggested by differentially expressed genes (DEG) of the MFP or PAR. The model is based on results from the current study and functions and pathways reported are relative to the functional analysis performed by Ingenuity Pathway Analysis® of interactive networks shown in Figure 5 and 6. See discussion for description of the model.

Microarray data verification through qPCR

Twenty-four out of 25 genes selected were verified via qPCR (Tables S3 and S4 in Additional file 1). These genes were chosen based on their level of expression between tissues as well as their associated biological function. In contrast with microarray data, IL1B had greater expression in PAR vs. MFP which was probably due to presence of unspecific binding of the labeled cDNA probe to the oligo on the microarray. Genes confirmed to have greater expression in PAR vs. MFP were LTF, TNC, SPP1, PRLR, CSN3, ESR1, A2M, ACTB, CDH1, MYC, CTNNB1, and TP53 (Additional file 1, Table S3 and S4).

Discussion

General considerations

The potential crosstalk between MFP and PAR during bovine mammary development has been recognized for some time (e.g., [2,11,12]). In addition to ovarian steroids [1] and other mammogenic hormones [4], it is believed that various proteins secreted by the MFP (e.g., cytokines, growth factors) act as local regulators of mammary gland growth and morphogenesis [2,11]. In this regard, greater expression of IGF1 in MFP and ESR1 in PAR tissue has been associated with greater rates of cell proliferation in mammary gland of early prepubertal animals (body weight = 100 kg) than older animals (body weight >100 kg) [13]. In support of this, regardless of nutritional management during the pre-weaning period, Ellis and Capuco [20] demonstrated that there was more epithelial cell proliferation at 2 mo than at 5 or 8 mo age. Those data indicated that mammary gland development and growth is substantial during the pre-weaning period.

Despite the large volume of research regarding development of mammary gland, the functional characterization of mammary PAR and MFP including the underlying gene networks and pathways, as well as putative interactions between both tissues during early development remain to be fully elucidated. The use of transcriptomics analysis via microarrays in combination with bioinformatics analysis can lead to an explosion of biological findings and help uncover potential crosstalk between the PAR and MFP which might affect long-term development of the mammary gland. Only a few studies have attempted to use transcriptomics to study PAR and MFP. For example, Li and colleagues [16] have characterized the transcriptome of MFP and PAR tissue in 5 mo old ovariectomized dairy heifers treated with estrogen and they indentified several genes, most novel, thought to be regulated by estrogen in MFP and PAR, underscoring a crucial role for estrogen in mammary development at this post-weaning stage. We are not aware of similar transcriptomics characterizations in pre-weaned heifers.

In the present study, we have attempted to provide a molecular signature of PAR and MFP using microarrays to uncover genes more highly expressed in one tissue vs. the other. We have corroborated those data through bioinformatics analysis in order to uncover potential molecules involved in the putative crosstalk between PAR and MFP in pre-weaned heifers. Our approach is novel but has the "limitation" common to most bioinformatics tools, as well as transcript annotation, in that the results are based on human and rodent data, as opposed to bovine. However, the large evolutionary similarity between mammalian species (including bovine; [21]) provides reasonable confidence that most of the results in the present study will be suitable for future hypothesis-driven molecular biology experiments.

Mechanisms of tissue development inferred from functional genomics and pathway analysis

The majority of genes that could be classified as preferentially-expressed in MFP or PAR were associated with cellular growth and proliferation, cellular or tissue development, metabolism, and cellular movement (Tables 2 and 3 and Additional file 3). It was evident from IPA analysis, however, that functional findings encompassed cells other than epithelial cells or adipocytes, e.g., immune and blood/endothelial cells, and fibroblasts (Additional file 3). These responses are not entirely unexpected given our tissue collection protocol, i.e., non-homogeneous cell types within PAR and MFP.

Overall, it was noteworthy that DEG more highly expressed in PAR vs. MFP were enriched significantly and in a greater number (e.g. >230 DEG in cellular growth and proliferation) in the above categories. Furthermore, proliferation and cell cycle were enriched in GO analysis of DEG that were more highly-expressed in PAR vs. MFP (Figure 1 and 2). The evidently larger growth and proliferation of PAR relative to MFP uncovered by transcriptomics agrees with findings by Meyer and colleagues [22] who showed that PAR DNA (mg) was affected more by age than weight of the animal; whereas, the MFP was more affected by nutrition than age. In addition, even though it was not discussed by those authors, calculation of data reported by Meyer and colleagues [22] indicated a greater relative increase (>2-fold from 50 to 350 kg of body weight) in PAR than MFP weight and also greater relative increase in mg PAR DNA than MFP DNA.

In support of greater relative growth in PAR vs. MFP, which is by definition an increase in size or number, the top signaling pathways in DEG more highly expressed in PAR vs. MFP (Table 4; even though considered non-significant at an FDR ≤ 0.05) were all related to cell cycle progression (Table 2). Among those, calcium signaling through PI3K has been associated with cell proliferation and stimulation of quiescent cells to re-enter the cell cycle, thus, initiating mitosis [23]. In addition, the Wnt/β-catenin signaling pathway has been related to control of cell proliferation and differentiation [24,25]. Besides cell growth and proliferation, significantly-enriched functions among the genes with greater expression in PAR vs. MFP were related to apoptosis, cell adhesion, cell organization and biogenesis, and development (Figure 1 and 2, Table 2, and Additional file 3). Those data clearly indicated that besides growth in number and dimension of cells, the mammary PAR underwent a larger degree of re-organization, both within each cell and between cells to form a highly-organized tissue compared to MFP, which is consistent with greater differentiation in PAR compared to MFP (Table 2). The larger degree of organization to support growth in PAR vs. MFP also is supported by the significant enrichment and, apparently, induced Actin Cytoskeleton Signaling (Table 2).

The DEG more highly expressed in MFP vs. PAR indicated that MFP had a more predominant "metabolism-associated" DEG profile compared with PAR, i.e., lipid metabolism, molecular transport (includes lipids and fatty acids), and carbohydrate metabolism were among the most-enriched functions (Table 3, Figure 2). Furthermore, LXR/RXR and TR/RXR activation were among enriched signaling pathways (Table 4). To some extent these findings are not unexpected and could partly be explained by the fact that this tissue is mainly composed of adipocytes, a feature reflected by the higher level of expression of RXRG (ca. 8-fold vs. PAR), PPARG (ca. 3-fold), NR1D1 (ca. 2-fold), and KLF4 (ca. 2.5-fold) all of which belong to classical pro-adipogenic pathways (Additional file 2) [26,27].

When we considered the effect on function among DEG more highly expressed in MFP vs. PAR, the large enrichment of tissue development (including connective tissue; Table 3) indicated morphological remodeling. Although the same analysis appeared to suggest that overall cellular growth and proliferation was lower in MFP vs. PAR, it should be noted that several pro-adipogenic factors (e.g., PPARG, NR1D1, and KLF4) were substantially enriched in MFP vs. PAR, some of those known to be associated with continued/sustained mitotic clonal expansion (e.g., NR1D1, KLF4) of committed pre-adipocytes as well as terminal differentiation, maturation, and hypertrophy of adipocytes (e.g., PPARG, ADIPOQ) [27]. Because this entire process also is regulated by hormones such as insulin, it is likely that nutritional intervention at this early age leading to altered insulin profiles over the long-term could alter the extent of adipogenesis. This point also is supported by the significant enrichment of response to insulin stimulus in the GO analysis of DEG more highly expressed in MFP vs. PAR (Figure 2).

The enrichment of morphological remodeling in DEG more highly expressed in MFP over PAR was probably due to the increase in size of the cells by accumulation of triacylglycerol, as suggested both by the enrichment of functions related to synthesis and transport of lipid as well as synthesis of acyl-CoA and carbohydrate utilization and also the MFP composition data reported by Daniels et al. for the same heifers [10] (Table 3, Figure 1 and 2). In the present experiment, the weight of MFP of the selected heifers at slaughter averaged 193 g and represented ca. 95% of total mammary gland weight [10], thus, indicating that a large proportion of the mammary tissue in these animals was fat pad. Sinha and Tucker [28] reported a considerable enlargement of the MFP during the pre-weaning phase (birth to 2-3 mo age), and more recent work showed similar results [22,29]. More recently, Meyer and colleagues [22] reported that the increase in MFP weight was more related to enhanced hypertrophy than proliferation, as also evidenced by ratio g of tissue/mg DNA from 50 to 350 kg body weight which was very similar for PAR but increased ca. 3-fold in MFP between 50 kg and 350 kg heifers. Our results suggesting that MFP development was partly related to adipocyte hypertrophy and remodeling appear to be supported by previous findings.

More highly-expressed genes in MFP vs. PAR were preferentially associated with cellular movement (e.g., migration and invasion of cells) and cell-to-cell signaling (e.g., activation and adhesion of cells), functions that are closely related to the immune system (Table 3, Additional file 3). In addition, top signaling pathways among DEG more highly expressed in MFP vs. PAR were primarily related to response to stimuli, e.g., acute phase response, complement system, LPS/IL-1 mediated inhibition of RXR function, IL-8 signaling, and oxidative stress response (Table 4). These data suggest that a significant amount of the DEG more highly expressed in MFP vs. PAR are immune-related genes or that the MFP vs. PAR compartment is more enriched with immune cells (e.g., macrophages) as it has been observed previously in adipose tissue from rodents and humans [30]. Our network analysis (see section below), and previous heifer mammary proteomics data [31], provide evidence that these immune-related pathways in MFP vs. PAR might be biologically relevant in the context of mammary gland development.

We observed a more significant enrichment of IGF-1 signaling among DEG more highly expressed in MFP vs. PAR than vice versa (Table 4 and Additional file 3). This pathway was previously related to estrogen and its receptor in 5 mo-old heifers, where increased IGF1 expression was observed after estrogen treatment [5]. In heifers and rodents [32], circulating estrogen seems to act through its receptor and induces MFP cells to secrete IGF-1, which will then act in a paracrine fashion on PAR cells. At least in mice, such a mitogenic effect is observed in spite of low plasma levels of IGF-1 [33]. Our findings, however, do not support an induction of IGF-1 via stimulation of estrogen in MFP compared with PAR. Instead, data suggest that PAR was probably more sensitive to estrogen due to the greater mRNA abundance of ESR1 which agrees with a previous study [5] where mRNA expression of ESR1 was more predominant (ca. 3-fold greater) in PAR than MFP. Furthermore, no difference in IGF1 expression was detected by the microarray analysis or qPCR [10] between the two tissues.

Pathway analysis underscored a more prominent role of IGF-1 signaling in MFP than PAR, which suggests greater sensitivity of MFP to IGF-1 signaling. This is supported by the detailed visualization of the pathway (Additional file 3) which indicates that several key factors in the IGF-1 signaling cascade had greater expression in MFP compared with PAR (e.g., AKT, IRS1, Ras). Because both IGFBP5 and 6 were more abundant in MFP vs. PAR (Additional files 2 and 3), it could be possible that IGF1 activity was reduced [34]. A similar result was found for IGFBP5 and 6 via qPCR in the study of Daniels et al. [10]. As mentioned above, PAR is thought to be affected by local stimulators (e.g., IGF-1) derived from MFP but in our study it appeared that IGF-1 was not a major player in this tissue as it was not among the significantly-enriched pathways (Table 4). Our findings leave open the possibility that other paracrine factors secreted by PAR cells, which can be released in response to circulating estrogen [1], or by the MFP can play a more prominent role during this stage of development in bovine mammary gland.

Transcriptional network analysis reveals a central role for several transcription regulators in PAR and MFP development

Mammary parenchyma

Gene network analysis among DEG that were more highly-expressed in PAR vs. MFP (Figure 3) revealed that the transcription regulators MYC (oncogene) and TP53 (tumor suppressor) could play central roles. MYC has been detected in prepubertal [35] and lactating bovine mammary tissue [36]. Expression of MYC mRNA is increased by IGF-1 [37], and a primary response of MYC is to enhance epithelial cell proliferation. Recent evidence showed that MYC is essential in mediating Wnt-signaling and subsequent cell proliferation and growth and it also appears to control protein expression through mRNA translation [38,39], all of which can be considered important in PAR development.

TP53 can induce cell apoptosis or cell cycle arrest in response to different types of stress, i.e., it is considered to be a central control point of cell transformation and tumorigenesis [40,41]. TP53 activation can induce a transient (cell cycle arrest) or a permanent block of cell proliferation (senescence), or can induce the activation of cell death pathways in response to genotoxic stress [42]. A classical feature of TP53 activation in response to physiological stress (e.g., oxidative stress) leading to DNA damage is the activation of signaling cascades leading to DNA repair, recombination, and the control of DNA replication ultimately protecting cells from endogenous DNA damage [42]. Therefore, on one hand our functional analysis (Figures 1 and 2 and Tables 2 and 4) does not support low cell cycle activity in PAR but, rather, a larger cell cycle activity compared with MFP. On the other hand the more highly expressed DEG in PAR vs. MFP were enriched significantly in functions associated with DNA replication, recombination, and repair (Table 2 and Figure 2). In addition to MYC and TP53, CTNNB1 also may have potentially controlled expression of several genes that were more highly expressed in PAR vs. MFP (Figure 3). The biological significance of CTNNB1 in the present experiment is not apparent, but this protein is known to have a pivotal role in alveologenesis during pubertal mammary development in mice [43].

Mammary fat pad

Transcriptional gene networks generated among DEG with ≥ 1.5-fold expression in MFP vs. PAR (Figure 4) encompassed a large number of metabolism-related genes, revealing a pattern that would be expected in a classical adipocyte network [26], i.e., several stages of adipogenesis were revealed by some of the overexpressed genes. For example, expression (Figure 4) of NR1D1 and KLF4 was suggestive of pre-adipocytes undergoing mitotic clonal expansion; whereas, expression of PPARG and its putative targets FASN, acetyl-Coenzyme A carboxylase alpha [ACACA], mitochondrial glycerol-3-phosphate acyltransferase [GPAM], elongation of very long chain fatty acids-like 6 [ELOVL6], LPL, CD36 molecule (thrombospondin receptor) [CD36], FABP4, acyl-CoA synthetase long-chain family member 1 [ACSL1], and SCD was indicative of lipid filling that is more characteristic of a mature adipocyte [27]. Not surprisingly, network analysis showed that PPARG is a central transcription factor among genes more highly expressed in MFP vs. PAR, and its expression at this stage of development could be essential for adipogenesis and lipid filling in MFP compared with PAR [44] but also for the production of adipokines such as ADIPOQ and LEP both of which could exert control over tissue inflammation (e.g., through regulation of prostaglandin-endoperoxide synthase 2 [PTGS2]).

PPARG interacts closely with RXRG, RXRB, and potentially with lipin 1 (LPIN1) [45]. RXR transcription factors participate in the regulation of cholesterol and fatty acid metabolism by interacting with other transcription regulators such as PPARG (i.e., via protein-protein interactions not shown in Figure 4) [44]. LPIN1 is one of 3 isoforms that is associated with triacylglycerol synthesis in rodent adipose tissue [46] and fatty acid oxidation in rodent liver [45]. All lipin isoforms were more expressed in MFP vs. PAR (Additional file 2), with LPIN1 being the only one having >1.5-fold larger mRNA abundance in MFP vs. PAR. Several other genes in addition to those discussed are considered to play central roles in adipose tissue, including DGAT2 (>11-fold higher expressed in MFP vs. PAR). Although not present in the network shown in Figure 4, DGAT2 is considered to be a PPARG-target gene in non-ruminants and is related to triglyceride synthesis [44]. THRSP (potential RXRB target gene; Figure 4) is a lipogenic transcription factor and it is synergistically regulated by thyroid hormone and insulin [47,48] as well as long-term via the transcription regulator carbohydrate responsive element binding protein (MLXIPL or ChREBP).

Other central transcription factors uncovered by IPA analysis among genes more highly-expressed in MFP vs. PAR were FOS, HNF4A, EPAS1, KLF2, and EGR2 (Figure 4, Additional file 2). FOS has been previously associated with eukaryotic cell proliferation and differentiation [49,50]. However, the overall transcriptomics functional analysis of DEG more expressed in MFP vs. PAR indicated a lower degree of proliferation/differentiation in MFP vs. PAR (Figures 1 and 2 and Table 3). Transcription of FOS is induced with different potency by epidermal growth factor, PDGF, transforming growth factor beta [TGFB], tumor necrosis factor, FGF, IL-1, cAMP, estrogen, and other growth factors that have diverse roles in the cell [49,50]. Thus, in developing bovine mammary tissue this transcription regulator might be essential for MFP development and would provide a functional link with PAR-derived cytokines/growth factors such as IL1B (Figure 5, discussed below). A recent study in mouse mammary cells showed that TGFB regulates the expression of HNF4A, which in turn could control cell proliferation [51]. Interestingly, MFP expressed greater mRNA of TGFB1 compared with PAR (Additional file 2), which might indicate that, at the moment of harvesting, the influence of PAR-derived TGFB was not relevant for MFP. KLF2, was among the TF which potentially could have controlled expression of several DEG more expressed in MFP vs. PAR (Figure 4). The KLF2 protein is an embryonic stem cell-related factor which appears to control pluripotency driven by MYC [52].

Also noteworthy was the higher expression of EPAS1 in MFP vs. PAR. EPAS1 is preferentially expressed in vascular endothelial cells and plays a pivotal role in the formation of mature vascular tissue [53]. As indicated by functional analysis, angiogenesis was highly-enriched among DEG with ≥ 1.5-fold greater expression in MFP vs. PAR (Table 4); however, the functional analysis did not indicate a greater degree of angiogenesis in MFP compared with PAR. The MFP is essential for ductal morphogenesis in developing mammary tissue, but also it is subject to marked angiogenesis upon stimulation by some factors released from PAR [54]. To our knowledge, EGR2 has not been studied in the context of mammary gland development. Our results suggest that this protein influenced the expression of several genes more highly-expressed in MFP vs. PAR (Figure 4). Among those putative targets, EGR2 could have control over the expression of growth hormone receptor [GHR] (Figure 4). The greater expression of GHR in MFP vs. PAR agrees with qPCR data from Daniels and colleagues [10].

Bidirectional crosstalk between tissues inferred from network analysis

The potential effects of MFP-derived growth factors on the developing mammary PAR have been previously studied by several groups with different species [11,55,56]. Similarly, an effect of paracrine factors from PAR on MFP has been suggested from studies of normal mammary epithelial cells [12] and breast tumor epithelial cells [57]. These types of interactions could occur in prepubertal mammary tissue, for example, through growth factors and cytokines such as molecules from the EGF and FGF families [11], which if secreted by either tissue could then act in a paracrine fashion. Thus, to uncover additional factors that could play a role in PAR and MFP development, we evaluated the presence of cytokines and growth factors that might be secreted preferentially by one tissue or the other. The criteria for this analysis was implemented using IPA and was based on the principle that the cytokines and growth factors more highly-expressed in one tissue and, very likely, secreted by it, could potentially affect those DEG that are more highly-expressed in the other tissue (i.e., higher sensitivity to the cytokine or growth factor release in greater amounts by the other tissue) via effects at the level of mRNA expression, functional activation, protein modification, protein-protein interaction, protein-RNA interaction, protein phosphorylation, and protein translocation (see legends in Figure 5 and 6). The multitude of relationships, mostly from rodent and human studies, underscores the importance of these molecules on mammary development.

Model of putative tissue interactions

Gene network and pathway analysis using IPA allowed us to develop a putative model that would allow both mammary compartments to interact and elicit large-scale changes in the transcriptome and, potentially, tissue function via the synthesis of cytokines and growth factors, many of which have not been studied previously in the context of bovine mammary development. Although our study deals with mRNA expression only, it is assumed that higher expression of cytokines and growth factors in one tissue likely results in greater amount of protein synthesized and very likely more protein secreted. The effect of those signaling molecules can be both at the paracrine, autocrine, as well as endocrine levels. The molecules uncovered by our analysis as being preferentially expressed in PAR compared with MFP or vice versa could represent paracrine factors which allow for crosstalk between the two tissues and play important roles in proliferation and development of both epithelial and MFP cells (Figure 7). More importantly, these could be considered an important starting point for future detailed molecular studies of the interaction between PAR and MFP.

Signaling molecules likely released in greater amounts by PAR compared with MFP, through more mRNA abundance in MFP vs. PAR and vice versa, and potentially acting on MFP appear to affect, in both tissues, similar functions including cell movement, development, growth and proliferation, cell-mediated immune response, and cell-to-cell signaling and interactions. However, those signaling molecules likely released in greater amounts by PAR compared with MFP can potentially affect cell death in MFP cells; whereas, those released in greater amounts by MFP compared with PAR can potentially affect DNA metabolism in PAR (Figure 7). The potential effects of cytokines and growth factors on cellular growth, development, proliferation, and cell signaling in mammary tissue have been reported for both tissues in a recent review of the literature [4].

Four signaling pathways could have been potentially stimulated reciprocally between the two tissues (Figure 7). Among those were signaling pathways related to acute phase response, axonal guidance, glucocorticoid receptor, and IGF-1 response. The biological significance of the first two is not apparent. Also for the glucocorticod receptor signaling a biological significance is not apparent in our study, nonetheless this pathway has been previously associated with mammary development. In fact, glucocorticoid receptor signaling was associated with the normal development of the virgin mouse mammary gland through stimulation of ductal epithelial cell proliferation [58] and with regulating lobuloalveolar development during pregnancy [59].

The importance of IGF-1 signaling in PAR and MFP development has been discussed above and in previous papers/reviews (e.g., [4,60]). Interestingly, our data suggested that MFP likely affected IGF-1 signaling in PAR not through higher synthesis of IGF-1 but rather by modulating the down-stream signaling. More importantly, our data suggested an active role of PAR in affecting IGF-1 signaling in MFP and based on the number of genes affected and the fact that IGF1 signaling was more significantly enriched in DEG with greater expression in MFP vs. PAR, the stimulation of this pathway appears to be as strong or stronger in MFP than PAR (Table 4 and Figure 7).

Several canonical pathways were uniquely affected in one tissue by signaling molecules released in greater amounts from the other tissue (Figure 7). The cytokines and growth factors preferentially released from PAR compared with MFP seem to influence several nuclear receptors related to lipid metabolism (FXR/RXR, LXR/RXR/PPAR, and VDR/RXR activation), but also maintenance of pluripotency and immune-related pathways (NRF2-mediated oxidative stress) in MFP. The signaling molecules preferentially released from MFP vs. PAR seem to influence PAR signaling pathways with distinctive roles during this stage of development, e.g., actin cytoskeletal signaling (Table 5).

Additional pathways not found to be highly-significant when all DEG were analyzed (Table 4) also seemed to be affected by signaling molecules likely released in greater amounts by MFP compared with PAR, including ERK/MAPK signaling, FGF signaling, and PDGF signaling (Figure 7). ERK/MAPK signaling is related to a broad range of intracellular functions [61] and its role in developing mammary gland is not apparent. FGF signaling is related to the induction of potent mitogenic and angiogenic cellular processes and it is involved in embryonic [62] and postnatal mammary development [63] as well as with tumors of diverse origin, including mammary tumors [64]. PDGF signaling promotes mammary cancer progression and can induce apoptosis of human and murine mammary carcinoma cells when inhibited [65], suggesting a fundamental role in normal mammary differentiation and development. In this regard, Orr Urtreger and Lonai [66] revealed a possible interaction during organogenesis between epithelial and mesenchymal tissue in mice, i.e., PDGFA was expressed in mammary epithelial tissue and its receptor in the surrounding mesenchymal tissue. We are not aware of research in the prepubertal mammary gland of heifer calves that studied these signaling pathways.

Conclusions

We uncovered specific transcriptomic signatures characterizing genes with large difference in expression between MFP and PAR tissues. Not surprisingly, most of the genes that were more highly expressed in MFP vs. PAR were characteristic of adipose tissue, and those more expressed in PAR vs. MFP were characteristic of an epithelial tissue undergoing expansion and remodeling. Overall, our analyses suggested a large degree of interaction between the two tissues and allowed envisaging a reciprocal influence between the two tissues during this stage of development. This was indicated by the potential effect that the signaling molecules preferentially expressed in PAR vs. MFP and, likely released, have on lipid metabolism-related functions/pathways, which from our data was what distinguished the most those genes more highly expressed in MFP compared with PAR. Similarly, the cytokines and growth factors more highly expressed in MFP compared with PAR potentially affected the functions/pathways related to cell cycle, development, and proliferation in PAR, which our data highlighted as the main functions represented among the genes more highly-expressed in PAR vs. MFP.

Recent efforts have largely focused on MFP as a source of paracrine factors but our study clearly showed that PAR cells could play the same role. Based on the current analysis, the number of cytokines and growth factors that potentially are secreted in greater amounts by each tissue and affect molecules in the other underscores the concept of crosstalk already postulated by several investigators [11,16,67]. Ultimately, these bidirectional interactions might be required to coordinate mammary tissue development under normal circumstances or in response to environmental stimuli, such as nutrition.

Overall, the model generated based on the results from the present experiment predicts a large degree of crosstalk between MFP and PAR with a reciprocal regulation. The main factors at play appear to encompass several cytokines and growth factors preferentially released by PAR including SPP1, CXCL10, PDGFA, DFF1, and NRG1 which are probably slowing down the proliferation of MFP and increasing its lipid accumulation. Concomitantly, cytokines and growth factors released preferentially by MFP such as ADIPOQ, FGF2, GRP, and NOV are probably inducing major re-organization and proliferation of the PAR.

Methods

Animals and sampling

Samples used in this study were a subset obtained from a larger experiment [10,68,69]. All animal procedures were conducted under protocols approved by the Virginia Tech Institutional Animal Care and Use Committee. Specific details on feeding, management, and sample collection have been reported previously [10,68,69]. For the current experiment, 19 PAR and 21 MFP samples from 21 Holstein heifer calves (65 d-old; 77.5 ± 2.6 kg BW) representing animals on all diets reported in [10,68,69] were used. Additionally, samples were used to conduct a direct transcriptomics comparison between parenchyma and fat pad. For the present experiment, MFP tissue was harvested from the mammary fat adjacent to the body wall, while PAR tissue was harvested from the macroscopic epithelial portion of the gland. Subsamples of PAR and MFP were snap-frozen in liquid- N2, shipped overnight to the University of Illinois, then stored in liquid- N2 until use.

Extraction of RNA, cDNA synthesis, microarrays, and real-time PCR

Details of these procedures are reported in Additional file 1, particularly in Tables S1-S4 and Figure S1. Briefly, PAR and MFP tissues were weighed (~0.5 g) and total RNA extracted using ice-cold Trizol (Invitrogen, Carlsbad, CA). The purity of RNA (A260/A280) was above 1.9. RNA quality was assessed using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). Samples had a median RNA integrity (RIN) value of 7.7 ± 0.7. cDNA synthesis for microarrays was carried out with a total of 10 μg of RNA (~1 μg/μL). Methods for cleanup and aminoallyl-labeling of cDNA were described previously [70]. Briefly, the aminoallyl-labeled cDNA sample was dried using a speed-vac (Eppendorf Vacufuge® Concentrator, Eppendorf, Westbury, NY) for ~1 h and then resuspended in 4.5 μl 0.1 M sodium carbonate buffer (pH = 9.0). Four and a half microliters of the appropriate Cy dye ester (Cy3 or Cy5; Amersham, Piscataway, NJ) was added to couple the aa-cDNA and incubated for at least 1 h at room temperature. Removal of uncoupled dye was done using the Qiagen PCR Purification Kit.

A bovine oligonucleotide microarray developed at the University of Illinois [71] with > 13,000 bovine oligonucleotides (70-mers) was used to identify large-scale changes in gene expression. Details on the development, annotation, hybridization protocol, and scanning of arrays have been reported previously [71]. In order to increase reliability of data, the following filtering criteria were applied: only slides with ≥ 20,000 (out of > 27,000) spots with a median signal intensity ≥ 3 SD above background in both Cy3 and Cy5 channels and a mean intensity ≥ 400 relative fluorescent units in both Cy3 and Cy5 channels were used. The microarray data have been deposited in NCBI's Gene Expression Omnibus [72] and are accessible through GEO Series accession number GSE20363 http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20363 webcite.

cDNA to be used in qPCR was synthesized starting from 100 ng total RNA mixed with 1 μg dT18 (Operon Biotechnologies, Huntsville, AL), 1 μL 10 mmol/L dNTP mix (Invitrogen, Carlsbad, CA), 1 μL Random Primers (3 μg/μL, Invitrogen, Carlsbad, CA), and 7 μL DNase/RNase free water. A total of 9 μL of Master Mix composed of 4.5 μL 5× First-Strand Buffer (Invitrogen, Carlsbad, CA), 1 μL 0.1 M DTT (Invitrogen, Carlsbad, CA), 0.25 μL (100 U) of SuperScript™ III RT (Invitrogen, Carlsbad, CA), 0.25 μL of RNase Inhibitor (Promega, Madison WI), 3 μL DNase/RNase free water was added. The reaction was performed in an Eppendorf Mastercycler® Gradient (Eppendorf, Westbury, NY) using the following temperature program: 25°C for 5 min, 50°C for 60 min and 70°C for 15 min. cDNA was then diluted 1:3 with DNase/RNase free water. Four μL of diluted cDNA mixed with 5 μL of SYBR green (Applied Biosystems, Foster City, CA), 0.4 μL of each 10 μM primers, and 0.1 mL of DNase/RNase free water. For real-time RT-PCR each sample was run in triplicate to control reproducibility of the essay and a 4 point relative standard curve (4-fold dilution of cDNA originate from a pool RNA of all samples) plus the non-template control were used. The reactions were performed in an ABI Prism 7900 HT SDS instrument (Applied Biosystems, Foster City, CA) using the following conditions: 2 min at 50°C, 10 min at 95°C, 40 cycles of 15 s at 95°C, and 1 min at 60°C. PPP1R11, MTG1, RPS15A were used as internal control genes to normalize qPCR data [73]. Additional details are reported in Additional file 1.

Data analyses

Data from a total of 82 microarrays (38 PAR and 44 MFP; 41 samples from 19 animals contributing both PAR and MFP and 3 animals contributing only MFP) were normalized for dye and array effects (i.e., Lowess normalization and array centering) and used for statistical analysis. All data were analyzed using the Proc MIXED procedure of SAS (SAS, SAS Inst. Inc., Cary, NC). To determine differences in mRNA expression between PAR and MFP, the statistical analysis had to be conducted with both PAR and MFP data together, i.e., fixed effects in the model were tissue and dye while random effects included calf and microarray. Raw P values for the tissue effect were adjusted using Benjamini and Hochberg's FDR [74]. Differences in relative expression between PAR and MFP were considered significant at an FDR-adjusted P = 0.05 for tissue. For a more stringent characterization between the two tissues, a ≥ 1.5-fold difference in mRNA expression was set as threshold among the DEG. Data from qPCR were analyzed using the same statistical model described above. Differences were considered significant at P ≤ 0.05. The complete statistical output of the microarray analysis is available in Additional file 2.

Data mining

Data mining was performed using IPA (Ingenuity Systems, Inc., http://www.ingenuity.com webcite) after uploading into the system the entire microarray and qPCR data set with associated FDR and fold differences between PAR and MFP. In IPA, thresholds of FDR = 0.05 and a ≥ 1.5-fold difference were applied to filter significantly affected genes for function, pathway, and network analyses effects. The significance of the association between the dataset filtered by these thresholds and the IPA functions was calculated by IPA using a Benjamini-Hochberg's FDR ≤ 0.01 using the mapped genes on our microarray as background. For canonical pathway analysis we used FDR ≤ 0.05 because the 1.5-fold DEG in PAR vs. MFP did not enrich any pathways at an FDR ≤ 0.01. For a full interpretation of the data generated by the functional analysis in IPA we used the "effect on function" feature in IPA, which allowed determining among those significantly-enriched functions, which specific sub-function/s, tissue, or cellular component/s affected the function and in which direction (i.e., increase/decrease the function). Details of the analysis are reported in Additional file 3 and a summary in Tables 2 and 3.

The gene ontology (GO) analysis was performed using DAVID [75] following the criteria suggested [76]. A Benjamini-Hochberg FDR correction of P-value ≤ 0.05 was set as significant for all categories and all GO terms. The GO analysis was performed in both combined and separate lists of DEG with ≥ 1.5-fold difference between the two tissues. Each separate list of genes with ≥ 1.5-fold difference between the two tissues was used for interpretation of all DEG list using Excel software. Details of the methods are reported in Additional File 1 and full results are reported in Additional file 4 and main findings in Figures 1 and 2.

Network analyses

Transcriptional networks among transcription factors and DEG with ≥ 1.5-fold difference between PAR and MFP and molecular relationships between cytokines and growth factors preferentially expressed in one tissue vs. the other with DEG more highly expressed in the other tissue were uncovered/built using IPA features and knowledge base.

Abbreviations

IPA: ingenuity pathway analysis; GO: Gene Ontology; MFP: mammary fat pad; PAR: parenchyma; qPCR: quantitative real-time RT-PCR.

Authors' contributions

PP performed microarray and qPCR analyses and wrote the manuscript. MB participated in data mining using IPA and GO and wrote the manuscript. DG extracted RNA from tissues. KMD coordinated and helped perform the animal experiment and tissue collection, and helped draft the manuscript. SLR performed microarray statistical analysis. REE and HAL contributed new reagents and tools. RMA helped conceive and design the animal study. WLH helped with data interpretation. JJL designed the microarray experiment and wrote the manuscript. All authors read and approved the final manuscript.

Acknowledgements

Funding for the gene expression and bioinformatics analysis was provided in part by the Cooperative State Research, Education, and Extension Service, USDA, under Hatch projects ILLU-538-307 and ILLU-538-391 (both to JJL).

References

  1. Connor EE, Meyer MJ, Li RW, Van Amburgh ME, Boisclair YR, Capuco AV: Regulation of gene expression in the bovine mammary gland by ovarian steroids.

    J Dairy Sci 2007, 90(Suppl 1):E55-65. PubMed Abstract | Publisher Full Text OpenURL

  2. Hovey RC, McFadden TB, Akers RM: Regulation of mammary gland growth and morphogenesis by the mammary fat pad: a species comparison.

    J Mammary Gland Biol Neoplasia 1999, 4(1):53-68. PubMed Abstract | Publisher Full Text OpenURL

  3. Radcliff RP, Vandehaar MJ, Chapin LT, Pilbeam TE, Beede DK, Stanisiewski EP, Tucker HA: Effects of diet and injection of bovine somatotropin on prepubertal growth and first-lactation milk yields of Holstein cows.

    J Dairy Sci 2000, 83(1):23-29. PubMed Abstract | Publisher Full Text OpenURL

  4. Akers RM: Major advances associated with hormone and growth factor regulation of mammary growth and lactation in dairy cows.

    J Dairy Sci 2006, 89(4):1222-1234. PubMed Abstract | Publisher Full Text OpenURL

  5. Meyer MJ, Capuco AV, Boisclair YR, Van Amburgh ME: Estrogen-dependent responses of the mammary fat pad in prepubertal dairy heifers.

    J Endocrinol 2006, 190(3):819-827. PubMed Abstract | Publisher Full Text OpenURL

  6. Purup S, Sejrsen K, Akers RM: Effect of bovine GH and ovariectomy on mammary tissue sensitivity to IGF-I in prepubertal heifers.

    J Endocrinol 1995, 144(1):153-158. PubMed Abstract | Publisher Full Text OpenURL

  7. Silva LF, Liesman JS, Etchebarne BE, Nielsen MS, Vandehaar MJ: Short communication: intramammary infusion of IGF-I increases bromodeoxyuridine labeling in mammary epithelial cells of prepubertal heifers.

    J Dairy Sci 2005, 88(8):2771-2773. PubMed Abstract | Publisher Full Text OpenURL

  8. Capuco AV, Ellis S, Wood DL, Akers RM, Garrett W: Postnatal mammary ductal growth: three-dimensional imaging of cell proliferation, effects of estrogen treatment, and expression of steroid receptors in prepubertal calves.

    Tissue Cell 2002, 34(3):143-154. PubMed Abstract | Publisher Full Text OpenURL

  9. Capuco AV: Identification of putative bovine mammary epithelial stem cells by their retention of labeled DNA strands.

    Exp Biol Med (Maywood) 2007, 232(10):1381-1390. PubMed Abstract | Publisher Full Text OpenURL

  10. Daniels KM, Capuco AV, McGilliard ML, James RE, Akers RM: Effects of milk replacer formulation on measures of mammary growth and composition in Holstein heifers.

    J Dairy Sci 2009, 92(12):5937-5950. PubMed Abstract | Publisher Full Text OpenURL

  11. Imagawa W, Pedchenko VK, Helber J, Zhang H: Hormone/growth factor interactions mediating epithelial/stromal communication in mammary gland development and carcinogenesis.

    J Steroid Biochem Mol Biol 2002, 80(2):213-230. PubMed Abstract | Publisher Full Text OpenURL

  12. Hovey RC, Davey HW, Mackenzie DD, McFadden TB: Ontogeny and epithelial-stromal interactions regulate IGF expression in the ovine mammary gland.

    Mol Cell Endocrinol 1998, 136(2):139-144. PubMed Abstract | Publisher Full Text OpenURL

  13. Meyer MJ, Rhoads RP, Capuco AV, Connor EE, Hummel A, Boisclair YR, Van Amburgh ME: Ontogenic and nutritional regulation of steroid receptor and IGF-I transcript abundance in the prepubertal heifer mammary gland.

    J Endocrinol 2007, 195(1):59-66. PubMed Abstract | Publisher Full Text OpenURL

  14. Meyer MJ, Capuco AV, Ross DA, Lintault LM, Van Amburgh ME: Developmental and nutritional regulation of the prepubertal bovine mammary gland: II. Epithelial cell proliferation, parenchymal accretion rate, and allometric growth.

    J Dairy Sci 2006, 89(11):4298-4304. PubMed Abstract | Publisher Full Text OpenURL

  15. Thorn SR, Purup S, Vestergaard M, Sejrsen K, Meyer MJ, Van Amburgh ME, Boisclair YR: Regulation of mammary parenchymal growth by the fat pad in prepubertal dairy heifers: role of inflammation-related proteins.

    J Endocrinol 2008, 196(3):539-546. PubMed Abstract | Publisher Full Text OpenURL

  16. Li RW, Meyer MJ, Van Tassell CP, Sonstegard TS, Connor EE, Van Amburgh ME, Boisclair YR, Capuco AV: Identification of estrogen-responsive genes in the parenchyma and fat pad of the bovine mammary gland by microarray analysis.

    Physiol Genomics 2006, 27(1):42-53. PubMed Abstract | Publisher Full Text OpenURL

  17. Akers RM, Ellis SE, Berry SD: Ovarian and IGF-I axis control of mammary development in prepubertal heifers.

    Domest Anim Endocrinol 2005, 29(2):259-267. PubMed Abstract | Publisher Full Text OpenURL

  18. Akers RM, McFadden TB, Purup S, Vestergaard M, Sejrsen K, Capuco AV: Local IGF-I axis in peripubertal ruminant mammary development.

    J Mammary Gland Biol Neoplasia 2000, 5(1):43-51. PubMed Abstract | Publisher Full Text OpenURL

  19. Siddall NA, Lin JI, Hime GR, Quinn LM: Myc--what we have learned from flies.

    Curr Drug Targets 2009, 10(7):590-601. PubMed Abstract | Publisher Full Text OpenURL

  20. Ellis S, Capuco AV: Cell proliferation in bovine mammary epithelium: identification of the primary proliferative cell population.

    Tissue Cell 2002, 34(3):155-163. PubMed Abstract | Publisher Full Text OpenURL

  21. Elsik CG, Tellam RL, Worley KC, Gibbs RA, Muzny DM, Weinstock GM, Adelson DL, Eichler EE, Elnitski L, Guigo R, et al.: The genome sequence of taurine cattle: a window to ruminant biology and evolution.

    Science 2009, 324(5926):522-528. PubMed Abstract | Publisher Full Text OpenURL

  22. Meyer MJ, Capuco AV, Ross DA, Lintault LM, Van Amburgh ME: Developmental and nutritional regulation of the prepubertal heifer mammary gland: I. Parenchyma and fat pad mass and composition.

    J Dairy Sci 2006, 89(11):4289-4297. PubMed Abstract | Publisher Full Text OpenURL

  23. Berridge MJ, Bootman MD, Roderick HL: Calcium signalling: dynamics, homeostasis and remodelling.

    Nat Rev Mol Cell Biol 2003, 4(7):517-529. PubMed Abstract | Publisher Full Text OpenURL

  24. Dravid G, Ye Z, Hammond H, Chen G, Pyle A, Donovan P, Yu X, Cheng L: Defining the role of Wnt/beta-catenin signaling in the survival, proliferation, and self-renewal of human embryonic stem cells.

    Stem Cells 2005, 23(10):1489-1501. PubMed Abstract | Publisher Full Text OpenURL

  25. Nusse R: Wnt signaling and stem cell control.

    Cell Res 2008, 18(5):523-527. PubMed Abstract | Publisher Full Text OpenURL

  26. Rosen ED, MacDougald OA: Adipocyte differentiation from the inside out.

    Nat Rev Mol Cell Biol 2006, 7(12):885-896. PubMed Abstract | Publisher Full Text OpenURL

  27. Lefterova MI, Lazar MA: New developments in adipogenesis.

    Trends Endocrinol Metab 2009, 20(3):107-114. PubMed Abstract | Publisher Full Text OpenURL

  28. Sinha YN, Tucker HA: Mammary development and pituitary prolactin level of heifers from birth through puberty and during the estrous cycle.

    J Dairy Sci 1969, 52(4):507-512. PubMed Abstract | Publisher Full Text OpenURL

  29. Brown EG, Vandehaar MJ, Daniels KM, Liesman JS, Chapin LT, Forrest JW, Akers RM, Pearson RE, Nielsen MS: Effect of increasing energy and protein intake on mammary development in heifer calves.

    J Dairy Sci 2005, 88(2):595-603. PubMed Abstract | Publisher Full Text OpenURL

  30. Hotamisligil GS: Inflammation and metabolic disorders.

    Nature 2006, 444(7121):860-867. PubMed Abstract | Publisher Full Text OpenURL

  31. Daniels KM, Webb KE Jr, McGilliard ML, Meyer MJ, Van Amburgh ME, Akers RM: Effects of body weight and nutrition on mammary protein expression profiles in Holstein heifers.

    J Dairy Sci 2006, 89(11):4276-4288. PubMed Abstract | Publisher Full Text OpenURL

  32. Walden PD, Ruan W, Feldman M, Kleinberg DL: Evidence that the mammary fat pad mediates the action of growth hormone in mammary gland development.

    Endocrinology 1998, 139(2):659-662. PubMed Abstract | Publisher Full Text OpenURL

  33. Richards RG, Klotz DM, Walker , Diaugustine RP: Mammary gland branching morphogenesis is diminished in mice with a deficiency of insulin-like growth factor-I (IGF-I), but not in mice with a liver-specific deletion of IGF-I.

    Endocrinology 2004, 145(7):3106-3110. PubMed Abstract | Publisher Full Text OpenURL

  34. Flint DJ, Boutinaud M, Tonner E, Wilde CJ, Hurley W, Accorsi PA, Kolb AF, Whitelaw CB, Beattie J, Allan GJ: Insulin-like growth factor binding proteins initiate cell death and extracellular matrix remodeling in the mammary gland.

    Domest Anim Endocrinol 2005, 29(2):274-282. PubMed Abstract | Publisher Full Text OpenURL

  35. Musters S, Coughlan K, McFadden T, Maple R, Mulvey T, Plaut K: Exogenous TGF-beta1 promotes stromal development in the heifer mammary gland.

    J Dairy Sci 2004, 87(4):896-904. PubMed Abstract | Publisher Full Text OpenURL

  36. Bionaz M, Loor JJ: Identification of reference genes for quantitative real-time PCR in the bovine mammary gland during the lactation cycle.

    Physiol Genomics 2007, 29(3):312-319. PubMed Abstract | Publisher Full Text OpenURL

  37. Conover CA, Bale LK: Insulin-like growth factor I induction of c-myc expression in bovine fibroblasts can be blocked by antecedent insulin receptor activation.

    Exp Cell Res 1998, 238(1):122-127. PubMed Abstract | Publisher Full Text OpenURL

  38. Cole MD, Cowling VH: Transcription-independent functions of MYC: regulation of translation and DNA replication.

    Nat Rev Mol Cell Biol 2008, 9(10):810-815. PubMed Abstract | Publisher Full Text OpenURL

  39. Cowling VH, Cole MD: Turning the tables: Myc activates Wnt in breast cancer.

    Cell Cycle 2007, 6(21):2625-2627. PubMed Abstract | Publisher Full Text OpenURL

  40. Finlay CA, Hinds PW, Levine AJ: The p53 proto-oncogene can act as a suppressor of transformation.

    Cell 1989, 57(7):1083-1093. PubMed Abstract | Publisher Full Text OpenURL

  41. Cadwell C, Zambetti GP: The effects of wild-type p53 tumor suppressor activity and mutant p53 gain-of-function on cell growth.

    Gene 2001, 277(1-2):15-30. PubMed Abstract | Publisher Full Text OpenURL

  42. Kim E, Giese A, Deppert W: Wild-type p53 in cancer cells: when a guardian turns into a blackguard.

    Biochem Pharmacol 2009, 77(1):11-20. PubMed Abstract | Publisher Full Text OpenURL

  43. Hatsell S, Rowlands T, Hiremath M, Cowin P: Beta-catenin and Tcfs in mammary development and cancer.

    J Mammary Gland Biol Neoplasia 2003, 8(2):145-158. PubMed Abstract | Publisher Full Text OpenURL

  44. Desvergne B, Michalik L, Wahli W: Transcriptional regulation of metabolism.

    Physiol Rev 2006, 86(2):465-514. PubMed Abstract | Publisher Full Text OpenURL

  45. Finck BN, Gropler MC, Chen Z, Leone TC, Croce MA, Harris TE, Lawrence JC Jr, Kelly DP: Lipin 1 is an inducible amplifier of the hepatic PGC-1alpha/PPARalpha regulatory pathway.

    Cell Metab 2006, 4(3):199-210. PubMed Abstract | Publisher Full Text OpenURL

  46. Reue K, Zhang P: The lipin protein family: dual roles in lipid biosynthesis and gene expression.

    FEBS Lett 2008, 582(1):90-96. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  47. Chou WY, Cheng YS, Ho CL, Liu ST, Liu PY, Kuo CC, Chang HP, Chen YH, Chang GG, Huang SM: Human spot 14 protein interacts physically and functionally with the thyroid receptor.

    Biochem Biophys Res Commun 2007, 357(1):133-138. PubMed Abstract | Publisher Full Text OpenURL

  48. Sudo Y, Goto Y, Mariash CN: Location of a glucose-dependent response region in the rat S14 promoter.

    Endocrinology 1993, 133(3):1221-1229. PubMed Abstract | Publisher Full Text OpenURL

  49. Weisz A, Rosales R: Identification of an estrogen response element upstream of the human c-fos gene that binds the estrogen receptor and the AP-1 transcription factor.

    Nucleic Acids Res 1990, 18(17):5097-5106. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  50. Angel P, Karin M: The role of Jun, Fos and the AP-1 complex in cell-proliferation and transformation.

    Biochim Biophys Acta 1991, 1072(2-3):129-157. PubMed Abstract | Publisher Full Text OpenURL

  51. Ishikawa F, Nose K, Shibanuma M: Downregulation of hepatocyte nuclear factor-4alpha and its role in regulation of gene expression by TGF-beta in mammary epithelial cells.

    Exp Cell Res 2008, 314(10):2131-2140. PubMed Abstract | Publisher Full Text OpenURL

  52. Cotterman R, Knoepfler PS: N-Myc regulates expression of pluripotency genes in neuroblastoma including lif, klf2, klf4, and lin28b.

    PLoS One 2009, 4(6):e5799. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  53. Takeda N, Maemura K, Imai Y, Harada T, Kawanami D, Nojiri T, Manabe I, Nagai R: Endothelial PAS domain protein 1 gene promotes angiogenesis through the transactivation of both vascular endothelial growth factor and its receptor, Flt-1.

    Circ Res 2004, 95(2):146-153. PubMed Abstract | Publisher Full Text OpenURL

  54. Neville MC, Medina D, Monks J, Hovey RC: The mammary fat pad.

    J Mammary Gland Biol Neoplasia 1998, 3(2):109-116. PubMed Abstract | Publisher Full Text OpenURL

  55. Cunha GR, Hom YK: Role of mesenchymal-epithelial interactions in mammary gland development.

    J Mammary Gland Biol Neoplasia 1996, 1(1):21-35. PubMed Abstract | Publisher Full Text OpenURL

  56. Cunha GR, Young P, Hamamoto S, Guzman R, Nandi S: Developmental response of adult mammary epithelial cells to various fetal and neonatal mesenchymes.

    Epithelial Cell Biol 1992, 1(3):105-118. PubMed Abstract OpenURL

  57. Singer C, Rasmussen A, Smith HS, Lippman ME, Lynch HT, Cullen KJ: Malignant breast epithelium selects for insulin-like growth factor II expression in breast stroma: evidence for paracrine function.

    Cancer Res 1995, 55(11):2448-2454. PubMed Abstract | Publisher Full Text OpenURL

  58. Reichardt HM, Horsch K, Grone HJ, Kolbus A, Beug H, Hynes N, Schutz G: Mammary gland development and lactation are controlled by different glucocorticoid receptor activities.

    Eur J Endocrinol 2001, 145(4):519-527. PubMed Abstract | Publisher Full Text OpenURL

  59. Wintermantel TM, Bock D, Fleig V, Greiner EF, Schutz G: The epithelial glucocorticoid receptor is required for the normal timing of cell proliferation during mammary lobuloalveolar development but is dispensable for milk production.

    Mol Endocrinol 2005, 19(2):340-349. PubMed Abstract | Publisher Full Text OpenURL

  60. Li RW, Capuco AV: Canonical pathways and networks regulated by estrogen in the bovine mammary gland.

    Funct Integr Genomics 2008, 8(1):55-68. PubMed Abstract | Publisher Full Text OpenURL

  61. Kolch W: Coordinating ERK/MAPK signalling through scaffolds and inhibitors.

    Nat Rev Mol Cell Biol 2005, 6(11):827-837. PubMed Abstract | Publisher Full Text OpenURL

  62. Dillon C, Spencer-Dene B, Dickson C: A crucial role for fibroblast growth factor signaling in embryonic mammary gland development.

    J Mammary Gland Biol Neoplasia 2004, 9(2):207-215. PubMed Abstract | Publisher Full Text OpenURL

  63. Hovey RC, Davey HW, Vonderhaar BK, Mackenzie DD, McFadden TB: Paracrine action of keratinocyte growth factor (KGF) during ruminant mammogenesis.

    Mol Cell Endocrinol 2001, 181(1-2):47-56. PubMed Abstract | Publisher Full Text OpenURL

  64. Grose R, Dickson C: Fibroblast growth factor signaling in tumorigenesis.

    Cytokine Growth Factor Rev 2005, 16(2):179-186. PubMed Abstract | Publisher Full Text OpenURL

  65. Jechlinger M, Sommer A, Moriggl R, Seither P, Kraut N, Capodiecci P, Donovan M, Cordon-Cardo C, Beug H, Grunert S: Autocrine PDGFR signaling promotes mammary cancer metastasis.

    J Clin Invest 2006, 116(6):1561-1570. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  66. Orr-Urtreger A, Lonai P: Platelet-derived growth factor-A and its receptor are expressed in separate, but adjacent cell layers of the mouse embryo.

    Development 1992, 115(4):1045-1058. PubMed Abstract | Publisher Full Text OpenURL

  67. Akers RM: Selection for milk production from a lactation biology viewpoint.

    J Dairy Sci 2000, 83(5):1151-1158. PubMed Abstract | Publisher Full Text OpenURL

  68. Daniels KM, Hill SR, Knowlton KF, James RE, McGilliard ML, Akers RM: Effects of milk replacer composition on selected blood metabolites and hormones in preweaned Holstein heifers.

    J Dairy Sci 2008, 91(7):2628-2640. PubMed Abstract | Publisher Full Text OpenURL

  69. Hill SR, Knowlton KF, Daniels KM, James RE, Pearson RE, Capuco AV, Akers RM: Effects of milk replacer composition on growth, body composition, and nutrient excretion in preweaned Holstein heifers.

    J Dairy Sci 2008, 91(8):3145-3155. PubMed Abstract | Publisher Full Text OpenURL

  70. Loor JJ, Dann HM, Everts RE, Oliveira R, Green CA, Guretzky NA, Rodriguez-Zas SL, Lewin HA, Drackley JK: Temporal gene expression profiling of liver from periparturient dairy cows reveals complex adaptive mechanisms in hepatic function.

    Physiol Genomics 2005, 23(2):217-226. PubMed Abstract | Publisher Full Text OpenURL

  71. Loor JJ, Everts RE, Bionaz M, Dann HM, Morin DE, Oliveira R, Rodriguez-Zas SL, Drackley JK, Lewin HA: Nutrition-induced ketosis alters metabolic and signaling gene networks in liver of periparturient dairy cows.

    Physiol Genomics 2007, 32(1):105-116. PubMed Abstract | Publisher Full Text OpenURL

  72. Edgar R, Domrachev M, Lash AE: Gene Expression Omnibus: NCBI gene expression and hybridization array data repository.

    Nucleic Acids Res 2002, 30(1):207-210. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  73. Piantoni P, Bionaz M, Graugnard DE, Daniels KM, Akers RM, Loor JJ: Gene expression ratio stability evaluation in prepubertal bovine mammary tissue from calves fed different milk replacers reveals novel internal controls for quantitative polymerase chain reaction.

    J Nutr 2008, 138(6):1158-1164. PubMed Abstract | Publisher Full Text OpenURL

  74. Reiner A, Yekutieli D, Benjamini Y: Identifying differentially expressed genes using false discovery rate controlling procedures.

    Bioinformatics 2003, 19(3):368-375. PubMed Abstract | Publisher Full Text OpenURL

  75. Dennis G Jr, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, Lempicki RA: DAVID: Database for Annotation, Visualization, and Integrated Discovery.

    Genome Biol 2003, 4(5):P3. PubMed Abstract | BioMed Central Full Text OpenURL

  76. Huang da W, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.

    Nat Protoc 2009, 4(1):44-57. PubMed Abstract | Publisher Full Text OpenURL