Co-Expression Network Analysis Identifies miRNA–mRNA Networks Potentially Regulating Milk Traits and Blood Metabolites

MicroRNAs (miRNA) regulate mRNA networks to coordinate cellular functions. In this study, we constructed gene co-expression networks to detect miRNA modules (clusters of miRNAs with similar expression patterns) and miRNA–mRNA pairs associated with blood (triacylglyceride and nonesterified fatty acids) and milk (milk yield, fat, protein, and lactose) components and milk fatty acid traits following dietary supplementation of cows’ diets with 5% linseed oil (LSO) (n = 6 cows) or 5% safflower oil (SFO) (n = 6 cows) for 28 days. Using miRNA transcriptome data from mammary tissues of cows for co-expression network analysis, we identified three consensus modules: blue, brown, and turquoise, composed of 70, 34, and 86 miRNA members, respectively. The hub miRNAs (miRNAs with the most connections with other miRNAs) were miR-30d, miR-484 and miR-16b for blue, brown, and turquoise modules, respectively. Cell cycle arrest, and p53 signaling and transforming growth factor–beta (TGF-β) signaling pathways were the common gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched for target genes of the three modules. Protein percent (p = 0.03) correlated with the turquoise module in LSO treatment while protein yield (p = 0.003) and milk yield (p = 7 × 10−04) correlated with the turquoise model, protein and milk yields and lactose percent (p < 0.05) correlated with the blue module and fat percent (p = 0.04) correlated with the brown module in SFO treatment. Several fatty acids correlated (p < 0.05) with the blue (CLA:9,11) and brown (C4:0, C12:0, C22:0, C18:1n9c and CLA:10,12) modules in LSO treatment and with the turquoise (C14:0, C18:3n3 and CLA:9,11), blue (C14:0 and C23:0) and brown (C6:0, C16:0, C22:0, C22:6n3 and CLA:10,12) modules in SFO treatment. Correlation of miRNA and mRNA data from the same animals identified the following miRNA–mRNA pairs: miR-183/RHBDD2 (p = 0.003), miR-484/EIF1AD (p = 0.011) and miR-130a/SBSPON (p = 0.004) with lowest p-values for the blue, brown, and turquoise modules, respectively. Milk yield, protein yield, and protein percentage correlated (p < 0.05) with 28, 31 and 5 miRNA–mRNA pairs, respectively. Our results suggest that, the blue, brown, and turquoise modules miRNAs, hub miRNAs, miRNA–mRNA networks, cell cycle arrest GO term, p53 signaling and TGF-β signaling pathways have considerable influence on milk and blood phenotypes following dietary supplementation of dairy cows’ diets with 5% LSO or 5% SFO.


Introduction
Bovine milk and its products constitute a rich source of proteins, energy, minerals (e.g., calcium), vitamins (A, B, D, E and K) and antioxidants in human nutrition. Milk supplies unsaturated fatty acids (USFA) which have been associated with decreased risk of cardiovascular diseases (stroke, high blood pressure, heart failure and coronary heart diseases), inflammatory diseases and some types of cancers [1][2][3]. Unsaturated fatty acids make up about 30% of the fatty acid content of milk, meanwhile it has been proposed that milk fat composition with potential positive effects on human health should contain about 70% USFA [4]. Nutrition is one of the factors that greatly impacts milk fat composition and the largest changes in milk fatty acid composition have been obtained either by changing the amounts and the nature of forages in the diets of cows, particularly pasture, or by adding plant or marine oils to the diet [5,6]. Plant products like linseed, soybeans, safflower and sunflower are the most effective sources of unsaturated plant lipids used to enhance the conjugated linoleic acid (CLA) and USFA contents of milk fat. Unsaturated fatty acids and other factors like physiological and metabolic state of the cow, breed and genetics are known to influence the concentration of blood metabolites like glucose, nonesterified fatty acids (NEFA), triacylglyceride (TAG) and β hydroxybutyric acid [7,8]. For instance, we reported significant increases in blood NEFA and TAG concentrations and significant reductions in milk fat and milk urea nitrogen contents in Holstein cows following dietary supplementation with USFA [8]. Moreover, the blood metabolic profile of dairy cows is used to assess the nutritional and health state of the dairy herd [9,10].
In our previous transcriptome studies of the bovine mammary gland, we identified mRNAs and miRNAs that were differentially expressed in response to diets rich in USFA [11,12]. Furthermore, we also examined the effect of diets rich in USFA on milk composition (fat, protein, milk yield and lactose) and blood metabolites (TAG and NEFA) of lactating Holstein cows [7]. The mRNA transcriptome analysis identified 1006 (460 up and 546 downregulated) and 199 (127 up and 72 down-regulated) genes that were significantly differentially regulated after linseed oil (LSO) and safflower oil (SFO) supplementation, respectively, meanwhile the miRNA transcriptome analysis detected 14 and 22 miRNAs significantly differentially regulated by LSO and SFO, respectively. Since a network of genes and regulatory factors work in concert to influence the phenotypic expression of traits, assessment of gene expression without taking into account the factors that regulate their activities may not adequately explain the complex biological mechanisms underlying the expression of traits. MiRNAs interact with mRNA(s) to regulate their (mRNA(s)) expression and consequently biological processes, so it is important to study their synergistic effects on the phenotypic expression of traits. Hence, an integrative approach in assessing gene expression in a network basis is necessary to unravel the molecular mechanism underlying milk fat traits.
Network approaches have proven to be powerful tools for exploring the biological mechanisms underlying complex traits [13][14][15]. It has been widely applied on output data from high throughput sequencing technology to identify key regulators and pathways in human diseases such as cancer and obesity [16], type 1 diabetes [17], Alzheimer's disease [18], livestock production traits [13,[19][20][21], and functional annotation of cattle genes [22]. Moreover, understanding gene networks also helps to better guide genomic selection in animal breeding programs [23]. In order to understand gene interaction, different methods have been developed to construct co-expression networks and to identify modules of highly connected genes. The weighted gene co-expression network analysis (WGCNA) is among the most established and widely used of such methods [24]. We and other authors have successfully used WGCNA to identify key genes and networks for various complex traits in livestock species such as meat quality traits in pigs [25], feed efficiency in cattle [26] and milk yield and component traits in cows [27]. Furthermore, integrative omics approaches have been applied on combined mRNA and miRNA expression data to detected major regulatory mechanisms in different phenotypes such as carcass and meat quality traits in porcine [25], abnormality in breast cancer patients [28], and colorectal cancer [29]. Moreover, the consensus module approach (finding common functions/processes) has proven to be a promising method for finding hub genes and regulators across different datasets [30][31][32][33][34][35]. Such hub genes and regulators may form targets for further functional validation of their roles in identified networks. Furthermore, the identification of key miRNAs, their networks, and their downstream target genes and pathways might also facilitate the use of genetic engineering technologies, such as RNA interference technologies or gene editing, to obtain desired phenotypes by controlling the expression of miRNAs and/or their target genes. Moreover, the miRNAs, genes, and network information might be useful for genomic predictions [36,37]. However, these approaches have not yet been applied to explore the regulatory mechanisms in the bovine mammary gland in response to diets rich in USFA (SFO or LSO). Therefore, this study aimed to (i) construct consensus modules across miRNA expression data from control, LSO and SFO treatments using the WGCNA approach; (ii) correlate important miRNA modules with milk and blood component phenotypes; (iii) enrich target genes (mRNA) of miRNAs from important modules to explore the possible biological processes, pathways and transcriptional regulators regulating milk and blood component phenotypes and (iv) identify miRNA-mRNA networks regulating milk and blood component phenotypes.

Phenotypic Data
A summary of the data on blood metabolites, milk and component yields including fatty acid profiles for the control and treatment periods used for co-expression and network analyses is shown in Table 1.

Identification of Consensus Modules and Module Trait Relationship
Using the WGCNA approach for miRNA read data described by Li et al. [12], we identified a total of three consensus modules (blue, brown, and turquoise) of co-expressed miRNAs during the control and treatment periods (Figures 1 and 2). The modules were made up of 70 (blue), 34 (brown), and 86 (turquoise) miRNA members (Figures 1 and 2). The grey module grouped miRNAs with no coherent co-expression patterns; therefore it was not further discussed. MiRNAs were further selected based on their intra modular connectivity or eigengene-based connectivity (k.ME). The k.ME is a measure of how a miRNA is correlated to module eigengene and miRNAs with high k.ME values (>0.6) are better representatives of module characteristics [38]. Therefore, miRNAs with k.ME > 0.6 were selected for downstream analyses. A total of 18, 12 and 19 miRNAs with k.ME > 0.6 in the blue, brown, and turquoise modules, respectively, were selected for downstream analyses (Table 2). Hub miRNAs or miRNAs with the most connections with other members of the module were bta-miR-30d, bta-miR-484 and bta-miR-16b for the blue, brown, and turquoise modules, respectively (Table 2).
In LSO treatment, the turquoise module correlated significantly with protein percent (p = 0.03) while protein yield (p = 0.003) and milk yield (p = 7 × 10 −4 ) correlated with the turquoise module, protein and milk yields and lactose percent (p < 0.05) correlated with the blue module and fat percent (p = 0.04) correlated with the brown module in SFO treatment ( Figure 1). Several fatty acids correlated (p < 0.05) with the blue (CLA:9,11) and brown (C4:0, C12:0, C22:0, C18:1n9c and CLA:10,12) modules in LSO treatment and with the turquoise (C14:0, C18:3n3 and CLA:9,11), blue (C14:0 and C23:0) and brown (C6:0, C16:0, C22:0, C22:6n3 and CLA:10,12) modules in SFO treatment ( Figure 2). Several measured parameters also correlated with the identified modules in the control samples (Figures 1 and 2).   The weighted gene co-expression network analysis (WGCNA) was used to group miRNAs into consensus modules based on their expression patterns. Four consensus modules were identified and each consensus module eigengene was tested for correlation with blood and milk parameters during (A) control period, (B) linseed oil and (C) safflower oil treatments. Correlation coefficients and corresponding p-values (in brackets) between turquoise, blue and brown modules in the y-axis, and blood and milk parameters in the x-axis. The module−trait relationship matrix is colored based on the intensity of the correlation: red is a strong positive correlation, while green is a strong negative correlation.
p-values (in brackets) between turquoise, blue and brown modules in the y-axis, and blood and milk parameters in the x-axis. The module−trait relationship matrix is colored based on the intensity of the correlation: red is a strong positive correlation, while green is a strong negative correlation.

miRNA Target Gene Prediction and Enrichment Analysis
Using TargetScan, we identified 3199, 3727, and 4045 target mRNAs for 18, 12, and 19 miRNAs in the blue, brown, and turquoise modules, respectively (Table S1a,c,e). Amongst them, 1311, 1533, and 1697 mRNAs from mRNA data of the same samples [11] had significant negative correlations (FDR < 0.05) with 18, 12 and 19 miRNAs in the blue, brown, and turquoise modules, respectively (Table S1b,d,e). These mRNAs (filtered target mRNAs) were used as input for enrichment analyses for GO, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and transcription factors.

Blue Module miRNAs and Their Potential Roles
The blue module miRNAs were significantly correlated with milk and protein yields during the control period, milk yield, protein percent and lactose during treatment with SFO, protein percent and protein yield with LSO treatment. The miRNA expression data of these animals revealed that several of the blue module miRNAs were differentially expressed in response to dietary supplementation with SFO (bta-miR-96, miR-99a-5p, miR-199b, miR-16a, miR-484 and miR-99b) and LSO (bta-miR-885, miR-23b-3p and miR-99a-5p) [12], thus supporting their relationship to milk yield and milk components.
Many of the blue module miRNAs have been reported to play diverse roles in many biological processes. For example, human homologue of bta-miR-191 has been found to be dysregulated in different types of tumors in humans including colorectal [46], breast and prostate cancers [47], while miR-199b is involved in acute myeloid leukemia [48,49] and breast cancer metastasis [50][51][52]. The implication of these miRNAs (bta-miR-191 and miR-199b) in the different types of cancers suggests roles in cellular functions in the mammary gland. The development of breast and prostate cancers is linked to cellular processes like cell death or apoptosis, therefore suggesting a link between these miRNAs and the milk phenotypes in this study.
Some blue module miRNAs like bta-miR-183 and bta-miR-2284b correlated negatively with many candidate genes for protein yield (CTDSP1, DGCR2, HLTF, ICA1, MTA1, SESN1, SPRY2, SRSF2, UTP6 and ZFAND5), milk yield (C14orf28, QARS, SLC20A1, HNRNPA1, MAFF, MGME1, PPP2R5C and RHPN2) and C17:0 fatty acid (ILF2 and SFT2D1) (Figure 4) thus suggesting potential roles in the regulation of these genes. In addition to the report of differential expression of some blue module miRNAs between lactation stages [39], the importance of blue module miRNAs for milk yield and components is further supported by results of pathways enrichments of their target mRNAs. For instance, ErbB, MAPK, Wnt, TGF β and Hippo signaling pathways are involved in the regulation of mammary gland development and lactation processes (reviewed in [53]). Furthermore, ErbB and TGF β signaling pathways have been associated with lactation persistency in Holsteins [54]. The p53 signaling pathway, the second most enriched pathway for the blue module, functions by preventing cancer formation and hence acts as a tumor suppressor [55]. The p53 gene has been nicknamed "guardian of the genome" due to its role in maintaining the stability of the genome [56]. In addition, this pathway regulates proper cellular differentiation and development, and it is also important for tissues undergoing postnatal development [57]. Furthermore, inappropriate expression of the p53 signaling pathway within the mammary epithelium of transgenic mice caused apoptotic cell death of the alveolar epithelium of the mammary gland [58]. In this study, enrichment of p53 signaling pathway supports the significant correlation of the blue module miRNAs with milk yield, protein yield and lactose percent.
Moreover, we also reported several enriched GO terms for the blue module filtered target mRNAs such as vesicle docking, negative regulation of transcription from RNA polymerase II promoter and proteasome−mediated ubiquitin-dependent protein catabolic process. However, it is not clear how these GO terms are linked to the studied phenotypes. We identified 22 important transcription factors which could mediate the functions of miRNAs in the regulation of the target mRNAs or phenotypes (Table 5). SMAD4, SP1 and EGR1 were the top three transcription factors enriched for the target mRNAs of the blue module miRNAs. SMAD4 is a tumor suppressor gene and it is essential for transforming growth factor beta (TGFβ) signaling [59], and plays important roles in cell differentiation, growth and apoptosis [60]. Other important transcription factors such as STAT3 and PPARG are well known to regulate milk and milk fat synthesis [61][62][63].

Brown Module miRNAs and Their Potential Roles
Bta-miR-484, a member of the brown module was previously reported as differentially expressed due to LSO and SFO treatments and it is also important in the regulation of lactation signaling [40]. Correlation analyses indicated that this miRNA negatively correlated with previously reported candidate genes for fat percentage (EIF1AD and NUDT16), C22:6n3 (CPPED1) and C16:0 (LY6E, DOLPP1 and QDPR). The effects of LSO and SFO supplementation reduced fat percentage in the studied animals by 30.38% and 32.42%, respectively [11] thus supporting the present observation. Additionally, bta-miR-484 has been observed to prevent cell proliferation and epithelial-mesenchymal transition process by targeting both ZEB1 and SMAD2 genes, thus functions like a tumor suppressor and may serve as a prospective biomarker for cervical cancer [64]. Other members (bta-miR-26b and bta-miR-107) of the brown module have functions related to cellular processes [40,53] which are essential for lactation processes or milk synthesis. The human homologue of bta-miR-26b was shown to play a protective role in the etiology of breast cancer by promoting apoptosis through targeting SLC7A11 [65] while bta-miR-107 is associated with mammary stem cell activities [66].
The most significantly enriched biological process GO term in the brown and turquoise modules, GDP binding, is involved in cell proliferation, signal transduction, protein synthesis, and protein targeting [67] (Table 3). Other enriched pathways in the brown module like MAPK signaling pathway plays an important part in numerous cellular processes such as apoptosis, proliferation and differentiation [68], stress responses, and immune defense [69,70] and have been noted to be important for mammary gland development and milk secretion in caprine [71].
The most enriched transcription factor for the brown module was specific protein 1 (SP1) known to regulate the expression of numerous genes involved in cell proliferation, apoptosis and differentiation, and increase in its transcriptional activities is associated with tumorigenesis [72].

Turquoise Module miRNAs and Their Potential Roles
The hub miRNA (bta-miR-16a) for the turquoise module has been reported to be differentially expressed in response to SFO treatment [12]. Bta-miR-16a being differentially regulated by SFO and also having the highest intra modular connectivity suggests involvement in the regulation of the traits (C14:0, C18:3n3 and 9, 11-CLA) that were significantly correlated with the turquoise module in SFO treatment. Although a direct role for this miRNA in mammary gland functions has not been demonstrated, a previous study suggests its involvement in tumor suppression through inhibition of cell cycle progression [73]. Some turquoise module miRNAs like bta-miR-130a and bta-miR-142-5p have been linked to milk fat synthesis [40,74] and disease parthenogenesis [75][76][77]. Additionally, overexpression of bta-miR-130a affects cellular triacylglyceride synthesis in bovine mammary epithelial cells via regulation of PPAR-γ [74], fatty acid storage and glucose metabolism. Besides, some predicted target genes of miR-130a have been associated with neurodevelopmental disorders such as autism, schizophrenia and hereditary spastic paraplegia [78]. GO enrichment indicated that target mRNAs of turquoise module miRNAs participate in many different processes such as GTP binding, GDP binding, transcription coactivator activity, membrane organization, protein ubiquitination and regulation of apoptotic processes. Several KEGG pathways such as Ubiquitin mediated proteolysis, TGF-β signaling pathway and cell cycle, and transcription factors such as STAT3, SMAD4, SP1 and EGR1 with roles in many different processes involving milk production and related traits [39,63] were enriched for target mRNAs of turquoise module miRNAs. TGF-β signaling pathway and p53 signaling pathway were common pathways enriched by target mRNAs from all three modules and their roles in related traits have been discussed above. Ubiquitin mediated proteolysis, the most enriched pathway for turquoise module is important for protein degradation [77] and many processes including mediation of lactation signal in skeletal muscle of dairy cows [76]. Top transcription factors enriched for the target mRNAs of turquoise module were SMAD4, SP1, and EGR1.The early growth response protein 1 (EGR1) acts as a transcription regulator of target genes, hence play roles in the regulation of cell survival, proliferation, cell death response to growth factors, DNA damage, and ischemia [79]. A notable enriched transcription factor was STAT3 with known association with milk production [75], fertilization and embryonic survival rates [80] in dairy cows.

Association Between miRNA and mRNA with Expressed Phenotypes
The most interesting part of this analysis is the different pathways linking miRNAs to the actual phenotypes through their target genes (Table 6 and Figure 4). Interestingly, different miRNAs that shared the same target genes regulated different traits. For instance, by targeting TP53, bta-let-7b might influence milk yield but also, bta-miR-96 might have an influence on protein percentage. Since many genes (mRNAs) and miRNAs influencing milk yield have been characterized such as TP53, GRHL2, bta-miR-183, bta-let-7b and bta-miR-96, we will not discuss about the network of miRNAs influencing milk yield. Meanwhile, for the first time, the link between some milk fatty acids, genes (mRNAs) and miRNAs have been reported. Interestingly, we observed negative correlation between bta-miR-484 and three different genes (QDPR, LY6E and DOLPP1) and positive correlation with C16:0 concentrations in milk. The roles of bta-miR-484 have been reported above, while it is not clear how these three genes are involved in the metabolism of C16:0. However, DOLPP1 has a potential role in regulating subcutaneous fat [81] while LY6E might play a role in the regulation of glycosylphosphatidylinositol. Also, QDPR is important for fat traits in pigs [82]. Therefore, these genes might be interesting candidates for milk fat traits. The influence of bta-miR-183 on protein yield and protein percentage suggests that it can down regulate 8 different genes to influence protein yield. Nevertheless, these connections are based on correlations so it might not reflect true associations; therefore more studies are needed to validate the identified links and how they participate in mammary gland response to dietary USFA.

Animal Management and Sampling
Animal management and sampling procedures were according to the national codes of practice for the care and handling of dairy cows (http://www.nfacc.ca/codes-of-practice) and approved by the Animal Care and Ethics Committee of Agriculture and Agri-Food Canada (CIPA#402, 04 April 2012).
Detailed procedures for animal management and data collection have been reported previously [12]. Briefly, 12 Canadian Holstein cows in mid-lactation were randomly assigned to either LSO or SFO treatment (6 cows/treatment). These animals were fed a control diet of total mixed ration of corn and grass silages (50:50) and concentrates for 28 days (control period), after which the control diet was supplemented with 5% LSO or 5% SFO on a dry matter basis (treatment period) for another 28 days. Animals were milked every day at 8:00 am and 6:00 pm. Analysis of fat, protein and lactose contents in milk samples collected on days −14 and −1 (control period) and +7, +14, +21 and +28 (treatment period) were done using 80 mL of milk (pool of morning (40 mL) and evening (40 mL) milk) by a commercial laboratory (Valacta Laboratories Inc., Ste. Anne de Bellevue, QC, Canada). Daily milk yield for each cow was recorded with electronic milk meters (MU-480, De Laval Inc., Kansas City, MO, USA). Milk fat from milk samples (40 mL) was extracted by centrifugation at 4500× g for 20 min at 4 • C.
Blood samples were aseptically collected from animals on days −14 and −1 (control period) and +7, +14, +21 and +28 (treatment period) and centrifuged at 7500× g for 20 min at room temperature. The resulting plasma was used for the analysis of non-esterified fatty acid (Wako Chemicals, Kit HR series NEFA-HR, Richmond, VA, USA) and triacylglyceride (enzychrom TAG assay kit, Bioassay System, Hayward, CA, USA), following manufacturers' instructions.
Mammary biopsies were collected from animals in each group on day −14, day +7 and day +28 which corresponded to middle of control period, early treatment and end of treatment periods, respectively, following an established protocol [83]. Biopsies were snap frozen in liquid nitrogen and stored at −80 • C pending isolation of total RNA.

RNA Isolation
Procedures for RNA isolation have been reported previously [11]. In brief, 50 mg of mammary gland biopsy sample was used for total RNA isolation using miRNeasy Kit (Qiagen Inc., Toronto, ON, Canada) following manufacturer's instructions. The Turbo DNase Kit (Ambion Inc., Foster City, CA, USA) was used to remove contaminating DNA from isolated RNA. The Nanodrop ND-1000 instrument (NanoDrop Technologies, Wilmington, DE, USA) was used to measure RNA concentration and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) was used to accessed RNA quality. The RNA Integrity Numbers of all the samples were >8.

mRNA Sequencing and Data Processing
Procedures for mRNA sequencing and data processing have been reported previously [11]. Briefly, 250 ng of total RNA/sample was use for library generation using the TruSeq stranded mRNA Sample Preparation Kit (Illumina Inc.

miRNA Sequencing and Data Processing
Procedures for miRNA library preparing, sequencing and bioinformatics management of data have been reported previously [12]. Briefly, the procedure for miRNA library preparation and barcoding for sequencing was according to Vigneault, et al. [84] with slight modifications [12]. Total RNA was first ligated to a primer (adaptor) at the 3 by T4 RNA Ligase 22tr K227Q (New England Biolabs Inc., Canada) and 5 ends by T4 RNA Ligase 1 (Enzymatics Inc., a division of Qiagene Inc., Beverly, MA USA) followed by reverse transcription into cDNA using Superscript III Kit (Life Technologies, Carlsbad, CA, USA). Barcoding of the different libraries was done followed by size separation by polyacrylamide gel electrophoresis and finally the concentration of the The FastQC program version 0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/) was used to check sequencing quality. Then, the cutadapt v1.2.2 program (http://code. google.com/p/cutadapt/) was used to trim adaptor sequences. Clean reads were parsed into one and mapped to the bovine genome (Bta_4.6.1) using bowtie 1.0.0. [85]. Reads that mapped to 1 to 5 positions were further used. miRDeep2 v2.0.0.5 tool was used to identify known miRNAs and also in the discovery of novel miRNAs.

Fatty Acid Analysis
Fatty acid methyl ester preparation and quantification of fatty acid profiles have been reported previously [11]. Briefly, fatty acid methyl esters were prepared and quantified using the Hewlett Packard 6890 N gas chromatographic system (Agilent Technology, Wilmington, DE, USA). The carrier gas was hydrogen and the capillary column used was SLB-IL111 (100 m × 0.25 mm, 0.2 µm in thickness, Supelco, Bellefonte, PA, USA). The oven temperature was set at 40 • C for 1 min followed by 80 • C to

Construction of Gene Co-Expression Networks
The expression of miRNAs was normalized using Deseq2 package (v1.11.19) [86] and the final normalized matrix of 321 miRNAs were used as input for co-expression network analysis using the WGCNA R-package [24]. For WGCNA analysis, a signed co-expression measure for each pair of miRNAs was computed based on their co-expression level. Then weighted adjacency matrix was calculated from the signed co-expression measure using a power function. A topological overlap measure (TOM) was then calculated based on a combination of a value between the adjacency of two miRNAs and the connection strength that the two miRNAs share with other miRNAs. A TOM of 0 or 1 was assigned to each pair of miRNAs. When miRNAs shared the same neighbor, a TOM value of 1 is assigned while a TOM value of 0 indicates that they do not share any neighbor. To produce a clustering tree (dendrogram), the dynamic tree−cutting algorithm was used [87]. To construct the consensus module, the blockwiseConsensusModulesfunction was run [88,89] with option of soft-thresholding power for network construction of 9, and minimum module size of 20. Moreover, the medium threshold was also applied to control the sensitivity of module detection (deepSplit of 2) and to merge modules in the dynamic tree (mergeCutHeight of 0.25). The signed network option was chosen when constructing the consensus modules. In the gene network, a gene might interact with many others to perform its function, therefore, a minimum module size of 30 genes has been recommended for the gene network construction, by the software developers [24]. Since a lower number of miRNAs might interact with each other to form networks [90,91], we applied a lower threshold of 20 miRNAs for minimum module size. Each branch of a tree is a module and a module with at least 20 miRNAs was assigned to a color ( Figure S1). Details about WGCNA and its merits have been reported previously [24,89,92].

Module−Trait Relationship
Module−trait relationships were computed based on Pearson's correlation between the module eigengene and blood and milk components data. The eigengene is defined as the first principal component of a given module and it represents a measure of miRNA expression profiles in the module. A module was chosen for further analysis if it presented a module−trait relationship > |0.5| and a p-value < 0.05. Potential biologically interesting (significant) modules were selected for downstream analysis. Furthermore, miRNAs in selected modules were used for functional enrichment analysis if eigengene−based connectivity (k.ME), a measure of how the miRNA is correlated to module eigengene, was > 0.6 (k.ME > 0.6) [26]. A k.ME > 0.6 indicates higher connectivity and thus higher representation of modular functions.

Predicted Target mRNAs of miRNAs
In order to investigate the function of miRNAs in module significantly correlated with traits, we first predicted their target mRNAs. The perl scripts from the TargetScan website (http://targetscan.org) were used to predict target mRNAs (targetscan_60.pl) and also to calculate their context scores (targetscan_61_context_scores.pl). TargetScan computes the context++ score for a specific site as the sum of the contribution of 14 features of the miRNA, miRNA site, or mRNA (including the mRNA surrounding sequence) (http://www.targetscan.org/vert_70/docs/context_score.html) to define sites on mRNAs most effectively targeted by miRNAs [93]. Predicted target mRNAs with context ++ scores above 95th percentile were further used [12,27,39]. The predicted target mRNAs were then filtered against the mRNA expression data from the same animals [11]. Only target mRNAs that were present in the mRNA expression data were retained for further analysis.

Co-Expression Analysis of miRNA-mRNA Expression
For miRNA-mRNA co-expression, the Pearson correlation coefficient between target mRNAs and miRNAs were calculated. A miRNA-mRNA pair was considered co-expressed if it had a negative and significant correlation value at FDR < 0.05. To further explore how miRNAs contributed to particular traits, we examined the correlation between miRNAs and their mRNA targets with the phenotypes in each significantly correlated module. Important interactions between miRNA and mRNAs were visualized using Cytoscape [94].

Gene Ontologies, Pathways and Transcription Factors Enrichment
Functional enrichment of GO terms of target mRNAs was performed for each selected module using EnrichR [95,96]. EnrichR presents results according to hierarchy and relationship between terms which facilitates the interpretation of results. In this enrichment, the p-values for each term were adjusted using Benjamini-Hochberg (BH) correction [95]. Gene ontology terms, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and transcription factors were considered significantly enriched at adjusted p < 0.05 (FDR < 0.05).

Conclusions
In this study, three consensus modules (blue, brown, and turquoise) composed of 70 (blue), 34 (brown) and 86 (turquoise) miRNA members were identified. We also demonstrated how miRNAs in these modules interacted with mRNAs to influence blood and milk phenotypes following dietary supplementation with USFA. Hub miRNAs for the blue, brown, and turquoise modules were bta-miR-30d, bta-miR-484 and bta-miR-16b, respectively. Turquoise module had the most significant correlations with several traits including protein percentage in LSO treatment, protein yield, milk yield, C14:0, C18:3n3n and 9, 11-CLA in SFO treatment. The association of miRNA modules with milk and blood phenotypes has provided information about miRNA modules, hub miRNAs, GO terms, transcription factors and pathways that are involved in the regulation of blood and milk parameters following dietary supplementation with diets rich in USFA. This study will contribute to the molecular understandings of the co-expression patterns of miRNAs, miRNA-mRNA, and regulatory activities in the bovine mammary gland following dietary supplementation with USFA.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/19/ 9/2500/s1. Table S1a-g: Predicted target genes (mRNAs) for miRNA members of the (a) blue, (c) brown, and (e) turquoise modules. Predicted filtered target genes for miRNA members of the (b) blue, (d) brown and (f) turquoise modules (these are predicted target mRNAs of the miRNAs that were present in the mRNA transcriptome data of the same animals and negatively correlated with any of the miRNAs in the blue, brown, and turquoise modules) (predicted target mRNAs that were not present in the mRNA transcriptome of the same data were not further used). (g) Specific mRNA-miRNA pairs in the blue, brown, and turquoise modules (each pair contains a miRNA that is negatively correlated with its target gene. Figure S1: The dynamic cut tree obtained from weighted co-expression network analyses. Funding: This study was funded by Agriculture and Agri-Food Canada (grant number J-000733).