Xiaogang Cui
Shengli Zhang
Qin Zhang
Xiangyu Guo3
Changxin Wu
Dongxiao Sun- 1Key Laboratory of Animal Genetics, Breeding and Reproduction of Ministry of Agriculture and Rural Affairs, National Engineering Laboratory of Animal Breeding, College of Animal Science and Technology, China Agricultural University, Beijing, China
- 2Key Lab of Medical Molecular Cell Biology of Shanxi Province, Institutes of Biomedical Sciences, Shanxi University, Taiyuan, China
- 3Center for Quantitative Genetics and Genomics, Aarhus University, Tjele, Denmark
A total of 31 differentially expressed genes in the mammary glands were identified in our previous study using RNA sequencing (RNA-Seq), for lactating cows with extremely high and low milk protein and fat percentages. To determine the regulation of milk composition traits, we herein investigated the expression profiles of microRNA (miRNA) using small RNA sequencing based on the same samples as in the previous RNA-Seq experiment. A total of 497 known miRNAs (miRBase, release 22.1) and 49 novel miRNAs among the reads were identified. Among these miRNAs, 71 were found differentially expressed between the high and low groups (p < 0.05, q < 0.05). Furthermore, 21 of the differentially expressed genes reported in our previous RNA-Seq study were predicted as target genes for some of the 71 miRNAs. Gene ontology and KEGG pathway analyses showed that these targets were enriched for functions such as metabolism of protein and fat, and development of mammary gland, which indicating the critical role of these miRNAs in regulating the formation of milk protein and fat. With dual luciferase report assay, we further validated the regulatory role of 7 differentially expressed miRNAs through interaction with the specific sequences in 3′UTR of the targets. In conclusion, the current study investigated the complexity of the mammary gland transcriptome in dairy cattle using small RNA-seq. Comprehensive analysis of differential miRNAs expression and the data from previous study RNA-seq provided the opportunity to identify the key candidate genes for milk composition traits.
Background
MicroRNAs (miRNAs), which are a class of non-coding small RNA (sRNA) molecules with the length of 18-24 nucleotides, are important regulators of gene expression. They can play important roles in a wide range of biological processes, including animal and plant development, cell differentiation, proliferation, apoptosis, and metabolism (Martello et al., 2010; Chen et al., 2012; Rottiers and Näär, 2012; Almughlliq et al., 2019; Barbu et al., 2020). In animal cells, miRNAs interact with a specific sequence in mRNA of the target gene and post-transcriptionally negatively regulate the expression of target genes by inhibiting their translation or inducing degradation of the target mRNAs (Huntzinger and Izaurralde, 2011; Barbu et al., 2020). MiRNAs have emerged as new potential biomarkers for miRNA-gene interactions and gene networks responsible for human diseases and economically important traits in livestock. Several diseases and conditions have been reported to be linked with the abnormal expression in miRNAs relating with differentiation, apoptosis and development (Lewis et al., 2005; Berezikov et al., 2006b; Lee et al., 2007). Many experimental techniques and computational methods have been developed to identify miRNAs (Aravin and Tuschl, 2005; Berezikov et al., 2006a; Landgraf et al., 2007), and large number of miRNAs have been identified in primates, rodents, birds, fish, and plants (Lagos-Quintana et al., 2003; Chen et al., 2005; Finucane et al., 2008; Glazov et al., 2008).
The bovine mammary gland is a complex organ which grows and develops after calving and is able to produce more than 30,000 kg of milk in a complete lactation cycle (Hennighausen and Robinson, 2005; Muroya et al., 2019). Because of its important functions, the mammary gland, especially mammary epithelial cells, has been used as the target tissue for gene expression profiling in order to identify key genes underlying milk production traits in dairy cattle (Silveri et al., 2006; Bionaz and Loor, 2008, 2011; Bionaz et al., 2012; Zhang et al., 2016; Pu et al., 2017; Cai et al., 2018; Ju et al., 2018; Yang et al., 2018; Billa et al., 2019; Li et al., 2020). However, only a few studies have been reported related to the miRNAs in the bovine mammary gland. A total of 798 mature bovine miRNAs have been deposited in miRBase (Luoreng et al., 2018), Release 22.1 (October 2018) and 55 of them were detected in the mammary gland. Li et al. (2012b) reported 283 known miRNAs and 74 novel miRNAs in the mammary gland of Holstein cows, among which 56 miRNAs were differentially expressed between lactating and non-lactating cows and might be involved in regulating lactation. Shen et al. (2016) identified 292 known miRNAs and 116 novel miRNAs in the bovine mammary epithelial cells, and three of them (bta-miR-33a, bta-miR-152 and bta-miR-224) might be involved in milk fat metabolism. Li et al. (2015) detected 370 known and 341 novel miRNAs in the bovine mammary gland infected with Staphylococcus aureus, and 358 known and 232 novel miRNAs in control group, 77 of which were differentially expressed between infected and healthy Holstein cows. In addition, Le Guillou et al. (2012) found that the overexpression of miR-30b caused a defect in lactation and delayed involution in mouse mammary gland.
In a previous study from our lab (Cui et al., 2014), 31 differentially expressed genes was identified by using RNA sequencing (RNA-Seq) to investigate the mammary gland epithelial tissues of four lactating Holstein cows with extremely high and low milk protein (PP) and fat percentages (FP). The objectives of the present study were to investigate the miRNA expression profiles in the same mammary gland samples that were used in the previous RNA-Seq study to identify known and novel miRNAs, and to perform an analysis of the differentially expressed miRNAs and previously identified genes. Some candidate miRNAs and their target genes that may be involved in milk protein and fat metabolism were identified.
Materials and Methods
Animals and Mammary Gland Tissue Samples
In the current study, the mammary gland epithelium samples of four lactating Chinese Holstein cows (high group vs. low group) same as our previous RNA-Seq experiment (Cui et al., 2014) were used. These four cows were selected from 30,000 Holstein cows in Beijing Sanyuanlvhe Dairy Farming Center, and the average PP and FP were 3.1% (2.7–3.8%) and 3.6% (3.1–4.5%) in this population. In order to keep the environmental factors identical, these four cows in almost the same period of lactation (353, 341, 377, and 325 days) were collected from the same farm possessing a total of 800 Holstein cows. Selected cows were divided into two groups according to the phenotypic values for PP and FP: two cows (high group) had high PP (3.6% and 3.8%) and FP (3.9% and 4.5%); the other two cows (low group) showed low PP (3.0%, 2.9%) and FP (3.2%, 3.1%).
The cows were killed by electroshock, and then they were bled, skinned, and dismembered in the same slaughterhouse. The rear mammary gland from each individual was harvested within 30 min after slaughtered. White mammary ducts and pink epithelium tissue were clearly observed when the right rear quarter of the mammary gland was cut in half lengthways from the teat and some milk were flowed out. Five pieces of epithelium tissue samples per cow were carefully collected and placed into a clean RNAse-free Eppendorf tube, and then stored in liquid nitrogen for subsequent RNA isolation. All procedures of collecting samples were carried out in strict accordance with the protocol approved by the Animal Welfare Committee of China Agricultural University (Permit Number: DK996). Total RNA was extracted from one piece of mammary gland epithelium samples from each cow and quality was controlled according to the protocols described by Cui et al. (2014). The value of RNA integrity number (RIN) from each sample was above 8.0.
Small RNA Sample Preparation and Sequencing
The preparation of small RNA library, including quality control and sequencing, was performed by Novogene (Beijing, China). The preparation of library was performed on 3 μg total RNA per sample using an IlluminaTruSeq™Small RNA Sample Preparation Kit (Illumina, San Diego, CA, United States). The samples were indexed using four codes in order to facilitate sequencing of these samples on one flow cell channel. Quality control in library preparation showed that adapter-adapter contamination was <5% and 85% of the sequences were miRNAs. The samples were subsequently sequenced on the Illumina Hiseq2000 platform and 50-bp single-end reads were obtained.
Sequencing Data Analysis
The sequencing data were obtained in the format of Illumina FASTQ (Illumina). The procedure of data filtering included removing low quality reads, reads containing poly-N stretches, reads with 5′primer contaminants, reads with 3′primers or the insert tag, and reads with poly-A, T, G, or C stretches. Thereafter, the sRNA tags within a certain range (18-30 nt) were retained for the successive steps. The Q20, Q30, and GC-content of the cleaned reads were calculated to evaluate the quality of data. Then, the sRNA tags were mapped to the bovine genome assembly (UMD3.1.66) using Bowtie (Langmead et al., 2009), no mismatches were allowed and the “seed” region size was set at 8 (Gupta et al., 2012; Giurato et al., 2013; Aggarwal et al., 2014; Kuksa et al., 2018). The mapped sRNA tags were aligned to the 798 bovine miRNA precursor sequences in miRBase (Release 22.1) to identify the known miRNA in the sRNA libraries allowing one mismatch. The sRNA tags that matched known miRNAs from species other than bovine may be novel bovine miRNAs, and were predicted the secondary structure, the Dicer cleavage site, and the minimum free energy of the mapped sRNA sequences using the miREvo (Wen et al., 2012) and miRDeep2 (Friedländer et al., 2012) software packages.
The expression of miRNA was measured as counts per million (CPM) using the following formula: normalized expression = mapped read count/total reads × 1000000 (Zhou et al., 2010), and DESeq2 R package (1.8.3) (Anders and Huber, 2010; Trapnell et al., 2013) was used to identify significantly differentially expressed miRNAs between high and low groups of cows. The threshold for differential expression was—log2 (FC)— > 1 and FDR p < 0.05 when using DESeq2 R package for differential expression miRNA analysis so that miRNAs with—log2 (FC)— > 1 and adjusted FDR p < 0.05 were designated as differentially expressed.
Furthermore, two cows in the same group were used to eliminate the background noise of individual-specific transcription by applying a pairwise approach, which enabled acquisition of more relevant data from the two groups.
Target Prediction, Pathway, and Annotation Analysis
TargetScan 6.2 and MiRanda (Enright et al., 2003) were used to predict putative target genes with the established miRNA seed database and the bovine genome sequence (UMD3.1.66). TargetScan 6.2 predicts targets by searching for the presence of conserved 8mer, 7mer, and 6mer sites that match the seed region of each miRNA. MiRanda predicts targets based on a development of the miRanda algorithm which incorporates current biological knowledge on target rules and on the use of an up-to-date compendium of mammalian miRNAs.
Gene ontology (GO) functional enrichment analysis was used for the candidate target genes of the miRNAs. GOseq with the Wallenius non-central hyper-geometric distribution (Young et al., 2010), which can adjust for the bias in gene length, was implemented for the GO enrichment analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) pathways analysis was performed using KOBAS 2.0 (Mao et al., 2005) software to test the statistical enrichment of the candidate target genes in the KEGG pathways.
Quantitative Real Time PCR
Expression levels of selected miRNAs were confirmed by quantitative real-time PCR (qRT-PCR) using the DyNAmo SYBR Green PCR kit (Applied Biosystems, Foster City, CA, United States) on a LightCycler480 (Roche Applied Science, Penzberg, Germany). qRT-PCR of target mRNAs was performed using specific miRNA stem-loop primers (Supplementary Table 8) and all reactions were run in triplicate. Relative quantification of miRNA was quantified using the 2–Δ Δ CT method and normalized against the U6 gene (ssD0904071006: Guangzhou RiboBio, Guangzhou, China) for each sample.
Plasmid Construction and Site-Directed Mutagenesis of 3′UTR in Predicted Target Genes
The 3′un-translated region (UTR) of four predicted target genes for the identified miRNAs, TRIB3, M-SAA3.2, PTHLH, and VEGFA, were PCR amplified using DNA collected from the bovine mammary gland samples applied for sequencing in this study as a template, and connected into pmirGLO Dual-Luciferase miRNA Target Expression Vector (pmirGLO, Promega) (Figure 1), respectively. The primers were listed in the Table 1. Afterward, the connected products were transfected into Escherichia coli, and then verified the correct sequence and orientation by sequencing. The QuikChange site-directed mutagenesis kit (Stratagene, La Jolla, CA, United States) was used to generate the 3′UTR variants of TRIB3, M-SAA3.2, PTHLH, and VEGFA where seed sequences recognized by microRNAs were deleted (Figures 2, 3). After the point mutation, same way was applied in order to find the correct mutant sequences for such four genes.
Figure 1. The pmirGLO vectors with the predicted 3’UTR target sequences of the 4 differentially expressed genes (A) pmirGLO-TRIB3-3′UTR; (B) pmirGLO- M-SAA3.2-3′UTR; (C) pmirGLO-PTHLH-3′UTR; (D) pmirGLO-VEGFA-3′UTR.
Table 1. PCR primers for TRIB3, PTHLH, VEGFA and M-SAA3.2.
Figure 2. Domain structures of the 4 differentially expressed genes showing the locations of the seed sequence of the miRNAs within the 3’UTR of theirs.
Figure 3. Locations and sequences of the miRNAs target sites in the 3’UTR of the 4 differentially expressed genes. The sequences of the miRNAs are indicated, along with mutations introduced in the target sites (underlined nucleotides) for generating the mutated reporter constructs.
Luciferase Reporter Assays
To further explore the repressing mechanism of miRNAs on the expression of 4 target genes (TRIB3, M-SAA3.2, PTHLH, and VEGFA) expression, the full-length TRIB3, M-SAA3.2, PTHLH and VEGFA 3′UTRs and the corresponding mutant version (the seed sequences were deleted) were transfected into human embryonic kidney HEK293 cells (GM-070001H: Shanghai, China), respectively. These cells were cultured at 37°C with 5% CO2 in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 4.5 g/liter glucose, 5% fetal bovine serum (Invitrogen), 2 mmol/liter glutamine, and antibiotics. Before transfection, HEK293 cells were plated into 24-well plates at 1.0 × 105 cells/well 24 h. 30 ng empty pmirGLO vector, pmirGLO-TRIB3/M-SAA3.2/PTHLH/VEGFA-3′UTR with 50 μl opti-MEM (Invitrogen) and 30 nM (final concentration) mimic miRNA, inhibitor miRNA, control miRNA (GenePharma) were co-transfected into each well with 1 μl Lipofectamine 2000 (Invitrogen). 30 ng mutants of the TRIB3/M-SAA3.2/PTHLH/VEGFA 3′UTR with 50 μl opti-MEM (Invitrogen) and 30 nM (final concentration) mimic miRNA, control miRNA (GenePharma) were co-transfected into each well with 1 μl Lipofectamine 2000 (Invitrogen). Relative firefly luciferase activities (normalized to Renilla luciferase activities) were measured 24 h after transfection with the Dual-Luciferase Reporter Assay Kit (Promega) on TECAN Infinite 200 multifunctional microplate reader (TECAN). All experiments were performed in triplicate so that data averaged from three independent experiments.
Results
Sequencing and Mapping of the sRNA Tags
Four new miRNA libraries were constructed using sRNA isolated from bovine mammary glands and sequenced using Illumina next-generation sequencing. A total of 10,538,878 (high milk PP and FP), 12,745,512 (high milk PP and FP), 9,744,027 (low milk PP and FP), and 9,682,136 (low milk PP and FP) high-quality cleaned reads were obtained from the four sRNA libraries (Supplementary Table 1; NCBI SRA accession numbers: SRR3631014, SRR3631016, SRR3631053, and SRR3631054). Distribution of the length for reads showed that most of the generated reads had 21 (>24%), 22 (>30%), and 23 (>13%) nucleotides (Supplementary Figure 1), which is the size of most known mature miRNAs. When aligning the sequenced reads against the bovine genome assembly (UMD3.1.66), it was found that 77.57%, 76.93%, 80.88%, and 78.15% of them uniquely aligned from the four libraries, respectively (Supplementary Table 1); 55-57% of them were aligned in the same direction as the reference genome sequence, and 20-25% were aligned in the opposite direction (Supplementary Table 2). The correlation coefficient (R2) between the two individuals within the high and low groups for milk PP and FP was calculated based on the CPM mapped fragment of each cow and was shown to be 0.988 and 0.980, respectively. This indicated that the similarity of the two biological replicates within each group was sufficiently high (Supplementary Figure 2).
MicroRNAs Identification and Target Prediction
Among the uniquely aligned reads across the four samples and six downloaded miRNA libraries (Cai et al., 2018), 24,320,809 (54.4%) matched known miRNAs in miRBase (Release21.0), which resulted in 497 known bovine miRNAs and 49 novel bovine miRNAs were identified (Supplementary Tables 3, 4). Subsequently, two well-established target prediction tools, TargetScan and miRanda, were used to predict target mRNAs of the miRNAs, and a total of 12,202 target genes were commonly predicted for the known and novel miRNAs (Supplementary Table 5). It is noteworthy that some well-known genes associated with milk composition traits were included such as β-casein (CSN2), κ-casein (CSN3), α-lactalbumin (LALBA), diacylglycerol O-acyltransferase 2 (DGAT2), growth hormone receptor (GHR), signal transducer and activator of transcription 5B (STAT5B), and stearoyl-coenzyme A desaturase (SCD) etc. This finding implied that the identified mammary miRNAs in this study were involved in metabolism of milk protein and lipid through the regulation of key genes affecting these traits.
Differentially Expressed miRNAs Between the High and Low Groups for Milk PP and FP and Target Prediction
The miRNAs that differed between the high and low PP and FP groups were determined in this study. A total of 71 top half miRNAs displayed significantly differential expression between the high and low groups using the DEseq2 algorithm (p < 0.05, FDR q < 0.05), with 35 were up-regulated and 36 were down-regulated in the high milk PP and FP group compared with the low group (Table 2). Subsequently, a total of 5,634 target genes were commonly obtained for these differentially expressed miRNAs by TargetScan and miRanda (Supplementary Table 7).
Table 2. Seventy-one differentially expressed miRNAs between the high and low milk protein and fat percentages groups.
Afterward, the results of the sequencing were validated with an independent method of real-time PCR assay. By using the same four mammary gland samples as used for sequencing, eight known miRNAs and seven novel miRNAs identified in the present study were randomly chosen for validation. It was found that the expression levels of miR-125a, miR-2904, miR-345-5p, miR-378c and Novel-18 were significantly higher in the high milk PP and FP group than in the low group (p < 0.05), and the expression levels of miR-21-3p, miR-29c, miR-106b and miR-190a were lower in the high group than in the low group (p < 0.05). Whereas, Novel-13, Novel-2, Novel-22, Novel-32, Novel-4 and Novel-42 did not display significant differences on miRNA levels between the two groups (p > 0.05) (Figure 4). Such expression patterns were exactly consistent with those shown by small-RNA sequencing data.
Figure 4. mRNA expression levels of the 15 randomly selected miRNAs validated with qRT-PCR. *indicates p < 0.05. Blue columns represent the relative miRNA expression levels by qRT-PCR normalized by U6 in the high group and red columns represent the relative miRNA expression levels by qRT-PCR normalized by U6 in the low group.
Gene Ontology Enrichment and Pathway Analysis
To further investigate the functional associations of the target genes, gene ontology (GO) annotation analysis was performed. It was found that these targets have a wide range of diverse functions, among which most were involved in protein and lipid metabolism, mammary gland development and differentiation, and immune functions (p < 0.01, FDR q < 0.01). Under the GO biological process category, the enriched terms related to lipid and protein metabolisms and cell growth were included such as protein binding, protein localization, protein transport, protein complex, regulation of protein metabolic process, lipid biosynthetic process, programmed cell death, protein targeting, lipid metabolic process, amino acid transport, regulation of protein kinase activity, and cellular response to mechanical stimulus (Supplementary Table 6).
A KEGG metabolic pathway analysis was also performed to identify functions that associate with the predicted target genes using KOBAS. Targets were enriched for functions such as mitogen-activated protein kinases (MAPK), adipocytokine, mammalian target of rapamycin (mTOR), glycosphingolipid biosynthesis, glycerophospholipid metabolism, hypoxia inducible factor-1 (HIF-1), and phosphatidylinositol 3 kinase-protein kinase B (PI3K-Akt) signaling pathways (Table 3).
Table 3. KEGG pathways assigned to the predicted target genes of the 497 known and 49 novel miRNAs identified in this study.
For the 71 top half differentially expressed miRNAs, 5,634 target genes were obtained and the targets were highly enriched in biological process consisting of synthesis and metabolism of protein and energy metabolism, as well as pathways mainly related to synthesis and metabolism of lipid and protein including glutathione metabolism, NF-kappa B signaling pathway, mTOR signaling pathway, fatty acid degradation, fatty acid metabolism and protein processing in endoplasmic reticulum (Table 4).
Table 4. KEGG pathways assigned to the predicted target genes of the 71 differentially expressed miRNAs identified in this study.
Comparison of the Target Genes of the Differentially Expressed miRNAs and the Differentially Expressed Genes Reported Previously
In our previous study (Cui et al., 2014), 21 of the target genes, which are listed in Table 5, were found to be differentially expressed between the high and low groups using the same four mammary gland samples in the current study. Among the 21 differentially expressed target genes, the expressions of only six down-regulated genes and one up-regulated gene matched the expression profiles of the differentially expressed miRNAs that targeted them. While 5 down-regulated genes were targeted by at least one up-regulated miRNA each, and 10 genes were targeted by both up-regulated and down-regulated miRNAs. Especially, 7 of the 21 differentially expressed target genes were the most promising candidate genes affecting milk protein and fat percentage identified by integrated analysis of differential gene expression, previously reported quantitative trait loci (QTLs) and genome-wide association studies (GWAS) (Cui et al., 2014), including tribbles homolog 3 (TRIB3), serum amyloid A1 (SAA1), serum amyloid A3 (SAA3), mammary serum amyloid A3 (M-SAA3.2), vascular endothelial growth factor A (VEGFA), parathyroid hormone-like hormone (PTHLH) and ribosomal protein L23A (RPL23A). In addition, KEGG pathway analysis using KOBAS, showed that two of the 21 target genes, DNA-damage-inducible transcript 3 (DDIT3) and nuclear receptor subfamily 4, group A, member 1 (NR4A1), were involved in the MAPK signaling pathway that plays critical role in protein synthesis and metabolism and fatty acid metabolism pathway (p < 0.05; Tables 3, 4), and 2 other genes, vascular endothelial growth factor A (VEGFA) and cyclin-dependent kinase inhibitor 1A (CDKN1A), were involved in the mTOR, HIF-1, PI3K-Akt, p53 and duct acid secretion signaling pathways, which are mostly related to synthesis and metabolism of protein and fat (p < 0.05; Tables 3, 4).
Table 5. Twenty-one differentially expressed target genes from our previous study for the 71 differentially expressed miRNAs.
MicroRNAs Repress the Expression of Target Genes Through the Binding of a Specific Target Sequence in Their mRNA 3′UTR
To study the regulatory functions of the identified miRNAs, four differentially expressed genes were chosen including TRIB3, M-SAA3.2, PTHLH and VEGFA, which having the expression pattern negatively correlated with their targeting miRNAs. Using the dual luciferase reporter assays, whether miR-2904, miR-339b/miR-146b/miR-339a, miR-29c/miR-106b/miR-190a, and miR-2904/miR-106b/miR-21-3p regulated the expression of the TRIB3, M-SAA3.2, PTHLH and VEGFA, respectively, were detected. Consequently, it was found that the luciferase level in HEK293 cells with mimics of miR-2904 decreased 40% relative to those with the empty vector, respectively (p < 0.05), while the inhibitor of miR-2904 yielded the same luciferase level as negative control (p > 0.05) (Figure 5). However, when the predicted binding sites of such miRNA seed sequences were mutated, luciferase activity was efficiently restored to the control levels (p > 0.05; Figure 5). Such results clearly indicated the notable regulatory role of the miR-2904 on the expression of TRIB3 by directly targeting its 3′UTR. Similarly, with regard to M-SAA3.2, it was also found that the overexpression of miR-146b, miR-339a and miR-339b decreased the luciferase levels in HEK293 cells by 80%, 72% and 74% after transfecting these mimics compared with the negative controls, respectively (p < 0.05), and the depressed expression of such miRNAs did not change the luciferase level in HEK293 cells transfected with their inhibitors (p > 0.05), respectively (Figure 6). When the mutant 3′UTR of M-SAA3.2 and mimics of the 3 miRNAs were co-transfected, the luciferase activity was same as the control level (p > 0.05; Figure 6). For PTHLH, the luciferase level in HEK293 cells transfected with the mimics of miR-29c, miR-106b and miR-190a was decreased by 37%, 49%, and 50% relative to the negative control, respectively (p < 0.05; Figure 7), however, the same level was kept by transfecting the inhibitors of such miRNAs (p > 0.05), respectively (Figure 7). When the mutant version of the PTHLH 3′UTR and mimics of miR-29c, miR-106b and miR-190a were co-transfected, respectively, the luciferase activities were same as the control levels (p > 0.05; Figure 7). Whereas, the expression of VEGFA was not affected by miR-2904, miR-106b and miR-21-3p (p > 0.05, Figure 8).
Figure 5. MicroRNAs represses the expression of TRIB3 via binding the 3’UTR target sequence. Luciferase activity in HEK293 cells co-transfected with miRNA mimic, miRNA inhibitor, miRNA control and empty vector for the TRIB3 3’UTR. Luciferase activity was assayed 24 h after transfection. All luciferase values were normalized to Renilla luciferase. Blue columns represent the luciferase activity co-transfected with miRNA mimic control; Red columns represent the luciferase activity co-transfected with miRNA inhibitor control; Green columns represent the luciferase activity co-transfected with miRNA mimic; Pink columns represent the luciferase activity co-transfected with miRNA inhibitor. (A) Represents the luciferase activity of TRIB3 after over- or down-expressed miR-2904 compared with controls. (B) Represents the luciferase activity of TRIB3 after transfecting mutant vector of miR-2904 compared with control. *Significant difference between the control and the treatment; **Very significant difference between the control and the treatment.
Figure 6. MicroRNAs represses the expression of M-SAA3.2 via binding the 3’UTR target sequence. Luciferase activity in HEK293 cells co-transfected with miRNA mimic, miRNA inhibitor, miRNA control and empty vector for the M-SAA3.2 3’UTR. Luciferase activity was assayed 24 h after transfection. All luciferase values were normalized to Renilla luciferase. The meanings of different colors are consistent with Figure 5. (A,C,E) Represents the luciferase activity of M-SAA3.2 after over- or down-expressed miR-146b, miR-339a and miR-339b compared with controls, respectively. (B,D,F) Represents the luciferase activity of M-SAA3.2 after transfecting mutant vector of miR-146b, miR-339a and miR-339b compared with control, respectively. **Very significant difference between the control and the treatment.
Figure 7. MicroRNAs represses the expression of PTHLH via binding the 3’UTR target sequence. Luciferase activity in HEK293 cells co-transfected with miRNA mimic, miRNA inhibitor, miRNA control and empty vector for the PTHLH 3’UTR. Luciferase activity was assayed 24 h after transfection. All luciferase values were normalized to Renilla luciferase. The meanings of different colors are consistent with Figure 5. (A,C,E) Represents the luciferase activity of PTHLH after over- or down-expressed miR-29c, miR-106b and miR-190a compared with controls, respectively. (B,D,F) Represents the luciferase activity of PTHLH after transfecting mutant vector of miR-29c, miR-106b and miR-190a compared with control, respectively. **Very significant difference between the control and the treatment.
Figure 8. MicroRNAs did not repress the expression of VEGFA via binding the 3’UTR target sequence. Luciferase activity in HEK293 cells co-transfected with miRNA mimic, miRNA inhibitor, miRNA control and empty vector for the VEGFA 3’UTR. Luciferase activity was assayed 24 h after transfection. All luciferase values were normalized to Renilla luciferase. The meanings of different colors are consistent with Figure 5. (A,B,C) Represents the luciferase activity of VEGFA after over- or down-expressed miR-2904, miR-106b and miR-21-3p compared with controls, respectively.
Discussion
The current study is the first comparative profiles of the mRNA and miRNA transcriptome in the mammary gland epithelium of dairy cows to the best of our knowledge. In this study, we generated an extensive miRNA expression profile of the mammary glands from lactating cows with extremely high and low milk PP and FP, and identified a total of 497 known bovine miRNAs and 49 novel bovine miRNAs. In previous studies, bovine miRNAs were identified using computational and direct cloning approaches (Coutinho et al., 2007; Gu et al., 2007; Jin et al., 2009, 2010; Long and Chen, 2009; Li et al., 2012b, 2015; Shen et al., 2016). Li et al. (2012b) identified 298 known miRNAs in lactating and non-lactating mammary gland of Holstein cows using miRNA-seq; 204 of them were among the 497 known miRNAs identified in the current study. Furthermore, 9 of the 71 differentially expressed miRNAs (miR-100, miR-10a, miR-133a, miR-1, miR-146b, miR-148a, miR-221, miR-30f, and miR-339b) identified in the current study were also reported by Li et al. (2012b) as differentially expressed between lactating and non-lactating bovine mammary glands. Gu et al. (2007) identified 31 distinct miRNAs in the mammary glands of Holstein cows, and all of these miRNAs was detected in the present study except miR-142b. Shen et al. (2016) identified 292 known miRNAs in the bovine primary mammary cells, among which 217 miRNAs and 38 differentially expressed miRNAs were also identified in the current study. For the 30 differentially expressed miRNAs in the lactating goat mammary gland fed ad libitum or deprived of food affecting milk composition reported by Mobuchon et al. (2015), only 6 miRNAs, including miR-660-5p, miR-451-5p, miR-125b, miR-196a, miR-223-3p, and miR-223-5p were detected as well in the current study. miR-30b related to lactation in mouse (Le Guillou et al., 2012) was also detected in this study, but did not show differential expression between high and low groups. The reason could be due to the mammary gland tissues were collected from different time points of lactation between the previous (Le Guillou et al., 2012) and the current studies.
miR-15a has been reported to be critical in cell development (Bonci et al., 2008), cell cycle (Bandi et al., 2009), and death (Cimmino et al., 2005; Aqeilan et al., 2010). Li et al. (2012a) found that miR-15a can inhibit the viability of mammary epithelial cells as well as the mRNA and protein expression of GHR, which is a major gene for milk composition traits (Bonci et al., 2008). In the current study, we also detected miR-15a and predicted that it may target GHR as well as candidate genes for milk PP and FP identified in our previous study, namely activating transcription factor 3 (ATF3), VEGFA, parathyroid hormone-like hormone (PTHLH), cation transport regulator homolog 1 (CHAC1), and NR4A1. Therefore, miR-15a was considered may affect milk composition by regulating the expression of these genes, although miR-15a was not one of the differentially expressed miRNA identified in this study. It was reported that miR-23b inhibited the expression of the transforming growth factor-beta (TGF-β) signaling (Finnerty et al., 2010). In the current study, miR-23b and 5 other miRNAs (miR-2454-3p, miR-496, miR-503-3p, miR-6520, and novel-6) were predicted to regulate STAT5B, which is known to be involved in TGF-β signaling (Passerini et al., 2008; Hosui et al., 2009). In addition, genes that are known to affect milk traits (CSN3, CSN2, LALBA, DGAT2, STAT5B, and SCD) were predicted to be targets of some of the identified miRNAs, which implied that they may play critical regulatory roles in mammary gland development and milk composition.
It was found that 21 of 31 differentially expressed genes detected in our previous study (Cui et al., 2014) were the predicted targets for some of the 71 differentially expressed miRNAs detected in the present study. Serum amyloid A1 (SAA1), serum amyloid A1 (SAA3), and mammary serum amyloid A3.2 (M-SAA3.2) were predicted to be regulated by miR-146b (SAA1 was also regulated by miR-125a and miR-125b); VEGFA was regulated by miR-125a, miR-125b, miR-106b, and miR-2904; and ribosomal protein L23a(RPL23A), tribbles homolog 3 (TRIB3), and PTHLH were regulated by miR-378c, miR-2904, and miR-106b, respectively. Moreover, Cai et al., performed RNA sequencing with mammary gland tissue samples from six Chinese Holstein cows with three extremely high and three low milk protein percentage phenotypes and miR-2904, miR-339b, miR-146b, miR-339a, miR-29c, miR-106b, miR-190a, miR-21-3p, miR-15a, miR-486, miR-135, miR-101a, miR-152 and miR-139 were found differentially expressed, which were also identified in our study and targeted on four differentially expressed genes (TRIB3, PTHLH, VEGFA, and M-SAA3.2). These seven genes represent the most promising candidates may affect milk PP and FP in dairy cattle (Cui et al., 2014). Specifically, miR-146b was reported to be involved mainly in leukemia, epidermal growth factor receptor (EGFR), MAPK, and nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) signaling pathways (Mathews et al., 2004; Taganov et al., 2006; Xiang et al., 2014). The EGFR and MAPK signaling pathways have been demonstrated to be related to adipocyte differentiation (Devaraj et al., 2009; Gao and Bing, 2011) and the NF-κB pathway controls the DNA transcription protein complexes. In human study, miR-146b was shown to regulate the NF-κB signaling pathway in which breast cancer metastasis suppressor 1 (BRMS1) has already been implicated, and inhibited both migration and invasion related to metastasis (Taganov et al., 2006; Xiang et al., 2014). Members of the miR-125 family were reported to be implicated in a variety of carcinomas and other diseases as either repressors or promoters. Sun et al. (2013) found that up-regulated miR-125 significantly inhibited the expression of VEGFA both in vitro and in vivo (Jiang et al., 1997). The miR-125 family was found to be a NF-κB-dependent gene in the study by Kim et al. (2012). miR-378c was shown to be involved in the regulation of RPL23A, which plays a critical role in translation and participates in apoptosis, cell division, and differentiation (Wool, 1996; Fang et al., 2012; Knezevic et al., 2012). This is consistent with previous reported study where miR-378c was found associated with apoptosis (Lee et al., 2007; Fang et al., 2012; Knezevic et al., 2012; Wang et al., 2014).
The GO and KEGG pathway analyses indicated that VEGFA, NR4A1, DDIT3, and CDKN1A were involved in the MAPK, mTOR, HIF-1, and PI3K-Akt signaling pathways, respectively. These four genes were predicted as target genes for miR-106b, miR-2904, miR-125a(b), miR-21-3p, miR-224, miR-31, miR-345-5p, and miR-3431. mTOR signaling is known as playing a fundamental role in adipogenesis (Laplante and Sabatini, 2009), which is the process that leads to the formation of adipose tissue and the most important energy storage site in mammals. It has been demonstrated that mTORC1 positively regulates the activity of sterol regulatory element binding protein 1 (SREBP1) and peroxisome proliferator-activated receptor gamma (PPARG) (Benmoussa et al., 2020), which are two transcription factors that control the expression of genes encoding proteins involved in lipid and cholesterol homeostasis (Kim and Chen, 2004; Porstmann et al., 2008; Kim et al., 2012). HIF-1 is a heterodimeric transcription factor that increases the phosphorylation of signal transducer and activator of transcription 5A (STAT5A) in mammary epithelial cells, and the phosphorylation of STAT5 is known to play important roles in the regulation of milk protein gene expression and mammary development (Shao and Zhao, 2014; Benmoussa et al., 2020). Several studies have shown that hypoxia causes mammary epithelial disorganization and induces a cancer cell-like phenotype in human mammary epithelial cells (MECs) (Whelan et al., 2010; Whelan and Reginato, 2011; Vaapil et al., 2012). The PI3K-Akt pathway has important functions in mammary gland development and function (Wickenden and Watson, 2010). One of the most important functions of Akt is the regulation of glucose homeostasis and metabolism, particularly in muscle and fat tissues (Enright et al., 2003). Therefore, these miRNAs could play critical roles in regulating formation of milk composition trait.
Considering that microRNAs regulate gene expression by targeting specific sequences in the 3′UTR of their cognate genes (Lewis et al., 2005; Friedman et al., 2009), the regulatory roles of some miRNAs on their predicted targets were verified using dual luciferase report assay transfected with mimics, inhibitors and mutants of seed sequences. The results demonstrated that miR-2904, miR-29c/miR-146b/miR-339a, miR-339b/miR-106b/miR-190a indeed down-regulated the expression of the TRIB3, M-SAA3.2 and PTHLH, respectively. The molecular mechanisms of how these miRNAs regulate their targets will be further validated through RNAi and over-expression in bovine mammary epithelial cell lines. In addition, it is generally recognized that miRNAs regulate the expression of target genes by inhibiting their translation or inducing degradation of the target mRNAs in animal cells. However, several predicted target genes were regulated in the same direction of expression as those of the corresponding miRNAs between high and low groups. The reason could be due to either target prediction error of the current commonly used prediction softwares (TargetScan 6.2 and MiRanda) or some unknown biological mechanisms. Actually, target prediction was only the first step for studies on interaction between miRNA and their targets. The miRNAs and targets with reverse expression patterns will be considered as the key components for further validation.
In this study, only two biological replicates, which were the same as in our previous RNA-seq investigation (Cui et al., 2014), were used for each condition due to the availability of mammary gland sample from lactating cows, especially high production ones. In order to minimize false-positive errors and ensure substantial detection power and accuracy, two strategies were applied to detect the differentially expressed miRNAs between milking Holstein cows with high PP and FP and cows with low PP and FP, by controlling the critical influencing factors. Small RNA transcripts were deeply sequenced (9-10G data per transcriptome), and only those differentially expressed miRNAs ranked in the top half of the expressed miRNAs were considered, as suggested by Rapaport et al. (2013), Trapnell et al. (2013). Rapaport et al. (2013) investigated the impact of different sequencing depths and number of replicates on the identification of differentially expressed genes, where the authors demonstrated that with most methods, over 90% of differently expressed genes at the top expression levels could be detected with using two replicates and 5% of the reads (Rapaport et al., 2013; Trapnell et al., 2013). The differentially expressed miRNAs expressed in the bottom half level were eliminated to ensure the power in detection. Although mRNA sequencing data was used in this study, detection of differentially expressed miRNAs is based on same statistical theory and software (Anders and Huber, 2010). However, more biological replications are still preferred and recommended in order to provide broader application (Rapaport et al., 2013; Trapnell et al., 2013). The more replicates are performed, the more the detection power is improved. The potential regulatory roles on target genes from such differentially expressed miRNAs will be validated further by performing more in-depth investigation.
Conclusion
Using sRNA sequencing, 497 known bovine miRNAs and 49 novel bovine miRNAs were identified in the mammary glands of lactating dairy cows. Among all these miRNAs, 71 were differentially expressed between cows with the high and low milk PP and FP. Combined with our previous RNA-Seq data, 21 differentially expressed genes were predicted as the targets for some of the 71 differentially expressed miRNAs. Biological processes related to protein metabolism, fat metabolism, and mammary gland development were enriched for some of the identified miRNAs, which indicated that they may play critical roles in regulating of milk protein and fat traits in dairy cattle. Expression of the TRIB3, M-SAA3.2 and PTHLH were significantly down-regulated by miR-2904, miR-29c/miR-146b/miR-339a and miR-339b/miR-106b/miR-190a through binding the specific target sequences in 3′UTR of these genes, respectively.
Data Availability Statement
The datasets GENERATED for this study can be found in NCBI SRA Accession Numbers: (SRR3631014, SRR3631016, SRR3631053, and SRR3631054).
Ethics Statement
All experiments, including all protocols for collection of the mammary gland tissues of experimental individuals and phenotypic observations, were performed in accordance with relevant guidelines and regulations of the Institutional Animal Care and Use Committee (IACUC) at China Agricultural University who have reviewed and approved the experiments. Samples were collected specifically for this study following standard procedures with the full agreement of the Beijing Sanyuanlvhe Dairy Farming Center who owned the animals.
Author Contributions
XC and MY performed the RNA-related experiments, data analysis, and wrote the manuscript. SZ, CW, and QZ participated in the experimental design and helped interpret the results. XG participated in the interpretation of the results, writing, and revision of the manuscript. DS conceived and designed the experiments, and wrote the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National Natural Science Foundation (31872330 and 31802041), National Science and Technology Programs of China (2013AA102504), Shanxi Province Science Foundation for Youths (201801D221285), and Program for Changjiang Scholar and Innovation Research Team in University (IRT_15R62).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We appreciate the Beijing Dairy Cattle Center and the Beijing Sanyuanlvhe Dairy Farming Center, for providing the mammary gland samples, milk traits records and pedigree data for the Chinese Holstein cows used in this study.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.548268/full#supplementary-material
Supplementary Figure 1 | The length distribution of the reads.
Supplementary Figure 2 | Correlation plots of the reads for two groups.
Supplementary Table 1 | Statistics of the small RNA sequencing data from four mammary gland samples.
Supplementary Table 2 | Annotations of the unique mapped sRNA reads from four mammary gland samples.
Supplementary Table 3 | Profile of the known bovine miRNAs.
Supplementary Table 4 | Profile of the bovine novel miRNAs.
Supplementary Table 5 | Prediction results of target genes for the known and novel miRNAs.
Supplementary Table 6 | GO terms assigned to the predicted target genes of the known and novel miRNAs.
Supplementary Table 7 | Prediction results of target genes for the 71 differentially expressed miRNAs.
Supplementary Table 8 | Fifteen randomly selected miRNAs for q-PCR experiments and their primers.
Abbreviations
PP, protein percentages; FP, fat percentages; MAPK, mitogen-activated protein kinases; mTOR, mammalian target of rapamycin; HIF-1, hypoxia inducible factor-1; PI3K-Akt, phosphatidylinositol 3 kinase-protein kinase B; IACUC, Institutional Animal Care and Use Committee; RIN, RNA integrity number; CPM, Counts per million; UTR, un-translated region; HEK, human embryonic kidney; DMEM, Dulbecco’s modified Eagle’s medium; FDR, false discovery rate; CSN2, β -casein; CSN3, κ -casein; LALBA, α -lactalbumin; DGAT2, diacylglycerol O-acyltransferase 2; GHR, growth hormone receptor; STAT5B, signal transducer and activator of transcription 5B; SCD, stearoyl-coenzyme A desaturase; QTL, quantitative trait loci; GWAS, genome-wide association studies; TRIB3, tribbles homolog 3; SAA1, serum amyloid A1; SAA3, serum amyloid A3 (SAA3); M-SAA 3.2, mammary serum amyloid A3; VEGFA, vascular endothelial growth factor A; PTHLH, parathyroid hormone-like hormone; RPL23A, ribosomal protein L23A; DDIT3, DNA-damage-inducible transcript 3; NR4A1, nuclear receptor subfamily 4; group A, member 1; VEGFA, vascular endothelial growth factor A; CDKN1A, cyclin-dependent kinase inhibitor 1A; ATF3, activating transcription factor 3; CHAC1, cation transport regulator homolog 1; TGF- β, transforming growth factor-beta; EGFR, epidermal growth factor receptor; NF- κ B, nuclear factor kappa-light-chain-enhancer of activated B cells; BRMS1, breast cancer metastasis suppressor 1; SREBP1, sterol regulatory element binding protein 1; PPARG, peroxisome proliferator-activated receptor gamma; STAT5A, signal transducer and activator of transcription 5A; MEC, mammary epithelial cells.
Footnotes
- ^ 1http://www.mirbase.org/
- ^ 2http://www.targetscan.org/vert_60/
- ^ 3http://www.microrna.org/microrna/home.do
- ^ 4http://geneontology.org/
- ^ 5http://www.genome.jp/kegg/
- ^ 6http://kobas.cbi.pku.edu.cn/home.do/
References
Aggarwal, P., Turner, A., Matter, A., Kattman, S. J., Stoddard, A., Lorier, R., et al. (2014). Rna expression profiling of human ipsc-derived cardiomyocytes in a cardiac hypertrophy model. PLoS One 9:e108051. doi: 10.1371/journal.pone.0108051
Almughlliq, F. B., Koh, Y. Q., Peiris, H. N., Vaswani, K., Holland, O., Meier, S., et al. (2019). Circulating exosomes may identify biomarkers for cows at risk for metabolic dysfunction. Sci. Rep. 9:13879.
Anders, S., and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol. 11:R106.
Aqeilan, R. I., Calin, G. A., and Croce, C. M. (2010). miR-15a and miR-16-1 in cancer: discovery, function and future perspectives. Cell Death Differ. 17, 215–220. doi: 10.1038/cdd.2009.69
Aravin, A., and Tuschl, T. (2005). Identification and characterization of small RNAs involved in RNA silencing. FEBS Lett. 579, 5830–5840. doi: 10.1016/j.febslet.2005.08.009
Bandi, N., Zbinden, S., Gugger, M., Arnold, M., Kocher, V., and Hasan, L. (2009). miR-15a and miR-16 are implicated in cell cycle regulation in a Rb-dependent manner and are frequently deleted or down-regulated in non-small cell lung cancer. Cancer Res. 69, 5553–5559. doi: 10.1158/0008-5472.can-08-4277
Barbu, M. G., Condrat, C. E., Thompson, D. C., Bugnar, O. L., Cretoiu, D., Toader, O. D., et al. (2020). MicroRNA Involvement in Signaling Pathways During Viral Infection. Front. Cell Dev. Biol. 8:143. doi: 10.3389/fcell.2020.00143
Benmoussa, A., Laugier, J., Beauparlant, C. J., Lambert, M., Droit, A., and Provost, P. (2020). Complexity of the microRNA transcriptome of cow milk and milk-derived extracellular vesicles isolated via differential ultracentrifugation. J. Dairy Sci. 103, 16–29. doi: 10.3168/jds.2019-16880
Berezikov, E., Cuppen, E., and Plasterk, R. H. (2006a). Approaches to microRNA discovery. Nat. Genet. 38(Suppl.), S2–S7.
Berezikov, E., Van Tetering, G., Verheul, M., Van De Belt, J., Van Laake, L., Vos, J., et al. (2006b). Many novel mammalian microRNA candidates identified by extensive cloning and RAKE analysis. Genome Res. 16, 1289–1298. doi: 10.1101/gr.5159906
Billa, P. A., Faulconnier, Y., Ye, T., Chervet, M., Le Provost, F., Pires, J. A. A., et al. (2019). Deep rna-seq reveals mirnome differences in mammary tissue of lactating holstein and montbéliarde cows. BMC Genomics 20:621. doi: 10.1186/s12864-019-5987-4
Bionaz, M., and Loor, J. J. (2008). Gene networks driving bovine milk fat synthesis during the lactation cycle. BMC Genomics 9:366. doi: 10.1186/1471-2164-9-366
Bionaz, M., and Loor, J. J. (2011). Gene networks driving bovine mammary protein synthesis during the lactation cycle. Bioinform. Biol. Insights 5, 83–98.
Bionaz, M., Periasamy, K., Rodriguez-Zas, S. L., Everts, R. E., Lewin, H. A., Hurley, W. L., et al. (2012). Old and new stories: revelations from functional analysis of the bovine mammary transcriptome during the lactation cycle. PLoS One 7:e33268. doi: 10.1371/journal.pone.0033268
Bonci, D., Coppola, V., Musumeci, M., Addario, A., Giuffrida, R., and Memeo, L. (2008). The miR-15a-miR-16-1 cluster controls prostate cancer by targeting multiple oncogenic activities. Nat. Med. 14, 1271–1277. doi: 10.1038/nm.1880
Cai, W., Li, C., Liu, S., Zhou, C., Yin, H., and Song, J. (2018). Genome Wide Identification of Novel Long Non-coding RNAs and Their Potential Associations With Milk Proteins in Chinese Holstein Cows. Front. Genet. 9:281. doi: 10.3389/fgene.2018.00281
Chen, B., Li, H., Zeng, X., Yang, P., Liu, X., and Zhao, X. (2012). Roles of microRNA on cancer cell metabolism. J. Transl. Med. 10:228. doi: 10.1186/1479-5876-10-228
Chen, P. Y., Manninga, H., Slanchev, K., Chien, M., Russo, J. J., Ju, J., et al. (2005). The developmental miRNA profiles of zebrafish as determined by small RNA cloning. Genes Dev. 19, 1288–1293. doi: 10.1101/gad.1310605
Cimmino, A., Calin, G. A., Fabbri, M., Iorio, M. V., Ferracin, M., Shimizu, M., et al. (2005). miR-15 and miR-16 induce apoptosis by targeting BCL2. Proc. Natl. Acad. Sci. U S A. 102, 13944–13949. doi: 10.1073/pnas.0506654102
Coutinho, L. L., Matukumalli, L. K., Sonstegard, T. S., Van Tassell, C. P., Gasbarre, L. C., Capuco, A. V., et al. (2007). Discovery and profiling of bovine microRNAs from immune-related and embryonic tissues. Physiol. Genomics 29, 35–43. doi: 10.1152/physiolgenomics.00081.2006
Cui, X., Hou, Y., Yang, S., Xie, Y., Zhang, S., Zhang, Y., et al. (2014). Transcriptional profiling of mammary gland in Holstein cows with extremely different milk protein and fat percentage using RNA sequencing. BMC Genomics 15:226. doi: 10.1186/1471-2164-15-226
Devaraj, S., Yun, J. M., Duncan-Staley, C., and Jialal, I. (2009). C-reactive protein induces M-CSF release and macrophage proliferation. J. Leukoc Biol. 85, 262–267. doi: 10.1189/jlb.0808458
Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA targets in Drosophila. Genome Biol. 5:R1.
Fang, J., Song, X. W., Tian, J., Chen, H. Y., Li, D. F., Wang, J. F., et al. (2012). Overexpression of microRNA-378 attenuates ischemia-induced apoptosis by inhibiting caspase-3 expression in cardiac myocytes. Apoptosis 17, 410–423. doi: 10.1007/s10495-011-0683-0
Finnerty, J. R., Wang, W. X., Hebert, S. S., Wilfred, B. R., Mao, G., and Nelson, P. T. (2010). The miR-15/107 group of microRNA genes: evolutionary biology, cellular functions, and roles in human diseases. J. Mol. Biol. 402, 491–509. doi: 10.1016/j.jmb.2010.07.051
Finucane, K. A., Mcfadden, T. B., Bond, J. P., Kennelly, J. J., and Zhao, F. Q. (2008). Onset of lactation in the bovine mammary gland: gene expression profiling indicates a strong inhibition of gene expression in cell proliferation. Funct. Integr. Genomics 8, 251–264. doi: 10.1007/s10142-008-0074-y
Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucl. Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688
Friedman, R. C., Farh, K. K., Burge, C. B., and Bartel, D. P. (2009). Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 19, 92–105. doi: 10.1101/gr.082701.108
Gao, D., and Bing, C. (2011). Macrophage-induced expression and release of matrix metalloproteinase 1 and 3 by human preadipocytes is mediated by IL-1β via activation of MAPK signaling. J. Cell Physiol. 226, 2869–2880. doi: 10.1002/jcp.22630
Giurato, G., De Filippo, M. R., Rinaldi, A., Hashim, A., Nassa, G., Ravo, M., et al. (2013). Imir: An integrated pipeline for high-throughput analysis of small non-coding rna data obtained by smallrna-seq. BMC Bioinformatics 14:362. doi: 10.1186/1471-2105-14-362
Glazov, E. A., Cottee, P. A., Barris, W. C., Moore, R. J., Dalrymple, B. P., and Tizard, M. L. (2008). A microRNA catalog of the developing chicken embryo identified by a deep sequencing approach. Genome Res. 18, 957–964. doi: 10.1101/gr.074740.107
Gu, Z., Eleswarapu, S., and Jiang, H. (2007). Identification and characterization of microRNAs from the bovine adipose tissue and mammary gland. FEBS Lett. 581, 981–988. doi: 10.1016/j.febslet.2007.01.081
Gupta, V., Markmann, K., Pedersen, C. N., Stougaard, J., and Andersen, S. U. (2012). Shortran: A pipeline for small rna-seq data analysis. Bioinformatics 28, 2698–2700. doi: 10.1093/bioinformatics/bts496
Hennighausen, L., and Robinson, G. W. (2005). Information networks in the mammary gland. Nat. Rev. Mol. Cell Biol. 6, 715–725.
Hosui, A., Kimura, A., Yamaji, D., Zhu, B. M., Na, R., and Hennighausen, L. (2009). Loss of STAT5 causes liver fibrosis and cancer development through increased TGF-{beta} and STAT3 activation. J. Exp. Med. 206, 819–831. doi: 10.1084/jem.20080003
Huntzinger, E., and Izaurralde, E. (2011). Gene silencing by microRNAs: contributions of translational repression and mRNA decay. Nat. Rev. Genet. 12, 99–110. doi: 10.1038/nrg2936
Jiang, H., Lin, J. J., Tao, J., and Fisher, P. B. (1997). Suppression of human ribosomal protein L23A expression during cell growth inhibition by interferon-beta. Oncogene 14, 473–480. doi: 10.1038/sj.onc.1200858
Jin, W., Dodson, M. V., Moore, S. S., Basarab, J. A., and Guan, L. L. (2010). Characterization of microRNA expression in bovine adipose tissues: a potential regulatory mechanism of subcutaneous adipose tissue development. BMC Mol. Biol. 11:29. doi: 10.1186/1471-2199-11-29
Jin, W., Grant, J. R., Stothard, P., Moore, S. S., and Guan, L. L. (2009). Characterization of bovine miRNAs by sequencing and bioinformatics analysis. BMC Mol. Biol. 10:90. doi: 10.1186/1471-2199-10-90
Ju, Z., Jiang, Q., Liu, G., Wang, X., Luo, G., Zhang, Y., et al. (2018). Solexa sequencing and custom microrna chip reveal repertoire of micrornas in mammary gland of bovine suffering from natural infectious mastitis. Anim. Genet. 49, 3–18. doi: 10.1111/age.12628
Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucl. Acids Res. 36, D480–D484.
Kim, J. E., and Chen, J. (2004). regulation of peroxisome proliferator-activated receptor-gamma activity by mammalian target of rapamycin and amino acids in adipogenesis. Diabetes 53, 2748–2756. doi: 10.2337/diabetes.53.11.2748
Kim, S. W., Ramasamy, K., Bouamar, H., Lin, A. P., Jiang, D., and Aguiar, R. C. (2012). MicroRNAs miR-125a and miR-125b constitutively activate the NF-κB pathway by targeting the tumor necrosis factor alpha-induced protein 3 (TNFAIP3. A20). Proc. Natl. Acad. Sci. U S A. 109, 7865–7870. doi: 10.1073/pnas.1200081109
Knezevic, I., Patel, A., Sundaresan, N. R., Gupta, M. P., Solaro, R. J., Nagalingam, R. S., et al. (2012). A novel cardiomyocyte-enriched microRNA, miR-378, targets insulin-like growth factor 1 receptor: implications in postnatal cardiac remodeling and cell survival. J. Biol. Chem. 287, 12913–12926. doi: 10.1074/jbc.m111.331751
Kuksa, P. P., Amlie-Wolf, A., Katanic, Ž, Valladares, O., Wang, L. S., and Leung, Y. Y. (2018). Spar: Small rna-seq portal for analysis of sequencing experiments. Nucl. Acids Res. 46, W36–W42.
Lagos-Quintana, M., Rauhut, R., Meyer, J., Borkhardt, A., and Tuschl, T. (2003). New microRNAs from mouse and human. Rna 9, 175–179. doi: 10.1261/rna.2146903
Landgraf, P., Rusu, M., Sheridan, R., Sewer, A., Iovino, N., Aravin, A., et al. (2007). A mammalian microRNA expression atlas based on small RNA library sequencing. Cell 129, 1401–1414.
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25.
Laplante, M., and Sabatini, D. M. (2009). An emerging role of mTOR in lipid biosynthesis. Curr. Biol. 19, R1046–R1052.
Le Guillou, S., Sdassi, N., Laubier, J., Passet, B., Vilotte, M., Castille, J., et al. (2012). Overexpression of miR-30b in the developing mouse mammary gland causes a lactation defect and delays involution. PLoS One 7:e45727. doi: 10.1371/journal.pone.0045727
Lee, D. Y., Deng, Z., Wang, C. H., and Yang, B. B. (2007). MicroRNA-378 promotes cell survival, tumor growth, and angiogenesis by targeting SuFu and Fus-1 expression. Proc. Natl. Acad. Sci. U S A. 104, 20350–20355. doi: 10.1073/pnas.0706901104
Lewis, B. P., Burge, C. B., and Bartel, D. P. (2005). Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell 120, 15–20. doi: 10.1016/j.cell.2004.12.035
Li, H. M., Wang, C. M., Li, Q. Z., and Gao, X. J. (2012a). MiR-15a decreases bovine mammary epithelial cell viability and lactation and regulates growth hormone receptor expression. Molecules 17, 12037–12048. doi: 10.3390/molecules171012037
Li, Q., Qiao, J., Zhang, Z., Shang, X., Chu, Z., Fu, Y., et al. (2020). Identification and analysis of differentially expressed long non-coding rnas of chinese holstein cattle responses to heat stress. Anim. Biotechnol. 31, 9–16. doi: 10.1080/10495398.2018.1521337
Li, R., Zhang, C. L., Liao, X. X., Chen, D., Wang, W. Q., and Zhu, Y. H. (2015). Transcriptome microRNA profiling of bovine mammary glands infected with Staphylococcus aureus. Int. J. Mol. Sci. 16, 4997–5013. doi: 10.3390/ijms16034997
Li, Z., Liu, H., Jin, X., Lo, L., and Liu, J. (2012b). Expression profiles of microRNAs from lactating and non-lactating bovine mammary glands and identification of miRNA related to lactation. BMC Genomics 13, 731. doi: 10.1186/1471-2164-13-731
Long, J. E., and Chen, H. X. (2009). Identification and characteristics of cattle microRNAs by homology searching and small RNA cloning. Biochem. Genet. 47, 329–343. doi: 10.1007/s10528-009-9234-6
Luoreng, Z. M., Wang, X. P., Mei, C. G., and Zan, L. S. (2018). Comparison of microrna profiles between bovine mammary glands infected with staphylococcus aureus and escherichia coli. Int. J. Biol. Sci. 14, 87–99. doi: 10.7150/ijbs.22498
Mao, X., Cai, T., Olyarchuk, J. G., and Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787–3793. doi: 10.1093/bioinformatics/bti430
Martello, G., Rosato, A., Ferrari, F., Manfrin, A., Cordenonsi, M., Dupont, S., et al. (2010). A MicroRNA targeting dicer for metastasis control. Cell 141, 1195–1207. doi: 10.1016/j.cell.2010.05.017
Mathews, D. H., Disney, M. D., Childs, J. L., Schroeder, S. J., Zuker, M., and Turner, D. H. (2004). Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure. Proc. Natl. Acad. Sci. U S A. 101, 7287–7292. doi: 10.1073/pnas.0401799101
Mobuchon, L., Marthey, S., Le Guillou, S., Laloe, D., Le Provost, F., and Leroux, C. (2015). Food Deprivation Affects the miRNome in the Lactating Goat Mammary Gland. PLoS One 10:e0140111. doi: 10.1371/journal.pone.0140111
Muroya, S., Ogasawara, H., Nohara, K., Oe, M., Ojima, K., and Hojito, M. (2019). Coordinate alteration of mRNA-microRNA transcriptomes associated with exosomes and fatty acid metabolism in adipose tissue and skeletal muscle in grazing cattle. Asian Australas J. Anim. Sci. 33, 1824–1836. doi: 10.5713/ajas.19.0682
Passerini, L., Allan, S. E., Battaglia, M., Di Nunzio, S., Alstad, A. N., Levings, M. K., et al. (2008). STAT5-signaling cytokines regulate the expression of FOXP3 in CD4+CD25+ regulatory T cells and CD4+CD25- effector T cells. Int. Immunol. 20, 421–431. doi: 10.1093/intimm/dxn002
Porstmann, T., Santos, C. R., Griffiths, B., Cully, M., Wu, M., and Leevers, S. (2008). SREBP activity is regulated by mTORC1 and contributes to Akt-dependent cell growth. Cell Metab. 8, 224–236. doi: 10.1016/j.cmet.2008.07.007
Pu, J., Li, R., Zhang, C., Chen, D., Liao, X., and Zhu, Y. (2017). Expression profiles of mirnas from bovine mammary glands in response to streptococcus agalactiae-induced mastitis. J. Dairy Res. 84, 300–308. doi: 10.1017/s0022029917000437
Rapaport, F., Khanin, R., Liang, Y., Pirun, M., Krek, A., Zumbo, P., et al. (2013). Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data. Genome Biol. 14:R95.
Rottiers, V., and Näär, A. M. (2012). MicroRNAs in metabolism and metabolic disorders. Nat. Rev. Mol. Cell Biol. 13, 239–250. doi: 10.1038/nrm3313
Shao, Y., and Zhao, F. Q. (2014). Emerging evidence of the physiological role of hypoxia in mammary development and lactation. J. Anim. Sci. Biotechnol. 5:9.
Shen, B., Zhang, L., Lian, C., Lu, C., Zhang, Y., Pan, Q., et al. (2016). Deep Sequencing and Screening of Differentially Expressed MicroRNAs Related to Milk Fat Metabolism in Bovine Primary Mammary Epithelial Cells. Int. J. Mol. Sci. 17:200. doi: 10.3390/ijms17020200
Silveri, L., Tilly, G., Vilotte, J. L., and Le Provost, F. (2006). MicroRNA involvement in mammary gland development and breast cancer. Reprod. Nutr. Dev. 46, 549–556. doi: 10.1051/rnd:2006026
Sun, Y. M., Lin, K. Y., and Chen, Y. Q. (2013). Diverse functions of miR-125 family in different cell contexts. J. Hematol. Oncol. 6:6. doi: 10.1186/1756-8722-6-6
Taganov, K. D., Boldin, M. P., Chang, K. J., and Baltimore, D. (2006). NF-kappaB-dependent induction of microRNA miR-146, an inhibitor targeted to signaling proteins of innate immune responses. Proc. Natl. Acad. Sci. U S A. 103, 12481–12486. doi: 10.1073/pnas.0605298103
Trapnell, C., Hendrickson, D. G., Sauvageau, M., Goff, L., Rinn, J. L., and Pachter, L. (2013). Differential analysis of gene regulation at transcript resolution with rna-seq. Nat. Biotechnol. 31, 46–53. doi: 10.1038/nbt.2450
Vaapil, M., Helczynska, K., Villadsen, R., Petersen, O. W., Johansson, E., Beckman, S., et al. (2012). Hypoxic conditions induce a cancer-like phenotype in human breast epithelial cells. PLoS One 7:e46543. doi: 10.1371/journal.pone.0046543
Wang, K. Y., Ma, J., Zhang, F. X., Yu, M. J., Xue, J. S., and Zhao, J. S. (2014). MicroRNA-378 inhibits cell growth and enhances L-OHP-induced apoptosis in human colorectal cancer. IUBMB Life 66, 645–654. doi: 10.1002/iub.1317
Wen, M., Shen, Y., Shi, S., and Tang, T. (2012). miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics 13:140. doi: 10.1186/1471-2105-13-140
Whelan, K. A., and Reginato, M. J. (2011). Surviving without oxygen: hypoxia regulation of mammary morphogenesis and anoikis. Cell Cycle 10, 2287–2294. doi: 10.4161/cc.10.14.16532
Whelan, K. A., Caldwell, S. A., Shahriari, K. S., Jackson, S. R., Franchetti, L. D., Johannes, G. J., et al. (2010). Hypoxia suppression of Bim and Bmf blocks anoikis and luminal clearing during mammary morphogenesis. Mol. Biol. Cell 21, 3829–3837. doi: 10.1091/mbc.e10-04-0353
Wickenden, J. A., and Watson, C. J. (2010). Key signalling nodes in mammary gland development and cancer. Signalling downstream of PI3 kinase in mammary epithelium: a play in 3 Akts. Breast Cancer Res. 12:202.
Wool, I. G. (1996). Extraribosomal functions of ribosomal proteins. Trends Biochem. Sci. 21, 164–165. doi: 10.1016/s0968-0004(96)20011-8
Xiang, M., Birkbak, N. J., Vafaizadeh, V., Walker, S. R., Yeh, J. E., Liu, S., et al. (2014). STAT3 induction of miR-146b forms a feedback loop to inhibit the NF-κB to IL-6 signaling axis and STAT3-driven cancer phenotypes. Sci. Signal 7:ra11. doi: 10.1126/scisignal.2004497
Yang, B., Jiao, B., Ge, W., Zhang, X., Wang, S., Zhao, H., et al. (2018). Transcriptome sequencing to detect the potential role of long non-coding rnas in bovine mammary gland during the dry and lactation period. BMC Genomics 19:605. doi: 10.1186/s12864-018-4974-5
Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11:R14.
Zhang, C., Wu, H., Wang, Y., Zhu, S., Liu, J., Fang, X., et al. (2016). Circular rna of cattle casein genes are highly expressed in bovine mammary gland. J. Dairy Sci. 99, 4750–4760. doi: 10.3168/jds.2015-10381
Keywords: mammary gland, mRNA, miRNA, RNA-seq, dairy cattle
Citation: Cui X, Zhang S, Zhang Q, Guo X, Wu C, Yao M and Sun D (2020) Comprehensive MicroRNA Expression Profile of the Mammary Gland in Lactating Dairy Cows With Extremely Different Milk Protein and Fat Percentages. Front. Genet. 11:548268. doi: 10.3389/fgene.2020.548268
Received: 02 April 2020; Accepted: 05 November 2020;
Published: 03 December 2020.
Edited by:
Shu-Hong Zhao, Huazhong Agricultural University, ChinaReviewed by:
Jun Luo, Northwest A&F University, ChinaZhibin Ji, Shandong Agricultural University, China
Copyright © 2020 Cui, Zhang, Zhang, Guo, Wu, Yao and Sun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Dongxiao Sun, sundx@cau.edu.cn; Xiaogang Cui, helloxiaogangcui@sxu.edu.cn; 398304310@qq.com