Transcriptome analysis uncovers the key pathways and candidate genes related to the treatment of Echinococcus granulosus protoscoleces with the repurposed drug pyronaridine

Background Cystic echinococcosis (CE) is a life-threatening zoonosis caused by the larval form of Echinococcus granulosus tapeworm. Our previous study showed that an approved drug pyronaridine (PND) is highly effective against CE, both in vitro and in an animal model. To identify possible target genes, transcriptome analysis was performed with E. granulosus sensu stricto protoscoleces treated with PND. Results A total of 1,321 genes were differentially expressed in protoscoleces treated with PND, including 541 upregulated and 780 downregulated genes. Gene ontology and KEGG analyses revealed that the spliceosome, mitogen-activated protein kinase (MAPK) pathway and ATP-binding cassette (ABC) transporters were the top three enriched pathways. Western blot analysis showed that PND treatment resulted in a dose-dependent increase in protein expression levels of EgMKK1 (MKK3/6-like) and EgMKK2 (MEK1/2-like), two members of MAPK cascades. Interestingly, several heat shock protein (HSP) genes were greatly downregulated including stress-inducible HSPs and their constitutive cognates, and some of them belong to Echinococcus-specific expansion of HSP70. Conclusions PND has a great impact on the spliceosome, MAPK pathway and ABC transporters, which may underline the mechanisms by which PND kills E. granulosus protoscoleces. In addition, PND downregulates HSPs expression, suggesting a close relationship between the drug and HSPs. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07875-w.


Background
Cystic echinococcosis (CE) is a parasitic zoonosis caused by the larval stage of the dog tapeworm Echinococcus granulosus. Chronically infected humans or domestic animals remain asymptomatic for a long time. Infection with E. granulosus leads to the development of one or multiple cysts located mostly in the liver and lungs, which triggers clinical signs in the late stage, including abdominal and chest pain, chronic cough, vomiting, even death [1,2]. CE is a globally distributed disease highly endemic in South America, Northern Africa, and Central Asia (especially Western China) [3]. The global health burden of CE is estimated at over 1 million disability-adjusted life-years (DALYs) each year [4]. The disease also affects the local livestock industry's economic benefits, with an estimated yearly loss of at least US$3 billion [5].
Among current CE treatment options, anti-parasitic drug therapy is widely used for most clinical cases [6,7]. Recently, we repurposed an approved anti-malarial drug pyronaridine (PND) as a promising candidate for CE treatment [8]. Oral administration of PND showed high concentrations of the drug in the liver and lungs, which are the most affected organs in CE. Oral administration, intraperitoneal injection or microinjection procedure (which mimics the clinical percutaneous techniques) significantly reduced the parasite burden in mice. However, the anti-CE mechanism of action of PND is not clear. Previous studies showed that the primary anti-malarial mode of action of PND is inhibition of β-hematin formation, enhancement of hematin-induced red blood cells lysis, and inducing the formation of abnormal vesicles in the food vacuole of plasmodium [9,10]. A transcriptome profiling of Plasmodium falciparum in response to PND reveals a striking abundance of genes encoding hostexported proteins [11]. In addition, PND has been characterized as a potential anticancer agent, which reverses the multi-drug-resistance (MDR) phenotype in MDR cancer cell lines by inhibiting the function of the efflux P-glycoprotein (Pgp) [12,13]. In this study, to obtain a comprehensive understanding of the anti-CE mechanism of PND, the global gene expression in E. granulosus protoscoleces (PSCs) following treatment with PND was analyzed using RNA-seq.

RNA sequencing data analysis
Global gene expression of PND-treated E. granulosus PSCs was analyzed using an Illumina platform. The obtained sequences were aligned against E. granulosus genome sequences. A total of 60.3 and 58.2 million clean reads were obtained from control and PND groups, respectively ( Table 1). The clean reads were mapped to the E. granulosus genome scaffold (https://www.ncbi. nlm.nih.gov/genome/?term=Echinococcus+granulosus) reported by Zheng et al. [14]. As shown in a volcano plot (Fig. 1), a total of 1,321 genes were significantly differentially expressed in the PND-treated group compared to the control group, including 541 upregulated and 780 downregulated E. granulosus genes.

Gene ontology (GO) classification
There were 47 GO terms significantly enriched for differentially expressed genes (DEGs) in PND-treated groups, including 19 biological process terms, 17 cellular component terms, and 11 molecular function terms ( Table 2). The most enriched GO terms were associated with cell part, voltage-gated calcium channel complex, binding, voltage-gated calcium channel activity and response to stress.

Biochemical pathway
The DEGs were also mapped to five KEGG subsystems, including environmental information processing, genetic information processing, organismal systems, metabolism, and cellular processes. The most significantly enriched 10 pathways are shown in Fig. 2. Spliceosome, mitogenactivated protein kinase (MAPK) signaling pathway and ATP-binding cassette (ABC) transporters were the top three significantly enriched pathways (Figs. S1, S2, S3). PND treatment resulted in significant changes in seven ABC transporter genes (Table 3), and three of them (EGR_00511, EGR_00512 and EGR_01347) are involved in MDR, coding MDR transporters and MDR-associated proteins. In addition, we found that PND treatment induced a significant downregulation in a large number of HSP genes including heat shock proteins (HSPs) and heat shock constitutive cognates (HSCs) ( Table 4). Most of the genes are of HSP70 family and involved in the top two enriched pathways. HSP70 is a highly expanded gene family in Echinococcus spp. Among the 19 differentially expressed HSP genes, six genes belong to Echinococcus-specific expansion of HSP70 (Table 4).

Validation of key DEGs by qRT-PCR
To validate the results of transcriptome sequencing, quantitative PCR (qRT-PCR) was used for confirmation of eleven representative genes selected from the top three most enriched KEGG pathways, including HSP72, SRF, ECSIT, PKC, MP3K and PTP from the MAPK signaling pathway, SYF, LSM4 and U2AF from spliceosome, and ABCB1 and ABCG2 from ABC transporters (Fig. 3). The qRT-PCR expression patterns of nine out of eleven DEGs were in agreement with the results of the transcriptome analysis, despite the variation of drug concentration or treatment time.

The protein levels of EgMKK1 and EgMKK2
To further evaluate the changes of key members of MAPK cascades, we determined the protein levels of EgMKK1 (MKK3/6-like) and EgMKK2 (MEK1/2-like) by Western blotting. As shown in Fig. 4, generally, PND treatment upregulated the protein levels of both EgMKK1 and EgMKK2 in a dose-dependent pattern. For EgMKK1, a significant elevation was observed in PSCs following the treatment of PND at the concentrations of LC 30 and LC 50 compared with the control group (p < 0.05), while in the case of EgMKK2, it was the treatment of LC 50 (p < 0.01).

Discussion
CE is a neglected disease that has remained "unattractive" to pharmaceutical companies, considering that, in the last few decades, no alternative drugs have been approved for treating it, although some efforts have been made [15][16][17]. In our previous study, PND, an approved antimalarial drug, was repurposed as an anti-CE candidate. PND killed 100 % of the cysts in a mouse infection model by intraperitoneal injection at 57 mg/kg/day for three days. When administered orally with a regimen of 57 mg/kg/day × 30 days, it produced 90.7 % cyst mortality, showing that PND is much more effective than albendazole (22.2 % cyst mortality at 50 mg/kg/day), the only anti-CE drug recommended by WHO [8]. However, the anti-parasitic mechanism of PND remains unclear.
In this study, RNA-seq technology was used to explore the genes affected by PND on E. granulosus. Using a suitable low dose, our study showed that PND treatment induced changes in the expression of a large number of genes, including 541 E. granulosus genes upregulated and 780 downregulated, which demonstrates that PND effectively targets E. granulosus PSCs. GO enrichment and KEGG analyses revealed that the significantly PND-altered biological processes and pathways were associated with a wide range of cellular components, biological processes, and metabolic pathways, including cellular structures and signaling pathways.
The MAPK cascade is an evolutionarily conserved signal transduction pathway that transmits and converts Fig. 1 Volcano plot for the differentially expressed genes (DEGs) between control and pyronaridine (PND)-treated groups. The X-axis shows the differences in gene expression (FDR: adjusted p-value), while the Y-axis indicates expression changes (log2 fold change) of the genes in different groups. Splashes were for different genes. Black dots, red and yellow dots, and blue and dark blue dots represent genes with no significant discrepancy, significantly upregulated genes, and significantly downregulated genes, respectively many extracellular signals by three consecutive phosphorylation events. MAPK pathways are implicated in several cellular processes, including proliferation, differentiation, apoptosis, inflammation, and stress response [18][19][20]. According to KEGG enrichment analysis, the MAPK pathway comes in second in the top-changed pathways affected by PND. In the last decade, a few components of the MAPK pathway have been identified in E. granulosus, including Egp38, EgERK, EgMKK1 and EgMKK2 [21][22][23]. Meanwhile, some MAPK inhibitors (e.g., sorafenib, U0126, SB202190) were found to effectively kill E. granulosus in vitro and/or in vivo [21,24,25], proving that the key kinases could be used as potential targets for anti-CE drug development. Following the exposure of PND, a direct effect on the gene levels of the key nodes of the MAPK pathway was not observed. We speculate that, rather than specifically targeting one key node, PND likely had a general impact on the whole pathway, which was demonstrated by the significantly elevated protein levels of EgMKK1 and EgMKK2.
ABC transporters are transmembrane proteins that actively mediate the translocation of a wide variety of molecules across the cell membrane, including drugs. A subset of ABC transporters is closely linked to MDR, e.g. the best-characterized multidrug transporter Pgp (ABCB1/MDR1). ABC multidrug transporters have been implicated in drug resistance in several parasites [26,27]. In the genome of E. granulosus, 22 putative ABC transporters were identified and could be classified into six subfamilies [28]. In this study, PND treatment  (Table 3), and three of them (EGR_00511, EGR_ 00512 and EGR_01347) are involved in MDR, coding MDR transporters and MDR-associated proteins. Usually, anti-parasitic drug treatment would result in increased/over expression of ABC transporters, especially MDR transporters, to remove or exclude xenotoxins from cells to guard the normal cellular physiology [29][30][31]. While in this study, after the PND treatment, an ABCB gene (EGR_00511) was significantly downregulated (validated by qRT-PCR, Fig. 3). This indicates that, besides the strong protoscolecidal ability, PND could also negatively regulate the expression of E. granulosus MDR transporter to favor its retention in PSC tissues as an add-on effect. It inspires us that future research efforts could be geared towards the combination of MDR modulators and current anthelmintics to enhance drug susceptibility.
In addition, we found that PND downregulated several HSP genes in E. granulosus PSC. HSPs are originally identified because of their roles in response to heat shock (or other stressors) and these molecules are also molecular chaperones involved in protein folding and maturation [32]. Some HSPs (such as heat shock cognate proteins, HSCs) are constitutively expressed in cells, and serve vital functions in cell metabolism maintenance. We showed that PND treatment induced a significant downregulation in a large number of HSP genes, indicating a close relationship between the drug and HSPs. The differentially expressed HSP genes (Table 4) included five downregulated and three upregulated heat shock proteins and also eleven downregulated constitutive cognates, indicating that not all the DEGs observed in transcriptome analysis were necessarily induced by stress (e.g. an external stimulus caused by PND drug treatment). In addition, most of the HSP genes are of the   HSP70 family. HSP70 is a highly expanded gene family in E. granulosus [14,33], and it has been found that some HSP70s may be non-functional transcribed pseudogenes [34]. Through orthology search, six differentially expressed HSP genes were identified to belong to Echinococcus-specific expansion of HSP70 (Table 4). HSPs are implicated in the cause and progression of various diseases, such as infections [35], cancer [36,37], and neurodegeneration [38,39]. In parasites, such as Plasmodium spp. [40], Leishmania spp. [41], and Trypanosoma spp. [42], HSPs have been investigated as potential drug targets. Some of the HSP genes have already been identified and studied in E. granulosus [43,44]. The results here reported motivate us to study the relationships of PND and HSPs further.

Conclusions
In this study, the transcriptome landscape of E. granulosus PSCs treated with PND was characterized. This allowed the identification of 1,321 DEGs, some of which were found to exhibit great influence on various life processes of E. granulosus, including MAPK pathway, ABC transporters and HSPs. These findings provide valuable genetic data to facilitate future studies toward understanding the anti-parasitic mechanism of PND.

RNA-seq bioinformatics analysis
The raw reads were subjected to adapter trimming and low-quality filtering using SeqPrep (https://github.com/ jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/ sickle) with default parameters. The high-quality clean reads were aligned to the reference genome using the HIASAT (https://ccb.jhu.edu/software/hisat2/index. shtml) software. The DEGs between control and PNDtreated PSCs were identified based on fragments per kilobases per million reads (FPKM) using the Transcripts Per Million reads (TPM) method. DESeq2 (http:// bioconductor.org/packages/stats/bioc/DESeq2/) was used to identify differential expression analysis. Gene expression with log2 fold change ≥ 1 or ≤ − 1, and differences in expression with an adjusted p-value < 0.05 were considered to be significant. In addition, the GO and KEGG databases were explored to identify which DEGs were significantly enriched in GO terms and KEGG pathways. GO functional enrichment was carried out by Goatools (https://github.com/tanghaibao/Goatools) [45].

qRT-PCR assay
Eleven representative genes (four upregulated genes: PKC, MP3K, PTP and ABCG2; nine downregulated genes: HSP72, ECSIT, SRF, SYF, LSM4, U2AF, and ABCB1) were selected from the top three enriched pathways, and their gene expression levels in the control and PND-treated groups were evaluated. GAPDH was used as an endogenous control. Gene expression was  Table S1.

Western blot analysis
The rabbit anti-EgMKK1/anti-EgMKK2 serums are generous gifts from Dr. Chuanshan Zhang, the First Affiliated Hospital of Xinjiang Medical University. Western blot analyses of EgMKK1 and EgMKK2 were performed as previously described [21]. β-Actin served as a loading control.