Gene expression in the mammary gland of the tammar wallaby during the lactation cycle reveals conserved mechanisms regulating mammalian lactation.

CSIRO PUBLISHING

Reproduction, Fertility and Development http://dx.doi.org/10.1071/RD14210

Gene expression in the mammary gland of the tammar wallaby during the lactation cycle reveals conserved mechanisms regulating mammalian lactation C. J. Vander Jagt A,B,D, J. C. Whitley B, B. G. Cocks B,C and M. E. Goddard A,B A

Department of Agriculture and Food Systems, Melbourne School of Land and Environment, The University of Melbourne, Parkville, Vic. 3010, Australia. B Computational Biology, Department of Environment and Primary Industries, AgriBio, Centre for AgriBioscience, 5 Ring Road, Bundoora, Vic. 3083, Australia. C School of Applied Systems Biology, La Trobe University, Bundoora, Vic. 3083, Australia. D Corresponding author. Email: [email protected]

Abstract. The tammar wallaby (Macropus eugenii), an Australian marsupial, has evolved a different lactation strategy compared with eutherian mammals, making it a valuable comparative model for lactation studies. The tammar mammary gland was investigated for changes in gene expression during key stages of the lactation cycle using microarrays. Differentially regulated genes were identified, annotated and subsequent gene ontologies, pathways and molecular networks analysed. Major milk-protein gene expression changes during lactation were in accord with changes in milkprotein secretion. However, other gene expression changes included changes in genes affecting mRNA stability, hormone and cytokine signalling and genes for transport and metabolism of amino acids and lipids. Some genes with large changes in expression have poorly known roles in lactation. For instance, SIM2 was upregulated at lactation initiation and may inhibit proliferation and involution of mammary epithelial cells, while FUT8 was upregulated in Phase 3 of lactation and may support the large increase in milk volume that occurs at this point in the lactation cycle. This pattern of regulation has not previously been reported and suggests that these genes may play a crucial regulatory role in marsupial milk production and are likely to play a related role in other mammals. Additional keywords: DAVID, EST, gene ontology, IPA, microarray, milk, network.

Received 16 June 2014, accepted 21 December 2014, published online 23 February 2015

Introduction Lactation is the key feature distinguishing mammals from other animals. Milk is essential for the early development of all mammalian newborn and infants. It not only provides the perfect natural source of nutrition for the young to grow and thrive, numerous protective factors in milk offer immunologic protection from infectious disease and may play a role in immune system development (Hasselbalch et al. 1996; Thompson et al. 2000). Despite being separated by at least 130 million years of evolution (Luo et al. 2003; Bininda-Emonds et al. 2007; Nilsson et al. 2010), the process of lactation in both eutherians (such as the cow) and marsupials (such as the tammar wallaby) produces the same outcome: milk. In eutherians, lactation is divided into two phases: preparation of the mammary gland during pregnancy (mammogenesis) and milk synthesis and secretion after birth (lactogenesis). Eutherians have adopted a strategy with a relatively short lactation, producing milk more or less uniform in composition after the birth of relatively mature young (Jenness 1986; Akers 2002). In contrast, marsupials have evolved a Journal compilation Ó CSIRO 2015

reproductive approach with a short gestation period, giving birth to fetus-like young that require milk for an extended period of time. Milk produced by the lactating mother alters considerably both in terms of volume and composition throughout the lactation cycle in order to provide the appropriate nutrition for the developing young (Green 1984; Nicholas 1988). This change in milk volume and composition is mediated by local factors in the mammary gland because it is possible for a tammar wallaby, for instance, to produce milk of different compositions from each gland simultaneously (Nicholas 1988). This suggests that the changes in milk secretion over the lactation cycle are mediated, at least in part, by changes in gene expression within the mammary gland. The aim of this paper is to report the changes in gene expression in the mammary gland of the tammar wallaby over the lactation cycle and to consider how these drive the observed changes in milk volume and composition. It is now recognised that marsupials have three lactation phases. In Phase 1, like eutherians, the mammary gland develops the capacity for milk synthesis and secretion. Phase 2, unique to www.publish.csiro.au/journals/rfd

B

Reproduction, Fertility and Development

marsupials, is the period from birth until approximately the time when the young vacates the maternal pouch. In the tammar wallaby, Phase 2 encompasses the first 200 lactation days and is further subdivided into Phases 2A and 2B (thus it is often said that there are four phases in tammar lactation). Milk secreted during this time is high in carbohydrate, low in fat and the rate of milk production is limited. Phase 3, equivalent to the whole of lactation in eutherians, is the period of extensive growth in the young. In the tammar, Phase 3 is characterised by copious secretion of lipid- and protein-rich milk and typically lasts till Day 300 (Messer and Green 1979; Nicholas 1988). Therefore, it is possible to observe larger changes in milk synthesis in the wallaby than in eutherian mammals and this increases our power to discover the factors driving these changes. To examine gene expression in the mammary gland during the lactation cycle, microarray technology can be employed. While this technique has been used extensively to investigate transcriptional regulation in the mammary gland of the mouse (Master et al. 2002; Rudolph et al. 2003; Clarkson et al. 2004; Stein et al. 2004; Oakes et al. 2006; Lemay et al. 2007; Ron et al. 2007), to date only two studies using microarrays to examine gene expression changes in the tammar wallaby mammary gland during lactation have been reported (Menzies et al. 2009; Khalil et al. 2011). To investigate the genetics driving the increase in milk protein seen in Phase 3 of tammar wallaby lactation, Menzies et al. (2009) analysed changes in mammary gland gene expression between Phase 2A and Phase 3 of lactation. Using a cDNA microarray (Lefe`vre et al. 2007), they observed a 2.4-fold upregulation in the expressed sequence tag (EST) representing folate receptor a (FOLR1). They hypothesised that FOLR1 may be a key regulatory point of folate metabolism for milk-protein synthesis within mammary epithelial cells. More recently, Khalil et al. (2011) identified a number of genes encoding novel putative milk proteins upregulated during involution (when compared with Phase 2B and 3 of lactation). It was hypothesised that these milk proteins secreted in involution may play a protective role against mammary gland infection. They also reported the upregulation of the pro-apoptotic tumour necrosis factor receptor superfamily 21 (TNFRSF21) gene, whose expression in the mammary gland had not previously been reported. Both these studies demonstrated that the marsupial reproductive strategy can be successfully exploited to identify novel genes implicated in the regulation of lactation. In this study we used microarray technology to identify differentially expressed genes in the tammar wallaby mammary gland during the lactation cycle. We have focussed on gene expression changes in the tammar mammary gland as it shifts from one lactation phase to the next. Transcriptome data was subjected to bioinformatic analyses for further understanding of how milk production is regulated through the three tammar lactation phases (including between Phase 2A and 2B). Transcriptional trends were explored via a non-hypothesis-driven approach, with potential key regulatory genes identified and discussed. Gene products and their interactions were visualised as networks and findings compared with public mouse and cow lactation expression data. To our knowledge, this is the first microarray study to investigate gene expression across Phases 1–3 of the tammar lactation cycle.

C. J. Vander Jagt et al.

E1 P18

L2

E2 P26

L22

E3

L62

L151

L130

L216

L87

L171

L151

L243

L110

L193

L171

L266

Fig. 1. Tammar wallaby lactation time-course microarray experimental designs. Three separate time-course microarray experiments were performed (E1, E2 and E3). Each experiment involved the comparison of various time points across two phases of tammar lactation (represented by adjoining lines). RNA used was derived from 14 time points spanning the length of the tammar lactation cycle; from Day 18 of pregnancy (time-point P18), to Day 266 of lactation (time-point L266). Each time-point sample came from a single animal and each microarray experiment was performed with a dye-swap replicate.

Materials and methods Tammar mammary gland samples Mammary gland samples were harvested from tammar wallabies housed in a purpose-built facility at the Department of Primary Industries (DPI) Victoria, Attwood. All procedures involving the use of the tammar wallabies were approved by the DPI Animal Ethics Committee (approval number 2566) and all animals were euthanased prior to tissue collection. Mammary gland sample time points were selected that would best represent the progressive changes occurring in the mammary gland throughout the lactation cycle and included samples from: early pregnancy (proliferative stage, Day 18 of pregnancy), late pregnancy (secretory pre-parturition, Day 26 of pregnancy), secretory post-parturition (Days 2 and 22 of lactation), established lactation (Days 62, 87, 110, 130, 151, 171 and 193) and late lactation (Days 216, 243 and 266 of lactation). Total RNA was extracted from all mammary samples and the subsequent mRNA was isolated. In total, mRNA from 14 time points was produced. Microarray experiments Three separate time-series microarray experiments were conducted (E1, E2 and E3) using microarray slides with 13 440 tammar mammary-derived ESTs custom-made at the DPI microarray facility using a Biorobotics MicroGrid II Microarrayer. The experimental design of each experiment is outlined in Fig. 1. Raw spot signal intensity data was quantified from the microarray images using ImaGene 5.5 Standard Edition (Biodiscovery Inc., Marina del Rey, CA, USA). For each of the three experiments, separate channel analysis of two-coloured data was performed in the R statistical environment (R Core Team 2013) using Linear Models for Microarray Data (LIMMA; Smyth and Speed 2003; Smyth 2004, 2005) as outlined in the LIMMA user’s guide (Smyth et al. 2005). For a detailed description of the pre-processing of quantified microarray data and the subsequent analysis procedures, including assessment of the most appropriate methods to use, refer to Data File S9, available as Supplementary Material to this paper. Briefly, microarray data was background corrected using the normexp method. An offset value of 50 was used during background

Tammar wallaby mammary gland gene expression

correction to add a constant to the intensities before logtransforming so that log-ratios are shrunk towards zero at the lower intensities. Print-tip loess normalisation was performed for each microarray slide before Aquantile normalisation was applied to allow for comparisons between slides. Following normalisation, separate channel analysis was performed by fitting the linear model: E(Yg) ¼ Xag, where Yg is the vector of normalised log ratios from the arrays, E(Yg) is the expected value of Yg, X is a design matrix indicating which RNA samples have been applied to each microarray and ag is the vector of log ratios to estimate (Smyth 2004; Smyth et al. 2005). A contrast matrix was then constructed. The contrasts were linear combinations of parameters from the linear fit model: bg ¼ CTag, where bg is a vector of contrasts for gene g, C is the contrast matrix and ag is a vector of coefficients (estimated log-fold changes) obtained from the linear model fit (Smyth 2004; Smyth et al. 2005). For the assessment of differential expression an empirical Bayes method was used to moderate the standard errors for the estimated log-fold changes. This results in a more stable inference and improved power especially for experiments with a small number of arrays (such as experiment E1; Smyth 2004). ESTs having a P value (adjusted for multiple testing using the Holm’s step-down Bonferroni method) of ,0.05 were considered for further analysis. Confirmation of microarray results Two types of confirmation of results were carried out – technical replication of individual micro-array results and biological replication using additional animals. Independent experimental validation of the microarray expression results (i.e. technical replication) was performed via northern analysis for a subset of differentially expressed ESTs and for one EST representing a housekeeping gene. In addition, where a gene was represented by more than one EST on the microarray, we checked that the different ESTs gave the same pattern of differential expression. There was only one animal at each time point. The statistical analysis judges the significance of the change in expression between two time points by comparing a target EST to the variation shown by ESTs in general. However, the resulting significance tests are still only based on one animal per time point. We confirmed results in additional animals (i.e. biological replication) in two ways. Firstly, some of the northern analyses were done on different animals from those used for the microarrays. Secondly, we showed that microarray results from different time points within the same phase of lactation showed similar patterns of expression as set out below. Confirmation of gene expression differences between stages of lactation To assess the similarity of differential expression between different time points within the same phase of lactation (Table 1) two analyses were performed: an ordered list comparison and hierarchical clustering. The ‘‘OrderedList’’ Bioconductor compliant package written in the language R was utilised for the ordered list comparison (Yang et al. 2008). EST lists from each comparison outlined in Table 1 were ordered highest to lowest based on their normalised M-values (log2(fold change)).


C

Table 1. Comparison of time points within each microarray experiment

Time-point comparisons

E1: Phase 1–2A

E2: Phase 2A–2B

E3: Phase 2B–3

P18 vs L2 P18 vs L22 P26 vs L2 P26 vs L22

L62 vs L151 L62 vs L193 L87 vs L151 L87 vs L193 L110 vs L151 L110 vs L193

L130 vs L216 L130 vs L243 L130 vs L266 L151 vs L216 L151 vs L243 L151 vs L266

All possible pairs of lists from within each experiment were interrogated using OrderedList (algorithm described in Yang et al. 2006). Similarity scores were calculated by the number of shared elements Sn in the first n elements of the lists for each n. The final scores were a weighted sum over Sn where the ends of the lists receive larger weights, thus ensuring that the more strongly induced genes dominate the score. Significance of detected list similarities were estimated by OrderedList via randomly perturbing the input data (via shuffling) to compute null distributions of the similarity score from which empirical P values were deduced. In addition, OrderedLists detected how far into the lists striking similarities occur (in terms of number of genes at the top and bottom of the lists). Unsupervised hierarchical clustering was also performed using the normalised fold-change values for ESTs with significant differential expression from P26 versus L2, L110 versus L151 and L151 versus L216 comparisons. Hierarchical clustering was implemented using the R-function ‘‘hclust’’ included in the package ‘‘stats’’ (R Core Team 2014). This clustering method is based on a set of dissimilarities between the samples (using Euclidean distance). Similarity between clusters was computed using Ward’s method (Ward 1963). The branch lengths of a cluster dendrogram are measures of the distance between the samples and thus a visual representation of the variability between samples within experimental contrasts and between individual experiments. Annotation of differentially expressed ESTs Annotating ESTs typically consists of a series of repetitive tasks. Consequently, due to the large number of significantly differentially regulated ESTs, an automated approach to annotation was taken. Custom Perl, Python and shell scripts were written to format the data and obtain tammar, human, opossum, bovine and mouse Ensembl transcript identifiers (ENSTs) for each EST (via a BLASTN sequence similarity search). Because of the reasonably large evolutionary time gap between the tammar and human, cow and mouse, a relatively lenient e-value of 0.01 was used as a match cut-off threshold. The corresponding ENST identifiers were then used as input for BioMart MartView (http://www.biomart.org, accessed December 2012) and ENSG (ensemble gene) and Entrez identifiers were obtained. Gene ontology analysis The Database for Annotation, Visualisation and Integrated Discovery (DAVID; Dennis et al. 2003; Huang et al. 2009), was

D


used to identify gene ontology (GO) terms that were overrepresented by significantly differentially regulated ESTs, when compared with the complete set of ESTs found on the microarray. Gene identifiers corresponding to the significant ESTs were imported into DAVID’s functional annotation tool and the biological process GO ontology terms were obtained using DAVID’s GO-Fat algorithm. The GO-Fat algorithm defines the term specificity based on the number of associated child terms. Broadest terms (terms with the most number of child terms) are filtered out so that they do not overshadow the more specific terms. Statistically enriched GO terms were identified based on the Expression Analysis Systematic Explorer (EASE), a variant of one-tailed Fisher exact probability (Hosack et al. 2003). For all DAVID analyses a Bonferroni-corrected P value less than or equal to 0.05 was used as a standard cut-off threshold. Pathway and network analysis Further functional interpretation of each gene set in context of relevance to canonical pathways, biological function and molecular networks was generated using the bioinformatics tool Ingenuity Pathway Analysis (IPA, Version 8.7; IngenuityÒ Systems, www.ingenuity.com). Gene sets were compared with the IPA library of canonical pathways, which includes 80 metabolic and 270 signalling pathways that have been incorporated from various resources and hand-curated. IPA categorised gene sets into the following biological functions: diseases and disorders, molecular and cellular functions and physiological system development and function and ranked them according to P values. P values less than 0.05 indicate a statistically significant, non-random association between a set of significant differentially regulated genes and the set of all genes related to a given function in Ingenuity’s knowledge base. The IPA analysis determined the subcategories within each category supplied with an appropriate P value and the number of genes identified. Gene sets were used to find possible connections between genes or gene products and other genes based on interactions previously reported in the literature. Significantly differentially expressed genes were mapped onto a global molecular network developed from information contained in the Ingenuity Pathways Knowledge Base. Networks of these genes were then algorithmically generated, based on their connectivity. Since the size of the created networks could be potentially enormous, the IPA software limited the maximum number of molecules in the network to 35, leaving only the most important ones based on the number of connections for each focus gene (focus genes ¼ a subset of uploaded significant genes having direct interactions with other genes in the database) to other significant genes. A network is scored based on the number of focus molecules in the network, its size, the total number of focus molecules analysed and the total number of molecules in the knowledge database that could potentially be included in the network. The biological functions that were most significant to these networks were determined and Fischer’s exact test was used to calculate P values determining the probability that each biological function assigned to a network was due to chance alone.


Results Validation of microarray expression Independent experimental validation of the microarray results was performed by northern analysis. The relative abundance of a-casein, early lactation protein (ELP), ferritin and Proteasome 26S Subunit 6 (PSMC4) were investigated across the perinatal period and TATA-binding protein-associated factor 172 (BTAF1) was investigated over the entire time-course. With the exception of PSMC4 (a housekeeping gene), probes were selected based on their expression profiles, with the expression of a-casein, ELP and ferritin during this time period already well established in the literature. BTAF1 was selected based on the sharp decrease in expression during the transition from Phase 2B to 3. It is during this time that there is a substantial increase in the amount of milk being secreted. Sizing from the RNA gel electrophoresis photo to scanned phospho-images corroborated EST size in all cases (data not shown). In each of the columns, northern hybridisation detected a single band. With the exception of the lactation Day 171 (L171) sample from E2 and E3, the correlation between the normalised northern results and the microarray profiles was exceptionally high (see Data File S10). In silico validation involved checking that ESTs representing the same gene were consistently expressed and comparing their expression profiles to those that are known from the literature. Fig. 2 depicts the expression profiles of ESTs representing major tammar milk-protein genes. It can be seen that the pattern of expression for ESTs representing the same milk-protein gene were consistent and in accord with the literature (Nicholas et al. 1997). Both the experimental and in silico approaches to validation cast doubt on the validity of the L171 sample. It was suspected that, at the time of mammary gland sampling, either the wrong gland was taken or the L171 sample containers were incorrectly labelled. Consequently, all data from the L171 time point was eliminated from this study. Cross-validation of significant differential expression between time points This study focuses on gene expression changes between the time points that are closest to each phase boundary. The expression data analysed was primarily derived from the following linear model contrasts: P26 versus L2 (Phase 1 vs 2A), L110 versus L151 (Phase 2A vs 2B) and L151 versus L216 (Phase 2B vs 3). Because each time point originated from one animal only, we need to ensure that differences in expression between time points were not due to animal variation alone. Similarity of ordered gene lists (SOGL) has been used successfully to measure and score the level of similarity between the results of two data sets (Yang et al. 2006; Yang and Sun 2007; Schneckener et al. 2011). We used the Bioconductor package ‘‘OrderedList’’ to measure the level of similarity between the normalised timepoint comparisons in each of the three lactation phase experiments. For each comparison, OrderedList-derived rankings using regularised t-scores were generated. We observed a significant similarity (P , 0.0001) between each pair of time points. While every comparison was significant, examining the weighted overlap scores indicated that the contrasts in E2



α-casein

5

1.5

3 2

M

M

β-casein

2.0

4

1

1.0 0.5

0 ⫺1 ⫺10

⫺5

0 0

5

10

15

20

⫺10

⫺5

0

5

Day

20

3

2

M

M

15

β-lactoglobulin

4

3

1

2 1

0 ⫺10

⫺5

0

5

10

15

0 ⫺10

20

⫺5

0

5

Day WAP 8

M

6 4 2 0 62

82

102

122

10

15

20

Day

10

M

10

Day

α-lactalbumin

4

142

162

182

3 2 1 0 ⫺1 ⫺2 ⫺3 ⫺4 130

WAP

150

170

Day

190

210

230

250

Day

ELP 2

8

0

6

⫺2

4

M

M

E

⫺4 ⫺6

LLP-A

2 0

⫺8 62

82

102

122

142

162

182

Day

⫺2 130

150

170

190

210

230

250

Day

LLP-B

8

M

6 4 2 0 130

150

170

190

210

230

250

Day Fig. 2. Normalised microarray expression profiles for milk-protein ESTs. Normalised microarray expression profiles for all ESTs representing the major tammar milk proteins: a-casein (seven ESTs), b-casein (four ESTs), a-lactalbumin (three ESTs), b-lactoglobulin (eight ESTs), whey acidic protein (WAP, three ESTs), early lactation protein (ELP, four ESTs), late lactation protein A (LLP-A, two ESTs) and LLP-B (one EST).

(Phase 2A–2B) were less similar than the contrasts compared in E1 and E3. E1 contrasts, examining expression changes between pregnancy and lactation initiation, were the most similar (see Data File S1).

Hierarchical clustering of significantly differentially expressed ESTs from the P26–L2, L110–L151 and L151–L216 data sets across all experimental comparisons outlined in Table 1 was also used to examine the similarity, or ‘related-ness’, between

F



Height 0

50

100

150

200

L266_vs_L151 L216_vs_L151 L243_vs_L151 L266_vs_L130

human, mouse and cow genomes via BLASTN sequence similarity searches (publicly available through Ensembl). Surprisingly, the genome returning the least number of ‘hits’ was the tammar wallaby. We believe that this is a reflection of the lack of annotation in the current version of the tammar genome (which was only completed to 2-fold coverage) and is not a true indication of the number of genes represented by the significantly differentially regulated ESTs. The genome returning the most ‘hits’ was the human genome. Therefore, all analysis from this point used annotations obtained from the human genome.

L216_vs_L130 L243_vs_L130 P18_vs_L22 P26_vs_L22 P18_vs_L2 P26_vs_L2 L151_vs_L110 L193_vs_L110 L151_vs_L62 L151_vs_L87 L193_vs_L62 L193_vs_L87

Fig. 3. Hierarchical clustering of tammar wallaby lactation microarray experiments. Hierarchical clustering dendrogram (using Euclidean distance and Ward’s method) of the three microarray experiments and the contrasts within each experiment spanning the lactation phase boundaries. Three major cluster branches are evident, each corresponding to an individual experiment.

time-point comparisons (Fig. 3). In every instance, expression sets examining the same lactation phase comparison clustered together. The variability between lactation phase comparisons (measured by the height of the dendrogram branches in Fig. 3) demonstrates that gene expression at lactation initiation is more similar to that occurring between Phase 2A and 2B than it is to gene expression changes between Phase 2B and 3. These results provide a level of confidence that the differential expression between time points is more likely indicative of changes in the mammary gland due to stage of lactation rather than animal differences. EST annotation ESTs displaying differential regulation (P , 0.05) between the phase boundaries were aligned to the tammar wallaby, opossum,

Phase 1–2A Examining expression changes between time-points P26 (2 days before parturition) and L2 (2 days after parturition) using a cutoff P value of 0.05, a total of 2191 ESTs displayed significant differential regulation (SDR). Of these ESTs, 932 (42.5%) were upregulated and 1259 (57.5%) were downregulated. These ESTs mapped to 1161 human genes (392 upregulated and 769 downregulated). Table 2 lists the top 10 up- and downregulated annotated ESTs (full lists available in Data Files S2–4). Singleminded homolog 2 (SIM2), a potential transcription repressor, headed the upregulated ESTs, increasing over 11-fold from late pregnancy to early lactation. Collagen, type VI, alpha-3 (COL6A3) headed the downregulated list decreasing over 5-fold. These genes and overall patterns of regulation are examined further in the discussion. Gene ontology (GO) analysis was used to screen for transcriptionally regulated biological processes. The top 10 most significant GO terms associated with both up- and downregulated genes are listed in Table 3. GO analysis revealed that the main upregulated genes were those associated with transport. There was an extensive list of GO terms associated with downregulated genes (refer to Data File S5 for complete list); the major theme, however, related to translation. To identify which of the known metabolic and signalling pathways are transcriptionally regulated, pathway analysis was applied to the significant Phase 1–2A gene set. Of the 350 canonical pathways available for evaluation in the Ingenuity Pathways Knowledge Base, 37 were significant (P , 0.05; see Data File S6). IPA also provided a summary of the top biological functions represented by the gene set (Data File S7). Cellular growth and proliferation had the top score for the molecular and cellular functions category, with 276 molecules having associations. Tissue development topped the physiological system development and function category with 79 molecules being associated and cancer topped the diseases and disorders category with 333 genes of the gene set having an association. In addition to identifying significant pathways and biological functions, IPA was used to identify gene networks. A total of 25 regulatory networks were linked to the Phase 1–2A data set (Data File S8). These networks had an IPA score of $21, indicating a less than 1021 chance that the genes in the network are associated together purely due to random events. The network score is a numerical value used to rank networks according to how relevant they are to the genes in the input dataset and allows the networks to be prioritised for further study. Fig. 4 depicts the highest scoring network with an IPA

Ensembl gene identifier

Phase 1–2A ENSG00000159263 ENSG00000160211 ENSG00000167531 ENSG00000104853 ENSG00000177463 ENSG00000101230 ENSG00000153406 ENSG00000091140 ENSG00000182253 ENSG00000170890 Phase 2A–2B ENSG00000090382 ENSG00000196544 ENSG00000118705 ENSG00000001629 ENSG00000132881 ENSG00000130158 ENSG00000020129 ENSG00000169180 ENSG00000145824 ENSG00000072501 Phase 2B–3 ENSG00000033170 ENSG00000110851 ENSG00000154305 ENSG00000226567 ENSG00000105723 ENSG00000107719 ENSG00000233360 ENSG00000136270 ENSG00000211664 ENSG00000099194

Direction of regulation

Upregulated SIM2 G6PD LALBA CLPTM1 NR2C2 C20orf82 NMRAL1 DLD SYNM PLA2G1B LYZ C17orf59 RPN2 ANKIB1 C1orf89 DOCK6 NCDN XPO6 CXCL14 SMC1A FUT8 PRDM4 MIA3 AC018495.2 GSK3A KIAA1274 Z83844.1 TBRG4 IGLV2-18 SCD

Lysozyme Uncharacterised protein C17orf59 Ribophorin II Ankyrin repeat and IBR-domain containing 1 Miro domain-containing protein C1orf89 Dedicator of cytokinesis 6 Neurochondrin Exportin 6 Chemokine (C–X–C motif) ligand 14 Structural maintenance of chromosomes 1A Fucosyltransferase 8 (alpha (1,6) fucosyltransferase) PR-domain containing 4 Melanoma inhibitory activity family, member 3 No description available Glycogen synthase kinase 3 alpha KIAA1274 No description available Transforming growth factor beta regulator 4 Immunoglobulin lambda variable 2-18 Stearoyl-CoA desaturase (delta-9-desaturase)

Gene name

Single-minded homolog 2 (Drosophila) Glucose-6-phosphate dehydrogenase Lactalbumin, alphaCleft lip and palate-associated transmembrane protein 1 Nuclear receptor subfamily 2, group C, member 2 Isthmin 1 homolog (zebrafish) NmrA-like family domain containing 1 Dihydrolipoamide dehydrogenase Synemin, intermediate filament protein Phospholipase A2, group IB (pancreas)

Gene description

Table 2. Top 10 significantly differentially regulated ESTs between phase comparisons (ranked based on fold-change)

170.51 42.88 42.81 41.18 31.11 10.49 9.50 7.52 7.39 7.32

4.79 4.64 4.32 4.08 3.63 3.01 2.37 2.21 2.09 2.00

11.59 7.58 7.14 6.91 6.59 6.28 6.10 5.81 5.63 5.40

Fold-change

(Continued )

5.57E-18 2.69E-15 1.51E-14 2.36E-16 1.35E-16 1.20E-10 2.04E-10 3.72E-06 3.90E-16 6.95E-11

1.81E-02 3.03E-07 4.75E-07 6.77E-04 3.47E-06 1.51E-04 1.02E-04 1.25E-02 3.47E-02 1.48E-04

3.49E-06 1.88E-05 4.50E-05 6.83E-05 1.49E-05 2.21E-04 1.47E-06 5.89E-04 1.82E-05 9.21E-05

P value

Tammar wallaby mammary gland gene expression Reproduction, Fertility and Development G

SPINT4 SPINT4 SLC30A2 SPINT4 AAGAB IPO13 RP11-169K16.7 SCTG2 YY1 SPINT4 COL6A3 FTH1 BTAF1 FTH1 FTH1 VASP AC098850.1 CLPB SDF2 L1

Kunitz-type protease inhibitor 4 precursor Kunitz-type protease inhibitor 4 precursor Solute carrier family 30 (zinc transporter), member 2 Kunitz-type protease inhibitor 4 precursor Alpha- and gamma-adaptin binding protein Importin 13 No description available Actin, gamma 2, smooth muscle, enteric YY1 transcription factor Kunitz-type protease inhibitor 4 precursor Collagen, type VI, alpha 3 Ferritin, heavy polypeptide 1 BTAF1 RNA polymerase II Ferritin, heavy polypeptide 1 Ferritin, heavy polypeptide 1 Vasodilator-stimulated phosphoprotein No description available ClpB caseinolytic peptidase B homologue (E. coli) Stromal cell-derived factor 2-like 1 Eukaryotic translation initiation factor 3, subunit B

ENSG00000106263

7.56

25.78 19.67 16.82 14.01 13.39 11.05 9.04 8.45 8.03

9.42 8.74 8.49 7.46 6.88 5.21 4.88 4.57 4.43 3.92

5.29 5.11 4.94 4.85 4.37 4.27 3.92 3.91 3.89 3.77

Fold-change

5.83E-11

1.25E-12 2.36E-16 1.56E-14 2.46E-16 5.12E-10 2.11E-12 2.24E-11 3.00E-15 3.04E-13

2.78E-04 8.55E-04 9.31E-05 5.42E-04 1.48E-04 1.13E-03 1.08E-02 5.42E-04 3.39E-04 4.97E-02

1.08E-04 5.59E-07 3.42E-05 9.20E-06 3.70E-06 7.04E-06 8.30E-05 2.43E-04 3.11E-05 2.53E-04

P value


EIF3B

COL6A3 MME NF1 RPL4 PSAP EEF2 RPL13AP5 EEF1A1 RPS3A RPS18

Gene name

Collagen, type VI, alpha 3 Membrane metallo-endopeptidase Neurofibromin 1 Ribosomal protein L4 Prosaposin Eukaryotic translation elongation factor 2 Ribosomal protein L13a pseudogene 5 Eukaryotic translation elongation factor 1 alpha 1 pseudogene 9 Small nucleolar RNA, C/D box 73A Ribosomal protein S18

Downregulated

Gene description

Phase 1–2A ENSG00000163359 ENSG00000196549 ENSG00000196712 ENSG00000174444 ENSG00000197746 ENSG00000167658 ENSG00000236552 ENSG00000156508 ENSG00000145425 ENSG00000226225 Phase 2A–2B ENSG00000149651 ENSG00000149651 ENSG00000158014 ENSG00000149651 ENSG00000103591 ENSG00000117408 ENSG00000233954 ENSG00000163017 ENSG00000100811 ENSG00000149651 Phase 2B–3 ENSG00000163359 ENSG00000167996 ENSG00000095564 ENSG00000167996 ENSG00000167996 ENSG00000125753 ENSG00000231458 ENSG00000162129 ENSG00000128228

Ensembl gene identifier

Direction of regulation

Table 2. (Continued)

H C. J. Vander Jagt et al.



I

Table 3. Significantly differentially regulated gene ontology (GO) biological process terms during lactation phase transitions Phase comparison

Phase 1–2A

Phase 2A–2B Phase 2B–3

Upregulated GO ID

GO term

GO:0046942 GO:0015849 GO:0008015 GO:0003013 GO:0015758 GO:0008645 GO:0015749 GO:0015837 GO:0006865 GO:0008643 GO:0008104 GO:0006955 GO:0042592 GO:0051384 GO:0031960 GO:0010033 GO:0009725 GO:0032868 GO:0009719 GO:0048878 GO:0030970

Carboxylic-acid transport Organic-acid transport Blood circulation Circulatory system process Glucose transport Hexose transport Monosaccharide transport Amine transport Amino-acid transport Carbohydrate transport Protein localisation Immune response Homeostatic process Response to glucocorticoid stimulus Response to corticosteroid stimulus Response to organic substance Response to hormone stimulus Response to insulin stimulus Response to endogenous stimulus Chemical homeostasis Retrograde protein transport, ER to cytosol

score of 44. The top network was assembled using 32 genes and encompasses functions associated with post-translational modification, amino-acid metabolism and nucleic-acid metabolism. Phase 2A–2B In the progression from Phase 2A (L110) to Phase 2B (L151) a total of 180 ESTs were found to be significantly differentially regulated (P # 0.05). Of these ESTs, 81 (45%) were upregulated and 99 (55%) were downregulated. These ESTs mapped to 90 human genes; the top 10 up- and downregulated genes are listed in Table 2. Heading the upregulated gene list and increasing close to 5-fold was lysozyme, one of the antimicrobial agents found in milk. Kunitz-type protease inhibitor 4 precursor headed the downregulated list, decreasing more than 9-fold. In fact, three of the top four downregulated ESTs matched Kunitz-type protease inhibitor 4 precursor. When subjected to an additional blast sequence similarity search against tammar wallaby genomic data, they all returned ‘hits’ to the tammar wallaby early lactation protein (ELP). ELP has been found to be a marsupialspecific whey protein that appears to be a member of the Kunitz proteinase inhibitor family. The Phase 2A–2B gene set was analysed for enriched biological function GO terms (Table 3). In comparison to the Phase 1–2A gene set, this gene set had only one significant associated GO term; protein localisation (associated with the upregulated genes). As with the Phase 1–2A data set, the Ingenuity Pathways Knowledge Base was used to identify significant transcriptionally regulated metabolic and signalling pathways for the Phase 2A–2B data set. A total of 15 pathways were significantly (P , 0.05) represented by the gene set.

Downregulated P value 1.18E-04 1.18E-04 2.69E-04 2.69E-04 1.14E-03 1.74E-03 1.74E-03 3.59E-03 4.37E-03 6.52E-03 3.46E-02 5.19E-03 1.16E-02 1.40E-02 1.67E-02 1.96E-02 2.11E-02 2.43E-02 2.92E-02 2.92E-02 3.45E-02

GO ID

GO term

P value

GO:0006414 GO:0006412 GO:0042254 GO:0022613 GO:0016072 GO:0006364 GO:0042274 GO:0042273 GO:0034470 GO:0006096

Translational elongation Translation Ribosome biogenesis Ribonucleoprotein complex biogenesis rRNA metabolic process rRNA processing Ribosomal small-subunit biogenesis Ribosomal large-subunit biogenesis ncRNA processing Glycolysis

6.62E-28 3.34E-14 5.92E-08 9.67E-07 1.22E-05 1.22E-05 2.61E-04 4.33E-04 8.90E-04 5.22E-03

GO:0030199 GO:0030029 GO:0009628 GO:0009612 GO:0030036 GO:0010033 GO:0006928 GO:0035295 GO:0032964 GO:0010638

Collagen fibril organisation Actin filament-based process Response to abiotic stimulus Response to mechanical stimulus Actin cytoskeleton organisation Response to organic substance Cell motion Tube development Collagen biosynthetic process Positive regulation of organelle organisation

7.57E-03 7.99E-03 1.21E-02 2.68E-02 2.91E-02 3.15E-02 3.70E-02 4.82E-02 4.83E-02 4.86E-02

The most significant molecular and cellular function identified in IPA was amino-acid metabolism. Connective-tissue development and function headed the physiological system development and function category and as with the Phase 1–2A gene set, cancer was the most represented disease. IPA network analysis identified four regulatory networks being linked to the Phase 2A–2B data set. These networks had an IPA score of $20, indicating a less than 1020 chance that the genes in each network are linked by chance. Fig. 5 depicts the highest scoring network, which was associated with cellular assembly and organisation, developmental disorder and genetic disorder. Phase 2B–3 ESTs displaying significant differential regulation (P # 0.05) between Phases 2B (L151) and 3 (L216) numbered 237. Of these ESTs, 115 (48.5%) were upregulated and 122 (51.5%) were downregulated. Table 2 lists the top 10 up- and downregulated ESTs that mapped to the Ensembl gene collection. Fucosyltransferase 8 (FUT8) topped the upregulated gene list, increasing a huge 170-fold. FUT8 belongs to the family of fucosyltransferases and catalyses the transfer of fucose from GDP-fucose to N-linked type complex glycopeptides. The gene heading the downregulated list (.25-fold downregulation) was the same gene that headed the downregulated list for the Phase 1–2A gene set; collagen, type VI, alpha-3 (COL6A3). This gene has been found to have many functions, one of which is organising matrix components. Gene ontology analysis revealed 20 significantly associated biological functions for the Phase 2B–3 data set, 10 upregulated and 10 downregulated (Table 3). Upregulated biological

J



Network shapes Cytokine Growth factor Chemical/Drug/Toxicant Enzyme G-protein coupled receptor Ion channel Kinase Ligand-dependent nuclear receptor Peptidase Phosphatase Transcription regulator Translation regulator Transmembrane receptor Transporter microRNA Complex/Group Other

Fig. 4. Top network for Phase 1–2A comparison. Top network of significantly differentially regulated genes between Phase 1 and Phase 2A. Function of this network is associated with post-translational modification, amino-acid metabolism and nucleic-acid metabolism. Node (gene) and edge (gene relationship) symbols are described below. Red nodes represent upregulated genes, green nodes represent downregulated genes, uncoloured nodes were not identified as significantly differentially expressed but were deemed relevant and integrated into the networks by the Ingenuity Pathway Systems software. The major interactions found between these proteins are indicated by different arrows. Solid thick line, binding; solid thick arrow, acts on; dashed thin line, indirect interaction; solid thin line, direct interaction.

functions focussed around the immune system, hormone stimulus and homeostasis, while the downregulated functions were associated with cytoskeleton and matrix organisation and cell movement. Pathway analysis of the Phase 2B–3 data set identified 23 pathways as being significantly represented (P , 0.05). The implications of these pathways and biological functions are discussed further below. In addition, IPA analysis identified a total of six significant gene networks (IPA score $13). The highest-scoring network is depicted in Fig. 6. The functions associated with the top network are cellular assembly and organisation, developmental disorder and genetic disorder. A full set of results for all of the phase comparisons can be found in Data Files S2–8. In addition, the complete set of microarray data (both raw and normalised) generated for this publication has been deposited in NCBI’s Gene Expression

Omnibus (GEO; Edgar et al. 2002) and is accessible through GEO Series accession number GSE63654 (http://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=gse63654). Discussion In this study, microarray analysis was carried out to profile geneexpression changes in the tammar wallaby mammary gland across the four phases of the lactation cycle. To do this, a custom-made cDNA microarray chip containing 13 440 mammaryderived tammar wallaby expressed sequence tags (ESTs) was utilised. Although it yielded an enormous amount of data, the experimental design has some limitations, which we have overcome as much as possible. Technical validation of the microarray measurements of expression change has been carried out by comparing the results from different ESTs representing



K

Fig. 5. Top network for Phase 2A–2B comparison. Top network of significantly differentially regulated genes between Phase 2A and Phase 2B. Function of top network is associated with tissue development, carbohydrate metabolism and molecular transport. For description of figure connectors see Fig. 4.

the same gene and by northern analysis of a sample of ESTs. Biological replication has been carried out by northern analysis and by showing that additional animals from the same stage of lactation have similar patterns of gene expression. Annotation of these ESTs using marsupial genomic resources proved difficult and we found that more ESTs were annotated when they were compared with the human Ensembl gene set than to marsupial gene sequences in Ensembl. Over-representation analysis (ORA) was used in this study to statistically evaluate the fraction of differentially regulated genes for each lactation phase comparison involved in a particular Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway or associated with a particular GO term. Despite the widespread usage of ORA for the biological interpretation of microarray experimental results, the statistical tests commonly available to assess the significance of enrichment in the GO and pathway analyses are anticonservative. In addition, we also performed a series of network analyses in this study using Ingenuity Pathway Analysis (IPA).

This type of analysis is limited by the integrity of the underlying data. It is also highly probable that protein interactions important to mammary development are incompletely represented by experiments with other cell types. Despite these limitations we believe the results show many changes in gene expression that are important to the changes in function of the mammary gland during the lactation cycle and we now discuss some of these changes. Secretory activation The interpretation of high-throughput expression data from complex tissues or an organ, such as the mammary gland, is complicated by the relative changes in different cell populations, especially during temporal analysis. This phenomenon could provide insight into the nature of the changes in the cell population of the mammary gland, or may serve to confound interpretation of the data generated in this investigation.

L



Fig. 6. Top network for Phase 2B–3 comparison. Top network of significantly differentially regulated genes between Phase 2B and Phase 3. Function of top network is associated with cellular assembly and organisation, developmental disorder and genetic disorder. For description of figure connectors see Fig. 4.

Our microarray data indicated that almost double the annotated genes were significantly downregulated compared with the number of upregulated ESTs at secretory activation (Phase 1–2A comparison). A large portion of these downregulated ESTs belonged to functional groups associated with the mammary stromal compartment (e.g. immune processes and cell adhesion). This somewhat surprising result is consistent with findings in both the cow (Finucane et al. 2008) and mouse (Rudolph et al. 2003; Lemay et al. 2007). In a study of gene expression in the mouse mammary gland over an equivalent period of lactation, Rudolph et al. (2003) postulated that this did not represent a real decrease in the expression of these genes involved in the stromal compartment, but simply reflected their dilution due to the massive expansion of the epithelial compartment, as occurs with rapid aveolarisation. Unexpectedly, we found that genes associated with translation and the ribosomal machinery were also displaying a pattern of downregulation as the mammary gland shifts from a pregnant state into secretory activation. This implies that both mRNA

and protein processing are arrested, or at the very least slowed, during secretory activation. Lemay et al. (2007) observed a similar finding in the further analysis of the Rudolph et al. (2003) mouse lactation data. Lemay et al. (2007) proposed that the functions of packaging and exporting proteins in late pregnancy by the traditional secretory pathway are attenuated during lactation. They also go on to suggest that the global downregulation of genes during secretory activation does not actually reflect a change in relative cell populations, as first suggested by Rudolph et al. (2003), but rather that milk production occurs as a result of widespread transcriptional suppression of functions such as protein degradation and cellular communication. Our data further supports these claims and suggests that this is not just a eutherian lactation mechanism, but may actually be a universal mechanism for mammalian lactation. It is interesting to note that, as with the Phase 1–2A comparison, we observed a decrease in expression of ribosomalrelated genes and in biological functions associated with cell


growth, proliferation, communication and movement during the transition from Phase 2B–3. This is in contrast to what might be expected, as it is at the onset of Phase 3 when there is a significant increase in the rate of milk synthesis. Milk composition changes from having elevated levels of carbohydrate and low levels of fat and protein, to being high in fat and protein and low in carbohydrates. This further supports the theory that both mRNA and protein processing are essentially slowed during secretion. Sim2 SIM2 was the most upregulated EST in our data between Phase 1 and 2A (Table 2). In the mouse, SIM2 expression peaks during lactation. SIM2 is a transcription factor that decreases expression of signal transducer and activator of transcription 3 (STAT3) and Nuclear Factor Of Kappa Light Polypeptide Gene Enhancer In B-Cells (NFKB) and hence opposes apoptosis, which occurs at involution of the mammary gland (Scribner et al. 2011). SIM2 is also necessary for mammary gland development, where it inhibits proliferation but increases ductal development (Laffin et al. 2008). Consistent with this, it suppresses tumours in the mammary gland and prostate (Lu et al. 2011; Scribner et al. 2011). Therefore, its role in the tammar wallaby may be to inhibit proliferation of mammary cells, to support lactation and inhibit involution. Milk lipid synthesis and secretion In the tammar wallaby, fat is present in milk predominantly as triglyceride, forming lipid droplets secreted directly from the cell. At the initiation of lactation, milk lipid is low in concentration, but by Phase 3 it becomes the major energy source for the young (Green et al. 1983; Green 1984; Kwek et al. 2007). In our study we found that the expression of several key enzymes for fatty-acid synthesis is altered at the mRNA level during secretory activation (for example phospholipase A2, group IB (PLA2G1B), lipoprotein lipase (LPL) and fatty-acid-binding protein 3 (FABP3)). Indeed, upon examining the biological functions associated with the transition from pregnancy to secretory activation, it can be seen that several functions associated with lipid synthesis are upregulated (lipid biosynthetic process, sterol metabolic process, steroid metabolic process and cholesterol metabolic process in Data File S7). In addition, lipid metabolism is one of the top functions associated with five of the 25 significant Phase 1–2A networks (Data File S8). This suggests that lipid synthesis is upregulated coincidently with secretory activation, where milk-fat droplets move from epithelial cells to the lumen just prior to lactation. This finding once again supports what has already been established in the mouse mammary gland (as reviewed in Anderson et al. 2007). In their review, Anderson et al. (2007) imply that sterol regulatory element binding protein (SREBP)-1c plays a major role in secretory activation with regard to lipid biosynthesis. The SREBP family of transcription factors is recognised as regulating fatty-acid and cholesterol biosynthesis. In our study we found that the SREBP-1c gene was one of the significantly upregulated genes during secretory activation (refer to Data File S2).


M

At the transition of Phase 2B to 3, milk production increases significantly and the gross milk composition changes markedly. Milk lipid levels increase dramatically as lactation progresses. This is not reflected in the significant biological functions or network analyses for the Phase 2B–3 dataset. It is interesting to note, however, that lipoprotein lipase (LPL) and fatty-acidbinding protein 3 (FABP3, a mammary-derived growth inhibitor) both showed upregulation during secretory activation and then again in the Phase 2B to 3 transition. LPL’s primary function is the hydrolysis of triglycerides of circulating chylomicrons and very low-density lipoproteins (VLDL) into fattyacid and glycerol molecules (as reviewed in Mead et al. 2002). LPL has also previously been found to be upregulated at secretory activation in the cow (Finucane et al. 2008). FABP3 is thought to participate in the uptake, intracellular metabolism and transport of long-chain fatty acids. It may also be responsible for the modulation of cell growth and proliferation (Bo¨rchers et al. 1997). Lactation and immunity The phases of tammar lactation are broadly correlated with major periods of growth and maturation of the young, including the development of the immune system. The tammar young is born in a very immature state. Immediately after birth it climbs from the birth canal and into the maternal pouch where it attaches (for the length of Phase 2A) permanently to one of the four maternal teats (Tyndale-Biscoe and Janssens 1988). This transition from Phase 1 to Phase 2A is arguably the time of greatest exposure to potential microbial pathogens and the transfer of maternal antibodies and other protective factors through colostrum is vital for the immuno-protection of the young. It was therefore surprising to find in the Phase 1–2A comparison that the majority of significantly differentially regulated ESTs associated with immune genes were downregulated, a finding that is not consistent with what is currently known to be true at the protein level (Joss et al. 2009). This downregulation of the immune system was also reflected in the significant GO terms. Genes associated with lymphocyte-mediated immunity and humoral immune response were statistically downregulated during early lactation relative to late pregnancy and negative regulation of B-cell activation was statistically upregulated. Similar results have been observed in mice (Lemay et al. 2007). The reason suggested for this phenomenon is that perhaps these genes, or a subset of them, are unique to colostrum and thus their transcription is switched off by the time lactation is truly established. However, Lemay et al. (2007) were looking at lactation Day 9 in mice (of an 18–21-day lactation cycle). We observed this at Day 2 of lactation and would expect that colostrum is still being made in the mammary gland during this period. Although colostrum is generally thought to be secreted around parturition and up to 48 h postpartum (Kruse 1983; Neville et al. 1991), it is synthesised prior to parturition and therefore initiation of full lactation is accompanied by a decrease in expression of genes needed only for colostrum. During Phase 2A, the young’s immune system undergoes rapid development but does not achieve a mature form until nearing the end of the phase and passing into Phase 2B. At the

N


start of Phase 2B the young relinquishes the teat, suckling less frequently but still permanently inhabiting the pouch. We see this increase in immuno-competence reflected by a further decrease in the expression of immune-related genes as lactation moves from Phase 2A to Phase 2B (refer to Data File S3). However, we see a significant increase in immune-related genes, functions and pathways during the transition from Phase 2B to 3. This observation is expected, as during Phase 3 the young has mature body systems and is moving towards independence from the mother. The young vacates the pouch for the first time and starts grazing on herbage, exposing it to the outside environment and coincidently to an array of pathogens. This increase in immune processes has been documented previously (Deane and Cooper 1988; Deane et al. 1990). Regulation of milk-protein genes All ESTs representing milk-protein genes increased in expression between late pregnancy and early lactation (a-casein increased 6–10.55-fold, b-casein increased 1.5–2.5-fold, a-lactalbumin increased 2–7.1-fold, b-lactoglobulin increased 1.2–5-fold and ELP increased 3.9–4.4-fold; for expression profiles of individual ESTs see Fig. 2; for normalised foldchange values of ESTs see Data Files S2–4). Some milk-protein genes stayed at high levels of expression throughout lactation (e.g. caseins) but ELP declined at about Day 100, whey acidic protein (WAP) increased at Day 100 and declined at Day 170 while late lactation protein A (LLP A) and B increased at Day 150. These changes are in agreement with previous reports (Nicholas et al. 1997). Although these changes in gene expression explain part of the change in milk-protein synthesis, Lemay et al. (2007) suggest that changes in translation due to changes in mRNA stability also contribute. This concept is not new and is supported experimentally in several published studies involving the casein genes (Guyette et al. 1979; Choi et al. 2004; Moshel et al. 2006). Nadin-Davis and Mezl (1985) observed that casein and WAP (the major milk whey protein) mRNA have longer poly-A tails, protecting the mRNA from degradation by endoribonucleases, when milk-protein production is high. In support of this hypothesis, we found that several genes associated with poly(A) tail shortening or degradation were significantly downregulated at secretory activation in our data. These genes include: decapping enzyme, scavenger (DCPS) (down 2.1-fold), terminal uridylyl transferase 1, U6 snRNA-specific (TUT1) (down 3.7-fold), Poly(A) Binding Protein, Cytoplasmic 1 (PABPC1) (down 2.8-fold) and DEAD (Asp-Glu-Ala-Asp) box helicase 5 (DDX5) (down 2.4-fold). Of these genes, DCPS actually displayed significant differential regulation in all three data sets, decreasing further as lactation moved from Phase 2A to 2B and then increasing again between Phases 2B and 3. This pattern of regulation is the opposite to the major milk protein WAP, perhaps lending weight to the theory that mRNA stability plays a major role in the regulation of milk-protein gene translation. The rate of casein translation has also been found to be reduced by amino-acid deprivation (Moshel et al. 2006). Our experimental findings lend weight to the suggestion that aminoacid availability is implicated in casein translation, with genes


belonging to the biological function group amino-acid and amine transport being significantly upregulated with the onset of lactation. Additionally, the canonical pathways of phenylalanine metabolism, lysine biosynthesis and methionine metabolism were all significantly represented by the Phase 1–2A data set, with the majority of genes within these pathways being upregulated (refer to Data Files S2 and S6). This potentially indicates that there is a cellular increase in the availability of these amino acids during the lactation initiation phase. WAP Whey acidic protein (WAP) is the principal whey protein found in a number of eutherian species, including mouse (Hennighausen and Sippel 1982), rat (Campbell et al. 1984), pig (Simpson et al. 1998b) and rabbit (Devinoy et al. 1988). It has also been found to be a major component of tammar milk in Phase 2B of lactation (Nicholas et al. 1995). In mice, rats and pigs, maximal induction of the WAP gene requires the presence of insulin, glucocorticoid and prolactin (Pittius et al. 1988; Schoenenberger et al. 1990; Puissant and Houdebine 1991; Simpson et al. 1998b). In contrast, maximal induction of WAP in mammary explants from Day 24 pregnant tammars requires insulin, glucocorticoid, prolactin as well as triiodothyronine (T3), which differs from the endocrine requirements of eutherians (Simpson et al. 1998a; Simpson and Nicholas 2002). In our study, we found that ESTs representing WAP increased expression significantly between Phase 2A and 2B, before decreasing again between Phases 2B and 3 (Fig. 2). This finding is supported by the literature (Simpson et al. 1998a; Simpson and Nicholas 2002). It has also been found that the level of brushtail possum (Trichosurus vulpecula) WAP mRNA closely correlates to the level of circulating prolactin in the lactating mother (Demmer et al. 2001). We found that prolactin signalling was one of the significant pathways represented by the Phase 2B–3 data set and that two (FBJ murine osteosarcoma viral oncogene homolog (FOS) and Nuclear Receptor Subfamily 3, Group C, Member 1 (Glucocorticoid Receptor) (NR3C1)) of the three significant differentially regulated genes (refer to Data Files S4 and S6) were downregulated following the same pattern of expression as WAP (despite the increase of milk production during this time). FUT8 a1,6-Fucosyltransferase (FUT8) catalyses the transfer of a fucose residue to N-linked oligosaccharides on glycoproteins (Wilson et al. 1976). In our data, we observed an EST representing FUT8 being significantly upregulated (over 170-fold) between Phase 2B and 3 (Table 2). We did not see any significant differential expression in either of the other phase comparisons. Studies of FUT8 double-knockout mice lacking this enzyme have observed growth retardation and early death during post-natal development (Wang et al. 2005). Core fucosylation has been shown to be crucial in the binding of epidermal growth factor receptor (EGFR) to its ligand, which subsequently stimulates cell differentiation and growth (Wang et al. 2005). Members of the EGFR system are essential local regulators of mammary gland development and function (Dahlhoff et al. 2011).


In the tammar wallaby, a recent study found that, during lactation, 30% of N-linked glycosylated structures contained a core (a1-6) fucose and that these structures exhibit statistically significant temporal changes over the lactation cycle (Wongtrakul-Kish et al. 2012). Glycans are involved in a number of different functions and cellular processes, such as cell adhesion and signalling (Ohtsubo and Marth 2006; Yago et al. 2010). This involvement mediates different downstream processes, one of which is immune system development and protection. Wongtrakul-Kish et al. (2012) speculated that the changes in structures indicated their developmental and functional significance in the maternal milk protection of the infant. The extensive upregulation of FUT8 we observed at the Phase 2B–3 boundary corresponds to the time period when the joey vacates the maternal pouch for the first time. It is possible that the increase of FUT8 in the mammary gland at this time is reflective of the increased immunological protection required by the young at this time (Deane and Cooper 1988; Deane et al. 1990). Conclusions High-throughput technology, such as microarray experimentation, often produces more data than can be feasibly investigated, as was evident in this study. However, despite recent advancements in transcriptomics, microarrays still provide a powerful (and relatively cheap) methodology to highlight known and novel genes for further investigation. Acknowledgements The authors would like to acknowledge the contributions of Jenny Lee, Jo Argento, Jenny Wadeson and Shannon Simpson who assisted with the animal husbandry and wet laboratory experiments. Thank you for your hard work.

References Akers, M. R. (2002). ‘Lactation and the Mammary Gland’. (Iowa State Press: Ames, Iowa.) Anderson, S. M., Rudolph, M. C., Mcmanaman, J. L., and Neville, M. C. (2007). Key stages in mammary gland development. Secretory activation in the mammary gland: it’s not just about milk-protein synthesis! Breast Cancer Res. 9, 204. doi:10.1186/BCR1653 Bininda-Emonds, O. R. P., Cardillo, M., Jones, K. E., Macphee, R. D. E., Beck, R. M. D., Grenyer, R., Price, S. A., Vos, R. A., Gittleman, J. L., and Purvis, A. (2007). The delayed rise of present-day mammals. Nature 446, 507–512. doi:10.1038/NATURE05634 Bo¨rchers, T., Hohoff, C., Buhlmann, C., and Spener, F. (1997). Heart-type fatty-acid-binding protein – involvement in growth inhibition and differentiation. Prostaglandins Leukot. Essent. Fatty Acids 57, 77–84. doi:10.1016/S0952-3278(97)90496-8 Campbell, S. M., Rosen, J. M., Hennighausen, L. G., Strech-Jurk, U., and Sippel, A. E. (1984). Comparison of the whey acidic protein gene of the rat and mouse. Nucleic Acids Res. 12, 8685–8697. doi:10.1093/NAR/12. 22.8685 Choi, K. M., Barash, I., and Rhoads, R. E. (2004). Insulin and prolactin synergistically stimulate b-casein messenger ribonucleic acid translation by cytoplasmic polyadenylation. Mol. Endocrinol. 18, 1670–1686. doi:10.1210/ME.2003-0483 Clarkson, R. W. E., Wayland, M. T., Lee, J., Freeman, T., and Watson, C. J. (2004). Gene expression profiling of mammary-gland development


O

reveals putative roles for death receptors and immune mediators in post-lactational regression. Breast Cancer Res. 6, R92–R109. doi:10.1186/ BCR754 Dahlhoff, M., Blutke, A., Wanke, R. D., Wolf, E., and Schneider, M. R. (2011). In vivo evidence for epidermal growth factor receptor (EGFR)mediated release of prolactin from the pituitary gland. J. Biol. Chem. 286, 39 297–39 306. doi:10.1074/JBC.M111.243493 Deane, E. M., and Cooper, D. W. (1988). Immunological development. In ‘The Developing Marsupial: Models for Biomedical Research’. (Eds C. H. Tyndale-Biscoe and P. A. Janssens.) pp. 190–199. (Springer: Berlin.) Deane, E. M., Cooper, D. W., and Renfree, M. B. (1990). Immunoglobulin G levels in fetal and newborn tammar wallabies (Macropus eugenii). Reprod. Fertil. Dev. 2, 369–375. doi:10.1071/RD9900369 Demmer, J., Stasiuk, S. J., Grigor, M. R., Simpson, K. J., and Nicholas, K. R. (2001). Differential expression of the whey acidic protein gene during lactation in the brushtail possum (Trichosurus vulpecula). Biochim. Biophys. Acta 1522, 187–194. doi:10.1016/S0167-4781(01)00334-7 Dennis, G., Sherman, B., Hosack, D., Yang, J., Gao, W., Lane, H. C., and Lempicki, R. (2003). DAVID: Database for Annotation, Visualisation and Integrated Discovery. Genome Biol. 4, 3. doi:10.1186/GB-20034-5-P3 Devinoy, E., Hubert, C., Jolivet, G., Thepot, D., Clergue, N., Desaleux, M., Dion, M., Servely, J. L., and Houdebine, L. M. (1988). Recent data on the structure of rabbit milk-protein genes and on the mechanism of the hormonal control of their expression. Reprod. Nutr. Dev. 28, 1145–1164. doi:10.1051/RND:19880711 Edgar, R., Domrachev, M., and Lash, A. E. (2002). Gene Expression Omnibus: NCBI gene expression and hybridisation array data repository. Nucleic Acids Res. 30(1), 207–210. doi:10.1093/NAR/30.1.207 Finucane, K. A., Mcfadden, T. B., Bond, J. P., Kennelly, J. J., and Zhao, F. Q. (2008). Onset of lactation in the bovine mammary gland: gene expression profiling indicates a strong inhibition of gene expression in cell proliferation. Funct. Integr. Genomics 8, 251–264. doi:10.1007/S10142008-0074-Y Green, B. (1984). Composition of milk and energetics of growth. In ‘Physiological strategies in lactation: the Proceedings of a Symposium held at the Zoological Society of London’. pp. 369–387. (Academic Press: London.) Green, B., Griffiths, M., and Leckie, R. (1983). Qualitative and quantitative changes in milk fat during lactation in the tammar wallaby (Macropus eugenii). Aust. J. Biol. Sci. 36, 455–461. Guyette, W. A., Matusik, R. J., and Rosen, J. M. (1979). Prolactin-mediated transcriptional and post-transcriptional control of casein gene expression. Cell 17, 1013–1023. doi:10.1016/0092-8674(79)90340-4 Hasselbalch, H., Jeppesen, D., Engelmann, M., Michaelsen, K., and Nielsen, M. (1996). Decreased thymus size in formula-fed infants compared with breastfed infants. Acta Paediatr. 85, 1029–1032. doi:10.1111/J.16512227.1996.TB14211.X Hennighausen, L. G., and Sippel, A. E. (1982). Mouse whey acidic protein is a novel member of the family of ‘four-disulfide core’ proteins. Nucleic Acids Res. 10, 2677–2684. doi:10.1093/NAR/10.8.2677 Hosack, D. A., Dennis, G., Jr., Sherman, B. T., Lane, H. C., and Lempicki, R. A. (2003). Identifying biological themes within lists of genes with EASE. Genome Biol. 4, R70. doi:10.1186/GB-2003-4-10-R70 Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi:10.1038/NPROT.2008.211 Jenness, R. (1986). Lactational performance of various mammalian species. J. Dairy Sci. 69, 869–885. doi:10.3168/jds.S0022-0302(86) 80478-7 Joss, J. L., Molloy, M. P., Hinds, L., and Deane, E. (2009). A longitudinal study of the protein components of marsupial milk from birth to weaning

P


in the tammar wallaby (Macropus eugenii). Dev. Comp. Immunol. 33, 152–161. doi:10.1016/J.DCI.2008.08.002 Khalil, E., Digby, M. R., Thomson, P. C., Lefe`vre, C., Mailer, S. L., Pooley, C., and Nicholas, K. R. (2011). Acute involution in the tammar wallaby: identification of genes and putative novel milk proteins implicated in mammary gland function. Genomics 97, 372–378. doi:10.1016/ J.YGENO.2011.03.003 Kruse, P. E. (1983). The importance of colostral immunoglobulins and their absorption from the intestine of the newborn animals. Ann. Rech. Vet. 14, 349–353. Kwek, J. H. L., Wijesundera, C., Digby, M. R., and Nicholas, K. R. (2007). The endocrine regulation of milk-lipid synthesis and secretion in tammar wallaby (Macropus eugenii). Biochim. Biophys. Acta 1770, 48–54. doi:10.1016/J.BBAGEN.2006.06.021 Laffin, B., Wellberg, E., Kwak, H. I., Burghardt, R. C., Metz, R. P., Gustafson, T., Schedin, P., and Porter, W. W. (2008). Loss of singleminded-2s in the mouse mammary gland induces an epithelial– mesenchymal transition associated with upregulation of slug and matrix metalloprotease 2. Mol. Cell. Biol. 28, 1936–1946. doi:10.1128/MCB. 01701-07 Lefe`vre, C. M., Digby, M. R., Whitley, J. C., Strahm, Y., and Nicholas, K. R. (2007). Lactation transcriptomics in the Australian marsupial, Macropus eugenii: transcript sequencing and quantification. BMC Genomics 8, 417. doi:10.1186/1471-2164-8-417 Lemay, D. G., Neville, M. C., Rudolph, M. C., Pollard, K. S., and German, J. B. (2007). Gene regulatory networks in lactation: identification of global principles using bioinformatics. BMC Syst. Biol. 1, 56. doi:10.1186/1752-0509-1-56 Lu, B., Asara, J. M., Sanda, M. G., and Arredouani, M. S. (2011). The role of the transcription factor SIM2 in prostate cancer. PLoS ONE 6, e28837. doi:10.1371/JOURNAL.PONE.0028837 Luo, Z.-X., Ji, Q., Wible, J. R., and Yuan, C.-X. (2003). An early Cretaceous tribosphenic mammal and metatherian evolution. Science 302, 1934– 1940. doi:10.1126/SCIENCE.1090718 Master, S. R., Hartman, J. L., D’Cruz, C. M., Moody, S. E., Keiper, E. A., Ha, S. I., Cox, J. D., Belka, G. K., and Chodosh, L. A. (2002). Functional microarray analysis of mammary organogenesis reveals a developmental role in adaptive thermogenesis. Mol. Endocrinol. 16, 1185–1203. doi:10.1210/MEND.16.6.0865 Mead, J. R., Irvine, S. A., and Ramji, D. P. (2002). Lipoprotein lipase: structure, function, regulation and role in disease. J. Mol. Med. 80, 753–769. doi:10.1007/S00109-002-0384-9 Menzies, K. K., Lefe`vre, C., Sharp, J. A., Macmillan, K. L., Sheehy, P. A., and Nicholas, K. R. (2009). A novel approach identified the FOLR1 gene, a putative regulator of milk-protein synthesis. Mamm. Genome 20, 498–503. doi:10.1007/S00335-009-9207-4 Messer, M., and Green, B. (1979). Milk carbohydrates of marsupials. II. Quantitative and qualitative changes in milk carbohydrates during lactation in the tammar wallaby (Macropus eugenii). Aust. J. Biol. Sci. 32, 519–531. Moshel, Y., Rhoads, R. E., and Barash, I. (2006). Role of amino acids in translational mechanisms governing milk-protein synthesis in murine and ruminant mammary epithelial cells. J. Cell. Biochem. 98, 685–700. doi:10.1002/JCB.20825 Nadin-Davis, S. A., and Mezl, V. A. (1985). Variation in the lack of polyadenylation of the rat milk-protein mRNAs during the lactation cycle. Int. J. Biochem. 17, 1067–1075. doi:10.1016/0020-711X(85) 90038-2 Neville, M. C., Allen, J. C., Archer, P. C., Casey, C. E., Seacat, J., Keller, R. P., Lutes, V., Rasbach, J., and Neifert, M. (1991). Studies in human lactation: milk volume and nutrient composition during weaning and lactogenesis. Am. J. Clin. Nutr. 54, 81–92.


Nicholas, K. R. (1988). Asynchronous dual lactation in a marsupial, the tammar wallaby (Macropus eugenii). Biochem. Biophys. Res. Commun. 154, 529–536. doi:10.1016/0006-291X(88)90172-6 Nicholas, K. R., Wilde, C. J., Bird, P. H., Hendry, K. A. K., Tregenza, K., and Warner, B. 1995. Asynchronous concurrent secretion of milk proteins in the tammar wallaby (Macropus eugenii). In ‘Intercellular Signalling in the Mammary Gland’. (Ed. C. J. W. E. Al.) pp. 153–170. (Plenum Press: New York.) Nicholas, K., Simpson, K., Wilson, M., Trott, J., and Shaw, D. (1997). The tammar wallaby: a model to study putative autocrine-induced changes in milk composition. J. Mammary Gland Biol. Neoplasia 2, 299–310. doi:10.1023/A:1026392623090 Nilsson, M. A., Churakov, G., Sommer, M., Tran, N. V., Zemann, A., Brosius, J., and Schmitz, J. (2010). Tracking marsupial evolution using archaic genomic retroposon insertions. PLoS Biol. 8, e1000436. doi:10.1371/JOURNAL.PBIO.1000436 Oakes, S. R., Hilton, H. N., and Ormandy, C. J. (2006). Key stages in mammary gland development. The alveolar switch: coordinating the proliferative cues and cell-fate decisions that drive the formation of lobuloalveoli from ductal epithelium. Breast Cancer Res. 8, 207. doi:10.1186/BCR1411 Ohtsubo, K., and Marth, J. D. (2006). Glycosylation in cellular mechanisms of health and disease. Cell 126, 855–867. doi:10.1016/J.CELL.2006. 08.019 Pittius, C. W., Sankaran, L., Topper, Y. J., and Hennighausen, L. (1988). Comparison of the regulation of the whey acidic protein gene with that of a hybrid gene containing the whey acidic protein gene promoter in transgenic mice. Mol. Endocrinol. 2, 1027–1032. doi:10.1210/MEND2-11-1027 Puissant, C., and Houdebine, L. M. (1991). Cortisol induces rapid accumulation of whey acid protein mRNA but not of asl and b-casein mRNA in rabbit mammary explants. Cell Biol. Int. Rep. 15, 121–129. doi:10.1016/ 0309-1651(91)90103-P R Core Team (2013). ‘R: A language and environment for statistical computing’. (R Foundation for Statistical Computing, Vienna, Austria.) Available from http://www.R-project.org/ [Verified 13 January 2015] Ron, M., Israeli, G., Seroussi, E., Weller, J., Gregg, J., Shani, M., and Medrano, J. (2007). Combining mouse mammary gland gene expression and comparative mapping for the identification of candidate genes for QTL of milk production traits in cattle. BMC Genomics 8, 183. doi:10.1186/1471-2164-8-183 Rudolph, M. C., Mcmanaman, J. L., Hunter, L., Phang, T., and Neville, M. C. (2003). Functional development of the mammary gland: use of expression profiling and trajectory clustering to reveal changes in gene expression during pregnancy, lactation and involution. J. Mammary Gland Biol. Neoplasia 8, 287–307. doi:10.1023/B:JOMG.0000010030. 73983.57 Schneckener, S., Arden, N., and Schuppert, A. (2011). Quantifying stability in gene-list ranking across microarray-derived clinical biomarkers. BMC Med. Genomics 4, 73. doi:10.1186/1755-8794-4-73 Schoenenberger, C.-A., Zuk, A., Groner, B., Jones, W., and Andres, A.-C. (1990). Induction of the endogenous whey acidic protein (Wap) gene and a Wap-myc hybrid gene in primary murine mammary organoids. Dev. Biol. 139, 327–337. doi:10.1016/0012-1606(90)90302-Y Scribner, K. C., Wellberg, E. A., Metz, R. P., and Porter, W. W. (2011). Singleminded-2s (Sim2s) promotes delayed involution of the mouse mammary gland through suppression of Stat3 and NFKB. Mol. Endocrinol. 25, 635–644. doi:10.1210/ME.2010-0423 Simpson, K. J., and Nicholas, K. R. (2002). The comparative biology of whey proteins. J. Mammary Gland Biol. Neoplasia 7, 313–326. doi:10.1023/A:1022856801175



Simpson, K., Shaw, D., and Nicholas, K. (1998a). Developmentally regulated expression of a putative protease inhibitor gene in the lactating mammary gland of the tammar wallaby, Macropus eugenii. Comp. Biochem. Physiol. B, Biochem. Mol. Biol. 120, 535–541. doi:10.1016/ S0305-0491(98)10040-8 Simpson, K. J., Bird, P., Shaw, D., and Nicholas, K. (1998b). Molecular characterisation and hormone-dependent expression of the porcine whey acidic protein gene. J. Mol. Endocrinol. 20, 27–35. doi:10.1677/JME. 0.0200027 Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3, Article 3. doi:10.2202/1544-6115.1027 Smyth, G. K. (2005). Limma: linear models for microarray data. In ‘Bioinformatics and Computational Biology Solutions using R and Bioconductor’. (Eds R. Gentleman, V. Carey, S. Dudoit, R. Irizarry, W. Huber) pp. 397–420. (Springer: New York) Smyth, G. K., and Speed, T. (2003). Normalisation of cDNA microarray data. Methods 31, 265–273. doi:10.1016/S1046-2023(03)00155-5 Smyth, G. K., Ritchie, M. E., Thorne, N., and Wettenhall, J. (2005). ‘LIMMA: Linear Models for Microarray Data User’s Guide’. Available from http://www.bioconductor.org/packages/release/bioc/vignettes/ limma/inst/doc/usersguide.pdf [Verified 13 January 2015] Stein, T., Morris, J., Davies, C., Weber-Hall, S., Duffy, M.-A., Heath, V., Bell, A., Ferrier, R., Sandilands, G., and Gusterson, B. (2004). Involution of the mouse mammary gland is associated with an immune cascade and an acute-phase response, involving LBP, CD14 and STAT3. Breast Cancer Res. 6, R75–R91. doi:10.1186/BCR753 Thompson, J. M. D., Becroft, D. M. O., and Mitchell, E. A. (2000). Previous breastfeeding does not alter thymic size in infants dying of sudden infant death syndrome. Acta Paediatr. 89, 112–114. doi:10.1111/J.1651-2227. 2000.TB01198.X Tyndale-Biscoe, C., and Janssens, P. (1988). ‘The Developing Marsupial: Models for Biomedical Research’. (Springer-Verlag: Heidelberg.) Wang, X., Inoue, S., Gu, J., Miyoshi, E., Noda, K., Li, W., MizunoHorikawa, Y., Nakano, M., Asahi, M., Takahashi, M., Uozumi, N.,

Q

Ihara, S., Lee, S. H., Ikeda, Y., Yamaguchi, Y., Aze, Y., Tomiyama, Y., Fujii, J., Suzuki, K., Kondo, A., Shapiro, S. D., Lopez-Otin, C., Kuwaki, T., Okabe, M., Honke, K., and Taniguchi, N. (2005). Dysregulation of TGF-b1 receptor activation leads to abnormal lung development and emphysema-like phenotype in core fucose-deficient mice. Proc. Natl. Acad. Sci. USA 102, 15 791–15 796. doi:10.1073/PNAS. 0507375102 Ward, J. H. (1963). Hierarchical grouping to optimise an objective function. J. Am. Stat. Assoc. 58, 236–244. doi:10.1080/01621459.1963.10500845 Wilson, J. R., Williams, D., and Schachter, H. (1976). The control of glycoprotein synthesis: N-acetylglucosamine linkage to a mannose residue as a signal for the attachment of L-fucose to the asparaginelinked N-acetylglucosamine residue of glycopeptide from alpha1-acid glycoprotein. Biochem. Biophys. Res. Commun. 72, 909–916. doi:10.1016/S0006-291X(76)80218-5 Wongtrakul-Kish, K., Kolarich, D., Pascovici, D., Joss, J. L., Deane, E., and Packer, N. H. (2012). Characterisation of N- and O-linked glycosylation changes in milk of the tammar wallaby (Macropus eugenii) over lactation. Glycoconj. J. 30(5), 523–536. doi:10.1007/S10719-0129452-8 Yago, T., Fu, J., Mcdaniel, J. M., Miner, J. J., Mcever, R. P., and Xia, L. (2010). Core 1-derived O-glycans are essential E-selectin ligands on neutrophils. Proc. Natl. Acad. Sci. USA 107, 9204–9209. doi:10.1073/ PNAS.1003110107 Yang, X., and Sun, X. (2007). Meta-analysis of several gene lists for distinct types of cancer: a simple way to reveal common prognostic markers. BMC Bioinformatics 8, 118. doi:10.1186/1471-2105-8-118 Yang, X., Bentink, S., Scheid, S., and Spang, R. (2006). Similarities of ordered gene lists. J. Bioinform. Comput. Biol. 4, 693–708. doi:10.1142/ S0219720006002120 Yang, X., Scheid, S., and Lottaz, C. (2008). ‘OrderedList: Similarities of Ordered Gene Lists. R package version 1.38.0’. Available from http://compdiag.molgen.mpg.de/software/index.shtml [Verified 13 January 2015]

www.publish.csiro.au/journals/rfd

Expression of the mammary gland-specific tammar wallaby early lactation protein gene is maintained in vitro in the absence of prolactin.

Alkaline ribonuclease and ribonuclease inhibitor in mammary gland during the lactation cycle and in the R3230AC mammary tumour.

Prolactin receptor gene expression in rat mammary gland and liver during pregnancy and lactation.

Prolactin binding in rat mammary gland during pregnancy and lactation.

Relaxin gene expression in the sow corpus luteum during the cycle, pregnancy, and lactation.

Developmental changes in mitochondria during the transition into lactation in the mouse mammary gland.

Mammary gland specific knockdown of the physiological surge in Cx26 during lactation retains normal mammary gland development and function.

Lactation during cholestasis: Role of ABC proteins in bile acid traffic across the mammary gland.

Miz1 deficiency in the mammary gland causes a lactation defect by attenuated Stat5 expression and phosphorylation.

Transcriptomic Analysis of the Mouse Mammary Gland Reveals New Insights for the Role of Serotonin in Lactation.

Pyrimidine nucleotide synthesis in the rat mammary gland: changes in the lactation cycle and effects of diabetes.

Pregnancy and lactation affect the microvasculature of the mammary gland in mice.

Evolution of prolactin receptors in rabbit mammary gland during pregnancy and lactation.

Hormonal regulation of beta-adrenergic receptors in the rat mammary gland during the estrous cycle and lactation: role of sex steroids and prolactin.

Impact of diet deprivation and subsequent overallowance during gestation on mammary gland development and lactation performance.

Production of parathyroid hormone-related protein by the rat mammary gland in pregnancy and lactation.

Expression of Human NSAID Activated Gene 1 in Mice Leads to Altered Mammary Gland Differentiation and Impaired Lactation.

Changes in Harderian gland activity in the female golden hamster during the oestrous cycle, pregnancy and lactation.

Regulation of mouse mammary-gland gamma-glutamyltranspeptidase mRNA during pregnancy, lactation and weaning.

Clock circadian regulator (CLOCK) gene network expression patterns in bovine adipose, liver, and mammary gland at 3 time points during the transition from pregnancy into lactation.

Dehydrogenases, nucleic acids and soluble proteins in mammary gland of the Mongolian gerbil during pregnancy, lactation and involution.

carbohydrate ratio regulate the expression of SNAT2 amino acid transporter in the mammary gland during lactation.

Evaluation of Suitable Internal Control Genes for RT-qPCR in Yak Mammary Tissue during the Lactation Cycle.

Effect of hypothyroidism on the expression of nuclear receptors and their co-regulators in mammary gland during lactation in the rat.