Bioinformatics and Gene Network Analyses of the Swine Mammary Gland Transcriptome during Late Gestation.

We used the newly-developed Dynamic Impact Approach (DIA) and gene network analysis to study the sow mammary transcriptome at 80, 100, and 110 days of pregnancy. A swine oligoarray with 13,290 inserts was used for transcriptome profiling. An ANOVA with false discovery rate (FDR < 0.15) correction resulted in 1,409 genes with a significant time effect across time comparisons. The DIA uncovered that Fatty acid biosynthesis, Interleukin-4 receptor binding, Galactose metabolism, and mTOR signaling were among the most-impacted pathways. IL-4 receptor binding, ABC transporters, cytokine-cytokine receptor interaction, and Jak-STAT signaling were markedly activated at 110 days compared with 80 and 100 days. Epigenetic and transcription factor regulatory mechanisms appear important in coordinating the final stages of mammary development during pregnancy. Network analysis revealed a crucial role for TP53, ARNT2, E2F4, and PPARG. The bioinformatics analyses revealed a number of pathways and functions that perform an irreplaceable role during late gestation to farrowing.

Introduction

Milk production by the lactating sow is a major limiting factor for the development of the neonatal pig.1 The mammary gland of the sow undergoes significant development during gestation2,3 and lactation.4 Development of the mammary gland is particularly dramatic during the final 1/3 of the gestation period when the gland undergoes rapid morphogenesis and expansion of parenchymal tissue.3,5 as well as initiation of the early stage of lactogenesis and colostrogenesis just prior to parturition.6 The physiological adaptations marking this period are a consequence of shifting profiles of gene expression. Characterizing the role of genes involved in metabolic and signaling pathways during this period will provide a more detailed understanding of how the mammary gland undergoes this structural and functional development.

Transcriptomic studies provide highly-effective approaches to unravel complex biological behaviors at a genome-wide level.7 Longitudinal studies of transcriptome profiling data are of particular value in providing a detailed view of physiological adaptations in tissues.8 High-throughput technologies play a key role in understanding physiological states by analyzing the data obtained through pre-designed time course experiments.9 In addition, functional genomics studies provide a means of identifying new connections and interactions among the pathways, regulatory networks and structural organizations that underlie the function of tissues.10 Although previous studies using transcriptional profiling techniques have provided some preliminary insights into the differential expression of genes in sow mammary glands during the peripartum period,11 our understanding of differentially expressed genes (DEG) and metabolic or signaling pathways in mammary tissue in that species is limited.

KEGG is a biological database (http://www.genome.jp/kegg/) comprised of systematic, genome and chemical information using a combination of 19 other databases.12 The bioinformatics framework of DAVID (http://david.abcc.ncifcrf.gov/) relies principally on the overrepresented approach (ORA), and uses large gene lists for biological interpretation. Together, both tools and the Dynamic Impact Approach (DIA)13 provide a better option as compared to other traditional analysis approaches, which mainly rely on a reductionist approach and are inadequate for functional analysis of time-course experiments.14 Thus, compared with ORA, the DIA analysis captures the well-established and important biological features of the mammary gland during pregnancy.13

In this study, we explored the role of biological pathways during late gestation in the sow by analyzing transcriptomics data using the DIA previously developed by our group.13 The DIA uses manually-curated reference pathways from Kyoto Encyclopedia of Genes and Genomes (KEGG)15 and Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics resources.16 The objective of this study was to undertake an integrated functional genomics study of the swine mammary gland transcriptome during late gestation by means of DIA and to highlight the role of its essential biological pathways.

Materials and Methods

Animal sampling and RNA extraction

All procedures for animal treatments were performed by following the guidelines approved by the University of Illinois Institutional Animal Care and Use Committee. 15 cross-bred primiparous gilts in late gestation were used for this study.5,17 Gilts were killed by exsanguination at the University of Illinois Meat Science Laboratory on day 80 (n = 5), 100 (n = 5) or 110 (n = 5) of gestation. Only tissue from gilts with at least 4 normal fetuses was used. Mammary glands were removed at slaughter and parenchymal tissue dissected from the central portion of the fourth anterior gland (previously observed to be the most developed gland during late pregnancy in primiparous gilts).5

After sample collection and storage at −80 °C, the ribonucleic acid (RNA) was extracted using ice-cold Trizol (Invitrogen Corp., San Diego, CA). Agarose gel electrophoretic analysis of 28S and 18S rRNA was performed to assess the integrity of RNA. The purity and concentration of RNA were assessed using the NanoDrop ND-1000 spectrophotometer and deemed to be high enough to carry out microarray hybridization. A swine oligonucleotide microarray (70 mer) containing 13,290 functional hits was used for transcript profiling.18 A dye-swap reference design was conducted using three time points (80, 100, and 110) during late pregnancy according to procedures published previously.19 The array annotation was based on similarity searches using BLASTN and BLASTX20 against human, mouse and pig UniGene databases, the human genome, and pig TIGR.21

Statistical analysis

Data were normalized for dye and array effects before statistical analysis as previously described.19 Microarray spots with median intensity standard deviation above the median of the background were applied as filters in order to ensure high-quality data. The data were analyzed using a PROC MIXED model of SAS.19 The resulting P-values were adjusted according to the Benjamini and Hochberg’s correction method.22 An ANOVA with a false discovery rate (FDR) correction resulted in 10.6% (1,409) DEG with a time effect (FDR < 0.15). The fold-change (FC) values were calculated using relative time points. Figure 1 contains the DEG at the 3 time point comparisons (namely, 100 versus (vs.) 80, 110 vs. 80, and 110 vs. 100 d). Tables 1–6 contain information on DEG that had a fold-change greater or lower than 2-fold at the 3 time point comparisons. The entire list of DEG can be found in Additional File 1.

Figure 1

Differentially expressed genes in pig mammary gland across different time point comparisons during late pregnancy (ie, 100 vs. 80, 110 vs. 80, and 110 vs. 100 d).

Table 1

List of DEG (FDR < 0.15 and raw P < 0.10) with 2-fold or greater upregulation on day 100 vs. 80. The Entrez Gene ID for each gene can be found in Additional File 1.

Gene symbol	Description	Fold
TF	Transferrin	4.51
HPR	Haptoglobin-related protein	4.03
TTN	Titin	3.61
UGP2	UDP-glucose pyrophosphorylase 2	3.52
OGFOD3	2-oxoglutarate and iron-dependent oxygenase domain containing 3	3.47
VIPR1	Vasoactive intestinal peptide receptor 1	3.17
DDX42	DEAD (Asp-Glu-Ala-Asp) box polypeptide 42	3.13
KIAA1755	KIAA1755	3.10
PROC	Protein C (inactivator of coagulation factors Va and VIIIa)	2.93
KLKB1	Kallikrein B, plasma (Fletcher factor) 1	2.89
TGDS	TDP-glucose 4,6-dehydratase	2.83
PIGR	Polymeric immunoglobulin receptor	2.77
PAX8	Paired box gene 8	2.75
C7	Complement component 7	2.73
BTN1A1	Butyrophilin, subfamily 1, member A1	2.70
NUDT7	Nucleoside diphosphate-linked moiety X motif 7	2.69
TXNDC17	Thioredoxin-like 5	2.67
RARB	Retinoic acid receptor, beta	2.66
SALL2	Sal-like 2 (Drosophila)	2.63
LSP1	Lymphocyte-specific protein 1	2.62
ABCG2	ATP-binding cassette, sub-family G (WHITE), member 2	2.62
MST1	Macrophage stimulating 1 (hepatocyte growth factor-like)	2.56
CPS1	Carbamoyl-phosphate synthetase 1, mitochondrial	2.52
ARRDC3	Arrestin domain containing 3	2.49
S100G	S100 calcium binding protein G	2.47
ACY1	Aminoacylase 1	2.45
LHX9	LIM homeobox 9	2.41
PDZD3	PDZ domain containing 3	2.37
KRT17	Keratin 17	2.36
SPTBN1	Spectrin, beta, non-erythrocytic 1	2.35
MYH4	Myosin, heavy chain 4, skeletal muscle	2.34
NNAT	Neuronatin	2.32
RNF11	Ring finger protein 11	2.32
COL8A1	Collagen, type VIII, alpha 1	2.31
B3GALTL	Beta 1,3-galactosyltransferase-like	2.30
KLC4	Kinesin light chain 4	2.28
MECP2	Methyl CpG binding protein 2 (Rett syndrome)	2.28
FTCD	Formiminotransferase cyclodeaminase	2.27
WHSC1L1	Wolf-Hirschhorn syndrome candidate 1-like 1	2.26
CYP2E1	Cytochrome P450, family 2, subfamily E, polypeptide 1	2.25
UVRAG	UV radiation resistance associated gene	2.25
TMCO5A	Transmembrane and coiled-coil domains 5A	2.20
IGJ	Immunoglobulin j polypeptide	2.19
TNK1	Tyrosine kinase, non-receptor, 1	2.17
PRRC2C	Proline-rich coiled-coil 2C	2.17
EEPD1	Endonuclease/exonuclease/phosphatase family domain containing 1	2.15
FBP1	Fructose-1,6-bisphosphatase 1	2.15
FMO1	Flavin containing monooxygenase 1	2.14
ACACA	Acetyl-Coenzyme A carboxylase alpha	2.13
IER5	Immediate early response 5	2.12
MT1A	Metallothionein 1A (functional)	2.10
FLRT2	Fibronectin leucine rich transmembrane protein 2	2.09
TESK1	Testis-specific kinase 1	2.08
IDH2	Isocitrate dehydrogenase 2 (NADP+), mitochondrial	2.08
HHIPL2	HHIP-like 2	2.07
CLEC18C	C-type lectin domain family 18, member C	2.06
GREB1	GREB1 protein	2.06
ITGA2B	Integrin, alpha 2b	2.05
PLOD3	Procollagen-lysine, 2-oxoglutarate 5-dioxygenase 3	2.05
PDIA4	Protein disulfide isomerase family A, member 4	2.04
AKR1B1	Aldo-keto reductase family 1, member B1 (aldose reductase)	2.02
ERLIN2	SPFH domain family, member 2	2.02
ATF5	Activating transcription factor 5	2.01
WDR26	WD repeat domain 26	2.01
CELF1	CUGBP, Elav-like family member 1	2.00

Table 6

List of DEG (FDR < 0.15 and raw P < 0.10) with 2-fold or greater downregulation on day 110 vs. 100. The Entrez Gene ID for each gene can be found in Additional File 1.

Gene symbol	Description	Fold
ARL6IP6	ADP-ribosylation-like factor 6 interacting protein 6	−3.26
OGFOD3	2-oxoglutarate and iron-dependent oxygenase domain containing 3	−2.60
KLKB1	Kallikrein B, plasma (Fletcher factor) 1	−2.57
NDRG1	N-myc downstream regulated gene 1	−2.26
TBC1D13	TBC1 domain family, member 13	−2.24
TDP2	Tyrosyl-DNA phosphodiesterase 2	−2.22
ELK3	ELK3, ETS-domain protein (SRF accessory protein 2)	−2.20
UGP2	UDP-glucose pyrophosphorylase 2	−2.03

Dynamic impact approach and enrichment analysis (DIA)

The newly developed DIA and the well-established DAVID bioinformatics resource were used to analyze the gene expression data. DIA uses the systems information from the KEGG database and ranks pathways based on their overall impact values across different time points. The detailed methodology of DIA analysis is described in.23 The DIA calculates the overall impact and flux (direction of impact) of biological pathways. For this purpose, the whole dataset (Additional File 1: sheet ‘Dataset’) with Entrez gene IDs, FDR, FC, and raw P-values of each time comparison were uploaded to calculate the impact and flux. An overall cut-off FDR < 0.15 and P-value < 0.05 between time points were applied as thresholds. The resultant summary of the DIA KEGG pathways (encompassing categories and subcategories) is shown in Figure 2. The DIA resulted in 195 biologically-impacted KEGG pathways (Additional File 2: sheet ‘DIA KEGG’). For enrichment analysis, the KEGG pathways were downloaded from the DAVID database with enrichment P-values and Benjamini-Hochberg FDR values (Additional File 4).

Figure 2

Overall summary of the KEGG (http://www.genome.jp/kegg/pathway.html) categories and sub-categories resulting from DIA analysis. Three time comparisons are shown viz., 100 vs. 80 d, 110 vs. 80 d and 110 vs. 100 d. Each time comparison consists of impact and flux. The impact values are shown in purple ranging from 0 (not impacted) to 15 (highly impacted), while flux values are represented from −10 (green = highly inhibited) to 10 (red = highly activated).

Gene ontology (GO) term analysis

GO terms such as biological process (BP), cellular components (CC), and molecular functions (MF) were downloaded from the DAVID database and were analyzed using DIA. For this purpose, a multigene list of both up and down regulated Entrez gene IDs was uploaded to the DAVID database.24 For background information, we uploaded Homo sapiens Entrez gene IDs corresponding to our customized oligo array.

KEGG pathways visualization

KegArray tool25 was used to analyze the KEGG pathways among DEG. For this purpose, a list of Entrez gene ID’s with their respective FC values was uploaded for each time comparison. An FC cut-off value of 1 was applied as a threshold criterion. The results are listed in Additional File 3 and are discussed in the sections below.

Clustering results

The ‘Functional Annotation Clustering’ results were also downloaded from the DAVID database. To select the most significant results, we applied an enrichment score cutoff criterion of ≥1.3. The clustering results can be found in Additional File 4.

Gene network analysis

This analysis was performed as describe previously by Piantoni et al.26 Briefly, we used Ingenuity Pathway Analysis® (IPA) to identify transcription regulators and their networks with other genes within the list of significant DEG at each of the 3 time comparisons. Software features and the IPA knowledge base were used for the analysis.

Results

Coverage of microarray elements and DEG

Out of 13,290 unique spots on the microarray, a total of 12,313 oligo IDs were annotated with human Entrez IDs which constituted approximately 92.65% of the total coverage. To identify the longitudinal transcriptional gene response during late pregnancy we relied on DEG between 100 vs. 80, 110 vs. 80, and 110 vs. 100 days. The 100 vs. 80 days comparison represents the difference in gene expression patterns between the relatively immature mammary gland at approximately 34 days prior to parturition and the substantially more developed gland at about 14 days prepartum in the primigravid gilt.3,5 The 110 vs. 80 days comparison represents the difference in gene expression patterns between the immature gland and the gland that is near full term at about 4 days prepartum. In contrast, the 110 vs. 100 days comparison represents the difference in mammary tissue between a time of rapid growth and limited functional differentiation at 100 days of gestation and the time at 110 days when the gland continues growing but also is entering into the functional differentiation stages associated with lactogenesis.

At P ≤ 0.05 and FDR ≤ 0.15, a total of 951 genes (348 up- and 603 down-regulated genes) at 100 vs. 80 days, 878 genes (336 up- and 542 down-regulated genes), and 221 genes (132 up- and 89 down-regulated genes) genes were differentially expressed (Fig. 1). Among these, only 701, 528 and 207 genes, respectively, were annotated with human Entrez gene IDs. For the 100 vs. 80 days and 110 vs. 80 days comparisons more genes were down-regulated as compared to up-regulated.

Overall summary of KEGG categories

Figure 2 provides a summary of KEGG categories (also see Additional File 2: Sheet ‘KEGG Summary’). The ‘Metabolism’ (M) category was up- regulated, whereas the ‘Genetic Information Processing’ category (GIP), ‘Environmental Information Processing’ category (EIP), ‘Cellular Processes’ category (CP) and ‘Organismal systems’ category (OS) were down-regulated in the 100 vs. 80 days comparison and the CP and OS categories were down- regulated in the 110 vs. 80 days comparison. In contrast, all categories and sub-categories were up-regulated in the 110 vs. 100 days comparison (Fig. 2).

Of the Metabolism sub-categories, metabolism of carbohydrate, energy, lipid, amino acids, cofactors and vitamins, and terpenoids and polyketides were up-regulated in the 100 vs. 80 days comparison. This up-regulated state was also observed in the 110 vs. 80 days comparison for metabolism of carbohydrate, energy, lipid, cofactors and vitamins, and terpenoids and polyketides, but not for amino acid metabolism. Expression of pathways associated with metabolism of carbohydrate, energy, lipid, and cofactors and vitamins continued to be up-regulated in the 110 vs. 100 days comparison, although with a more limited impact value. In contrast, expression pathways associated with metabolism of amino acids, terpenoids and polyketides was unchanged in that latter comparison. Carbohydrate metabolism consistently had the greatest impact factor in each day comparison.

Transcription, translation and replication and repair sub-categories of GIP and the signal transduction sub-category of EIP each were down-regulated in the 100 vs. 80 days and 110 vs. 80 days comparisons, but were up-regulated in the 110 vs. 100 days comparison. In contrast, the folding, sorting and degradation sub-category of GIP was consistently upregulated in each comparison, as was the membrane transport and signaling molecules and interaction sub-categories of EIP.

All sub-categories of CP and OS (except environmental adaptation) were strongly down-regulated in the 100 vs. 80 days comparison. This expression relationship was also seen in the 110 vs. 80 days comparison, with the exception of the immune system sub-category of the OS category, which was up- regulated. All sub-categories of the CP and OS categories were up-regulated in the 110 vs. 100 days comparison except for sensory system, which was not different in that comparison.

Most impacted KEGG pathways

The DIA analysis provides an evaluation of the impact of pathways within sub-categories of KEGG (Additional File 2: Sheet ‘DIA KEGG’). The overall most-impacted pathway was fatty acid biosynthesis. This impact was particularly strong in the 100 vs. 80 days comparison, continued to be a high impact in the 110 vs. 80 days comparison, and then remained unchanged in the 110 vs. 100 days comparison. Galactose metabolism, one carbon pool by folate, glyoxalate and dicarboxylate metabolism, and linoleic acid metabolism were among the 25 overall most impacted pathways among the metabolism subcategories.

The intestinal immune network for the IgA production pathway of the Immune System was found to be highly impacted in the overall analysis, with the highest impacts in 110 vs. 80 days and 110 vs. 100 days comparisons. In addition, 5 pathways of the Immune System Diseases sub-category were among the 25 most impacted pathways in the DIA analysis (Additional File 2: Sheet ‘DIA KEGG’).

Several other pathways were identified as overall highly impacted by the DIA analysis. The ECM- receptor interaction pathway of the Signaling Molecules and Interaction sub-category seemed to switch from up-regulation in the 100 vs. 80 days comparison, to down-regulation in the 110 vs. 80 days comparison, and back to up-regulation in the 110 vs. 100 days comparison. The vascular smooth muscle contraction pathway (circulatory System sub- category) and gap junction pathway (cell communication sub-category) were consistently down-regulated in the 100 vs. 80 days and 110 vs. 80 days comparisons, and then up-regulated in the 110 vs. 100 days comparison. In contrast, the ATP-binding cassette (ABC) transporters pathway (membrane transport sub-category) and the protein export pathway (folding, sorting and degradation sub-category) were up-regulated or not different in the 100 vs. 80 days comparison, strongly up-regulated in the 110 vs. 80 days comparison and up-regulated in the 110 vs. 100 days comparison.

Interestingly, 4 pathways of the cardiovascular diseases sub-category were among the 25 most impacted sub-categories in the DIA analysis, including arrhythmogenic right ventricular cardiomyopathy (ARVC), hypertropic cardiomyopathy (HCM), dilated cardiomyopathy (DCM) and viral myocarditis. In each of the latter pathways, the highest impact factors were observed in the 100 vs. 80 days or 110 vs. 80 days comparisons.

From the KEGG pathways analyzed by the overrepresentation approach (ORA) as shown in Additional File 4 (at an FDR enrichment P-value of 0.05) it could be discerned that the vascular smooth muscle contraction and long-term potentiation pathways both were most significantly enriched and were down- regulated at 100 vs. 80 days. The ECM-receptor interaction pathway was up-regulated in the 100 vs. 80 days and 110 vs. 80 days comparisons (Additional File 2, DIA KEGG sheet; Additional File 3); gap junction and intestinal immune network for IgA production pathways both were up-regulated in the 110 vs. 100 days comparison (Additional File 3; DIA KEGG sheet). The cell cycle and neurotrophin signaling pathways both were down-regulated in the 100 vs. 80 days and 110 vs. 80 days comparisons (Additional File 2, DIA KEGG sheet), while the Jak-STAT signaling pathway was up-regulated in the 110 vs. 80 days comparison (Fig. 3). The complement and coagulation cascades pathway was up-regulated in the 100 vs. 80 days comparison. These pathways can also be explored in Additional File 3 which lists the pathways obtained from the KegArray tool.

Figure 3

Changes in flux during late-pregnancy of the most-impacted KEGG signaling pathways in mammary tissue.

Gene ontology (GO) and cluster analysis

The DIA analysis was used to evaluate GO terms under the 3 main GO classes, including Biological Processes (BP), Cellular Components (CC) and Molecular Function (MF), as summarized in Figure 4 and Additional File 2 (Sheets ‘GOTERM BP’, ‘GOTERM CC’ and ‘GOTERM MF’, respectively). Among the ten most-enriched BP, DNA double- strand break processing involved in repair via single-strand annealing, connective tissue growth factor production, connective tissue growth factor biosynthetic process, and double-strand break repair via single-strand annealing each were only mildly changed in the 100 vs. 80 days comparison, but then strongly impacted and up-regulated in the 110 vs. 80 days and 110 vs. 100 days comparisons. Negative regulation of mRNA 3′-end processing was strongly impacted and up-regulated only in the 110 vs. 100 days comparison. Positive regulation of isotype switching to IgE isotypes was most strongly impacted and up-regulated in the 110 vs. 80 days comparison. Response to anoxia and ERK1 and ERK2 cascade were not changed in the 100 vs. 80 days comparison, however both became highly impacted and down-regulated in the 110 vs. 80 and 110 vs. 100 days comparisons.

Figure 4

The top 10 most enriched GO terms within Biological Processes (BP), Cellular Components (CC) and Molecular Function (MF). The ‘overall’ column of impact and flux was calculated by taking into account the values for each GO term at each of the three time comparisons. The range of impact and flux values is shown in the legends at the top of the Figure.

Among the 10 most-enriched CC, the Six1-Six4 complex were not changed in the 100 vs. 80 days comparison, and then were highly impacted and up-regulated in the 110 vs. 80 and 110 vs. 100 days comparisons. Collagen type VIII and short-chain collagen were strongly impacted and up-regulated in the 100 vs. 80 and 110 vs. 80 days comparisons, with no differences observed in the 110 vs. 100 days comparison. In contrast, collagen type III was only impacted in the 110 vs. 80 days comparison when it was down-regulated. Glycerol-3-phosphate dehydrogenase complex and intercellular canaliculus both were not different in the 100 vs. 80 days comparison, then up-regulated in the 110 vs. 80 and 110 vs. 100 days comparisons. In contrast to the above, PCNA complex was highly impacted and down-regulated in the 100 vs. 80 days comparison, but not observed to be different in the other comparisons.

Of the ten most-enriched MF, 3 patterns emerged (Fig. 4). Crossover junction endodeoxyribonuclease activity, interleukin 4 receptor, methylenetetrahydrofolate reductase (NADPH) activity, sodium-exporting ATPase activity, phosphorylative mechanism, bile acid-exporting ATPase activity, and flap endonuclease activity each were not observed to be different in the 100 vs. 80 days comparison, then were highly impacted and up-regulated in the 110 vs. 80 and 110 vs. 100 days comparisons. In contrast, UTP:glucose-1-phosphate uridylyltransferase activity, UTP-monosaccharide-1-phosphate uridylyltransferase activity, and dTDP-glucose 4,6-dehydratase activity each were highly impacted and up-regulated in the 100 vs. 80 days comparison, minimally impacted in the 110 vs. 80 days comparison, and down-regulated in the 110 vs. 110 days comparison. The glycoprotein-N-acetylgalactosamine 3-beta-galactosyltransferase activity was not different in the 100 vs. 80 days comparison, however, was highly impacted and down-regulated in the 110 vs. 80 and 110 vs. 100 days comparisons.

Cluster analysis27 was used to develop an integrated understanding of cellular processes in which genes involved in most pathways tend to share similar roles and functions. Also, the clustering analysis denoting 80, 100 and 110 days of late gestation using diverse annotation sources will provide a different representation of functional groups of GO terms. Given the conspicuous differences during late gestation in the clustering analysis, we centered on comparisons of possible roles of annotated pathways.

The Gene Set Enrichment Analysis (GSEA) generated many pathways that were enriched in either the 80, 100 or 110 days samples (Additional File 4). From these results, cell cycle, vascular smooth muscle contraction, and long-term potentiation were inhibited and enriched in the 100 vs. 80 days comparison (Additional File 4). Particularly, the expression of genes down-regulated in the 100 vs. 80 days comparison was enriched in GnRH signaling pathway (Additional File 2: Sheet ‘DIA KEGG’; Additional File 3). In addition, the expression of genes that were downregulated in the 110 vs. 80 days comparison may be associated with activities such as ECM-receptor interaction, focal adhesion, and spliceosome (Additional File 2: Sheet ‘DIA KEGG’).

This analysis focused on uncovering relationships between transcription factors and all the DEGs by means of the IPA knowledgebase, and specifically on published evidence at the level of gene expression (E), transcription (T), and protein-DNA interactions (PD) (Figs. 5–7). In the comparison of 100 vs. 80 days (Fig. 5), this analysis identified transcriptional networks encompassing 148 DEG and 70 transcription factors (67 transcription regulators and 3 ligand-dependent nuclear receptors). Similarly, the comparison of 110 vs. 80 days (Fig. 6) revealed 105 DEG and 74 transcription factors (no ligand-dependent nuclear receptors). The comparison of 110 vs. 100 days (Fig. 7) uncovered 33 DEG and 58 transcription factors (54 transcription regulator and 4 ligand-dependent nuclear receptors).

Figure 5

Transcriptional networks encompassing 148 DEG and 70 transcription factors (67 transcription regulators and 3 ligand-dependent nuclear receptors) generated by Ingenuity Pathway Analysis® using the list of DEG in the comparison of 100 vs. 80 d. The analytical approach was set to uncover relationships between transcription factors and all the DEGs in the comparison combined (with expression [denoted by E in the edge of the line], transcription [T], and protein-DNA interactions [PD]). The transcription factors are denoted by large bold fonts. The shapes are filled with red when consistently up-regulated and green when consistently down-regulated at 100 vs. 80 d. Arrows denote the direction of the interaction and solid lines denote direct and dashed lines indirect interactions.

Figure 7

Transcriptional networks encompassing 33 DEG and 58 transcription factors (67 transcription regulators and 3 ligand-dependent nuclear receptors) generated by Ingenuity Pathway Analysis® using the list of DEG in the comparison of 110 vs. 100 d. The analytical approach was set to uncover relationships between transcription factors and all the DEGs in the comparison combined (with expression [denoted by E in the edge of the line], transcription [T], and protein-DNA interactions [PD]). The transcription factors are denoted by large bold fonts. The shapes are filled with red when consistently up-regulated and green when consistently down-regulated at 110 vs. 100 d. Arrows denote the direction of the interaction and solid lines denote direct and dashed lines indirect interactions.

Figure 6

Transcriptional networks encompassing 105 DEG and 74 transcription factors (67 transcription regulators and 3 ligand-dependent nuclear receptors) generated by Ingenuity Pathway Analysis® using the list of DEG in the comparison of 110 vs. 80 d. The analytical approach was set to uncover relationships between transcription factors and all the DEGs in the comparison combined (with expression [denoted by E in the edge of the line], transcription [T], and protein-DNA interactions [PD]). The transcription factors are denoted by large bold fonts. The shapes are filled with red when consistently up-regulated and green when consistently down-regulated at 110 vs. 80 d. Arrows denote the direction of the interaction and solid lines denote direct and dashed lines indirect interactions.

Discussion

Lactogenesis: general considerations

The primiparous gilts used in this study had an average expected gestation length of approximately 114 days. As such, 80 days of gestation represents a time when the gland is starting to undergo a dramatic increase in the rate of growth.2,3,6 At 100 days of gestation, or 14 days prepartum, the gland has increased in mass significantly, increased in total gland protein and decreased in total gland adipose relative to 80 days.3 The 110 days tissue represents about 4 days prepartum when the gland is still growing, but also has shifted substantially into lactogenesis, or the initiation of lactation, leading up to parturition. The histological and molecular change of the gland during pregnancy can be appreciated by the increase in cross-section area and total DNA, as well as the breadth of pathways affected (Fig. 8).

Figure 8

Summary figure of the most relevant biological functions in pig mammary gland during late gestation as uncovered by the Dynamic Impact Approach (DIA) using KEGG pathways and Gene Ontology (GO) terms. Shown at the center are the changes in mammary gland cross-section area, DNA content, and morphometry during late gestation. Red titles are the most relevant metabolic pathways as revealed by the KEGG and GO analysis. Other color titles denote main functions impacted during late gestation. The X-axis depicts the day relative to parturition and the Y-axis denote the flux (the direction of the impact) as calculated by DIA.

Abbreviations: AA, amino acid; Arg, arginine; Ca, calcium; Cu, copper; Cys, cysteine; ER, endoplasmic reticulum; FA, fatty acids; GnRH, Gonadotropin-releasing hormone pathway; His, histidine; IL-4, interleukin-4; LCFA, long-chain fatty acids; Lys, lysine; MAPK, mitogen-activated protein kinase; Met, methionine; mTOR, mammalian target of rapamycin signaling pathway; Na, sodium; Nod, Nucleotide Oligomerization Domain; Notch, Notch signaling pathway; Phe, phenylalanine; Pro, proline; Ser, serine; TCA cycle, tricarboxylic acid cycle; TF, transcription factor; Trp, tryptophan; Val, valine; VEGF, vascular endothelial growth factor signaling pathway.

Lactogenesis is described as occurring in 2 stages.28 Stage I involves the initial cellular and enzymatic differentiation of the mammary epithelial cells resulting in the development of the capacity to synthesize and secrete milk during late pregnancy. The initial stage of lactogenesis also coincides in part with the period of colostrum accumulation. During stage I there is an increase in the cellular organelles associated with synthesis of milk components (eg, Golgi apparatus, ribosomes, endoplasmic reticulum; Fig. 8). The presence of these cellular organelles in epithelial cells from the porcine mammary gland at 105 days of gestation, but not 90 days, suggests that the start of the stage I occurs between those times.29 This is followed by a further increase in lipid droplets and endoplasmic reticulum in the epithelial cells between 105 and 112 days of gestation.

The concentration of total RNA in mammary tissue of pregnant gilts, while remaining relatively constant up to 90 days of gestation, then increases rapidly by 105 days and continues to increase at least through parturition.6 Furthermore, the rate of incorporation of acetate into lipid in mammary tissue cultures, a reflection of fatty acid synthesis, increases between 90 and 105 days, and then continues to increase at least to 112 days.6 Collectively, these observations are consistent with stage I of lactogenesis starting between 90 and 105 days of gestation in the sow and continuing up to parturition. Stage II of lactogenesis, referred to as copious milk secretion, occurs starting at or immediately after parturition in the sow. Changes in composition of mammary secretions from a few hours prior to parturition through several days postpartum confirm that stage II of lactogenesis does not occur until parturition in this species.30,31

In the present study, gene expression profiles in mammary tissue collected at 80 days of gestation reflect a growing gland that has not yet started stage I of lactogenesis, while at 100 days the tissue has initiated stage I, and at 110 days of gestation the tissue has fully entered into stage I. The precise timing of the start of stage I of lactogenesis in swine remains to be determined. Mechanistically, our analysis suggests that the marked increase in lipid droplet synthesis was driven by the sustained activation of ‘fatty acid biosynthesis’, ‘Tricarboxylic acid cycle (TCA) cycle’, and ‘glyoxylate and dicarboxylate’ flux, with a modest activation of ‘oxidative phosphorylation’ (Additional Files 1–2). Those results agree with quantitative measurements of acetate and glucose use as energy sources and also for synthesis of lipid during 75 through 105 days of pregnancy,6 ie, decrease in oxidation of both compounds coupled with increased use for lipogenesis particularly acetate.

Immune system

Both hormonal and local mammary gland factors are involved in control of immunoglobulin transport in the ruminant mammary gland.32 The induction of protein export observed at 110 vs. 80 and 110 vs. 100 days (Fig. 8) suggests that genes responsive to those hormones also might be part of the mechanism that enhance immune-related protein accumulation before farrowing. In the context of immune function and mammary development, the dramatic activation of ‘IL-4 receptor binding’ and ‘Positive regulation of isotype switching to IgE isotypes’ at 110 vs. 80 compared with 100 vs. 80 days (Fig. 8) was driven by the marked upregulation of IL-4 observed at 110 vs. 80 and 110 vs. 100 days (Tables 2 and 3). This cytokine is produced primarily by lymphocytes,33 suggesting that those cells are part of the mammary tissue “microenvironment” during pregnancy.

Table 2

List of DEG (FDR < 0.15 and raw P < 0.10) with 2-fold or greater upregulation on day 110 vs. 80. The Entrez Gene ID for each gene can be found in Additional File 1.

Gene symbol	Description	Fold
ABCB11	ATP-binding cassette, sub-family B (MDR/TAP), member 11	9.75
SLX1B	SLX1 structure-specific endonuclease subunit homolog B (S. cerevisiae)	6.02
DEAF1	Deformed epidermal autoregulatory factor 1 (Drosophila)	5.82
ZNF362	FLJ25476 protein	5.63
NDUFS5	NADH dehydrogenase (ubiquinone) Fe-S protein 5	5.59
MTHFR	5,10-methylenetetrahydrofolate reductase (NADPH)	5.56
CHERP	Calcium homeostasis endoplasmic reticulum protein	5.37
LALBA	Lactalbumin, alpha	5.31
KLHDC8A	Kelch domain containing 8A	5.09
CXXC1	CXXC finger 1 (PHD domain)	4.85
SEPT9	Septin 9	4.60
RWDD2A	RWD domain containing 2	4.54
IL4	Interleukin 4	4.39
ACO2	Aconitase 2, mitochondrial	4.37
APBB3	Amyloid beta (A4) precursor protein-binding, family B, member 3	4.21
FBXW9	F-box and WD-40 domain protein 9	4.18
ACSBG2	Acyl-CoA synthetase bubblegum family member 2	4.08
KLHDC8B	Kelch domain containing 8B	3.96
HPR	Haptoglobin-related protein	3.94
ZNF474	Zinc finger protein 474	3.85
GPD1	Glycerol-3-phosphate dehydrogenase 1 (soluble)	3.64
PIGR	Polymeric immunoglobulin receptor	3.47
LRRC8A	Leucine rich repeat containing 8 family, member A	3.38
MBD5	Methyl-CpG binding domain protein 5	3.36
RAD54L	RAD54-like (S. cerevisiae)	3.36
LSP1	Lymphocyte-specific protein 1	3.33
SIX6	Sine oculis homeobox homolog 6 (Drosophila)	3.22
RAB2A	RAB2, member RAS oncogene family	3.06
HDAC6	Histone deacetylase 6	3.04
KRT17	Keratin 17	3.00
CRELD2	Cysteine-rich with EGF-like domains 2	2.99
CSNK1A1	Casein kinase 1, alpha 1	2.94
CCDC159	Coiled-coil domain containing 159	2.90
HHIPL2	HHIP-like 2	2.86
MARCH3	Membrane-associated ring finger (C3HC4) 3	2.85
RNLS	Renalase, FAD-dependent amine oxidase	2.81
AP2B1	Adaptor-related protein complex 2, beta 1 subunit	2.73
SRPRB	Signal recognition particle receptor, B subunit	2.69
SSNA1	Sjogren’s syndrome nuclear autoantigen 1	2.66
BTN1A1	Butyrophilin, subfamily 1, member A1	2.65
RNF11	Ring finger protein 11	2.63
ITGB4	Integrin, beta 4	2.61
MAP2K5	Mitogen-activated protein kinase kinase 5	2.57
FAM13C	Family with sequence similarity 13, member C	2.56
KLF9	Kruppel-like factor 9	2.45
C20orf24	Chromosome 20 open reading frame 24	2.39
SPTBN1	Spectrin, beta, non-erythrocytic 1	2.37
RARB	Retinoic acid receptor, beta	2.36
COL8A1	Collagen, type VIII, alpha 1	2.28
TARP	TCR gamma alternate reading frame protein	2.26
ZNF197	Zinc finger protein 197	2.24
ARRDC3	Arrestin domain containing 3	2.24
EPS15	Epidermal growth factor receptor pathway substrate 15	2.24
MRPL11	Mitochondrial ribosomal protein L11	2.23
PAX8	Paired box 8	2.22
IFI44L	Interferon-induced protein 44-like	2.20
C1orf210	Chromosome 1 open reading frame 210	2.16
ZSWIM5	Zinc finger, SWIM-type containing 5	2.11
FBXO41	F-box protein 410	2.10
SNX29	Sorting nexin 29	2.10
ACY1	Aminoacylase 1	2.08
TMEM194A	Transmembrane protein 194A	2.07
ABCG2	ATP-binding cassette, sub-family G (WHITE), member 2	2.06
PIGF	Phosphatidylinositol glycan anchor biosynthesis, class F	2.05
VIPR1	Vasoactive intestinal peptide receptor 1	2.05
PLK5	Polo-like kinase 5	2.03
CD86	CD86 molecule	2.03
RBM6	RNA binding motif protein 6	2.02
FBP1	Fructose-1,6-bisphosphatase 1	2.00

Table 3

List of DEG (FDR < 0.15 and raw P < 0.10) with 2-fold or greater upregulation on day 110 vs. 100. The Entrez Gene ID for each gene can be found in Additional File 1.

Gene symbol	Description	Fold
BARD1	BRCA1 associated RING domain 1	5.35
SLX1B	SLX1 structure-specific endonuclease subunit homolog B	5.26
STRADA	STE20-related kinase adaptor alpha	5.12
KLHDC8A	Kelch domain containing 8A	5.12
CHERP	Calcium homeostasis endoplasmic reticulum protein	4.26
DEAF1	Deformed epidermal autoregulatory factor 1 (Drosophila)	4.22
UROC1	Urocanase domain containing 1	4.02
MTHFR	5,10-methylenetetrahydrofolate reductase (NADPH)	3.70
FLYWCH1	FLYWCH-type zinc finger 1	3.66
ABCB11	ATP-binding cassette, sub-family B (MDR/TAP), member 11	3.51
GPD1	Glycerol-3-phosphate dehydrogenase 1 (soluble)	3.37
KLHDC8B	Kelch domain containing 8B	3.31
LALBA	Lactalbumin, alpha-	3.01
SHARPIN	SHANK-associated RH domain interactor	2.98
ZNF474	Zinc finger protein 474	2.93
MRI1	Methylthioribose-1-phosphate isomerase homolog	2.86
ACO2	Aconitase 2, mitochondrial	2.77
ITGB4	Integrin, beta 4	2.72
MST1	Macrophage stimulating 1 (hepatocyte growth factor-like)	2.64
TGM3	Transglutaminase 3 (protein-glutamine-gamma-glutamyltransferase)	2.61
IL4	Interleukin 4	2.58
CDC37	Cell division cycle 37 homolog (S. cerevisiae)	2.52
MBD5	Methyl-CpG binding domain protein 5	2.50
MAP2K5	Mitogen-activated protein kinase kinase 5	2.45
RAB2A	RAB2, member RAS oncogene family	2.38
SPPL3	Signal peptide peptidase like 3	2.36
EPS15	Epidermal growth factor receptor pathway substrate 15	2.34
ARHGAP17	Rho GTPase activating protein 17	2.34
SEZ6L	Seizure related 6 homolog (mouse)-like	2.33
SNORD88C	Small nucleolar RNA, C/D box 88C	2.28
ZNF197	Zinc finger protein 197	2.25
DEF6	Differentially expressed in FDCP 6 homolog (mouse)	2.22
SRPRB	Signal recognition particle receptor, B subunit	2.20
NOMO1	NODAL modulator 1	2.17
SSNA1	Sjogren’s syndrome nuclear autoantigen 1	2.13
CSNK1A1	Casein kinase 1, alpha 1	2.03
AP2B1	Adaptor-related protein complex 2, beta 1 subunit	2.03
DNAJB12	DnaJ (Hsp40) homolog, subfamily B, member 12	2.02

Although we are unaware of similar measurements in livestock species, a recent study revealed that IL-4 concentration in lactating mice was very low relative to other cytokines such as TNF, several interleukins, and several chemokines.34 Because it has been established that IL-4 induces epithelial cell proliferation in vitro via activating MAPK (including JNK),35 the high expression of this cytokine as well as the activation of the MAPK pathway (Fig. 3) in mammary tissue during late-pregnancy suggest that it could play a role in the late stages of mammary development.

Although the concentration of IgE in milk is very low,36 another potential biological role for the upregulation of IL4 might be related with the synthesis of IgE by mononuclear cells37 that may form part of the mammary microenvironment. Synthesis of IgE requires T and B cell contact. Naïve B cells undergo class switching leading to the development of memory B cells and plasma cells that produce IgE.38 Antibodies for IgE are involved in the host defense against some types of infection.39 The observed responses seem to suggest the infiltration of IgE-producing memory B cells into the mammary gland at a point in which the functional differentiation of the mammary gland seems to be nearly complete. It remains to be determined if the marked activation of IgE isotype switching has a functional role, eg, to expedite the re-encounter with antigen.40

From a functional standpoint the high impact and marked activation of ‘IL-4 receptor binding’ seems to be mechanistically related with IgE switching, because IL-4 is a key regulator of that process and also of B-cell proliferation.41 It might be possible that the marked activation of the immune system observed close to farrowing serves to provide protection to the mammary gland from potential pathogens.

Metabolism

Our finding of ‘lipid’ metabolism as being highly-impacted on day 100 vs. 80 and 110 vs. 80 (Fig. 2) and the marked upregulation (Additional File 1, Tables 1 and 2) of VLDLR, FABP1, SLC2A6, BTN1A1, ACACA, ABC G2, IDH2, ACO2, and FBP1 underscored the degree of coordination of several key pathways that would allow for use of long-chain fatty acids from the circulation (VLDLR, FABP1) but also lipogenic substrates (SLC2A6) by mammary cells. The combined data suggested a greater flux via the TCA cycle (ACO2, IDH2), lipogenesis via ACACA, glyceroneogenesis to provide glycerol-3-phosphate (FBP1), lipid-droplet associated proteins (BTN1A1) to complete the formation of milk fat, and lastly the transport of lipid droplets into the gland cistern (ABCG2).42–45 Although not all these genes are targets of PPARG (eg, ACACA, FABP1, ABCG2), its upregulation (Figs. 5 and 6; Additional File 1) at 100 and 110 days likely elicited a degree of control over the milk fat synthesis process.46 Together, these data agree with biochemical studies performed several years ago and outlined in the introductory section.

Although we are unaware of quantitative measures of amino acid (AA) metabolism in mammary tissue during late pregnancy, our data indicated that the inhibition of AA metabolism observed at 110 vs. 80 relative to 100 vs. 80 days was associated with reduced degradation of Essential Amino Acids (EAA) such as Lys, Val, Leu, and Ile, and also reduced ‘Glutathione’, ‘Propanoate’, ‘Selenoamino acid’, ‘Arginine and Proline’, and ‘Glycine, Serine, and Threonine’ metabolism (see Additional File 2: Sheet ‘DIA KEGG’ for details; see also Fig. 8). From a functional standpoint, an inhibition of ‘Glutathione metabolism’ would allow for sparing use of NADPH for fatty acid synthesis and re-synthesis of methionine from 5,10-methylene tetrahydrofolate via MTHFR (Tables 2 and 3), ie, both pathways use NADPH.47,48 More importantly, because AA can feed carbon into the TCA cycle or can be used for tissue (parenchyma) and colostrum protein synthesis, the observed adaptations could indicate that the swine mammary gland shifts the priority of AA utilization from greater use/catabolism in early pregnancy, eg, to sustain parenchymal growth, to sparing AA catabolism near parturition either for colostrum synthesis or fetal utilization.3,6 The inhibition of ‘Propanoate metabolism’ at 110 vs. 100 days relative to other time comparisons, in particular, supports the view of a gradual decrease in use of AA for catabolism. This pathway encompasses the oxidative deamination of glutamate to alpha-ketoglutarate and ammonia via the key enzyme GLUD1, leading to an increased flux through the TCA cycle and an overall increase in ATP production.47

The estimated daily gain of fetal crude protein is 4.6 g/fetus/day between day 70 of gestation and farrowing, during which the rate of crude protein deposition in each mammary gland is 3.4 g/day.4 In fact, the use of EAA for both fetal and mammary tissue crude protein deposition increases dramatically between day 80 of pregnancy to farrowing. Such response is particularly evident for Leu, Lys, Val, Thr, Ile, and Arg,49 and could partly explain the activation of ‘Valine, Leucine, and Isoleucine’ and ‘Lysine’ degradation at 100 vs. 80 days relative to 110 vs. 100 days (see Additional File 2: Sheet ‘DIA KEGG’ for details; also Additional File 3 for pathway details). Thus, in the context of our study, the marked transcriptional changes observed between 100 and 110 vs. 80 days provide a window into the mechanisms that drive the functional differentiation of mammary gland from the non-lactating to the lactating state.

A novel finding of our study was the marked activation of ‘Galactose metabolism’ at 100 vs. 80 days, a response driven primarily by the marked upregulation of UGP2 (UDP-glucose pyrophosphorylase) and AKR1B1 (aldo-reductase) (Tables 1 and 2) coupled with the downregulation of LCT (lactase) and B4GALT1 (Additional Files 1 and 3). Based on the evaluation of the KEGG pathway our data would suggest greater UDP-glucose and D-galactose production but lower D-lactose (Additional File 3). The ultimate fate of D-galactose, particularly at 100 vs. 80 days, could have been its metabolism through glycolysis as a source of energy. However, the marked upregulation of LALBA (Tables 2 and 3) at 110 vs. 80 and 100 days and the downregulation of UGP2 (namely at 110 vs. 100 days) likely resulted in greater lactose synthesis. Although not measured in earlier studies, it is likely that the observed accumulation of secretory products in alveoli of mammary tissue at 105 days of gestation also is driven by lactose synthesis such that by ~110 to ~112 days of parturition expression of genes such as LALBA is up-regulated. The inhibition of ‘Galactose’, ‘Starch and sucrose’, and ‘Glyoxylate and dicarboxylate’ metabolism at 110 vs. 100 days compared with 110 vs. 80 days might have been associated with the insulin-resistant state of pregnancy to support the marked rate of fetal growth at the latter stages of pregnancy.50

The fact that there was a decrease in the impact and a lack of change in flux for most metabolic pathways at 110 vs. 100 days (Fig. 1 and Additional File 2: Sheet ‘KEGG Summary’) seems to suggest that other molecular adaptations take precedence during the immediate time preceding parturition, ie, help “complete” the transition of the tissue to “fully-functional” gland or from stage I to stage 2 of lactogenesis. That was particularly evident when exploring the overall impact and flux of most pathways associated with GIP, EIP, CP, and OS in the KEGG database.

Cellular communication

The mammary gland experiences repeated cycles of pregnancy, lactation, and involution. During these cycles, the interactions of mammary epithelial cells with their surrounding extracellular matrix (ECM), particularly the basement membrane and the connective tissue, play a pivotal role in building a communication bridge between epithelial cells and their surrounding environment.51,52 In rat cervical stromal cells it was demonstrated that ITGA11 and ITGA2 expression increases through most of the pregnancy and decreases along with progesterone as the time of parturition nears.53 The functions of ITGA11 include fluid homeostasis54 by maintaining ECM tension and inhibition of cell migration;55 hence, upregulation of ITGA11 expression during pregnancy may help prevent excessive tissue remodeling, whereas the decrease close to parturition may serve to enhance influx of inflammatory cells and fluid.53 The upregulation of collagen genes (Table 1), vitronectin (VTN), and ITGA11 (Additional File 1) and the lack of change in laminin at 100 vs. 80 days (Additional File 3) seem to agree with the data from pregnant rats. However, at 110 vs. 80 days the downregulation of collagen (Tables 4 and 5; Additional File 3) coupled with upregulation of ITGA11 and ITGB4 (Additional File 1) suggests that these integrins may play a different role in the mammary gland. For instance, they may protect against excessive tissue remodeling via the upregulation of matrix metallopeptidases (eg, MMP1; Additional File 1), which are involved in collagen digestion and breakdown of the extracellular matrix.56

Table 4

List of DEG (FDR < 0.15 and raw P < 0.10) with 2-fold or greater downregulation on day 100 vs. 80.

Gene symbol	Description	Fold
ZNF692	Zinc finger protein 692	−2.82
DSP	Desmoplakin	−2.47
SOX4	SRY (sex determining region Y)-box 4	−2.45
STMN1	Stathmin 1/oncoprotein 18	−2.35
RANBP3	RAN binding protein 3	−2.35
YIF1A	Yip1 interacting factor homolog A (S. cerevisiae)	−2.33
COL5A2	Collagen, type V, alpha 2	−2.31
GNAQ	Guanine nucleotide binding protein (G protein), q polypeptide	−2.30
METAP2	Methionyl aminopeptidase 2	−2.30
CENPF	Centromere protein F, 350/400ka (mitosin)	−2.27
POSTN	Periostin, osteoblast specific factor	−2.27
TRGJP	TRGJP T cell receptor gamma joining P	−2.25
H2AFZ	H2A histone family, member Z	−2.25
MDK	Midkine (neurite growth-promoting factor 2)	−2.20
PTPN11	Protein tyrosine phosphatase, non-receptor type 11	−2.18
PADI2	Peptidyl arginine deiminase, type II	−2.18
RCOR1	REST corepressor 1	−2.15
MXRA7	Matrix-remodelling associated 7	−2.15
CYP2C18	Cytochrome P450, family 2, subfamily C, polypeptide 18	−2.11
PRRC2B	Proline-rich coiled-coil 2B	−2.10
PDLIM5	PDZ and LIM domain 5	−2.09
TUBA1C	Tubulin, alpha 6	−2.08
ZCCHC3	Zinc finger, CCHC domain containing 3	−2.08
NAB1	NGFI-A binding protein 1 (EGR1 binding protein 1)	−2.07
NPM1	Nucleophosmin (nucleolar phosphoprotein B23, numatrin)	−2.05
G2E3	G2/M-phase specific E3 ubiquitin protein ligase	−2.05
EGR1	Early growth response 1	−2.04
CDC45	Cell division cycle 45	−2.03
RPL15	Ribosomal protein L15	−2.03
MAP1B	Transcribed locus	−2.02
SS18	Synovial sarcoma translocation, chromosome 18	−2.02
ATRX	Alpha thalassemia/mental retardation syndrome X-linked	−2.01
ZNF136	Zinc finger protein 44	−2.01
PPP1R12A	Protein phosphatase 1, regulatory (inhibitor) subunit 12A	−2.01
TUBA1C	Tubulin, alpha 6	−2.01
NREP	Neuronal regeneration related protein	−2.00

Table 5

List of DEG (FDR < 0.15 and raw P < 0.10) with 2-fold or greater downregulation on day 110 vs. 80. The Entrez Gene ID for each gene can be found in Additional File 1.

Gene symbol	Description	Fold
COL1A2	Collagen, type I, alpha 2	−3.57
COL1A1	Collagen, type I, alpha 1	−3.55
COL1A2	Collagen, type I, alpha 2	−3.15
RANBP3	RAN binding protein 3	−3.09
POSTN	Periostin, osteoblast specific factor	−2.93
COL3A1	Collagen, type III, alpha 1	−2.90
ZNF692	Zinc finger protein 692	−2.86
DSP	Desmoplakin	−2.67
SGSM3	RUN and TBC1 domain containing 3	−2.63
SOX4	SRY (sex determining region Y)-box 4	−2.58
H2AFZ	H2A histone family, member Z	−2.49
EGR1	Early growth response 1	−2.49
GBP1P1	Guanylate binding protein 1, interferon-inducible pseudogene 1	−2.45
H3F3B	H3 histone, family 3B (H3.3B)	−2.45
ZFP36L1	Zinc finger protein 36, C3H type-like 1	−2.45
OR1M1	Olfactory receptor, family 1, subfamily M, member 1	−2.44
C1orf55	Chromosome 1 open reading frame 55	−2.40
UNC13C	Unc-13 homolog C (C. elegans)	−2.36
NXNL2	Chromosome 9 open reading frame 121	−2.35
TUBA1C	Tubulin, alpha 6	−2.34
OXTR	Oxytocin receptor	−2.33
C1GALT1	Core 1 synthase, glycoprotein-N-acetylgalactosamine 3-beta-galactosyltransferase, 1	−2.29
RPL15	Ribosomal protein L15	−2.28
ELK3	ELK3, ETS-domain protein	−2.28
LAPTM4B	Lysosomal associated protein transmembrane 4 beta	−2.26
PADI2	Peptidyl arginine deiminase, type II	−2.25
PPP4R1	Protein phosphatase 4, regulatory subunit 1	−2.23
COL5A2	Collagen, type V, alpha 2	−2.23
HNRNPH2	Heterogeneous nuclear ribonucleoprotein H2	−2.21
ZFP36L1	Zinc finger protein 36, C3H type-like 1	−2.20
METAP2	Methionyl aminopeptidase 2	−2.18
ARL6IP6	ADP-ribosylation-like factor 6 interacting protein 6	−2.16
HNRNPR	Heterogeneous nuclear ribonucleoprotein R	−2.14
RPL39	Ribosomal protein L39	−2.13
TUBA1C	Tubulin, alpha 6	−2.13
FKBP3	FK506 binding protein 3, 25 kDa	−2.12
SMTNL2	Smoothelin-like 2	−2.12
CALM2	Calmodulin 2 (phosphorylase kinase, delta)	−2.09
PTPRF	Protein tyrosine phosphatase, receptor type, F	−2.08
C7orf60	Chromosome 7 open reading frame 60	−2.08
RPL17	Ribosomal protein L17	−2.06
TK2	Thymidine kinase 2, mitochondrial	−2.05
CEP135	Centrosomal protein 135 kDa	−2.04
PARP3	Poly (ADP-ribose) polymerase family, member 3	−2.03
HNRPDL	Heterogeneous nuclear ribonucleoprotein D-like	−2.01
HLA-A	Major histocompatibility complex, class I, A	−2.00

Junctional complexes including tight junctions (TJs), adherens junctions (AJs), and gap junctions serve an important role in mediating cell–cell interactions.57 The importance of gap junctions in the proper development, differentiation and proliferation of the mammary gland has been known for some time.58–61 Within gap junctions there are intercellular channels that allow direct communication between cells, allowing for direct transport of small molecules such as ions, amino acids, and other metabolites.62 These constitute normal physiological processes between adjacent cells. Thus, the marked increase in activation of gap junction (Fig. 8; Additional File 1) is in accordance to well-established knowledge of the underlying physiological events.63

Mammary gland tight junctions are permeable but experience closure around farrowing in the pregnant animal.64 In turn, the closure of tight junctions around parturition is often accompanied by a decrease in some ions such as sodium and chloride, and an increase in lactose and potassium in milk.65 The drop in progesterone prior to parturition is also responsible for the closure of mammary tight junctions in the late-pregnant mouse.66 Thus, the observation of a modest activation of genes associated with tight junctions at 110 vs. 100 days (Fig. 8) supports the role of progesterone during late-pregnancy.

Signaling regulation

The induction of Gonadotropin-releasing hormone (GnRH) signaling during pregnancy (Fig. 3) is similar to what was observed recently in bovine mammary gland, and suggests this novel mechanism might be conserved among mammals. Expression of GnRH receptor (GnRHR), which allows for GnRH signaling and synthesis of gonadotropins in the anterior pituitary, also occurs in the mammary gland.67 The function of GnRHR in extrapituitary tissues is not clear, but it could play a role in reducing cell proliferation.67 Many hormones have been detected in milk but not Leutenizing hormone (LH) or follicle stimulating hormone (FSH).68 Further studies are required to clarify these findings.

Previous in vivo data from the bovine mammary gland23,69 clearly indicated that mTOR plays an important role in milk protein synthesis. The marked induction of the mTOR signaling pathway (Fig. 3) was unexpected because this pathway changed little during late-pregnancy in the bovine mammary gland23 and, in fact, there was a small degree of inhibition of this pathway during early lactation. In a recent microarray study of the swine mammary gland, a pathway analysis indicated that few genes composing the mTOR pathway were upregulated on day 17 compared with day −17 around parturition.11 Although at a different stage of the lactation cycle (ie, non-lactating mammary), the present data on ribosome, Golgi apparatus, amino acid metabolism (eg, His, Lys), and mTOR signaling (Fig. 3) suggest an important role of this pathway in preparing the mammary gland for active milk protein synthesis.69

Blood flow regulation

An increase in mammary blood flow is observed concurrently with the fall in progesterone close to parturition, and is thought to be a physiological adaptation to increase volume of secretions within the gland.70 At a more basic level, the control of blood volume through vessels is closely regulated via contraction and relaxation of vascular smooth muscle cells (VSMC) at least in part via signaling through adrenergic receptors or intracellular Ca2− concentration.71,72 The fact that vascular smooth muscle contraction was activated at 100 and 110 vs. 80 days (Supplementary File 2, sheet DIA KEGG) seems to argue against a role of Ca2− signaling. It is also possible that estrogen and relaxin might play an inhibitory role on VSMC function (ie, proliferation and contraction).73,74

There was little change observed in the angiogenesis pathway (Fig. 8), and the increase in activation of Ca2− signaling observed at 110 vs. 100 days was modest compared with the marked increase in VEGF and MAPK signaling (Fig. 3). In that context, the fact that ‘positive regulation of endothelial cell differentiation’ was substantially inhibited (Additional File 2, sheet GOTERM BP) and ‘endothelial cell migration’ was markedly activated at 110 vs. 80 days compared with 110 vs. 100 days seems to suggest that the effect of VEGF signaling on VSMC of the mammary gland during the late stages of pre-partum mammary development is to promote endothelial cell migration rather than formation of new blood vessels.75 The marked increase in the activation of MAPK signaling76 also supports a functional role of VEGF in the formation and migration of endothelial cells within the mammary gland.

Epigenetic regulation and transcription

As opposed to what was recently observed in mature bovine mammary gland,23 our data revealed a marked and sustained activation of the components of epigenetic modification (euchromatin and nuclear euchromatin; Fig. 8) during late-pregnancy. Those data, however, contrast with the lower-impact and activation of the GO terms ‘chromatin assembly or disassembly’ and ‘chromatin remodeling’ on day 100 and 110 along with the marked inhibition of ‘chromatin-mediated maintenance of transcription’ on day 110 vs. 80 (Additional File 2, sheet GOTERM BP). The pattern of activation of euchromatin at 100 and 110 vs. 80 days was in line with the small change in the total number of DEG (Fig. 1) but the large increase in tissue DNA content (Fig. 8) at both time points. Interestingly, there were fewer DEG at 110 vs. 100 days, which corresponds to the lack of change in tissue DNA content between both times, and potentially with the fact that the GO term ‘chromatin-mediated maintenance of transcription’ was markedly inhibited between 100 to 110 vs. 80 days and was not impacted at all on day 110 vs. 100.

The above changes, along with the marked and progressive activation of the GO term ‘positive regulation of transcription factor activity’, suggest an active degree of transcription through at least 100 days of pregnancy. By 110 days, however, the continued development of the mammary gland as discerned by the increase in cross sectional area (Fig. 8) and previous physiological data6 would have been driven at least in part by activation of genes associated with DNA replication, cell cycle, and cell growth (Figs. 3 and 8). Overall, results highlight a seemingly crucial role of epigenetic regulation and transcription factor activity in swine mammary development during pregnancy. Such a role has been reported to exist in the lactating bovine mammary gland.77

Gene networks

The consistent downregulation of the transcription regulators TP53, ARNT2, EGR1, HDAC1, and E2F4 coupled with the upregulation of PPARG, EPAS1, and KLF9 on day 100 and/or 110 compared with day 80 underscored their potential importance in the local control of swine mammary development at the latter stages of pregnancy (Figs. 5 and 6). Furthermore, from a functional standpoint the upregulation at day 110 vs. 100 of the transcription factors GABPA, which controls expression of oxytocin receptor (OXTR) and prolactin (PRL), and STAT5B, which controls milk protein synthesis, appear quite important at the closest point to parturition that we studied (Fig. 6). Although we are unaware of similar data in pig mammary gland, previous work in bovine indicated that TP53 and HDAC1 are primarily expressed in parenchymal than fat pad tissue, but PPARG, KLF9, EGR1, EPAS1 are primarily expressed in fat pad.26 It also can be deduced from mammary tissue studies in dairy cattle26,78 that ARNT2, NR3C1, and ETV1 are likely preferentially expressed in parenchyma as opposed to the fat pad.

The apparent role of TP53 networks in mammary gland development is not entirely surprising because this transcription regulator has been well-studied in the context of cancer, where its inhibition leads to cell cycle progression and eventually neo-plastic transformation.79 Furthermore, E2F4 plays an important role in the control of cell cycle and the action of tumor suppressor proteins, causing an overall decrease in expression of proliferation-associated genes.80 It was also recently demonstrated that nuclear E2F4, when unregulated, plays a key role in the downregulation of several mitotic genes, hence, promoting cell cycle arrest.81 Therefore, in the context of mammary development, the downregulation of TP53 and E2F4 along with the observed levels of expression of its target genes could be functionally related to cellular proliferation (eg, greater activation of the cell cycle pathway; Figs. 5 and 6) and the morphological changes that were observed between day 80 and 110 (Fig. 8).

A recent study clearly established that the role of TP53 in enhancing cell-cycle arrest is mediated at the transcriptional level.82 Whereas the role of TP53 appears clear at both 100 and 110 days, the activation of several signaling pathways known to promote cell differentiation (eg, Jak-STAT, Hedgehog signaling; Fig. 3) and cell growth (mTOR signaling) were markedly activated at 110 vs. 100 days suggesting they may play a more important role in the latter stages of pregnancy and in the case of mTOR signaling may even be important for the complete activation of the milk protein synthesis machinery.69 Further support for the functional role of TP53 at the transcription level is the recent demonstration that its activation results in the inhibition of the mTOR pathway and consequently a repression of protein translation.82

Although much less is known about the functional role in mammary development of several transcription regulators that we uncovered, eg, DEAF1, MECP2, RARB, EGR1, ETV1, and EPAS1, the fact that they experienced marked upregulation or downregulation depending on day of pregnancy suggests that they are important. For instance, recent studies revealed that 17β-estradiol increased the expression of DEAF1 in vitro,83 which in turn could downregulate serotonin signaling84 leading to a potential impairment in milk secretion.85 The fact that DEAF1 was markedly upregulated (~6-fold) at 110 vs. 80 days supports its role in the process of lactogenesis. EPAS1 is a transcription factor that induces genes regulated by oxygen, ie, it is induced when cellular oxygen levels fall, and, for instance, it has been demonstrated to play a role in the proliferation-promoting effect of HIF2A, partly through enhancing mTOR signaling.86 The role of RARB in limiting cell growth by regulating gene expression is well-established;87 thus, it could be envisioned that it forms part of the network controlling mammary cell proliferation.

A novel and seemingly important finding in this study was the marked upregulation of MECP2 on day 100 and 110 vs. 80, which encodes a nuclear protein that possesses a methyl-CpG binding domain; thus, it is capable of binding specifically to methylated DNA.88 This protein functions primarily as a transcriptional repressor. Although the 5-hydroxylmethylcytosine (5 hmC) and 5-methylcytosine content in the mammary gland genome is unknown, it has been shown in neuronal cells containing high levels of 5 hmC that MECP2 is the major methyl-CpG-binding protein.89 Because mutations that inhibit 5 mhC binding are known to cause disease, it is tempting to speculate that in the pregnant swine mammary gland of primiparous animals this nuclear protein may serve as one of the epigenetic mechanisms for regulating euchromatin structure and gene expression. A recent review discussed the role of epigenetics in regulation of bovine mammary function.90

Conclusions

The present bioinformatics and gene network functional analysis of the swine mammary transcriptome strengthens the understanding of the biological adaptations occurring in this organ during late gestation. The data underscored a novel role for IL-4 in mammary development during pregnancy. Several elements across a wide number of pathways were identified as potential regulatory factors of the local adaptations in mammary gland. Further research will help confirm the functional relevance of the pathways uncovered, and how they can be manipulated to control sow milk yield and quality.

Additional File 1

Additional File 2

Additional File 3

Additional File 4