1. Introduction
Udders are vital functional organs in cows and provide a significant source of protein (milk) for humans. Milk is produced by the mammary glands and is collected into the milk ducts through physical methods such as squeezing and excreting[
1]. Tyoically, after the mammary gland is infected with pathogenic microorganisms, or stimulated by physical, chemical, and other factors, inflammatory changes can occur in the mammary plasma or parenchymal tissue, resulting in mastitis. This disease is primarily caused by various pathogenic bacteria or microorganisms in the environment infecting the mammary glands via milk ducts, leading to inflammation[
2]. Based on clinical symptoms, mastitis can be divided into clinical and subclinical forms [
3,
4]. Clinical mastitis is characterized by pronounced signs of inflammation in the udder, microbial presence, and chemical changes in milk[
5]. Cows with clinical mastitis are affected by pathogenic bacteria that cause fibrosis and atrophy of the mammary gland, leading to premature culling[
6]. However, cows with subclinical mastitis exhibit only mild inflammation after pathogenic bacteria enter teat duct. The somatic cell count (SSC) in milk was increases, with no other clinical manifestations[
7]. Although subclinical mastitis does not pose an immediate risk, if left untreated, it will eventually develop into clinical mastitis and cause huge economic losses.
The study by Paramanandham et al. indicated that
Staphylococcus aureus is the most common mastitis pathogen,
Escherichia coli is the main pathogen causing clinical mastitis and
Streptococcus spp. can cause both subclinical and clinical mastitis worldwide[
8]. In many cases of subclinical mastitis, the inflammatory response is predominantly caused by
Streptococcus agalactiae[
9,
10,
11].
S. agalactiae is classified as group B according to the Lancefield bacterial taxonomy. The bacterium is mostly β-hemolytic, although some non-hemolytic and CAMP-positive strains have been observed[
12]. Due to the variety of virulence factors,
S. agalactiae can adhere to and invade host cells, inducing an inflammatory response in the mammary gland. Therefore, it is essential to analyze the transcriptional regulation of inflammatory genes during
S. agalactiae invasion of the mammary gland. The mammary epithelium in first line of defense in the mammary gland. It can effectively initiate an immune response, eliminating pathogens before any abnormal changes occur in the mammary gland, which is crucial for resistance to mastitis and affects susceptibility[
13]. When the pathogen successfully breaches the host's physical defenses, mammary gland epithelial detect bacteria through specific pattern recognition receptors and initiate a series of immune responses. Currently, researchers have conducted high-throughput sequencing studies on mammary infected with
S. agalactiae and identified several immune-related receptors or pathways, but the findings are not consistent. However, the mechanisms underlying
S. agalactiae infection in subclinical mastitis are not well understood.
Zhang et al. found that, compared with healthy cow mammary glands, a total of 129 differentially expressed genes and 144 differentially expressed proteins were identified in mammary glands infected with
S. agalactiae[
14]. Intramammary infection with
S. agalactiae triggered a complex host innate immune response involving complement and coagulation cascades, ECM-receptor interaction, focal adhesion, and phagosome and bacterial invasion of epithelial cells pathways. Tong et al. used proteomics to discover that differentially expressed proteins included enzymes and proteins associated with various metabolic processes and cellular immunity in
S. agalactiae-infected bovine mammary epithelial cells[
15]. Subsequently, they used the ubiquitinome analysis to find that ubiquitinated proteins are associated with regulating cell junctions in the host[
16]. Sbardella et al. infected dairy cow mammary glands with
S. agalactiae via manual intervention, and screened 122 differential genes in the sequencing data using three different statistical methods, but only platelet activation showed a significant enrichment pathway[
17]. Additionally, it has been found that differential genes in mammary alveolar tissue infected with
S. agalactiae are mainly involved in innate immune response, inflammatory response, chemokine signaling, Wnt signaling, complement and coagulation cascades compared with normal tissues[
18]. Richards et al. concluded that lactose metabolism is an important metabolic pathway for
S. agalactiae to adapt to the bovine mammary environment, as determined by sequencing analysis of isolated
S. agalactiae[
19]. In mammary glands infected with
S. agalactiae, differentially expressed miRNAs were mainly involved in the RIG-I-like receptor signaling pathway, cytosolic DNA sensing pathway, and Notch signaling pathway[
20].
Most of the studies are in vivo, and research on immune changes following mammary epithelial cell infection is lacking. Therefore, studying the immune regulation of mammary epithelial cells infected with
S. agalactiae is essential. The transcriptome serves as a powerful indicator of the physiological state of a cell, whether healthy or diseased. Consequently, transcriptome analysis has become a crucial tool for understanding the molecular changes that occur during bacterial infections of eukaryotic cells. Previously, transcriptomic studies were technically limited to analyzing mRNA expression changes in either the bacterial pathogen or the infected eukaryotic host cell. However, the increasing sensitivity of high-throughput RNA sequencing now enables "Dual RNA-seq" studies, simultaneously capturing all classes of coding and noncoding transcripts in both the pathogen and the host[
21].
To the best of our knowledge, this study represents the first report of Dual RNA-seq on mammary epithelial cells infected with S. agalactiae. Identifying key differentially expressed genes and pathways in infected mammary epithelial cells provides a foundation for a better understanding of the central mechanisms of host defense during subclinical infections like mastitis.
4. Discussion
S. agalactiae, commonly referred to as Group B
Streptococcus (GBS), is a zoonotic pathogen and a highly infectious Gram-positive bacterium[
28]. This pathogen can infect the mammary glands of dairy cows, leading to subclinical mastitis, which can progress to clinical mastitis if left unchecked[
28]. Research has demonstrated that
S. agalactiae infection in the mammary gland initiates a series of innate immune responses in the host. Concurrently, the bacteria invade mammary epithelial cells to evade host defenses and antibiotics[
14]. During this infection process, various molecular interactions occur, highlighting the importance of studying the molecular changes that mediate mastitis.
Zhang et al. utilized transcriptomics and proteomics to analyze mammary tissue infected by
S. agalactiae, investigating the host’s immune response to the pathogen[
14]. Tong et al. conducted ubiquitination sequencing and analysis on mammary epithelial cells infected with
S. agalactiae[
16]. Mayara et al. perfused mammary tissue with
S. agalactiae to observe transcriptional changes, focusing on the most affected biological functions and pathways during this process[
18]. These studies primarily examined the resistance mechanisms of mammary tissue or the ubiquitination processes in mammary epithelial cells following
S. agalactiae infection, each emphasizing different aspects. This study utilized in vitro cultured mammary epithelial cells infected with varying times of
S. agalactiae and performed interaction transcriptome analysis to investigate the molecular mechanisms of host-pathogen interactions during transcriptional regulation.
The GO enrichment analysis of DEmRNAs showed enrichment in processed related to inflammation, disease occurrence, damage repair, and apoptosis regulation. Upregulated genes were enriched in exogenous apoptosis regulation processes, while downregulated genes were enriched in endogenous apoptosis signaling pathways. This could be due to the induction of exogenous apoptotic signals and the inhibition of endogenous apoptotic processes following pathogen invasion. Upregulated genes were significantly enriched in small GTPase-mediated signal transduction and regulation of the ERK1 and ERK2 cascades.
Rho-GTPase acts as a molecular switch during inflammatory cell migration by cycling between the inactive Rho-GDP and active Rho-GTP forms. It plays a crucial role in actin cytoskeleton dynamics and precise regulation of leukocyte immune functions. Existing reports indicate that dysregulation of Rho-GTPase signaling is associated with various inflammatory diseases[
29]. Small GTPase can mediate ERK1/2 to enter the signal nucleus through the MAPK signaling pathway, leading to apoptosis, inflammatory stress responses, and reactive oxygen species (ROS) production[
30]. ROS acts as a double-edged sword, potentially playing roles in both pro-inflammatory and anti-inflammatory processes.
Recent studies have revealed the physiological importance of ROS as signaling molecules, crucial for maintaining cellular function and homeostasis[
31]. Previous research indicated that
S. agalactiae infection in human endothelial cells induces ROS[
32], which can persist for a week compared to
Staphylococcus aureus and
Escherichia coli[
33]. In this study, we found that genes related to the regulation of reactive oxygen species metabolic processes were upregulated as mammary epithelial cells were increasingly infected by the pathogen. This indicates that mammary epithelial cells employ ROS signaling as a defensive measure during infection. However, ROS production further activates the ERK1/2 pathway, triggering various immune processes to eliminate cells, but also leading to apoptosis[
34,
35]. ROS can also induce autophagy, allowing the pathogen to survive in mammary epithelial cells[
36]. These findings suggest that during
S. agalactiae infection in mammary epithelial cells, the ERK1/2 pathway induces inflammatory responses through MAPK signaling, accompanied by increased ROS production.
KEGG analysis identified several crucial pathways, including transcription misregulation in cancer, IL-17 signaling pathway, P53 signaling pathway, TGF-beta signaling pathway, MAPK signaling pathway, pathways in cancer, and apoptosis. Transcription misregulation in cancer pathways is triggered by pathogen invasion, interferes with the regulation of cancer-related transcription factors. The upregulated
RUNX1 can inhibit the invasiveness of most breast cancer subtypes, especially in the early stages of tumorigenesis, and prevent the epithelial-to-mesenchymal transition in breast cancer cells[
32]. In time-series analysis,
RUNX1 expression gradually increases with the severity of bacterial infection, likely due to bacterial interference and disruption of the cell's regulatory system. Recent studies have shown that the IL-17 signaling pathway is also involved in the occurrence of mastitis[
37], particularly related to the inflammatory response in mammary epithelial cells[
38]. The p53 signaling pathway and MAPK signaling pathway are enriched in exosomes from bacterial-infected cells[
39]. The p53 signaling pathway is a complex cellular stress response network with various inputs and downstream outputs, related to its role as a tumor suppressor pathway[
40]. The MAPK signaling pathway is involved in tumor formation, invasion, metastasis, and apoptosis[
41], and its activation is also implicated in mastitis[
42]. The TGF-beta signaling pathway induces apoptosis in mammary epithelial cells[
43]. Moreover, TGF-beta1 can cooperate with the ERK1/2 pathway to promote Gram-positive bacterial adhesion and infection of mammary epithelial cells[
44]. Additionally, the RhoA/Rho kinase signaling cascade aids in TGF-beta induced changes in cytoskeletal organization and cell permeability[
45].
BCL2L11, involved in the apoptosis pathway, primarily participates in the extrinsic apoptotic signal pathway in the absence of ligand, positive regulation of cysteine-type endopeptidase activity involved in the apoptotic process, and protein kinase binding among upregulated genes. It is also involved in EGFR tyrosine kinase inhibitor resistance.
It is conceivable that upon pathogens interference, host cells initiate immune and inflammatory responses to combat bacterial infection, potentially leading to adverse effects such as progression toward cancer. Enriched signaling pathways in the upregulated clusters 1 and 3 suggest that cells may gradually transform towards a cancerous state. However, cells cannot evade the process of apoptosis. In Cluster 3, several disease-related pathways are upregulated, including the pathway in cancer. Concurrently, BCL2L11 is enriched in this pathway along with other genes.
Alternative splicing is a crucial mechanism of genetic regulation, enhancing the diversity and complexity of the transcriptome and proteome from a limited number of genes. Numerous studies suggest that alternative splicing events can lead to changes in protein expression or function during disease onset and progression[
46]. In this study, variable splicing showed that DEmRNAs were mainly involved in regulating GTPase activity, protein ubiquitination, ATP binding, and RNA binding. Some studies suggest that immune-related GTPase directly mediates pathogen membrane destruction by binding to the pathogen membrane, thereby exposing the pathogen to cytoplasmic defenses[
47]. Immune-related GTPase can also utilize the ubiquitination system to tag intracellular pathogens[
48]. In this study, differentially spliced genes were involved in ubiquitin-dependent protein catabolic processes and protein polyubiquitination. This indicates that
S. agalactiae infection in mammary epithelial cells activates the host cell's immune defense mechanisms, leading to the overactivation of the GTPase system, which collaborates with the intracellular ubiquitination system to combat intracellular bacterial damage. However, the influence of bacterial intracellular molecular regulation disrupts normal transcriptional regulation in host cells, particularly the normal alternative splicing process, potentially causing dysregulation of GTPase signaling. Differentially spliced genes were enriched in ATP binding and RNA binding functions, indicating that
S. agalactiae interferes with the splicing forms of ATP binding proteins and RNA binding proteins, indirectly affecting cellular energy metabolism and transcriptional regulation.
Among the three groups of DEmRNAs, the upregulated gene
RUNX1 exhibited MXE, A3SS, and SE events in the S_M group (all not significant), A3SS and SE events in the H_M group (not significant), and significant A3SS and SE events in the H_S group. KEGG analysis of the H_S group indicated that
RUNX1 is primarily involved in cancer-related pathways, including cancer transcription misregulation pathways and tight junctions. GO analysis revealed involvement in nucleoplasm localization, ATP binding, and protein-containing complex processes. This suggests that as bacterial infection intensifies, the splicing process of
RUNX1 in host cells is disrupted, altering its expression pattern and potentially increasing the risk of cellular transformation. Studies have demonstrated that
RUNX1 is associated with RNA Pol II-transcribed protein and lncRNA genes, as well as RNA Pol I-transcribed ribosomal genes, which are crucial for the growth and phenotype maintenance of mammary epithelial cells[
49]. This further implies that bacterial interference with the transcription process indirectly induce morphological changes in mammary epithelial cells, thereby contributing to carcinogenesis. Recent research indicates that RUNX1 plays a role in breast cancer migration and invasion[
50].
In this study, through the analysis of DElncRNA target gene prediction, alternatively spliced genes, and DEmRNA, which identified eight candidate genes at the intersection, including
RUNX1 and
BCL2L11.
RUNX1 is a crucial transcription factor that induces the expression of numerous genes. Among the upregulated genes,
BCL2L11 is one that
RUNX1 induces[
51].
BCL2L11 plays a crucial dual role in disease mechanisms by inhibiting autophagy and initiating apoptosis[
52]. Under normal conditions,
BCL2L11 undergoes alternative splicing, forming at least 18 different isoforms[
51]. In this study, however,
BCL2L11 underwent a MXE event between the H and M groups. MXE results in different exon combinations that may maintain protein folding but alter the specificity and selectivity of protein functions[
53]. This indicates that under continuous bacterial infection, the alternative splicing of
BCL2L11 is affected, altering its splicing form and resulting in changes to
BCL2L11 conformation. This splicing pattern reduces proteome diversity, leading to protein dysfunction and altered biological functions.
In the time series analysis, genes involved in BP of DNA binding transcription activator activity are upregulated in clusters 1 and 6. However, the transcription factors of the host cell are downregulated in clusters 5 and 8. This indicates that some factors from the pathogens have replaced the host's transcription factors, thereby affecting the host's transcriptional regulation. It is hypothesized that during pathogen infection of host cells, the host's spliceosome is disrupted and consequently downregulated. The host's cell nucleus is similarly affected, which interferes with the host's transcriptional regulation, particularly the activation of nucleic acid transcription factors and the specificity of RNA polymerase.
Interestingly, in the GO enrichment analysis of DEmRNAs in the host cell transcriptome, both spliceosome and mRNA splicing processes were downregulated, indicating that bacterial interference affects the normal gene splicing process of the mammary epithelial cell. Conversely, some proteins secreted by
S. agalactiae are involved. Combined with the differential expression analysis of the pathogenic bacteria, multiple genes are upregulated during infection, suggesting an influence on the host cell molecular transcriptional regulation and their interaction during infection. For example,
nusB,
rimN,
yhbY,
infC,
prfB, and
tsf are involved. Studies have shown that the
nusB is transcribed in opposition to the eukaryotic system and can be a potential antibacterial target[
54]. Additionally, it utilizes biological coupling of transcription and translation to control downstream gene expression[
55]. The protein YrdC translated from
rimN is involved in tRNA modification and preferentially binds to RNA[
56,
57].
yhbY is involved in ribosome assembly and exhibits RNA binding activity[
58].
infC guides transcriptional regulation and is upregulated in bacterial infection in the animal liver[
59,
60]. The RF2 protein encoded by
prfB is required for recognizing stop codons during bacterial translation termination[
61]. This suggests that bacteria proliferate extensively and express proteins capable of invading host cells.
tsf can contribute to the production of biologically active bacterial keratinase[
62], and this site can also confer strong antibiotic resistance[
63], making it a potential therapeutic target. To adapting to the host system, bacteria employ various strategies, including the production of virulence factors and the formation of biofilms, to evade the host immune system and resist antibiotics[
64].
Bacteria can manipulate the host signaling pathways by regulating the host lncRNAs to escape immune clearance. Therefore, bacteria also can induce significant alteration in the cell transcriptome and develop various strategies to modify immune signaling for its survival[
65]. Currently, lncRNAs have been shown to play a crucial role in regulating alternative splicing in response to various stimuli or diseases[
66]. Additionally, increasing evidence indicates that lncRNAs are important in the regulatory circuits controlling innate and adaptive immune responses to bacterial pathogens[
67]. In the co-expression network analysis, these three bacterial genes (
tsf,
prfB, and
infC) have the most significant effect on the differential targeting of lncRNAs to
RUNX1 and
BCL2L11. Both
RUNX1 and
BCL2L11 undergo abnormal alternative splicing during the infection process. It is hypothesized that during the infection of bovine mammary epithelial cells by
S. agalactiae, genes such as
tsf,
prfB, and
infC, which are involved in RNA binding, infiltrate host cells and disrupt the targeting of lncRNAs to
RUNX1 and
BCL2L11. This affects the host's normal alternative splicing process, disrupting normal cell proliferation regulation and apoptosis.
Figure 1.
Figure 1 Screening differently expressed mRNAs (DEmRNAs) in
S. agalactiae infected mammary epithelial cells among the control (n= 5),
S. agalactiae – S groups (n= 5), and
S. agalactiae - H groups (n= 3). A. Gene expression level analysis in the S_M, H_M, H_S. The X-axis of the box plot represents the sample name, while the Y-axis represents log10 (RPKM). The box plots for each region correspond to five statistical measures (maximum, upper quartile, median, lower quartile, and minimum values, respectively). B. Cluster analysis of DEmRNAs in mammary epithelial cells between the control group (M1, M2, M3, M4, and M5), S groups (S1, S2, S3, S4, and S5) and H treated groups (H1, H2, and H3). 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. 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 is and the trend of the more recent gene expression. C. D. E. Volcano plot of global DEmRNAs in S_M, H_M and H_S, respectively. Red dots (Up) represent significantly up-regulated genes (p-values< 0.05, log2(fold-change) >1); blue dots (Down) represent significantly down-regulated genes (p-values< 0.05, log2(fold-change) <-1); grey dots represent insignificantly differential expressed genes. F. Upset map analysis of S_M, H_M and H_S. The origin and connecting lines of the X-axis represent intersections, while the black bars represent the number of differentially expressed genes in each group; The number of differentially expressed genes at the intersection of each group on the Y-axis.
Figure 1.
Figure 1 Screening differently expressed mRNAs (DEmRNAs) in
S. agalactiae infected mammary epithelial cells among the control (n= 5),
S. agalactiae – S groups (n= 5), and
S. agalactiae - H groups (n= 3). A. Gene expression level analysis in the S_M, H_M, H_S. The X-axis of the box plot represents the sample name, while the Y-axis represents log10 (RPKM). The box plots for each region correspond to five statistical measures (maximum, upper quartile, median, lower quartile, and minimum values, respectively). B. Cluster analysis of DEmRNAs in mammary epithelial cells between the control group (M1, M2, M3, M4, and M5), S groups (S1, S2, S3, S4, and S5) and H treated groups (H1, H2, and H3). 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. 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 is and the trend of the more recent gene expression. C. D. E. Volcano plot of global DEmRNAs in S_M, H_M and H_S, respectively. Red dots (Up) represent significantly up-regulated genes (p-values< 0.05, log2(fold-change) >1); blue dots (Down) represent significantly down-regulated genes (p-values< 0.05, log2(fold-change) <-1); grey dots represent insignificantly differential expressed genes. F. Upset map analysis of S_M, H_M and H_S. The origin and connecting lines of the X-axis represent intersections, while the black bars represent the number of differentially expressed genes in each group; The number of differentially expressed genes at the intersection of each group on the Y-axis.
Figure 2.
GO and KEGG analysis of DEmRNAs in the S_M, H_M, H_S. A. The Y- axis on the left represents GO terms of up-regulated genes, including biological process (BP), cellular component (CP), and molecular function (MF), the X- axis indicates different comparison groups. The area of a circle represents DEGs number. Low p- values are shown in the red circle, and high p-values are shown in the blue circle. B. The Y- axis on the left represents GO terms of down-regulated genes, including biological process (BP), cellular component (CP), and molecular function (MF), the X- axis indicates different comparison groups. The area of a circle represents DEGs number. Low p-values are shown in the red circle, and high p-values are shown in the blue circle. C. The Y- axis on the left represents KEGG pathways, the X- axis indicates genes enrichment of each term. The shapes represent different groups. The area of shapes represents DEmRNAs numbers.
Figure 2.
GO and KEGG analysis of DEmRNAs in the S_M, H_M, H_S. A. The Y- axis on the left represents GO terms of up-regulated genes, including biological process (BP), cellular component (CP), and molecular function (MF), the X- axis indicates different comparison groups. The area of a circle represents DEGs number. Low p- values are shown in the red circle, and high p-values are shown in the blue circle. B. The Y- axis on the left represents GO terms of down-regulated genes, including biological process (BP), cellular component (CP), and molecular function (MF), the X- axis indicates different comparison groups. The area of a circle represents DEGs number. Low p-values are shown in the red circle, and high p-values are shown in the blue circle. C. The Y- axis on the left represents KEGG pathways, the X- axis indicates genes enrichment of each term. The shapes represent different groups. The area of shapes represents DEmRNAs numbers.
Figure 3.
GO enrichment result of differentially spliced genes in the S_M, H_M, H_S. The Y-axis on the left represents GO terms, including biological process (BP), cellular component (CP), and molecular function (MF), and the X-axis indicates gene enrichment of each term. Low p-values are shown in the red circle, and high p-values are shown in the green circle. The area of a circle represents of DEmRNAs number.
Figure 3.
GO enrichment result of differentially spliced genes in the S_M, H_M, H_S. The Y-axis on the left represents GO terms, including biological process (BP), cellular component (CP), and molecular function (MF), and the X-axis indicates gene enrichment of each term. Low p-values are shown in the red circle, and high p-values are shown in the green circle. The area of a circle represents of DEmRNAs number.
Figure 4.
The pathways of spliced genes in the S_M, H_M, H_S. The Y-axis on the left represents KEGG pathways, the Y-axis on the right represents the major category to which each pathway belongs, the X-axis indicates DEmRNAs numbers of each pathway.
Figure 4.
The pathways of spliced genes in the S_M, H_M, H_S. The Y-axis on the left represents KEGG pathways, the Y-axis on the right represents the major category to which each pathway belongs, the X-axis indicates DEmRNAs numbers of each pathway.
Figure 5.
Screening DElncRNAs compared between M group, S group, and H group. A. Distribution of DElncRNAs in each sample, with the Y-axis representing the number of genes and the X-axis representing different samples; B. Venn analysis of novel DElncRNAs obtained from four software programs: CNKI, COC, Pfam, and CPAT. C. Cluster analysis of DElncRNAs in mammary epithelial cells between the control group (M1, M2, M3, M4, and M5), normally treated groups (S1, S2, S3, S4, and S5) and deeply treated groups (H1, H2, and H3). Red indicates highly expressed genes, and green indicates low expressed genes. Each column represents a sample. D E F. Volcano plot of global DElncRNAs in S_M, H_M and H_S, respectively. Gradient red dots represent significantly regulated genes (p<0.05, |log2(fold-change)| >1); dark green dots represent significantly differential expressed genes.
Figure 5.
Screening DElncRNAs compared between M group, S group, and H group. A. Distribution of DElncRNAs in each sample, with the Y-axis representing the number of genes and the X-axis representing different samples; B. Venn analysis of novel DElncRNAs obtained from four software programs: CNKI, COC, Pfam, and CPAT. C. Cluster analysis of DElncRNAs in mammary epithelial cells between the control group (M1, M2, M3, M4, and M5), normally treated groups (S1, S2, S3, S4, and S5) and deeply treated groups (H1, H2, and H3). Red indicates highly expressed genes, and green indicates low expressed genes. Each column represents a sample. D E F. Volcano plot of global DElncRNAs in S_M, H_M and H_S, respectively. Gradient red dots represent significantly regulated genes (p<0.05, |log2(fold-change)| >1); dark green dots represent significantly differential expressed genes.
Figure 6.
Analysis of DElncRNA target genes. The red circle represents lncRNA, and the purple circle represents mRNA.
Figure 6.
Analysis of DElncRNA target genes. The red circle represents lncRNA, and the purple circle represents mRNA.
Figure 7.
Venn map analysis of DElncRNA target genes, DEmRNA and AS genes. DEG means differential expression genes; AS means Alternative splicing; S_M, H_S and H_M represent three comparison groups; Target lncRNA means three comparative groups targeting mRNA.
Figure 7.
Venn map analysis of DElncRNA target genes, DEmRNA and AS genes. DEG means differential expression genes; AS means Alternative splicing; S_M, H_S and H_M represent three comparison groups; Target lncRNA means three comparative groups targeting mRNA.
Figure 8.
Trend analysis of S. agalactiae infection in breast epithelial cells. This series of charts uses Mfuzzy to illustrate the dynamic changes in DEmRNAs during pathogen infection. The yellow and green lines represent genes with small differences in expression changes, while the red and purple lines indicate genes with large differences in expression changes.
Figure 8.
Trend analysis of S. agalactiae infection in breast epithelial cells. This series of charts uses Mfuzzy to illustrate the dynamic changes in DEmRNAs during pathogen infection. The yellow and green lines represent genes with small differences in expression changes, while the red and purple lines indicate genes with large differences in expression changes.
Figure 9.
Screening and enrichment analysis of pDEmRNAs of pathogens in S. agalactiae- normally treated groups (n= 5) compared with S. agalactiae- deeply treated groups. A. Cluster analysis of pDEmRNAs in S. agalactiae between normal treated group (S1, S2, S3, S4, and S5) and deeply treated group (H1, H2, and H3). 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. B. Volcano plot of global pDEmRNAs in S. agalactiae between normal treated group and deeply treated 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); grey dots represent insignificantly differential expressed genes. C. KEGG pathway classified annotation of pDEmRNAs in S. agalactiae. The pathway is exhibited on the left axis, and the area of the circle represents the number of genes listed on the right axis. D. Annotation of pDEmRNAs using Gene Ontology (GO) in S. agalactiae. The rich factor of mRNAs for each GO annotation is exhibited in the left axis.
Figure 9.
Screening and enrichment analysis of pDEmRNAs of pathogens in S. agalactiae- normally treated groups (n= 5) compared with S. agalactiae- deeply treated groups. A. Cluster analysis of pDEmRNAs in S. agalactiae between normal treated group (S1, S2, S3, S4, and S5) and deeply treated group (H1, H2, and H3). 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. B. Volcano plot of global pDEmRNAs in S. agalactiae between normal treated group and deeply treated 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); grey dots represent insignificantly differential expressed genes. C. KEGG pathway classified annotation of pDEmRNAs in S. agalactiae. The pathway is exhibited on the left axis, and the area of the circle represents the number of genes listed on the right axis. D. Annotation of pDEmRNAs using Gene Ontology (GO) in S. agalactiae. The rich factor of mRNAs for each GO annotation is exhibited in the left axis.
Figure 10.
Co-expression network of host cell DElncRNA, DEmRNA, and pathogen pDEmRNA. Red indicates upregulation, blue indicates downregulation, and the color intensity represents strength.
Figure 10.
Co-expression network of host cell DElncRNA, DEmRNA, and pathogen pDEmRNA. Red indicates upregulation, blue indicates downregulation, and the color intensity represents strength.
Figure 11.
The expression level of candidate gene.
Figure 11.
The expression level of candidate gene.
Table 1.
mRNA sequence quality.
Table 1.
mRNA sequence quality.
Sample |
Group |
Total raw reads |
Total clean reads |
Total clean base (G) |
Effective Rate (%) |
Reads with UIDs |
Dedup Reads |
M1 |
Control (M Group) |
81298832 |
70041960 |
10.38 |
86.15 |
64794252(92.51%) |
61112168(87.25%) |
M2 |
80917920 |
70388474 |
10.46 |
86.99 |
65098992(92.49%) |
60322418(85.70%) |
M3 |
92091998 |
79660982 |
11.81 |
86.5 |
73798986(92.64%) |
69202752(86.87%) |
M4 |
82288064 |
70524828 |
10.46 |
85.7 |
65226928(92.49%) |
61014132(86.51%) |
M5 |
91369374 |
79413762 |
11.81 |
86.92 |
73703106(92.81%) |
67255788(84.69%) |
S1 |
Treat1 (S Group) |
102946650 |
92613194 |
13.63 |
89.96 |
86948336(93.88%) |
79863190(86.23%) |
S2 |
71198348 |
63027196 |
9.25 |
88.52 |
59120464(93.80%) |
56144004(89.08%) |
S3 |
86815548 |
76934788 |
11.25 |
88.62 |
72256950(93.92%) |
68365948(88.86%) |
S4 |
92741594 |
82904518 |
12.12 |
89.39 |
77832968(93.88%) |
73052626(88.12%) |
S5 |
104732530 |
93522804 |
13.67 |
89.3 |
87853102(93.94%) |
82741910(88.47%) |
H1 |
Treat2 (H Group) |
60552640 |
43447378 |
6.3 |
72.94 |
41488876(95.49%) |
40567390(93.37%) |
H2 |
72552826 |
51713480 |
7.61 |
71.28 |
49400638(95.53%) |
46882558(90.66%) |
H3 |
83335490 |
62442120 |
9.2 |
74.93 |
59544572(95.36%) |
55991706(89.67%) |
Table 2.
Type of alternative splicing and statistics of differential alternative splicing events.
Table 2.
Type of alternative splicing and statistics of differential alternative splicing events.
EventType. |
NumEvents. JC. only |
SigEvents. JC. only (up: down) |
NumEvents. JC+ Reads On Target |
SigEvents. JC+ reads On Target (up: down) |
S_ M |
H_ M |
H_ S |
S_ M |
H_M |
H_ S |
S_ M |
H_ M |
H_ S |
S_ M |
H_ M |
H_ S |
SE |
36032 |
33110 |
35535 |
428:497 |
557:1034 |
434:818 |
36038 |
33113 |
35536 |
450:533 |
593:1091 |
462:851 |
MXE |
7816 |
6661 |
7424 |
993:1090 |
1398:1132 |
1159:794 |
7816 |
6661 |
7424 |
982:1083 |
1375:1130 |
1144:787 |
A5SS |
348 |
324 |
309 |
17:16 |
19:21 |
9:12 |
349 |
324 |
309 |
17:15 |
21:23 |
13:13 |
A3SS |
418 |
405 |
393 |
12:12 |
18:13 |
12:10 |
418 |
405 |
393 |
12:13 |
19:13 |
12:09 |
RI |
494 |
455 |
426 |
7:21 |
13:37 |
8:17 |
503 |
457 |
432 |
6:18 |
13:32 |
7:13 |