Gene expression profiling of mammary gland development reveals putative roles for death receptors and immune mediators in post-lactational regression

Introduction In order to gain a better understanding of the molecular processes that underlie apoptosis and tissue regression in mammary gland, we undertook a large-scale analysis of transcriptional changes during the mouse mammary pregnancy cycle, with emphasis on the transition from lactation to involution. Method Affymetrix microarrays, representing 8618 genes, were used to compare mammary tissue from 12 time points (one virgin, three gestation, three lactation and five involution stages). Six animals were used for each time point. Common patterns of gene expression across all time points were identified and related to biological function. Results The majority of significantly induced genes in involution were also differentially regulated at earlier stages in the pregnancy cycle. This included a marked increase in inflammatory mediators during involution and at parturition, which correlated with leukaemia inhibitory factor–Stat3 (signal transducer and activator of signalling-3) signalling. Before involution, expected increases in cell proliferation, biosynthesis and metabolism-related genes were observed. During involution, the first 24 hours after weaning was characterized by a transient increase in expression of components of the death receptor pathways of apoptosis, inflammatory cytokines and acute phase response genes. After 24 hours, regulators of intrinsic apoptosis were induced in conjunction with markers of phagocyte activity, matrix proteases, suppressors of neutrophils and soluble components of specific and innate immunity. Conclusion We provide a resource of mouse mammary gene expression data for download or online analysis. Here we highlight the sequential induction of distinct apoptosis pathways in involution and the stimulation of immunomodulatory signals, which probably suppress the potentially damaging effects of a cellular inflammatory response while maintaining an appropriate antimicrobial and phagocytic environment.


Introduction
Mammary gland development during the pregnancy cycle is characterized by successive phases of cell growth, differentiation, high metabolic activity and apoptosis. At the ultrastructural level this includes dramatic changes in tissue architecture, involving ductal epithelial branching and morphogenesis, invasion of tissue compartments, vas-cularization and subsequent organized remodelling. These events are highly reproducible and strictly controlled at the transcriptional level by circulating hormones and locally derived factors [1]. Thus, many transcription factors have been shown either to directly affect this developmental program or to exhibit altered activity at specific stages in the pregnancy cycle [1].
Our laboratory has focused on the role of transcription factors in postlactational regression of the gland after weaning. Involution can be divided into at least two phases [2][3][4], broadly comprising the following: an initial reversible phase whereby the gland maintains its gross morphology but undergoes a substantial increase in the rate of epithelial cell apoptosis [5]; and a secondary irreversible phase, which involves the destruction of basement membrane by matrix metalloproteinases, phagocytic clearance of milk, and apoptotic bodies and alveolar collapse [6]. Immune cells are present at all stages of mammary development, including involution [7], but the precise role of the immune system during postlactational regression has yet to be fully established. Many genes have been shown to be differentially regulated during involution [5,8,9]. However, there has been no comprehensive analysis of gene expression in mammary involution.
Microarray analysis has had a major impact on our understanding of the transcriptional basis of complex biological systems. Normal tissue development and homeostasis has been studied in a variety of mouse tissues, including retina [10], liver [11], pancreas [12], uterus [13] and mammary gland [14].
The few microarray studies of normal mouse mammary gland described in the literature have either focused on early stages in the developmental cycle [14], or have used mammary data as a tool for illustrating methods of data analysis [15,16].
In the present study we applied a microarray approach to study the transcriptional expression (mRNA levels) of 8618 mouse genes in mammary gland during the pregnancy cycle. We focused on five time points in involution but we also included earlier time points in the pregnancy cycle to establish the specificity of these involution-related genetic changes. We provide the raw microarray data files and processed genelists for online analysis or download (see the Mammary Apoptosis and Development Group Home Page: www.path.cam.ac.uk/~madgroup). Here we highlight specific aspects of the microarray data that are pertinent to involution by identifying common patterns of gene expression and relating them to specific biological functions. Thus, we demonstrate the sequential activation of death receptor genes followed by components of the mitochondrial (intrinsic) pathway after weaning, suggesting that different apoptosis mechanisms are employed at different phases of involution. We also highlight a prominent role for immune related genes in involution, equating this with phagocytic clearance of apoptotic cells and the maintenance of an antimicrobial environment during milk stasis. We provide evidence that the proapoptotic, acute phase transcription factor Stat3 (signal transducer and activator of transcription-3) [17,18] may simultaneously suppress a potentially damaging cellular infiltrate. A similar study conducted in Balb/c mice and described in the accompanying paper in this issue [19] confirms proinflammatory gene expression in involution. Together, these data provide a useful online resource for studying the activity of specific genes in mammary involution.

Method Mouse husbandry
Six-week-old C57/Bl/6 mice were obtained from Harlan Laboratories (Bicester, UK). At 8 weeks females were either mated or culled for virgin mammary glands. Postpartum, a minimum of six pups per suckling female was ensured by cross-fostering where appropriate. Females used in involution time points were force weaned at 10 days of lactation. Lactating and involuting mammary glands were monitored daily for signs of localized inflammation or mastitis. Two out of 90 animals were excluded from the study for this reason. The study was approved by the local ethics committee and complies with the Helsinki Declaration.

Harvesting mammary gland RNA
A total of 12 stages of adult mouse mammary gland development were selected for this study: 8-week-old virgin; 5 days, 10 days and 15 days of gestation (in which day 1 was the first day postcoitum); 0 days (first day postpartum), 5 days and 10 days of lactation; and 12, 24, 48, 72 and 96 hours after forced weaning. All mammary glands were harvested between 11:00 and 13:00 to limit circadian effects. A total of six animals were used per time point (three per hybridization). In each case a single abdominal gland was removed following excision of lymph nodes and immediately frozen in liquid nitrogen. RNA was extracted from frozen tissue using Trizol reagent (Invitrogen, Paisley, UK) followed by additional column purification (Rneasy; Qiagen, Crawley, UK). Briefly, three samples of mammary tissue from the same time point were ground under liquid nitrogen, and 20 mg of each was pooled into 1.5 ml Trizol reagent. Resuspended RNA samples were then passed through purification columns according to the manufacturer's instructions and RNA integrity was monitored with Lab-on-a-chip (Agilent, West Lothian, UK) before use in microarray analysis.

Microarray hybridization and data analysis
RNA 5 µg was labelled according to manufacturers protocols (Affymetrix, High Wycombe, UK) and hybridized to Affymetrix MGU74ver2a chips representing 12 488 transcripts or 8618 genes. Computation of expression values was performed using the perfect match/mismatch model implemented in dChip [20].
Twenty-four Affymetrix hybridisations (two for each mammary time point) were selected following removal of outliers identified by dChip. All subsequent analyses were performed with Genespring analysis software (Silicon Genetics, Redwood City, CA, USA). Two filters were applied to the complete data set of 12 488 transcripts. First, 755 transcripts (6%) were removed following elimination of genes whose mean signal intensities were below 20 normalized intensity units in at least 11 of the 12 time point conditions. Second, 5149 transcripts (44%) were removed on the basis that they exhibited no significant change in mean signal intensity across all time points (Welch ANOVA P = 0.05, using the Benjamini and Hochberg False Discovery Rate controlling procedure). A separate lactation/involution specific time course gene set was established by applying the same cutoff filters to all 12 488 transcripts using 5 day and 10 day lactation time points, and 12, 24, 48, 72 and 96 hour involution time points.
Both data sets were clustered (K means) according to expression pattern similarity across all 12 conditions (complete time course) or seven conditions (involution time course) using standard correlation. Nine K mean classifications were performed for each data set, starting with 10 clusters and increasing by five incrementally to 50 clusters. The optimal number of clusters was determined empirically based on highest observed variability and redundancy between similar clusters.
All genes were annotated according to known function using the Gene Ontology Consortium categories [21]: biological process, cellular component and molecular function. Onto-Express [22] was used to determine whether clusters of genes with similar expression profiles were enriched in specific GO functional categories. Based on the genes present on the GeneChip, Onto-Express calculated the expected number of occurrences of each functional category in each cluster. The probability that each functional category was over-represented in a cluster was derived using a binomial model and the P values were corrected for multiplicity using the Bonferroni method. Thus, Onto-Express provided information about the statistical significance of each of the pathways and categories represented by the genes in each cluster (these data are available online as Additional files 1 and 2).

Archived data
All genelists, including expression profile clusters, ontological definitions and all raw data .cel files, fulfilling MIAME criteria, were uploaded to ArrayExpress (http://www.ebi. ac.uk/arrayexpress/) and are also available through GeNet (http://genet.hgmp.mrc.ac.uk:8080/servlet/GeNet) for download and online analysis. (For information on accessing data through GeNet or from Breast Cancer Research see our website: www.path.cam.ac.uk/~madgroup/). Gene lists use Affymetrix probe_ID and accession number as transcript identifiers. In the body of the text and figures, gene/transcripts are named according to the Mouse Genome Database [23].

Results
We analyzed the expression of over 12 000 transcripts, or approximately one-third of the protein encoding capacity of the mouse genome, at 12 different time points in the mammary gland pregnancy cycle. We focused our analysis on the changes in gene expression that occur during the first 4 days of postlactational regression, harvesting at five time points during this period. We also included three time points in lactation, three in gestation and a single time point in nulliparous animals in order to distinguish between involution-specific genes and genes that were also highly expressed at other times in the pregnancy cycle. All raw and normalized data files were submitted to GeNet and made available for download or online analysis through our website (www.path.cam.ac.uk/~madgroup/). Expression of 6796 transcripts (54% of the transcripts assayed) significantly changed across the 12-point time course. This suggested that approximately half of the protein encoding content of the mouse genome was regulated in mammary gland at some point during the pregnancy cycle. Among the genes that did not significantly change (nonsignificant profiles genelist in GeNet) were a number of recognized 'housekeeping' genes, including 18s ribosomal and α-actin genes. However, other genes that are ubiquitously expressed in some cell types (e.g. GAPDH and cyclophilins) exhibited differential expression in the mammary time course. This reflects a combination of factors, including changes in tissue composition and dramatic shifts in cellular metabolism and macromolecular synthesis throughout the mammary cycle, in particular during lactation.
In order to verify that the microarray time course reflected expected global expression patterns in vivo, we compared our results with the known mammary expression patterns of 32 genes described in the literature, all of which exhibited the expected expression profiles (see Additional file 3).

Expression profile analysis of the pregnancy cycle (12-point time course)
It was expected that a proportion of the 6797 differentially expressed transcripts would be specific to involution, whereas others may have additional roles at earlier stages in the pregnancy cycle. In order to identify these subsets of involution related genes, the expression patterns of all differentially expressed transcripts were compared and subsequently grouped according to similarity into 35 gene profile clusters. Some clusters exhibited superficially similar patterns (Fig. 1). In order to simplify the classification of expression patterns observed in the time course, these similar profiles were pooled, resulting in 11 broadly distinct expression profiles (Fig. 1a).
These 11 clusters represented maximal gene expression at different stages during the pregnancy cycle. Thus, one cluster of 651 transcripts (5% of all transcripts on the array) were maximally expressed in virgin mammary gland (V SL ). Two clusters with a total of 2317 transcripts (19%) were maximally expressed in gestation (G and G SL ) and three clusters (529 transcripts, 4%) were maximally induced during lactation (L T , L and L G ). Four clusters (3113 transcripts, 25%) exhibited maximum expression during involution (I L , I T , I P and I G ) whereas one cluster Available online http://breast-cancer-research.com/content/6/2/R92

Figure 1
Clusters of gene expression profiles from a 12-point time course of adult mammary gland development. Thirty-five K means clusters were obtained from 6796 differentially expressed genes across a 12-point time course of the adult mammary gland pregnancy cycle. In order to identify basic trends in expression, these clusters were arbitrarily organized into groups based on the similarity between their mean profiles (b). These clusters were pooled and displayed as combined average expression profiles (a). (a) Combined expression profiles named according to developmental stage exhibiting maximal expression (capital) and their associated expression pattern (subscript). V SL , virgin and suppressed lactation; G, gestation; G SL , gestation and suppressed lactation; L T , lactation, transient; L, lactation; L G , lactation and gestation; I L , involution and lactation; I T , involution, transient; I P , involution and parturition; I G , involution and gestation; PC, postcoitum. Y axis is the mean normalized signal intensity. Dotted lines indicate boundaries between major phases of development: virgin to gestation; parturition; and 10-day lactation (forced weaning). Total number of genes in each combined cluster shown to the right. (b) Thirty-five K means clusters. Total number of genes (brackets) and the ID for each cluster is shown. Each differentially expressed transcript is represented only once in the 35 clusters. Error bars represent standard error of the mean. X axis is the 12 developmental time points used in microarray analysis: v, virgin (8 week); G, gestation (days 5, 10 and 15 postcoitum); L, lactation (days 0, 5 and 10 postpartum); I, involution (days 0.5, 1, 2, 3 and 4 after forced weaning).
included 187 transcripts (1%) that were uniformly expressed throughout the pregnancy cycle (PC). Further details on the nomenclature of these combined clusters are provided in the legend to Fig. 1. Genelists for each combined cluster and the 35 primary clusters are available online (see Method section, above).
Six clusters consisted of transcripts that were either maximally expressed during involution (I L , I T , I P and I G ) or specifically suppressed following weaning (L and L G ), suggesting that a proportion of these genes may play important roles in regression of the mammary gland. These six clusters consisted of 3463 transcripts, or 28% of all transcripts on the array. Only one group of 518 genes exhibited involution-specific expression, all of which were transiently activated within 24 hours after weaning (I T ). It was expected that a proportion of these represent genes that may contribute to the initial phase of cell death characteristic of the reversible first stage of involution. All other involution related clusters exhibited elevated expression at earlier stages in the pregnancy cycle. For example, one group of 236 involution related genes (I P ) exhibited an additional peak of expression specifically at parturition. This suggested that many transcriptionally regulated processes in involution were common to gestation (I G ), parturition (I P ) or lactation (I L ).

Correlation of gene function with gene expression pattern
Of the 6747 differentially expressed transcripts, 3273 (49%) encoded proteins with known function. In order to establish what proportion of the transcripts in each cluster from Fig. 1a shared common biological functions, transcripts with known function were assigned to categories as defined by the Gene Ontology Consortium (www.geneontology.org), and statistically significant associations between gene function and each cluster were identified using Onto-Express [22]. The significant associations (corrected P < 0.05) are listed in Table 1. Further details on the distribution of functionally related genes between these clusters are available in Additional file 4.
Unexpectedly, similar proportions of apoptosis related genes were expressed in each of the 11 clusters; the only significant correlation being with the I L cluster, in which maximal gene expression occurred at 12 hours of involution. However, a number of other biological processes were significantly linked to particular stages of mammary development, and were consistent with our current understanding of the physiology of the mammary gland. Thus, during gestation there was a significant increase in the expression of cell cycle regulatory genes (G, G SL and I G ) and a concomitant increase in nucleic acid and macromolecular synthesis ('biogenesis' and 'metabolism' in Table 1). By lactation, the proportion of cell cycle genes had diminished, to be replaced by an increase in the Breast Cancer Research Vol 6 No 2 Clarkson et al. Table 1 Relationship between gene expression pattern and gene function in the 12 Transcripts representing genes with known function were classified by molecular function and biological process, according to the ontological definitions of the Gene Ontology Consortium. The total numbers of known transcripts classified in this way are shown for each cluster. A statistical test for nonrandom distribution of these functional groups (taking into account the total number of transcripts in each category) was performed on each cluster using Onto-Express. All statistically significant associations (corrected P < 0.05) are listed for each cluster using broad ontological definitions. The genes making up these broad categories are provided on GeNet (see Method section). Additional information on the details of these significant associations and a graphical representation of all ontologies are available as Additional files 1 and 4. ECM, extracellular matrix. number of genes involved in development and differentiated function, including fatty acid biosynthesis and other metabolic processes (L G , L and I L ). These three clusters should therefore contain genes that are important for terminal differentiation and transition to a secretory phenotype (genelists for each cluster are available online; see Method section, above).

R96
Morphogenesis and tissue remodelling are characteristic features of the mammary gland both during pregnancy and late in postlactational regression. We observed significant correlations between morphogenesis and three clusters with broadly similar expression profiles (L G , L and I L ), each being characterized by increased gene expression during gestation and lactation followed by a marked and progressive decline at the lactation-involution boundary. A significant proportion of remodelling-related genes were also observed in the V SL cluster, which exhibited an inverse pattern to the three profiles described above. Because the transcripts in the V SL cluster were, by definition, distinct from those in the L G , L and I L groups, it is possible that the former may contain remodelling related genes that contribute to mammary regression whereas the latter group of three clusters contains genes that are involved in invasion and remodelling of the fat pad during pregnancy.
A strong statistical relationship existed between involution and immune-related genes ( Table 1 and Fig. 2). Of the four clusters exhibiting maximal gene expression during involution, three were significantly associated with inflammation (I T ), the acute phase response (I P and I L ), or humoral immunity (I P ). Of the remaining eight clusters, V SL , whose transcripts exhibited a modest increase during involution (and were linked with remodelling), was associated with innate cellular defence (Fig. 2). Immune genes were most prevalent in the I P cluster, in which transcripts were transiently expressed at parturition and maximally expressed by 48 hours of involution (Fig. 2). Of the 57 immune-related genes present in the I P cluster (24% of known genes in the cluster), 43 (18%) were immunoglobulin genes, seven (3%) were acute phase protein genes and seven (3%) were inflammatory mediators.

Acute phase gene expression correlates with Stat3 activity
We have previously shown that the IL-6 family cytokine leukaemia inhibitory factor (LIF) induced phosphorylation (and activation) of Stat3 in mammary gland during involution, and that this transcriptional signalling pathway was important for the induction of apoptosis in the early involuting gland [17,24]. Stat3 is also a recognised mediator of Available online http://breast-cancer-research.com/content/6/2/R92

Figure 2
Relative expression of immune related genes in the 12-point timecourse of mammary development. Immune-related transcripts differentially expressed in the 12-point developmental time course were grouped into five categories that represent different immunological functions. The number of transcripts in each category is shown in brackets. Genelists of all transcripts assigned to these immunological categories are available online. (a) The frequency of transcripts in each expression profile cluster for each immunological category. Statistically significant associations, as determined by Onto-Express, are indicated with asterisks. Y axes indicate the percentage for each cluster. X axis indicates the combined expression profile clusters. (b) The mean expression profiles for all transcripts in each immunological category. Y axis indicates the mean normalized signal intensity, and the X axis indicates the 12 developmental time points (see Fig. 1). the acute phase response, and because acute phase related genes exhibited a parturition/involution expression profile (Fig. 2) we assessed the relationship between Stat3 activity and acute phase gene expression.
Acute phase proteins are classified according to the transcriptional signalling pathways that regulate their expression [25]. Thus, class I acute phase genes are activated by cooperative binding of CAAT-enhancer binding protein (C/ebp)β, C/ebpδ and/or nuclear factor (NF)-κB (and in some cases also Stat3) to specific sites in their promoters. Class II genes are activated by Stat3 alone. Class I genes were induced at parturition and involution (similar to the I P profile). Class II genes were transiently, and specifically, expressed at the onset of involution (similar to the I T profile). We previously described Stat3 and NF-κB activities in the mammary pregnancy cycle [24,26], whereas Cebp gene expression in mammary gland has also been described [27]. We compared the expression of these transcription factors and found that Stat3 and both Cebp genes in addition to phosphorylated Stat3 DNA binding activity [17] correlated with class I acute phase gene expression (Fig. 3).

Expression profile analysis of lactation/involution (seven-point time course)
In order to define patterns of gene expression pertinent to the lactation-involution transition, irrespective of expression in pregnancy, differentially expressed genes were reclustered using a truncated time course involving lactation and involution time points only (Fig. 4). A total of 4103 transcripts (33% of all transcripts) exhibited significant changes in expression across the seven time points, namely 5 days and 10 days of lactation, and 12, 24, 48, 72 and 96 hours of involution. These profiles were sorted according to similarity into 15 clusters (Fig. 4b). Similar clusters were then combined to produce seven basic profiles that differed according to the kinetics of gene induction or suppression during lactation and involution (Fig. 4a). Thus, four patterns of activation were observed: a rapid but transient increase in gene expression 12 hours after weaning (Inv1); a rapid activation, which was maximal by 12 hours and sustained for up to 4 days (Inv2); a slightly delayed induction, maximal by 24/48 hours, and sustained for at least 96 hours after weaning; and a gradual increase in expression, which peaked at or beyond 96 hours (Inv4). Similarly, three patterns of suppression were observed: a transient decrease in transcript levels that mirrored the pattern observed in Inv1 (Lac1); a rapid loss of expression, reaching baseline levels at 12/24 hours and remaining low for up to 96 hours (Lac2); and a delayed reduction in expression, originating from a peak in expression at 12 hours involution and reaching a minimum beyond 96 hours (Lac3).
Consistent with our analysis of the complete 12-point time course, we found that involution related gene expression was associated with an overall increase in immune gene expression (Table 2 and Fig. 5; also see Additional file 5). There was an initial transient increase in proinflammatory genes (Inv1). However, the association was most pronounced with transcripts exhibiting delayed kinetics (Inv4), in which known genes involved in innate immunity, antimicrobial defence and inflammation were all preferentially expressed. The delayed immunological changes also coincided with a significant increase in apoptosis markers (see below) and factors that affect tissue architecture, including genes encoding structural proteins (decorin, laminins, Breast Cancer Research Vol 6 No 2 Clarkson et al.

Figure 3
Mean expression profiles of acute phase genes. (a) Expression profile of the signal transducer and activator of transcription (Stat)3 transcript (red), as compared with CAAT-enhancer binding protein (c/ebp)β (dark blue) and c/ebpδ (light blue) from the microarray. (b) The mean expression profiles of 12 class I acute phase response genes (cI APR; red) and five class II acute phase response genes (cII APR; blue). Y axes indicate mean normalized intensity, and X axes indicate the developmental time points in days. Error bars represent standard error of the mean. Gest, gestation; Lac, Lactation; Inv, involution.
Nid1, collagen 1α1, 1α2, 3α1 and 4α2) and matrix associated proteins (Mmp3, Mmp12, Gelsolin, uPa and Adamts1), thus corroborating previous reports of the role of matrix proteases in remodelling of the regressing mammary gland [2]. Lactation specific expression, as expected, was associated with downregulation of differentiation and metabolic processes after weaning (Lac2 and Lac3; Table 2).

Cytokine expression highlights two immune phases during involution
In order to gain further insight into the immune responses associated with mammary involution, we listed the cytokines, chemokines and immune-related genes associated with each involution profile from Fig. 4 (Table 3). This confirmed the two phases of immune-related gene expression illustrated by Fig. 5: an initial transient expression of Available online http://breast-cancer-research.com/content/6/2/R92

Figure 4
Clusters of gene expression profiles from a seven-point time course of mammary lactation and involution. Fifteen K means clusters of 4103 differentially expressed transcripts from seven stages of lactation and involution. Fourteen primary clusters (b) were grouped according to similarity into (a) seven basic profiles: Inv1, Inv2, Inv3, Inv4, Lac1, Lac2 and Lac3. These similar clusters were pooled and plotted as seven combined expression profiles. The 15th cluster represented differential expression between 5-day and 10-day lactation only (not shown). The number of genes for each pooled cluster is indicated to the right of panel A. The number of genes (brackets) and the cluster ID is indicated for each primary cluster in panel B. Each differentially expressed transcript is represented only once in the 15 primary clusters. The number of genes in each combined profile (panel A) is the sum of its primary profiles (panel B). Y axis is the mean of normalized intensities. Error bars represent the standard error of the mean for all transcripts in each cluster. X axis = 7 point timecourse across lactation (Lac) and involution (Inv). Dotted lines highlight proinflammatory cytokines and TNF superfamily genes followed by delayed expression of monocyte and lymphoid chemokines and immunoglobulin genes.
Thus, an initial burst of cytokine and cytokine receptor expression was observed within 12 hours of forced weaning (Inv1; Table 3), which included the acute-phase mediators LIF receptor (Lifr), IL-11 (Il11) and the neutrophil chemoattractant Gro1 (Cxcl1), and the proinflammatory mediators IL-1a, IL-1b and IL-13 (Il1a, Il1b and Il13). This was accompanied by a transient decrease in other cytokine receptors (including Il11r) and an inhibitor of IL-1 receptor signalling (Il1rn), suggesting that specific proinflammatory signals may be suppressed at the transcriptional level. R100 Table 2 The relationship between gene expression pattern and gene function in the seven-point lactation/involution time course Transcripts representing genes with known function were classified by molecular function and biological process, according to the ontological definitions of the Gene Ontology Consortium. The total numbers of known transcripts classified in this way are shown for each cluster. A statistical test for nonrandom distribution of these functional groups (taking into account the total number of transcripts in each category) was performed on each cluster using Onto-Express. All statistically significant associations (corrected P < 0.05) are listed for each cluster using the same ontological definitions used in Table 1 and available on GeNet. Additional information on the details of these significant associations and a graphical representation of all ontologies are available as Additional files 2 and 5. ECM, extracellular matrix.   Inv1  Inv2  Inv3  Inv4  Lac 1  Lac2  Lac3   Cytokines  Lifr  Il11ra2  Il12rb1  Osmr  Il1rn  Il1rl1l  Il1r2  Il1a  Ltbr  Il1r1  Il4ra  Il1b  Il10rb  Il7r  Il1r2 Csf1r The identity of immune-related transcripts was determined for each of the combined profile clusters illustrated in Fig. 4. Thus, transcripts belonging to the ontological category Immunity from each cluster were listed according to the immunological subcategories shown in the first column of the table. All gene names conform to the mouse genome informatics database [23], synonyms are included in brackets where appropriate. The number of transcripts of the same immunoglobulin chain genes are indicated. Complete genelists for these immunological categories are available on GeNet (see Method section). marked increase in the monocyte/macrophage markers Lrp, Csf1r and Cd14 (Table 3 and Fig. 6a), but no evidence for an accompanying increase in neutrophils or eosinophils (polymorphonuclear leucocyte markers Fut4, Fcgr3, Caecam1, Il8rb and Tnfrsf5 were not induced during involution; data not shown). This is consistent with previous observations of the immune cell complement in regressing murine mammary gland [28] and supports the hypothesis that neutrophil numbers are suppressed in involuting tissues in order to prevent subsequent tissue damage [28,29]. Interestingly, the gene for oncostatin M receptor (Osmr) was upregulated in involution (Fig. 6b). In models of inflammation oncostatin M is produced locally to prevent proinflammatory cytokines such as IL-8 from recruiting neutrophils, thus attenuating the inflammatory response [30]. We also confirm a previous report that uterocalin (Lcn2) and Sgp2, a subunit of clusterin, are upregulated during involution and parturition [19], because their expression correlates with the I P profile (Fig. 6b). Uterocalin, an acute-phase response gene that is regulated by Stat3 [31,32], prevents neutrophil accumulation possibly through targeted destruction of neutrophils [29], whereas clusterin protects cells from the damaging effects of neutrophils [33,34]. Furthermore, secretory leukocyte protease inhibitor (Slpi), a potent inhibitor of neutrophil proteases [35], and a Stat3-induced gene (Boland M and coworkers, unpublished data) was also markedly upregulated during involution (Fig. 6b). Thus, although the transcriptional effects of the acute-phase mediators LIF-Stat3 and C/ebpβ/δ are clearly evident (Fig. 3), the immunological outcome at the cellular level is likely to be minimised by coexpression of anti-inflammatory factors such as these.

Identity of immune-related genes in lactation and Involution clusters
Following the initial 'spike' in inflammatory cytokine expression, six different chemokine genes were slowly induced over the 4-day period of involution (Inv4; Table 3). These included chemoattractants for monocytes and macrophages (Ccl6, Ccl7, Ccl8 and Cxcl14) and lymphocytes (Ccl6, Ccl21a, Cxcl9 and Cxcl14). This chemokine induction was accompanied by a marked increase in immunoglobulin gene expression (Fig. 5, Table 3). Consistent with these observations, the monocyte/macrophage markers Lrp1, Csf1r and Cd14, and the lymphocyte antigens Tnfrsf9, H2-M1 and Tnfrsf4 were increased during involution (Table 3 and Fig. 6a).

Soluble innate defence factors are induced in involution
In addition to the immunoglobulins, genes that encode a number of innate soluble defence factors were highly expressed in involution (Table 3 and Fig. 6c). These soluble factors possess antimicrobial properties that are secreted to protect against mastitis during the vulnerable period of milk stasis [36]. Thus, the iron chelator lactoferrin (Ltf), which is found at high levels in the milk of humans and cattle [37,38], is induced with kinetics similar to those of the antibody transcripts described above (Fig. 5). Comple-Breast Cancer Research Vol 6 No 2 Clarkson et al.

Figure 6
Involution related expression profiles: apoptosis and immune related genes.  (Table 3). These factors are found at high levels in mastitic milk and in normal bovine mammary glands during involution, and are thought to play a dual role in opsonization and lysis of bacteria and as a proinflammatory mediator [39,40]. Cd14 was significantly upregulated in involution and was also highly expressed during parturition (Fig. 6c). This antigen is expressed on bovine mammary polymorphonuclear neutrophils and macrophages [41] and plays an integral role in suppressing bacterial infections in response to endotoxin. Furthermore, Lbp [42] and Ly86 [43], which in conjunction with Cd14 augment responses to lipopolysaccharide (LPS), were also upregulated. The significant regulation of these genes in uninfected mammary gland confirms a role for soluble defence factors in normal development of mouse mammary gland.

Involution is associated with a transient induction of death receptor genes
Involution is characterized by extensive apoptosis of epithelial cells, which occurs both before and during the tissue remodelling process. By subdividing the first 4 days of involution into three discrete phases of gene expression, we demonstrated a link between the expression of apoptosis-related genes and the delayed onset of extracellular matrix degradation and alveolar collapse at 48 hours (Inv3 and Inv4; Table 2). However, we were unable to identify a significant increase in the number of apoptosis genes expressed at the onset of involution. One explanation for this was that specific combinations of cell death genes, rather than an increase in the number of genes activated, influenced apoptosis at the onset of involution. In order to address this, all cell death related genes found to be regulated during involution were grouped according to their lactation/involution expression pattern (Table 4). Regulators (Bax, BclX, Mcl1, mXIAP and Apaf-1) and effectors (caspases) of classical apoptosis and immune cell-mediated killing (granzyme A, granzyme B, Fadd, Fas, Fas ligand) were regulated with different kinetics during involution and were relatively evenly distributed across the four involution clusters (Table 4).
However, we observed a strong bias in the expression of extrinsic (death receptor) and intrinsic (mitochondrial) apoptosis genes, in which extrinsic-related components were generally induced within the first 24 hours of involution whereas intrinsic genes, with the exception of Hrk and Diva, were not maximally transcribed until at least 96 hours.

Inv1
Inv2 The identity of cell death-related transcripts was determined for each of the profile clusters illustrated in Fig. 4. Thus, transcripts belonging to the ontological category apoptosis or cell death from each cluster were listed according to their expected role as proapoptotic or survival factors. All gene names conform to the mouse genome informatics database [23], and synonyms are included in brackets where appropriate. A complete genelists for this apoptosis category is available on GeNet (see Method section).
with previous studies in mammary gland of two of these ligand-receptor pairs [44,45]. NF-κB activity also correlated with the rapid activation of these TNF superfamily ligands. DNA binding activity of this transcription factor is markedly upregulated within 3 hours of forced involution and is suggested to promote survival of a subpopulation of mammary epithelial cells [26]. Thus Nfkb2 and activators of NF-κB (Tnfsf12, Bcl10, Tnfrsf11a, Traf6 and Tnfsf10) were upregulated, whereas an inhibitor of NF-κB function (Tnip1) was inhibited during the initial 24 hours following forced weaning, when NF-κB activity was induced [26]. Previous studies demonstrating the direct activation of NF-κB by TNF superfamily ligands in mammary epithelial cells [26,46,47] therefore suggests that the observed transient increase in these ligands may be functionally relevant.
In contrast to the rapid response of extrinsic factors, intrinsic apoptosis genes were predominantly (but not exclusively) expressed at least 24 hours after forced weaning. Proapoptotic and antiapoptotic members of the Bcl-2 family (Bax, Mcl1 and Bcl2l) were induced together with the downstream effector of the intrinsic pathway Apaf1.
Although the precise mechanism(s) of cell death in mammary gland have yet to be established, the trends in apoptosis gene expression identified here hint that extrinsic (death receptor) and intrinsic death programmes may predominate at different times. A marked induction of Igfbp5 complemented the delayed induction of intrinsic apoptosis genes (Fig. 6d). The kinetics of this potent inducer of mammary apoptosis and involution [48] suggests that it probably contributes to the stimulation of intrinsic cell death at this time.
Apoptotic bodies initially shed into the lumen and the residual dying cells left in the epithelium are believed to be cleared from the mammary gland by nonprofessional phagocytes and interstitial macrophages, respectively [28]. Receptors and ligands that mediate recognition and engulfment of apoptotic cells include members that are common to both types of phagocytic cell (e.g. Cd36, Lrp1 and Mfge8) and those specific to professional macrophages (e.g. Cd14 and Cd68). The genes for each of these factors (except Cd36, which exhibits an initial downregulation before a modest increase by 96 hours) are induced with similar kinetics (Fig. 6e). This is indicative of the increasing demands placed on the mammary gland to clear dead/dying cells during regression, and also probably reflects the increasing numbers of macrophages occupying the gland during the first 4 days of involution.

Discussion
Microarray analysis of large sequence verified cDNA probe sets provides a powerful tool for analyzing complex biological systems because it can extract patterns of gene expression from a significant proportion of the total genomic content of an organism. In this study we simulta-neously analyzed approximately 12 000 transcripts, representing 8618 unigene clusters, or about one-third of the protein encoding capacity of the mouse, to determine basic expression patterns in adult mammary gland during the pregnancy cycle.
Although a previous study of approximately 5000 transcripts described gene expression in early phases of the pregnancy cycle [14], we focused our analysis on the first 4 days after forced weaning, using time points from earlier stages in development to establish the specificity of these involution-related genetic changes. In the accompanying paper, Stein and coworkers [19] describe a similar study in Balb/c animals. Despite differences in the methods of analysis, their findings on the activation of inflammatory mediators and monocyte/lymphocyte markers in involution correlate with the results described here, further verifying our microarray data. Thus, strain differences did not appreciably influence the immune responses associated with involution.
The principal objective of the present study was to provide a comprehensive set of gene expression data from late stages of mouse mammary development to be used as an online resource for further analysis and/or download. The raw microarray data files and processed genelists from the study are therefore available for download or online analysis (see www.path.cam.ac.uk/~madgroup/).
Here we have outlined the common patterns of gene expression across the developmental time course and have correlated these patterns with biological function. Thus, as expected, pregnancy was associated with the induction of genes that are involved in cell proliferation, differentiation and tissue morphogenesis. This was superceded in lactation by genes involved in biosynthetic (fatty acid synthesis) and metabolic pathways, correlating with the findings of an earlier microarray study of the pregnancy cycle [14]. We have not attempted to provide a detailed account of the genetic changes that occur early in the pregnancy cycle; however, the data are available online for further scrutiny.
By concentrating six of our 12 time points around the first 96 hours of involution, we were able to determine the general kinetics of gene expression during this period. Thus, according to the cluster analysis presented in Fig. 4, most differentially expressed transcripts in involution were initially induced (or suppressed) within 24 hours of weaning, whereas maximum expression generally occurred at 12 (Inv1, Inv2), 48 (Inv3) or 96 (Inv4) hours. This implies three phases of gene expression in the first 4 days of involution -phases that are characterized both by their time of initiation and their rate of induction. At the cellular level this period of involution has previously been divided into two stages [2]: an initial reversible period of approximately 48 hours characterized by increased apoptosis and loss of differentiated function, followed by matrix degradation and alveolar collapse. The pattern of gene expression appeared to correlate with this cellular model. We observed a swift increase in specific apoptosis transcripts (Inv1; Table 4) and a subsequent decline in differentiated function by 48 hours (Inv3; Table 2), followed by an increase in factors affecting tissue architecture (Inv4; Table 2) and phagocyte function (Fig. 6). Further detailed analysis of the Inv3 cluster, which exhibits maximal expression at 48 hours, may elucidate underlying mechanisms that are responsible for the transition to phase two involution.

Transcriptional regulation of apoptosis in involution: two phases of death?
Our data also provide provisional evidence for the existence of a two-phase apoptotic process in involution and which also coincides with the two phases previously described [2]. While the number of apoptosis genes was similar in each involution cluster, there was a pronounced bias in the distribution of extrinsic and intrinsic apoptosis genes. The initial spike of activity consisted of four death receptor ligands, namely Tnf (TNF-α) Tnfsf6 (Fas ligand), Tnfsf10 (TRAIL) and Tnfsf12 (TWEAK). These ligands bind to the principal death receptors Fas, TNFR-1, TNFR-2, DR3 and DR4, which initiate apoptosis in cell mediated immune cell death and in a variety of physiological and pathological situations. These data confirm previous observations that Fas and Fas ligand are upregulated within a day of involution in mouse mammary gland, and that the deletion of these genes results in delayed regression and loss of apoptosis [44]. Deregulation of Fas signalling through aberrant Akt and NF-κB activity is believed to play a role in the immune escape of breast cancers and in resistance to chemotherapy [49]. TNF-α also induces NF-κB signalling to prevent apoptosis in a subpopulation of cultured mammary epithelial cells [26,47], and has been shown to be downregulated late in involution (day 7); however, it has not been studied at earlier time points [50]. Furthermore, elevated IL-10 dependent TRAIL expression (and its receptor DR4) has previously been demonstrated in involution [45].
In contrast to the transient activation of death receptor ligands, components of mitochondrial (or intrinsic) apoptosis exhibit sustained, often delayed induction. Thus the caspase 9 effector, Apaf1, and Bax were induced with delayed kinetics. Costimulation of the survival factors Mcl1 and Bcl2l (BclXl) suggests that the transcriptional regulation of Bax may be counterbalanced, and therefore controlled, by expression of Bcl-2 homologues. Suppression of BclX, for example, promotes death during mammary involution [51], whereas delayed induction of Bax has previously been demonstrated [52]. Regulation of Bcl homologues was not exclusively restricted to late stage involution. Hrk and Boo/Diva were differentially expressed within 24 hours. Hrk is a proapoptotic factor that is transcriptionally suppressed by progesterone signals in breast cancer cell lines [53]. The function of Boo/Diva, which has only previously been described in mouse ovary, is poorly defined [54]. Cross-talk between death receptor pathways and Bcl-2 regulatory factors is well established [55] and may explain the early activation of these two Bcl homologues. Probable upstream initiators of a putative delayed apoptotic mechanism include the loss of survival signals, such as insulin-like growth factor (IGF)-1 [48,56], and the destruction of extracellular matrix-integrin contacts leading to an anoikis-type death [6,57]. The dramatic induction of the proapoptotic mammary protein IGFbinding protein 5 proposed previously [48] coincides with the onset of second phase involution (Fig. 6d, Table 4) and the upregulation of metalloproteinases [2,6]. Both IGF-1 and integrins signal through Akt, which was also downregulated in involution ( Table 4). Loss of this signal could then induce Bax [57]. Our data suggest that apoptosis is induced in involution following the specific transcriptional activation of a subset of proapoptotic genes encoding death receptors, and that delayed expression of inhibitors of cell survival (e.g. IGF-binding protein 5) induce alternate apoptosis pathways during the second phase of involution. This hypothesis warrants further investigation.
Integral to the apoptotic process is the phagocytic clearance of apoptotic bodies by macrophages. It has been suggested that monocytic macrophages are recruited to the mammary gland to supercede phagocytic epithelial cells in the clearance of milk and dead cells [28]. We have demonstrated monocyte specific (Cd14 and Cd68) and nonspecific phagocyte receptors during involution, supporting a role for monocytic phagocytes in mammary involution.

Multiple roles for immune cells in mammary involution
Measurements of the cellular content in milk from involuting glands of sheep, pigs and cattle indicate an early rise in polymorphonuclear neutrophils, a subsequent increase in macrophages and low numbers of eosinophils, B cells and T cells [7,28]. However, the immune cell complement of human and rodent mammary gland is poorly defined. It is believed that these phagocytes are recruited to the gland to perform several functions in phagocytic clearance and cellular immunity.
We have described the sequential activation of cytokine and chemokine genes in involution. The effect of these would be to attract monocytes and lymphocytes to the gland, and we identified an increase in a number of immune cell markers, supporting this hypothesis. For example, increased expression of the Csf-1 receptor, c-fms (Csf1r; Fig. 6a), which is restricted to macrophages of the monocyte lineage in mammary gland [58] and is required for appropriate branching morphogenesis in pregnancy [59], confirmed that macrophages are critical for normal development [58]. A role for these cells in phagocytic clearance during involution is emphasized by the concomitant increase in two receptors (Cd14 and Cd91) and the milk factor Mgfe8, which mediate recognition of apoptotic cells by macrophages (Fig. 6c,e).
Recruitment of cells of the monocytic lineage comprises an inflammatory response. However, in the absence of a persistent bacterial infection, this does not occur in the regressing mammary gland, even though inflammatory cytokines and downstream mediators are activated, presumably to perform additional roles during involution.
During chronic inflammatory diseases such as rheumatoid arthritis or inflammatory bowel disease, uncontrolled accumulation of neutrophils yield aberrantly elevated levels of inflammatory compounds, which become major contributors to tissue damage [60,61] Thus, regulation of neutrophil recruitment into inflammatory sites and their clearance are critical processes that ensure effective host defence without tissue injury. In mammary gland, the LIF-Stat3 signaling pathway is known to be crucial for appropriate regression during early involution by initiating epithelial apoptosis [17,24]. With the concomitant downregulation of IL-11 receptor ( Table 3) and lack of a significant increase in IL-6 or its receptor (data not shown), it is also likely to be the predominant protagonist of acute phase signalling. This is confirmed by previous observations that C/ebpβ and C/ebpδ are both transcriptional targets of Stat3 in vivo [24] (Boland M, unpublished data) and that the profiles of Class I acute phase genes correlated precisely with phosphorylated Stat3 activity [24] and Stat3, C/ebpβ and C/ebpδ expression (Fig. 3). Class II acute-phase response genes exhibit specific and transient expression in involution. It is thought that these genes respond transiently to proinflammatory signals [25] and therefore may respond in a similar manner to the rapid activation of Stat3 at the onset of involution. Why they are not induced at parturition is unclear, although a molecular switch that determines the specificity of signalling pathways downstream of LIF has been proposed [24].
Thus, inhibitory mechanism(s) must exist to prevent the accumulation of potentially damaging cellular infiltrates in response to LIF. In a previous study of involution in uterus and mammary gland, the acute phase protein uterocalin (Lcn2) was shown to be induced [29]. It was proposed that one role of this protein in involuting mammary gland and uterus was to eliminate infiltrating neutrophils by inducing apoptosis [62] and therefore to prevent subsequent tissue damage [63]. We are able to confirm significant activation of uterocalin expression in mammary involution (Fig. 6b) and the concomitant absence of neutrophil antigens (Table 3). In conjunction with the induction of uterocalin, clusterin (Sgp2), oncostatin M receptor (Osmr) and secretory leukocyte inhibitory protein (Slpi) were also increased (Fig. 6b). Clusterin protects cells against oxidative and chemical damage [33,34] and has been suggested to be coexpressed with uterocalin as part of a localized acute-phase response in order to limit the damaging effects of residual neutrophils in involuting mammary gland [29,63]. SLPi, an antichymotrypsin, also contributes to the elimination of neutrophil function, whereas Osmr may possess similar anti-inflammatory effects through the sequestration of proinflammatory cytokines. Osmr, Slpi and Sgp2 are all directly induced by Stat3 in mammary epithelial cells (Boland M, unpublished data) [64], whereas uterocalin expression conformed with Stat3 activity (Fig. 3 and 6b). A putative anti-inflammatory role for Stat3 was suggested in our previous study of Stat3-deficient mammary glands [17], in which we observed an increase in mastitis in the absence of Stat3. A similar protective mechanism may exist during parturition, when the LIF-Stat3 signal is also active [24,65]. We are currently investigating the possibility that LIF-Stat3 signalling is responsible for the induction of these four genes in mammary gland, therefore confirming a dual role for LIF in involuting mammary gland. In addition to the direct effects of these anti-inflammatory acute phase genes, immunosuppressive factors including transforming growth factor-β (Inv3; Table 3) expressed by phagocytic macrophages are likely also to contribute to the suppression of an inflammatory response [64,66].

Soluble defence factors are induced during involution
Soluble defence factors and cells that are involved in humoral immunity are thought to be induced to prevent mastitis during this extended period of milk stasis. For example, the iron chelator lactoferrin is induced in mammary involution (Fig. 6c), which plays a role in mammary defence by protecting against coliform infection [36]. Indeed, low levels of lactoferrin is a predisposing factor for mastitis in humans [38]. Immunoglobulins not only have an important role in supplementing milk but are also required to facilitate opsonization by macrophages and neutrophils (IgM and IgG 2 ) [67] and the direct sequestration of bacteria (IgA).
Cd14 exhibited a dramatic pattern of expression throughout the complete time course (Fig. 6c). Cd14 is a core component of innate immunity, which binds LPS and augments LPS responses [42]. It is found bound to intramammary monocytes or neutrophils [41] in cattle subjected to endotoxin, and its soluble form is a known constituent of milk, in which it helps Cd14-negative cells respond to infection [68] and contributes to inhibition of gastrointestinal infections in human infants [69]. It may also play a role in recognition and phagocytosis of apoptotic cells by macrophages expressing cell surface Cd14 [70]. Thus, the role of elevated Cd14 may be both to prevent infection and potentiate apoptotic cell clearance. Binding of soluble Cd14 to LPS is augmented by the acute phase protein Lbp, which also exhibited an I P expression profile (Fig. 6c).

Conclusion
We highlighted two general trends in gene expression during mouse mammary involution: an increase in immunerelated genes and a shift in the identity of apoptosis components. The increase in immune genes correlated with previously described roles for immune cells in involuting mammary gland, namely antimicrobial defence, cytokine signalling and phagocytic clearance. The proposed change in the cell death machinery based on extrinsic and intrinsic apoptosis gene expression has not previously been observed. This putative mechanism clearly needs further investigation. One common link between apoptosis and immune responses in mammary gland is the previously defined role for LIF and Stat3. We are currently investigating the possibility that Stat3 mediates the antiinflammatory effects vital for the protection of the mammary gland in addition to its recognized role as a proapoptotic factor in mammary gland.

Competing interests
None declared.

R107
The following Additional files are available online: