Co-Expression Network and Pathway Analyses Reveal Important Modules of miRNAs Regulating Milk Yield and Component Traits

Co-expression network analyses provide insights into the molecular interactions underlying complex traits and diseases. In this study, co-expression network analysis was performed to detect expression patterns (modules or clusters) of microRNAs (miRNAs) during lactation, and to identify miRNA regulatory mechanisms for milk yield and component traits (fat, protein, somatic cell count (SCC), lactose, and milk urea nitrogen (MUN)) via miRNA target gene enrichment analysis. miRNA expression (713 miRNAs), and milk yield and components (Fat%, Protein%, lactose, SCC, MUN) data of nine cows at each of six different time points (day 30 (D30), D70, D130, D170, D230 and D290) of an entire lactation curve were used. Four modules or clusters (GREEN, BLUE, RED and TURQUOISE) of miRNAs were identified as important for milk yield and component traits. The GREEN and BLUE modules were significantly correlated (|r| > 0.5) with milk yield and lactose, respectively. The RED and TURQUOISE modules were significantly correlated (|r| > 0.5) with both SCC and lactose. In the GREEN module, three abundantly expressed miRNAs (miR-148a, miR-186 and miR-200a) were most significantly correlated to milk yield, and are probably the most important miRNAs for this trait. DDR1 and DDHX1 are hub genes for miRNA regulatory networks controlling milk yield, while HHEX is an important transcription regulator for these networks. miR-18a, miR-221/222 cluster, and transcription factors HOXA7, and NOTCH 3 and 4, are important for the regulation of lactose. miR-142, miR-146a, and miR-EIA17-14144 (a novel miRNA), and transcription factors in the SMAD family and MYB, are important for the regulation of SCC. Important signaling pathways enriched for target genes of miRNAs of significant modules, included protein kinase A and PTEN signaling for milk yield, eNOS and Noth signaling for lactose, and TGF β, HIPPO, Wnt/β-catenin and cell cycle signaling for SCC. Relevant enriched gene ontology (GO)-terms related to milk and mammary gland traits included cell differentiation, G-protein coupled receptor activity, and intracellular signaling transduction. Overall, this study uncovered regulatory networks in which miRNAs interacted with each other to regulate lactation traits.


Introduction
MicroRNAs (miRNAs) are small noncoding RNA molecules about 22 nucleotides long which regulate gene expression post-transcriptionally, and play key roles in a wide range of biological processes. Many lines of evidence indicate that several miRNAs can work together to affect target genes in the same or different biological pathways [1][2][3]. Complex relationships exist between miRNAs, since they (1) might be clustered together both by sequence similarity and genomic location, (2) might be clustered into the same miRNA family, (3) may regulate or are regulated by the same transcription factor, and (4) might share target genes in certain biological processes [3,4]. Indeed, several approaches have been proposed to explore these relationships, such as miRNA-miRNA synergistic network [5] and co-expression analyses [6,7]. The first approach is based on the downstream study of miRNA target genes through construction of networks based on different weighted methods on common target genes, as well as the interactions among them [3,5]. Meanwhile, weighted co-expression network is based on construction of interaction networks (modules) of miRNAs with similar expression patterns, whereby miRNAs in the same module interact with one another to regulate the same or similar biological processes [3,8]. Moreover, every module is assigned an eigenvalue, which enables determination of the relationship (correlation) between modules and traits of interest. The hub genes of each module points to the most active miRNAs in each network, which are potentially the most important miRNAs regulating the transcriptomic mechanisms underlying the traits.
Lactation is a complex process known to be controlled by various regulatory mechanisms, including miRNAs [9,10]. The roles of miRNAs in dairy lactation or mammary gland development have been investigated by several authors [11][12][13][14][15][16][17][18]. Recently, we reported that 58 miRNAs were dynamically and differentially expressed across lactation stages, and that 19 miRNAs were differentially expressed throughout lactation in a significant and time-dependent manner [19]. These results suggest that miRNAs might interact with each other to regulate gene expression throughout lactation. However, it was not clear if these miRNAs could interact with each other to regulate lactation phenotypes, and what mechanisms underlie possible interactions. This study therefore aimed to (i) characterize groups (modules) of high interacting miRNAs during an entire bovine lactation curve using weighted gene co-expression network analysis (WGCNA), (ii) correlate important modules with milk yield and component traits, and (iii) enrich target genes of miRNAs from important modules, for exploration of the signaling pathways and upstream transcriptional regulators of milk yield and component traits.

Milk Component Yield Trend during a Lactation Curve
Quarter and interaction between quarter and day had no significant effects on milk yield and milk components, while day significantly affected (p < 0.05) milk yield and milk components (except Fat%, p = 0.132) ( Table 1). Milk production was similar from D30 up to D170, and decreased significantly by D230 and D290 (p < 0.001). Except for D70, Protein% on other days numerically increased slightly with significantly (p < 0.05) higher values by D230 (3.22 ± 0.29%) and D290 (3.57 ± 0.40%) as compared to D30 (2.98 ± 0.29%). Similarly, somatic cell count (SCC) increased numerically being significant (p < 0.05) on D230 and D290 as compared to D30. On the contrary, Lactose% decreased significantly (p < 0.01) on D230 (3.16 ± 1.21%) and D290 (3.33 ± 0.93%) as compared to D30 (4.36 ± 0.17%). There was no clear trend for milk urea nitrogen (MUN) content (mg/dL), since it had highest numerical increase on D170 (p < 0.05) but remained similar across all other time points throughout lactation.

Important miRNA Modules for Milk and Component Yields
Using the WGCNA approach, we identified eight different modules of co-expressed miRNAs, which were assigned different colors ( Figure 1). miRNA membership in the modules ranged from 32 (BLACK module) to 164 (TURQUOISE module) ( Figure 1). Results of the correlations between module eigengene values with traits are shown in Figure 1. Four modules (BLUE, GREEN, RED and TURQUOISE) were each significantly correlated with at least one trait (|r| > 0.5). A positive correlation was found between the GREEN module and milk yield (r = 0.57 and p = 2 × 10 −7 ). Two modules, RED and TURQUOISE, were significantly negatively (r = −0.57, p = 2 × 10 −7 and r = −0.57, p = 3 × 10 −7 ) and positively (r = 0.58, p = 2 × 10 −7 , r = 0.54, p = 1 × 10 −6 ) correlated with lactose and SCC, respectively ( Figure 1). The BLUE module was significantly negatively (r = −0.53, p = 2 × 10 −6 ) correlated with lactose ( Figure 1). No module was significantly correlated with Fat%, Protein% and MUN. Shared target genes of members (11 miRNAs) of the GREEN module are shown in Figure 2. Numbers of miRNAs with module membership values >0.8 were 34, 11, 39 and 15 for the BLUE, GREEN, TURQUOISE and RED modules, respectively. Details of significant module memberships and correlations of modules with traits are shown in Table 2 (GREEN module), Table 3 (BLUE module), Table 4 (RED module) and

Important miRNA Modules for Milk and Component Yields
Using the WGCNA approach, we identified eight different modules of co-expressed miRNAs, which were assigned different colors ( Figure 1). miRNA membership in the modules ranged from 32 (BLACK module) to 164 (TURQUOISE module) ( Figure 1). Results of the correlations between module eigengene values with traits are shown in Figure 1. Four modules (BLUE, GREEN, RED and TURQUOISE) were each significantly correlated with at least one trait (|r| > 0.5). A positive correlation was found between the GREEN module and milk yield (r = 0.57 and p = 2 × 10 −7 ). Two modules, RED and TURQUOISE, were significantly negatively (r = −0.57, p = 2 × 10 −7 and r = −0.57, p = 3 × 10 −7 ) and positively (r = 0.58, p = 2 × 10 −7 , r = 0.54, p = 1 × 10 −6 ) correlated with lactose and SCC, respectively ( Figure 1). The BLUE module was significantly negatively (r = −0.53, p = 2 × 10 −6 ) correlated with lactose ( Figure 1). No module was significantly correlated with Fat%, Protein% and MUN. Shared target genes of members (11 miRNAs) of the GREEN module are shown in Figure 2. Numbers of miRNAs with module membership values >0.8 were 34, 11, 39 and 15 for the BLUE, GREEN, TURQUOISE and RED modules, respectively. Details of significant module memberships and correlations of modules with traits are shown in Tables 2 (GREEN module), 3 (BLUE module), 4 (RED module) and 5 (TURQUOISE module).

Target Genes of miRNA Members in BLUE, GREEN, TURQUOISE and RED Modules
The miRNA members of the BLUE, GREEN, TURQUOISE and RED modules were predicted to target 3361, 1232, 4241 and 979 unique genes, respectively (Table S1a-d). Many miRNAs shared the same target genes, especially when they had the same seed region (Tables 2-5). The common target genes shared among members of the GREEN module are shown in Figure 2, and of the other three modules in Figure 3a (BLUE module), Figure 3b (RED module) and Figure 3c (TURQUOISE module). DDHD1 and DR1 genes were the most common targets for the GREEN module members, since they were each targeted by six miRNAs (Figure 2). MIER1, RNPEP and YME1L1 genes were also targeted by five different miRNAs in the GREEN module. Meanwhile, ENTPD3, RSBN1 and NAA15 genes were the most common targets in the BLUE (targeted by six miRNAs), TURQUOISE (targeted by seven miRNAs) and RED (targeted by six miRNAs) modules, respectively (Figure 3a-c and Table S1).

Target Genes of miRNA Members in BLUE, GREEN, TURQUOISE and RED Modules
The miRNA members of the BLUE, GREEN, TURQUOISE and RED modules were predicted to target 3361, 1232, 4241 and 979 unique genes, respectively (Table S1a-d). Many miRNAs shared the same target genes, especially when they had the same seed region (Tables 2-5). The common target genes shared among members of the GREEN module are shown in Figure 2, and of the other three modules in Figure 3a (BLUE module), 3b (RED module) and 3c (TURQUOISE module). DDHD1 and DR1 genes were the most common targets for the GREEN module members, since they were each targeted by six miRNAs (Figure 2). MIER1, RNPEP and YME1L1 genes were also targeted by five different miRNAs in the GREEN module. Meanwhile, ENTPD3, RSBN1 and NAA15 genes were the most common targets in the BLUE (targeted by six miRNAs), TURQUOISE (targeted by seven miRNAs) and RED (targeted by six miRNAs) modules, respectively (Figure 3a-c and Table S1). Predicted target genes of miRNAs in the GREEN module, and important hub genes and transcription factors. The red square nodes are miRNA members of the GREEN module, the blue round nodes are common predicted target genes for these miRNAs, the yellow round nodes are most highly predicted common (hub) target genes for these miRNAs, and the black V shape is the transcription regulator targeted by miRNAs, which also targets other predicted target genes in the networks. Predicted target genes of miRNAs in the GREEN module, and important hub genes and transcription factors. The red square nodes are miRNA members of the GREEN module, the blue round nodes are common predicted target genes for these miRNAs, the yellow round nodes are most highly predicted common (hub) target genes for these miRNAs, and the black V shape is the transcription regulator targeted by miRNAs, which also targets other predicted target genes in the networks.     The yellow V-like shape represents miRNAs, red round shape represents target genes, and the hub gene (most common target) is in the center, and is represented by a hexagon shape (yellow).

Enriched Gene Ontologies for Target Genes of miRNA Members of the BLUE, GREEN, TURQUOISE and RED Modules
The number of gene ontology (GO)-terms enriched for the BLUE, GREEN, TURQUOISE and RED modules were 158, 141, 174, and 65, respectively (Table S2a-d). Among them, 11 GO-terms were common to all the groups ( Figure S1). The scatter plots of enriched GO-terms for the GREEN, BLUE, RED and TURQUOISE modules are shown in Figures 4 and S2-S4, respectively. In the GREEN module, cell differentiation (p = 3.2 × 10 −4 ), intracellular components (p = 9.6 × 10 −10 ), and G-protein coupled receptor activity (p = 3.7 × 10 −5 ) were the most significantly enriched biological process, The yellow V-like shape represents miRNAs, red round shape represents target genes, and the hub gene (most common target) is in the center, and is represented by a hexagon shape (yellow).

Enriched Gene Ontologies for Target Genes of miRNA Members of the BLUE, GREEN, TURQUOISE and RED Modules
The number of gene ontology (GO)-terms enriched for the BLUE, GREEN, TURQUOISE and RED modules were 158, 141, 174, and 65, respectively (Table S2a-d). Among them, 11 GO-terms were common to all the groups ( Figure S1). The scatter plots of enriched GO-terms for the GREEN, BLUE, RED and TURQUOISE modules are shown in Figure 4 and Figures S2-S4, respectively. In the GREEN module, cell differentiation (p = 3.2 × 10 −4 ), intracellular components (p = 9.6 × 10 −10 ), and G-protein coupled receptor activity (p = 3.7 × 10 −5 ) were the most significantly enriched biological process, cellular component, and molecular function GO-terms, respectively ( Figure 4, Table S2b). The BLUE and TURQUOISE modules shared three most significantly enriched GO-terms, which were sensory perception of chemical stimulus, intracellular part and G-protein coupled receptor activity for biological process, cellular component and molecular function, respectively ( Figures S2 and S4, and Table S2a,c). For the RED module, cellular component morphogenesis (p = 3 × 10 −3 ), cytoplasm (p = 1.1 × 10 −5 ) and N-acyltransferase activity (p = 3.7 × 10 −2 ) were the most significantly enriched biological process, and cellular component and molecular function GO-terms, respectively ( Figure S3 and Table S2d).  Table S2b). The BLUE and TURQUOISE modules shared three most significantly enriched GO-terms, which were sensory perception of chemical stimulus, intracellular part and G-protein coupled receptor activity for biological process, cellular component and molecular function, respectively ( Figures S2 and S4, and Table S2a,c). For the RED module, cellular component morphogenesis (p = 3 × 10 −3 ), cytoplasm (p = 1.1 × 10 −5 ) and N-acyltransferase activity (p = 3.7 × 10 −2 ) were the most significantly enriched biological process, and cellular component and molecular function GO-terms, respectively ( Figure S3 and Table  S2d).

Signaling Pathways and Transcription Factors Enriched for miRNA Members of the BLUE, GREEN, TURQUOISE and RED Modules
A total of 18, 15, 21 and 11 canonical signaling pathways were significantly enriched for the target genes predicted for 34, 11, 39 and 15 miRNA members of the BLUE, GREEN, TURQUOISE and RED modules, respectively, using IPA core analysis ( Figure 5 and Tables S3a-d). Rac, PTEN and protein kinase A signaling pathways were the most significantly enriched for target genes of miRNA members of the GREEN module (Table S3a), and consequently are the most significant pathways for milk yield. Other enriched relevant pathways in the module for milk yield included 3-phosphoinositide biosynthesis, RAR activation, and signaling pathways of PTEN, ErbB, DNA methylation and transcriptional repression, eNOS, and Neuregulin ( Figure 5, Table S3a). Meanwhile, eNOS and Endothelin-1 signaling were the most significantly enriched for the BLUE module (Table S3b). Other relevant pathways for the BLUE module included NF-κB signaling, Notch signaling, phototransduction, and TR/RXR activation pathways. In the TURQUOISE module, several pathways related to cell cycle checkpoint were enriched (Table S3c); meanwhile, several pathways involved in nucleic acid metabolism were significantly enriched for the RED module (Table S3d).
A total of 35, 13, 99 and 9 transcription factors were enriched for target genes of miRNA members of the BLUE, GREEN, TURQUOISE and RED modules, respectively (Tables 6 and S4). The transcription factor HHEX was targeted by miR-148a. Moreover, HHEX also regulates genes (VEGFA, NRP1 and

Signaling Pathways and Transcription Factors Enriched for miRNA Members of the BLUE, GREEN, TURQUOISE and RED Modules
A total of 18, 15, 21 and 11 canonical signaling pathways were significantly enriched for the target genes predicted for 34, 11, 39 and 15 miRNA members of the BLUE, GREEN, TURQUOISE and RED modules, respectively, using IPA core analysis ( Figure 5 and Table S3a-d). Rac, PTEN and protein kinase A signaling pathways were the most significantly enriched for target genes of miRNA members of the GREEN module (Table S3a), and consequently are the most significant pathways for milk yield. Other enriched relevant pathways in the module for milk yield included 3-phosphoinositide biosynthesis, RAR activation, and signaling pathways of PTEN, ErbB, DNA methylation and transcriptional repression, eNOS, and Neuregulin ( Figure 5, Table S3a). Meanwhile, eNOS and Endothelin-1 signaling were the most significantly enriched for the BLUE module (Table S3b). Other relevant pathways for the BLUE module included NF-κB signaling, Notch signaling, phototransduction, and TR/RXR activation pathways. In the TURQUOISE module, several pathways related to cell cycle checkpoint were enriched (Table S3c); meanwhile, several pathways involved in nucleic acid metabolism were significantly enriched for the RED module (Table S3d).
A total of 35, 13, 99 and 9 transcription factors were enriched for target genes of miRNA members of the BLUE, GREEN, TURQUOISE and RED modules, respectively (Table 6 and Table S4).

Milk Yield and Components during Lactation
The characteristics of milk yield during a bovine lactation curve are well documented. Milk yield increases from the first day of lactation, to peak milk (around day 60), followed by a gradual decrease until the end of lactation [20,21]. Previously, we have shown that milk yield of cows in this study followed a similar pattern expected of a standard lactation curve [19], such as in Wood's model [22]. Significant variations by day recorded for milk components (except for Fat%) is supported by previous studies, which reported significant variations for milk component traits by lactation day, and also by periods of lactation [23][24][25][26]. Nevertheless, the number of phenotypic records in our data was small, which might explain the non-significant effect of day observed for milk Fat% in this study.

miRNAs, Hub Target Genes, Gene Ontologies, Pathways and Transcription Factors Involved in Milk Yield
A notable result in this study was the significant positive correlation (r = 0.57, p = 2 ×10 −7 ) between the GREEN module and milk yield. Some miRNA members of the GREEN module are known to have potential roles in various aspects related to milk yield, including mammary gland development, such as miR-200a [27,28], immune response and epigenetic upregulation, such as miR-148a [29][30][31], or nutrient response, such as miR-141 [18]. In fact, three miRNAs (miR-148a, miR-186 and miR-200a) that displayed the highest positive correlations with milk yield (r = 0.52-0.54) ( Table 2) are among the most abundantly expressed miRNAs in mammary gland tissues or milk [18,19,32,33], suggesting important roles for these miRNAs during lactation [34]. For instance, both miR-148a and miR-200a are known to be involved in the regulation of cell proliferation and death [28,[35][36][37], crucial processes for the regulation of milk yield. Recently, Melnik et al. [30] reviewed the epigenetic regulatory roles of miR-148a in lactation, and suggested that bovine miR-148a could be a critical factor for human health, since it is a component of milk fat globule and milk exosomes and highly expressed in milk, as well as having an identical seed region with human miR-148a. Moreover, 8 of the 11 miRNAs with module membership values >0.8 in the GREEN module were predicted to target two hub genes (DR1 and DDHD1) (Figure 2). Interestingly, DR1 (encodes a TATA box-binding protein associated phosphoprotein that represses RNA polymerase II [38]) was also predicted to be targeted by miR-2285e, miR-2285c, miR-141, and miR-200a. Although no direct association has been reported between DR1 and milk yield up to now, its role in gene transcription regulation suggest that DR1 is an important hub gene for miRNA network regulation of milk yield. DDHD1 gene, a member of the intracellular phospholipase family, encodes phospholipase enzyme, and is involved in catalyzing the degradation of phosphatidic acid, as well as attenuating cell activation [39]. The DDHD1 gene might play roles in milk yield via its roles in various biological processes, such as phospholipid metabolism and related signaling pathways [40]. Further supporting evidence for the role of GREEN module miRNAs in milk yield came from the enrichment results of their target genes. Gene ontology enrichment showed the importance of cell differentiation (p = 3.2 × 10 −4 ) in milk yield, since it was the most significantly enriched biological process GO-term ( Figure 4). Moreover, intracellular signaling transduction and embryo development were also among the most significant enriched biological processes. This supports the relationship between reproduction and milk yield in cows [41][42][43], since most lactating dairy cows are also gestating. Similarly, signaling pathways important for mammary gland development, as well as for lactation, such as GO-term intracellular signaling transduction, were highly enriched for milk yield [44]. Several known and notable signaling pathways for milk yield were enriched including PTEN, prolactin, Rac, protein kinase A, neuregulin, ErbB, and TGF-β signaling ( Figure 5). The roles of these pathways have been discussed in our previous study on the same animals [19]. For instance, prolactin and PTEN signaling are important pathways for hormone regulation of milk secretion [44,45], while TGF-β might have roles in mammary gland health [46,47] which is important for the maintenance of milk yield.
It is well documented that miRNAs and transcription factors can co-regulate gene expression [48][49][50]. However, little is known about the interaction of transcription factors and miRNAs in lactation regulation in general, and in milk yield regulation in particular. Enrichment of transcription factors for target genes of miRNAs in the GREEN module identified significantly enriched transcription factors, such as EHMT2, ZNF350, HHEX, MITF, CCND1, TP53, SP1, RYBP, TAL1, and members of SMAD family, including SMAD 2, 3, and 7. Interestingly, miR-148a regulates HHEX, while HHEX regulates several genes, including VEGEA (targeted by miR-186), NRP1 (targeted by miR-148a), or MYH10 (targeted by miR-141 and miR-200a) in the GREEN module ( Figure 2). HHEX plays an important role in the development of multiple organs, including liver, thyroid, and forebrain, via interaction with other signaling molecules [51]. Using cell culture experiments, Puppin et al. [52] observed that modification of HHEX protein localization occurs during lactation and tumorigenesis, and further suggested that HHEX may play a role in differentiation of mammary epithelial cells. However, there has been no functional study on the role of HHEX in lactation and milk production. Other transcription factors like CCND1, RYBP, and TP53 have been reported to play important roles in cellular cycle regulation and related processes [53][54][55], so they can either function as mediators or regulators of miRNA networks controlling milk yield. Furthermore, the regulation of milk fat synthesis might be an interplay between transcription factors and miRNAs (as well as other non-coding RNAs), and since these relationships are temporal and spatial, each type of molecule may have dominant roles at certain time points. Further studies, on the expression profile of transcription factors and milk at different time points throughout a lactation curve, may provide more clues about these relationships.

miRNAs, Hub Genes, Gene Ontologies, Pathways, and Transcription Factors Regulating Milk Components
Fat content did not change significantly during the course of lactation (Table 1), and no miRNA module was significantly correlated with Fat% and MUN. Meanwhile, lactose was significantly correlated with three different modules (RED, BLUE, and TURQUOISE). miR-18a in the BLUE module was the most significantly correlated to lactose (r = −0.65) (Table 3). miR-18a is a member of miR-17/92 cluster, and plays a role in tumor progression [56,57]. For example, miR-18a impairs cancer cell growth via inhibition of the expression of CDC42 gene [56], or negatively regulates the expression of PIAS3, an inhibitor of STAT3 [58]. Wu et al. [58] reported that the expression of miR-18a was positively correlated with the expression of STAT3 gene in gastric adenocarcinoma tissues. These authors further indicated that the correlated expression was via a down-regulation of PIAS3 gene by miR-18a. Moreover, Cai et al. [59] reported that promoter binding of STAT3 is required for transactivation of miR-148a (and other members of miR-17/92 cluster) in human macrophage cells following Toxoplasma infection. According to Brock et al. [60], a promoter region of miR-18a and other members of miR-17/92 contain a functional binding site for STAT3, which transcriptionally activates these miRNAs. However, it is not clear how miR-18a regulates milk lactose metabolism, but we speculate that its regulatory role is via targeting important molecules involved in lactation, such as STAT3.
Some notable members of the BLUE module, such as members of miR-221/222 cluster, are involved in breast cancer cell regulation by targeting multiple pathways [61], such as those promoting epithelial-to-mesenchymal transition [62]. miR-221 also targets PTEN [63], an important transcription factor in lactation regulation [44]. PTEN is known to regulate mammary cell growth, proliferation, and survival, by down-regulating important pathways of lactation, such as PI3K-AKT and mitogen-activated protein kinase pathways [64]. For instance, PTEN can down-regulate PTEN-AKT pathway, which is required for the initiation of lactation through the induction of autocrine prolactin [65]. Overexpression of PTEN was shown to down-regulate the expression of MAPK, CCND1, AKT, MTOR, S6K1, STAT5, SREBP1, PPARγ, PRLR, and GLUT1 genes, or up-regulate EIF4EBP1 gene, which are important for lactation related signaling pathways [44]. Moreover, Wang et al. [32] observed an increase in the expression of miR-221 during lactation, and suggested a role in the control of endothelial cell proliferation or angiogenesis. The most significantly enriched biological process GO-term for target genes of BLUE module miRNAs, was "sensory perception of chemical stimulus", required for an organism to receive a sensory chemical stimulus, convert it to a molecular signal, and recognize and characterize the signal [66]. Interestingly, the most significant molecular function GO-term was "G-protein coupled receptor activity", which involves a combination of extracellular signaling and signal transduction across membranes by activation of associated G-protein [66]. Wickramasinghe et al. [67] reported that "G-protein coupled receptor activity" was significantly enriched for differentially expressed genes at the peak of lactation. G-protein coupled receptors play roles in mediating oxytocin hormone [68], an important hormone for mammary gland cell differentiation and lactation [69][70][71]. Interestingly, both signaling pathways and upstream transcription regulator enrichments identified Notch gene family ( Table 6) and Notch signaling pathway (Figure 4) as important for milk lactose metabolism. Notch signaling is also important for mammary gland development and lactation [72], as it controls mammary epithelial cell fate [73]. Moreover, Notch 3 and 4 were also targeted by miR-874 and miR-2331-3p in the BLUE module, respectively. HOXA7, the most enriched upstream transcription regulator in the BLUE module, regulates the transcription of several genes including CD93, EDNRA, EGFR, IL7R, KCNA3, MSI2, PGR, TSC22D1, and SOX4 (a master regulator of epithelial-mesenchymal transition [74]). These results support the roles of BLUE module miRNAs in epithelial-mesenchymal transition.
The RED module showed significant correlation with two traits, SCC and lactose. However, since all members of the RED module, except miR-2285v, are novel miRNAs, no documented functions are available for these miRNAs. The most important enriched biological processes GO-term for target genes of miRNAs in the RED module was "cellular component morphogenesis". This term is defined as the process in which cellular structure is generated and organized, and so supports an important role of the RED module in cell regulation. Furthermore, biological processes GO-terms enriched for the RED module are involved in cell fate, such as "cell differentiation", "cell mobility", "cell development", and "cell migration". Further supporting evidence of the relationship between the RED module miRNAs and SCC, was through the significant enrichment (using IPA core analysis) of the planar cell polarity (PCP) pathway ( Figure 5). The PCP pathway involves a set of core molecular regulators coordinating the orientation of cells within a tissue sheet [75,76]. Interestingly, this pathway also contributes to glucose homeostasis, which is important for lactose metabolism [77]. Activation of PCP signaling locally activates Rac signaling to regulate cell fate specification events and cellular movements [78]. Indeed, Rac signaling pathway was also significantly enriched for the RED module. Notably, MYB was the most significant transcription factor enriched for target genes (including BCL2, CD4, COL4A1, FUT8, GATA3, ITPR1, KIT, KITLG, MMP1, MMP3, and SOCS2) of miRNAs in the RED module. The MYB gene is well documented to play important roles in cell growth, differentiation, and apoptosis [79].
The TURQUOISE module, just like the RED module, was significantly correlated with SCC (r = 0.54) and lactose (r = 0.57). The TURQUOISE module has two hub miRNAs in the duplex form of miR-1423p and miR-1425p. It is well documented that miRNAs can form miRNA-miRNA duplexes through reverse complementary binding events if they have completely or partially complementary structures [3,80,81]. Functional analysis of miR-142-3p indicated its importance in regulating signal to balance cell cycle processes such as balance of cell proliferation and differentiation of mesenchymal cells [82]. Interestingly, miR-223 and miR-142 showed strong interaction in the TURQUOISE module. miR-223 has been shown to up-regulate miR-142 expression through transcription factors (LMO2-L/-S isoforms and CEBP-β) [83]. Moreover, miR-223 was down-regulated after lactation peak, and might play a role in the mammary response to pathogens after parturition [32]. Similar to the BLUE module, "sensory perception of chemical stimulus" and "G-protein coupled receptor activity" GO-terms were the most significant enriched biological processes and molecular function terms for the TURQUOISE module. The role of miRNAs in the TURQUOISE module in the regulation of SCC is reflected by several cell cycle pathways enriched for its miRNAs target genes, such as cell cycle G2/M DNA damage checkpoint regulation, cell cycle regulation by BTG family proteins, and the role of CHK proteins in cell cycle checkpoint control pathways ( Figure 5). Moreover, many other relevant signaling pathways for mammary gland and lactation were significantly enriched, such as TGF-β, HIPPO, Wnt/β-catenin, epithelial adherent junction, NF-κB, integrin, CDK5, BMP and prolactin, as well as the STAT3 pathway. These pathways were also significantly enriched by target genes of differentially expressed miRNAs in our previous study [19]. Notably, STAT3 and STAT5A (key transcription factors during lactation and mammary gland involution [84][85][86][87][88][89]) were also enriched for the target genes of miRNAs in the TURQUOISE module. Interestingly, STAT3 is also directly targeted by miR-27a-5p of the TURQUOISE module.

Animal Management and Milk Sampling
Animal use and experimental procedures were according to the national codes of practice for the care and handling of dairy cattle (www.nfacc.ca) and approved by the Animal Care and Ethics Committee of Agriculture and Agri-Food Canada.
Procedures for animal management and data collection have been reported previously [19]. Briefly, nine Canadian Holstein cows, first to third parity, were used. Cows were fed a mixed ration of corn and grass silages (50:50) and concentrate, and managed following standard procedures. Animals were housed in individual tie stalls and allowed ad libitum access to feed and water at all times. Milk samples, about 50 mL composite of left and right back quarters of morning and evening milking, were collected on different days (D) in milk (D30, D70, D130, D170, D230, and D290) to represent an entire lactation curve for the measurement of milk components (fat (Fat%) and protein percentages (Protein%), lactose, milk urea nitrogen (MUN) and somatic cell count (SCC). Daily milk yield for each cow was recorded with electronic milk meters (MU-480, De Laval Inc. Kansas City, USA). Milk samples for RNA isolation (50 mL) were collected two hours after the morning milking at different times throughout the lactation curve; day 1 in milk [D1], D7, D30, D70, D130, D170, D230 and D290. Samples were placed on ice, transferred to the laboratory, and immediately processed to reduce potential RNA degradation. In the laboratory, milk was mixed well and centrifuged at 1900 g for 15 min at 4 • C. The fat layer (upper phase) was transferred to a 50 mL RNase free falcon tube and~7.5 mL Qiazol lysis reagent (Qiagen Inc., Hilden, Germany) was added, and vigorously mixed by vortexing until the fat was well dispersed. The homogenized fat was stored at −80 • C until used.

Statistical Analysis
Statistical analysis was performed with SAS version 9.3 software (SAS Institute Inc., Cary, NC, USA). The effect of day on milk yield and components was analyzed using a two-way analysis of variance in a completely randomized design with repeated measures. The effect of day was tested with multiple comparisons to D30 using a Dunnett correction. Quarter (cow) was used as the subject of the repeated analysis, and the best variance-covariance matrix was selected using the AICC (Corrected Akaike Information Criteria) criteria for each variable analysed. ARH (I) variance-covariance matrix was used for all variables, except that CSH was used for lactose.
The statistical model included the fixed effects of quarter and day.
where Y ijk = observation for quarter i sampled on day j for animal k, µ = general mean, q i = fixed effect of quarter i = left, right, cow (q) k (i) = random effect (subject of the repeated measure), d j = fixed effect of time of sampling (day) j, qd ij = interaction between quarter and day of sampling, and e ijk = residual error. Probabilities lower than 0.05 (p < 0.05) were declared significant.

Total RNA Isolation, miRNA Sequencing, and Bioinformatics Management of Data
Procedures for RNA isolation, miRNA sequencing, and bioinformatics management of data have been reported previously [19]. In brief, total RNA was extracted from fat homogenate using miRVana miRNA isolation kit (Life Technologies, Carlsbad, CA, USA). RNA was digested with Turbo DNase (Ambion Inc., Austin, TX, USA) and purified using Zymo RNA clean and concentrator-25 kit (Zymo Research, Irvine, CA, USA). The integrity of RNA (RIN) was determined on an Agilent 2100 Bioanalyzer using RNA 6000 Pico kit (both from Agilent Technologies, Richardson, TX, USA). Isolated total RNA with RIN values from 2.3 to 8.5 and having an intact small RNA fraction on an Agilent 2100 Bioanalyzer electropherogram were used for library preparation [19]. A total of 71 libraries were prepared and subjected to 50 bp single end sequencing on an Illumina HiSeq 2000 system (Illumina Inc., San Diego, CA, USA) using TruSeq v3 reagents. The identification of known miRNAs and discovery of novel miRNAs were performed using miRDeep2 v2.0.0.8 [90], which uses a probabilistic algorithm based on the miRNA biogenesis model and designed to detect miRNAs from deep sequence reads. The expression of miRNAs (count table) was normalized using Deseq2 package (v1.11.19) [91]. Deseq2 models read count data with a negative binomial distribution. The normalized data was extracted using count () function and a final normalized matrix containing 713 miRNAs (475 known and 238 novel) [19] was used for network construction.

Gene Co-Expression Network Analysis
The weighted gene correlation network analysis (WGCNA) R-package [92] was used for miRNA co-expression network analysis. The input for co-expression analysis was normalized data of 713 miRNAs obtained from miRNA expression data as described above and in details in Do et al. [19]. To compute the co-expression network for whole lactation data, an adjacency matrix was generated by calculating Pearson's correlation between all miRNAs, and raising it to a power β (soft threshold) of 7. The power value was chosen using a scale-free topology criterion (r 2 = 0.95). Then, miRNAs were clustered using degree of overlap in shared neighbors between them, a topological overlap measure (TOM) was calculated and a value between 0 and 1 was assigned to each miRNA pair. When miRNAs share the same neighbors, a TOM value of 1 is assigned, and when they do not share any neighbor, a TOM value of 0 is assigned. An average linkage hierarchical clustering using the dynamic tree-cutting algorithm to produce a clustering tree (dendrogram) [93] was performed. Each branch of the tree is a module, and modules with at least 20 miRNAs were assigned a color. Details about methodology and the relative merits of WGCNA have been provided previously [92][93][94].
The module-trait relationship was used to select potential biologically interesting modules for downstream analysis. Module-trait relationship was computed based on Pearson's correlation between the module eigengene values and milk yield and components. The eigengene value is defined as the first principal component of a given module, and it represents a measure of gene expression profiles in the module. A module was chosen for further analysis if it had a value of module-trait relationship >|0.5| for any phenotype. Furthermore, miRNAs in selected modules were used for functional enrichment analysis only if their intra-modular connectivity with the module was >0.8. Intra-modular connectivity measures how co-expressed a given gene or miRNA is with the other miRNAs (members) within the module, and can also be called module membership.

Function Enrichment of Target Genes of miRNAs in Significant Modules
TargetScan was used to predict the target genes of miRNAs (known and novel) in modules significantly associated with milk traits. Perl scripts from the TargetScan website (www. targetscan.org) were used to predict (targetscan_60.pl), and to calculate, the context scores (targetscan_61_context_scores.pl) of target genes. Predicted target genes with context + scores above 95th percentile were further used [18]. The target genes of module specific miRNAs were enriched for gene ontology (GO) terms using ClueGO [95], and pathways and upstream transcription regulators using Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems Inc., Redwood City, CA, USA). For ClueGO analysis, a hypergeometric test was used for GO-term enrichment and Benjamini-Hochberg [23] correction for multiple testing-controlled p-values. For IPA core analysis, pathways were considered significantly over-represented at a Benjamini-Hochberg corrected p-value ≤ 0.05, and contained at least two target genes from data input.

Conclusions
Overall, this study identified regulatory networks and related mechanisms by which miRNAs interact with each other to regulate lactation phenotypes. Important hub miRNAs, transcription factors and regulatory networks for milk traits were uncovered. Moreover, the enrichment of important signaling pathways for milk yield and milk components enhances our knowledge about the regulation of milk composition at the molecular level. In addition, miRNAs demonstrated various ways of interaction, including shared common target genes, enriched transcription factors and regulatory pathways as well as formation of miRNA-miRNA duplexes in the regulation of lactation phenotypes. Further functional characterization of important module miRNAs and hub genes will further understanding of their roles, as well as inform of their potential use as biomarkers in selection programs for high milk yield or milk with specific requirements.