MicroRNA expression profiling of porcine mammary epithelial cells after challenge with Escherichia coli in vitro

Coliform mastitis is a symptom of postpartum dysgalactia syndrome (PDS), a multifactorial infectious disease of sows. Our previous study showed gene expression profile change after bacterial challenge of porcine mammary epithelial cells (PMECs). These mRNA expression changes may be regulated through microRNAs (miRNAs) which play critical roles in biological processes. Therefore, miRNA expression profile was investigated in PMECs. PMECs were isolated from three lactating sows and challenged with heat-inactivated potential mastitis-causing pathogen Escherichia coli (E. coli) for 3 h and 24 h, in vitro. At 3 h post-challenge with E. coli, target gene prediction identified a critical role of miRNAs in regulation of host immune responses and homeostasis of PMECs mediated by affecting pathways including cytokine binding (miR-202, miR-3277, miR-4903); IL-10/PPAR signaling (miR-3277, miR-4317, miR-548); and NF-ĸB/TNFR2 signaling (miR-202, miR-2262, miR-885-3p). Target genes of miRNAs in PMECs at 24 h were significantly enriched in pathways associated with interferon signaling (miR-210, miR-23a, miR-1736) and protein ubiquitination (miR-125, miR-128, miR-1280). This study provides first large-scale miRNA expression profiles and their predicted target genes in PMECs after contact with a potential mastitis-causing E. coli strain. Both, highly conserved miRNAs known from other species as well as novel miRNAs were identified in PMECs, representing candidate predictive biomarkers for PDS. Time-dependent pathogen clearance suggests an important role of PMECs in inflammatory response of the first cellular barrier of the porcine mammary gland.


Background
Coliform mastitis (CM) represents one of the most important cardinal symptoms of postpartum dysgalactia syndrome (PDS) in sows [1]. This multifactorial infectious disease has high economic relevance for the pig industry, because both postpartum sows and their piglets can be critically affected by mastitis and lactation failure [2]. Husbandry and management conditions, pathogen contamination, and genetic predisposition are suspected etiological factors that contribute to PDS [3]. Gram-negative pathogens such as Escherichia coli (E. coli) are most commonly isolated from milk of PDS-positive sows, but are found in non-affected sows as well [4,5]. Therefore, it is unclear which determinants result in development of subclinical or clinical signs of infection within 12-48 h postpartum, or otherwise preserve clinical health.
Lipopolysaccharide (LPS) endotoxin is a pathogenic factor of E. coli and can induce acute and severe inflammation in sows and other animal species [6]. Our previous study showed that E. coli induce a fast and strong inflammatory response in porcine mammary epithelial cells (PMECs) [7]. In these cells, bacterial challenge strongly upregulates mRNA expression of genes encoding cytokines, chemokines, and cell adhesion molecules. Furthermore, these changes suggest rapid activation of other immune-competent cells mediating pathogen clearance and host homeostasis through interleukin 1 beta (IL1B) and tumor necrosis factor (TNF) expression [7]. Therefore, cellular as well as molecular factors may be highly relevant to animals' susceptibility to intramammary infection and mastitis.
MiRNAs are small non-coding RNA molecules (∼22 nucleotides long) and are highly conserved across higher eukaryotic species and can modulate host immune responses [8] as well as are involved in a variety of signaling pathways and biological processes, such as development, cell proliferation, cell differentiation, cell death, metabolism, and disease [9].
In contrast, the function of miRNAs in innate and adaptive immunity of porcine mammary glands has not yet been described. Therefore, our study focused on miRNA expression changes in PMECs induced by a potential mastitis-causing E. coli strain isolated from milk of PDS-positive sows. We characterized miRNA expression differences in PMECs at 3 h and 24 h after challenge with heat-inactivated E. coli to investigate and to compare early and late phase of immune response since specific kinetics and extents of global changes in the transcriptome of bovine mammary epithelial cells were found in a study with comparable experimental design [17]. The response of PMECs was determined by comparing miRNA expression profiles between pathogenchallenged and unchallenged control groups. Together with our previous mRNA expression profiles, miRNA expression data were then integrated based on pairwise correlations and computational target predictions to identify potentially affected signaling pathways. Finally, we performed selective analysis of the most and strongest affected biological processes, molecular functions, and Kyoto Encyclopedia of Genes and Genomes (KEGG)/ canonical pathways to investigate the role of miRNAs in pathogen clearance of PMECs.
Statistical analysis identified 143 significantly differentially expressed (DE) probesets (p < 0.05; fold change >1.5 or < −1.5) corresponding to 102 mature miRNA sequences (41 up-and 61 downregulated) at 3 hpc with E. coli when compared with unchallenged controls ( Fig. 1 with E. coli ( Fig. 1; Additional file 1). Some of these miR-NAs are not yet annotated in Sus scrofa. As concern to this point, we compared these miRNA sequences with Sus scrofa miRNA in miRBase (http://www.mirbase.org) and alignment with pig genome (Sus scrofa 10.2 download from NCBI: http://www.ncbi.nlm.nih.gov/ on 1.9.2015) and added in Additional file 1 in the column of pig chromosome and position. On the basis of the top 100 significantly DE miRNAs after pathogen challenge, a heat map was generated to identify distinguishable miRNA expression profiles between short-term (3 hpc) and long-term (24 hpc) challenge with E. coli compared to unchallenged controls (0 hpc) (Fig. 2). Hierarchical clustering analysis determined similar patterns in miRNA expression profiles of both pathogen-challenged groups compared to unchallenged control group (Fig. 2).
Prediction of miRNA targets identifies biological processes and molecular functions that differ at 3 hpc and 24 hpc with E. coli Pairwise correlation coefficient analysis was performed to evaluate associations of expression levels between 442 significantly DE miRNA probesets and 1345 significantly DE mRNA probesets (previously identified by Jaeger et al. [7]). In total, 36,610 miRNA-mRNA pairs were significantly negatively correlated (p < 0.05) and in silico predicted as functionally linked. TargetScan and RNAhybrid software were used to further filter predicted miRNA-mRNA pairs by setting the threshold for microRNA:target duplex energy to −25 kcal/mol. Finally, we identified 135 miRNA-mRNA pairs corresponding to 34 mature miRNA sequences and 71 genes at 3 hpc with E. coli (Table 1; Additional file 2). At 24 hpc with E. coli, 531 miRNA-mRNA pairs corresponding to 67 mature miRNA sequences and 214 genes were filtered ( Table 2; Additional file 3).

Validation of selected miRNA expression by real time quantitative PCR
Seven miRNAs (miR-210, miR-423-5p, miR-17-3p, miR-193a-5p, miR-320, miR-339-5p, and miR-362) that were DE in miRNA microarray analyses were selected for validation by real time quantitative PCR (RT-qPCR). MiRNA expression analyzed by RT-qPCR was in good accordance with microarray data, as shown by similar fold changes (FC) (Additional file 6). This suggested that our microarray data were reliable, although we used pools of samples for microarrays compared to individual samples for qPCR.

Discussion
The role of miRNAs in innate and adaptive immunity of the mammary gland is increasingly well described in different species, e. g. human, mice, and cattle. In pigs, the miRNA role regulating inflammatory defense mechanisms during microbial infection of the mammary gland has not yet been described. Therefore, we used the previously described PMEC model [7] to analyze specific pathogen defense mechanisms under standardized experimental conditions in vitro. To the best of our knowledge, this is the first study describing specific miRNAs and their predicted target genes in PMECs after contact with a potential mastitis-causing E. coli strain, indicating DE miRNAs that were identified in PMECs at 3 hpc and 24 hpc with E. coli were analyzed together with our previous mRNA expression profile and then integrated based on pairwise correlations and computational target predictions.
Pathogen challenge interfered with many host cellular processes, including modulating miRNA production, which influences host cell physiology and defense. As shown in this study, many miRNAs were significantly upregulated in PMECs at 3 hpc and 24 hpc with E. coli. Expression data showed that long incubation of PMECs with heat-inactivated E. coli resulted in more DE miR-NAs than short incubation. Early response of PMECs at 3 hpc with E. coli was followed by a late, more intensive host response at 24 hpc, as indicated by an approximately two-fold increase of significantly DE miRNAs. This is in accordance with our mRNA microarray analysis, which also identified a more intense host response after long-term challenge with E. coli [7].
Approximately 46% and 17% of identified mRNAs are predicted to be miRNA-regulated in PMECs at 3 hpc and 24 hpc with E. coli, respectively. Some of the identified DE miRNAs have not yet been described in the current literature. But based on our target prediction analysis and existing literature, most of the identified DE miRNAs likely play important roles in regulation of innate immunity and cellular homeostasis of PMECs. Among these, miR-101, miR-24-star, miR-30b-star, miR-210, miR-34c, and miR-17-3p were significantly DE in PMECs after E. coli challenge and are also affected during host responses of human, mice, or cattle to LPS Fig. 5 Top KEGG pathways and canonical pathways regulated by predicted target genes of significantly DE miRNAs at 3 hpc and 24 hpc with E. coli. a Predicted target genes of predominantly downregulated (black bars), but also upregulated (grey bars), miRNAs from PMECs at 3 hpc and 24 hpc with E. coli were significantly enriched in KEGG pathways such as cytokine-cytokine receptor interactions, apoptosis, cancer pathways, NODlike receptor signaling pathways, and Jak-STAT signaling pathways. b Predicted target genes of downregulated miRNAs were significantly enriched in canonical pathways such as IL10 signaling, NF-kB signaling, PPAR signaling, TNFR2 signaling, and the role of osteoblasts, osteoclasts, and chondrocytes in rheumatoid arthritis in PMECs at 3 hpc with E. coli. Predicted target genes of up-and downregulated miRNAs were significantly enriched in canonical pathways such as interferon signaling, protein ubiquitination, glioma invasiveness signaling, myc-mediated apoptosis signaling, and Gq signaling in PMECs at 24 hpc with E. coli stimulus or E. coli infection [12,[20][21][22][23][24]. In addition, some DE miRNAs in our analysis play pivotal roles in regulation of host defenses against other bacterial infections (miR-23a-star, miR-125-star, and miR-128) [12,25,26]; viral infections (miR-29a, miR-505, and miR-423-5p) [27][28][29]; parasitization (miR-8-star) [30]; and inflammatory bowel disease (miR-10a and miR-5128) [31,32]. Furthermore, miR-101, miR-29b, miR-10a, miR-17-3p, and miR-92b were DE expressed in PMECs after pathogen challenge and are important regulators in bovine (endo-) metritis [33]. In addition, miR-484, miR-660, miR-34c, miR-339, miR-193a-5p, and miR-362 have regulatory functions in breast cancer progression [34][35][36][37][38][39] and were also deregulated in PMECs after pathogen challenge. Likewise, some miRNAs that were DE in PMECs at 3 hpc and 24 hpc with E. coli are potential regulators of genes that were enriched in KEGG pathways, such as cancer-related or apoptosis pathways. Therefore, a deregulation of growth regulation or an alteration of cellular physiological processes such as cell proliferation, cell death and repair of genes in PMECs by E. coli infection is assumed, which might play a pivotal role in the genesis of PDS.
Staedel and Darfeuille [25] reported that recognition of intruding pathogens by PAMPssuch as membrane-associated TLRs and cytoplasmic NOD-like receptorsis essential for subsequent activation of innate and adaptive immunity, mostly stimulated by NF-ĸB signaling. Accordingly, we identified miR-202, miR-2262, miR-4031-5p, miR-4317, miR-548u, and miR-885-3p as potential regulators of target genes that are enriched in NOD-like receptor signaling pathways in PMECs, referring to recognition of E. coli. In contrast, the TLR signaling pathway was not identified as one of the top enriched canonical pathways in PMECs after pathogen challenge as discuss in our previously study [7]. It is possible that the preparation of heat kill particles somehow removed TLR agonists from the particles from pMECs but not from bMECs [40]. This is still lack of knowledge in this field of experiment. However, miR-101, miR-125-star, miR-8-star, and miR-29a were DE in PMECs after pathogen challenge and are involved in TLR signaling [20,25,30].
MiRNAs can interact with several signaling molecules, such as cytokines, chemokines, and transcription factors, and therefore influence various signaling pathways and cellular processes [25]. In our study, miR-202, miR-3277, miR-4903, miR-648, miR-8-5p, miR-885-3p, and miR-890 were identified to play critical roles in regulation of genes that control cytokine receptor activity and cytokine binding. Additionally, some miRNAs that were DE in PMECs after pathogen challenge are known to regulate genes that are involved in NF-ĸB signaling (miR-210 and miR-23a-star) [41,42] or production of cytokines IL-6, TNF-alpha, IL-1 beta, and IL-10 (miR-30b-star, miR-210, miR-17-3p, miR-23a-star, miR-24-star, miR-29b, miR-125-star, miR-1934-star, and miR-320) [21,22,24,[42][43][44][45][46][47]. Furthermore, our analysis identified predicted target genes of DE miRNAs in PMECs at 3 hpc with E. coli that are enriched in top canonical pathways such as IL-10 signaling, PPAR signaling, NF-ĸB signaling, TNFR2 signaling, or arthritis-associated roles. These signaling pathways are important for fast pathogen recognition and effective nonspecific removal of infectious agent as well as moderate the inflammatory response to limit tissue damage of the host. At 24 hpc of PMECs with E. coli, predicted target genes of DE miRNAs were enriched in top canonical pathways such as interferon signaling and protein ubiquitination, which were also significantly enriched in our previous mRNA analysis of pathogen-challenged PMECs [7]. In this late immune response phase these signaling pathways can increase host-defense by activation of immune cells and regulation of a variety of fundamental cellular processes.
Altogether, the number of target genes predicted for DE miRNAs and involved in regulation of biological processes or molecular functions in PMECs was two-to tenfold higher at 24 hpc than at 3 hpc with E. coli.

Conclusions
Combined with current literature, our data show that many functions of several miRNAs seem to be highly conserved across species. However, while many miRNAs can have similar functions, specific miRNAs also can have different functions in PMECs. We suppose the identified miRNAs have critical roles in pathogen recognition and initiation of local and systemic immune responses in PMECs. Altogether, our results strongly suggest that DE miRNAs in PMECs contribute to pathogenesis of porcine mastitis induced by E. coli. The detection of several known and novel miRNAs in context with predicted target genes in PMECs is another indicator for the immune competence of these cells and provides a strong impulse to further examine porcine mammary gland immune defense mechanisms.

Cell culture and pathogen challenge
PMECs are originated from our previous study [7]. Therefore, cell cultures were established from mammary glands of three lactating sows of commercial herds. Animal care and tissue collection was performed in compliance with the German Law of Animal Protection. The experimental protocol was approved by the Animal Care Committee of the Leibniz Institute for Farm Animal Biology, Dummerstorf, Germany. Animals were exsanguinated following electronarcosis and dissected. Cryopreserved mammary epithelial cells were thawed and resuspended in complete medium consisting of Dulbecco's Modified Eagle Medium/Nutrient Mixture F-12 (DMEM/ F12, PAN Biotech, Aidenbach Germany), 10% fetal bovine serum (FBS, PAA, Cölbe, Germany), 1% Antibiotic/Antimycotic Solution (APS, 10,000 U/ml penicillin, 10 mg/ml streptomycin sulphate, 25 μg/ml amphotericin B, PAA), 10 μg/ml insulin (Sigma-Aldrich, Taufkirchen, Germany) and 1 μg/ml hydrocortisone (Sigma-Aldrich).
E. coli strain (gMEc240, sequence type 101, phylogroup B1, C+) used for this study is an isolate from milk of PDS-positive sows and was provided by the Institute of Microbiology and Epizootics, Department of Veterinary at the Freie Universität Berlin, Berlin, Germany. E. coli was grown in brain-heart-infusionbroth (BHB, Oxoid, Wesel, Germany) at 37°C to the logarithmic phase of culture growth (Optical Density at 600 nm [OD 600 ] 0.5~5 × 10 7 /ml). Bacteria were heatinactivated at 80°C for 1 h, spun down at 3000 rpm for 15 min, washed twice with Dulbecco's Modified Eagle Medium/Nutrient Mixture F-12 (DMEM/F12) medium and resuspended herein at a density of 1 × 10 8 /ml.
The experimental design was identical with the procedure as described before [7]. Briefly, approximately 4.4 × 10 5 cells from each 3 individual (3 biological replicates) were seeded and cultured in collagen-coated 6-well plates (1:10 collagen R in destilled water, Menal, Emmendingen, Germany) in complete medium without APS (three technical replicates per individual and treatment condition). After 24 h, medium was changed. Within the next 48 h after seeding, cells reached about 90% confluency. Then PMECs were challenged with 10 7 / ml heat-inactivated E. coli for 3 h and for 24 h. Unchallenged cells were used as a negative control. After incubation periods, medium was discarded, and cells were washed three times with phosphate buffered saline (PBS, PAA) to remove the bacteria. Afterwards, cells were collected for isolation of total RNA and small RNA.

RNA isolation
Total RNA was isolated from each of the 27 samples using the TRI® reagent (Sigma-Aldrich) and extracted with phenol-chloroform according to the manufacturer's instructions. Isolated RNA was purified using RNeasy Mini Kit (Qiagen, Hilden, Germany), and contaminating DNA was removed by DNase I digestion (Qiagen). Isolation and enrichment of small RNA was performed by using a miRNeasy kit and an RNeasy MinElute Cleanup kit (both from Qiagen) according to the manufacturer's protocols. RNA integrity and quantity were assessed with an Agilent 2100 Bioanalyzer using an RNA 6000 Nano kit for total RNA and a separate kit for small RNA (both from Agilent Technologies Inc., Santa Clara, CA, USA). The mean RIN and standard deviation for the samples used in the analysis was 9.72 ± 0.24. microRNA microarray analysis and statistics microRNA expression profiling was performed using hybridization to GeneChip miRNA 3.0 Arrays (Affymetrix Inc., Santa Clara, CA, USA) according to the manufacturer's recommendations. The array has a very high sensitivity and specificity and offers 100% coverage of miRBase 17 (www.mirbase.org) by a one-color approach. In addition, the array contains 16,772 entries representing hairpin precursor, total probe set 19,724 for detection of 1733 mature miRNA, 2216 human snoRNA and scaRNA products in 153 species, and provides a greater than 3-log dynamic range with higher than 95% reproducibility and 85% transcript detection at 1.0 amol for a total RNA input of 130-500 ng.
In total 27 small RNA samples were isolated. Equivalent amount of small RNA from three technical replicates per individual was pooled. In total 9 pooled samples of small RNA were used for further experiment. An amount of 200 ng of each small RNA sample pool was labeled with the FlashTag Biotin HSR RNA Labeling kit (Genisphere, Hatfield, PA). Labeled RNA was hybridized for 16 h to the miRNA arrays according to the manufacturer's recommendations. A total of nine microarrays were obtained. Afterwards, arrays were washed and stained in a Fluidics Station 450 (Affymetrix), and scanned on the G3000 GeneArray Scanner (Affymetrix). Pre-processing of raw expression data was done using the Expression Console™ software (Affymetrix). Robust multi-array average (RMA) background correction, log-2 transformation, quantile normalization and statistical analysis were performed using SAS 9.4 containing JMP Genomics 7 (SAS Institute, Cary, NC, USA). The ANOVA procedure in JMP genomics 7 was performed to determine relative changes in miRNA levels, including effects mediated by 3 experimental group (before (0 h) and after challenged (3 h and 24 h) with each 3 replicates of each experimental group. Alterations in transcript abundances were considered to be statistically significant at p < 0.05. Subsequently, data were filtered by fold change >1.5 or < −1.5. The miRNA microarray data were deposited in the Gene Expression Omnibus (GEO) public repository (GEO accession number GSE88870: GSM2350436 -GSM2350444).

Prediction and functional annotation of miRNA target genes
We used our previous microarray-based mRNA expression data to integrate with the DE miRNAs and scan for potential target genes [7]. TargetScan was used to detect predicted target genes based on seed complementarity on both 3′-, 5′-UTR and coding sequences of the porcine mRNA sequences (Sus scrofa 10.2 download from NCBI: http://www.ncbi.nlm.nih.gov/ on 1.9.2015) and miRNA sequences [86]. RNAhybrid software (http:// bibiserv.techfak.uni-bielefeld.de/rnahybrid) was used for direct prediction of multiple, energetically most favorable potential binding sites of our differentially expressed miR-NAs [87,88]. microRNAs were tested due to the following default parameters: (i) number of hits per target, 1; (ii) energy cutoff, −25 kcal/mol; (iii) maximal internal or bulge loop size per side, 4. Pearson correlation of miRNA and mRNA expression levels was calculated based on the predicted miRNA-mRNA relationships. We selected only negative correlation pairs of miRNA-mRNA with p < 0.05. The most affected biological processes and molecular functions as well as enriched KEGG pathways of predicted target genes were determined using the Database for Annotation, Visualization and Integrated Discovery 6.7 (DAVID) [89]. In addition, the top canonical pathways, miRNA targets are involved in, were identified using the Ingenuity Pathway Analysis software (p ≤ 0.05, Fisher's exact test; IPA, Ingenuity Systems, Redwood City, CA, USA).

Validation of microRNA microarray results by Fluidigm real-time quantitative PCR
Transcript levels of seven miRNAs were measured by RT-qPCR using the Fluidigm BioMark HD System (Fluidigm Corporation, San Francisco, CA, USA). First, 200 ng of miRNA was converted to cDNA using the Megaplex RT Primers, Human Pool Kit v2.1/v3.0 (Thermo Fisher Scientific, Waltham, MA, USA). The pre-amplification sample mixtures were prepared containing 1.25 μl of diluted cDNA, 2.5 μl of 2× Taqman pre-amplification master mix (Applied Biosystems, Waltham, MA, USA), and 0.5 μl of the primer pool (200 nmol each primer/μl; Promega, Mannheim, Germany; Additional file 6). The pre-amplification PCR program consisted of a 10 min 95°C denaturation step and 14 cycles of 15 s at 95°C and 4 min at 60°C. The preamplified reactions were diluted 10 times in Tris-EDTA (TE) buffer. Five microliters from sample mix containing pre-amplified cDNA, pooled primers and EvaGreen master mix (Fluidigm BioMark) were combined as outlined in the manufacturer's protocol. The 48.48 dynamic array chip (Fluidigm Corporation) was first primed with control line fluid and then loaded with the sample and assay mixtures via the appropriate inlets using an IFC controller. The array chip was placed in the Fluidigm BioMark Instrument for PCR at 95°C for 10 min, followed by 30 cycles at 95°C for 15 s and 60°C for 1 min. The samples were investigated in duplicate, and data were analyzed with RT-PCR analysis software in the BioMark HD instrument. The mean Ct values of the different miRNAs were normalized to U6 snRNA and miR-24-3p (internal controls), and then subjected to analysis using the 2 −ΔΔCt method [90]. Both internal controls were selected, because their transcript levels were most stable over time and treatment in PMECs. In addition, U6 snRNA is known as one of the most commonly used internal control for the normalization of miRNA qPCR data.