Regulation of milk synthesis and secretion is controlled mostly through local (intramammary) mechanisms. To gain insight into the molecular pathways comprising this response, an analysis of mammary gene expression was conducted in 12 lactating cows shifted from twice daily to once daily milking. Tissues were sampled by biopsy from adjacent mammary quarters of these animals during the two milking frequencies, allowing changes in gene expression to be assessed within each animal. Using bovine-specific, oligonucleotide arrays representing 21,495 unique transcripts, a range of differentially expressed genes were found as a result of less frequent milk removal, constituting transcripts and pathways related to apoptotic signaling (NF-κB, JUN, ATF3, IGFBP5, TNFSF12A) mechanical stress and epithelial tight junction synthesis (CYR61, CTGF, THBS1, CLDN4, CLDN8), and downregulated milk synthesis (LALBA, B4GALT1, UGP2, CSN2, GPAM, LPL). Quantitative real-time PCR was used to assess the expression of 13 genes in the study, and all 13 of these were correlated (P < 0.05) with values derived from array analysis. It can be concluded that the physiological changes that occur in the bovine mammary gland as a result of reduced milk removal frequency likely comprise the earliest stages of the involution response and that mechano-signal transduction cascades associated with udder distension may play a role in triggering these events.
- mammary involution
- mechanical stress
in mammals, reduced suckling frequency initiates the downregulation of milk synthesis and the induction of apoptotic pathways and structural remodeling of mammary tissue. Milk removal frequency appears to regulate milk synthesis through local intramammary pathways, since, if one of two opposing mammary glands is milked less frequently compared with the other, there is a corresponding reduction in secretory activity of the less frequently milked gland only (65, 51). Furthermore, when milk is removed at an increased frequency, milk production increases relative to the control gland (27).
The physiological triggers and discrete signaling pathways that comprise these intramammary responses to milking frequency have not been well defined, although several hypotheses attempt to explain local regulation of milk synthesis (10). The accumulation of compounds in milk that serve as feedback inhibitors of secretion (40) has been suggested. Recent evidence suggests that the neurotransmitter serotonin may actively regulate milk secretion in the mouse (54) and cow (19), and a casein-derived peptide has also been hypothesized as a putative regulator of lactation in the cow (46). Another hypothesis proposes that apoptotic signaling cascades are initiated by cell stretch-sensing mechanisms resulting from the physical distortion of the engorged mammary gland (39).
In cattle, distension of mammary alveoli is essentially complete 16–18 h postmilking (10), and at this time a number of physiological events occur including a progressive reduction in the rate of milk secretion (10). This is manifested in a daily yield reduction of 20–30% in once daily milked cows relative to twice daily milked animals (10), and this acute downregulation of milk secretion is compounded over time through an acceleration of udder involution (5, 31). In an involution time series experiment, histological analysis revealed apoptosis as a relatively early event in response to milk accumulation, with a numerical increase in mammary epithelial and luminal cell apoptosis shown 36 h after cessation of milking and significant increases following 72 h of milk accumulation (47). A significant increase in milk lactoferrin concentration is also seen during less frequent milking (13), which is also a characteristic of early mammary involution (64, 29). These observations suggest that very early involution processes are activated in the glands of cows milked once daily, with mammary tissue largely “rescued” from progression of involution through the daily removal of milk.
To gain insight into early-stage involution mechanisms and to understand the key molecular changes accompanying milk yield reduction, the current study examined changes in mammary gene expression in response to 26 h of milk accumulation, sampled on a background of 5 days of once daily milking. In particular, we aimed to test the hypothesis that a reduction in milking frequency would result in gene expression changes indicative of apoptosis and the mechanical and chemical feedback inhibition mechanisms that have been proposed previously as the earliest signals of milk secretion inhibition.
Animals and Treatments
Approximately 300 Holstein Friesian-Jersey crossbred dairy cows in midlactation were milked once daily for 5 days to identify animals that showed either a minimal or high level of milk production loss under the once daily milking system. Milk yield was recorded before, during, and after the challenge on a daily basis (see Fig. 1).
Twenty cows were selected from this herd for the microarray study (∼4 wk later), representing animals that lost <15% of milk yield (n = 10) or >30% of milk yield (n = 10) during the initial once daily milking period. Further criteria for selection of cows included body condition score [>4.5 on a 1–10 scale (43)], milk yield when milked twice per day (>7.0 kg/day), pregnancy status, uniform udder (no blind teats, dry quarters), and lack of apparent intramammary infections in the current lactation. Animals then underwent an additional once daily milking challenge, with tissue biopsies collected as described below.
Milk yields of the cows selected for biopsy are shown in Fig. 1. Of the 20 cows sampled, milk yield data were missing for a single cow and mRNA was of marginal quality for two other cows. For the remainder, microarray analyses were conducted on a random selection of 12 of these cows.
The trial was approved by the Ruakura Ethics Committee, Hamilton, New Zealand. All animals were handled in accordance with the principles set out in the Declaration of Helsinki and the American Physiological Society's Guiding Principles in the Care and Use of Animals. Selected cows were grazed together on rye-grass/white clover-dominant pasture and milked twice daily, with milk yield data recorded for 6 days prior to the initial mammary biopsy. Milking times were 06:30 and 15:30 during twice-daily milking and 07:30 during the once-daily period. On the day of the first biopsy, cows were milked as usual and the biopsy was taken from a rear quarter of each animal between 09:30 and 11:00. Cows were milked that afternoon to remove residual blood clots, and then the once daily regime was imposed the following day for 5 days. On the sixth day a repeat biopsy was taken from the contralateral gland of each animal between 09:30 and 11:00, ∼26–27 h after the last milking. Twice daily milking resumed that afternoon.
The biopsy site was clipped, cleaned, and scrubbed with iodine antiseptic. Local anesthetic (2% lignocaine, 4 ml) was injected subcutaneously at the biopsy site. A small (2–3 mm) stab incision was made through the skin with a scalpel blade, and the BARD biopsy needle was then inserted into the mammary tissue to a depth of ∼5 cm. The biopsy gun was then fired to take the tissue sample (30 mg), and biopsies were visually inspected. On occasions, repeat samples were taken from the same site. Tissue samples were placed in Nalgene 2.0 ml cryogenic vials (Nalgene, Rochester, NY) and snap-frozen in liquid nitrogen. Firm hand pressure was then applied to the biopsy site to limit bleeding, and the site closed with a michel clip.
Prophylactic treatment with systemic antibiotic was then applied (Excenel, Pharmacia). All cows were hand stripped before and after the subsequent three milkings to ensure removal of blood clots from the gland and teat cisterns and were monitored for any signs of mastitis (none were evident).
RNA Extraction and mRNA Analysis
Isolation of RNA.
Mammary tissue biopsies were homogenized in Qiagen buffer RLT (QIAGEN, Hilden, Germany) using Fastprep Lysing Matrix D tubes in a FastPrep instrument (MP Biomedicals, Solon, OH). Total RNA was extracted using a Qiagen RNeasy kit (Qiagen). RNA quantity was determined by spectrophotometry in a Nanodrop ND-1000 (Nanodrop Technologies, Wilmington, DE). RNA integrity was checked with the Agilent 2100 Bioanalyzer with a RNA 6000 Nano LabChip kit (Agilent Technologies, Palo Alto, CA). RNA samples that had a poor 18S:28S ribosomal RNA ratio indicating degradation were not used for microarray analysis.
For bovine microarray development, 325,000 expressed sequence tags (EST), mRNA, and CDS region sequences for Bos taurus genes were retrieved from GenBank and supplemented with paired 5′- and 3′-reads from 100,000 ESTs obtained from normalized bovine fetal brain and lactating and nonlactating mammary gland cDNA libraries. After repeat masking, sequences were clustered by pairwise comparison using Paracel Transcript Assembler (Paracel, Pasadena, CA) to reduce redundancy, and the resulting clusters were mapped to genes predicted from the 2.1 assembly of the bovine genome. Gene-specific cluster sequences were then supplemented with novel transcript sequences predicted by the NCBI RefSeq project (http://www.ncbi.nlm.nih.gov/RefSeq/). Probe sequences (60-mer) targeted toward the 3′-regions of the transcripts were designed by Agilent Technologies. We selected 21,495 unique probe sequences for array manufacture by in-situ synthesis. To facilitate regulatory network and pathway discovery, array probes were annotated with identifiers for the human orthologs of bovine genes represented on the array. Overall, the arrays represent the bovine orthologs of ∼19,000 unique human genes. Detailed information about the array probes is available on the Agilent Technologies website (http://www.chem.agilent.com).
Production and Labeling of Amplified RNA and Microarray Hybridization
Amplified RNA (aRNA) for microarray hybridization was prepared using the Amino Allyl MessageAmp aRNA Kit (Ambion, Austin, TX) in accordance with the protocol provided with the kit. We used 1 μg of total RNA to generate the amino allyl modified aRNA. Modified aRNA (5 μg) was then vacuum dried and labeled with Cy3 and Cy5 NHS ester dyes (Amersham Cy3 and Cy5 Mono-Reactive Dye Packs; GE Healthcare UK, Little Chalfont, Buckinghamshire, UK). Purified, dye-incorporated aRNA was then quantified and assessed for dye coupling efficiency using the Nanodrop (Nanodrop Technologies) instrument. Microarray hybridization was performed using the Agilent Gene Expression Hybridization Kit in accordance with the Agilent 60-mer oligo microarray processing protocol, version 2.0 (Agilent Technologies). Fragmented, dye-labeled aRNA was added to each Agilent 22K 60-mer oligonucleotide microarray, hybridized overnight (17 h), and washed and dried the following day according to the same protocol guidelines. Arrays were scanned immediately after this drying step using an Agilent DNA microarray scanner.
A total of 12 microarrays were used in this study. Each array was hybridized with RNA derived from a mammary sample taken at each biopsy time point, such that the twice-a-day milked sample and once-a-day milked sample were contained on the same slide for each animal. Six of the microarrays represented animals previously shown to demonstrate a “large” reduction in milk yield, whereas the other six arrays represented animals that previously displayed a “small” reduction in milk yield when milked once per day (Fig. 1).
Data Analysis and Statistics
Data analyses were performed with Genespring GX 7.3.1. (Agilent Technologies). Array data were imported into Genespring using Agilent's two-color “Enhanced FE” import scenario, which included “Per Spot: Divide by control channel” and “Per Chip: Normalize to 50th percentile” normalization steps. A dye swap normalization was also carried out on the six arrays that used reverse signal and control channel dye labeling. Filters were then applied to the data to improve the quality of the normalized dataset. First, data were filtered “on flags” such that any probes that were not deemed “present” or “marginal” (according to feature extraction spot quality guidelines) in at least six of the 12 arrays were omitted from analysis. Secondly, probes that did not have a minimum threshold of 100 raw intensity units in at least six of the 12 arrays were also omitted from analysis. These two quality control filtering criteria reduced the number of probes to 20,416, and this list was then used to generate subsequent probe lists of interest. Probe lists were derived from volcano plots to give one list representing up- and downregulated genes across all animals, and one list representing genes associated with milk yield loss. Animals were grouped into large and small yield loss groups according to their milk yield performance recorded over the biopsy period (Fig. 1) (as opposed to performance during the selection trial), each consisting of seven and five animals respectively. Volcano plot filtering criteria consisted of a P value threshold of 0.05 and an expression value cut-off of 1.5-fold greater or less than control. This expression value was an average across all 12 arrays in the case of the “all arrays” upregulated and downregulated probe list, and an average across the five and seven small and large yield loss groups respectively for the yield loss associated list. The data discussed in this publication have been deposited in the National Center for Biotechnology Information's Gene Expression Omnibus (12) and are accessible through GEO Series accession number GSE16876 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE16876).
Gene Function, Network and Pathway Identification and Analysis
Differentially regulated genes were annotated with broad-range biological functions using the PANTHER (Protein ANalysis THrough Evolutionary Relationships) Classification System (http://www.pantherdb.org/) (55).
Biological network generation and functional analyses were carried out through the use of Ingenuity Pathways Analysis (IPA; Ingenuity Systems, www.ingenuity.com). The two gene lists generated by volcano plot filtering were combined to create a single dataset of sufficient size to allow meaningful pathway assignment. In this way, a master list comprising genes representing expression changes across all animals, as well as differences between the two yield loss groups, was generated. The list of human homologs corresponding to these bovine genes was then uploaded into the application. Each gene identifier was mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base (IPKB), with a subset of these used for network generation based on the information present in the Ingenuity molecular network knowledge base (127 of 176 mapped identifiers eligible for network construction). Global molecular networks of these “focus genes” were then generated based on their connectivity. A score, which is based on a P value calculation and calculates the likelihood that these genes are part of a network and do not associate due to random chance alone, was assigned to each network. This numerical network score was used to rank the networks. Although IPA has the provision to assemble networks and associations via indirect associations, only molecules that had been demonstrated to associate directly with one another were used in the analysis.
IPA Functional Analysis was also used to identify the biological functions that were most significant to this data set. Fisher's exact test was used to calculate a P value determining the probability that each biological function assigned to the data set was due to chance alone.
Canonical pathways were generated through the use of IPA. This analysis identified the pathways from the IPA library that were most significant to the data set. Genes from the dataset that were associated with a canonical pathway in the IPKB were considered for the analysis. The significance of the association between the data set and the canonical pathway was measured in two ways: 1) A ratio of the number of genes from the data set that map to the pathway divided by the total number of genes that map to the canonical pathway is displayed. 2) Fisher's exact test was used to calculate a P value determining the probability that the association between the genes in the dataset and the canonical pathway is explained by chance alone.
Quantitative Real-Time Reverse Transcription Polymerase Chain Reaction
Thirteen genes spanning both gene lists were subjected to quantitative real-time PCR analysis. RNA samples were treated with DNase I Amplification Grade (Invitrogen, Carlsbad, CA) and purified with Qiagen RNeasy columns (Qiagen). Superscript III reverse transcriptase and random hexamer primers (Invitrogen) were used to generate cDNA from 5 μg of total RNA as per the Invitrogen Superscript III protocol (Invitrogen part no. 18080.pps). FAM-labeled Taqman probes were designed using the ABI Filebuilder program (www.appliedbiosystems.com/filebuilder), with Refseq mRNA sequences aligned with genomic sequence derived from the version 3 Bos taurus genome build (http://www.ncbi.nlm.nih.gov/spidey/) to identify intron-spanning probe target regions. Real-time PCR reactions (20 μl) consisted of the following: 10 μl TaqMan Gene Expression Master Mix (ABI, Singapore), 1 μl 20× Assay Mix [premix custom primers and TaqMan probe (ABI, Foster City, CA)], 5 μl of diluted cDNA (optimum concentration determined empirically for each assay, ranging from between 1/10 and 1/100), and 4 μl millipure H2O. Custom oligonucleotide and Taqman probe sequences are indicated in Table 4. Levels of 18S ribosomal RNA (rRNA) were determined using the Eukaryotic 18S rRNA TaqMan probe (ABI) and served to normalize the amount of input RNA and differences in cDNA synthesis efficiencies for relative real-time PCR calculations. The expression of 18S rRNA was found to be sufficiently stable for endogenous control of experimental data with standard deviations of <0.5 cycles and a data range spanning <2 cycles. Quantification of the endogenous control 18S assay and two other experimental assays were done on each 96-well plate, such that there were standard curves and eight experimental samples repeated in triplicate for each assay on the plate. PCR and fluorescence detection were carried out on an ABI7000 sequence detection system, using the following cycling parameters: initial incubation of 50°C for 2 min, followed by a polymerase activation step of 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Real-time PCR data were then assessed using the ABI PRISM 7000 SDS Software, and results were exported to Microsoft Excel for analysis. Statistical analysis of results was conducted using the relative standard curve method as outlined in the ABI technical bulletin “Guide to Performing Relative Quantitation of Gene Expression Using Real-Time Quantitative PCR” (http://www3.appliedbiosystems.com/cms/groups/mcb_support/documents/generaldocuments/cms_042380.pdf).
Expression Changes Revealed by Microarray Profiling
Volcano plot filtering (P < 0.05, fold change >1.5×) of microarray data was used to generate a gene list representing the most up- and downregulated genes across the 12 animals, as well as a list containing genes found to differ in expression between the large and small yield loss animals in the experiment. These lists of up- and downregulated genes consisted of 125 probes (115 discrete genes) and 51 probes (47 discrete genes), respectively, and are shown in Supplementary Tables S1 and S2.1
The two gene lists generated by volcano plot filtering (Supplementary Tables S1 and S2) were combined to generate a single list of sufficient size to allow molecular network generation and pathway analysis using IPA. Of these genes, 127 were deemed “network eligible” genes and were used to generate networks, pathways, and biological functions with relevance to the dataset. The top three networks built by IPA were: Network 1 - cell death, cancer, cellular movement (score = 51); Network 2 - cell death, cancer, gene expression (score = 26); and Network 3 - metabolic disease, lipid metabolism, molecular transport (score = 20). Network 1 had by far the largest “score” of these three networks and is illustrated in Fig. 2.
The top 10 biological functions identified by IPA analysis are shown in Fig. 3A. The genes (and associated P values and fold changes) specific to each of these biological functions are indicated in Table 1. P values for each of these biological functions were: cellular movement P = 2.20 × 10−9−7.23 × 10−3; cancer P = 2.07 × 10−8−7.23 × 10−3; cell death P = 2.07 × 10−8−7.23 × 10−3; tissue development P = 3.91 × 10−7−7.23 × 10−3; hematological disease P = 1.04 × 10−6−7.23 × 10−3; cell morphology P = 1.58 × 10−6−7.23 × 10−3; cellular growth and proliferation P = 2.77 × 10−6−7.23 × 10−3; cell-to-cell signaling and interaction P = 3.19 × 10−6−7.23 × 10−3; tumor morphology P = 5.91 × 10−6−7.23 × 10−3; cellular development P = 7.80 × 10−6−7.23 × 10−3.
The top five canonical pathways (and associated P values/ratios) generated by IPA are shown in Fig. 3B. Genes (and associated P values and fold changes) found within each of these canonical pathways are indicated in Table 1. P values for each of these canonical pathways were: PPARa/RXRa activation (P = 3.47 × 10−3, 5/145); p53 signaling (P = 3.63 × 10−3, 4/87); IGF-1 signaling (P = 3.78 × 10−3, 4/90); galactose metabolism (P = 4.26 × 10−3, 4/112); and Toll-like receptor (TLR) signaling (P = 6.62 × 10−3, 3/53).
Quantitative real-time PCR was used to assess the expression changes of 13 genes that showed differential expression by microarray analysis. Genes were chosen from the array dataset from each of the gene lists, such that a range of fold changes and signal intensities were examined. In all cases real-time PCR data corroborated the differential expression measured on the microarrays. This was true of relative intensity as well as the magnitude of change indicated by microarray analysis. When the fold changes of individual transcripts across the 12 arrays were plotted alongside the corresponding fold changes derived from the relative qRT-PCR results, significant (P < 0.05) Pearson product moment correlation coefficients were yielded (Table 2).
The primary objectives of this study were to gain insight into the molecular signaling networks driven by reduced milk removal frequency and to assess whether these networks reflect early apoptosis and milk synthesis inhibition pathways. The design of the experiment aimed to capture these changes by reducing milking frequency from twice per day to once per day for 5 days, giving insight into gene expression changes associated with a chronic reduction in milking frequency. Biopsies from once daily milked animals were taken from “full” glands to also include changes in gene expression that could reflect shorter-term responses to mammary engorgement, since, as discussed elsewhere, yield reduction in cattle begins 16–18 h postmilking (10), with the earliest morphological evidence of apoptosis occurring after 36–72 h of milk accumulation (47). While changes in milk yield during reduced milking frequency in ruminants are induced through local mechanisms (51), it is possible that some of the genes found to be differentially expressed in the current study vary due to the hormonal differences that result from reduced milking frequency and/or the improved energy balance associated with reduced milk yield (10).
Analysis of the dataset with IPA software supports the hypothesis that reducing milking frequency initiates apoptosis and regression processes in the mammary gland, with the top network assigned being “cell death, cancer, and cellular movement” (Fig. 2). The three functions identifying this network were also ranked as the top biological functions across the whole dataset (a ranking that encompasses many other molecules outside those reported in the functional network). This network, therefore, gives an indication of the possible direct interactions of these molecules (Fig. 2), with the biological functions of cell death, cancer, and cell movement describing broad-range functions relevant to this network, as well as to the entire dataset. Within each of these crude biological classifications, relevant subclassifications are also ranked, and it is noteworthy that within the functional class of cancer, the highest P values are paired to apoptosis and cell death subclassifications (P = 2.07 × 10−8 and P = 1.23 × 10−7, respectively). Other proliferative cancer subclassifications may also be relevant to the dataset, since recent evidence shows many similarities between advanced stage mammary involution and tumorigenic tissue remodeling (4, 45).
IPA Canonical Pathways
The top five canonical pathways identified by IPA were PPARa/RXRa activation, p53 signaling, IGF-1 signaling, galactose metabolism, and TLR signaling. Despite P values <0.05 being reported by Ingenuity for all five of these pathways, the size of the composite gene list entered into the program was relatively small, and so the numbers of genes represented in each of the canonical pathways is correspondingly small (∼3.5–5.5% of all molecules comprising the complete canonical pathway, Fig. 3B). Conservative interpretation as to the relevance of these pathways is therefore required. However, three of these pathways can be related to previously demonstrated facets of mammary apoptosis and involution.
The inference from IPA that galactose metabolism is altered in the study animals is consistent with one of the most obvious phenotypes of reduced milking frequency, a decrease in lactose synthesis associated with the reduction in milk volume (10, 18). The assignment of this canonical pathway derives from the appearance of several downregulated genes, namely B4GALT1, LALBA, and UGP2, all of which are key genes involved in the lactose synthesis pathway (28, 11, 41). Mitochondrial glycerol-3-phosphate acyltransferase gene (GPAM) and the major milk protein gene casein beta 2 (CSN2) are two other genes of particular relevance to milk synthesis (10, 53) that were also downregulated in the current study.
The IGF-1 canonical pathway also has obvious relevance to this study. The dataset contains three IGF binding proteins (IGFBP5, CTGF, and CYR61) that fall under the single node of IGFBPs in the canonical pathway. IGFBPs can modulate cell survival signaling elicited by IGFs and also have additional IGF-independent functions (see reviews in Refs. 2 and 24). Apoptosis and cell survival functions have been ascribed for all three of these molecules during mammary involution, and each of these is discussed further, in the section IGFBPs, Apoptotic Signaling, and Mechanical Stress Pathways below.
The p53 pathway has been shown to be necessary for normal mammary involution in mice (25). This too is consistent with the proposition that the glands of once daily milked animals share similarities with signaling pathways active in involuting mammary glands, particularly since the proapoptotic effects of p53 signaling are restricted to the earliest stage of mammary involution (25).
TLR signaling promotes apoptosis and activation of innate immunity (for reviews see Refs. 44, 55). However, the principal gene associating this IPA canonical pathway with our dataset is the TLR2, which was actually downregulated when the animals were milked less frequently. This is the reverse of what would be expected, since a decrease in TLR2 expression would likely reduce apoptotic and immune signaling during once daily milking. An explanation to reconcile this observation is that this change may simply be a reflection of the milk yield reduction phenotype, given that TLR2 has recently been reported to be a component of the milk fat globule membrane (42). The same may be true of assignment of the PPARa/RXRa activation pathway to our dataset since two of its key molecules CD36 and lipoprotein lipase (LPL), are also expressed in milk (42, 20). Due to these observations, further putative roles for these pathways have not been assigned.
Transcription Factor Activation and Apoptotic Signaling Networks
At the center of the top IPA-generated network (Fig. 2), and a gene found in four of the five canonical pathways assigned to the dataset, is the transcription factor JUN. JUN is strongly upregulated in nearly all of the once daily milked animals, with three separate probes independently placing this molecule in the upregulated transcripts list. The extent of JUN overexpression also appears to be related to milk yield reduction, as one probe indicated higher expression in the large yield loss animal group compared with the small yield loss group. JUN is a viral oncogene homolog (63) that can bind several other partner molecules to form the transcription factor complex AP1. These partner molecules include members of the FOS and CREB/ATF families, and notably, differently partnered AP1 complexes have already been demonstrated to be upregulated during involution of the murine mammary gland (35, 36) and during prostate involution (34). Also indicated at the center of the Ingenuity network is ATF3 (Fig. 2). This is the single most differentially expressed gene between the two yield loss groups in the study and is a potential AP1 complex binding partner with demonstrated affinity for JUN (23, 37). ATF3 has been associated with apoptosis regulation and involved with modulating the response to various different cell stressors including mechanical stretch and injury (6). It is unknown whether JUN and ATF3 are directly interacting in the mammary gland in the current study, but their presence indicates mammary cells from animals milked less frequently may be embarking on routes of apoptotic signaling that are also found during more advanced stage mammary involution.
The RELA gene forms part of another transcription complex, NF-κB, which is also associated with apoptotic and cell stress signaling (14). RELA is also upregulated in the study animals and is another key molecule occupying a central position in the highest scoring IPA-generated network (Fig. 2). NF-κB plays a central role in a myriad of inflammatory pathways (22) and has been reported to be upregulated in involuting mouse mammary tissue (8).
IGFBPs, Apoptotic Signaling, and Mechanical Stress Pathways
A member of the IGF-1 signaling canonical pathway, IGFBP5, was upregulated in response to reduced milk removal frequency in all of the animals in the current study and has been assigned developmental and apoptotic roles in a variety of tissues (3). IGFBP5 expression has been associated with apoptosis and early involution in several cell culture and rodent studies (33, 59, 60, 61), making its upregulation in once daily milked animals consistent with the premise that these cows display an early involution phenotype.
CYR61, a member of the CCN gene family, was the single most overexpressed transcript in the data set and is another member of the IGF-1 signaling canonical pathway. This ECM-associated protein is reported to be involved with the physiological processes of cell adhesion, angiogenesis, and tumorigenesis (17, 1, 26, 62). Through integrin-mediated adhesion, CYR61 is reportedly able to selectively modulate cell survival, with the molecule promoting apoptosis in fibroblasts (58) and transducing cell survival signals in other cell types (30, 31, 66). Another CCN family member, connective tissue growth factor (CTGF), which is also upregulated in the study animals, has been demonstrated to induce apoptosis in the breast cancer cell line MCF-7 (21). It is unknown whether CYR61 or CTGF promotes apoptosis of mammary epithelial cells within the context of the current study; however, their presence in the dataset may also point to the mechanisms underlying lactation inhibition more broadly. Both CCN molecules have been ascribed key roles in mechanically stressed tissues (7, 16), and CYR61 has been shown to be directly involved in transducing cell stretch signals through a number of molecular pathways, including mechanisms that sense actin structural dynamics (56). Thrombospondin (THBS1), the second most upregulated gene in our dataset, has also been demonstrated to be involved in mechano-signal transduction and participates specifically at the level of apoptosis induction (15). These findings indicate that lactation inhibition in less frequently milked animals, which culminates in the apoptotic induction of secretory epithelia, may occur as a result of signals tied directly to udder distension and proceed through pathways encompassing these three molecules.
Another finding that may reflect the mechanical stress exerted on mammary tissue as a result of reduced milk removal frequency is the observation that both claudin 4 (CLDN4) and claudin 8 (CLDN8) transcript levels were increased. Claudins are key tight junction proteins, and it has been proposed that mammary engorgement disrupts epithelial tight junction integrity in the mammary gland (10, 38). Increased tight junction permeability has been demonstrated as a result of less frequent milking in cows (52), and it is reasonable to suggest that disruption of these junctions could drive an increase in expression of specific claudin genes. In support of this hypothesis is the finding that several CLDN variants (CLDN1, CLDN3, and CLDN4) are upregulated on the first day postweaning in the mouse mammary gland (4).
Distention and mechanical stress in the engorged mammary gland have been proposed as a physiological trigger of mammary involution (39) and are concordant with the observation that milk secretion is regulated locally within the mammary gland, as opposed to systemically (51). The accumulation of feedback inhibitors in milk also fits this model, and recent work implicating serotonin as a molecule capable of regulating milk synthesis is noteworthy (54, 19). No obvious sign of serotonin-associated signaling was found in the current study except for the appearance of a single putative tryptophan transporter SLC6A14 (49). Tryptophan is the amino acid precursor of serotonin, and although the appearance of this single molecule is far from indicative of a serotonin-derived response to milking frequency, it is noteworthy that SLC6A14 is the second most differentially expressed molecule between the two yield loss groups in our dataset.
Overlap With Mammary Involution Gene Expression Analyses
Microarrays have been used by other investigators to examine the expression profiles of genes in the involuting tissues of mice (50, 9, 4) and the cow (48). Singh et al. (2008) (48) used cDNA microarrays in an involution time series experiment in cows and found SAT1, LALBA, and LTF differentially expressed in accordance with the changes observed in the current study. Stein et al. (50) conducted an involution time series in mice, and three of the 16 genes reportedly downregulated over this time span (CSN2, LPL, and CD36) were also downregulated in response to less frequent milk removal in the current study. Stein et al. also reported upregulation of JUN, ATF3, CTGF, SAT1, and CLDN4, corroborating changes in the current study that, as already discussed, are likely important to the physiological changes that occur as a result of udder distension. Blanchard et al. (4) examined gene expression changes that occur on the first day of mouse mammary involution, and in line with the current dataset as well as the Stein et al. (50) study, an increase in the expression of CTGF and CLDN4 was observed. Other molecules shared by these two datasets include TNFSF12A, GADD45A, GADD45B, MCL1, THBS1, KIT, and SLC34A2. With the exception of GADD45B, the expression direction of all of these molecules agreed with the expression changes exhibited in the animals in the current study.
In conclusion, a range of genes was found to be up- and downregulated in the mammary glands of cows shifted from twice daily to once daily milking, as well as genes whose expression correlated with the extent of milk yield reduction measured in the trial. These changes represent a reduction in milk protein synthesis as well as an induction of the apoptotic signaling networks characteristic of involution. This study demonstrates that apoptotic signaling is an early event in response to a change in milking frequency and, importantly, that initiation of these pathways may occur through sensing of mechanical stress present in the engorged mammary gland.
This research was supported in part by Dairy InSight (Wellington, New Zealand).
No conflicts of interest are declared by the authors.
We thank Gabrielle Bracefield for involvement in preparation of animal ethics applications and Brian Walsh for animal handling and management.
↵1 The online version of this article contains supplemental material.
- Copyright © 2010 the American Physiological Society