A total of 285 unique pig-encoded miRNAs, including 260 known and 25 novel mature miRNAs, were found to be shared amongst the six libraries. Of these mature miRNAs, 252 miRNAs (235 known and 17 novel miRNAs) were successfully divided into 22 subclusters by K-means clustering (named by K1~K22), 29 subclusters by SOM clustering (named by SOM1~SOM29) and 6 subclusters by Hierarchical clustering (named by H1~H6), each characterized by a unique expression pattern. The patterns corresponding to individual clusters are visualized as line plots (Fig. 1–3). A total of 106 known and 3 novel miRNAs were significantly differentially expressed between the infected and control samples at the three time points (q < 0.001). A total of 4908 predicted target genes of 224 miRNAs were in agreement with the overlap of two prediction methods (miRanda and RNAhybrid), and these constituted our conservative set. However, twenty-eight miRNAs have no target genes in agreement with the overlap of these two methods.
K-means clustering
For K-means clustering, the number of miRNAs and their expression trends at different infection times in different subclusters were shown in Fig. 1. Subclusters K1~K22 had the different number of subjects, and subclusters K9 and K11 had only one subject, with miR-216 and miR-105-1, respectively, while subcluster K20 had the greatest number of subjects, with 27 miRNAs, followed by subcluster K2 that had 26 miRNAs. Additionally, 109 identified differential expression miRNAs (DEMs) were found to be distributed in 15 subclusters, namely K1~5, 7, 12~15, 17, 18, and 20~22. Seven subclusters (K6, K8~11, K16 and K19) did not exhibit the differential miRNA expression patterns specific to any infection stage.
In GO enrichment, 513, 58 and 63 significantly enriched GO terms (q < 0.05) were respectively identified as biological processes (BP), cellular components (CC) and molecular functions (MF) for the target genes of the miRNAs in 9 K-subclusters, including K1~3, 12~15, 17 and 18. Eight BP terms, including regulation of cellular process (GO:0050794), cell communication (GO:0007154), biological regulation (GO:0065007), metabolic process (GO:0008152), signaling (GO:0023052), regulation of biological process (GO:0050789), response to stimulus (GO:0050896) and cellular process (GO:0009987), were shared amongst these nine K-subclusters. Subcluster K13, containing 19 miRNAs, was found to enrich into the least number of GO terms, with 32 terms, including 27 BP, 2 CC and 3 MF (Table 1). Subcluster K15, containing 20 miRNAs, was found to enrich into the greatest number of GO terms, with 443 terms, including 361 BP, 38 CC and 44 MF (Table 1). The enriched GO terms were related to signal transduction, response to stimulus/stress, immune process, apoptotic process, metabolic process and biosynthetic process for BP, and catalytic activity, GTPase activity and nucleotide binding for MF (Additional file 1: Table S1). And the GO terms related to defense response, innate immune response, NF-κB transcription, tumor necrosis factor-mediated signaling pathway and inflammatory response involving to immune response-regulating cell surface receptor signaling pathway, cytokine secretion and leukocyte activation, aggregation and migration were only enriched by the K15 subcluster. Additionally, GO terms related to angiogenesis, interleukin-1-mediated signaling pathway, Notch signaling pathway and vasculature development were found to be subcluster K1-sepcific (Additional file 1: Table S1). And amino acid transport-related terms were the most frequently represented objects in subcluster K2-sepcific GO terms, and cell cycle process and DNA replication-related terms for K12, protein localization and transport for K18 (Additional file 1: Table S1).
In the KEGG enrichment, 15 KEGG pathways involving to Cellular Processes (Peroxisome, Regulation of actin cytoskeleton, Focal adhesion, Phagosome), Human Diseases [Dilated cardiomyopathy (DCM), Toxoplasmosis, Hepatitis B, Chagas disease (American trypanosomiasis), Pertussis, Proteoglycans in cancer], Signal transduction (MAPK signaling pathway, TNF signaling pathway), Immune system (Th1 and Th2 cell differentiation, Hematopoietic cell lineage), Metabolism (Fatty acid degradation), were significantly enriched by 6 subclusters (q < 0.05), including K1, 2, 12, 14, 15 and 17 (Table 1, Additional file 2: Table S2). Of these, Peroxisome (ssc04146) was shared amongst K-subclusters K1, K15 and K17. Subcluster K2, containing 26 miRNAs, was found to enrich into the greatest number (6) of KEGG pathways. Of these, Th1 and Th2 cell differentiation (ssc04658) and Toxoplasmosis (ssc05145) was subcluster K2-specific (Additional file 2: Table S2).
In the Reactome pathways enrichment, 16 Reactome pathways involving to Metabolism, Immune system, Cell-Cell communication and Signal transduction were significantly enriched (q < 0.05) by subclusters K1, K2, K3, K12, K15 and K17 (Table 1, Additional file 3: Table S3). Interestingly, all of Metabolism-related pathways, including biosynthesis of phosphatidylcholine and Glycerophospholipid (R-SSC-1483191, R-SSC-1483206), and Metabolism of Phospholipid and Cobalamin (R-SSC-1483257, R-SSC-196741), were found only in sublcuster K1. All of Cell connection-related pathways, including Cell junction organization (R-SSC-446728), Cell-Cell communication (R-SSC-1500931) and Cell-extracellular matrix interactions (R-SSC-446353), were found only in sublcuster K2. The Immune System-related pathways were found to be the most frequently represented pathways, which were found to distribute in subclusters K2, K3, K12 and K15 (Additional file 3: Table S3). Additionally, no significantly enriched GO, KEGG and Reactome terms were found in subclusters K4~11, 16 and 19~22.
SOM clustering
For SOM clustering, the number of miRNAs and their expression trends at different infection times in different subclusters were shown in Fig. 2. Subclusters SOM1~SOM29 had the different number of subjects, and subcluster SOM26 had only one subject, with miR-326, while subcluster SOM1 had the greatest number of subjects, with 32 miRNAs, followed by subcluster SOM29 that had 19 miRNAs. Additionally, 109 identified DEMs were found to be distributed in 26 subclusters, namely SOM1~5, 7~14, 16~25, 27~29. Three subclusters (SOM6, SOM15 and SOM26) did not exhibit the differential miRNA expression patterns specific to any infection stage.
In GO enrichment, a total of 660 GO terms, including 552 BP, 46 CC and 62 MF, were significantly enriched (q < 0.05) by the target genes of the miRNAs from 11 SOM-subclusters, including SOM1, 8~10, 12~14, 17, 18, 25 and 29, respectively (Table 2). Response to stimulus (GO:0050896) was shared amongst 9 SOM-subclusters (SOM1, 8~10, 13, 14, 17, 25 and 29), and was found to be the most frequently represented terms. Subcluster SOM25, containing 17 miRNAs, were found to be the most number of GO terms, with 457 terms, including 386 BP, 29 CC and 42 MF (Table 2), of which 110 terms were only found in this subcluster and were mostly associated with calcium ion transport, stress-activated responses, MAPK cascade signaling, and cell apoptosis (Additional file 4: Table S4). Interestingly, the specific GO terms were found to be enriched by some subclusters. Ion binding and ion transport-related terms were only found in subclusters SOM8 and SOM12, respectively. And all of the four terms for subcluster SOM9-specific were related to leukocyte migration and chemotaxis, including lymphocyte migration (GO:0072676), leukocyte migration (GO:0050900), chemokine activity (GO:0008009) and chemokine receptor binding (GO:0042379). Furthermore, GO terms related to cytokine production and secretion, interleukin-1 production, GTPase activity, and NF-κB signaling were found to be the most frequently represented objects in the subcluster SOM17-specific terms (Additional file 4: Table S4).
In the KEGG enrichment, 26 KEGG pathways involving to Cellular Processes, Signal transduction, Human Diseases, and Organismal Systems were significantly enriched by 8 subclusters (q < 0.05), including SOM1, SOM8, SOM12~14, SOM17, SOM25 and SOM29 (Table 2, Additional file 5: Table S5). Of these, 19 KEGG pathways were classied into Human Diseases (11) and Signal transduction (8). MAPK signaling pathway (ssc04010) and TNF signaling pathway (ssc04668) were showed to be the most frequently pathways enriched amongst SOM subclusters (Additional file 5: Table S5). MAPK signaling pathway was shared amongst subclusters SOM8, SOM25 and SOM29, and TNF signaling pathway was common to subclusters SOM14, SOM25 and SOM29. Noteworthy, subcluster SOM29, containing 19 miRNAs, was found to enrich into the greatest number (10) of KEGG pathways, followed by SOM25 with enrichment into 9 KEGG pathways. Of these, 8 and 5 KEGG pathways were found to be subcluster SOM29- and SOM25-specific, respectively (Additional file 5: Table S5).
In the Reactome pathways enrichment, 15 Reactome pathways involving to Cell-Cell communication, Gene expression, Immune system, Signal transduction, Developmental Biology, DNA Repair, and Hemostasis were significantly enriched by subclusters SOM8, SOM9, SOM17, SOM25 and SOM29 (q < 0.05). SOM25 was found to enrich into the greatest number (9) of Reactome pathways (Table 2). Notably, all of the Gene expression-related pathways including RNA Polymerase III Transcription Initiation From Type 1 Promoter (R-SSC-76061), RNA Polymerase III Transcription Initiation From Type 2 Promoter (R-SSC-76066), RNA Polymerase III Transcription (R-SSC-74158) and RNA Polymerase III Transcription Initiation (R-SSC-76046), and DNA repair-related pathways including Gap-filling DNA repair synthesis and ligation in TC-NER (R-SSC-6782210) and Dual incision in TC-NER (R-SSC-6782135), were enriched only by subcluster SOM 25. Cytokine Signaling in Immune system (R-SSC-1280215) was the most frequently represented pathways, which were found to distribute in subclusters SOM17, SOM25 and SOM29 (Additional file 6: Table S6). Additionally, no significantly enriched GO, KEGG and Reactome terms were found in subcluster SOM2~7, 11, 15, 16, 19~24, 26 and 27.
Hierarchical clustering
For Hierarchical clustering, the number of miRNAs and their expression trends at different infection times in different subclusters were shown in Fig. 3. Subclusters H1~H6 had the different number of subjects, and subcluster H5 had the greatest number of subjects, with 67 miRNAs, while subcluster H6 had the least number of subjects, with 8 miRNAs.
In GO enrichment, a total of 477 GO terms, including 371 BP, 46 CC and 60 MF, were significantly enriched (q < 0.05) by the target genes of the miRNAs from 5 H-subclusters, including H1~H5, respectively (Table 3). Of these, 41 GO terms were shared amongst 5 H-subclusters (H1~H5), which the most were associated with cell communication, cellular process and response to stimulus. Subcluster H5, containing 67 miRNAs, were found to be the most number of GO terms, with 376 terms, including 304 BP, 32 CC and 40 MF (Table 3), of which 94 terms were only found in this subcluster and were mostly associated with apoptotic process, DNA-binding transcription factor activity, and ion homeostasis and transport (Additional file 7: Table S7). Additionally, further analysis revealed that some GO terms were subcluster-specific. For example, RNA splicing-related terms were the most frequently represented objects in the subcluster H1-specific terms, and lymphocyte activation-related terms were common in the subcluster H3-specific terms (Additional file 7: Table S7).
In the KEGG enrichment, 14 KEGG pathways involving to Signal transduction, Human Diseases and Immune system were significantly enriched by 4 subclusters (q < 0.05), including H1~H3 and H5 (Table 3, Additional file 8: Table S8). Of these, 7 KEGG pathways were classified into Human Diseases involving to Cancer (3), Cardiovascular disease (1), Immune disease (1) and Infectious disease (2) (Additional file 8: Table S8). Subcluster H3, containing 39 miRNAs, were found to be the most number of KEGG pathways (6), followed by H2 with 4 KEGG pathways (Table 3).
In the Reactome pathways enrichment, 7 Reactome pathways involving to immune system and autophagy were significantly enriched (q < 0.05) by subclusters H1~H5 (Table 3, Additional file 9: Table S9). Of these, Immune system (R-SSC-168256) pathway was shared amongst subclusters H2~H5. H5 was found to enrich into the greatest number (4) of Reactome pathways. It is worth mentioning that Macroautophagy (R-SSC-1632852) was enriched only by sublcuster H1. And Innate immune system-related pathways including Innate Immune System (R-SSC-168249), Antimicrobial peptides (R-SSC-6803157) and Toll Like Receptor 4 (TLR4) Cascade (R-SSC-166016), were enriched by the subclusters H2 and H5. Furthermore, two Cytokine Signaling in Immune system-related pathways, including Cytokine Signaling in Immune system (R-SSC-1280215) and Signaling by Interleukins (R-SSC-449147), were found only in subcluster H5 (Additional file 9: Table S9). Additionally, no significantly enriched GO, KEGG and Reactome terms were found in subcluster H6.
The overlap of miRNAs among K-means clustering, SOM clustering and Hierarchical clustering
Considerable bias was found among K-means, SOM and Hierarchical clustering in regards to the representation of clustered genes for the three function enrichments, which mainly due to discrepancy in calculation methods for miRNA expression patterns. The clusters obtained by three algorithms were compared to find out similar expression patterns, up to 22 miRNAs mainly showing the downregulated expression at 50 DPI were identified into one subcluster (namely subcluster H3-K17-SOM1) through the three algorithms in this study, and were included in the subcluster K17 for K-means clustering, subcluster SOM1 for SOM clustering and subcluster H3 for Hierarchical clustering, respectively. In addition to novel-201, the remaining 21 miRNAs were classified into 12 miRNA families (let-7, mir-1, mir-10, mir-128, mir-129, mir-139, mir-1468, mir-27, mir-28, mir-30, mir-328, mir-455), of which 6 miRNA members, including let-7a, let-7c, let-7e, let-7f, let-7g and let-7i, belonged to the let-7 family, and 4 miRNAs, including miR-30a-3p, miR-30b-3p, miR-30c-3p and miR-30e-3p, belonged to the mir-30 family. A total of 860 target genes corresponding to 20 miRNAs (no containing novel-201 and miR-99a) in subcluster H3-K17-SOM1 were further used to analyze their cellular functions.
GO enrichment analysis found that a total of 400 GO terms were significantly enriched (p < 0.05) by the subcluster H3-K17-SOM1, and the main enriched GO terms is shown in Fig. 4. And the target genes of miRNAs in this subcluster were involved in apoptotic process, cell adhesion, GTPase activity, protein deubiquitination, and especially immune response including B cell differentiation, immunoglobulin production, inflammatory response, innate immune response, interleukin production, natural killer cell mediated cytotoxicity, neutrophil chemotaxis, T cell cytokine production, T cell differentiation and T cell proliferation (Additional file 10: Table S10). In addition, enrichment results showed that the most of GO terms (14 terms) of cellular components were related to mitochondrion, suggesting many target genes of miRNAs in subcluster H3-K17-SOM1 located in mitochondrion may play important roles in mitochondrial function during T. gondii infection.
KEGG enrichment analysis revealed that a total of 121 pathways were potentially regulated by the miRNAs in subcluster H3-K17-SOM1 during T. gondii infection, and the main enriched KEGG pathways is shown in Fig. 5. The 28 KEGG pathways involving to human diseases, signal transduction, immune system, cellular processes and metabolism were significantly enriched (q < 0.05) by the subcluster (Additional file 11: Table S11). Of these, three immune system-related KEGG pathways, including Chemokine signaling pathway (ssc04062), B cell receptor (BCR) signaling pathway (ssc04662), and C-type lectin receptor (CLR) signaling pathway (ssc04625), were not found to be significantly enriched by any subcluster of K-means clustering, SOM clustering or Hierarchical clustering.
Reactome enrichment analysis showed that 12 pathways were significantly enriched (q < 0.05) by the subcluster H3-K17-SOM1, of which five pathways were classed into immune system, especially cytokine Signaling in Immune system, following by three metabolism- and two cell cycle-related pathways (Additional file 11: Table S11).
Regulatory network was constructed for the immune-related target genes and miRNAs in subcluster H3-K17-SOM1 (Fig. 6). A total of 174 target genes encoding the different cytokines involving in immune regulation of the host, including chemokines, interleukins (ILs), tumor necrosis factors (TNF), interferons (IFN), Toll-like receptors, C-type lectin, Immunoglobulin, etc., were potentially regulated by 16 miRNAs. Among these miRNAs, ssc-miR-328 regulated the greatest number of target genes, with 70 targets, following by ssc-miR-30b-3p regulating 31 target genes. The target genes regulated by the greatest number of miRNAs were PTPN6 and TMEM173, each regulated by three miRNAs (Fig. 6). These results indicate that these DEMs and cytokines might play central roles in the host response to T. gondii infection.