Differentially expressed homologous genes reveal interspecies differences of Paragonimus proliferus based on transcriptome analysis

Freshwater crabs were collected from a natural stream in Mengla County, Xishuangbanna Dai Autonomous Prefecture, Yunnan Province, China. The shells of the crabs were removed, and the limbs, muscles, viscera, and other soft tissues were torn or triturated into tiny pieces, which were then repeatedly ground in a wire mesh. Samples were then rinsed by passing a large amount of water through the wire mesh into a 1000 ml tapered cylinder. The turbid liquid was allowed to settle for 30 minutes and the supernatant was carefully discarded, leaving approximately one third of the cylinder filled with sediment. After three to five washes, the remaining turbid liquid was shaken and poured into several petri dishes on a black backdrop. Enough suspension was poured into the dishes for tiny worms to easily be detected by visual inspection or under a light microscope. Using a pipette, the detected metacercariae were moved into normal saline and kept at 4°C until further use. The obtained metacercariae were used to infect Sprague-Dawley rats (eight metacercariae per rat) via subcutaneous injection (Lin & Xueming, 2001) into the abdominal wall. The rats were maintained at the Animal Lab Center of Kunming Medical University and provided with food and water ad libitum. Eight weeks post-infection (PI), the rats were euthanized under anesthesia by intraperitoneal injection of pentobarbital. Adult P. proliferus worms were isolated from the lungs of the rats, and a sample of three adult worms, whose viscera were removed under a light microscope, was stored at -80 °C until further use.

Transcriptome data of P. proliferus and other Paragonimus species Three P. proliferus adults were homogenized in 1 ml TRIzol reagent using a microcentrifuge pestle, and total RNA was extracted from the homogenate using a TRIzol Plus RNA Purification Kit (Thermo Fisher Scientific). RNA was treated with DNase I, and oligo(dT) was used to isolate mRNA. mRNA was fragmented by mixing with fragmentation buffer, and cDNA was synthesized using the mRNA fragments as templates. Short fragments were purified and resolved with elution buffer for end reparation and single nucleotide A (adenine) addition, followed by the addition of adapters. Suitable fragments were then selected for PCR amplification. An Agilent 2100 Bioanaylzer and an ABI StepOnePlus Real-Time PCR System (Applied Biosystems) were used to assess the quantity and quality of the sample library. The library was then sequenced using an Illumina Hi-Seq 4000 RNA-Seq platform with a depth of 12 Gb. The four available Paragonimus transcriptomic datasets were downloaded from the following NCBI database websites: P. skrjabini (https://www.ncbi.nlm.nih.gov/sra/SRX1507709)P. kellicotti (https://www.ncbi.nlm.nih.gov/sra/SRX530756)P. miyazakii (https://www.ncbi.nlm.nih.gov/sra/SRX1100074) and P. westermani (https://www.ncbi.nlm.nih.gov/sra/SRX1507710)

After obtaining raw reads of the five Paragonimus species, namely P. proliferus, P. skrjabini, P. kellicotti, P. miyazakii, and P. westermani, low-quality and adaptor-polluted reads and reads with a high content of unknown bases (N) were removed to obtain clean reads, which were then used to perform de novo assembly of unigenes using Trinity (version 2.0.6).

Transcript annotation: We used BLAST (version 2.2.23) to align unigenes to Non-Redundant Protein Sequence (NR), Nucleotide Sequence (NT), Clusters of Orthologous Groups (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Swiss-Prot annotations. Blast2GO (version 2.5.0) was used to obtain NR and Gene Ontology (GO) annotations, and InterProScan5 (version 5.11-51.0) was used to obtain InterPro annotations.

Four bilateral comparisons of unigenes were performed among the five species (P. proliferus versus P. skrjabini, P. kellicotti, P. miyazakii, and P. westermani) using BLAST. Transcripts that were the best match between two species (e-value=10e-5) were defined as homologous genes. Pairwise comparisons were also performed using NOISeq to identify differentially expressed homologous genes with |log2FC|≥1 and probability > 0.8. Mutually differentially expressed genes among the four pairs of comparisons were considered core genes.

Their extremely thin and fragile cystic wall, spindle or scaphoid shape, and huge size of 2.253 ± 0.364 mm × 0.668 ± 0.071 mm allowed metacercariae to easily and rapidly excyst after isolation from crabs and exposure to the environment. The ventral sucker, which was approximately three times the size of the oral sucker, was located in the anterior third end of the body. The two intestinal branches that extended alongside the body wall and towards the tail end were thin and curved before the ventral sucker level, and became thick and smooth thereafter (Fig. 1a).

Fig. 1

Eight weeks PI, adult worms (identified by mature reproductive organs that could be clearly visualized under a light microscope) of 7.238 ± 0.704 mm × 3.571 ± 0.655 mm in size were detected in the lung tissue or thorax of rats. An abundance of eggs filled an eiloid uterus, which was located right at the ventral sucker level. One ovary was located at the opposite side of the uterus, and the two testicles, located in the middle-posterior part of the body, were flamboyancy with 4 – 6 lobules. The uterus and testes were also extremely huge (Fig. 1b).

The morphological features of both metacercariae and adults were identical to those of P. proliferus described by Ben-jiang ( 2004) and Doanh et al. (2008; 2013).

As shown in Fig. 2, RNA-Seq was used to sequence the transcriptome of adult P. proliferus worms. The sequences were then assembled de novo and filtered to include only high-confidence transcript sequences. Of 89.34 Mb clean reads obtained from 119.20 Mb raw reads, a total of 47,959 transcripts corresponding to 29,967 unigenes were generated after removal of low-quality and adaptor-polluted reads and reads with a high content of unknown bases (N). The total length, mean length, N50, and GC content were 17,089,849 bp, 570 bp, 826 bp, and 47.36 %, respectively. Overall, a total of 20,669 (68.97 %) unigenes were annotated to seven functional databases (NR, NT, GO, COG, KEGG, Swiss-Prot, and InterPro). Fig. 3 shows that the annotated unigenes belonged to Opisthorchis viverrini (26.25 %), Culex quinquefasciatus (22.36 %), Clonorchis sinensis (20.83 %), Halyomorpha halys (3.7 %), and other species (28.86 %).

Fig. 2

Fig. 3

As shown in the Venn diagram in Fig. 4, a total of 10,629 (19.13 % of P. proliferus genes), 12,074 (21.74 %), 13,558 (24.41 %), and 8,051 (14.49 %) homologous genes were identified between P. proliferus and P. skrjabini, P. kellicotti, P. miyazakii, and P. westermani, respectively.

Fig. 4

To characterize the interspecies differences in homologous gene expression, pairwise comparisons were performed using NOISeq, which identified 8192 differentially expressed homologous genes among the five species, 7393 of which had |log2FC| ≥ 1 and probability > 0.8. Surprisingly, as shown in Fig. 5, 3950/5622 (70.26 %), 1049/1084 (96.77 %), 388/473 (82.03 %), and 189/214 (88.32 %) genes were expressed at lower levels in P. proliferus compared to P. kellicotti, P. skrjabini, P. westermani, and P. miyazakii, respectively. The biological characteristics of P. proliferus related to lower expression of these genes requires further exploration, but there is indeed a possibility that these genes influence growth, development, invasion, reproduction, virulence, and motility.

Fig. 5

Fig. 6

There were 527, 693, and 757 genes annotated as cellular components (CC; 207 terms, 18 of which were significantly enriched with P < 0.05), biological processes (BP; 930 terms, 10 with P < 0.05), and molecular functions (MF; 756 terms, 14 with P < 0.05), respectively. The top 42 enriched GO terms (P < 0.05) are shown in Fig. 6.

Of the 10 significantly enriched GO terms belonging to BP, four participate in “genetic central dogma”, such as DNA replication/ synthesis of RNA primers, regulation of translational initiation, rRNA modification, and DNA-dependent DNA replication. The remaining six terms (monovalent inorganic cation transport, potassium ion transport, hydrogen ion transmembrane transport, ornithine metabolic process, and fructose metabolic process) are related to energy metabolism. These findings are in line with the enrichment of MF; at least five of the 14 terms (P < 0.05) take part in “duplication, transcription, or translation”, such as endonuclease activity, DNA primase activity, ribonuclease activity, and transcriptional repressor activity. The majority of the remaining terms were related to the oxidation respiratory chain and the activity of some energy production-related enzymes. Therefore, we inferred that changes in the processes of “duplication, transcription, or translation” and “energy metabolism” in P. proliferus may lead to key interspecies differences in development or biological behavior. This was further verified by the significantly enriched GO terms belonging to CC; seven of 18 terms (P < 0.05) belonged to chondriosome-related cellular components, four of 18 to ATPase, and three of 18 to ribosome. As is well known, the chondriosome is the site of the oxidation respiratory chain, and the ribosome is critical for protein synthesis.

A total of 1412 genes were annotated to 290 KEGG pathways (as shown in Table 1, 13 pathways were significantly enriched with P < 0.05). Of the 13 pathways, four were assigned to signal transduction, including forkhead box, class O (FoxO), transforming growth factor (TGF) β, vascular endothelial growth factor (VEGF), and mechanistic target of rapamycin (mTOR) signaling, all of which participate in important biological processes; two were related to metabolism, mainly pyruvate metabolism and terpenoid backbone biosynthesis; and two (aldosterone-regulated sodium reabsorption and bile secretion) were assigned to the excretory and digestive systems. However, these genes also intrinsically belong to ion or lipid metabolism as well as pathways related to DNA replication and repair.

In mammals, FoxO family members are involved in cell metabolism, growth, differentiation, oxidative stress, senescence, autophagy, and aging (Lee & Dong, 2017), and we suspect that they have similar biological functions in Paragonimus. TGF-β or its parasite mimics may regulate the immune response of the host (Musah-Eroje & Flynn, 2018), and mTOR signaling is related to Leishmania proliferation (Jaramillo et al., 2011). Energy metabolism is closely related to parasitic growth or adaptation to the host environment (Caddigan et al., 2017), as well as replication and virulence (Mancio-Silva et al., 2017). Recently, it was reported that pyruvate homeostasis is a determinant of parasite growth and metabolic plasticity in Toxoplasma gondii (Xia et al., 2019).

The processes of signal transduction, metabolism, and DNA replication and repair may affect the biological characteristics of P. proliferus; further functional verifications should be performed to determine the effect of genes on the key targets of these signaling pathways.

As shown in the Venn diagram in Fig. 7, of the 8192 differentially expressed homologous genes, 49 mutually differentially expressed genes were identified among the four pairs of comparisons and defined as core genes that may play key roles in the biological differences between P. proliferus and other Paragonimus species. Interestingly, a majority (48/49) of these core genes were expressed at lower levels in P. proliferus than in the four other species. A total of 16 of the 49 core genes were annotated by functional databases (partial annotation information is listed in Table 2), whereas the remaining 33 were not.

Fig. 7

Fig. 8

Six core genes were annotated to 11, 21, and 26 GO terms belonging to CC (DNAH1, NOP56, and PI3K), BP (SSP7, PFK-2/ FBPase2, DNAH1, Adk, PI3K, and NOP56), and MF (SSP7, PFK-2/FBPase2, DNAH1, Adk, PI3K, and NOP56), respectively, as shown in Fig. 8. Of the 21 terms belonging to BP, 17 were involved in metabolism, and at least 12 of the 26 nucleotide or ribonucleotide-related terms belonging to MF may be involved in DNA replication and repair. Of all the enriched GO terms, four belonging to BP (phosphate-containing compound metabolic process, organophosphate metabolic process, phosphorus metabolic process, and carbohydrate derivative metabolic process) and three belonging to MF (phosphotransferase activity/alcohol group as acceptor, kinase activity, and transferase activity/transferring phosphorus-containing groups) were significantly enriched with P < 0.05, which seems to indicate that phosphorus and carbohydrate metabolism may play an important role in the interspecies differences between P. proliferus and other Paragonimus species.

The six genes annotated to the GO terms mentioned above, as well as four others annotated to the KEGG Orthology (MRP, Tuba, Wasf2, and Ppp1r7), were also annotated to 87 KEGG pathways, 28 of which had a P value < 0.05 (Table 3) and all of which were enriched with four genes (PI3K, WASF2, MRP, and PFK-2/FB-Pase2). Of the 28 pathways, four were assigned to signal transduction (including 5’ adenosine monophosphate-activated protein kinase [AMPK], Janus kinase/signal transduces and activators of transcription [JAK-STAT], VEGF, and mTOR signaling), four participated in immune-related signaling (Fc gamma R-mediated phagocytosis, and Toll-like receptor, Fc epsilon RI, and B cell receptor signaling pathways), and seven were related to digestive or endocrine systems, as well as carbohydrate metabolism, and signaling pathways that are always involved directly or indirectly in metabolism. It is worth noting that there were also nine pathways related to level 2 cancers according to KEGG enrichment, which suggests that parasite infections may be a risk factor for tumorigenesis, a view that is supported by other reports (Feng & Cheng, 2017; van Tong et al., 2017).

Phosphoinositide 3-kinase (PI3K), which is critical for parasite virulence, was involved in 26 significantly enriched pathways. In Entamoeba histolytica, PI3K signaling plays an important role in motility, phagocytosis, host cell adhesion, and proteolysis of the extracellular matrix (Koushik et al., 2014). In Trypanosoma cruzi, PI3K-like activity plays an important role in host cell invasion, and parasite entry is impaired in trypomastigotes treated with a PI3K inhibitor prior to infection (Wilkowsky et al., 2001). Recently, PIKs have been proposed to be targets for the treatment of T. cruzi infection (Gimenez et al., 2015). Phosphofructokinase (PFK-2)/ fructose 1,6-bisphospatase (FBPase), which is also found in Trypanosoma brucei, Leishmania major, and T. cruzi (Chevalier et al., 2005), was identified as being involved in fructose and mannose metabolism and thyroid hormone signaling in our study.

According to KEGG enrichment of the analyzed transcriptome, Wiskott-Aldrich syndrome protein family member 2 (WASF2), which was annotated to signaling of choline metabolism in cancer, Fc gamma R-mediated phagocytosis, and bacterial invasion of epithelial cells, may be involved in parasite metabolism, invasion, and host immunity. Mitochondrial RNA processing (MRP) belongs to signaling of microRNAs in cancer and ATP-binding cassette (ABC) transporter classes. Reports have revealed that disruption of MRP in Plasmodium falciparum was related to a reduction in parasitemia, removal of toxic metabolites, and transport of antimalarial drugs (Raj et al., 2009). Additionally, drug efflux transporters of the ABC gene superfamily can confer drug resistance to pathogens (Carmona-Antonanzas et al., 2015).

Tubulin alpha-8 chain (TUBA) is also essential for host cell invasion and parasite survival in T. gondii (Morrissette, 2015). Adenosine kinase, encoded by the Adk gene, is a key enzyme in the purine-salvage pathway, adenosine phosphorylation, and activation of cytotoxic analogues (Romanello et al., 2013; Timm et al., 2014). In T. brucei, inhibition or downregulation of adenosine kinase results in resistance to cordycepin, demonstrating its role in the activation of adenosine antimetabolites (Luscher et al., 2007) and silencing of nucleolar protein (NOP) 1, which contributes to defects in rRNA-processing and causes lethal modifications (Barth et al., 2008). However, it is unclear whether NOP56 has a similar role. Protein phosphatases (PPs) are also related to reproduction and development. For example, in T. cruzi protein phosphatase 2A is important for the complete transformation of trypomastigotes into amastigotes (Gonzalez et al., 2003), and recent research in Toxocara canis has shown that PP1 regulates kinetochore-microtubule interactions during spermatogenesis via PP1cα-PP1r7 mechanisms (Ma et al., 2015). The functions of dynein axonemal heavy chain 1 (DNAH1) and Ssp7 are still unknown in parasites, especially in P. proliferus, and further research on the molecular biology of these genes in parasites is necessary.

Serine/threonine-protein phosphatase 2A regulatory subunit B (STP) may play a functional role in parasite reproduction (Boag et al., 2003; Ma et al., 2014). In T. canis, STP is not detected in adult females, but is expressed at high levels in the spermary, vas deferens, and musculature of the adult male, suggesting that STP may be involved in spermatogenesis and mating behavior and may represent a potential molecular target for controlling the reproduction of T. canis (Ma et al., 2014). The outer dense fiber (ODF), a coiled-coil protein, is the main cytoskeletal structural component of the sperm tail in animals and is potentially involved in building the cellular cytoskeleton (Petersen et al., 2002). Some ODF genes were found to be aberrantly expressed in tumors, such as prostate adenocarcinomas (Ghafouri-Fard et al., 2012), but the functions of parasite ODF, retrovirus-related Pol polyprotein from transposon opus, STP, and the cell wall protein Awa1p, require further research.