Transcriptome analysis of mammary epithelial cell gene expression reveals novel roles of the extracellular matrix.

BACKGROUND: The unique lactation strategy of the tammar wallaby (Macropus eugeni) has been invaluable in evaluating the role of lactogenic hormones and the extracellular matrix (ECM) in the local control of mammary gland function. However molecular pathways through which hormones and ECM exert their effect on wallaby mammary gland function remain unclear. This study undertakes transcriptome analysis of wallaby mammary epithelial cells (WallMEC) following treatment with mammary ECM from two distinct stages of lactation.
METHODS: WallMEC from MID lactation mammary glands were cultured on ECM from MID or LATE lactation and treated for 5 days with 1 μg/ml cortisol, 1 μg/ml insulin, 0.2 µg/ml prolactin, 650 pg/ml triodothyronine and 1 pg/ml estradiol to induce lactation. WallMEC RNA from triplicate ECM treatments was used to perform RNAseq.
RESULTS: ECM from MID and LATE lactation differentially regulated key genes in sugar and lipid metabolism. Seven pathways including galactose metabolism, lysosome, cell adhesion molecules (CAM), p53 signaling, the complement and coagulation and Nod-like receptor signaling pathways were only significantly responsive to ECM in the presence of hormones. The raw RNA-seq data for this project are available on the NCBI Gene Expression Omnibus (GEO) browser (accession number GSE81210).
CONCLUSIONS: A potential role of ECM in regulation of the caloric content of milk, among other functions including apoptosis, cell proliferation and differentiation has been identified. GENERAL SIGNIFICANCE: This study has used a non-eutherian lactation model to demonstrate the synergy between ECM and hormones in the local regulation of mammary function.

Background

The mammary extracellular matrix (ECM) as a regulator of mammary epithelial cell (MEC) function continues to attract the interest of researchers not only because of its role in breast homeostasis but also because of downstream implications in milk-driven neonatal development. In the latter context, understanding how MEC's are conditioned to secrete the complex mix of molecules that constitute milk would, in principal, enable the engineering of custom milks to suit specific developmental needs of individual babies. Studies have indeed demonstrated that the mix of milk proteins secreted by MEC's can be predetermined by varying the composition of ECM on which the MEC's are cultured [1], [2], [3]

However, the vast majority of studies linking the ECM to mammary epithelial function and milk composition have been conducted in eutherian models in which the milk composition remains largely constant throughout the lactation period [4]. The constancy of milk composition in these models makes it difficult to understand the relationship between the components of the ECM and their potential role in regulating individual milk protein genes. Non-eutherian models in which milk composition progressively changes during lactation have the potential to reveal novel molecular relationships associated with the ECM.

In marsupials, including tammar wallaby the lactating mother profoundly alters the composition of milk during the course of lactation [5]. Further, the progressive changes in milk composition in the individual mammary glands is marked by phase-specific induction and silencing of specific genes, providing what has been called natural gene knockout [6].

Tammar lactation is divided into four phases. The gestation phase (P1) lasts approximately 26 days culminating in the birth of a fetus-like young [7], [8], [9]. Phase 2 A (EARLY lactation) starts at parturition when the neonate, although extremely altricial, has motile capability and crawls from the birth canal to the pouch and attaches permanently to one of the four teats for approximately 100 days post-partum [7], [10], [11], [12], [13]. During EARLY lactation the mother produces relatively small volumes of dilute milk that is rich in carbohydrates but low in protein and lipid [14], [15]. Phase 2B (MID lactation) commences 100–120 days post-partum and continues for approximately 100 days during which the neonate remains in the pouch but relinquishes the teat and only re-attaches to suck [13], [14], [15]. Similar to EARLY lactation, the milk produced during MID lactation is dilute and high in carbohydrates but low in protein and lipids.

During Phase 3 (LATE lactation) the mammary gland enlarges significantly [16], producing large amounts of concentrated milk that is rich in protein and lipid but low in carbohydrates [5], [7], [14], [15]. LATE lactation also represents the period of most dramatic change in morphology and growth of the young including the switch from ectothermic to endothermic regulation of body temperature [5], [13], [17].

In an earlier study we demonstrated that treatment of MID lactation MEC's with ECM extracted from a late lactation mammary gland switches the MEC's to a late lactation phenotype [3]. In a separate study we showed that ECM regulates the temporally regulated expression of the MaeuCath1 gene, a cathelicidin that may play a role in MEC proliferation and antimicrobial protection of both the mother's mammary gland and the neonate [18]. The current study used the Illumina platform to perform RNA-seq of MID lactation MEC's cultured on MID and LATE lactation ECM in the presence or absence of lactogenic hormones. The intention was to explore the molecular pathways which are collectively impacted by changes in ECM resulting in altered milk protein gene expression during the tammar lactation cycle.

Methods

Experimentation with animals

Wallabies were kept in an open yard with adequate vegetation and water. Tissues were obtained following animal ethics approval by the Deakin University Animal Welfare Committee (AWC) in adherence to the “Australian code for the care and use of animals for scientific purposes (NHMRC).

Preparation of ECM

Prior to extraction of the ECM, the lactation phase of the mammary glands was determined by measuring the expression of tELP, tWAP and LLP-B (representing EARLY, MID and LATE lactation respectively) using quantitative PCR. The extraction of ECM from mammary tissue was performed as previously described [3], [19]. Briefly, snap-frozen mammary tissue (3 animals per group) was added to a ceramic mortar containing liquid nitrogen and pulverised to a fine powder followed by homogenisation in a high salt buffer (1 g of tissue in 2 ml) containing 2 mM n-ethylmaleimide (NEM), 200 μg/ml phenylmethylsulfonyl fluoride (PMSF) and 2x protease inhibitor (PI) cocktail (Sigma, P8340). The homogenate was enriched for insoluble ECM proteins by performing a second high salt wash and urea extraction, including a high speed centrifugation step (1 h at 100,000g) between the washes and the ECM serially dialysed to remove protease inhibitors and urea. The final dialysis was performed in Hams F12/ M199 media (Gibco®) in the absence of serum and the ECM aliquoted for freezing at −20 °C until required. ECM isolated using this procedure is enriched for ECM proteins by up to 60 fold and is depleted of soluble proteins by up to 90 fold but contains matrix-associated growth factors [20]

Isolation of tammar mammary epithelial cells (WallMEC)

Wallabies at day 136 of lactation representing MID lactation were euthanised and mammary tissue surgically removed and immediately placed in 1× Hanks’ Balanced Salt Solution (HBSS) (Sigma 55021C) containing 10 μg/ml penicillin/streptomycin (Gibco, USA) and 2.5 μg/ml Fungizone (Gibco). After removal of fat, the tissue was weighed and sliced repeatedly into tiny fragments followed by incubation at 37 °C for up to 4 h in 400 units/ml Collagenase Type 3 and 100 units/ml Hyaluronidase (Worthington) observing a ratio of 25 g tissue per 100 ml digest media and cells harvested by filtration, first through 200 µm mesh and then 53 µm mesh Nalgene filter. The filtrates were centrifuged at 80g for 5 min and cell pellets washed twice by suspending in HBSS and centrifuging at 80×g for 5 min. After the final wash, cells were resuspended in freezing media (90% FCS (Invitro Technologies)/10%DMSO) (Sigma, D8418), at a density 2 × 107 cells/ml and frozen in liquid nitrogen.

WallMEC culture and isolation of RNA

Culturing of WallMECs and isolation of total RNA was performed as earlier described by Wanyonyi et al. [3], [18]. Briefly, 200 µl of undiluted ECM was added to each well of a 6 well tissue culture plate (Costar; Corning Incorporated) and incubated at 37 oC for 1 h to allow coating of the wells. WallMECs at passage 4 or earlier were seeded at 5×104 cells per well in growth media (M199/Hams/HEPES media supplemented with 1 μg/ml cortisol (F), 10 ng/ml EGF, 1 μg/ml insulin (I), 1 mM glutamine, 20% horse serum, 5% foetal bovine serum and penicillin/streptomycin) into the ECM-coated wells. The cells were incubated at 37 oC/5% CO2 until they formed mature acini as earlier described [21], typically 13 days after seeding. Acini were considered mature after attaining a diameter of 200 µm or wider and appeared dark under bright field. Mature acini were washed once with 2 ml differentiation media (IF) consisting of 2% FBS in growth media without EGF and incubated for 5 days in IF containing lactogenic hormones (0.2 µg/ml prolactin, 650 pg/ml triodothyronine and 1 pg/ml estradiol). Control wells were treated with media containing IF. Representative images of acini were captured at 40x magnification using the bright field setting of the Leica TCS SP5 confocal microscope (Leica Microsystems) after estimating the average width of acini in five fields per replicate well (Fig. 2B-E). RNA was purified from each of the three replicate WallMEC cultures using the RNeasy Mini Kit (Qiagen). In order to ensure that there was no difference in the gene profiles of replicates, reverse transcription-PCR (RT-PCR) was performed on individual replicates using β-casein, tELP, tWAP, tLLP-B and GAPDH primers and the SsoFast EvaGreen method (BioRad). The final RNA libraries were generated by pooling the three replicates upon confirming that there was no difference in the level of expression of the lactation phase specific genes between them. The RNA libraries were designated MID_NH, MID_H, LATE_NH and LATE_H representing WallMECs cultured on MID ECM with no hormones, MID ECM with hormones, LATE ECM with no hormones and LATE ECM with hormones respectively. The sequencing service was contracted to Beijing Genomics Institute (BGI).

Fig. 2

: Acini morphology after culturing WallMEC on ECM extracted from MID and LATE ECM. A: Acini after culturing WallMEC on MID ECM for 13 days (before treatment with hormones). A represents mature acini. WM represents WallMEC cells in a monolayer. B: Acini cultured on MID ECM after 5 days of treatment with hormones. C: Acini after culturing WallMEC on LATE ECM for 13 days (before treatment with hormones). D: Acini cultured on LATE ECM after 5 days of treatment with hormones. Acini images were captured at 40x magnification in a bright field setting of the Leica TCS SP5 confocal microscope (Leica Microsystems).

RNA-seq data analysis

The sequence output in FastQ format [22] was screened for quality using FastQC software provided by Babraham Bioinformatics (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and the Macropus eugenii genome (Ensembl reference, Macropus_eugenii. Meug_1) used as a reference for annotation of the sequence reads. The pipeline for sequence data analysis consisted of the software described in the Tuxedo protocol [23]. Briefly, sequence reads were mapped to the M.eugenii genome using the Bowtie2 function of TopHat 2.0.5 [24], [25] and transcripts assembled using Cufflinks [26]. The final transcriptome assembly was performed by merging assembled transcripts using Cuffmerge and SAMtools [27] and the differential expression of the mapped reads determined using Cuffdiff. In order to visualise the differential expression of transcripts, volcano plots were generated using CummeRbund (R Package Package) and the location where the reads aligned confirmed using SeqMonk (version 0.22.0). In order to visualise the differential expression of genes between treatments, volcano plots were constructed comparing MID ECM versus LATE ECM with and without hormones. The plots are shown in Fig. 3A–D.

Fig. 3

Volcano plots comparing the fold change in gene expression of WallMEC cultured on MID and LATE ECM. A: MID ECM without hormones (NH) versus MID ECM in the presence of hormones (H). B: LATE ECM without hormones (NH) versus LATE ECM in the presence of hormones (H). C: MID ECM without hormones (NH) versus LATE ECM without hormones (NH). D: MID ECM in the presence of hormones (H) versus LATE ECM in the presence of hormones (H). The dots on the left side of zero represent down regulated genes while the dots on the right side represent up-regulated genes. Blue dots represent genes whose expression level is significantly different. The statistical significance of the difference between treatments expressed as –log10(p value) is shown on the y-axis. In order to enhance the level of confidence only genes with –log10(p value) higher than 2.5 were considered to be differentially expressed between the treatments.

Functional clustering and pathway analysis

In order to perform functional annotation clustering, the gene ID's for the differentially expressed genes (both up-regulated and down-regulated) were uploaded into DAVID v6.7 (the Database for Annotation, Visualization and Integrated Discovery) [28], [29] and the gene ID's converted to symbols which are compatible with the DAVID database. Only gene clusters with an enrichment value equal to or higher than 1.3, following authors’ recommendations [28] were considered representative of the treatments. Pathway analysis was performed using the BioCarta and KEGG databases which are available at DAVID (http://david.abcc.ncifcrf.gov/). A summary of the experimental protocol is presented in Fig. 1.

Fig. 1

Summary of the experimental protocol.

Statistical analysis

Statistical analysis of qPCR data from replicate samples was performed using one way ANOVA. In order to rank the abundance (enrichment) of genes, functional clusters and functional pathways, an enrichment score was calculated as the geometric mean of all the enrichment p-values (EASE scores) for each annotation term associated with the gene members in the cluster or pathway [28]. Since the geometric mean is a relative score and not an absolute p-value, the mean p-values were expressed as minus log. Higher enrichment scores are indicative of greater response of the gene members to the specific treatment and in the case of functional clusters, a cut-off of 1.3 was applied as recommended by Huang et al. [28]. The p-value (EASE score) was determined using a modified Fisher's exact test with a default cut-off of 0.1. Smaller p-values indicated greater significance of the enrichment of individual gene-terms.

Results and discussion

ECM – dependent changes in WallMEC acini morphology

Consistent with an earlier study, culturing of WallMEC on mammary ECM induced the formation of alveoli-like acini capable of expressing milk protein genes in response to lactogenic hormones [3]. Additionally, the size and morphology of WallMEC-derived acini was modifiable by culturing WallMEC on ECM from different phases of lactation, LATE ECM produced significantly smaller (210 µm in diameter versus 320 µm) but more numerous acini than MID ECM (Fig. 2) and modified the milk protein gene profile from a mid-lactation phenotype to a late-lactation phenotype. Since the expression profile of milk protein genes replicated previously published data [3] and was only used in this study to confirm the switch of WallMEC from a MID to LATE phenotype, this data has not been included.

Transcript abundance and differential regulation of gene expression by ECM and hormones

On average 6.6 million reads, 49 bp long were obtained for each of the four treatment groups, MID-NH, MID-H, LATE-NH and LATE-H, representing WallMEC treated with MID ECM without hormones, MID ECM with hormones, LATE ECM without hormones and LATE ECM with hormones respectively. Up to 78.8% of reads successfully mapped to transcripts in the wallaby reference genome (Supplementary data, Table 1). A summary of the counts of genes which were differentially regulated by ECM and hormones is provided in Supplementary data, Table 2. Splice indices based on the relative abundances of the splice variants did not reveal any difference between the treatment groups. Comparison of the amount of isoform switching between isoforms originating from common transcription start sites did not reveal any difference between the treatment groups. Isoform switching was determined based on the square root of the Jensen-Shannon divergence (sqrt(JS)) computed on the relative abundances of the splice variants (Cuffdiff)”.

Functional clusters of differentially expressed genes

The relative abundance of functional clusters for each dataset are shown in Supplementary data, Tables 3A and B. No clustering was done for the MID_NH versus MID_ H because only one gene was differentially expressed. Although functional annotation was attempted for the 5 differentially expressed genes in the LATE_NH versus LATE_H dataset, no functional cluster was obtained. In all the datasets, glycoproteins constituted the highest percentage of differentially expressed genes, averaging 16% (Supplementary data, Tables 3A and B).

Post-translational glycosylation of proteins regulates the function of a large number of proteins and in the context of milk, glycosylation protects proteins from proteolysis [30], modulates interaction of milk antimicrobials with target pathogens [31] and plays a role in neonatal development. Given the ubiquity of glycoproteins, it is difficult to ascertain from this study the role of ECM in their differential expression and further investigation is required.

Also prominent in all the datasets were the ECM and cell adhesion clusters (3%), dominated by matrix metalloproteases (MMP's) and ADAM metallopeptidase with thrombospondin (ADAMTS's) and they appeared to be regulated both by the ECM and hormones. For example while MMP1, MMP7, MMP9, MMP16, ADAMTS1, ADAMTS8 and ADAMTS15 were up-regulated in the MID_H versus LATE_H dataset, only MMP1, MMP16 and ADAMTS15 were up-regulated in the MID_NH versus LATE_NH dataset (higher in LATE_NH; Tables 4 and 5 in Supplementary data). The ECM related cluster also consisted of structural ECM proteins including collagen VII, a basal lamina protein secreted by epithelial cells and collagen XV which is expressed mainly in fibroblasts but also in renal, pancreatic, lung and placental epithelial cells [32]. Tenascin C which is expressed in fibroblasts but also in mammary epithelial cells [33] in a lactation phase-dependent manner [34], [35] was also differentially regulated by the ECM (Supplementary data, Tables 4 and 5).

Similar to the ECM related cluster, genes in the hormone response cluster were abundant (3%) in both NH and H datasets (Supplementary data, Tables 3A and B) but surprisingly absent in the down-regulated gene category of both clusters (Supplementary data, Tables 6 and 7). This lop-sided expression of hormone response genes requires further investigation because it can help dissect the hormone effect from the ECM effect.

Mammary luminal epithelial cells have been shown to express various MMP's [36], [37] and ADAMTS's [38] constitutively and in an inducible manner although these proteases are mainly expressed in the stroma and myoepithelium [38]. Importantly, the overarching cause for change in WallMEC phenotype was evidently the synergistic effect of hormones and the ECM since there were marked differences in gene expression between MID ECM and LATE ECM treated cells.

Differentially regulated functional pathways

Since only one gene (unknown) was differentially expressed between no hormone (NH) and hormone-treated (H) WallMEC cultured on MID ECM, pathway analysis was not performed for the MID-NH versus MID_H dataset. Similarly, although pathway analysis was performed for the 5 differentially expressed genes in the LATE_NH versus LATE_H dataset, the query did not return any pathway. The comparison between WallMEC cultured on MID ECM and LATE ECM in the presence of hormones (MID_H versus LATE_H) identified 12 pathways in which genes were differentially up-regulated (higher in LATE_H; Table 1). Nine out of these pathways were affected in the comparison between MID ECM and LATE ECM in the absence of hormones (MID_NH versus LATE_NH). Therefore although there appeared to be little or no influence of hormones in the change of WallMEC phenotype (Table 1A), the hormones exerted their effect in the presence of ECM. For example, while the entire p53 signaling pathway was not affected in the MID-NH versus LATE_NH dataset, five genes in the pathway were up-regulated in the MID-H versus LATE_H dataset (higher in LATE_H; Table 1B). Similarly, although the Galactose metabolism and Nod-like receptor signaling pathways were up-regulated in the MID_NH versus LATE_NH dataset (higher in LATE_NH), their enrichment score was significantly lower than the MID_H versus LATE_H dataset and below the cut-off of 1.3.

Table 1A

: Genes down-regulated by ECM grouped according to functional pathways.

MID_ NH versus LATE_NH	Regulated genes
No pathway	–
MID_H versus LATE_H	Regulated genes
Lysosome	cathepsins, ACP5, ASAH1, LAPTM
Complement and coagulation pathway	C1S, C2, C3,
Cell adhesion molecules	PVRL1, CNTN1, CDH3,

Table 1B

: Genes up-regulated by ECM grouped according to functional pathways.

MID_ NH versus LATE_NH	Regulated genes
Jak-STAT signaling pathway	IL6, CytokineR, PI3K, CycD
Cytokine-cytokine receptor interaction	IL6, IL11, LIF, BSF3, CSF3, GHR, CSF2RB, SF11B, SF21, IL1R2, HGF, FLT, IL18
Complement and coagulation cascades	F3, THBD, serpinE1
Hematopoietic cell lineage	G-SCF, IL−6, IL−11,CD121
TGF-beta signaling pathway	THBS1, BMP, FST, Smad1/5/8
Cell cycle	Bub1, BubR1, cd20, cyclinB, Plk1, cyclinD, MCM5
oocyte meiosis	IGF−1, Bub1, PI3K, cdc20, Plk1, CycB1, AC
Neurotrophin signaling pathway	TrkC, PI3K, IkB, PKCδ, RIMID, IRAK
Wnt signaling pathway	FRP, Wnt, Frizzled, Nkd, NFAT
MID_H versus LATE_H	Regulated genes
Jak-STAT signaling pathway	IL6, CytokineR, PI3K, CyclinD
Cytokine-cytokine receptor interaction	IL6, IL11, LIF, BSF3, CSF3, GHR, CSF2RB, SF11B, SF21, IL1R2
Complement and coagulation cascades	F3, THBD, serpinE1, PLAU, DF
Hematopoietic cell lineage	G-SCF, IL−6, IL−11,CD121
TGF-beta signaling pathway	Smad2,3, THBS1, BMP, FST, Smad1/5/8
p53 signaling pathway	CyclinD, CyclinB, IGF, PAI, p53R2
Cell cycle	Bub1, BubR1, Smad2,3, cd20, cyclinB, Plk1, cyclinD,
Oocyte meiosis	IGF−1, Bub1, PI3K, cdc20, Plk1, CycB1, AC
Neurotrophin signaling pathway	TrkC, PI3K, IkB
NOD-like receptor signaling pathway	RIMID, IkB, IL−6
Wnt signaling pathway	FRP, Wnt, Frizzled, Nkd, Smad3, NFAT, fra1, CycD, uterine
Galactose metabolism	Glucose phosphomutase, Glucose−1-phosphate, Hexokinase, uridylyltransferase UDP-galactose 4-epimerase

Equally interesting, three down-stream member genes of the canonical Wnt signaling pathway including Fra1, CyclinD and Uterine which link Wnt signaling to cell cycle and Smad3 which links it to TGF-β signaling were only up-regulated in the presence of hormones (Table 1B). The effect of hormones on WallMEC phenotype is further confirmed by the down-regulation of genes in three functional clusters including the lysosome pathway, complement and coagulation pathway and cell adhesion molecules (Table 1A) only in the presence of hormones. These observations confirm a key hypothesis in this study that the ECM sensitizes mammary epithelial cells to lactogenic hormones.

Galactose metabolism pathway

During LATE phase, the milk is rich in protein and lipid but low in carbohydrate and whereas the main carbohydrates in MID phase milk are complex oligosaccharides, LATE milk contains mainly monosaccharides [39]. Previous studies have attributed the decrease in oligosaccharides during LATE phase to the down-regulation of 3βGalT and 4βGalT and the increase in the expression of UDP-galactose hydrolase. 3βGalT and 4βGalT collectively catalyse the synthesis of oligosaccharides while UDP-galactose hydrolase breaks down UDP-galactose, a building block of oligosaccharides [40], [41].

Since changes in enzymes that catalyse galactose metabolism (β−1,3 galactosyltransferase and β−1,4 galactosyltransferase) are dependent on the phase of lactation, it has been suggested that these enzymes are regulated by the endocrine system and unknown factors within the mammary gland [39]. Prior to this study no association had been shown between ECM and the regulation of milk carbohydrate content during wallaby lactation.

This study reports for the first time the increase in the expression of the genes for four enzymes (EC 2.7.1.1, EC 2.7.7.9, EC 5.1.3.2 and EC 5.4.2.2) involved in galactose metabolism after culturing WallMEC on LATE ECM relative to MID ECM (Fig. 4A; Table 1B). UDP-glucose 6-dehydrogenase (EC 2.7.7.9) catalyses the reversible conversion of UDP-glucose into α-D-Glucose-1-P while glucose phosphomutase (EC 5.4.2.2) catalyses the conversion of α-D-Glucose-1-P to α-D-Glucose-6P both of which are downstream steps in the incorporation of D-galactose into the glycolysis pathway [42]. UDP-galactose 4-epimerase (EC 5.1.3.2) reversibly converts UDP-galactose to UDP-glucose [43] while hexokinase (EC 2.7.1.1) phosphorylates α-D-Glucose to form α-D-Glucose-6P [42]. This data suggests that although the reactions catalysed by these enzymes are reversible, changes in the ECM may progressively shift the equilibrium towards catabolism of carbohydrates leading to the marked decrease in carbohydrates during the LATE phase.

Fig. 4

A: Galactose metabolism pathway. Genes that were up-regulated after culturing WallMEC on LATE ECM relative to MID ECM are marked in blue and underlined. An enrichment score of greater than 1.3 (calculated as the geometric mean of all the enrichment p-values (EASE scores), was used to determine differentially expressed genes in the pathway. The p-value (EASE score) was determined using a modified Fisher's exact test with a default cut-off of 0.1. 4B:Lysosome pathway: Genes that were down-regulated following treatment of WallMEC with LATE_H relative to MID_H are marked in blue and underlined. Criteria for determining differentially expressed genes in the pathway are explained in the legend for Fig. 3A. 4C:p53 signaling pathway: Genes that were up-regulated following treatment of WallMEC with LATE_H relative to MID_H are marked in blue ovals. Criteria for determining differentially expressed genes in the pathway are explained in the legend for Fig. 3A. 4D:Wnt signaling pathway. Genes that were up-regulated following treatment of WallMEC with LATE_H relative to MID_H are marked in blue. Criteria for determining differentially expressed genes in the pathway are explained in the legend for Fig. 3A.

Lysosome pathway

The lysosome functional cluster was among the three differentially down-regulated pathways in the MID_H versus LATE_H dataset (lower in LATE_H), with four genes including cathepsin S, ACP 5 (Acid phosphatase 5), ASAH1 (N-acylsphingosine amidohydrolase) and LAPTM (Lysosomal-associated transmembrane protein 4B) being down-regulated (Table 1A; Fig. 4B). ACP5 catalyses the release of mannose 6-phosphate, the lysosome-targeting label, from proteins which are marked for proteolytic degradation in the lysosome [44] while ASAH1 (N-acylsphingosine amidohydrolase) an acid ceramidase digests ceramides [45].

The down-regulation of lysosomal proteases, manose 6-phosphatase and ceramidases upon culturing WallMEC on LATE ECM suggests decreased proteolysis and lipid catabolism which is consistent with high protein and lipid content in milk during LATE phase of lactation in the wallaby [5], [13], [46]. Whereas the ceramide content of wallaby milk is not well documented, in bovine milk, ceramides represent close to 30% of total phospholipid [47]. It is possible that during LATE lactation the ECM regulates milk lipid content by down-regulating lipid catabolic enzymes. Consequently, this study is the first to report a possible association between the ECM and the elevated milk lipid content during LATE lactation.

p53 signaling pathway

Four members of the p53 signaling pathway including Cyclin B, Cyclin D, IGF (insulin-like growth factor) and p53 R2 (the p53-regulated R2 subunit of ribonucleotide reductase) were up-regulated in LATE_H treated WallMEC relative to the MID_H (Table 1B; Fig. 4C). The relative up-regulation of cyclins B and D by treating cells with LATE-H (LATE ECM in the presence of hormones) suggests that LATE ECM may promote cell division [48] since both the cyclins stimulate the progression of cells through the cell cycle. This is consistent with the significant increase in the size of the mammary gland during phase 3 of wallaby lactation [16]. Similarly the increased expression of IGF implies that LATE ECM promotes the proliferative effect of growth hormone (GH) since GH exerts its growth stimulus through the induction of IGF expression.

However the up-regulation of p53R2 was somewhat discordant with the proliferative function of LATE ECM since p53R2 has been shown to arrest cell cycle by suppressing cyclin D1 expression [49]. In the context of wallaby lactation, the balance between the proliferative activity of cyclins and the cell-cycle-arresting activity of p53R2 may be needed to enable enhanced mammary growth during LATE lactation while allowing for gradual involution which is evident in late lactation [50].

Wnt signaling pathways

Seven genes in the canonical Wnt signaling pathway and one gene in the Wnt/Ca2+ signaling pathway (Table 1B; Fig. 4D), were up-regulated in the MID-H versus LATE_H dataset (higher in LATE_H). Four of these genes including Wnt, Frizzled, Nkd and NFAT were also up-regulated in the MID_NH versus LATE_NH dataset (higher in LATE_NH). Notably, three genes, Fra1, CycD and Uterine which link the Wnt signaling pathway to cell cycle progression [51] were only up-regulated in the hormone treated samples. Similarly, Smad3 which links Wnt signaling to TGF-β signaling pathway [51], [52] was expressed exclusively in the hormone treated samples.

Apart from Smad3, three genes in the TGF-β (Transforming growth factor beta) signaling pathway, including BMP, Smad8, THBS1 were up-regulated in the MID-NH versus LATE_NH and MID-H versus LATE_H datasets (higher in LATE). TGF-β negatively regulates epithelial cell proliferation [53] and has been known to induce apoptosis in several cells through the Smad pathway [54]. Although Smad8 is mainly regulated through the BMP pathway [9], studies in diverse cell types have shown that TGF-β signaling can be conveyed through activin-like kinase 1(ALK1) which induces phosphorylation of Smad8 [5], [6], [7], [8]. In endothelial cells, TGF-β signaling through ALK1 stimulates proliferation, consistent with increase in Smad8 expression in late lactation when cell proliferation is required for enlargement of the lactating gland. Studies have shown that TGF-β is down–regulated when MEC are cultured on ECM [55] and it represses β-casein expression in a Smad-dependent manner in MEC [56] implying that TGF-β may inhibit MEC differentiation. Given the pleiotropic nature of the wnt and TGF-β signaling pathways, these changes in gene expression may highlight the central role of ECM in cell signaling rather than specific regulation of mammary epithelial cell function.

JAK/STAT signaling pathway

The response of the mammary gland to prolactin, leading to milk protein gene expression and lactation is mediated by the JAK/STAT signaling pathway through Jak2 and Stat5 [57]. In this study four members of the JAK/STAT signaling cascade including IL-6, CytokineR PIK3 and CycD were up-regulated in both the MID_NH versus LATE_NH and MID_H versus LATE_H datasets (higher in LATE) (Table 1B). Therefore it appears that the differential expression of these signaling molecules was not hormone-dependent and was mediated by the ECM. This is consistent with earlier studies in the wallaby which showed that although there is increased prolactin expression during LATE lactation [13], the change in milk protein gene expression pattern between MID and LATE lactation was not dependent on hormones [5].

The complement cascade

The complement system is a multi-tasking family of more than 30 proteins that function in a cascade to effect immunological reaction to pathogens by opsonising pathogens, activation of immune cells, chemotaxis, and clearance of immune cells among other functions [58], [59]. In this study two of the five up-regulated genes (PLAU and DF) in the complement and coagulation system were up-regulated only in the MID_H versus LATE_H dataset (Table 1B). PLAU (Urokinase-type plasminogen activator) is a serine protease that plays a role in the breakdown of the ECM by cleaving plasminogen to form plasmin. PLAU therefore may provide a link between a possible hormone-ECM mediated remodelling of the extracellular matrix during LATE lactation in the wallaby and is worth further investigation. DF (Complement factor D) plays a role in the alternative complement pathway by cleaving Bf [60]

Of particular interest is the down-regulation of the three serially connected members (C1s, C2 and C3) of the classical complement pathway in the MID_H versus LATE_H dataset (Table 1A). C1s, C2 and C3 are part of a multi-step process involving the mu chain of immunoglobulin M (IgM) that generates C3b, a well known opsonin for macrophages and neutrophils [61], [62], [63]. Whether this down-regulation of the complement molecules is aimed at preventing autophagy of immune cells during LATE lactation is not clear from this study. However it is apparent that the ECM and hormones act synergistically to regulate immune activity during LATE lactation through the complement cascade.

Conclusion

This study has used RNAseq to suggest novel roles of ECM in the control of mammary epithelial function. The ground-breaking association between ECM and the expression of key enzymes in carbohydrate and lipid metabolism provides a possible mechanism for the control of the calorie content of milk. Further, the study has demonstrated the differential regulation of p53, TGF-β, Wnt and JAK-STAT signaling and the complement and coagulation cascade by lactation phase-specific ECM.

Author contributions

SSW, KRN and CL conceived the study and formulated its design. SSW performed experiments including isolation of mammary ECM and generation of WallMEC primary cultures and RNA libraries. AK and CL performed read quality control, annotation and transcriptome assembly while RDP performed final pathway analysis. SSW wrote the first draft of the manuscript and co-ordinated input by all co-authors. All the authors have reviewed and approved the submission of this manuscript.

Conflict of interest

The authors individually and collectively declare that they have no conflict of interest.

Funding

This work was supported by a Ph.D. scholarship awarded to SSW by Deakin University, Australia.