Mastitis affects one-third of the dairy cattle in the United States and causes economic losses exceeding $2 billion a year in dairy farming (Sordillo & Streicher, Reference Sordillo and Streicher2002) mostly through its effects on milk production, an average milk yield loss per cow of 0–9% during the first lactation and 0–11% during the second lactation (Heikkilä et al. Reference Heikkilä, Nousiainen and Pyörälä2012). Reduced milk production is the largest indirect cost of mastitis. The immune response related to mastitis is a very complex biological process that involves immune cells, mammary epithelial cells and endothelial cells. Therefore, understanding the pathophysiological process and the host immune response to mastitis in mammary tissue on a molecular level is very important for the development of new strategies to control mastitis (Rinaldi et al. Reference Rinaldi, Li and Capuco2010). Streptococcus agalactiae (S. agalactiae) is one common and major cause of bovine mastitis (Trigo et al. Reference Trigo, Ferreira, Ribeiro, Dinis, Andrade, Melo-Cristino, Ramirez and Tavares2008; Zadoks et al. Reference Zadoks, Middleton, McDougall, Katholm and Schukken2011). S. agalactiae is a Gram-positive, obligate pathogen of the mammary gland that causes subclinical mastitis in 11–43% of infections (Keefe, Reference Keefe1997). The pathogen generally causes a low-grade persistent of infection and does not have a high self-cure rate (Farnsworth, Reference Farnsworth1987).
MicroRNAs (miRNAs) are small, non-coding RNA molecules approximately 21–23 nucleotides (nt) in length (Tang, Reference Tang2010). The mature miRNA is loaded onto the RNA-induced silencing complex, which combines with its target mRNA 3′ untranslated region (UTR) to inhibit translation of the mRNA or directly degrade the mRNA (Bartel, Reference Bartel2004; Filip, Reference Filip2007). miRNAs have been shown to control a wide range of biological processes, including cellular growth, differentiation, proliferation, apoptosis, metabolism and immune response (O'Connell et al. Reference O'Connell, Rao, Chaudhuri and Baltimore2010; Taft et al. Reference Taft, Pang, Mercer, Dinger and Mattick2010; Dorn, Reference Dorn2011). Recent studies have demonstrated that miRNAs help regulate the innate and adaptive immune responses. For example, miR-155, miR-146a, miR-21, and miR-9 regulate innate immune cell activation (Gantier, Reference Gantier2010). Of these, miR-145 silences the metastasis gene mucin 1 (MUC1) to inhibit cell invasion and metastasis in breast cancer cell lines (Sachdeva & Mo, Reference Sachdeva and Mo2010) whilst miR-146 mediates the congenital immune response to bacterial infections by targeting TNF receptor associated factor 6 (TRAF6) and interleukin receptor associated enzyme I (IRAK1) in alveolar A549 epithelial cells (Williams et al. Reference Williams, Perry, Moschos, Larner-Svensson and Lindsay2008). These examples suggest that certain miRNAs may be important in the modulation of the immune response.
The roles of miRNAs in the host immune response to infection are much less studied in livestock than in humans and mice (Sun et al. Reference Sun, Aswath, Schroeder, Lippolis, Reinhardt and Sonstegard2015). Several studies have reported a variety of differentially expressed miRNAs in cattle during bacterial infections, suggesting roles for these miRNAs in the host immune response. For example, five miRNAs (miR-184, miR-24-3p, miR-148, miR-486, and let-7a-5p) were up-regulated in bovine mammary epithelial cells challenged with Escherichia coli (E. coli), whereas four other miRNAs (miR-2339, miR-499, miR-23a and miR-99b) were up-regulated in the same cells when challenged with Staphylococcus aureus (Jin et al. Reference Jin, Ibeagha-Awemu, Liang, Beaudoin, Zhao and Guan2014). Additionally, five inflammation-related miRNAs (miR-9, miR-125b, miR-155, miR-146a and miR-223) were up-regulated in bovine monocytes stimulated with lipopolysaccharide and S. aureus enterotoxin B (Dilda et al. Reference Dilda, Gioia, Pisani, Restelli, Lecchi, Albonico, Bronzo, Mortarino and Ceciliani2012). These in vitro cell culture studies suggest that at least some miRNAs are involved in bovine mammary gland immunity in response to infection. However, the functions of miRNAs in the mammary gland immune response in vivo remain unknown. More studies on the bovine mammary gland may help establish the involvement of miRNAs in immune modulation in infectious diseases, such as mastitis. This knowledge will facilitate the development of biomarkers for mastitis diagnosis and therapeutic agents for treatment.
The aim of this study is to identify differentially expressed miRNAs in the bovine mammary gland following challenge with S. agalactiae and to discuss possible functions of the identified miRNAs. These findings will help provide an understanding of the roles of miRNAs in the process of S. agalactiae-type mastitis and the mechanisms of target gene expression and regulation.
Materials and methods
Ethics statement
This project was approved by the Animal Care and Use Committee of the College of Animal Science and Technology, Yangzhou University, China. All experiments were conducted according to the Regulations for the Administration of Affairs Concerning Experimental Animals published by the Ministry of Science and Technology, China, in 2004.
Animals
Primiparous Chinese Holstein cows (n = 3) were selected from a cattle farm on the outskirts of Yangzhou, Jiangsu, China. The cows were mid-lactation (175 ± 31 d) without a history of mastitis or treatments for mastitis. Milk samples were collected to determine somatic cell count (SCC, <100 000 cells/ml) and the presence of bacteria to confirm that the cows did not have subclinical mastitis or an intramammary infections. The cows were fed ad libitum and individually.
Bacterial strains
S. agalactiae (ATCC13813) was provided by the Veterinary College of Yangzhou University. S. agalactiae was suspended in 100 ml Todd-Hewitt Broth and cultured for 6 h at 37 °C in 5% CO2. Broth cultures were diluted in sterile mammal Ringer's solution (Oxoid, Basingstoke, England) to generate 1 × 106 CFU/ml S. agalactiae suspensions. The presence of S. agalactiae in milk samples was tested according to the Microbiological Examination of Food Hygiene – Examination of Streptococcus hemolyticus (GB/T 4789.11-2003).
Inducing mastitis in dairy cows
Two mammary quarters were assigned to the test and to the control group (online Supplementary Table S1). A 5-ml 1 × 106 CFU/ml S. agalactiae suspension or 5 ml of PBS was injected into the test or control mammary glands, respectively, through the milk duct via a tip needle.
Tissue sample collection
Mammary tissues were collected by surgical operation 24 h after challenge. Biopsy sites were selected to avoid obvious subcutaneous blood vessels and cisternae. The biopsy sites were clipped of hair and sterilised with an iodine surgical scrub and 70% ethanol three times. Next, 5 ml of lidocaine HCl was injected subcutaneously as local anaesthetic and a 6–7-cm incision made in the skin and fascia where the mammary gland capsules were visible. A sample of mammary gland tissue (≥1 g) was excised and wiped to remove any visible milk secretions and connective tissue. The sample was divided into two equal parts; the first part was placed into liquid nitrogen immediately, and the second part was placed into fixative solution to allow production of paraffin sections.
Histological examination
The paraffin-fixed blocks were serially sectioned into 6–7-μM coronal slices and stained with haematoxylin and eosin using traditional paraffin section methods. The haematoxylin and eosin-stained sections were analysed by light microscopy using a Nikon fluorescence microscope (Nikon, Tokyo, Japan).
Extraction of total RNA and Solexa sequencing
Total RNA was extracted from the cryo-frozen mammary gland tissues using a mirVana™ RNA Isolation Kit (Applied Biosystems AM1556, Carlsbad, CA, USA). Then, the RNA was purified using a QIAGEN RNeasy® Kit (QIAGEN, Mainz, Germany) and analysed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). The RIN value of the sample must be greater than 7 for high-throughput RNA sequencing. Subsequently, 5 µg of total RNA was sent to Shanghai Oebiotech Co. (Shanghai, China) for small RNA library construction and high-throughput sequencing.
Bioinformatics analysis
Sequencing reads from the small RNA libraries were pre-processed by removing adaptors, contaminants, inserts, polyA signals and low quality reads (with more than four bases whose quality was less than ten and more than six bases whose quality was less than thirteen). The length distribution of the small RNAs was summarised. Then, common and specific sequences among all of the samples were analysed. rRNAs, tRNAs, and other RNAs were identified by alignment using the GenBank database (http://www.ncbi.nlm.nih.gov) and Rfam database 10.1 (http://rfam.xfam.org/). The small RNAs were divided into several categories based on priority. Known miRNAs were identified by their alignment to the designated part of miRBase18.0 (http://www.mirbase.org/). The detailed criteria were: (1) Align the tags to the miRNA precursor in miRBase with no mismatches, and (2) Based on the first criterion, align the tags to the mature miRNA in miRBase with at least a 16-nt overlap to allow offsets. Most of the miRNA genes are located predominantly in the introns of protein-coding genes or intergenic regions (Starega-Roslan et al. Reference Starega-Roslan, Koscianska, Kozlowski and Krzyzosiak2011). The characteristic hairpin structure of a miRNA precursor can be used to predict novel miRNAs (Bizuayehu et al. Reference Bizuayehu, Fernandes, Johansen and Babiak2013). Novel miRNAs and their secondary structures were predicted using non-annotated small RNAs with the Mireap software (http://sourceforge.net/projects/mireap).
Differential expression analysis of miRNAs
Expression of miRNAs in the samples from the two groups (S. agalactiae and Control) was normalised to obtain the expression of transcripts per million (TPM). If the normalised expression value of a given miRNA was zero, its expression value was modified to 0·01. If the normalised expression (NE) of a given miRNA was less than 1 in both libraries, the miRNA was removed from the differential expression analysis. The fold change and P-value were calculated with the normalised expression using the following formulas (Audic & Claverie, Reference Audic and Claverie1997).
Normalised expression = (Actual miRNA read count/Total clean read count) × 106.
Fold change = Log2 (S. agalactiae-NE/Control-NE)
$$\eqalign{P(x\left \vert {y{\rm )}} \right. = & \left( {\displaystyle{{{\rm N}_{\rm 2}} \over {N_1}}} \right)\displaystyle{{(x + y)} \over {x!y!{(1 + (N_2/N_1))}^{(x + y + 1)}}} \cr C(y \le y_{\min} \left \vert {x)} \right. = & \sum\limits_{y = 0}^{y \le y\min} {\,p(y\left \vert {x)} \right.} \cr D(y \ge y_{\max} \left \vert {x)} \right. = & \sum\limits_{y \ge y_{\max}} ^\infty {\,p(} y\left \vert {x)} \right.} $$
N 1 and x represent the total count of clean reads and the normalised expression level of a given miRNA in the control sRNA library, respectively. N 2 and y represent the total count of clean reads and the normalised expression level of a given miRNA in the S. agalactiae RNA library, respectively.
Verification of miRNA expression by qPCR
Eight miRNAs were chosen from the differentially expressed miRNAs between the two groups to validate the accuracy of the Solexa sequencing. A miRcute miRNA First-Strand cDNA Synthesis Kit (TIANGEN, Beijing, China) was used to perform reverse transcription of the miRNA, and the SYBR green method was used to detect the expression level of the chosen miRNA. The miRNA primers consisted of a specific primer and a universal primer. The eight miRNA-specific primers (listed in Table 1) were synthesised by Sangon Biotech Co., Ltd. (Shanghai, China). The universal primers were supplied in the miRcute miRNA qPCR detection kit (TIANGEN). Bta-S18 was used as an internal reference. All operations were performed according to the kit's instructions. All reactions were run in triplicate. The fold changes of miRNA expression were calculated using the 2−ΔΔCT method (Livak & Schmittgen, Reference Livak and Schmittgen2001; Chao, Reference Chao2008; Chen et al. Reference Chen, Gelfond, Mcmanus and Shireman2009).
Primer sequences used for real-time PCR
Prediction and analysis of target genes
The RNAhybrid software was used to predict the target genes of the miRNAs that complied with the following seed region criteria: (1) no mismatches can occur among the 1–9 nt region on the 5′ end and (2) G-U is permitted, but the numbers cannot be more than 3. No other criteria were defined because miRNAs were combined with mRNAs (target genes) by permitting mismatches with the exception of the seed region. KEGG pathway analysis (Kanehisa et al. Reference Kanehisa, Araki, Goto, Hattori, Hirakawa, Itoh, Katayama, Kawashima, Okuda, Tokimatsu and Yamanishi2008) was conducted to identify significantly enriched metabolic pathways or signal transduction pathways in target gene candidates by comparison with the whole reference gene background. The formula below was used to calculate the gene numbers for each pathway and then applied the hypergeometric test to find significantly enriched pathways in the target gene candidates.
$$P = \sum\limits_{i = 0}^{m - 1} {\displaystyle{{\left( \matrix{M \hfill \cr i \hfill} \right)\left( \matrix{N - M \hfill \cr {\rm} n - i \hfill} \right)} \over {\left( \matrix{N \hfill \cr n \hfill} \right)}}} $$
Here, N is the number of all genes with KEGG annotation, n is the number of target gene candidates in N, M is the number of all genes annotated to a certain pathway, and m is the number of target gene candidates in M. A Bonferroni correction of the P-value was used to obtain a corrected P-value (q-value).
Results
Induction of S. agalactiae-type mastitis
Twenty-four hours following the S. agalactiae challenge, the test udders demonstrated clinical signs of mastitis, including redness, swelling, pyrexia, hardness and pain. The milk SCC from the test quarters increased rapidly (>2 000 000 cells/ml). In contrast, the SCC from the control udders was normal. The only bacteria detected in the test milk samples were S. agalactiae. Furthermore, the biopsy of the infected mammary gland tissues (Fig. 1a) stained with haematoxylin and eosin demonstrated pathological changes compared with the control (Fig. 1b). Obvious lesions appeared in the infected udders, the mesenchyme was swollen, parts of the lumen were atrophied, the mammary epithelial cells were loosely connected, and the intercellular gap was increased. Desquamated mammary epithelial cells, polymorphonuclear neutrophils, macrophages, lymphocytes and other inflammatory cells were concentrated in the lumen. Based on these signs, mastitis was induced successfully.
Haematoxylin and eosin staining of mammary gland tissue 400×. (a) Control mammary gland tissue: normal mammary with an integrated structure; (b) S. agalactiae-infected mammary gland tissue with loosely connected mammary epithelial cells and increased intercellular gaps.
Length distribution and miRNA composition of the small RNAs
A total of 20·85 and 18·58 M clean reads (accounting for 98·09 and 97·08% of the high-quality sequences, respectively) were found in the control and S. agalactiae groups, respectively (online Supplementary Table S2). Figure 2 shows the length distributions of the clean reads. The majority of the clean sequences had lengths of 21–23 nt in the two groups (accounting for 67·36 and 52·47% of the control and S. agalactiae group reads, respectively). The lengths of the reads peaked at 21 nucleotides (49·32% for the control and 35·96% for the S. agalactiae group), indicating that the clean reads included a large number of miRNA sequences.
Small RNA length distributions in the two groups. (a) control group: approximately 67·36% of the small RNAs had 21–23 nt lengths; (b) S. agalactiae group: approximately 52·47% of the small RNA had 21–23 nt lengths.
Summary of the small RNA composition
The common and specific tags for the species (unique) and the numbers (total) of the samples from the two groups are summarised. A total of 11·86% unique reads were common between the two groups of samples (Fig.3a), whereas the number of common reads between the samples of the two groups were concentrated (12 156 904; 90·62%) (Fig. 3b). These results demonstrated that the uniformity of the different samples was good at the whole sequencing level. In the test group, matching of the small RNA tags to the genome using the Short Oligonucleotide Analysis Package (SOAP) demonstrated that the unique and total small RNAs matched approximately 25·95 and 72·48% of the genome, respectively, whereas approximately 29·64 and 70·32% of the total small RNAs matched the genome in the control group (online Supplementary Table S3). To obtain unique small RNAs mapped to only one annotation, the rRNAetc priority rule was used (GenBank > Rfam > known miRNA) (online Supplementary Tables S4–S6). If a small RNA was not annotated, it was marked with ‘unann’. Most of the unique small RNAs were unann (accounting for approximately 89·39% of the unique RNAs in the control group and 90·62% in the S. agalactiae group) (online Supplementary Table S7). In the control group, miRNAs constituted the majority of the total small RNAs (total sRNAs; 70·31%), although only a small proportion (1·71%) of the small RNAs unique to the control group were miRNAs. In contrast, in the S. agalactiae challenged group, miRNAs accounted for 66·59% of the total small RNAs and 1·07% of the small RNAs unique to the group.
Venn diagram showing the unique sRNAs (a) and total sRNAs (b).
Known miRNA Analysis and novel miRNA prediction
A total of 517 miRNA precursors, 373 miRNAs, 59 miRNA-5ps and 59 miRNA-3ps were found in the test group, whereas 488 miRNA precursors, 358 miRNAs, 53 miRNA-5ps and 51 miRNA-3ps were found in the control group. In the test group, 10 376 unique sRNAs and 17 703 700 total sRNAs were matched to the miRNA precursors, whereas 8066 unique sRNAs and 14 380 650 total sRNAs were matched to the miRNA precursors in the control group (online Supplementary Table S6).
In total, 10 063 novel miRNAs were predicted in the control group, of which 232 were unique, and 14 127 novel miRNAs were predicted in the S. agalactiae group, of which 399 were unique (online Supplementary data file SD1).
Differential expression of known miRNAs
The expression of miRNAs in the two groups (S. agalactiae and Control) was demonstrated in a log2 ratio scatter plot (Fig. 4). A total of 35 differentially expressed known miRNAs that responded to S. agalactiae infection were identified (Table 2).
Differential expression of known miRNAs in the control and S. agalactiae groups (red points represent miRNAs with a ratio >2, blue points represent miRNAs with 1/2< ratio ≤2; and green points represent miRNAs with a ratio ≤1/2; Ratio = normalised expression of the S. agalactiae group/normalised expression of the control).
Differentially expressed known miRNAs in the Control and S. agalactiae groups
† std: Normalised expression level of the miRNA in a sample.
‡ Fold change (log2*/*): Fold change in the miRNA expression in the pair of samples.
Verification of miRNA expression by qPCR
The qPCR results for the eight chosen miRNAs (including four miRNAs with up-regulated and four miRNAs with down-regulated expression) are shown in Fig. 5. The expression levels of bta-miR-136, −378, −3431, and −3660 were decreased, and the expression levels of bta-miR-16a, −451, −383, and −451 were increased. All of the miRNAs had the same expression trends as those observed in the Solexa sequencing results. These results verified that the Solexa sequencing was reliable.
Comparison of the qPCR results with the Solexa sequencing results.
Prediction and analysis of target genes
miRNAs regulate mRNAs (target genes) incompletely in animals. The RNAhybrid software was used to predict the miRNA target genes and obtained 18 801 target genes and 82 832 sites for the 35 differentially expressed miRNAs (online Supplementary data file SD2).
KEGG analysis showed that the predicted target genes were related to immune system pathways, including the RIG-I-like receptor signalling pathway and cytosolic DNA-sensing pathway, signal transduction pathways, including the Notch signalling pathway, and cellular processes, included endocytosis, tight junctions and adherens junctions (Table 3). Therefore, we hypothesised that the differentially expressed miRNAs regulated the host immune response to S. agalactiae infection by targeting these predicted genes.
Pathway annotation
Discussion
miRNAs play important roles in sustaining cellular homeostasis. Abnormal miRNA expression is associated with specific diseases (Vimalraj & Selvamurugan, Reference Vimalraj and Selvamurugan2013). Studying the adaptive changes in miRNAs in mammary tissues in response to mastitis can elucidate the mechanisms by which miRNAs modulate the immune response and provide appropriate strategies to control mastitis. In this study, we identified 35 differentially expressed miRNAs, including 10 up-regulated and 25 down-regulated miRNAs, in bovine mammary gland tissues with S. agalactiae-type mastitis. Of these, miR-223 exhibited the highest level of up-regulation. In previous investigations of bovine mastitis, miR-223 was also up-regulated by 2–5 fold (Naeem et al. Reference Naeem, Zhong, Moisá, Drackley, Moyes and Loor2012; Li et al. Reference Li, Zhang, Liao, Chen, Wang, Zhu, Geng, Ji, Mao, Gong and Yang2015). miR-223 plays pivotal roles in inflammation; for example, this miRNA negatively regulates neutrophil granulocyte proliferation and differentiation and reduces extended factor E2F1 expression to limit the cell cycle process (Fazi et al. Reference Fazi, Rosa, Fatica, Gelmetti, Marchis, Nervi and Bozzoni2005; Johnnidis et al. Reference Johnnidis, Harris, Wheeler, Stehling-Sun, Lam, Kirak, Brummelkamp, Fleming and Camargo2008). miR-223 also inhibits several signalling pathways by targeting IGR1R (Naeem et al. Reference Naeem, Zhong, Moisá, Drackley, Moyes and Loor2012). The validated targets of miR-223, including GZMB, IKKa, RC3H1 and STAT3, have effects on inflammation and infection (Haneklaus et al. Reference Haneklaus, Gerlic, O'Neill and Masters2013). This study showed that some of the predicted target genes of miR–223 participated in the immune pathway during mammary tissue infection with S. agalactiae. Thus, this miRNA could be useful as a candidate biomarker of mastitis in future research.
Similar to miR-223, miR-16 was also related to the immune response. Previous studies indicated that one important function of miR-16a was maintenance and regulation of the levels of inflammatory mediators in cells (Jing et al. Reference Jing, Huang, Guth, Zarubin, Motoyama, Chen, Padova, Lin, Gram and Han2005). miR-16a up-regulated interleukins IL-6, IL-8, and IL-10 following challenge with S. uberis in cows (Naeem et al. Reference Naeem, Zhong, Moisá, Drackley, Moyes and Loor2012). miR-16 also modulated macrophage polarisation, the inflammasome and NF-kB signalling (Chen et al. Reference Chen, Liu, Li, Wang, Liu, He, Yang, Niu, Wang and Guo2014).
miR-136 has been reported to function as a tumour suppressor in recent studies. For example, miR-136 suppressed tumour invasion and metastasis by targeting RASAL2 in triple-negative breast cancer (Yan et al. Reference Yan, Li, Tong, Han, Zhao, He and Jin2016). Additionally, miR-136 regulated two identified anti-apoptotic genes (AEG-1 and Bcl-2) in human glioma (Yang et al. Reference Yang, Wu, Guan, Cai, Fang, Li and Li2012). However, the function of miR-136 in bovine mastitis is unclear. In this study, CD93 was one of the predicted target genes of miR-136. CD93, which is a membrane-associated glycoprotein on the surface of cells, mediates phagocytosis, inflammation and cell adhesion (Yanaba et al. Reference Yanaba, Asano, Noda, Akamata, Aozasa, Taniguchi, Takahashi, Ichimura, Toyama, Sumida, Kuwano, Tada, Sugaya, Kadono and Sato2012). These functions provide a research direction to probe the relationship between miR-136 and CD93.
miR-3660 expression was down-regulated in this experiment, but the function of miR-3660 was unclear. Interestingly, VWF was the predicted target gene of miR-3660. VWF is a compound that non-covalently combines with glycoprotein and coagulation factor VIII in the plasma. Recent studies have suggested that VWF exerts an additional antitumor effect by negatively modulating angiogenesis and apoptosis (Franchini et al. Reference Franchini, Frattini, Crestani, Bonfanti and Lippi2013). KEGG analysis suggested that VWF took part in the complement and coagulation cascade pathways, ECM-receptor interactions and focal adhesion. These pathways were related to the immune system, signalling molecule interactions and the cellular community.
Other miRNAs (i.e., miR-335 and miR-378) have been shown to be related to fat metabolism and adipocyte differentiation (Gerin et al. Reference Gerin, Bommer, McCoin, Sousa, Krishnan and MacDougald2010; Iliopoulos et al. Reference Iliopoulos, Drosatos, Hiyama, Goldberg and Zannis2010). Due to the potential involvement of host immune responses to mastitis with other metabolic changes (i.e., fatty acids and proteins), the predicted targets suggested the involvement of complex immune responses and the synergy of biological processes.
Conclusion
Two small RNA libraries corresponding to the control and S. agalactiae-type mastitis mammary gland tissues were constructed, and 35 differential miRNAs were identified. The bioinformatics analysis revealed that these miRNAs had a total of 18 801 predicted target genes. The target genes of these miRNAs were enriched in pathways related to the immune response. We should continue to study the functions of these miRNAs, especially miR-223, miR-16, miR-3660 and miR-136, and verify the modulated relationships between these miRNAs and their predicted target genes based on the bioinformatics analysis results provided in this study.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022029917000437.
This research was supported by the National High Technology Research and Development Programme of China (2013AA102505-5), National Natural Science Foundation of China (Grant No. 31472067, 31272407, 31372286), Project Funded by the Priority Academic Programme Development of Jiangsu Higher Education Institutions, and Jiangsu province agriculture ‘San Xin’ project (SXGC [2015] 328).
Authors’ contributions
Zhangping Yang and Yongjiang Mao conceived the study. Rui Li, Xiangxiang Liao and Dan Chen constructed the small RNA libraries. Junhua Pu conducted the qRT-PCR, performed the bioinformatics analysis and wrote the manuscript. Chenglong Zhang, Wenqiang Wang, and Xiaohan Geng collected samples and performed the SCCs and bacterial testing. Yihui Zhu performed the mammary gland tissue sample haematoxylin and eosin staining. Yunchen Gong reviewed the manuscript. Dejun Ji and Yongjiang Mao read the manuscript. All authors approved the final manuscript.
