Expression profiling of peripheral blood miRNA using RNAseq technology in dairy cows with Escherichia coli-induced mastitis

E. coli is the main causative agent of mastitis in dairy cows, but the mechanism of molecular regulation underlying the occurrence and development of mastitis has not yet been fully elucidated. In this study, an E. coli-induced mastitis model was created and RNASeq technology was used to measure the miRNA expression profiles at different times post-infection (0, 1, 3, 5, 7 dpi), as well as to screen for differentially expressed miRNA. The results show detection of 2416 miRNAs, including 628 known miRNAs and 1788 newly discovered miRNAs. A total of 200 differentially expressed miRNAs were found at different time points. Bioinformatics analysis showed that these differentially expressed miRNAs may regulate the occurrence and development of mastitis in dairy cows through seven signal transduction pathways, namely cytokine-cytokine receptor interaction, MAPK signaling pathway, chemokine signaling pathway, leukocyte transendothelial migration, T cell receptor signaling pathway, Toll-like receptor signaling pathway, and cell adhesion molecules. In addition, bta-miR-200a, bta-miR-205, bta-miR-122, bta-miR-182 and the newly discovered conservative_15_7229 might be involved in immune process in late stage of E. coli-induced mastitis. The results of this study lay the foundation for molecular network analysis of mastitis and molecular breeding of dairy cows.

miRNA expression overall distribution. Based on our deep sequencing data, we calculated the TPM (transcripts per million) values of the miRNAs from the five sample groups. The overall expression pattern of miRNAs is shown in Fig. 3. As shown in Fig. 3, there are differences in the expression patterns of miRNAs between the five samples, suggesting that these differentially expressed miRNAs may play different roles in the development of mastitis. In addition, the abundance of expressed miRNAs (TPM >1000) (Supplementary Table S3) shows that the top 10 highly expressed miRNAs in all samples are: bta-miR-486, bta-miR-451, bta-miR-92a, bta-let-7f, bta-miR-25, bta-let-7i, bta-let-7g, bta-miR-26a, bta-miR-21-5p and bta-miR-191. In addition, out of all the new miRNAs, unconservative_3_22065 had the highest level of expression in all five samples with a TPM of over 1200.    Table S5 for sequence information). Pairwise comparisons revealed that, compared to the control group (BE0), there were fewer miRNAs differentially expressed in BE1 and BE5, while more miRNAs were significantly differentially expressed in BE3 and BE7 (Fig. 4). It is noteworthy that the expression of bta-miR-200a and bta-miR-205 in BE1, BE5 and BE7 groups were significantly increased (bta-miR-200a log 2 FC between 4.36 and 7.13, bta-miR-205 log 2 FC between 1.93 and 4.70). At the same time, the expression of bta-miR-122 was significantly downregulated in BE1 and BE3 (log 2 FC was -1.82), while it was significantly upregulated in BE5 and BE7 (log 2 FC 0.91 and 1.82, respectively). This indicates that bta-miR-122 plays complex roles in mastitis. The expression level of bta-miR-200a and bta-miR-205 in blood significantly increased on the first day after infection, and decreased to pre-infection level on the 3rd day after infection, while on the 5th and 7th day it was significantly increased (Fig. 5). In addition, our study found that compared with the BE0 group, the expression of the novel miRNA conservative_15_7229 in BE5 and BE7 groups was significantly decreased. The above results indicate that these 4 miRNAs (bta-miR-122, bta-miR-200a, Our results showed that, on days 1, 3, 5, and 7 following induction, there were many differentially expressed miRNAs in dairy cow blood following E. coli-induced mastitis. The differences in expression were very pronounced (Supplementary Table S4). Notably, bta-miR-182 expression significantly increased in BE1_vs_BE7, BE3_vs_BE5, BE3_vs_BE7 and BE5_vs_BE7 comparisons (log 2 FC of 2.65, 2.30, 3.95 and 1.64, respectively) (Supplementary Table S4). This indicates that bta-miR-182 is significantly upregulated at 5 or 7 dpi and may also participate in the immune process in the late stage of mastitis.
Differentially expressed miRNA qPCR validation. Ten differentially expressed miRNAs (five of which were upregulated and five of which were downregulated) were randomly selected to perform qPCR experiments and the qPCR results were compared with the high-throughput sequencing results as shown in Fig. 6. The results show that qPCR results are in complete agreement with the results of high-throughput sequencing, which indicates that our findings of differential expression in miRNAs through high-throughput sequencing are accurate and reliable.

DIE-miRNA target gene prediction and functional annotation and KEGG signaling pathway analysis.
To explore the function of DIE-miRNA in the development of mastitis, we predicted the target genes of DIE-miRNAs, and our results show that the 200 DIE-miRNAs may have 18,369 target genes. The results of KEGG pathway enrichment analysis showed that these target genes may be involved in the following 6 aspects: disease, environmental information processing, cellular processes, metabolism, organismal systems, and genetic information processing, and encompass 50 significant signal pathways, of which 7 signal pathways (cytokine-cytokine receptor interaction, chemokine signaling pathway, leukocyte transendothelial migration, T cell receptor signaling pathway, Toll-like receptor signaling pathway, and cell adhesion molecules) are associated with immunity ( Fig. 7). Therefore, we conclude that these DIE-miRNAs may participate in the disease development process by regulating their target genes during E. coli-induced mastitis.  . qPCR validation of differentially expressed miRNAs. qPCR expression was normalized with cel-miR-39. There were 10 differentially expressed miRNAs for qPCR validation, with bta-miR-200a, bta-miR-205 and bta-miR-342 from the BE0 vs BE7 comparison; conservative_15_7229 from BE0 vs BE5; bta-miR-326 and conservative_X_30752 from BE3 vs BE5; bta-miR-145 and bta-miR-331-3p from BE3 vs BE7; bta-miR-122 from BE1 vs BE5; bta-miR-214 from BE1 vs BE7 comparison.

Discussion
Dairy cow mastitis is an inflammatory disease caused by pathogenic microorganisms. Studies have shown that the expression level of immunity-related genes rapidly changes after pathogens infect tissues or cells, and the process of disease occurrence and development can be finely regulated [18][19][20] . Therefore, screening differentially expressed genes may play a very important role in the research of disease molecular mechanism or rapid diagnostic technology. In recent years, with the development of sequencing technology, high-throughput sequencing techniques are often used to mine new miRNAs related to immunity and disease and to screen for differentially expressed miR-NAs. In the case of dairy cow mastitis, screening of differentially expressed miRNAs based on high-throughput sequencing has focused on both mammary epithelial cells and mammary tissue. For example, Jin et al. 21 found 113 novel miRNAs in E. coli and S. aureus-infected bMECs and found 10 DIE-miRNAs. Lawless et al. 13 screened 21 DIE-miRNAs from S. uberis infected bMECs. In addition, other researchers have also used RNASeq technology to measure the miRNA expression profiles of E. coli and S. aureus-infected mammary tissues and the results show that there are multiple differentially expressed miRNAs in infected mammary tissues 10,22,23 .
Although changes in the expression of blood miRNAs have been shown to be putative biomarkers for early diagnosis of disease 7,9 , the differential expression of miRNAs can also be used as biomarkers of animal performance 24 . However, little has been done so far on miRNAs in the blood of dairy cows affected by mastitis. Lawless et al. 16 detected three DIE-miRNAs (bta-miR-146b, -451, and -411a) in blood monocytes in dairy cows infected with S. uberis in vivo. Li et al. 14 and Chen et al. 15 used qRT-PCR and Illumina's Genome Analyzer IIe to detect DIE-miRNAs in the peripheral blood of mastitis-affected dairy cows in the clinic. Their results showed that, compared with healthy cows, there were 123 mastitis-related miRNAs that experience significant upregulation (such as bta-miR-16a, bta-miR-125b, bta-miR-15a, bta-miR-29b, bta-miR-301a, bta-miR-21-3p, bta-miR-181a and bta-miR-181b, etc.) and 50 miRNAs that undergo significant downregulation (such as bta-miR-223, bta-let-7f, bta-miR-375, bta-miR-148a, etc.). However, it is surprising that the results of the above two scholars vary greatly, and some findings are even in opposition to one another. For example, their results show that bta-miR-146a and bta-miR-146b are significantly downregulated in mastitis cow blood in one reference 15 , but are significantly upregulated in another reference 14 . Thus it is clear that the specific dairy cow population, feeding environment and the type of infection and time of infection (or the degree of mastitis) have a significant effect on blood miRNA expression profiles. Therefore, we felt that an in-depth study of E. coli-induced mastitis miRNA expression profile was essential to clarify the molecular mechanism of E. coli-induced mastitis. Therefore, we used a small dose of E. coli for a long-term infection of mammary tissue, and used RNASeq technology to detect miRNA expression in dairy cow blood at different times after infection. Our results show that 2416 miRNAs were detected, including 628 known miRNAs and 1788 newly predicted miRNA ( Fig. 2 and Supplementary Table S2). The result of comparison of blood miRNA expression at different times revealed that 200 DIE-miRNAs (including 76 known miRNAs and 124 novel miRNAs) were differentially expressed, of which bta-miR-200a, bta-miR-205, bta-miR-122, bta-miR-182 expression were significantly upregulated in late mastitis, while the newly found con-servative_15_7229 was significantly downregulated. As for the function of the differentially expressed conservative_15_7229, thus far, there have been no reports on this miRNA. A closely related gene found in recent years is miR-200a. The expression of miR-200a is significantly increased in many types of disease, consistent with the results of this study. Zhao et al. 25 showed that the expression of miR-200a-3p is increased in hepatocytes of alcoholic liver disease, and may target ZEB2 to induce apoptosis. Wang et al. 26 have shown that the expression of miR-200a in CD4 + T cells correlates with the expression of Th17/Treg cells and related cytokines in psoriasis vulgaris. It is noteworthy that, Zarate-Neira et al. 27 found that mouse miR-200a-3p expression was significantly increased in a systemic lupus erythematosus (SLE) model and can inhibit the TLR4/MyD88 signalling pathway and promote interferon production, thereby regulating the immune response in mice and proving its usefulness as a biomarker of SLE. In addition, the expression of miR-200a is significantly increased in fish lymphocytes 28 and micro-vascular endothelial cells 29 induced by LPS (cell surface component of E. coli), targeting the expression TLR1 28 and ZEB1 30 genes to regulate the development and progression of disease 30,31 . However, the current relationship between miR-200a and mastitis in dairy cows is unknown. Therefore, in order to explore the function of bovine miR-200a, we predicted its target gene, and the results indicated that ZEB1 gene may be the target gene of bovine miR-200a. In addition, we found that the 3′-UTR of the ZEB1 gene of other species, such as bovine and human, are highly homologous and that the sequences of bovine and human miR-200a are substantially identical (Supplementary Fig. S1). Therefore, we conclude that bovine miR-200a may also target the bovine ZEB1 gene to regulate the development of E. coli-induced dairy mastitis disease. The specifics remain to be further studied.
As for the expression of miR-122, the results of this study show that the expression level of bta-miR-122 in dairy cow blood is significantly increased in the late stage of E. coli infection (5 dpi and 7 dpi). Consistent with the trend of our results, miR-122 also significantly increases in the bloodstream late in the murine liver injury model and has been shown to be a biomarker of liver injury 32 . In addition, it has been shown that miR-122 is involved in the regulation of Toll-like receptor signaling pathway after Vibrio anguillarum infection of fish (miiuy croaker) by targeting TLR14 33 . Therefore, we hypothesize that bta-miR-122 may also be involved in inflammatory response in the late stage of E. coli-induced mastitis.
There are only a few studies on the roles of miR-205 and miR-182 in immune regulation. Studies have shown that the expression of miR-205 is significantly upregulated in patients with allergic rhinitis 34 , and that miR-205 can also regulate the expression of erbB2/erbB3 in breast cancer cells to promote apoptosis 35 . IL-2 cytokine stimulation of Th cells can result in a significant increase in the expression of miR-182, reaching its maximum three days after stimulation. Furthermore, in vitro and in vivo studies using the mouse arthritis model showed that IL-2 induced mouse miR-182 expression shows a dose-dependent increase, while inhibiting miR-182 may be beneficial in the treatment of arthritis 36 . Therefore, we hypothesized that the expression of IL-2 and other inflammatory cytokines in E. coli-induced dairy cow mastitis was low in early stages, but that during later stages IL-2 and other inflammatory cytokines continued to accumulate and induced a significant increase of miR-182 expression in the late stage of E. coli-induced mastitis.
In addition, in order to explore the function of DIE-miRNAs, we predicted the target genes of 200 DIE-miRNAs in E. coli-induced mastitis and analyzed their KEGG signaling pathway. The results indicated that these DIE-miRNAs may be involved in seven animal immunity-related signalling pathways, of which the Toll-like receptor signaling pathway 37 and chemokine signaling pathway 38 are closely related to innate immunity and inflammatory response, while the T cell receptor signaling pathway 39 is related to adaptive immunity. Therefore, we concluded that these DIE-miRNAs may be involved in the innate or adaptive immunity in dairy cow mastitis and may thus regulate the development of mastitis.
In summary, this study used RNASeq technology to detect the miRNA expression profile in the peripheral blood of dairy cows infected with E. coli-induced mastitis at five different post-induction times. The results show that a total of 2416 miRNAs were detected, of which 628 miRNAs were known, and 1788 miRNAs were unknown. A total of 200 DIE-miRNAs were found to be involved in cytokine-cytokine receptor interaction, MAPK signaling pathway, chemokine signaling pathway, leukocyte transendothelial migration, T cell receptor signaling pathway, Toll-like receptor signaling pathway, and cell adhesion molecules in a total of seven signal transduction pathways that regulate the development of mastitis. In addition, bta-miR-200a, bta-miR-205, bta-miR-122, bta-miR-182 and the newly discovered conservative_15_7229 might be involved in immunity in late stage of dairy cow mastitis caused by E. coli. The results of this study have laid the foundation for molecular network analysis and molecular breeding of mastitis, and also provide new molecular materials for the biotherapy of mastitis.

Materials and Methods
Ethics Statement. All animal procedures were conducted with attention to animal welfare and humanitarian procedures according to guidelines laid down by the Institutional Animal Care and Use Committee of Northwest A&F University, and the protocols were approved by the Experimental Animal Manage Committee of Northwest A&F University.
Experimental animals and mastitis model induction. Experimental animals and methods used for the induction and verification of mastitis models have been previously described 22 . In brief, healthy 2-year-old half-sibling dairy cows were divided into control and experimental groups. Animals in the experimental group were injected with 5 mL of 10 5 CFU/mL E. coli (ATCC 25922) via the nipple. The control group was injected with an equivalent amount of PBS. We adopted a painless surgical method to collect mammary tissue and used hematoxylin-eosin staining and pathology analysis to ensure the reliability of the model.
Small RNA library construction and high-throughput sequencing. Based on protocols from previous publications 22,40,41 , 1.5 μg of qualified RNA was taken and a small RNA library was constructed using the method described in the NEBNext ® Ultra ™ Small RNA Sample Library Prep Kit for Illumina ® (NEB, USA). After the library was constructed, the library concentration was tested using Qubit 2.0 and was adjusted to 1 ng/μl. Insert Size was measured with an Agilent 2100 bioanalyzer, and qPCR was used to accurately determine the effective concentration of the library to ensure library quality. After qualification of the library, high-throughput sequencing was performed with the Illumina HiSeq2500 of BioMerker Biotech Co., Ltd. (Beijing).
High-throughput sequencing quality control. Sequencing reads were single-end (SE) 50 nt. In order to ensure accurate and credible analysis results, the raw data from each sample after sequencing was screened as follows: (1) low quality sequences were removed; (2) sequences with unknown base N (unidentifiable base) content greater than or equal to 10% of reads were removed; (3) reads without 3′ linker sequences and insertions were removed; (4) reads contaminated with 5′ linkers were removed; (5) cleavage of 3′ adapter sequences and T/C/G reads were removed; (7) sequences shorter than 18 nt or longer than 30 nt were removed. This process yielded high-quality, clean, reliable data. sRNA classification annotation. Clean reads were aligned using Silva, GtRNAdb, Rfam, and Repbase databases using Bowtie software 17 to filter out non-coding RNAs such as rRNA, tRNA, snRNA and snoRNA and repeated sequences to obtain the unannotated reads.
UMD3.1 (http://www.ensembl.org/Bos_taurus/Info/Index) was used as the reference genome for sequence alignment and subsequent analysis. Using miRDeep2 software 42 , unannotated reads were aligned with the reference genome for positional information on the reference genome, these are subsequently referred to as mapped reads. miRNA identification. To determine miRNA sequences in the unannotated reads, known miRNA identification and novel miRNA prediction were performed using the miRDeep2 software package 42 , i.e. reads and genomic sequences were aligned to obtain possible miRNA precursor sequences using RNAfold and randfold to calculate the energy of the precursor structure and to predict the secondary structure. Combined with information on subsets of sequences (e.g. miRNA production characteristics, mature, star and loop sequences), RNAfold and randfold were scored using the Bayesian model to finally identify miRNAs. MiRNAs with scores greater than 80% were considered to be reliable. miRNA expression and differential expression analysis. Expression of miRNA in each sample was statistically analyzed and was normalized using the TPM algorithm 43 . TPM normalization formula is as follows: = × TPM (specific miRNA reads 1, 000, 000)/total miRNA mapped reads The miRNA differential expression analysis between two samples was performed using IDEG6 software 44 . To remove statistical false positives, Benjamini Hochberg's correction method was used to correct the p-value to the q-value during analysis 45 . Differentially expressed miRNAs were screened using |log 2 (Fold Change) | ≥ 1, FDR ≤ 0.01, and q-value < 0.005 as standard. miRNA expression qPCR verification. Ten differentially expressed miRNAs (five upregulated and five downregulated miRNAs) were randomly selected from the sequencing data, and qPCR was used to verify the expression results from RNASeq. The specific method used was as follows: using miRcute miRNA isolation kit (Tiangen Biotech Co., Ltd. Beijing, China), miRNA extraction was performed according to the manufacturer's instructions, and the purity and concentration of the miRNA was measured by a trace UV spectrophotometer and diluted to 50 ng/μl. In addition, at the time of miRNA extraction, exogenous cel-miR-39 (5′-UCACCGGGUGUAAAUCAGCUUG), which has no bovine homologue, was added to each tube for subsequent normalization of qPCR quantification results 46,47 . The reaction system consisted of 4 μl of miRNA template, 10 μl of 2× miRNA RT Reaction Buffer, 2 μl of miRNA RT Enzyme Mix, and 4 μl RNase-free H 2 O. The reaction conditions were 42 °C for 60 min and 95 °C for 3 min to produce the first strand cDNA, which was then diluted to a concentration of 80 ng/μl. The qPCR reaction was performed using the miRcute Plus miRNA qPCR Detection Kit (SYBR Green) (Tiangen Biotech Co., Ltd. Beijing, China). The reaction system consisted of: 1 μL cDNA template, 10 μL 2× miRcute miRNA Premix (with SYBR & ROX), 0.4 μL Reverse Primer (10 μM), 0.4 μL Forward Primer, and 8.2 μL RNase-Free ddH 2 O. The reaction conditions were as follows: 95 °C for 15 min, 5 cycles of 94 °C for 20 s, 64 °C for 30 s, 72 °C for 34 s, followed by 45 cycles of 94 °C for 20 s, then 60 °C for 34 s. Experiments using qPCR followed the MIQE guidelines 48 : NTC and NRC groups were set up for each qPCR experiment to ensure that the entire experiment was free of DNA or RNA contamination. The qPCR primer sequences are shown in Supplementary  Table S1. The qPCR amplification efficiency of each primer ranged from 99.8% to 103.9% (R 2 > 0.99). To ensure the specificity of primer amplification, primer melting curve detection was performed at the end of qPCR. All quantitative experiments were performed in three technical duplicates.
The miRNA expression level was calculated using the 2 -ΔΔCt method 49 , and the exogenous cel-miR-39 was used to normalize the quantitative results. The expression level was expressed as mean ± SD, and Student's t-test was used to test the significance of the difference. A P value of less than 0.05 indicated a significant difference.