Transcriptome microRNA profiling of bovine mammary epithelial cells challenged with Escherichia coli or Staphylococcus aureus bacteria reveals pathogen directed microRNA expression profiles

Background MicroRNAs (miRNAs) can post-transcriptionally regulate gene expression and have been shown to be critical regulators to the fine-tuning of epithelial immune responses. However, the role of miRNAs in bovine responses to E. coli and S. aureus, two mastitis causing pathogens, is not well understood. Results The global expression of miRNAs in bovine mammary epithelial cells (MAC-T cells) challenged with and without heat-inactivated Staphylococcus aureus (S. aureus) or Escherichia coli (E. coli) bacteria at 0, 6, 12, 24, and 48 hr was profiled using RNA-Seq. A total of 231 known bovine miRNAs were identified with more than 10 counts per million in at least one of 13 libraries and 5 miRNAs including bta-miR-21-5p, miR-27b, miR-22-3p, miR-184 and let-7f represented more than 50% of the abundance. One hundred and thirteen novel miRNAs were also identified and more than one third of them belong to the bta-miR-2284 family. Seventeen miRNAs were significantly (P < 0.05) differentially regulated by the presence of pathogens. E. coli initiated an earlier regulation of miRNAs (6 miRNAs differentially regulated within the first 6 hrs post challenge as compared to 1 miRNA for S. aureus) while S. aureus presented a delayed response. Five differentially expressed miRNAs (bta-miR-184, miR-24-3p, miR-148, miR-486 and let-7a-5p) were unique to E. coli while four (bta-miR-2339, miR-499, miR-23a and miR-99b) were unique to S. aureus. In addition, our study revealed a temporal differential regulation of five miRNAs (bta-miR-193a-3p, miR-423-5p, miR-30b-5p, miR-29c and miR-un116) in unchallenged cells. Target gene predictions of pathogen differentially expressed miRNAs indicate a significant enrichment in gene ontology functional categories in development/cellular processes, biological regulation as well as cell growth and death. Furthermore, target genes were significantly enriched in several KEGG pathways including immune system, signal transduction, cellular process, nervous system, development and human diseases. Conclusion Using next-generation sequencing, our study identified a pathogen directed differential regulation of miRNAs in MAC-T cells with roles in immunity and development. Our study provides a further confirmation of the involvement of mammary epithelia cells in contributing to the immune response to infecting pathogens and suggests the potential of miRNAs to serve as biomarkers for diagnosis and development of control measures.


Background
Mastitis is an infectious disease of the mammary gland that has been considered as one of the most economically important diseases in the dairy industry for several decades [1]. Mastitis can be caused by a variety of bacteria including Gram-negative Escherichia coli (E. coli) and Gram-positive Staphyloccucus aureus (S. aureus) [2]. The mastitis caused by E. coli and other Gram-negative bacteria is often clinical with an acute and severe inflammation and the pathogens may be eventually cleared by the immune system within days or with antibiotic treatment [3,4]. In contrast, infection with Gram-positive bacteria like S. aureus often causes mild mastitis and the clearance of the pathogens by antibiotics is often ineffective [5].
Efforts at mastitis control include understanding host response mechanisms to infecting pathogens and development of appropriate control strategies. Transcriptomic as well as proteomic profilings have shown a marked difference in the response of host to E. coli and S. aureus bacteria. Transcriptomic studies on bovine mammary epithelial cells in vitro challenged with mastitis pathogens and/or mammary gland tissues collected after intramammary infection revealed very different mechanisms in host innate immune responses to pathogens [6][7][8]. Differential cytokines and chemokines (at protein level) and other immune response proteins were also observed in bovine mammary epithelial cells and tissues or milk, in response to E. coli and S. aureus [4,9,10].
Recently, an endogenous noncoding small RNA known as microRNA (miRNA) has been shown to be involved in a wide variety of biological processes such as development, differentiation, apoptosis and viral infection. The role of miRNA in the regulation of the innate and acquired immune response has been well reviewed [11][12][13][14]. It was shown that miR-155 promotes the production of TNF-α in human embryonic kidney cells (HEK-293), suggesting the positive role of miR-155 to modulate the release of inflammatory mediators [15]. Recent studies on human epithelial cells infected by a protozoan parasite (Cryptosporidium parvum) revealed the up-regulation of let-7i or miR-27b. While let-7i directly regulates the expression of toll-like receptor 4 (TLR4) [16], miR-27b targets KH-type splicing regulatory protein (KSRP) to coordinate TLR4-mediated epithelial immune responses to pathogens [17], indicating that miRNA may act as a critical regulator in epithelial immune responses. In addition, it was suggested that miRNA may modulate epithelial immune responses at every step of the innate immune network, including production and release of cytokines/ chemokines, expression of adhesion and costimulatory molecules [18]. It has been clearly demonstrated that bovine mammary epithelial cells mount a robust immune response to the presence of bacterial pathogens [19][20][21].
However, how miRNAs modulate this robust immune response by the bovine mammary epithelial cells to the presence of the mastitis pathogens, E. coli and S. aureus bacteria, is poorly understood. Furthermore, studies elucidating the regulatory roles of miRNA in bovine immunity and infection are few. A recent study using quantitative real time PCR technique revealed differential expression of five inflammation related miRNAs (miR-9, miR-125b, miR-155, miR-146a and miR-223) after stimulation of bovine monocytes with lipopolysaccharide (LPS) and S. aureus enterotoxin B [22]. Using the same technique, four miRNAs (bta-miR-181a, miR-16, miR-31 and miR-223) out of 14 miRNAs associated with regulation of innate immunity and mammary cell function were shown to be differentially regulated in bovine mammary tissue challenged with Streptococcus uberis (S. uberis) [23]. These studies suggest roles of miRNAs in bovine mammary gland immunity but the extent of information is limited by the technique used. Increasingly, more robust techniques like microarray and next-generation sequencing are finding applications in the elucidation of miRNAs in immunity [24,25]. In particular, next-generation sequencing has the ability to profile the expression of both known and novel miRNAs with high resolution and accuracy and to distinguish miRNAs that are very similar in sequence as well as isomiRs [26]. Only two studies to date have used next-generation RNA sequencing to study the regulatory roles of miRNAs upon bacterial or viral infections in bovine [25,27]. Using RNA sequencing, 21 miRNAs were detected as significantly differentially expressed upon infection of bovine primary epithelial cells with S. uberis, as well as a unique miRNA profile in response to a Gram-positive bacterial infection [25]. However, more information is needed to understand the roles of miRNAs in modulating bovine mammary gland infections for their effective application as biomarkers of mastitis or as therapeutic agents.
To investigate the role of miRNA in host defense to two mastitis pathogens, miRNA expression in bovine mammary epithelial cells (MAC-T) challenged with heat-inactivated Gram-negative E. coli stain P4 or Gram-positive S. aureus strain Smith CP bacteria were characterized by nextgeneration sequencing at different time points (6, 12, 24 or 48 hr) following challenge and at 0, 6, 12, 24 or 48 hr without challenge. The next-generation sequencing technology allowed simultaneous detection of known and novel bovine miRNAs, and global miRNAs expression as well as pathogen directed differential miRNA expression patterns.

miRNA sequencing
Thirteen small RNA libraries were constructed and sequenced simultaneously. A total of 15,315,312 high-quality reads were generated. Among them, 12,272,208 sequences ranging from 18 to 30 nucleotides were obtained after adaptor trimming, accounting for 80.1% of all small RNA (sRNA) sequences. Alignment with miRBase (Release 19) revealed that miRNAs were highly enriched in all libraries. Of the 18 to 30 nucleotide sRNA fraction, more than three quarters (76.02%) of them were identified as known bovine miRNAs, while only a small number (< 2%) aligned to bovine tRNAs, rRNAs and snoRNAs. The remaining reads were other sRNAs including novel bovine miRNAs (~0.2%), loop sequences of miRNA precursors and sequencing artifacts ( Figure 1A). The majority of reads from sRNAs and known miRNAs were 20 to 24 nucleotides in length (comprising 89.4% and 96.6% of their total number, respectively). Dominant reads of sRNAs or known miRNAs were 22 nucleotides in length (~50%), followed by 23, 24 or 20 nucleotides, and lastly 21 nucleotides ( Figure 1B, C).

Identification of novel miRNAs and miRNA candidates
With the high-throughput sequencing data, miRDeep2 [28] generated a score for each known miRNA and novel miRNA. The score of 5.0 yielded a signal-to-noise ratio of 16.2 [28] and was used as cut-off for novel miRNA prediction (Additional file 1: Table S1.1). A total of 114 novel miRNA hairpins were identified (Additional file 2: Table S2). Moreover, 186 hairpins with lower scores (0-4.9) but having more than 10 total counts were further analyzed and 5 or 153 of them were defined as precursors of novel miRNAs or miRNA candidates, respectively (Additional file 3: Table S3). Although almost all miRNAs have a variety of isoforms, most of the known miRNAs detected or undetected by miRDeep2 (Additional file 1: Table S1.2 and Additional file 1: Table S1.3) were within the size range of 20 to 24 nucleotides. Moreover, a false precursor of bta-miR-532 predicted by miRDeep2 (Additional file 1: Table S1.2) contained a star miRNA of 17 nucleotides in length, indicating that a precursor generating a short miRNA may have a high false positive rate. Therefore, precursors which form mature miRNAs < 20 nucleotides in length were not further analyzed.
In total, five homologs of known miRNA precursors from other animal species and ten reversely complementary to known bovine miRNA precursors were identified to generate 24 novel miRNAs ( Table 2, Additional file 2:  Table S2 and Additional file 3: Table S3). We also found an ortholog of bta-miR-222 (bta-miR-222b, Table 2, Additional file 3: Table S3). In addition, more than one third of novel or candidate miRNA precursors (41/114, 65/153, respectively) were found to be clustered together and have a high similarity to those of bta-mir-2284 family in sequences (Additional file 2: Table S2,  Additional file 3: Table S3 and Additional file 4: Figure S1). The remaining 62 novel miRNA precursors were found to produce 56 mature miRNAs which are less conserved to any known miRNA in miRBase (Release 19, Additional file 2: Table S2).

Time effect on miRNA expression in MAC-T cells
Considering that miRNAs play roles in almost all biological processes, we hypothesized that the expression of miRNAs   Figure S2), suggesting that only a limited number of miRNAs might be significantly regulated during the 48 hrs of cell challenge with mastitis pathogens. The lowest correlation was found between 0 hr and 48 hrs in control samples (R 2 = 0.961), indicating that the expression of miRNA might be significantly regulated over time. Therefore, differentially expressed miRNAs between different time points and 0 hr in control cells were analyzed using DEseq, an R/Bioconductor package [29] method to compare differential expression based on sequence counts. We observed that five miRNAs (bta-miR-193a-3p, miR-423-5p, miR-30b-5p, miR-29c and miR-un116) were differentially expressed (P < 0.05) in at least two time points in control cells as compared to 0 hr ( Figure 2). The expression of bta-miR-193a-3p was the most regulated over time, increasing lineally and reaching 8.03 fold increase (P < 0.0001) within 48 hrs of cell growth. The expression of bta-miR-30b-5p and bta-miR-29c also increased over time but to a lesser magnitude as compared to bta-miR193a-3p. The expression of bta-miR-423-5p was highest at 6 hrs (P < 0.001) and returned to 0 hr level by 24 hrs (Figure 2). Conversely, the expression of bta-miR-423-5p and bta-miR-un116 decreased with time. The expression of bta-miR-193a-3p and bta-miR-423-5p in control cells was further validated by qRT-PCR. Both miRNAs were significantly up-regulated and reached a peak at 12 hrs, although no change larger than 2-fold was found between any two time points (Additional file 6: Figure S3).  (Table 3) by the presence of pathogenic bacteria in MAC-T cells. Within 6 hrs of the presence of E. coli, the expression of 6 miRNAs in MAC-T cells was significantly altered (P < 0.05), three were down regulated (bta-miR-193a-3p, miR-30c and miR-30b-5p) while three were up-regulated (bta-miR-365-3p, miR-184 and miR-24-3p) ( Table 3). With time, the number of differentially expressed miRNAs in the presence of E. coli  Sequence here only shows the dominant one and the star sequence that was absent was predicted by miRDeep2. 5 miRNA with score less than 5 is identified from miRNA candidate precursor.    Table S4). Gene ontology (GO) functional annotation showed that the target genes of differentially regulated miRNAs were significantly (P < 0.05) enriched in different functional groups, namely; cellular process, developmental process, localization, biological regulation, cellular component organization, cell death, multicellular organismal process, metabolic process, establishment of localization and cell growth ( Table 5). The  highest involvement of target genes was of 14 miRNAs in the cellular and developmental process followed by target genes of 11 miRNAs in localization and of 10 miRNAs in biological regulation. Interestingly, several KEGG pathways were significantly (P < 0.05) enriched by target genes of 14 of the 17 differentially expressed miRNAs (Table 6 and Additional file 8: Table S5). Notably, pathways of the immune system, signal transduction, cellular processes, nervous system, development and pathways in human diseases, especially cancer were significantly enriched by target genes of at least 4 miRNAs ( Table 6).

Discussion
Next-generation sequencing technologies have not only enabled the detection of more lowly expressed miRNAs and novel miRNA discovery, but also accurate profiling of the expression of miRNAs at a high-throughput level.
Using deep-sequencing, we have identified a number of novel bovine miRNAs and revealed a differential regulation of miRNAs in MAC-T cells challenged with E. coli or S. aureus bacteria, thus suggesting the regulation of miRNAs in host response to these mastitis pathogens. The number of miRNAs deposited in miRBase has increased exponentially in recent years and high throughput sequencing has contributed to almost all sequences submitted to miRBase after 2007 [30]. However, high throughput sequencing generates millions of reads per sample and requires powerful and accurate computational methods to mine known and novel miRNAs. Several computational software packages such as miRDeep [28], MIReNA [31], miRanalyzer [32], and miRTRAP [33] have been developed to identify miRNAs from high-throughput data. Our study used a combination of miRDeep2 and Blastn and identified 113 novel miRNAs in MAC-T cells. It is notable that computational based approaches may over-predict novel miRNAs and our data remain to be validated by further studies.
It is necessary to give a consistent and unique gene name for each novel miRNA. The nomenclature of miRNA has been well described previously [30,34,35]. However, the annotation of miRNAs still presents misperceptions in some instances. Some recently identified miRNAs and their precursors are reversely complementary to previously known miRNAs and precursors. For example, hsa-mir-103a and hsa-mir-103b (miRBase Release 19) are paired (reversely complementary) but are usually thought to be distinguished by 1 or 2 nucleotide differences in their nucleotide sequences. Here we added a symbol "`" after miRNA to represent the antisense of a known miRNA precursor (Table 2), i.e. hsa-mir-103 and mir-103′instead of mir-103a, b, respectively. In addition and to the best of our knowledge, there is no established criterion that defines conserved miRNAs. We suggest that a novel miRNA is homologous or orthologous if: 1) it precursor aligns to a known miRNA precursor with > 70% identity and > 75% sequence coverage; or 2) it has > 90% identity and > 90% length coverage with a known miRNA as well as same arm location on the precursor. Previously, only sequence similarity was taken into account for the annotation of conserved miRNAs and disparities have been known to occur. For example, bta-miR-193a [36] (miRBase Release 19), generated from bta-mir-193a-2, is reversely complementary to bta-miR-193a-3p from 1 to 19 nucleotides but not correctly named. Using this criterion, one miRNA candidate identified in this study (un219, Additional file 3: Table S3) was aligned to bta-miR-222. Un219 is at the same arm as bta-miR-222 on the precursor and therefore named as bta-miR-222b ( Table 2, Additional file 3: Table S3).
Only a few miRNAs (Table 1) were highly expressed, accounting for 78.55% of all known bovine miRNAs detected. This observation is consistent with a recent study that detected similar miRNAs expression levels in bovine primary mammary epithelial cells infected with Streptococcus uberis and also, seven of the top expressed miRNAs in that study [25] are the same as in this study. Most of the detected top expressed miRNAs are conserved in human, mouse, and bovine, and belong to several miRNA families, vis, miR-31, miR-26a, miR-27a-3p/27b, let-7a-5p/7f/7i, miR-21-5p, miR-22-3p, miR-184, miR-186, miR-191, miR-205 and miR221/222. In particular, several of these highly expressed miRNAs are known to play roles in growth and immunity. The top expressed miRNA, bta-miR-21-5p in this study, was also the most highly expressed miRNA in human blood monocytes challenged with Mycobacterium leprae and also found to negatively regulate the vitamin D-dependent antimicrobial pathway in leprosy [37]. miR-21-5p has been shown to  and to promote the expression of interleukin 10 (IL-10) [38]. Further evidence of the immune capacity of miR-21 was from the study of Narducci et al. [39] who provided in vitro evidence for involvement of miR-21, miR-214 and miR-486 in cell survival in Sézary syndrome in humans. It has been recently shown that the level of immune and development related miRNAs, including miR-27b, were significantly higher in colostrum than in mature milk [40]. Up-regulation of miR-27b by LPS was found to destabilize proliferator-activated receptor γ1 (PPARγ1) mRNA which is often associated with chronic inflammatory diseases [41]. In addition, miR-27b was reported to target KSRP and increase iNOS mRNA stability for host's defense against cryptosporidial infection [17], further supporting a role for this miRNA in immunity. Therefore, the highly expressed miR-21-5p and miR-27b in mammary epithelial cells in this study might be associated with a function in immunity. Five miRNAs were found to be significantly regulated during the 48 hr cell incubation without pathogens. It is not surprising owing to their involvement in almost all biological processes as demonstrated by GO functional annotation of the target genes of three (bta-miR-193a-3p, miR-423-5p and miR-30b-5p) of these miRNAs. GO functional annotation of target genes of bta-miR-193a-3p and miR-30b-5p showed enriched genes related to cell growth and death, e.g. growth arrest specific gene 1 (GAS1), myeloid cell factor 1 (MCL1, also BCL2 related), BCL2-like 11 (apoptosis facilitator) and programmed cell death 10 (PDCD10). GAS1 has been shown to be involved in growth suppression through blocking of entry to S-phase and also prevents cycling of normal and transformed cells [42,43]. Also, involvement of MCL1 and PDCD10 in the regulation of apoptosis and cell survival has been demonstrated [44][45][46]. It is well known that cells enter the stationary phase after reaching confluence. Additionally, we used fetal bovine serum (FBS) free media, known to inhibit proliferation, during the 48 hr period of cell incubation. The up-regulation of miR-193a-3p and miR-30b-5p may play regulatory roles in cell death which need to be further confirmed in the context of mastitis.
Out of the 231 known and 113 novel miRNAs identified in this study, only 17 showed differential expression after challenge of MAC-T cells with E. coli or S. aureus bacteria, which may look small considering the large number of identified miRNAs. It is probable that observed miRNA expression profiles in response to heat inactivation of bacteria in this study may be different than expression elicited by live bacterial cells. In a recent study [47], a substantially different miRNA expression profile emerged from killed virulent bacilli as compared to live active bacteria thus suggesting an active influence of living bacteria on cellular miRNA metabolism. However, our results on response to heat inactivated S. aureus are in line with literature information whereby 21 miRNAs were shown to be differentially expressed after challenge of bovine primary mammary epithelia cells with live S. uberus (Gram-positive) bacteria [25]. Regardless, further studies on the miRNA expression patterns during infection with live and active bacteria is necessary to verify the miRNA expression variation identified between these two species. In addition, there were discrepancies in expression pattern of some low abundant miRNAs during the infection detected by sequencing as compared to RT-PCR. This difference may be due to the pooling approach used for library preparation and subsequent sequencing, as well as the relative low sequence reads which may limit the proper detection of the low abundant miRNAs by miRNA-Seq. Furthermore, the miRNA signature in different organs, tissues or cell types is different and also regulated differently. miRNA expression profile has been shown to vary between mature milk and colostrum [48], between lactating and non-lactating mammary glands [49] and, between organs (lung, brain, liver and spleen) and whole blood [50]. During the onset of mastitis and subsequent progression, the whole system reacts by recruiting immune cells into the mammary gland to help combat infection. These immune cells and other factors, including resident factors and cells in the mammary gland work in concert to ward off infection. It should be noted here that, the mammary epithelial cells used in this study only constitutes a small portion of that army of cell types that secrete defence factors against invading pathogens. Thus, the present findings portray only the role of epithelial cells derived miRNAs in the fight against mastitis pathogens, which may be different from the holistic picture that can be obtained after examination of all players. In profiling the transcriptome of S. uberis induced mastitis, fundamental differences were revealed between immune gene expression in the mammary gland and in a primary cell culture model thus demonstrating the complexity of the bovine mammary gland immune response to an infecting pathogen, and also indicating that a coordinated response exists between resident, recruited, and inducible immune factors [51]. Furthermore, transcriptomics comparisons between MAC-T cells and mammary tissue during late pregnancy and peak lactation showed a larger overall similarity in gene expression between MAC-T cells and lactating than non-lactating mammary tissue [52]. Despite these demonstrated differences [51,52] and or similarities [52], our results showed similarity between highly expressed miRNA profiles with primary bovine mammary epithelial cells challenged with Gram positive bacteria [25] and nonchallenged mammary gland tissues (data not shown) thus indicating that the MAC-T cell is still a valuable tool to study mammary gland biology. Importantly, the findings from our study show the positive involvement of mammary epithelia cells in contributing to the immune response to infecting pathogens.
Furthermore, our study has demonstrated that E. coli initiated a stronger response at the early hours of infection as demonstrated by the differential response of six miRNAs within the first six hours of cell challenge with E. coli bacteria as opposed to only one differentially expressed miRNA in the presence of S. aureus. By 48 hrs, two miRNAs were differentially regulated in E. coli challenged cells while in S. aureus cells, more miRNAs were implicated with time. Furthermore, three of the miRNAs that responded significantly to the presence of E. coli at 6 hrs post infection only showed a differential response in S. aureus cells by 24 or 48 hrs. Interestingly, our study shows that a different set of five miRNAs (miR-184, miR-24-3p, miR-148, miR-486 and let-7a-5p) were unique to E. coli bacteria while another set of four (miR-2339, miR-499, miR-23a and miR-99b) were unique to S. aureus bacteria. Our data therefore demonstrate a differential response pattern of MAC-T cells to Gram-negative E. coli, as compared to Gram-positive S. aureus. Mastitis caused by pathogenic ligands from Gram-positive or Gram-negative bacteria usually elicits different response patterns from the host thus leading to different degrees of mastitis. Our data supports the intense immune reaction usually caused by E. coli mastitis, which often leads to acute mastitis. On the contrary, the slower initial response of miRNAs to S. aureus bacteria may support the slow progression of mastitis caused by this type of bacteria. Differential response patterns of miRNAs to Gram-positive and Gramnegative bacteria have been demonstrated previously [22,49]. Furthermore, our data is supported by previous studies that have demonstrated pathogen-dependent differences of host (mammary gland) responses to other immune factors (genes and proteins) [4,[7][8][9][10].
It is notable that predicted gene targets of most of the differentially expressed miRNAs in response to pathogens in this study are significantly enriched for GO functionally annotated roles like developmental process, cellular process, localization, biological regulation, cellular component organization, multicellular organismal process, cell death, metabolic process, establishment of localization and growth, thus supporting the involvement of these miRNAs in numerous developmental and physiological processes, including disease development. This also goes a long way to portray the significant roles of miRNAs in supporting the mammary gland in its growth and productive capacities as well as defense capabilities. Interestingly, pathway analysis of miRNA predicted gene targets, showed significant enrichments in KEGG pathways already associated with host response to different diseases. Notably, many of the KEGG pathways significantly enriched by differentially expressed miRNA target genes in this study have been previously associated to mammary gland responses to bacterial pathogens in bovine [6,23,25]. For example, gene targets of five differentially expressed miRNAs (miR-365-3p, miR-30b-5p, miR-30c, let-7a-5p and miR-23a) were enriched for pathways in immune system (B-cell receptor signaling pathway, chemokine signaling, T-cell receptor signaling and Fc gamma R-mediated phagocytosis). These miRNAs are amongst several miRNAs with demonstrated roles in immunity. A recent study showed that the expression levels of three miRNAs including miR-23a were significantly higher in breast cancer with lymph node metastasis, compared with that from patients without lymph node metastasis or normal tissue and also that the expression of the miR-23a/24-2/27a cluster promoted mammary carcinoma cell migration, invasion, and hepatic metastasis, through targeting Sprouty2 (SPRY2) and consequent activation of p44/42 MAPK (mitogen-activated protein kinase) [53]. The MAPK pathway was highly activated in mammary cells and tissues challenged with S. uberis [23,25]. In addition, the role of Let-7 miRNA family in immunity is well established. The let-7 miRNA family was identified as the common denominator of Salmonellaregulated miRNAs in macrophages and epithelial cells, thus suggesting that repression of let-7 relieves cytokine IL-6 and IL-10 mRNAs from negative post-transcriptional control, thus establishing a paradigm of miRNA-mediated feed-forward activation of inflammatory factors when mammalian cells are targeted by bacterial pathogens [54]. Let-7b/d/e miRNAs, were recently shown to be differentially regulated in infected (S. uberis) bovine mammary cells [25].

Conclusion
Using deep sequencing, we characterized the miRNome of bovine mammary epithelia cells (MAC-T cells) challenged with heat inactivated E. coli or S. aureus bacteria and detected 231 known bovine miRNAs and 113 novel miRNAs, and also a pathogen dependent differential regulation of miRNAs. In particular, E. coli elicited an earlier differential regulation of miRNAs as opposed to a delayed regulation by S. aureus. Highly expressed and differentially regulated miRNAs have demonstrated roles in diverse biological processes as well as in immunity. GO and KEGG pathway analysis showed significant enrichments of predicted target genes of differentially regulated miRNAs in different functional groups (e.g. cellular process, developmental process, biological regulation, cell death, etc.) and KEGG pathways of the immune system, signal transduction, cellular process, nervous system, development and pathways in cancer and other human diseases. Importantly, our study provides further confirmation of the involvement of mammary epithelia cells in contributing to the immune response to infecting pathogens and a potential role for miRNAs as biomarkers in early diagnosis of mastitis and in development of control measures.

Cell culture and bacteria challenge
Cell culture and bacteria challenge were carried out as previously described with some modifications [4]. Briefly, MAC-T cells (a bovine mammary epithelial cell line) were seeded at a concentration of 1.5 × 10 5 cells in a 6 well cell culture plate (BD Biosciences, Mississauga, Ontario, Canada) and grown in a growth medium for 24 hours at 37°C in 5% CO 2 humidified incubator. The growth medium contained DMEM and RPMI 1640 (Wisent, St-Bruno, Quebec, Canada) at a concentration of 1:1, 10% fetal bovine serum (FBS) (Wisent), 10 μL/mL ITS (insulintransferrin-selenium solution) (Wisent), and 1% antibiotic antimycotic solution (100×) (Wisent). At 85% confluence, the growth medium was changed for infection medium (same as the growth medium without FBS) and cells were allowed to grow for another 24 hr period before infection with bacterial pathogens.
Escherichia coli (E. coli) strain P4 and Staphylococcus aureus (S. aureus) strain Smith CP were the infection agents. Bacteria were initially grown overnight on Luria Bertani (LB) agar (E. coli) or on tryptic soy agar (TSA) (S. aureus) aerobically in a humidified incubator at 37°C. A single colony of S. aureus was transferred to a 50 mL conical tube containing 20 ml of tryptic soy broth (TSB) and incubated at 37°C in an open air shaker at 225 RPM. Similarly, a single colony of E. coli was grown the same way in LB broth. The bacteria were grown until an OD 600nm of 0.6 was reached and then plated in triplicates on their respective media to confirm the number of bacteria per mL. Growth of bacteria and subsequent manipulations were carried out in the biosafety containment level II laboratory of the Dairy and Swine Research and Development Centre of Agriculture and Agri-Food Canada, Sherbrooke, following institutional safety procedures.
Prior to infection, cells in 6 wells were individually trypsinized and counted using the Countess® Automated Cell Counter (Life Technologies, Burlington, Ontario, Canada). The two bacterial strains were then diluted to achieve a concentration enabling a ratio of 10 bacteria to 1 MAC-T cell in the infection medium and then heat inactivated at 63°C for 30 minutes to prevent overgrowth during the period of cell challenge. After refreshing infection medium, each triplicate cell well representing different time points (6, 12, 24 and 48 hr) and each bacterial strain was challenged with the infection medium containing bacteria to achieve an infection rate of 1:10. Nonchallenged triplicates (control) for each time point (0, 6, 12, 24 and 48 hr) were also included. Uninfected cells were treated with the same volume of heat inactivated infection media without bacteria. After 0, 6, 12, 24 and 48 hr, the media were removed, cells washed once in Hank's balanced salt solution 1X (HBSS 1X; Wisent) without trypsinisation, harvested in 1 mL of lysis/binding buffer (mirVana miRNA Isolation Kit, Ambion Inc., Austin, TX, USA) and stored at −80°C until RNA extraction.

Total RNA isolation and library construction
Total RNA was extracted using the mirVana miRNA Isolation Kit (Ambion®, life technologies, USA) following manufacturer's protocol. The concentration and purity of the isolated RNA was checked by NanoDrop® ND-1000 Spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA). The quality of RNA was further assessed with the RNA 6000 Nano Labchip Kit using the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). All the 39 RNA samples had a RNA integrity number (RIN) of 9.7-10 and were stored at −80°C until further analysis.
Equal amounts (330 ng) of total RNA from triplicate samples at the same time point were pooled to construct small RNA library using the TruSeq Small RNA Sample Preparation kit (Illumina, San Diego, CA, USA) according to the manufacturer's instruction. PCR amplification was performed for 13 cycles. In total, 13 small RNA libraries were constructed and pooled in equal amounts for gel purification and sequencing. Sequencing was carried out on the HiScan SQ system (Illumina) using the TruSeq™ SBS kit v3 (50 cycles, Illumina). Real-time analysis and base calling was done using the HiSeq Control Software Version 1.4.8 (Illumina).

Small RNA sequence analysis
Low-quality reads were removed from the raw reads using CASAVA 1.8 based on chastity. After trimming the 3′ adaptor sequence by a perl script ("clip_adapters. pl") provided by miRDeep2 software (version 2.0.0.5), small RNA reads with 18-30 nts in length from all libraries were extracted and pooled to identify known and novel miRNAs using miRDeep2 with the default parameters [28]. All reads were aligned to the bovine genome (2009, UMD3.1) with Bowtie (allowing 1 mismatch) [55] and reads that mapped to bovine rRNAs, tRNAs, and snoRNAs in the Rfam RNA family database [56] were discarded. The remaining reads were scored by miRDeep2 as known or potential novel bovine miRNAs. Small RNAs that mapped only to genomic repeat loci were removed. The lowest score cut-off of 5 that yielded a signal-to-noise ratio of 16.2 [28] was used for novel miRNA predication and small RNAs with higher scores than threshold were regarded as potential novel miRNAs. Furthermore, small RNAs having lower scores than cut-off but more than 10 reads in total [30] were further analyzed as miRNA candidates. Small RNAs < 20 nucleotides in size were ignored because of higher false-positive discovery rates and the remaining novel miRNAs and miRNAs candidates were further annotated using Blastn [57] (word size = 6). Conserved novel miRNAs were classified as homologous known miRNAs or orthologous known bovine miRNAs and named according to the following guidelines: (1) it matches > 90% identity and > 90% of the length of a known animal miRNA (miRBase Release 19), or (2) it precursor aligns to any precursor of a known miRNA with > 70% identity and > 75% sequence coverage [58]. In addition, a huge bovine specific miRNA family (bta-miR-2284, miRBase Release 19) containing more than 80 members, shows very high similarity in their precursor sequences. Therefore, novel miRNAs were regarded as new members of bta-miR-2284 family if their precursors were clustered by CLUSTAL W [59] with bta-miR-2284 precursors. Non-conserved novel miRNAs were also named and inputted into miRBase, while the remaining miRNA candidates that didn't follow any of the categories above were not further analyzed.
The expression of miRNAs (known and novel miRNAs) for each library was calculated by miRDeep2 using the Quantifier module [28]. Reads that mapped equally well to the positions of two or more mature miRNAs were divided equally to the read counts of those mature miRNAs. To investigate the regulation of miRNAs in MAC-T cells in response to stimulation by pathogens, miRNA expression in challenged cells at different time points were compared to expression in unchallenged cells (controls) at the same time point using DEseq and raw reads as suggested [29]. DEseq allows a P-value to be determined in the absence of any available biological replicates, by treating the two conditions as replicates, under the assumption that only a small proportion of transcripts are differentially expressed [29]. Similarly, miRNA regulation in control cells over time was determined by comparing expression levels at different time points and 0 hr.

Target gene prediction and pathway analysis
The putative target genes of differentially expressed miRNAs were predicted using Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems Inc., Redwood City, CA, USA). IPA uses a miRNA target filter that depends on experimentally validated gene-miRNA interactions from TarBase database (http://diana.cslab.ece.ntua. gr/tarbase) and miRecords (http://mirecords.biolead.org/), and the presence of conserved 8mer and 7mer sites on 3 prime untranslated regions (3′UTR) of genes that match the seed region of each miRNA from TargetScan (http://www.targetscan.org). Additionally, manually curated findings within the IPA knowledgebase are also used. Predictions returned by IPA are either experimentally validated interactions, high or moderate interactions between the 8mer/7mer seed regions of miRNA that match sites on 3′UTR of genes. In this study, only genes showing high predictions and experimentally validated interactions were further analysed. Predicted gene targets were imputed into the Database for Annotation, Visualization and Integrated Discovery (DAVID version 6.7, http://david.abcc.ncifcrf.gov/) for gene ontology (GO) functional annotation and pathway analysis.

Quantitative PCR and data analysis
Five differentially expressed miRNAs (bta-miR-193a-3p, miR-423-5p, miR-21-3p, miR-365-3p and miR-486) were validated by quantitative RT-PCR using TaqMan® miRNA Assays following manufacturer's recommendations (Applied Biosystems, Foster City, CA, USA). Briefly, cDNA was reversely transcribed from 10 ng of total RNA using 5 × specific miRNA RT primer for each miRNA. PCR products were amplified from cDNA samples using 20 x TaqMan® miRNA Assays. Fluorescence signals were detected on an ABI StepOnePlus Real-time PCR System (Applied Biosystems). Triplicates for each reaction were performed. Bovine miR-192 was used as internal control since it was the best stable miRNA detected by sequencing (C.V =5.14%). The relative expression of target miRNA was calculated following the formula: ΔCt target miRNA = Ct target miRNA -Ct internal control . The control group at each time point was selected as a calibrator in each comparison group for the relative quantification. The relative expression of miRNA in a sample normalized to internal control and relative to the calibrator was calculated as follows: Relative expression target miRNA = 2 -ΔΔCt , where ΔΔCt = ΔCt target miRNA, sample -ΔCt target miRNA, calibrator . T-test (SAS version 9.2, 2008) was used to compare the difference between each comparison group. Significant differential expression was declared at P < 0.05.

Availability of supporting data
All the sequencing data sets supporting the results in this study have been deposited in the publicly available NCBI's Gene Expression Omnibus Database (http://www.ncbi. nlm.nih.gov/geo/). The data are accessible through GEO Series accession number GSE51979 (http://www.ncbi.nlm. nih.gov/geo/query/acc.cgi?acc=GSE51979).