- 1College of Veterinary Medicine, Shandong Agricultural University, Tai‘an, China
- 2Research Center for Animal Disease Control Engineering, Shandong Agricultural University, Tai‘an, China
- 3China Animal Health and Epidemiology Center, Qingdao, China
- 4College of Veterinary Medicine, China Agricultural University, Beijing, China
Staphylococcus aureus (S. aureus), a common mastitis pathogen widespread in the natural environment of dairy farms, is capable of invading mammary epithelial cells making treatment difficult. However, the mechanism of the response of bovine mammary epithelial cell to S. aureus invasion remains elusive. In this study, transcriptomic analysis and bioinformatics tools were applied to explore the differentially expressed RNAs in bovine mammary epithelial cells (bMECs) between the control and S. aureus-treated group. A total of 259 differentially expressed mRNAs (DEmRNAs), 27 differentially expressed microRNAs (DEmiRNAs), and 21 differentially expressed long non-coding RNAs (DElncRNAs) were found. These RNAs mainly enrich the inflammatory response, immune response, endocytosis, and cytokine-cytokine receptor interaction. qRT-PCR was used to analyze the quality of the RNA-seq results. In particular, to the defense mechanism of bovine mammary epithelial cells against intracellular S. aureus, the PPAR signaling pathway and the genes (ACOX2, CROT, and NUDT12) were found to be up-regulated to promote the production of peroxisomes and ROS, DRAM1 expression was also up-regulated to facilitate the activation of autophagy, indicating that the above mechanisms were involved in the elimination of intracellular S. aureus in bovine mammary epithelial cells.
Introduction
Staphylococcus aureus is a frequently isolated pathogen responsible for bovine mastitis (1). This bacterial disease is economically significant in dairy cows, causing considerable economic losses and a series of food safety concerns (2). Despite the tremendous progress in understanding the pathogenesis of bovine mastitis, the disease remains an important issue in the dairy industry worldwide (3, 4).
Related studies have shown that S. aureus invades and survives in bovine mammary epithelial cells, in bMECs protected from host defenses, and in the bactericidal effect of some antimicrobials, thus making the treatment of mastitis caused by S. aureus difficult and prone to recurrence (5–7). Moreover, although S. aureus is resistant to serval antibiotics (8–10). The administration of antibiotics is currently the most common method for treating bovine mastitis (11). Hence, understanding the potential molecular regulatory mechanisms of S. aureus invading bovine mammary epithelial cells is particularly important.
MicroRNAs (miRNAs) are endogenous small non-coding RNAs (22–25 nucleotides) universally expressed in higher eukaryotic cells that play a crucial role in post-transcriptional gene regulation (12). Previous research has shown microRNAs as a vital part of the mammalian host response to bacterial infection, involved in the host immune response (13–16). Long non-coding RNAs (lncRNAs) are a class of non-coding RNAs over 200 nucleotides in length that are essential regulators of the immune response. Several researchers have reported the specific involvement of lncRNAs in the response of host cells to bacterial infection (17, 18). Their role in regulating gene-encoding products involved in the immune response, including direct interactions with chromatin, RNA, and proteins has been one of the most recent discoveries (19, 20). In general, the characterization of RNA regulatory networks represents a new area in the field of host-pathogen interactions.
High-throughput transcriptome sequencing has effectively been used to explore candidate mRNAs, lncRNAs, and miRNAs with a differential expression that may be involved in specific biological processes (21–23). Thus, it laid the foundation for subsequent integration of whole transcriptome analysis. To date, several studies have focused on the bovine mammary tissue or epithelial cells transcriptional response to S. aureus showing significant changes in gene expression following S. aureus infection (24–28). However, complete transcriptome analysis of the response of bovine mammary epithelial cells to intracellular S. aureus has not been reported. This study, for the first time, details the interpretations of whole transcriptome bMECs infected with intracellular S. aureus.
The purpose of this study was to explore the transcriptional regulation of bovine mammary epithelial cells after S. aureus invasion and identify related candidate mRNAs, lncRNAs, and miRNAs. Finally, the possible functions of the identified RNAs were also discussed. These findings provide a base for the study on the pathogenic mechanism of intracellular S. aureus and offer several potential targets for the treatment of S. aureus-mastitis.
Materials and Methods
Bacterial Strains and Growth Conditions
Staphylococcus aureus strain ATCC25923 was inoculated on a Luria-Bertani (LB) Agar and incubated at 37°C for 24 h. A single colony was randomly selected and cultured in LB broth with agitation at 37°C for 12 h and the growth was monitored by measuring the OD600nm.
Cell Culture
MAC-T cells are sensitive to the regulation of the extracellular matrix and lactogen. MAC-T cells after differentiation can synthesize and secrete α- and β-casein. It represents an in vitro model of bovine lactation (29). Bovine mammary epithelial cell line MAC-T cells were cultured in a cell culture plate with a growth medium consisting of Dulbecco's modified eagle culture medium (DMEM) with 10%(v/v) FBS and maintained in 5% CO2 at 37°C. Cells that cultured to monolayer confluence were used for further experiments.
Intracellular Infection Model
The model of intracellular infection was established in our previous study (30). In brief, MAC-T cells were seeded on 6-well cell plates, and when cells were cultured to monolayers, S. aureus (OD600nm = 0.8–1.2) was inoculated with (treatment group) or without (control group) an MOI of 8:1. After 2 h of incubation, the cells were washed three times in PBS. The cells from the treatment and control groups were further placed into fresh medium supplemented with lysostaphin (10 μg/mL) and gentamicin (50 μg/mL) to kill extracellular bacteria. After 12 min, extracellular fluid was collected for plate culture experiments to verify that extracellular bacteria were eliminated. The cells were again washed three times with PBS to remove extracellular bacteria and dead cells (30), and incubated for 2 h in 10% FBS-DMEM.
RNA Extraction and cDNA Library Construction
The total RNA from the control and S. aureus-treated group (three samples per group) were extracted using RNAiso Plus (Takara, Dalian, China) according to the manufacturer's instructions. Qualitative and quantitative total RNA was quantified using Nano Drop and the Agilent 2100 Bioanalyzer (Thermo Fisher Scientific, MA, USA). These RNAs were divided into two parts for library construction of RNA or small RNA, respectively.
Total RNA was treated with DNase I which degraded double-stranded and single-stranded DNA and ribosomal RNA (rRNA) was removed using the Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, CA, USA). Fragmentation of the purified RNA was carried out at a specific temperature and ionic environment. The first strand cDNA was produced by PCR in the first strand reaction system, along with the second-strand cDNA. A-Tailing Mix and RNA Index Adapters were added and incubated to carry out end repair, and purification was performed using Ampure XP Beads. Distribution of fragment sizes in the sequencing database was examined using the Agilent 2100 Bioanalyzer, and the libraries were quantified using qRT-PCR (TaqMan Probe). Finally, the qualified libraries were pair-end sequenced on the Illumina Hiseq 4000 platform (BGI, Shenzhen, China).
Quality Control of Raw Data
Before alignment, FastQC was used to check the quality of the raw reads generated by the Illumina Hiseq 4000 platform that were filtered by removing the dirty raw reads. Reads containing adapter, an unknown base >10%, and low-quality reads (The base number of threshold mass ≤10 accounts for more than 50% of the total reading) were removed to obtain clean reads of mRNA. PolyA, adapter, low-quality, and length <18 nt tags were removed to get clean reads of miRNAs.
Reads Alignment and Differential Expression Analysis of RNA-Seq
The transcriptome data were submitted to the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/sra/), with the BioProject ID, PRJNA591729. In this study, gene model annotations and reference genomes (ARS-UCD1.2/bosTau9) were accessed from UCSC (31), lncRNA model annotations were accessed from NONCODE (32), and miRNA model annotations were obtained from miRbase (33).
The clean reads of each sample were mapped using HISAT2 for mRNA and lncRNA to the index from the reference genome. StringTie was used to assemble transcripts for each sample to obtain all assemblies. After initial assembly, assembled transcripts were merged by the StringTie merge module, creating a uniform set of transcripts for all samples. StringTie was run again by the merge function to re-calculate the abundances of the merged transcripts after merging all assemblies and generated read coverage tables. StringTie also provided a Python script to calculate the count matrix of mRNA or lncRNA directly from the file created in the previous step (34).
Conversely, Bowtie was applied for miRNA to map the clean reads of each sample to the index from the reference genome (35). miRDeep2 calculated the count matrix of each miRNA for each sample (36).
The DESeq2 package in R identified the count matrix between the control and S. aureus-treated group, and the results according to the P-values (37). The pheatmap in R were used for clustering. lncRNAs, miRNAs, and mRNAs with a P < 0.05 and | log2(fold-change) | > 1 were set as the threshold for being differentially expressed.
Protein–Protein Interaction (PPI) Analysis
Understanding all functional interactions between expressed proteins is essential for a system-wide understanding of cellular function. The STRING database aims to consolidate known and predicted protein–protein association data. Network analysis of DEmRNAs protein-protein interactions were obtained using the STRING database (38), and Cytoscape was used to draw the PPI network.
Prediction of the Target Gene of lncRNAs and miRNAs
The target genes of the lncRNAs were identified using a large-scale RNA-RNA interaction prediction tool, RIsearch2 (39). The target genes of the miRNAs were predicted using TargetScan and miRWalk (40, 41).
GO Annotation and KEGG Pathway Analysis of DEmRNAs and Target Genes of DEmiRNAs and DElncRNAs
Gene Ontology (GO) enrichment analysis was performed using the DAVID 6.8 functional annotation tool. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were analyzed by the clusterProfiler R package (42, 43).
Competitive Endogenous RNAs (ceRNAs) Analysis of DEmRNAs, DEmiRNAs, and DElncRNAs
DEmRNAs, DEmiRNAs, and DElncRNAs crosstalk networks were constructed based on the hypothesis of competitive endogenous RNA (ceRNA) (44). ceRNA act as a sponge for microRNA (miRNA) through its binding site, and changes in ceRNA abundance in individual genes can regulate miRNA activity (45). Thus, understanding this new RNA crosstalk could provide insight into gene regulatory networks. Correlation analysis was performed in R (P < 0.05 & cor < − 0.4) and the networks were constructed by cytoscape software (46).
qRT-PCR Identification of Differentially Expressed mRNAs
Randomly selected five (PIDD1, TNFRSF18, DEPDC5, TMEM102, TRIM32) and three (CD36, DRAM1, CXCR1) genes of biological interest from DEmRNAs were tested by qRT-PCR to verify the reproducibility and repeatability of the differentially expressed genes in RNA-Seq data. Total RNA was extracted using RNAiso Plus (Takara, Dalian, China) and reverse-transcribed using the PrimeScript® RT regent kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China) according to the manufacturer's instructions. Quantitative PCR analysis was performed using the LightCycler® 96 (Roche Diagnostics GmbH, Germany) with TB Green™ Premix Ex Taq™ II (Tli RNaseH Plus) Kit (Takara, Dalian, China) following the manufacturer's instructions. The reaction mixtures in a 96-well plate were run at 95°C for 120s, followed by 40 cycles at 95°C for 15 s and 65°C for 30 s. qRT-PCR was performed in triplicate for each cDNA sample to ensure repeatability. The glyceraldehyde-3-phosphate dehydrogenase (GAPDH), actin beta (ACTB), and peptidylprolyl isomerase A (PPIA) genes were set as endogenous reference genes. All primers used in this study are listed in Table 1.
Statistical Analysis
The relative fold change of target gene expression was calculated by the 2−ΔΔCt method (47), and a Students t-test examined significant differences between conditions. ΔCt value was determined by subtracting the target Ct of each sample from the GAPDH, ACTB, and PPIA Ct values.
Results
Transcriptome Assembly Profiles Evaluation
The quality of six mRNA and lncRNA transcriptome sequence expression profiles of mammary epithelial cells are summarized in Table 2. Clean data were obtained after filtration from downstream analysis.
Analysis of Differentially Expressed mRNAs
Cluster pattern analysis of DEmRNAs between the control (n = 3) and S. aureus-treated groups (n = 3) is illustrated in Figure 1A. A total of 259 DEmRNAs (Supplementary File 1) were filtered by the thresholds of P < 0.05 and | log2(fold-change) | > 1, out of which 124 were up-regulated and 135 were down-regulated (Figure 1B).
Figure 1. Screening and enrichment analysis of differently expressed mRNAs (DEmRNAs) in S. aureus infected mammary epithelial cells compared with non-infected mammary epithelial cells. (A) Cluster analysis of DEmRNAs in mammary epithelial cells between the control group (C1, C2, and C3) and treatment group (T1, T2, and T3). Red indicates highly expressed genes, and blue indicates low expressed genes. Each column represents a sample, and each row represents a gene. On the left is the tree diagram of mRNA clustering, and on the right is the name of each mRNA. The closer the two mRNA branches are, the closer their expression level is. The upper part is the tree diagram of sample clustering, and the bottom is the name of each sample. The closer the two-sample branches are to each other, the closer the expression pattern of all genes in the two samples are and the trend of the more recent gene expression. (B) Volcano plot of global DEmRNAs in mammary epithelial cells between the control group and treatment group. Red dots (Up) represent significantly up-regulated genes [P < 0.05, | log2(fold-change) | >1]; blue dots (Down) represent significantly down-regulated genes [P < 0.05, log2(fold-change) < −1]; black dots represent insignificantly differential expressed genes. (C) KEGG pathway classified annotation of DEmRNAs in mammary epithelial cells. The pathway is exhibited in the left axis, and the size of the circle represents the number of genes listed in the right axis. (D) The relationship between these pathways is illustrated.
Gene Ontology (GO) enrichment analysis showed that 259 DEmRNAs were distributed into 83 GO terms (Supplementary File 2) and divided into three categories based on molecular function, biological process, and cellular components (Figure 2). The top 10 enrichment GO terms, such as the metal ion binding, integral component of the membrane, extracellular exosome, intracellular, and extracellular space were also significantly enriched.
Figure 2. Annotation of DEmRNAs using Gene Ontology (GO) in mammary epithelial cells. The number of mRNAs for each GO annotation is exhibited in the right axis, and the proportion of genes for each GO annotation is listed on the left axis.
The enrichment analysis of the KEGG pathway in DEmRNAs is depicted in Figure 1C. The top 20 enriched pathways included “Endocytosis,” “Cytokine-cytokine receptor interaction,” and “PPAR signaling pathway,” which may be associated with inflammation development. The relationship between these pathways is reflected in Figure 1D.
The association in STRING includes direct (physical) and indirect (functional) interactions. The PPI network (Figure 3) showed that myeloid differentiation-factor 88 (MyD88) is a central protein and interacts with multiple proteins.
Figure 3. Protein–protein interaction (PPI) analysis of DEmRNAs. Map node size and color to degree, low values to small sizes and dark colors. Map edge size to the combined source, low values to small sizes.
Analysis of Differentially Expressed lncRNA
A total of 21 DElncRNAs were obtained (Supplementary File 3), of which 13 were up-regulated and 8 were down-regulated (P < 0.05 and | log2(fold-change) | > 1). Volcano plot showed the DElncRNAs between the control group and treatment groups (Figure 4A). The heatmap showed hierarchical clustering for DElncRNAs in Figure 4B. A total of 288 DElncRNA target genes were predicted. Figure 5 depicts 107 GO terms (Supplementary File 4) from the GO analysis. The “integral component of the membrane” term was the found to be most significant enrichment. Figure 4C depicts the enrichment analysis of the KEGG pathway in DElncRNA targets. Figure 8A shows two common target genes (Supplementary File 5) from DEmRNA and DElncRNA targets and KEGG enrichment analysis of these genes is shown in Figure 8C.
Figure 4. Screening and enrichment analysis of differently expressed lncRNAs (DElncRNAs) in S. aureus infected mammary epithelial cells compared with non-infected mammary epithelial cells. (A) Cluster analysis of DElncRNAs in mammary epithelial cells between the control group and treatment group. (B) Volcano plot of global DElncRNAs in mammary epithelial cells between control group and treatment group. (C) KEGG pathway classification for the annotation of DElncRNAs in mammary epithelial cells.
Summary of miRNA Sequencing
The quality of six miRNA transcriptome sequence expression profiles of mammary epithelial cells is presented in Table 3. The data used for downstream analysis was filtered.
A total of 27 DEmiRNAs were identified (P < 0.05 and | log2(fold-change) | > 1) (Supplementary File 6) and the volcano plot also shows these DEmiRNAs (Figure 6B). The heatmap of the DEGs shows clustering of DEmiRNAs (Figure 6A). Figure 7 shows the predicted 784 target genes of DEmiRNAs that enriched the 417 GO terms (Supplementary File 7). KEGG pathway enrichment analysis of DEmiRNAs is shown in Figure 6C, suggesting DEmiRNAs major role in the “Ras signaling pathway,” “Endocytosis,” and “PI3K-AKT signaling pathway”. There were 14 common target genes (BASP1, RAB11FIP4, PIP5K1C, VAMP4, DHRS12, LBH, MOB3B, ABO, FOXQ1, WDFY2, DYRK1B, NCKAP5L, TMEM59L, LIMK1) from DEmRNA and DEmiRNA targets Figure 8B. Figure 8D shows KEGG enrichment analysis of DEmRNA and DEmiRNA targets.
Figure 6. Screening and enrichment analysis of differently expressed miRNAs (DEmiRNAs) in S. aureus infected mammary epithelial cells compared with non-infected mammary epithelial cells. (A) Cluster analysis of DEmiRNAs in mammary epithelial cells between the control group and treatment group. (B) Volcano plot of global DEmiRNAs in mammary epithelial cells between the control group and treatment group. (C) KEGG pathway classified annotation of DEmiRNAs in mammary epithelial cells.
Figure 8. Analysis of DElncRNA target genes and DEmiRNA target genes. (A) Venn map analysis of DElncRNA target genes and DEmRNAs. (B) Venn map analysis of DEmiRNA target genes and DEmRNAs. (C) KEGG pathway classification for the annotation of two common target genes from the DEmRNAs and DElncRNA target genes in mammary epithelial cells. (D) KEGG pathway classification for the annotation of 14 common target genes from the DEmRNAs and DEmiRNA target genes in mammary epithelial cells.
Competitive Endogenous RNAs (ceRNAs) Analysis
The ceRNAs network is shown in Figure 9 which reflects the competitive endogenous relationships between DEmRNAs, DEmiRNAs, and DElncRNAs.
Confirmation of Gene Expression With qRT-PCR
Nine differential expressed genes (PIDD1, DRAM1, TNFRSF18, DEPDC5, TMEM102, TRIM32, CD36, CXCR1) were randomly selected and quantified using qRT-PCR to confirm differentially expressed genes in mammary epithelial cells obtained by RNA-seq between the control and treatment group. The gene symbol corresponding to all mRNAs is shown in Supplementary File 8. PIDD1, DRAM1, TNFRSF18, and DEPDC5 were the up-regulating genes while TMEM102, TRIM32, TNFRSF25, CD36, and CXCR1 were down-regulating in the high-throughput RNA-Seq. The expression of these genes was confirmed by qRT-PCR (Figure 10).
Figure 10. Expression levels of selected DEGs quantified by quantitative reverse transcription-PCR (qRT-PCR). GAPDH, ACTB, and PPIA were used as an internal control, and data are presented as fold change (N = 3 per group).
Discussion
Staphylococcus aureus often causes subclinical bovine mastitis and persistent intramammary infections. Ineffective pathogen clearance often leads to chronic and persistent infections (7). Studies have shown that S. aureus invades udder epithelial cells, which protects the pathogen from host defenses and antibiotics (6). Thus, it is particularly important to understand the response associated with bovine mastitis at a molecular level. In this study, we performed a comprehensive assessment of the whole transcriptomic profile of bMECs infected by S. aureus intracellularly, it contributed to a more embedded understanding of the transcriptome regulation behind this biological process.
Wang et al. investigated the transcriptional responses of primary bovine mammary epithelial cells against three S. aureus strains with different virulent factors using a tag-based high-throughput transcriptome sequencing technique (24). Li et al. identified functional miRNAs in bovine mammary glands infected with S. aureus on Chinese Holstein cows using Solexa sequencing and bioinformatics (25). Fang et al. compared the expression of mRNAs and miRNAs at after 24 h of intra-mammary infection (IMI) with high or low concentrations of S. aureus (26). The current study mainly focused on the defense mechanism of bMECs against intracellular S. aureus. Differences were observed from the previous research from the experimental designs to the conclusion.
Through the GO functional enrichment analysis for DEmRNAs, functional molecular processes were found to be related to the inflammatory response, immune response, peroxisome, regulation of autophagy, and positive regulation of ERK1 and ERK2 cascade. Related research showed that peroxisomes play an indispensable role in the generation of reactive oxygen species (ROS) (48). ROS is crucial in the signaling and defense of biological organisms and can also contribute to bacterial killing (49–51). Previous studies showed that S. aureus induces bMECs triggering immune responses and ROS production (52–54). Autophagy is a highly conserved mechanism that maintains homeostasis by removing damaged organelles and cytoplasmic proteins. It is also an immune response pathway that can eliminate intracellular bacteria (55–57). Bacterial infection induces autophagy and transports bacteria to the lysosome for degradation by autophagosomes (58, 59). The ERK1/2 signaling pathway has been shown to respond to the autophagy regulation of intracellular pathogens (60). KEGG pathway analysis identified several critical pathways, such as endocytosis, cytokine-cytokine receptor interaction, PPAR signaling pathway, and the inflammation mediator regulation of TRP channels. Endocytosis is essential for the entry of pathogens into cells and followed by its replication. Bacteria use this mechanism to utilize the osmotic network of the endocytic organelles to enter the cytoplasmic or replicating site and reach the relevant intracellular compartment (61, 62). Cytokine-cytokine receptor interaction is also reported to be involved in clinical mastitis (63). Peroxisome proliferator-activated receptor (PPAR) has anti-inflammatory effects inevitably linked to inflammation (64). Studies have shown that transient receptor potential (TRP) channels are associated with various factors and mechanisms that can activate/modulate inflammation through innate immunity (65, 66). They also play a key role in S. aureus elimination from bovine mammary epithelial cells.
In this study, some genes were up-regulated in the S. aureus-treated bMECs, such as ACOX2 (acyl-CoA oxidase 2), CROT (carnitine O-octanoyltransferase), NUDT12 (nudix hydrolase 12), and DRAM1 (DNA damage regulated autophagy modulator 1). Among these, three genes (ACOX2, CROT, and NUDT12) were enriched in peroxisomes. Additionally, peroxisomes play an integral role in the production of ROS. ROS levels acutely increased during cellular stress and the process of bacterial killing. Many studies supported the indispensable role of cellular stressors in regulating the innate immune responses (51). A new study reported that ROS can enhance macrophage antimicrobial activity against intracellular S. aureus (67). In the S. aureus treated group, DRAM1 expression was up-regulated. Moreover, in the case of DRAM1 over-expression, autophagosome production could be triggered by enriching DRAM1 on the Golgi membrane (68). Besides, DRAM1 can affect autophagy by the acidification of lysosomes, the fusion of lysosomes with autophagosomes, and autophagosome clearance (69). However, although autophagy has an antibacterial effect, there is evidence that pathogens have developed several ways to evade this mechanism (70). S. aureus can avoid autophagy clearance in bMECs by impairing lysosome function (30). Which may be an important mechanism of S. aureus to evade the host's immune response.
Conversely, certain genes including CD36 (CD36 molecule), CXCR1 (chemokine C-X-C motif receptor 1), BMP4 (bone morphogenetic protein 4), and TNFRSF25 (TNF receptor superfamily member 25) were down-regulated in S. aureus-treated bMECs. CD36 is a transmembrane glycoprotein receptor and contributes to the host's innate defense against S. aureus. It binds to TLR2 (toll like receptor 2) and recognizes the S. aureus cell wall diacylglycerol, inducing phagocytosis and cytokine production. In addition, CD36 expression on macrophages plays a significant role in host control of inflammation and skin damage during skin infection caused by S. aureus (71). Previous studies showed that blocking cell surface CXCR1 expression could be used as a beneficial treatment against S. aureus infection by disrupting the balance between inflammation and bacterial clearance (72). BMP4 facilitates leukocyte recruitment and inflammation improvement (73). The down-regulation of CD36, CXCR1, and BMP4 expression is not conducive to eliminating the bacteria but may be beneficial for S. aureus to evade its destruction.
Analysis of DElncRNA target profiles revealed that DElncRNAs participate in the regulation of the “TNF signaling pathway” and “NF-kappa B signaling pathway.” The pro-inflammatory cytokine TNF plays a significant role in apoptosis, cell proliferation, differentiation, necrosis, and cytokine induction (74). Activation of NF-kappa B is thought to guide the transcription of genes associated with inflammation and cell death or survival (75).
The integrated analysis identified DEmiRNA targets to be related to the “mTOR signaling pathway,” “Endocytosis,” and “PI3K-AKT signaling pathway”. PI3K regulates its downstream effector's activity through the AKT/mTOR/p70S6K signaling axis, thereby altering the production of cytokines, and thus could be a potential target for inflammatory diseases (76). Previous studies have also shown that Bta-miR-145 expression was down-regulated in udder tissues after S. aureus infection, consistent with our research. It was also found that overexpression of Bta-miR-145 significantly inhibited the proliferation of bMECs, and also reduced the secretion of IL-12 and TNF-α, but increased the secretion of IFN-γ (77). Therefore, the down-regulation of Bta-miR-145 is conducive to the production of IL-12 and TNF-α, thus inducing an immune response.
Conclusion
In the current study, we characterized the whole transcriptome profiles of bovine mammary epithelial cells infected intracellularly with S. aureus by RNA-seq. S. aureus invading into bovine mammary epithelial cells can trigger the immune responses, ROS production, and the expression of genes involved in autophagy. These differentially expressed RNAs may be critical in understanding the molecular mechanisms of S. aureus to survive in bovine mammary epithelial cells. It thus provides novel insights into the responses to bovine mammary epithelial cells with intracellular S. aureus.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.
Author Contributions
YL and BH: conceptualization. XW and NG: data curation and writing—original draft. XW and FS: formal analysis. XY: investigation. NG and XY: methodology. JL and YL: project administration, writing—review and editing, and supervision. RW: resources. XW: software. LL and RW: validation. MZ: visualization. All authors contributed to the article and approved the submitted version.
Funding
The project was supported by the National Science Foundation of China (31802259, 31872535), the Shandong Key R&D Program (2019GNC106141), the Shandong Natural Science Foundation of China (ZR2018MC027, ZR2016CQ25), the China Postdoctoral Science Foundation (2018M632704, 2019T120601), and the Funds of Shandong Double Tops Program.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2020.00642/full#supplementary-material
Supplementary File 1. Detailed information of DEmRNAs. Differently expressed mRNAs (DEmRNAs) in S. aureus infected mammary epithelial cells compared with non-infected mammary epithelial cells.
Supplementary File 2. Detailed information of GO enrichment of DEmRNAs. GO term function classification of DEmRNAs.
Supplementary File 3. Detailed information of DElncRNAs. Differently expressed lncRNAs (DElncRNAs) in S. aureus infected mammary epithelial cells compared with non-infected mammary epithelial cells.
Supplementary File 4. Detailed information of GO enrichment of DElncRNAs. GO term function classification of DElncRNAs.
Supplementary File 5. The two gene common target genes from DEmRNA and DElncRNA targets. Intersection of DEmRNA and DElncRNA.
Supplementary File 6. Detailed information of DEmiRNAs. Differently expressed miRNAs (DEmiRNAs) in S. aureus infected mammary epithelial cells compared with non-infected mammary epithelial cells.
Supplementary File 7. Detailed information of GO enrichment of DEmiRNAs. GO term function classification of DEmiRNAs.
Supplementary File 8. The gene symbol corresponding to all mRNAs.
References
1. Barkema H, Schukken Y, Lam T, Beiboer M, Wilmink H, Benedictus G, et al. Incidence of clinical mastitis in dairy herds grouped in three categories by bulk milk somatic cell counts. J Dairly Sci. (1998) 81:411–9. doi: 10.3168/jds.S0022-0302(98)75591-2
2. Halasa T, Huijps K, Østerås O, Hogeveen H. Economic effects of bovine mastitis and mastitis management: a review. Vet Q. (2007) 29:18–31. doi: 10.1080/01652176.2007.9695224
3. Zhao X, Lacasse P. Mammary tissue damage during bovine mastitis: causes and control. J Anim Sci. (2008) 86:57–65. doi: 10.2527/jas.2007-0302
4. Ruegg PL. A 100-year review: mastitis detection, management, and prevention. J Dairy Sci. (2017) 100:10381–97. doi: 10.3168/jds.2017-13023
5. Gruet P, Maincent P, Berthelot X, Kaltsatos V. Bovine mastitis and intramammary drug delivery: review and perspectives. Adv Drug Deliv Rev. (2001) 50:245–59. doi: 10.1016/S0169-409X(01)00160-0
6. Dego OK, Van Dijk J, Nederbragt H. Factors involved in the early pathogenesis of bovine Staphylococcus aureus mastitis with emphasis on bacterial adhesion and invasion. A review. Vet Q. (2002) 24:181–98. doi: 10.1080/01652176.2002.9695135
7. Anaya-López JL, Contreras-Guzmán OE, Cárabez-Trejo A, Baizabal-Aguirre VM, Lopez-Meza JE, Valdez-Alarcon JJ, et al. Invasive potential of bacterial isolates associated with subclinical bovine mastitis. Res Vet Sci. (2006) 81:358–61. doi: 10.1016/j.rvsc.2006.02.002
8. Gentilini E, Denamiel G, Llorente P, Godaly S, Rebuelto M, DeGregorio O. Antimicrobial susceptibility of Staphylococcus aureus isolated from bovine mastitis in Argentina. J Dairy Sci. (2000) 83:1224–7. doi: 10.3168/jds.S0022-0302(00)74988-5
9. Tenhagen B-A, Köster G, Wallmann J, Heuwieser W. Prevalence of mastitis pathogens and their resistance against antimicrobial agents in dairy cows in Brandenburg, Germany. J Dairy Sci. (2006) 89:2542–51. doi: 10.3168/jds.S0022-0302(06)72330-X
10. Idriss SE, Foltys V, Tančin V, Kirchnerová K, Tančinová D, Zaujec K. Mastitis pathogens and their resistance against antimicrobial agents in dairy cows in Nitra, Slovakia. Slovak J Anim Sci. (2014) 47:33–8. Available online at: https://sjas.ojs.sk/sjas/article/view/206
11. Gomes F, Henriques M. Control of bovine mastitis: old and recent therapeutic approaches. Curr Microbiol. (2016) 72:377–82. doi: 10.1007/s00284-015-0958-8
12. Bushati N, Cohen SM. microRNA functions. Annu Rev Cell Dev Biol. (2007) 23:175–205. doi: 10.1146/annurev.cellbio.23.090506.123406
13. Eulalio A, Schulte L, Vogel J. The mammalian microRNA response to bacterial infections. RNA Biol. (2012) 9:742–50. doi: 10.4161/rna.20018
14. Maudet C, Mano M, Eulalio A. MicroRNAs in the interaction between host and bacterial pathogens. FEBS Lett. (2014) 588:4140–7. doi: 10.1016/j.febslet.2014.08.002
15. Zhang P, Wang L, Li Y, Jiang P, Wang Y, Wang P, et al. Identification and characterization of microRNA in the lung tissue of pigs with different susceptibilities to PCV2 infection. Vet Res. (2018) 49:18. doi: 10.1186/s13567-018-0512-3
16. Aguilar C, Mano M, Eulalio A. MicroRNAs at the host–bacteria interface: host defense or bacterial offense. Trends Microbiol. (2019) 27:206–18. doi: 10.1016/j.tim.2018.10.011
17. Heward JA, Lindsay MA. Long non-coding RNAs in the regulation of the immune response. Trends Immunol. (2014) 35:408–19. doi: 10.1016/j.it.2014.07.005
18. Zur Bruegge J, Einspanier R, Sharbati S. A long journey ahead: long non-coding RNAs in bacterial infections. Front Cell Infect Microbiol. (2017) 7:95. doi: 10.3389/fcimb.2017.00095
19. Wang H, Liu H, Liu J, Zhao K, Wang C, Yang W. High-level exogenous trans10, cis12 conjugated linoleic acid plays an anti-lipogenesis role in bovine mammary epithelial cells. Anim Sci J. (2014) 85:744–50. doi: 10.1111/asj.12204
20. Chen YG, Satpathy AT, Chang HY. Gene regulation in the immune system by long noncoding RNAs. Nat Immunol. (2017) 18:962–72. doi: 10.1038/ni.3771
21. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. (2008) 5:621. doi: 10.1038/nmeth.1226
22. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. (2009) 10:57. doi: 10.1038/nrg2484
23. Soon WW, Hariharan M, Snyder M. High-throughput sequencing for biology and medicine. Mol Syst Biol. (2013) 9:640. doi: 10.1038/msb.2012.61
24. Wang X, Xiu L, Hu Q, Cui X, Liu B, Tao L, et al. Deep sequencing-based transcriptional analysis of bovine mammary epithelial cells gene expression in response to in vitro infection with Staphylococcus aureus stains. PLoS ONE. (2013) 8:e82117. doi: 10.1371/journal.pone.0082117
25. Li R, Zhang CL, Liao XX, Chen D, Wang WQ, Zhu YH, et al. Transcriptome microRNA profiling of bovine mammary glands infected with Staphylococcus aureus. Int J Mol Sci. (2015) 16:4997–5013. doi: 10.3390/ijms16034997
26. Fang L, Hou Y, An J, Li B, Song M, Wang X, et al. Genome-wide transcriptional and post-transcriptional regulation of innate immune and defense responses of bovine mammary gland to Staphylococcus aureus. Front Cell Infect Microbiol. (2016) 6:193. doi: 10.3389/fcimb.2016.00193
27. Wang XG, Ju ZH, Hou MH, Jiang Q, Yang CH, Zhang Y, et al. Deciphering transcriptome and complex alternative splicing transcripts in mammary gland tissues from cows naturally infected with Staphylococcus aureus mastitis. PLoS ONE. (2016) 11:e0159719. doi: 10.1371/journal.pone.0159719
28. Kosciuczuk EM, Lisowski P, Jarczak J, Majewska A, Rzewuska M, Zwierzchowski L, et al. Transcriptome profiling of staphylococci-infected cow mammary gland parenchyma. BMC Vet Res. (2017) 13:161. doi: 10.1186/s12917-017-1088-2
29. Rejman J, Oliver S, Muenchen R, Turner J. Proliferation of the MAC-T bovine mammary epithelial cell line in the presence of mammary secretion whey proteins. Cell Biol Int Rep. (1992) 16:993–1001. doi: 10.1016/S0309-1651(06)80052-4
30. Geng N, Wang X, Yu X, Wang R, Zhu Y, Zhang M, et al. Staphylococcus aureus avoids autophagy clearance of bovine mammary epithelial cells by impairing lysosomal function. Front Immunol. (2020) 11:746. doi: 10.3389/fimmu.2020.00746
31. Karolchik D, Baertsch R, Diekhans M, Furey TS, Hinrichs A, Lu YT, et al. The UCSC genome browser database. Nucleic Acids Res. (2003) 31:51–4. doi: 10.1093/nar/gkg129
32. Zhao Y, Li H, Fang S, Kang Y, Wu W, Hao Y, et al. NONCODE 2016: an informative and valuable data source of long non-coding RNAs. Nucleic Acids Res. (2016) 44:D203–8. doi: 10.1093/nar/gkv1252
33. Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. (2006) 34:D140–4. doi: 10.1093/nar/gkj112
34. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. (2016) 11:1650–67. doi: 10.1038/nprot.2016.095
35. Langmead B. Aligning short sequencing reads with Bowtie. Curr Protoc Bioinformatics. (2010) 32:11.7.1–14. doi: 10.1002/0471250953.bi1107s32
36. Mackowiak SD. Identification of novel and known miRNAs in deep-sequencing data with miRDeep2. Curr Protoc Bioinformatics. (2011) 36:12.10.11–15. doi: 10.1002/0471250953.bi1210s36
37. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
38. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, et al. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic Acids Res. (2016) 45:D362–8. doi: 10.1093/nar/gkw937
39. Alkan F, Wenzel A, Palasca O, Kerpedjiev P, Rudebeck AF, Stadler PF, et al. RIsearch2: suffix array-based large-scale prediction of RNA-RNA interactions and siRNA off-targets. Nucleic Acids Res. (2017) 45:e60. doi: 10.1093/nar/gkw1325
40. Ritchie W, Flamant S, Rasko JE. Predicting microRNA targets and functions: traps for the unwary. Nat Methods. (2009) 6:397–8. doi: 10.1038/nmeth0609-397
41. Sticht C, De La Torre C, Parveen A, Gretz N. miRWalk: an online resource for prediction of microRNA binding sites. PLoS ONE. (2018) 13:e0206239. doi: 10.1371/journal.pone.0206239
42. Dennis G Jr, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, et al. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. (2003) 4:P3. doi: 10.1186/gb-2003-4-5-p3
43. Yu GC, Wang LG, Han YY, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics J Integr Biol. (2012) 16:284–7. doi: 10.1089/omi.2011.0118
44. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell. (2011) 146:353–8. doi: 10.1016/j.cell.2011.07.014
45. Denzler R, Agarwal V, Stefano J, Bartel DP, Stoffel M. Assessing the ceRNA hypothesis with quantitative measurements of miRNA and target abundance. Mol Cell. (2014) 54:766–76. doi: 10.1016/j.molcel.2014.03.045
46. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
47. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2– ΔΔCT method. Methods. (2001) 25:402–8. doi: 10.1006/meth.2001.1262
48. Wang B, Apanasets O, Nordgren M, Fransen M. Dissecting peroxisome- mediated signaling pathways: a new and exciting research field. In: Brocard C, Hartig A, editors. Molecular, Machines Involved in Peroxisome Biogenesis and Maintenance. Vienna: Springer (2014). p. 255–73.
49. West AP, Brodsky IE, Rahner C, Woo DK, Erdjument-Bromage H, Tempst P, et al. TLR signalling augments macrophage bactericidal activity through mitochondrial ROS. Nature. (2011) 472:476. doi: 10.1038/nature09973
50. Fransen M, Nordgren M, Wang B, Apanasets O. Role of peroxisomes in ROS/RNS-metabolism: implications for human disease. Biochim Biophys Acta. (2012) 1822:1363–73. doi: 10.1016/j.bbadis.2011.12.001
51. Holmström KM, Finkel T. Cellular mechanisms and physiological consequences of redox-dependent signalling. Nat Rev Mol Cell. (2014) 15:411. doi: 10.1038/nrm3801
52. Rainard P, Riollet C. Innate immunity of the bovine mammary gland. Vet Res. (2006) 37:369–400. doi: 10.1051/vetres:2006007
53. Gilbert FB, Cunha P, Jensen K, Glass EJ, Foucras G, Robert-Granié C, et al. Differential response of bovine mammary epithelial cells to Staphylococcus aureus or Escherichia coli agonists of the innate immune system. Vet Res. (2013) 44:40. doi: 10.1186/1297-9716-44-40
54. Zheng L, Xu Y, Lu J, Liu M, Dai B, Miao J, et al. Variant innate immune responses of mammary epithelial cells to challenge by Staphylococcus aureus, Escherichia coli and the regulating effect of taurine on these bioprocesses. Free Radic Biol Med. (2016) 96:166–80. doi: 10.1016/j.freeradbiomed.2016.04.022
55. Klionsky DJ, Cuervo AM, Dunn J, William A, Levine B, van der Klei IJ, et al. How shall I eat thee?. Autophagy. (2007) 3:413–16. doi: 10.4161/auto.4377
56. Shahnazari S, Brumell J. Mechanisms and consequences of bacterial targeting by the autophagy pathway. Curr Opin Microbiol. (2011) 14:68–75. doi: 10.1016/j.mib.2010.11.001
57. Ktistakis NT, Tooze SA. Digesting the expanding mechanisms of autophagy. Trends Cell Biol. (2016) 26:624–35. doi: 10.1016/j.tcb.2016.03.006
58. He C, Klionsky DJ. Regulation mechanisms and signaling pathways of autophagy. Annu Rev Genet. (2009) 43:67–93. doi: 10.1146/annurev-genet-102808-114910
59. Bah A, Vergne I. Macrophage autophagy and bacterial infections. Front Immunol. (2017) 8:1483. doi: 10.3389/fimmu.2017.01483
60. Shinojima N, Yokoyama T, Kondo Y, Kondo S. Roles of the Akt/mTOR/p70S6K and ERK1/2 signaling pathways in curcumin-induced autophagy. Autophagy. (2007) 3:635–7. doi: 10.4161/auto.4916
61. Lonhienne TG, Sagulenko E, Webb RI, Lee K-C, Franke J, Devos DP, et al. Endocytosis-like protein uptake in the bacterium Gemmata obscuriglobus. Proc Natl Acad Sci USA. (2010) 107:12883–8. doi: 10.1073/pnas.1001085107
62. Cossart P, Helenius A. Endocytosis of viruses and bacteria. Cold Spring Harb Perspect Biol. (2014) 6:a016972. doi: 10.1101/cshperspect.a016972
63. Tiezzi F, Parker-Gaddis KL, Cole JB, Clay JS, Maltecca C. A genome-wide association study for clinical mastitis in first parity US Holstein cows using single-step approach and genomic matrix re-weighting procedure. PLoS ONE. (2015) 10:e0114919. doi: 10.1371/journal.pone.0114919
64. Murphy GJ, Holder JC. PPAR-γ agonists: therapeutic role in diabetes, inflammation and cancer. Trends Pharmacol. (2000) 21:469–74. doi: 10.1016/S0165-6147(00)01559-5
65. Parenti A, De Logu F, Geppetti P, Benemei SJ. What is the evidence for the role of TRP channels in inflammatory and immune cells? Br J Pharmacol. (2016) 173:953–69. doi: 10.1111/bph.13392
66. Khalil M, Alliger K, Weidinger C, Yerinde C, Wirtz S, Becker C, et al. Functional role of transient receptor potential channels in immune cells and epithelia. Front Immunol. (2018) 9:174. doi: 10.3389/fimmu.2018.00174
67. Abuaita BH, Schultz TL, O'Riordan MX. Mitochondria-derived vesicles deliver antimicrobial reactive oxygen species to control phagosome-localized Staphylococcus aureus. Cell Host Microbe. (2018) 24:625–36. e5. doi: 10.1016/j.chom.2018.10.005
68. Nagata M, Arakawa S, Yamaguchi H, Torii S, Endo H, Tsujioka M, et al. Dram1 regulates DNA damage-induced alternative autophagy. Cell Stress. (2018) 2:55. doi: 10.15698/cst2018.03.127
69. Zhang XD, Qi L, Wu JC, Qin ZH. DRAM1 regulates autophagy flux through lysosomes. PLoS ONE. (2013) 8:e63245. doi: 10.1371/journal.pone.0063245
70. Huang J, Brumell JH. Bacteria-autophagy interplay: a battle for survival. Nat Rev Microbiol. (2014) 12:101–14. doi: 10.1038/nrmicro3160
71. Castleman MJ, Febbraio M, Hall PR. CD36 Is Essential for regulation of the host innate response to Staphylococcus aureus α-toxin–mediated dermonecrosis. J Immunol. (2015) 195:2294–302. doi: 10.4049/jimmunol.1500500
72. Dutta P, Bishayi B. Neutralization of TNF-α and IL-1β regulates CXCL8 production through CXCL8/CXCR1 Axis in macrophages during Staphylococcus aureus Infection. Immunol Invest. (2020) 1–26. doi: 10.1080/08820139.2020.1787436
73. Vestweber D. Relevance of endothelial junctions in leukocyte extravasation and vascular permeability. Ann N Y Acad Sci. (2012) 1257:184–92. doi: 10.1111/j.1749-6632.2012.06558.x
74. Rothe J, Gehr G, Loetscher H, Lesslauer W. Tumor necrosis factor receptors–structure and function. Immunol Res. (1992) 11:81–90. doi: 10.1007/BF02918612
75. Wang T, Zhang X, Li JJ. The role of NF-kappaB in the regulation of cell stress responses. Int Immunopharmacol. (2002) 2:1509–20. doi: 10.1016/S1567-5769(02)00058-9
76. Xie S, Chen M, Yan B, He X, Chen X, Li D. Identification of a role for the PI3K/AKT/mTOR signaling pathway in innate immune cells. PLoS ONE. (2014) 9:e94496. doi: 10.1371/journal.pone.0094496
Keywords: Staphylococcus aureus, bovine mastitis, bovine mammary epithelial cells, transcriptome, microRNA, LncRNA 3
Citation: Wang X, Su F, Yu X, Geng N, Li L, Wang R, Zhang M, Liu J, Liu Y and Han B (2020) RNA-Seq Whole Transcriptome Analysis of Bovine Mammary Epithelial Cells in Response to Intracellular Staphylococcus aureus. Front. Vet. Sci. 7:642. doi: 10.3389/fvets.2020.00642
Received: 03 June 2020; Accepted: 07 August 2020;
Published: 07 December 2020.
Edited by:
Jasim Muhammad Uddin, Bangladesh Agricultural University, BangladeshReviewed by:
Md. Aminul Islam, Tohoku University, JapanViswas K. Nagaleekar, Indian Veterinary Research Institute (IVRI), India
Copyright © 2020 Wang, Su, Yu, Geng, Li, Wang, Zhang, Liu, Liu and Han. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yongxia Liu, liuyongxia@sdau.edu.cn; Jianzhu Liu, liujz@sdau.edu.cn
†These authors have contributed equally to this work