Characterization of miRNA profiles in the mammary tissue of dairy cattle in response to heat stress

MicroRNAs (miRNAs) are a class of small noncoding RNAs that play important roles in the regulation of gene expression. However, the role of miRNAs in bovine mammary gland responses to heat stress is not well understood. In the present study, we performed a deep RNA sequencing analysis to identify miRNAs associated with the heat stress potential of the bovine mammary gland. We identified 27 miRNAs that were differentially expressed significantly between the mammary tissue of Holstein cattle heat stress and normal conditions. Twenty miRNAs had higher expression in the mammary tissue of heat-stressed Holstein cattle. The seven highest differentially expressed candidate miRNAs (bta-miR-21-5p, bta-miR-99a-5p, bta-miR-146b, bta-miR-145, bta-miR-2285 t, bta-miR-133a, and bta-miR-29c) identified by deep RNA sequencing were additionally evaluated by stem-loop qPCR. Enrichment analyses for targeted genes revealed that the major differences between miRNAs expression in the mammary gland of heat-stressed versus control were associated with the regulation of Wnt, TGF-β, MAPK, Notch, and JAK-STAT. These data indicated that the differentially expressed miRNAs identified in this study may act as dominant regulators during heat stress. We might reduce heat stress damage of Holstein cows by up-regulating or down-regulating these differentially expressed miRNAs.


Background
Stress can be defined as an external condition that produces a "strain" in a biological system [1]. The environmental stress may be measured by changes in body temperature, metabolic rate, or productivity. Heat stress negatively impacts all features of dairy cattle production including milk composition and mammary gland pathogens. It substantially influences a cow's growth and development [2]. Reduction in reproductive performance of lactating cows during summer is associated with decreased thermoregulatory competence [3]. Heat stress causes cow metabolic disorders and a reduction in milk production [4,5], and it also decreases immunity and increases susceptibility to mastitis, endometritis disease and even death in severe cases [6][7][8]. Diminished milk yield and reproductive losses during summer months seriously affect the economic potential of the dairy industry. In addition, global warming may boost the occurrence of heat stress [9]. Thus, for the dairy industry, heat stress has been a bottleneck which limits the efficiency of the dairy supply throughout the year.
The heat stress response is a complex molecular process that involves the transcriptional and post-transcriptional regulation of stress-related genes. Acute environmental change initiates the heat stress response at the cellular level. Changes in gene expression are associated with a reaction to an environmental stressor as well as changes across a variety of organs and tissues associated with the acclimation response. Functional genomics establishes a verifiable link between gene expression and phenotype. Endogenous noncoding small RNAs known as microRNAs (miRNAs) are increasingly being recognized as important modulators of gene expression at the post-transcriptional level and have been shown to be involved in diverse biological processes such as differentiation, development, apoptosis, and viral infection. RNA-Seq, in particular, allows a global analysis of gene expression responses to environmental change.
Some miRNAs have been shown to be involved in plant stress responses by down-regulating the respective target genes encoding regulatory and functional proteins [10]. Differential expression of miRNAs is also associated with thermal stress in cattle. Down-regulation of miR-181a can reduce heat stress damage in peripheral blood mononuclear cells (PBMCs) of Holstein cows [11]. miRNA profiles in bovine mammary tissue infected with Staphylococcus [12,13] and miRNA profiles in serum or PBMC cell of heat-stressed cattle have been studied [14][15][16]; however, miRNA profiles in bovine mammary gland of stressed lactating cattle have not been compared to normal lactating cattle. Therefore, it was our objective to profile miRNA expression under heat stress using bovine mammary glands. To investigate the role of miRNA in the heat stress response, miRNA expression in bovine mammary glands were characterized by next-generation sequencing during summer and spring with and without a heat stress challenge.

Tissue collection
Samples were collected from eight mammary gland tissues of the four lactating Chinese Holstein cows from Yucheng dairy farm, China. The samples were collected in two different environmental seasons, viz. spring and summer, with the temperature ranges between 15 and 20°C (March; designated as non-heat stressed, NHS) and 30-38°C (July; designated as heat stressed, HS), respectively. The temperature humidity index (THI) was used as a heat stress indicator. The temperature probing procedure has been described in detail elsewhere [17]. Dry-bulb and wet-bulb temperatures were recorded using a dry and wet bulb thermometer consisting of two thermometers, a dry bulb thermometer and a wet bulb thermometer. THI was calculated as THI = 0.72 (Td + Tw) + 40.6 where Td is the dry-bulb temperature and Tw is the wet-bulb temperature. Stress response of the animals was characterized by recording physiological parameters. The rectal temperature (RT) and heat tolerance coefficient (HTC) were measured for both NHS and HS groups according to a method previously described, HTC = 100-10 (RT-38.3) [18]. Four cows were designated as biological replicates for NHS and HS. In each case, one complete mammary gland was removed after excision of the intramammary lymph node. The samples were frozen in liquid nitrogen until RNA extraction.

Total RNA extraction
Total RNA from the mammary gland was isolated using Trizol Reagent (Life Technologies, USA) according to manufacturer's instructions. RQ1 DNase (Promega, USA) was used to treat the total RNA to digest the DNA. The concentration and purity of the extracted RNA were determined with a NanoDrop 1000 Spectrophotometer (NanoDrop Technologies, USA). The quality of RNA was assessed through the Agilent 2100 Bioanalyzer (Agilent Technologies, USA). All eight RNA samples had an RNA integrity number (RIN) of ≥8 and were stored at − 80°C for further analysis.

Small RNA library construction
Equal amounts of total RNA (330 ng) from each mammary gland sample of both NHS and HS groups were used to construct a miRNA library through a TruSeq Small RNA Sample Preparation kit (Illumina, USA) following the manufacturer's protocol. PCR amplification was performed including 13 cycles. Eight small RNA libraries were constructed in equal amounts for gel purification. Quality and quantity of purified small RNA were estimated using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Sequencing was carried out on a HiSeq-2000 system (Illumina, USA). Real-Time Analysis (RTA) and base calling were performed by HiSeq Control Software Version 1.4.8 (Illumina, USA).

Small RNA sequence analysis
Low-quality reads were removed from the raw data. After trimming the 3′ adaptor sequence, small RNA reads with the length of 18-30 nt from all libraries were extracted using miRDeep2 (version 2.0.0.5) with the default parameters to identify known and novel miRNAs [19]. Each library was processed separately. Subsequently, all reads were mapped to the bovine genome with Bowtie [20]. The reads that mapped to bovine tRNAs, rRNAs, and snoRNAs in the Rfam RNA family database [21] were discarded. Small RNAs that only mapped to genomic repeat loci were removed. Novel miRNA and precursors were identified by the core module miRDeep2.pl. Datasets of novel miRNA and precursors were created through adding miRNA predicted with a miRDeep2 score > 0 to known miRNA. The expression of detected miRNAs for each library was estimated by the Quantifier module of miRDeep2 [19]. To investigate the regulation of miRNAs in mammary gland in response to heat stress, miRNA expression under heat stress (summer) were compared to expression in controls (spring) using DEseq [22].

Target gene prediction and pathway analysis
The target genes of differentially expressed miRNAs were predicted through TargetScan (http://www.targetscan.org) and miRanda (http://www.microrna.org/microrna/home.do) according to a previous study [23]. To comprehensively describe the properties of genes and gene products, we executed gene ontology (GO) annotation and enrichment analysis from three ontologies: molecular function, cellular component and biological process. Functional analysis of predicted gene targets was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID version 6.7, http://david.abcc.ncifcrf.gov/) for pathway analysis.

Stem-loop quantitative PCR and data analysis
Seven differentially expressed miRNAs (bta-miR-21-5, bta-miR-99a-5p, bta-miR-146b, bta-miR-145, bta-miR-2285 t, bta-miR-133a, and bta-miR-29c) were validated through a standardized and reliable stem-loop qRT-PCR procedure. The RT reaction mixture contained a 4 μg aliquot of total RNA and a mixture of 1 μL of each RT primers (10 μM) for all of the mature miRNAs and U6, which was chosen as a reference control (Additional file 1), 10 μL 2 × reaction mix, 1.5 μL RT enzyme mix (Sangon, China). The mixture was incubated at 16°C for 30 min, 37°C for 30 min, 85°C for 5 min, subsequently, frozen on ice for at least 5 min. SYBR® Green Realtime PCR Master Mix (Sangon, China) was used to detect miRNA expression by a Bio-Rad IQ5 System (USA). Briefly, cDNAs were diluted 10 times and a 1 μL diluted sample was used as a template in a 20 μL PCR reaction, which contained 10 μL 2 × SYBR Green Realtime PCR Master Mix, 0.25 μM of a miRNA-specific forward primer and universal reverse primer, respectively. The quantitative PCR was conducted in triplicate for 1 min at 95°C, followed by 40 cycles of 15 s at 95°C and 30 s at 60°C. For each PCR, dissociation curve analysis was carried out to discriminate specific products from primer dimers. The fold changes of miRNA in different samples were calculated by ΔCt method.
The relative expression of target miRNA was calculated by the following formula: ΔCt (target miRNA) = Ct (target miRNA)-Ct (internal reference). U6 was used as an internal reference. The NHS group at spring was selected as a calibrator for the relative quantification. The relative expression of miRNA normalized to internal control and relative to the calibrator was calculated as follows: Relative expression of target miRNA = 2 -ΔΔCt , where ΔΔCt = ΔCt (target miRNA, sample) -ΔCt (target miRNA, calibrator). The difference between the two groups was compared using T-test (SAS version 9.2, 2008). Statistical difference was at P < 0.05.

Characterization of the stress response
Rectal temperature and HTC at two different environmental temperature ranges are presented in Table 1. A THI of above 72 and a rectal temperature of above 39°C were regarded as an indication of heat stress. The physiological parameters of HS were significantly different from the NHS group of cows (P < 0.05). The average rectal temperature of HS group was higher, and a significant decrease of HTC as compared with the control. Rectal temperature and HTC data showed that the cattle were under heat stress.

miRNA sequencing
Eight small RNA libraries were constructed and sequenced. A total of 20,818,666 high-quality reads were generated. Among them, 19,859,218 sequences ranging from 18 to 30 nucleotides were obtained after adaptor trimming, accounting for 95.39% of all small RNA (sRNA) sequences. A summary of the data is provided in Table 2. Alignment with miRBase (Release 21) indicated that miRNAs were highly enriched in all libraries. To assess the efficiency of high-throughput sequencing for sRNA detection, all sequence reads were annotated and classified by alignment with Rfam databases. Of the 18 to 30 nucleotide sRNA fraction, four out of five (86.72%) were identified as other noncoding sRNA, while only a small number (13.28%) aligned to rRNAs, tRNAs, snRNAs and snoRNAs in bovine (Fig. 1). The major reads from sRNAs were 20 to 24 nucleotides in length (accounting for 91.14% of total number). Dominant reads of sRNAs were 22 nucleotides in length (41.38%), followed by 21 and 23 nucleotides, and lastly 20 and 24 nucleotides (Fig. 2).

Identification of novel miRNAs and miRNA candidates
To determine whether these small RNA sequences are genuine bovine miRNAs, we scanned the bovine genome for hairpin structures comprising the candidate miRNA with miRDeep2 software (version 2.0.0.5), which can be used to identify both known and novel miRNAs from deep sequenced sRNA libraries. In total, 483 loci possessed the typical stem-loop structures matching the known miRNA hairpins (808 miRNA, miRbase 21) in the mammary gland. A total of 139 novel miRNA hairpins were identified (Additional file 2).
Analyses of the first nucleotide bias of the 18-30 nt miRNAs candidates revealed that uridine (U) was the most common at the 5′ end of the 20-24 nt miRNAs in the mammary gland (Fig. 3). Moreover, miRNA nucleotide bias at each position also showed that U was the dominant first nucleotide (Fig. 4).

Effects of heat stress on miRNA expression in mammary gland
Given that miRNAs play important roles in many biological processes, we speculated that the expression of miRNAs might be regulated in the mammary glands of Holstein cattle under heat stress. The global expression of miRNAs under normal condition and heat-stressed condition were profiled and the correlations between libraries were performed with the normalized counts of all detected miRNAs. Only a small number of miRNAs were significantly regulated in mammary glands of Holstein cattle under heat stress. DEseq, an R/Bioconductor package [22] method was used to analyze differentially expressed miRNAs between different conditions based on sequence counts. We observed that 24 miRNAs were differentially expressed (P < 0.05) in the heat stress condition when compared to the control (Additional file 3). The highest differentially expressed miRNAs (bta-miR-21-5p, bta-miR-99a-5p, bta-miR-146b, bta-miR-145, bta-miR-2285 t, bta-miR-133a, and bta-miR-29c) are shown in Fig. 5. The expression of bta-miR-145, bta-miR-2285 t, bta-miR-133a, and bta-miR-29c was increased and reached a 3.50, 3.24, 4.30, and 4.03-fold increase (P < 0.05), respectively. Conversely, the expression of bta-miR-21-5p, bta-miR-99a-5p, and bta-miR-146b were decreased under heat stress condition (P < 0.05). Furthermore, we also verified the expression of these miRNAs using qRT-PCR (P < 0.05) and the expression pattern was consistent with sequencing results (Additional file 4).
Biological function enrichment analysis of predicted target genes of differentially expressed miRNAs in mammary glands of Holstein cattle under heat stress Target gene prediction of 24 miRNAs differentially regulated by heat stress indicates that about 26,824 genes that were enriched in 59 GO terms may be regulated by these miRNAs (Fig. 6). GO analysis showed that the predicted targets of differentially expressed miRNAs were significantly enriched (P < 0.05) in different functional groups, namely, biological regulation, cellular process, metabolic process, multicellular organismal process, regulation of biological process, response to stimulus and single-organism process (Fig. 6, Table 4). One of the GO terms included more than 2000 responses to stimulus-related genes. Interestingly, several KEGG (Kyoto encyclopedia of genes and genomes) pathways were significantly enriched (P < 0.05) by target genes of the 24 differentially expressed miRNAs (Table 5). Notably, pathways of the RNA degradation, mTOR signaling pathway, immune system and pathways in human diseases, especially cancer were significantly enriched by target genes. This indicated that the differentially expressed  Total reads, total sequenced reads; High quality reads, number of high quality reads with no N, no more than 4 bases whose quality score is lower than 10 and no more than 6 bases whose quality score is lower than 13; 3′ adapter null, number of reads with no 3'adaptor; Insert null, number of reads with no insertion; 5′ adapter contaminants, number of 5′ contaminants; Clean reads, number of clean reads after adaptors and contaminants are removed which are used in the following analysis miRNAs identified in this study might act as a dominant regulator during heat stress.

Discussion
Heat stress triggers a dramatic and complex program of altered gene expression in mammary glands similar to patterns investigated in other cell types exposed to thermal stress. As reported by Sonna et al. [24], these changes include inhibition of DNA synthesis, RNA transcription, and translation, disruption of cytoskeletal components, and alterations in metabolism. In the present study, thermal stress induced changes in the miRNA expression in dairy cattle mammary glands. Although the miRNA expression profiles in dairy cattle have been studied, there is still a limited understanding of their role in heat-stressed bovine mammary glands and the influence of this on dairy potential. In this study, 483 known miRNAs were detected in dairy cattle mammary glands by RNA-Seq. Several studies have uncovered the roles of the highly abundant miRNAs in mammary glands. For example, two of the eight miRNAs with the highest expression in the dairy cattle mammary gland in this study, bta-miR-21-5p, the other arm of bta-miR-21, and bta-miR-143, were shown to promote adipogenesis [25,26], which accounts for 12.70% of the miRNAs expressed in the dairy cattle mammary gland. The expression of miR-21-5p was strongly induced at 7 d postpartum compared with the dry period suggesting that it might promote mammary cell proliferation during early lactation [27]. MiR-21 decreases the expression of Tgfbr2 by targeting TGF b receptor II (Tgfbr2) and eventually enhances adipogenic differentiation [26]. Moreover, miR-21 is associated with thermal stress in Frieswal crossbred dairy cattle [15]. Our data showed that the expression of bta-miR-21-5p was lower in heat stressed group than that in control, which is consistent with bta-miR-21 expression in PBMC of Sahiwal cows [16]. The up-regulation of miR-143 decreased the expression of pleiotrophin and increased some adipocyte-important genes, enhancing the rate of adipocyte differentiation at early stages of adipogenesis [25]. Decreased expression of miR-143 by its antisense sequence suppressed differentiation of preadipocytes through repressing ERK5, suggesting this miRNA may play a key role in adipocyte differentiation [28].  One of the other eight most abundant miRNAs in this study is miR-148a, orthologs of miR-148. MiR-148 was highly abundant in the lactating mammary gland of mouse and goat [29,30]. MiR-148a and miR-26a show the similar expression patterns during the lactation period in cow milk [31]. Accumulating evidence indicates that miR-148a induce cell proliferation and differentiation [32,33]. MiR-148a promotes adipogenesis by repressing Wnt signaling [34]. In addition, the expression of bta-miR-145, bta-miR-133a, and bta-miR-29c was increased in the heatstressed group. MiR-145 is associated with cow mastitis caused by Staphylococcus aureus [12]. Moreover, the expression of miR-145 during differentiation can regulate the insulin receptor substrate 1 to inhibit adipogenesis [35]. Inhibition of miR-133 led to the expression of GLUT4 and insulin-mediated glucose uptake attenuation in cardiomyocytes [36]. Overexpression of miR-29a-c reduced the protein levels of PGC-1a and G6 Pase in primary hepatocytes and mouse livers [37]. These researches indicate that the highly expressed miRNAs may be related to the mammary gland biology, milk synthesis, and lactation process.
Bta-miR-2285 t was significantly increased in the mammary tissue of dairy cattle (Holstein-Friesian) compared with beef (Limousin) postpubertal heifers [38]. However, its role in the mammary gland remains unclear. The expression of bta-miR-2285 t and bta-miR-146b was decreased in the heat-stressed group compared with the control. Similarly, miR-146b has a lower expression in the serum of heatstressed Holstein cows [14]. It has been suggested that miR-146b can modulate Sirtuin1, suppressing the negative regulators of adipogenesis, and eventually promoting adipogenesis [39]. Up-regulation of miR-146b is found during pregnancy, especially in the luminal progenitors compared to the basal/stem cells, suggesting it is involved in the differentiation of mammary epithelial cells [40]. Furthermore, the expression of miR-146b was upregulated in the luminal progenitors in pregnant mice, which indicates that miR-146b is involved in the differentiation of the mammary stem cells [40].
Interestingly, we found that the expression of bta-miR-21-5p and bta-miR-146b tended to decrease, and bta-miR-145 tended to increase in the heat-stressed group   compared with control group. These results are consistent with the phenomenon that, acute heat stress affects the lipolysis and the rate-limiting enzyme of lipogenesis in bovine adipocytes [41]. Considering that adipocytes in the mammary gland can regulate the growth and biological function of the mammary epithelium [42], we speculated that these miRNAs might have a role in the regulation of milk fat synthesis. A total of 139 novel miRNA hairpins were detected. High-throughput sequencing and bioinformatics analysis have become the main methods to identify the potentially novel miRNAs since 2008 [43]. Given the number of known miRNAs (808 miRNAs) in the miRBase (Release 21), our data may enrich the miRNA resources in bovine.
To further understand and provide some molecular insight into the physiological functions and biological processes involving these miRNAs in heat-stressed mammary glands, target genes were predicted based on miRNA/ mRNA interactions. The predicted target genes were classified using KEGG function annotations to character the pathways that were actively regulated by miRNAs in the mammary gland. In the present study, pathway analysis of miRNA targets revealed that Wnt, TGF-β, MAPK, Notch, and JAK-STAT signaling pathways may play key roles in the mammary gland in the process of heat stress. Wnt, EGFR, TGF-β, and insulin signaling pathways are known to play a key role in normal development of the mammary gland [44][45][46][47]. The MAPK pathway is an important regulator of mammary epithelial cell differentiation and function [48]. These reports are supported by our findings.  In addition, MAPK signaling, PI3k-Akt signaling, and immune-regulatory are greatly influenced by miRNAmediated regulation in Frieswal cattle [15]. The function analysis of differentially expressed miRNAs and their target genes suggested the effects of heat stress on signaling mechanisms. Although many target gene candidates were predicted by bioinformatics methods, structural verification and signaling pathways analysis in vitro need to be further performed to validate the relationship between miRNAs and mRNA. Target genes with pathway annotation, The number of target genes which are associated with the pathway (ratio). All genes of the bos taurus with pathway annotation, The number of all reference genes which are associated with the pathway (ratio)

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author details