This is the first report of a Weighted gene co-expression network analysis in Echinococcus granulosus sensu stricto protoscolex stage. These networks are useful in understanding parasite basic biology and provide the baseline for further research. Immune evasion is a relevant topic in the host-parasite interaction, and the identified networks could provide new molecular targets in the medical treatment of this disease, such as the COP9 signalosome subunit.

Cystic Echinococcosis (CE) is a zoonotic larval disease, caused by the infection of the tapeworm Echinococcus granulosus sensu lato (s.l.), and affects both animals and humans globally. According to the World Health Organization (WHO), the global burden of trying to stop this zoonosis exceeds three billion US dollars every year [1]. CE is characterized by the development of cysts (formerly called hydatid cysts) in the viscera (mainly, liver and lungs) of intermediate hosts such as cattle and sheep, among other herbivores, whereas humans act as dead-end hosts [2]. CE cysts are comprised of three layers: the germinal and laminated layer of parasite origin, and an adventitial layer, which is the result of the host immune response against the parasite [3]. The germinal layer produces protoscoleces, the infective stage for the definitive host (dogs and other canids), who become infected when consuming viable protoscoleces [4]. E. granulosus s.l. is a cluster that groups different species of Echinococcus such as E. equinus, E. ortleppi, E. canadensis, E. felidis and E. granulosus sensu stricto (s.s.) [5]. The latter, which corresponds to the sheep-dog cycle, is the main E. granulosus s.l. species that infects humans [6], and it is found in Asia, Europe, Oceania and the Americas [7]. For causes that remain to be fully understood, some E. granulosus s.l. CE cysts are unable to produce viable protoscoleces and are termed non-fertile CE cysts [3]. Since cattle and sheep host can harbor simultaneously fertile and non-fertile CE cysts in the same viscera, from the same E. granulosus s.s. species [2], there must be microenvironment factors that affect parasite development and shape the host immune response.

Differences in morphology and fertility of CE cysts regarding infected organ have been previously described. In E. granulosus s.l. cysts from buffalo, the maturation rate of protoscoleces of liver origin was found to be different as compared to those of lung origin [8]. Also, a morphological study showed that larger protoscoleces were found in lungs CE cysts compared to those found in liver CE cysts, and presented size variability depending on the organ localization, among other variables [9]. So far, the identification of the molecular mechanisms that drive these differences remain unknown.

In a “classic approach” of gene expression profiles, the analyses focus on the individual genes, ignoring completely the interactions and expression minor variation among them, leaving a gap of information as genes play roles not by isolation but by interaction with each other. The gene-to-gene co-expression analysis emerged as an approach to solve the gene interaction problems [10]. Weighted gene co-expression network analysis (WGCNA) has been applied to many studies since late 2008, when the R package WGCNA was released [11], and then growing exponentially, being the 2012 the first year to reach 10 articles per year with this methodology and in the 2021 reaching more than 600 articles per year already published. Focusing on high complexity and multifactorial biological problems, WGCNA can be used as a data exploratory tool or as a gene screening method; having applications in gene expression and protein interaction data, among others [12]. Comparing gene expression networks (not only an individual gene) between set samples may identify modules of co-expressed genes and expression profiles of intramodular hub genes inside an interesting module, showing the Connectivity (the relative importance of a gene in a network). This approach may be performed to identifying this hub genes as possible therapeutic targets based on the constructed network [13, 14].

WGCNA has been already proposed to explore the parasite—host interaction. In human patients with chronic Chagas disease (Trypanosoma cruzi), hub genes were found associated with immune cell signaling pathways, T cell activation and B cell cellular immunity, although gene expression analysis of the parasite were not done [15]. In another instance, transcriptomic analyses of the freshwater snail Oncomelania hupensis, were performed after 3 different times of invasion with the parasite Schistosoma japonicum, finding a module related to ribosomes, translation, mRNA processing, among others, associated with Schistosoma infection [16]. Another work explored the presence of non-conserved modules and its hub genes between complicated and uncomplicated disease generated by protozoan Plasmodium falciparum, finding differences among those groups and identifying key genes for the development of the parasite [17]. Similar experiences have been performed in Trypanosoma brucei species, being those findings key in the potential identification of molecular targets to control the disease [18]. These findings show that WGCNA is a useful tool in the analysis of transcriptomic data in samples from parasite diseases, although most of them are focused on the host instead of the parasite.

Thus, it would not be unexpected that more complex statistical methods, such as the generation of co-expression networks, could help to elucidate subtle differences associated with parasite survival in different host organs. Therefore, to understand the differences in gene expression between E. granulosus s.s. protoscoleces found in liver and lung CE cysts, a WGCNA was performed. We found several modules with hub genes coding for hypothetical protein products and immune evasion related modules, revealing gene co-expression networks. Furthermore, we discuss possible molecular mechanisms performed by Echinococcus granulosus s.s. in response to adverse environments found in different organs that they can parasitize.

Despite morphological and pathophysiological differences in E. granulosus s.l. protoscoleces from lung and liver that have been previously documented [8, 24, 25], molecular processes associated with this issue has never been addressed based on next-generation sequence data. WGCNA is a method for the analysis of the gene expression patterns of multiple samples, clustering genes and form modules by similar gene expression patterns, creating co-expression networks and identifying intramodular hub genes (highly connected genes inside each module). The main differences among classical differential expression analysis and functional enrichment analysis, is that they cannot reveal connections and interactions among genes that are crucial in biological processes [26, 27]. A recent work involving 10 parasitic platyhelminths (strobilated and non-strobilated) species, identified a set of 34 evolutionary conserved cestode proteins, as possible components of developmental pathways required for strobilation, including E. granulosus s.l. inside the group of strobilated platyhelminths [28], This work is the only one to explore in this way some aspects of the E. granulosus s.s. development, but did not assess differences among parasites collected from different viscera. A previous transcriptome work did not identify differential gene expression among viscera in protoscoleces from sheep liver CE cyst and sheep lung CE cyst [23].

From the identified modules of this work, there are two interesting approaches for the discussion. The first one is to look for immune related hub genes in the modules, that can add information about the host-parasite relationship and/or immunoregulation mechanisms. The second one is the hypothetical proteins found as hub genes, as interesting targets for the study of specific and unknown mechanisms of parasitism in each organ.

On the first approach we may highlight the RoyalBlue, White and specially the MidnightBlue modules, in first two modules there is an immune related hub gene specific to E. granulosus s.l., which is composed of tegument protein and tegumentary antigen. The tegumentary antigen is an immunomodulatory molecule associated with chronic infection, as it inhibits chemotaxis, induces IL-4-positive T lymphocytes and non-complement fixing antibodies (IgG4) [29]. As a specific protein from E. granulosus s.l., this data suggest that suppression of a co-expressed gene may consequently also suppress the tegumentary antigen, but it has been demonstrated in our work that it will occur in liver CE cysts and could not happen in lung CE cysts.

In the MidnightBlue module, there are two genes related to protein degradation pathway, first, the ubiquitin hydrolase isozyme L3 (EGR_02941) has been reported before, highlighting that this protein can hydrolyze UBB(+1), (a form of ubiquitin associated with neurogenerative disorders) which is not effectively degraded by the proteasome [30]. On the other hand, we also found the COP9 signalosome complex subunit 3 (EGR_02738), an essential component of COP9 signalosome. This COP9 signalosome is a conserved molecule found in various disease-causing protozoans such as Leishmania spp., Trypanosoma spp., Toxoplasma spp., and Entamoeba histolytica. Moreover, in the latter, it was found to be essential for the normal functioning and a druggable parasite target, since disruption of the COP9 from the parasite, leads to a disrregulation of the ubiquitin proteasome pathway, impairing degradation of the proteins and leading to cell death [31]. Also, the COP9 complex has been involved in the immune response, regulating the NF-kB function by protecting the IκBα (an inhibitor of NF-kB) from degradation, leading to reduction of the NF-kB under TNF stimulation [32]. With those two Hub genes inside this module, we speculate that this non-conserved module represents an important part of the protein degradation pathway and could be useful for the parasite to suppress the immune response from the host. Interestingly, NF-κB signaling pathway has particular relevance to several liver diseases and hepatoprotective agents, which could explain why the MidnightBlue module would be specific to this organ [33].

In the same context of the immune response, inside the MidnightBlue module, it is important to note that the tetraspanin CD9 antigen (EGR_09512) acts also as a hub gene in this module. It is known that mononuclear phagocytes fuse to form multinucleated giant cells, a hallmark of the immune response in the adventitial layer of non-fertile CE cysts [34]., When these cells cannot eliminate pathogens, cell-mediated immunity is activated and mononuclear phagocytes coalesce to form Langhans giant cells [35]. A study investigated the role of this tetraspanin associated with CD81 tetraspanin, finding that CD9 and CD81 coordinately prevent the fusion of mononuclear phagocytes [36]. These results suggest a potential role of tetraspanins associations in the immune evasion of the host.

The methyl-CpG-binding protein Mbd2 (EGR_02905) is another important hub gene found in the MidnightBlue, in a study performed in bone marrow-derived dendritic cells (DC) from wild-type (WT) and Mbd2^−/− mice, the DC from Mbd2^−/− mice presented different mRNA expression, with 70 genes downregulated and 49 genes upregulated (compared with WT). Among these downregulated genes, many of them were related to immunological processes like antigen presentation, showing that this gene is able to induce control over other genes and also suggesting that the inhibition of Mbd2 may impair the optimal DC function and the initiation of a Th2 immune response (which is the main response against helminths [37]). It is curious that a parasite expresses a gene which is responsible of maintaining an adequate immune response against helminthes. We hypothesize that the parasite molecule can mimics host molecule and then evading the immune response, although more studies are needed to confirm this.

With these findings we may think that the MidnightBlue module is important in the immune evasion response and may represent possible therapeutic targets, as all hub genes inside the module are highly connected. This approach opens a window to new therapeutic strategies, as we found specific gene interaction occurring only in protoscoleces from liver but not in lung CE cysts. These differences should be taken into account as the generation of the cyst in the liver may involve different pathways than the occurring in the lung, so it is not misguided to think, as it seems to be a different process of cyst development, as different strategies to control the liver and lung echinococcosis.

Another interesting approach is that SkyBlue and Violet hub genes are exclusively Hypothetical proteins. It means that there are a group of unknown function proteins performing key roles in the parasitism of protoscolex in liver CE cysts, that are not co-expressed in lung CE cysts protoscoleces. This finding also suggests that the microenvironment promoted by the host induces changes in parasite gene expression, which could not be identified by more conventional transcriptomic analyses [23]. As no significative results of homology analysis, we may think that there are many genes of the unknown function that should be explored.

An important limitation of our work that is important to know, is that this approach needs mainly RNA-seq data that considers large number of replicates with low variation among them [28]. Despite this difficulties, choosing a topology close to a scale free-network seems to be a good model to work on, as despite universality of scale-free networks remains controversial, biological networks appear strongly scale free [38]. Conversely, while it is true that would be better to work with more samples, our data will be available for future exploration complemented with more samples obtained from other researchers.

A second approach is to find non-conserved modules in the protoscoleces from lung CE cysts compared to protoscoleces from liver CE cysts. Attempts to perform were unsuccessful, as it was not possible to get a scale-free network in protoscoleces from lung CE cysts, which could lead to a loss of information since the WGCNA analyzes are optimized for this type of network. However, the impossibility of obtaining the scale-free network only prevented an analysis of the conservation of the data of protoscoleces from liver CE cysts data for the identified modules in protoscoleces from lung CE cysts. Thus, data regarding the possible existence of new protoscoleces from lung CE cyst modules, which are not identified in protoscoleces from liver CE cysts, were not included in this work.

Nonetheless, the investigation into hub genes presents a promising avenue for translational research, potentially informing the development of novel therapeutic or control strategies against E. granulosus. A pertinent example is the work of Cancela et al. [39], who identified a novel subfamily of nuclear receptors with two DNA-binding domains (2DBDs)—a feature not yet reported in vertebrates. Their elucidation of the full-length 3D structure of the, Eg2DBDα.1 nuclear receptor in E. granulosus offers valuable insights into the receptor and structure-function relationship, highlighting its potential as a target for novel anthelmintic drugs.

This research contributes to the understanding of differential gene expression regulation in E. granulosus, underscoring the potential role of the organ-specific environment in exerting selective pressure on the parasite. Such environmental influences may have been pivotal in the evolution of certain E. granulosus s.l. lineages, enhancing their infective efficacy under specific conditions.

This work shows that there are differences in gene expression in protoscoleces from Cystic Echinococcosis cysts found in the liver and lungs. We found 34 modules in protoscoleces from liver CE cysts, 12 of which showed differential co-expression compared to protoscoleces from lung CE cysts. Several of these differentially co-expressed modules contain genes related to immune evasion as well as hypothetical proteins of unknown function that are co-expressed exclusively in protoscoleces from liver CE cysts. This suggests there are molecular differences in protoscoleces based on the organ environment of the CE cyst, which has implications for understanding and controlling CE infection.