MicroRNA-mRNA regulatory network related to lipid metabolism in bovine mammary epithelial cells

Milk fat percentage is an important factor of milk quality in dairy cattle. Functional microRNAs and genes can affect lipid synthesis and metabolism through differential expression in bovine mammary epithelial cells (BMECs). It is necessary to screening the crucial candidate gene and miRNA on milk fat percentage. In this study, we extract total RNA of BMECs isolated from Chinese Holstein cows with high and low milk fat percentages for the conjoint analysis of RNA-seq and Solexa sequencing data. 190 differentially expressed genes and 33 differentially expressed microRNAs (DERs) were enriched in 488 GO terms and 12 KEGG pathways significantly (p <0.05) based on the conjoint analysis. The detection of triglyceride production in BMECs showed that bta-miR-21-3p and bta-miR-148a promote triglyceride synthesis, whereas bta-miR-124a, bta-miR-877, bta-miR-2382-5p and bta-miR-2425-5p inhibit triglyceride synthesis. Meanwhile, the target relationships between PDE4D and bta-miR-148a, PEG10 and bta-miR-877, SOD3 and bta-miR-2382-5p, and ADAMTS1 and bta-miR-2425-5p were verified using luciferase reporter assays and quantitative RT-PCR. The conjoint analysis can more accurately screen candidate regulator related to milk fat percentage at the molecular level, which provided a scientific research method screening functional miRNA and gene for the breeding of new high-quality dairy cows.

research area of dairy cattle genetics and breeding [1,2]. The mammary gland is composed of lactation tissue and connective tissue, has the physiological function of synthesizing and secreting milk in dairy cattle [3]. Bovine mammary epithelial cells (BMECs) cultured in vitro can still maintain the ability to synthesize and secrete milk and can be used as an experimental cell model. Understanding the mechanism of fat biosynthesis and metabolism in BMECs can provide a theoretical basis for further studies on the genetic improvement and breeding of new dairy cattle.
MicroRNA (miRNA) is a class of single-stranded RNA molecules encoded by endogenous genes with a length of approximately 22 nucleotides [4]. MiRNAs cannot encode proteins themselves, but they can regulate a variety of biological processes by recognizing the 3'untranslated region (UTR) of downstream target genes to degrade the mRNA or obstruct protein translation [5,6]. The differential expression of miRNAs in the mammary gland at different developmental stages of different cattle species indicates that it may play an important role in regulating fat deposition in milk [7,8]. It has been found that most mRNAs in milk whey are present in exosomes, whereas miRNAs in milk whey are present in the supernatant as well as exosomes according to miRNA and mRNA microarray analyses of bovine raw milk and total RNAs purified from exosomes (prepared by ultracentrifugation) and ultracentrifuged supernatants [9]. To screen the regulatory factors that influence traits more accurately, the conjoint analysis of miRNAs and mRNAs can obtain candidate miRNAs and target genes more effectively to clarify the regulatory mechanism of traits more clearly [10]. For example, to understand the molecular mechanism underlying intramuscular fat (IMF) deposition after castration in beef cattle, there was an investigation of differentially expressed miRNAs (DERs) and mRNAs in the IMF of bulls and steers in Qinchuan cattle using next-generation sequencing, which showed that 580 upregulated genes were mainly associated with lipid metabolism, lipogenesis, and fatty acid transportation signaling pathways, and 52 DERs were identified [11].
At present, there are article on the miRNA-mRNA regulatory network fine-tuning the porcine muscle fiber, but no artical on the conjoint analysis of functional genes and miRNAs related to milk fat metabolism in BMECs. Therefore, we extracted total RNA from BMECs with high and low milk fat percentages for a conjoint analysis of RNA-seq and Solexa sequencing data to screen candidate DERs and their target differentially expressed genes (DEGs) related to milk fat percentage and verify the regulatory mechanism of DERs on milk lipid metabolism to provide a molecular basis for future gene and miRNA studies of milk fat percentage and the marker-assisted selection of milk fat percentage trait.

Functional DEGs on milk fat percentage screened by RNA-seq analysis
Among all expressed genes, a total of 17,771 genes were identified as differentially Functional enrichments of the DEGs by Gene Ontology (GO) analysis showed that 4132 GO terms were enriched, of which 2725, 397 and 1010 were enriched in biological process, cellular component and molecular function GO terms, respectively. A total of 531 biological process GO terms were significantly enriched (p<0.05), 54 cellular component GO terms were significantly enriched (p<0.05), and 134 molecular function GO terms were significantly enriched (p<0.05). The most significant GO terms are listed in ascending order of the p value ( Fig 1A). To verify the relationship and connection of significant GO terms, a GO-Tree (Fig 1B) was constructed by exploring the regulatory networks of significant biological process GO terms (p<0.01), which suggested that most GO terms were upregulated, and the GO terms of the upregulated genes were related to the unsaturated fatty acid biosynthetic process and the fatty acid biosynthetic process, which are related to the regulation of lipid metabolism.
Functional enrichments on the DEGs by Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that 228 KEGG pathways were enriched, of which 45 KEGG pathways were significantly enriched (p<0.05). The most significant KEGG pathways are listed in ascending order of the p value (Fig 2A). A Pathway-Act-Network (Fig 2B) was constructed by analyzing the regulatory network of significant pathways of both the up-and downregulated genes (p<0.05), which suggested that most KEGG pathways were downregulated, and the downregulated genes are related to fat digestion and absorption and the P-53 signaling pathway.
To systematically study the regulatory relationship between different genes, all DEGs were used to construct an intergene signal transduction network, and the intergene interaction network was integrated to obtain the interaction network between genes (Fig 3), such as secreted frizzled related protein 1 (SFRP1) gene inhibits the expression of Wnt family member 5B (WNT5B), and SH3 domain containing GRB2 like 3, endophilin A3 (SH3GL3) gene binds fibroblast growth factor receptor 3 (FGFR3) gene.

Functional DERs on milk fat percentage screened by Solexa sequencing
Among all expressed miRNAs, a total of 347 miRNAs were identified as differentially  (Table 3).
All the target genes regulated by the DERs were predicted by the combination of the miRbase and the TargetScan. Functional enrichments by GO analysis showed that 2240 GO terms were enriched, of which 1513, 213 and 514 were enriched in biological process, cellular component and molecular function GO terms, respectively. A total of 386 biological process GO terms were significantly enriched (p<0.05), 45 cellular component GO terms were significantly enriched (p<0.05), and 120 molecular function GO terms were significantly enriched (p<0.05). The most significant GO terms are listed in ascending order of the p value ( Fig 4A). Functional enrichments by KEGG analysis showed that 152 KEGG pathways were enriched, of which 21 KEGG pathways were significantly enriched (p<0.05). The most significant KEGG pathways are listed in ascending order of the p value ( Fig 4B).

Analysis
After the conjoint analysis of RNA-seq and Solexa sequencing data in BMECs with high and

The detection of triglyceride production regulated by miRNAs in BMECs
To confirm the effect of miRNAs on triglyceride content, we collected BMECs 48 h posttransfection with the mimics, inhibitor and negative control plasmids of 6 DERs. Compared with the negative control (Nc): triglyceride content was significantly increased in the miR-21-3p mimics group but significantly decreased in the miR-21-3p inhibitor group (p<0.05); there was no significant change in the miR-124a mimics group but a significant increase in the miR-124a inhibitor group (p<0.05); triglyceride content was significantly increased in the miR-148a mimics group but significantly decreased in the miR-148a inhibitor group (p<0.05); there was no significant change in the miR-877 mimics group but a significant increase in the miR-877 inhibitor group (p<0.05); triglyceride content was significantly decreased in the miR-2382-5p mimics group (p<0.05), but there was no significant change in the miR-2382-5p inhibitor group; and triglyceride content was significantly decreased in the miR-2425-5p mimics group but significantly increased in the miR-2425-5p inhibitor group (p<0.05) (Fig 7), suggesting that bta-miR-21-3p and bta-miR-148a promote triglyceride production, whereas bta-miR-124a, bta-miR-877, bta-miR-2382-5p and bta-miR-2425-5p inhibit triglyceride production.

Quantitative RT-PCR analysis of the DERs and their target DEGs
Compared with the negative control: the expression of the 6 DERs was significantly increased in the mimics group and significantly decreased in the inhibitor group (p<0.05); the expression of glutamic-oxaloacetic transaminase 2 (GOT2) and SH3 domain binding protein 5 (SH3BP5) was upregulated in both the miR-21-3p mimics and miR-21-3p inhibitor groups (p<0.05); the expression of peroxiredoxin 6 (PRDX6) was significantly upregulated in the miR-124a inhibitor group (p<0.05); the expression of immunoglobulin like domain containing receptor 2 (ILDR2) was upregulated in both the miR-124a mimics and miR-124a inhibitor groups; the expression of phosphodiesterase 4D (PDE4D) was significantly upregulated in the miR-148a inhibitor group but downregulated in the miR-148a mimics group (p<0.05); the expression of paternally expressed 10 (PEG10) was significantly downregulated in the miR-148a mimics group but upregulated in the miR-148a inhibitor group (p<0.05); the expression of superoxide dismutase 3 (SOD3) was significantly downregulated in the miR-2382-5p mimics group but upregulated in the miR-2382-5p inhibitor group (p<0.05); and the expression of ADAM metallopeptidase with thrombospondin type 1 motif 1 (ADAMTS1) was significantly downregulated in the miR-2425-5p mimics group but upregulated in the miR-2425-5p inhibitor group (p<0.05) ( Fig   8). The quantitative RT-PCR results showed that PRDX6, PDE4D, PEG10, SOD3 and ADAMTS1 were negatively regulated by DERs.

Target verification of miRNAs and target genes by luciferase reporter assay
Compared with the DERs mimics+pmiR-BR-REPORT (si) groups: luciferase activity was increased in both the miR-21-3p mimics+GOT2-wild type (WT) and miR-21-3p mimics+GOT2-mutant type (mut) groups; luciferase activity was increased in both the miR-21-3p mimics+SH3BP5-WT and miR-21-3p mimics+SH3BP5-mut groups; luciferase activity was increased in both the miR-124a mimics+PRDX6-WT and miR-124a mimics+PRDX6-mut groups; luciferase activity showed no difference between the miR-124a mimics+ILDR2-WT and miR-124a mimics+ILDR2-mut groups; luciferase activity was significantly decreased in the miR-148a mimics+PDE4D-mut group (p<0.05); luciferase activity was significantly decreased in the miR-148a mimics+PEG10-mut group (p<0.05); luciferase activity was significantly decreased in the miR-2382-5p mimics+SOD3-mut group (p<0.05); and luciferase activity was significantly decreased in the miR-2425-5p mimics+ADAMTS1-mut group (p<0.05) (Fig 9). The results revealed that no targeted binding sequence exists in the 3'-UTR of genes which we predicted to be located between GOT2, SH3BP5 and bta-miR-21-3p, PRDX6, ILDR2 and bta-miR-124a. However, the targeted binding sequence exists in the 3'-UTR of genes which we predicted to be located between PDE4D and bta-miR-148a, PEG10 and bta-miR-877, SOD3 and bta-miR-2382-5p, and ADAMTS1 and bta-miR-2425-5p. By analyzing the regulatory network of target genes and miRNAs, we found that the mechanism of miRNAs on milk fat percentage is very complex. A single miRNA regulates a single gene [12], a single miRNA regulates multiple target genes at the same time [13], or a single target gene is regulated by multiple miRNAs simultaneously [14,15], thus affecting the lipid metabolism of BMECs. In our experiments, 6 miRNAs (miR-21-3p, miR-124a, miR-148a, miR-877, miR-2382-5p, and miR-2425-5p) and 8 predicted target genes (GOT2, SH3BP5, PRDX6, ILDR2, PDE4D, PEG10, SOD3, and ADAMTS1) were selected to verify the target relationship. The following may explain why some did not have a significant regulatory relationship. The DEGs were screened according to the targeted regulatory relationships predicted by TargetScan and miRBase in the conjoint analysis of RNA-seq and Solexa sequencing data; however, the regulation of target genes by miRNAs in cells is complex. In our luciferase reporter assay, we inserted only approximately 60 bp of the DNA sequence before and after the predicted target DNA sites rather than the entire 3'-UTR of the gene. Therefore, there may be a site effect between the binding of miRNAs to the predicted target DNA sequence. We will verify this speculation in a future study.

Discussion
Triglyceride is a fatty molecule formed by the polymerization of three molecules of longchain fatty acids and one molecule of glycerol that can be transported in the blood in the form of lipoproteins [16]. The metabolic regulation of triglyceride in mammals plays an important role in the ability and balance of the mammalian body. Triglyceride has the highest lipid content in the human body. Excessive triglyceride can lead to obesity, cardiovascular, fatty liver, hyperlipidemia and other diseases, and the abnormal state of triglyceride synthesis can also lead to lipid metabolism disorders and other diseases [17,18]. Triglyceride is also an important part of milk, accounting for 90% of milk fat, and the synthesis of triglyceride in cattle mammary cells can affect milk fat percentage [19,20].
MiR-148a and miR-17-5p can promote the synthesis of triglyceride and fatty acid metabolism in goat mammary epithelial cells by negatively regulating peroxisome proliferator activated receptor alpha (PPARa) and PPARG coactivator 1 alpha (PPARGC1a), and interference with the expression of PPARa and PPARGC1a can promote the synthesis of milk fat by regulating the genes related to lipid metabolism [21]. In our study, conjoint analysis showed that bta-miR-148a was significantly downregulated in BMECs with high milk fat percentage, but miR-148a can significantly increase the content of triglyceride (p < 0.05). Studies have shown that miR-124a can inhibit the expression of patatin like phospholipase domain containing 2 (ATGL) and abhydrolase domain containing 5 (CGI-58) to reduce fat decomposition, leading to an increase in triglyceride content and the accumulation of lipids [22], which is contrary to our conjoint analysis data showing that bta-miR-124a was significantly downregulated in BMECs with high milk fat percentage and triglyceride content was significantly increased in the miR-124a inhibitor group (p < 0.05).
Previous research has shown that miR-21-3p plays an important role in triglyceride production by regulating the expression of the target gene ELOVL fatty acid elongase 5 (Elovl5) [23], miR-29b can promote triglyceride production and suppress apoptosis in BMECs by regulating the expression of the target genes lipoprotein lipase (LPL) and thymine DNA glycosylase (TDG), which are consistent with our conjoint analysis data showing that bta-miR-21-3p and bta-miR-29b are significantly upregulated in BMECs with high milk fat percentage. In summary, a large number of known and novel DERs were screened by conjoint analysis are mainly related to triglyceride production and milk fat metabolism, suggesting that DERs and target DEGs can assist in the candidate markerassisted selection of milk fat percentage.

Conclusion
Together, in the present study provides a scientific and effective research method screening functional miRNA and gene related to milk fat percentage by conjoint analysis of RNA-seq and Solexa sequencing of BMECs with high and low milk fat percentages. We identified the regulation of 6 DERs on lipid metabolism, which showed that bta-miR-21-3p and bta-miR-148a promote triglyceride synthesis, whereas bta-miR-124a, bta-miR-877, bta-miR-2382-5p and bta-miR-2425-5p inhibit triglyceride synthesis. Furthermore, the quantitative RT-PCR results showed that the expression level of PRDX6, PDE4D, PEG10, SOD3 and ADAMTS1 were negatively regulated by DERs, and luciferase reporter assay showed that the target relationships existed between PDE4D and bta-miR-148a, PEG10 and bta-miR-877, SOD3 and bta-miR-2382-5p, and ADAMTS1 and bta-miR-2425-5p, which means miRNA affect the expression of gene via not only target relationship but also other regulatory factors. Our findings identified miRNA-mRNA regulatory network involved in lipid metabolism in bovine mammary epithelial cells and provided a foundation for the study of the molecular mechanisms related to milk fat percentage.

Total RNA extraction
Total RNA from BMECs with high and low milk fat percentages was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. Total RNA was further treated with RNase-free DNase I (NEB, Beijing, China), and RNA longer than 18 nucleotides was recovered. RNA concentration was measured by NanoDrop 2000 (Thermo Scientific, USA). The quality of RNA was assessed by the Agilent2200 (Agilent, USA).

Conjoint analysis of RNA-seq and Solexa sequencing
(1) RNA-seq analysis: The sequencing library of each RNA sample was prepared using Ion Total RNA-Seq Kit v2 according to the protocol provided by manufacturer (Life technologies, USA). The transcript with poly (A)+ RNA of Cattle were analyzed using highthroughput Life technologies Ion Proton Sequencer. Reads sequenced were filtered and mapped to Cattle genome (download from NCBI) using Mapsplice software. The mapped reads was counted to achieve the expression of each gene based on the gene annotation information from NCBI database. In order to achieve the best quality of the RNA Sequencing, Novelbio applied the reads filtration to filter the reads with lower quality and short sequence under following criteria: read length > 50; over 30% base quality >13.
Using the internationally recognized algorithm EBSeq based on the fold change (FC) and false discovery rate (FDR) thresholds to screen the DEGs in the BMECs of Chinese Holstein cows with high and low milk fat percentages [25,26]. Genes were considered differentially expressed according to log 2 FC >1 or log 2 FC <-1and FDR <0.05.
(2) Solexa sequencing analysis: We performed a differential expression analysis according to the data of Solexa sequencing obtained by a previous study in our laboratory and screened the DERs in the BMECs of Chinese Holstein cows with high and low milk fat percentages using the internationally recognized EBSeq algorithm based on the FC and FDR thresholds [27]. MiRNAs were considered differentially expressed according to log 2 FC >1 or log 2 FC <-1and FDR <0.05.
(3) Conjoint analysis of RNA-seq and Solexa: All target genes regulated by DERs were predicted by the combination of miRbase (http://www.mirbase.org) and TargetScan (http://www.targetscan.org). All the predicted target genes were intersected with the DEGs obtained by RNA-seq, and then the DEGs and DERs were obtained by a combined analysis.

GO and KEGG pathway enrichment analyses
GO analysis was applied to analyze the main functions of the DEGs according to GO, which is the functional classification of NCBI (http://www.ncbi.nlm.nih.gov). The GO category was classified by Fisher's exact test and the χ 2 test, and the FDR was calculated to correct the P value: the smaller the FDR was, the smaller the error in judging the P value [28]. We computed the P values for the GO terms of all the DEGs. Enrichment provides a measure of the significance of the function: as the enrichment increases, the corresponding function is more specific, which helps us identify those GO terms with concrete functional descriptions in the experiment.
Pathway analysis was used to determine the significant pathways of the DEGs. Pathway annotations of microarray genes were downloaded from KEGG (http://www.genome.jp/Wegg/). Fisher's exact test was used to find the significant enrichment pathways. The resulting P values were adjusted using the BH FDR algorithm (Benjamini and Hochberg 1995). Pathway categories with an FDR < 0.05 were reported.
(2)Construction of luciferase reporter vectors: We predicted the DNA-binding sites of miRNAs to the target gene's 3'-UTR with TargetScan. For the luciferase reporter assay, the WT and mut DNA sequences were obtained by annealing performed in a 10 μL reaction with 0.5 μL forward primer, 0.5 μL reverse primer, 1 μL 10 × Ex Taq buffer and 8 μL ddH 2 O under the following procedure: 95°C for 5 minutes, then dropped at a rate of 1.5°C per minute to 22°C. The DNA fragment was then ligated with the linearized pmiR-BR-REPORT plasmid to construct the WT and mut vectors. All DNA oligo nucleotides used for annealing were synthesized by Sangon (Shanghai, China), and the sequence details are shown in Table 1.

Determination of triglyceride content
The Triglyceride Extraction Kit (APPLYGEN, China) was used to separate the cells following the manufacturer's instructions. Then, the cell lysates were collected in 1.5 mL centrifuge tubes, heated at 70°C for 10 minutes, centrifuged at 2000 g for 5 minutes at room temperature, and the supernatant was extracted to detect the contents of triglyceride and protein on a SpectraMax M5 Microplate Reader. The protein content was used as a reference to correct the variation in cell number, each sample was analyzed in triplicate.

Luciferase reporter assay
The luciferase reporter experiments were performed with BMECs at a density of 1 × 10 5 cells in 24-well plates (BIOFIL, Guangzhou, China). To validate the target relationship of the candidate miRNAs and genes, 1.5 µg miRNA mimic plasmid was cotransfected with 1.5 µg WT, 1.5 µg mut or 1.5 µg pmiR-BR-REPORT plasmid into BMECs. Forty-eight hours after transfection, the firefly (hluc+) and Renilla (hRluc) luciferase activities were measured with the Dual-Glo luciferase assay system (Promega) according to the manufacturer's instructions. The hluc+ gene was used as a reference gene to correct the variation in transfection efficiency, and the relative luciferase activity was calculated as hRluc/hluc+.
All experiments for the luciferase reporter assay were performed in triplicate.

Statistical analysis
The data were analyzed by GraphPad 6.0 statistical software with a completely random design of Student's t test for the comparison of two independent groups. The mean and standard deviation (SD) or standard error of mean (SEM) among triplicate samples were calculated, and p<0.05 was considered to indicate statistical significance.

Availability of data and materials
The data sets used and analyzed during the current study are available. The RNA-seq and Solexa sequencing data have been submitted to the GenBank databases under accession number GSE137488.

Competing interests
The authors declare that they have no competing interests        The detection of triglyceride production regulated by DERs in BMECs A Relative