In order to build, visualize, and analyze the DEG interaction networks in the context of innate immunity, we mapped our data to the InnateDB Network Analysis. Table 1 shows an example of InnateDB-curated molecular interactions generated from our data. In total, 2,458 InnateDB-curated interactions were generated from the 5,409 DEG at the acute stage of infection (FDR < 0.05). A smaller set comprising 224 InnateDB-curated interactions were generated from the 2,364 DEG at the chronic stage of infection (FDR < 0.05). The full molecular interactions from both DEG datasets are listed in Tables S1 and S2 in Supplementary Material, respectively.

Although curated interactions in InnateDB are annotated in a wide range of cell and tissue types, the majority of these interactions are from studies involving cell lines rather than primary cells (21). This is reflected in our results with 1,622 interactions from the acute stage and 159 interactions from the chronic stage of infection predicted from at least one source based on cell lines. Since the interactome is very much dependent on the context of a particular cell type under investigation (21), we focused on the interactions from studies involving PBMC-related cell types including B cells, T cells, NK cells, and monocytes. From samples representing the acute stage of infection, 458 PBMC-related interactions were identified from 2,458 InnateDB-curated interactions and 336 distinct DEG involved in these interactions (Table S1 in Supplementary Material). For the chronic stage of infection, 47 PBMC-related interactions were identified from 224 InnateDB-curated interactions, in which 55 distinct DEG were involved (Table S2 in Supplementary Material).

Among the 2,458 InnateDB-curated interactions representing the acute stage of infection, 92 interactions were identified, which were associated with evidence from fibrosis-associated cell types including fibroblast and stellate cell lines (the major cells involved in hepatic fibrosis) (Table S1 in Supplementary Material).

In order to gain a systems-level understanding of interactions of immunological interest, we used NetworkAnalyst to construct and visualize the molecular interaction networks (18). We first input the 336 DEG from the acute stage of infection that are involved in the PBMC-associated interactome with the corresponding fold-change values and generated a default network with 7,549 nodes and 24,354 edges (first-order interactors, including all identified proteins that directly interact with the input proteins, also known as seed nodes). It is notable that, apart from the seed nodes, there were a large number of directly interacting proteins identified, which made the default network too large and would lead to a “hairball” effect that rarely produces any informative outcome. For practical (visual, biological, and computational) reasons, a “Reduce” function was used to prune the network down to a more manageable size (200–2000 nodes) as recommended by Xia et al. (18), and a refined network with 331 nodes and 1,769 edges was retained for further analysis (zero-order interactions that refer to connections between two input seed nodes). The “Reduce” function retained only the seed nodes and the direct connections (edges) among them. Two widely used topological measures—“degree of connectivity” and “between centrality” were used to help identification of importance of nodes within the network (Figure S1 in Supplementary Material) (18). Degree of connectivity refers to the number of connections the node has to other nodes. Betweenness centrality refers to the number of shortest paths going through the node (22). In the current network (Figure S2 in Supplementary Material), the top five nodes with the highest node degree values are “reticuloendotheliosis viral oncogene homolog A” (RELA; degree value = 66), “tumor protein p53” (TP53; 65), “transcription factor Sp1” (SP1; 63), “amyloid precursor protein” (APP; 56), and “tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein” ζ (YWHAZ; 51). The top five nodes with the highest betweenness values are TP53 (betweenness value = 2,264.97), RELA (2,164.28), SP1 (1,505.93), “C-terminal Src kinase” (SRC; 1,271.48), and “breast cancer 1” (BRCA1; 1,269.96).

In the chronic stage of infection, 47 interactions are confirmed to stem from studies involving PBMC-related cell types, in which only 55 distinct DEG were involved. This figure was much lower than the recommended number of seed nodes (200–2,000) for network construction (18). Hence, we also performed network construction on the 186 distinct DEG involved in the 224 InnateDB-curated interactions using the same method for the dataset of the acute stage of infection. A preliminary network with 5,035 nodes and 11,422 edges was generated. After processing by the “Reduce” function built-in to NetworkAnalyst, a refined network with 158 nodes and 305 edges was retained for further analysis (only zero-order interactions among the input seed nodes were kept). Figure S3 in Supplementary Material shows the topology of the network generated. The top five nodes ranked by node degree are COPS5 (COP9 constitutive photomorphogenic homolog subunit 5), IKBKB (inhibitor of nuclear factor kappa-B kinase subunit β), YWHAE (14-3-3 protein epsilon, tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein epsilon), BRCA1, and JUN (Jun proto-oncogene, AP-1 transcription factor subunit) with their node degree values 21, 18, 17, 16, and 15 respectively. The top five nodes ranked by betweenness values are BRCA1, COPS5, LYL1 (lymphoblastic leukemia-associated hematopoiesis regulator 1), S-phase kinase-associated protein 2 (SKP2), and CSNK2B (Casein kinase 2, β polypeptide) with their betweenness values 271.25, 148.23, 139.17, 120.4, and 103.67, respectively.

We first performed function-first approach analysis on the networks representing acute and chronic stages of infection using the Reactome database.

In total, 121 and 63 Reactome functional pathways were enriched (P < 0.001) in networks representing acute and chronic stages of infection, respectively. Tables 2 and 3 show the list of the top 25 pathways enriched in both stages of infection, respectively.

In the network representing the acute stage of infection, five pathways associated with TLR3/TLR4 signaling were ranked in the top 10 overlapping pathways identified, including “TLR4 cascade,” “activated TLR4 signaling,” “TRIF-mediated TLR3/TLR4 signaling,” “MyD88-independent cascade,” and “TLR3 cascade” (Table 2). Since TLR3 and TLR4 share a specific TRIF (TIR-domain-containing adapter-inducing interferon-β)-dependent pathway (unlike other TLRs) (23, 24), we extracted a subnetwork that contained all nodes involved in these five pathways (Table S3 in Supplementary Material). Subsequent module exploration analysis generated four significantly overlapping modules within this subnetwork (Figure 1). Then, functional enrichment analysis was repeated on these modules. Among them, Module 1 contains six nodes [NOD2, mitogen-activated protein kinase (MAPK) 14, MAPK1, MAPK3, interleukin-1 receptor-associated kinase (IRAK1), and IRAK4], which are involved in a range of TLR cascades, including TLR3/4/7/8/10, indicating a cross talk between TLR pathways. Module 3 contains five nodes mainly involved in My88-independent TRIF-dependent TLR signaling, including interferon regulatory factor 7 (IRF7), interferon regulatory factor 3 (IRF3), conserved helix-loop-helix ubiquitous kinase (CHUK), mitogen-activated protein kinase kinase kinase 7 (MAP3K7), and Fas-associated protein with death domain (FADD).

Similar to the findings for the acute stage of infection, the five pathways related to TLR4/TLR3 also presented in the top 10 enriched pathways in the chronic stage of infection, including “TLR4 cascade,” “activated TLR4 signaling,” “TRIF-mediated TLR3/TLR4 signaling,” “MyD88-independent cascade,” and “TLR3 cascade.” We extracted a subnetwork that contained all nodes that are involved in these five pathways (Table S3 in Supplementary Material; Figure 2). There were no significant overlapping modules found within this subnetwork, which may be due to the small size of the subnetwork.

Almost all the top ranked pathways enriched by the function-first approach in both networks are associated with immune responses. It is, thus, clear that the pre-enrichment of DEG from the immune-related interactome (which have been confirmed in the context of PBMC-related cell types) is a reliable and effective way to extract immune-related networks and refine the network size for further analysis. In order to obtain tightly clustered modules, we conducted a complementary approach (connection-first) to re-analyze the two networks separately.

Ninety-two InnateDB-curated interactions from the acute stage of infection were identified, which included evidence from fibrosis-associated cell types including fibroblast and stellate cell lines. Some genes related to the TGFβ-mediated fibrogenic progresses (such as SMAD3, SMAD4, COL1A1, and TGFB1) were involved in the 92 interactions. In our previous IPA analysis on this dataset, we hypothesized that the increased expression of these genes may be possibly associated with the fibrosis in the ovine liver during F. hepatica infection (15). A previous study using hepatic lymph nodes has also proposed that TGFβ1 plays a central role in fibrosis during ovine F. hepatica infection (25). It is well known generally that persistent TGFβ signaling leads to excessive fibrosis and ultimately scarring of internal organs (26). TGFβ induces the transcription of the gene COL1A1, encoding collagen type I, which is the major fibrous collagen and plays a central role in wound-healing (26). Another study demonstrated that the transcript expression levels of TGFB1 and COL1A1 were significantly increased in ovine liver tissue at eight wpi of F. hepatica (27). Notably, the cells being analyzed in this study are PBMC, and it is, therefore, notable that 92 interactions also have been demonstrated in fibrosis-associated cell types that may directly interact with F. hepatica during the local infection in the liver. These data may indicate that in response to F. hepatica infection, blood mononuclear cells could migrate into the infected liver parenchyma and release TGFB1 and COLIA1 to activate fibrosis-associated cell types (e.g., stellate cells) and facilitate hepatic fibrosis. Or perhaps the ovine PBMC and fibrosis-associated cell types might perform same gene expression changes associated with TGFβ—collagen pathway in the context of F. hepatica infection. Further investigations into links between the systemic immune system and local infection in the liver will be valuable.

Node importance in network analysis can be estimated as degree of connectivity (the number of connections the node has to other nodes) and betweenness centrality (the number of shortest paths going through the node) (22). Nodes with a high degree of centrality are referred to as “hub nodes,” whereas nodes with a high betweenness centrality are referred as “bottleneck nodes” (28). Notably, in our analysis, some molecules performed as both hub and bottleneck nodes in the network, indicating their crucial roles in regulating the innate immune response to F. hepatica infection. For instance, RELA and SP1 are in the top of both hub and bottleneck nodes in the network representing the acute stage of infection. RELA encodes the transcription factor NF-kBp65, one of the active subunits of the NF-Kb complex involved in NF-Kb binding to DNA. Previous studies have suggested that NF-kBp65 participates in adaptive immunity and responses to invading pathogens via NF-kB activation (29). One study has demonstrated that F. hepatica tegument can target NF-kBp65 and suppress its protein expression in DC, which subsequently resulted in decreased production of pro-inflammatory cytokines in DC (14). It has also been shown that the recombinant form of a major F. hepatica-secreted molecule, sigma class glutathione transferase (rFhGST-si) could induce partial maturation of DC in vitro by enhancing the protein production of NF-kBp65 and signaling through it (12). These studies suggest that different components of F. hepatica can affect DC maturation and function through regulating RELA expression in different ways. SP1 is a zinc finger transcription factor that is a central mediator of IL-10 induction (30, 31). Previous studies proposed that a helminth immune modulator (cystatin from Acanthocheilonema viteae) can bind to macrophages and activate SP1 via p38 signaling pathway to initiate transcript expression of IL10 (32). Since IL10 is an important regulatory cytokine induced by F. hepatica infection (33), the high importance of SP1 within the network may indicate its potential role in regulating IL10 during the acute stage of F. hepatica infection.

Breast cancer 1 is one of the top five bottleneck nodes in both networks at the acute and chronic stages of infection. Much current research on BRCA1 focuses on its role in DNA damage repair (34), but its role in the context of immunology/infection is rarely characterized. One study suggested that BRCA1 can activate IFN-γ signaling to induce target gene IRF-7 expression and to subsequently stimulate the innate immune response and elicit apoptosis (35). In our DEG data, BRCA1 showed a significantly decreased transcript expression at the acute stage of F. hepatica infection (Table S1 in Supplementary Material), which dropped down to an even lower level at the chronic stage of infection (Table S2 in Supplementary Material). This provides a potential mechanism for the decreased transcription of IFNG in our DEG dataset (Table S1 in Supplementary Material) and for the downregulation of Th1 immune responsiveness seen in F. hepatica infection. The high betweenness value of BRCA1 in both networks indicates that F. hepatica infection may block IFN-γ signaling and further attenuate Th1 immune responses through suppressing the expression of BRCA1.

The initiation of TLR3/TLR4 signaling by F. hepatica infection and its downstream effects are not fully understood. The ligand of TLR3 is double-stranded RNA (dsRNA), which has been well-studied in the field of viral infection. One study has demonstrated that dsRNA from Schistosoma mansoni eggs is a natural ligand for TLR3, and this egg-derived dsRNA was capable of activating DC through a fully TLR3-dependent pathway (36). Hence, regulation of TLR3-dependent cellular pathways could be used by helminth parasites to modulate host innate immune responses. TLR3 suppression has been demonstrated to play a role in promoting Th2 immune responses to S. mansoni infection. In TLR3 gene-deficient mice, the absence of TLR3 resulted in significantly decreased Th1-associated factors (e.g., IL12) and significantly increased Th2 cytokines and chemokines in response to S. mansoni infection (37). In addition, macrophages from TLR3 deficient mice were alternatively activated (37). These data suggest that the absence of TLR3 leads to a strongly skewed Th2 response during S. mansoni infection. During the early stage of infection, the penetration and migration of F. hepatica can cause severe injury of liver tissue. To avoid massive immunopathology, F. hepatica migrating juveniles have to prevent cells of the innate immune system promoting a Th1-adaptive response, thus making them more likely to drive the development of Th2 immune responses (38). Regulation of TLR3 signaling has been considered to be involved in this polarization process. Previous studies have demonstrated that F. hepatica cathepsin L1 (FhCL1, a major component of the excretory-secretory proteins) could suppress the production of proinflammatory mediators from LPS-stimulated macrophages by inhibiting TLR3 signaling (38). In this study, Module 3 extracted from the TLR3/TLR4-associated subnetwork at the acute stage of infection contains five nodes that play crucial roles in My88-independent TRIF-dependent TLR signaling (Figure 1). Among them, IRF3 is an important transcription factor involved in the downstream regulation of TLR3-TRIF (39), which regulates type 1 IFNs (e.g., IFN-β), leading to subsequent induction of iNOS (inducible nitric oxide synthase, encoded by gene NOS2) expression (38, 39). IRF3 showed significantly decreased expression at the acute stage of infection (logFC = −7.07) (Table S1 in Supplementary Material), with similarly low expression levels until the chronic stage of infection (15). The decreased expression of IRF3 and the related Module 3 might bridge the suppression of TRIF-dependent TLR3 signaling to the extremely decreased transcript expression of NOS2 during F. hepatica infection observed in our previous study (15). iNOS catalyzes the conversion of arginine into citrulline and nitric oxide (NO), which is a defense mechanism in response to pathogen invasion (40). The suppression of NOS2 transcript expression in macrophages is an important characteristic of F. hepatica infection (41), which could be one mechanism of immune evasion by the parasite to avoid damage from host-generated nitric oxide. Previous studies have shown that FhCL1 can suppress murine macrophage activation by the degradation of TLR3 protein within endosomes, rather than downregulation of TRIF and TLR3 transcript expression (38). This is consistent with the fact that TRIF and TLR3 are absent in our DEG datasets. It would be useful to investigate the effect of FhCL1 on downstream TLR3 signaling, such as transcript expression of IRF3. Our study also highlighted an interaction between TLR3 and TLR4 TRIF-dependent pathways, which may be affected by FhCL1 (38). However, other as yet unidentified mechanisms may be involved in the blocking of TRIF-dependent TLR3/TLR4 pathways induced by F. hepatica at the transcription regulation level.

Toll-like receptor 4 is well known as the receptor for bacterial endotoxins, several viral proteins, and polysaccharide and plays a pivotal role in the induction of inflammatory responses (42). In comparison to TLR3, there are more studies on the association between TLR4 and F. hepatica infection. One previous study suggested that FhES was capable of inhibiting the activation of bovine macrophages by TLR4 antagonists, indicating F. hepatica may modulate the host innate immune response through either altering the expression of TLR4 or interfering with its subsequent signaling (43). A recent study suggested that TLR4 may be involved in F. hepatica fatty acid-binding protein (FhFABP) induced alternative-activation of human macrophages (44). Another study suggested that FhFABP can inhibit TLR4 activation in macrophages and suppresses the inflammatory cytokines induced following incubation of LPS (45). It has also been suggested that TLR4 plays a crucial role in mediating the activation of a suppressive DC phenotype by two major components of FhES, FhCL1 and FhGST-si (12). However, TLR4 does not always account for the modulation of host DC. One study has demonstrated that FhTeg suppresses DC maturation and function independently of the TLR4 pathway, since it still functioned in DC generated from TLR4-deficient mice (14). Hence, it can be concluded that different components of F. hepatica have diverse influence on TLR4 transcript expression and related signaling, and TLR4 may play various roles in different cell types during F. hepatica infection. In the TLR3/TLR4-associated subnetwork from the chronic stage of infection (Figure 2), JUN has the highest node degree and betweenness values. JUN encodes the transcription factor c-Jun, which can be activated and phosphorylated by c-Jun N-terminal kinases (JNK) (46). Previous studies suggested that JNK was implicated in the TLR4-mediated intracellular signaling pathway through which rFhGST-si (the recombinant form of one major F. hepatica secreted protein) signaled to suppress DC maturation and function (12). Downstream of the JNK signaling pathway, activated c-Jun, and c-Fos function as subunits of a group of dimeric transcription factors were collectively referred to as activating protein (AP-1) that regulate transcript expression of target genes (47). It has been proposed that AP-1 complexes containing c-Fos might repress genes that are implicated in DC maturation, since continued expression of c-Fos was detrimental for DC maturation (48). In our analysis, genes coding JNK (MAPK8), c-Jun (JUN), and c-Fos (FOS) were all implicated in the subnetwork from the chronic stage of infection (Figure 2), and all showed significantly increased transcript expression levels at this stage (Table S2 in Supplementary Material). Hence, we hypothesize that activation of the TLR4-JNK-AP-1 signaling pathway could be a potential mechanism used by F. hepatica to modulate maturation and function of host DC during chronic infection.

By comparing the TLR3/TLR4-related subnetworks generated from acute and chronic stages of infection, we found several nodes (including LY96, TLR1, TLR4, IRAK1, MAPK3, MAP3K7, BRCA1, NFKBIB, COPS5, and JUN) that were consistently involved in both subnetworks (Figure 2), suggesting that these DEG may play a consistent role in regulation of TLR3/TLR4 function during F. hepatica infection. Among them, LY96 encodes a protein named lymphocyte antigen 96 (also known as MD2), which has been shown to play a role in the binding between LPS and TLR4 (49). Our group previously proposed that the ability of FhES to produce M2 macrophages may be partially dependent on parasite glycan residues, and the interaction between FhES and TLR4 is possibly involved in this mechanism (43). Given the consistent expression regulation of LY96 in this study and its role in TLR4 function, we propose that LY96 is critical in TLR4 recognition of fluke-derived glycan residues and may play a role in induction of M2 macrophages in the presence of F. hepatica secretory glycoproteins.

In summary, we provide here the first investigation of gene networks associated with TLR3/TLR4 signaling/cascade within PBMC during F. hepatica infection, which provides useful information for further investigating the roles of TLR3/TLR4 in single cell types, during helminth infection.

Module 6 extracted from the network representing the acute stage of infection significantly enriched four Reactome pathways associated with intrinsic apoptosis and shows a general trend of increased gene expression (Figure 3). This may indicate an activated intrinsic pathway of apoptosis in PBMC at early stage of F. hepatica infection, which is in agreement with our previous observations from the IPA pathway “Death Receptor Signaling” and “Apoptosis Signaling” (15). BH3-only proteins, structurally characterized by the BH3 domain [Bcl-2 homology (BH) domain 3], are a main pro-apoptotic subgroup in the BCL-2 (B-cell CLL/lymphoma 2) family that are crucial regulators of apoptosis during the intrinsic pathway (50). To date, several BH3-only proteins have been discovered in human and mouse (51). In our DEG data, there were three BH3-only member genes showing significantly increased expression during the acute stage of infection, including BIK (BCL-2-interacting killer; also known as NBK; logFC = 5.04), BNIP1 (BCL2/adenovirus E1B 19 kDa interacting protein 1; logFC = 1.69), and BCL2L11 (Bcl-2-like protein 11; also known as BIM; logFC = 1.43). In addition, one gene encoding another BH3-only protein, PMAIP1 (phorbol-12-myristate-13-acetate-induced protein 1; also known as NOXA; logFC = 3.89), also showed an increased expression in chronic stage of infection (15). Previous studies have shown that BIM is likely to be involved in apoptosis of lymphocytes [summarized from Strasser (50)], and also has a crucial role in the termination of T-cell immune responses (52, 53). We hypothesize that these four BH3-only proteins may be associated with previously observed apoptosis of peripheral blood immune cells during F. hepatica infection (54, 55). Since the top three enriched pathways in Module 6 were associated with BH3-only proteins, the increased expression of these BH3-only genes in the DEG dataset could be potentially attributed to the activation of Module 6.

Module 1 extracted from the network representing the chronic stage of infection shows a general trend of activation and enriched three BCR-related pathways (Figure 4). Several nodes involved in all the three pathways were UBA52 (ubiquitin A-52 residue ribosomal protein fusion product 1), growth factor receptor-bound protein 2 (GRB2), IRS2 (insulin receptor substrate 2), MAP3K7 (also known as TAK1), IKBKB, and NFKBIB. BCR plays a critical role in initiation of B-cell responses and guides cell maturation, survival, energy, and the production of antibodies in plasma B-cells (summarized from Dal Porto et al. (56) and Liu et al. (57). As shown in Figure 4, GRB2 is the most important node in this module, which has the highest degree and betweenness values. GRB2 encodes an adapter protein that plays a role in recruiting and assembling the BCR signalosome and the amplification of its signal (58, 59). In another RNA-seq study on ovine PBMC following F. hepatica infection, the authors suggested that F. hepatica infection might induce activation of BCR through downregulating CD22, a negative modulator of BCR (60). Although CD22 was not detected in our DEG dataset, our data do suggest that the increased transcript expression of GRB2 and the associated module might be involved in activation of downstream signaling of BCR during F. hepatica infection. Another highlighted DEG of interest is MAP3K7, encoding an enzyme that is critical for B-cell maturation and BCR-mediated B-cell proliferation (61, 62). It has previously been demonstrated that mice with B-cell-specific MAP3K7 deficiency have a significantly reduced basal titers of IgG1, IgG2a/b, IgG3, and IgA in serum compared with those of control mice (61). In addition, the production of antigen-specific IgM and IgG1 in response to a T cell-dependent antigen (nitrophenol-chicken g-globulin) was considerably impaired in MAP3K7-deficient mice compared with that of controls. This study indicated that MAP3K7 is indispensable for in vivo induction of humoral immune responses as well as isotype switching (61). Notably, transcript expression levels of MAP3K7 in our study consistently decreased from the acute through to the chronic stage of F. hepatica infection (Tables S1 and S2 in Supplementary Material). It is well known that F. hepatica influences IgG1/IgG2 class-switching and antibody affinity maturation (63–65). Hence, we propose that the consistent suppression of MAP3K7 and its role in this module could be implicated in the F. hepatica-induced modulation of BCR signaling and humoral immune responses.

This study for the first time investigates the global gene expression changes associated with innate immunity to F. hepatica infection at a systems level. Our findings revealed novel insights into F. hepatica-mediated modulation of host immune responses, which facilitate parasite survival and chronic infection. The data from this study will benefit the rational design of vaccines which to be effective, will need to overcome the extreme parasite-mediated immunoregulation seen in fasciolosis. As well as their direct applicability in the design of vaccines for the control of this major disease of livestock, these insights also contribute to the general understanding of helminth-mediated immunoregulation [reviewed in Ref. (66)]. A substantial body of evidence now points to the multiple layers of immunoregulation and immunodiversion employed by helminths to limit immunopathology and survive long-term in their hosts. However, a complete mechanistic understanding of how such profound effects of the host’s immune system are orchestrated is still developing. Furthering this understanding is important not only for advancing helminth vaccinology, but also for harnessing similar immunoregulatory strategies for the treatment of autoimmune and immunopathological states. Our data were obtained using a complex PBMC population, and it is, therefore, not surprising that a complex pattern of gene expression changes was observed. Further verification of specific pathway changes in individual cell types is an obvious next step.