Next Article in Journal
Identification and Functional Characterization of Alternative Transcripts of LncRNA HNF1A-AS1 and Their Impacts on Cell Growth, Differentiation, Liver Diseases, and in Response to Drug Induction
Previous Article in Journal
lncRNA-mRNA Co-Expression and Regulation Analysis in Lung Fibroblasts from Idiopathic Pulmonary Fibrosis
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
ncRNA
Volume 10
Issue 2
10.3390/ncrna10020027
ncrna-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Long Non-Coding RNA Levels Are Modulated in Schistosoma mansoni following In Vivo Praziquantel Exposure

by
Pedro Jardim Poli
1,
Agatha Fischer-Carvalho
1,
Ana Carolina Tahira
1,
John D. Chan
2,
Sergio Verjovski-Almeida
1,3 and
Murilo Sena Amaral
1,*
1
Laboratório de Ciclo Celular, Instituto Butantan, São Paulo 05503-900, SP, Brazil
2
Global Health Institute, University of Wisconsin-Madison, Madison, WI 53792, USA
3
Instituto de Química, Universidade de São Paulo, São Paulo 05508-900, SP, Brazil
*
Author to whom correspondence should be addressed.
Non-Coding RNA 2024, 10(2), 27; https://doi.org/10.3390/ncrna10020027
Submission received: 30 January 2024 / Revised: 5 April 2024 / Accepted: 13 April 2024 / Published: 19 April 2024
(This article belongs to the Section Long Non-Coding RNA)
Browse Figures
Versions Notes

Abstract

:
Schistosomiasis is a disease caused by trematodes of the genus Schistosoma that affects over 200 million people worldwide. For decades, praziquantel (PZQ) has been the only available drug to treat the disease. Despite recent discoveries that identified a transient receptor ion channel as the target of PZQ, schistosome response to this drug remains incompletely understood, since effectiveness relies on other factors that may trigger a complex regulation of parasite gene expression. Long non-coding RNAs (lncRNAs) are transcripts longer than 200 nucleotides with low or no protein-coding potential that play important roles in S. mansoni homeostasis, reproduction, and fertility. Here, we show that in vivo PZQ treatment modulates lncRNA levels in S. mansoni. We re-analyzed public RNA-Seq data from mature and immature S. mansoni worms treated in vivo with PZQ and detected hundreds of lncRNAs differentially expressed following drug exposure, many of which are shared among mature and immature worms. Through RT-qPCR, seven out of ten selected lncRNAs were validated as differentially expressed; interestingly, we show that these lncRNAs are not adult worm stage-specific and are co-expressed with PZQ-modulated protein-coding genes. By demonstrating that parasite lncRNA expression levels alter in response to PZQ, this study unravels an important step toward elucidating the complex mechanisms of S. mansoni response to PZQ.

1. Introduction

Schistosomiasis is a neglected tropical disease caused by trematodes of the genus Schistosoma. Conservative estimates suggest that schistosomiasis affects over 200 million people worldwide; the disease has been reported in 78 countries [1,2,3]. Additionally, according to WHO’s Global Health Observatory data from 2020, schistosomiasis ranks among the top five leading causes of DALYs (disability-adjusted life years) and death by parasitic and vector diseases [4]. The three main species that infect humans are Schistosoma haematobium, S. japonicum, and S. mansoni. The latter is responsible for infection in the Americas and sub-Saharan Africa [1,5]. S. mansoni worms usually inhabit the mesenteric veins of their mammalian host, where males and females are paired. The females lay hundreds to thousands of fertilized eggs per day, which are carried by the circulation to the liver or pass through the intestine wall, being eliminated with the host’s feces.
Since the mid-1980s until today, treatment of schistosomiasis relies on praziquantel (PZQ), which is considered a safe and effective drug, showing only mild side effects in humans, that eliminates adult schistosome infections, often in a single dose [5,6,7,8]. For decades it has been known that PZQ affects the tegument and induces calcium influx into the whole parasite, resulting in muscular contraction and paralysis [6]. However, only very recently, the target of PZQ in S. mansoni has been successfully identified and characterized as a transient receptor potential melastatin ion channel, known as Sm.TRPMPZQ [9,10]. Nevertheless, worm response to this drug is still not fully understood since it depends on a complex disruption of the parasite homeostasis. This involves the host’s immune system [1,6,7] in addition to the regulation of gene expression in the parasite when exposed to PZQ, as shown by McCusker et al. [11].
Long non-coding RNAs (lncRNAs) are transcripts longer than 200 nucleotides with low or no protein-coding potential [12,13]. In mammals, lncRNAs play a well-established role in numerous biological processes, such as gene regulation [14,15] and drug resistance [16,17]. Remarkably, lncRNAs exhibit high tissue- and species-specificity [13,14,18], hinting towards their potential as therapeutic targets against parasites [19,20]. In S. mansoni, the presence of lncRNAs was initially described by our group in 2011 [21] and has since been reported by other groups [22] as reviewed by Silveira et al. [19]. Furthermore, we have previously demonstrated that S. mansoni lncRNA levels are modulated by 5-azacytidine [23], a drug used to treat human myelodysplastic syndrome, and more recently, we highlighted their indispensable role in adult worm homeostasis and fertility [24].
Given the intricate nature of schistosomes’ response to PZQ [6,11], which remains incompletely understood, and recognizing the significant functions and regulatory roles exhibited by S. mansoni lncRNAs [24], we hypothesized that S. mansoni lncRNAs could be involved in parasite response following in vivo PZQ exposure. In this study, we demonstrate that in vivo exposure to PZQ induces significant changes in the expression of the worm’s lncRNAs, even several hours after the treatment. In our re-analyses of the RNA-Seq data of McCusker et al. [11], we found hundreds of lncRNAs differentially expressed both in mature and immature worms. Through RT-qPCR we validated the differential expression of selected lncRNAs, underscoring the impact of PZQ on altering S. mansoni lncRNA expression. This represents an important step toward fully elucidating the complexity of drug response in schistosomes, a crucial effort considering the emerging concern regarding PZQ-resistant worms [1,6,8].

2. Results

2.1. Sets of LncRNAs Are Differentially Expressed in S. mansoni following In Vivo Praziquantel Exposure

We re-analyzed the RNA-Seq data generated by McCusker et al. [11] to search for long non-coding RNAs (lncRNAs) differentially expressed (DE) in S. mansoni after in vivo treatment of infected mice with praziquantel (PZQ) (Supplementary Table S1 shows the samples used and alignment statistics). Three RNA-Seq datasets were re-analyzed: (i) RNA-Seq from adult worm couples obtained from mice across ten time points at 0, 0.25, 1, 3, 6, 9, 12, 24, 48, and 96 h after treatment with a single, curative dose of PZQ (400 mg/kg) delivered by oral gavage at 7 weeks post-infection (“time-course experiment”); (ii) RNA-Seq from adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 4 weeks post-infection (“4-weeks experiment”); (iii) RNA-Seq from adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 7 weeks post-infection (“7-weeks experiment”).
As expected, principal component analysis (PCA) resulted in transcriptomes of the PZQ-treated and control groups segregating broadly into two distinct regions with replicates from the same conditions clustering together, both for their control or PZQ-treated groups (Supplementary Figure S1).
In that study, McCusker et al. [11] analyzed only the protein-coding genes differentially expressed in S. mansoni after PZQ treatments. As lncRNA expression levels have been shown to be modulated or potentially affected by drugs in other eukaryotes [25] and also in S. mansoni [23], we hypothesized that lncRNA levels could be also modulated by PZQ in S. mansoni. Indeed, re-mapping the RNA-Seq data to the genome using a S. mansoni reference transcriptome previously built by us that comprises lncRNAs in addition to protein-coding genes [26], we found lncRNAs DE in S. mansoni worms recovered from mice treated with PZQ in all three re-analyzed RNA-Seq experiments: 189, 181 and 348 lncRNAs were found DE in the “time-course experiment” (Figure 1), in the “4-weeks experiment” (Figure 2a) and in the “7-weeks experiment” (Figure 2b), respectively.
In the “time-course experiment” the following lncRNAs were found as DE (Figure 1 and Supplementary Table S2): 127 long intergenic non-coding RNAs (lincRNAs, being 101 upregulated and 26 downregulated), 58 antisense non-coding RNAs (SmLNCAs, being 41 upregulated and 17 downregulated) and 4 sense non-coding RNAs (SmLNCSs, being 2 upregulated and 2 downregulated). In the “4-weeks experiment”, a higher number of SmLNCSs was found as DE when compared with the “time-course experiment” (Figure 2a and Supplementary Table S3): 108 lincRNAs, being 61 upregulated and 47 downregulated, 57 SmLNCAs, being 44 upregulated and 13 downregulated and 16 SmLNCSs, being 13 upregulated and 3 downregulated. Interestingly, in the “7-weeks experiment”, the highest number of lncRNAs was found as DE (Figure 2b and Supplementary Table S4): 236 lincRNAs, being 163 upregulated and 73 downregulated, 95 SmLNCAs, being 51 upregulated and 44 downregulated and 17 SmLNCSs, being 13 upregulated and 4 downregulated.
Protein-coding genes (Smps) were also detected as DE in the “time-course experiment” (Supplementary Figure S2), in the “4-weeks experiment” (Supplementary Figure S3) and in the “7-weeks experiment” (Supplementary Figure S4). As expected, we found a good overlap in the protein-coding genes found as DE when we compared our RNA-Seq re-analyses with the analyses from McCusker et al. [11]. For the “time-course experiment”, we found in our re-analysis 760 (41%) out of the 1848 Smps previously found as DE in McCusker et al. [11]. For the “4-weeks” dataset, we found in our re-analysis 112 (83%) out of the 135 Smps previously found as DE. For the “7-weeks” dataset, we found in our re-analysis 356 (88%) out of the 405 Smps previously found as DE. These are reasonable proportions considering the differences in the reference transcriptomes used and the different read-mapping and counting software used in both analyses, with more stringent DE statistical cutoff criteria used here.

2.2. A Set of LncRNAs Is Concomitantly Differentially Expressed under Distinct PZQ Treatments

We evaluated whether the lncRNAs found DE are shared between the three different PZQ-treatment regimens. Interestingly, most of the lncRNAs found DE in the “time-course experiment” (101 out of 189 lncRNAs, 54%, Figure 3), are also DE in the “4-weeks” or in the “7-weeks” experiment. The overlap of 54% of lncRNAs DE in the time-course experiment with lncRNAs DE in the other treatment regimens point to a limited set of affected lncRNAs, given the total of 16,583 lncRNAs present in the reference transcriptome [26]. A higher number of DE lncRNAs is shared between the “time-course experiment” and the “7-weeks experiment” (95 out of 189, 50%) when compared with the ones shared between the “time-course experiment” and the “4-weeks experiment” (26 out of 189, 14%), and this is probably because both “time-course” and the “7-weeks” experiments were performed with adult worms.
On the other hand, the “4-weeks experiment” shares more DE lncRNAs with the “7-weeks experiment” (62 out 181 lncRNAs) than with the “time-course experiment” (26 out 181 lncRNAs). In addition, 20 lncRNAs were found DE in all three experiments concomitantly, and we selected six of these for differential expression validation by RT-qPCR, as further explained below.
We also investigated whether the DE lncRNAs identified in our re-analysis overlapped with those modulated by 5-azacytidine [23], a drug known to affect S. mansoni fecundity [27]. Since 5-azacytidine is an epigenetic drug known to prevent DNA and RNA methylation [28,29], impairing the S. mansoni females’ transcription, translation and stem cell activities [30], we did not expect to find considerable overlap. Indeed, only a small proportion of DE lncRNAs from in vivo PZQ treatment experiments overlapped with those from 5-azacytidine (Supplementary Figure S5): 69 (8%) from the “4-weeks experiment”, 90 (10%) from the “7-weeks experiment”, and 59 (6%) from the “time-course experiment” out of 912 lncRNAs listed as DE after 5-azacytidine treatment.

2.3. Validation by RT-qPCR of LncRNAs Differential Expression

We selected ten candidate lncRNAs for RT-qPCR differential expression validation based on the differentially expressed genes (DEGs) identified simultaneously in the three experiments (six lncRNAs) or in the “time-course experiment” exclusively (four lncRNAs). These lncRNAs were selected based on the following criteria: (i) only lincRNAs (long intergenic non-coding RNAs) were selected, to avoid quantification of the levels of protein-coding genes expressed at the same genomic loci of antisense and sense lncRNAs and (ii) they show a wide range of expression levels in the RNA-Seq (trimmed mean of M values (TMM) from 0.10 to 172 in at least one of the conditions).
First, six protein-coding genes were used as controls for RT-qPCR differential expression validation. In the RT-qPCR assays, samples from adult worms collected after PZQ treatment of mice (400 mg/kg) at time points 0, 6, 12, 24, and 48 h post-treatment were used. Four out of the six tested protein-coding genes were validated: Smp_000660.1 (Ornithine aminotransferase), Smp_008660.2 (Severin), Smp_307450.3 (U-actitoxin-Avd3s) and Smp_311670.1 (U-actitoxin-Avd3s) showed statistically significant differential expression in at least two of the time points tested by RT-qPCR (Supplementary Figure S6, black lines). Despite the RT-qPCR results exhibiting the same expression kinetic profile observed in the RNA-Seq (Figure S6, blue lines) for Smp_132670.1 (Myosin regulatory light chain 2 smooth muscle) and Smp_307020.1 (Actin-2), these two genes did not show a statistically significant differential expression at any of the time points tested.
We then tested by RT-qPCR, in the same samples, ten lincRNAs identified as differentially expressed in the RNA-Seq datasets. Seven of the tested lincRNAs were confirmed by RT-qPCR to show a statistically significant differential expression in at least one of the time points (Figure 4, black lines): SmLINC101519-IBu, SmLINC105115-IBu, SmLINC110492-IBu, SmLINC121232-IBu, SmLINC161393-IBu, SmLINC163938-IBu, and SmLINC172840-IBu, with changes in expression ranging from 0.003 to 81.4-fold (Figure 4, black lines). The kinetic expression profile found in the RT-qPCR assays (Figure 4f, black lines) for SmLINC142502-IBu mirrored those from the RNA-Seq (Figure 4f, blue lines), however it did not show a statistically significant differential expression in the RT-qPCR assays. For SmLINC133371-IBu (Figure 4e) and SmLINC159037-IBu (Figure 4g), distinct expression profiles were found in the RT-qPCR assays when compared with the RNA-Seq.
Overall, differential expression patterns observed in the RT-qPCR data mirrored those obtained in the RNA-Seq re-analysis: high correlations between the fold changes found in the RT-qPCR assays and in the RNA-Seq data were found for the time points 12, 24 and 48 h (Supplementary Figure S7), with Pearson correlation coefficients of 0.828 (p-value = 6.85 × 10−13), 0.829 (p-value = 6.08 × 10−13) and 0.830 (p-value = 5.46 × 10−13), respectively. To ensure the validity of our assays, non-differentially expressed lincRNAs were also tested by RT-qPCR as an additional control. As expected, the five tested lincRNAs showed no statistically significant change in expression, relative to the untreated control, throughout the time points (Supplementary Figure S8).

2.4. LncRNAs Modulated by Praziquantel Are Differentially Expressed along S. mansoni Life-Cycle Stages

To evaluate whether the lncRNAs differentially expressed after PZQ treatment and tested by RT-qPCR are also expressed in other S. mansoni life-cycle stages, we looked at the “lncRNA transcriptome” [31] for the expression patterns of the ten selected lincRNAs in control assays. We observed a heterogeneous expression pattern distribution (Figure 5): expression of SmLINC105115-IBu (Figure 5b), SmLINC133371-IBu (Figure 5e) and SmLINC161393-IBu (Figure 5h) are enriched in miracidia and sporocysts stages, whereas SmLINC163938-IBu (Figure 5i) and SmLINC172840-IBu (Figure 5j) have higher expression levels in schistosomula.
While SmLINC121232-IBu shows higher expression levels in adult males (Figure 5d), SmLINC101519-IBu (Figure 5a), SmLINC142502-IBu (Figure 5f) and SmLINC159037-IBu (Figure 5g) are enriched in adult females, with SmLINC142502-IBu being exclusively expressed in this life-cycle stage. In turn, SmLINC110492-IBu (Figure 5c) is broadly expressed in cercariae, schistosomula, adult males and adult females. These results show that most of the tested lincRNAs are not specific to a single stage and may play roles in different S. mansoni life-cycle stages. As expected, most of the protein-coding genes were found not to be stage-specific, with only Smp_132670.1 and Smp_307020.1 showing enrichment in adult males (Supplementary Figure S9).

2.5. Weighted Gene Co-Expression Network Analysis Shows Terms Related to Drug Response Mechanisms

We performed a weighted gene co-expression network analysis (WGCNA) [32] on the RNA-Seq data from all three experiments together to identify co-expression networks integrating gene expression differences following in vivo PZQ treatment. After counts normalization and filtering, 14,040 transcripts remained (11,796 protein-coding genes and 2244 lncRNAs), and one sample (SRR10947815; “PZQ 96hr”), out of 30, was removed from the analysis. A total of 15 co-expression modules were identified (Supplementary Figure S10 and Supplementary Table S5).
We found that, except for SmLINC142502-IBu, all DE lncRNAs tested in RT-qPCR are part of four different modules, namely “brown2” (SmLINC163938-IBu), “floralwhite” (SmLINC101519-IBu, SmLINC105115-IBu, SmLINC110492-IBu, SmLINC133371-IBu and SmLINC159037-IBu), “plum” (SmLINC161393-IBu), and “plum2” (SmLINC121232-IBu and SmLINC172840-IBu). These lincRNAs show a high and significant membership to their corresponding modules (Supplementary Figure S11).
In addition to containing the lincRNAs tested by RT-qPCR, these four modules also stand out for other reasons. Two of them show the highest number of DEGs among all modules (Figure 6a), making them the most enriched in DEGs: in the “brown2” module, 244 out of 294 genes are DE (83%; p-value = 7.32 × 10−65, Fisher’s Exact Test), and in the “floralwhite” module, 584 out of 909 genes are DE (64%; p-value = 1.97 × 10−75, Fisher’s Exact Test). Although “plum” and “plum2” modules also contain a considerable amount of DEGs, only “plum” is significantly enriched, with 837 out of 2241 genes being DE (37%; p-value = 1.40 × 10−2). On the other hand, “plum” and “plum2” modules exhibit, respectively, the highest association with PZQ treatment (R2 = 0.53, p-value = 7.34 × 10−6; and R2 = 0.52, p-value = 1.13 × 10−5). As a result, we observe a clear separation of PZQ-treated samples and control group samples into two distinct clusters, distributed in the plot area of these modules’ eigengenes (the first principal component, PC1, of the expression matrix of the corresponding module) (Supplementary Figure S12).
Functional enrichment analysis of the protein-coding genes inside these four modules reveals that they are enriched with gene ontology (GO) terms that might indicate drug response mechanisms (Figure 6b,c and Supplementary Figure S13). For example, the “brown2” module is most significantly enriched in biological processes (BPs) such as “response to stimulus”, “signaling”, “response to stress”, and “cell surface receptor signaling pathway” (Figure 6b). The “floralwhite” module is also enriched in “response to stimulus”, “signaling” and “response to stress”; in addition, BP terms such as “transmembrane transport”, “monoatomic ion (transmembrane) transport”, “small molecule metabolic process”, and “DNA repair” are noteworthy (Figure 6c). Although in “plum” and “plum2” modules most of the enriched terms are related to metabolism, they also show enrichment among the terms as found in “brown2” and “floralwhite” modules (Supplementary Figure S13). Regarding molecular function GO terms, these four modules show significant enrichment in relevant terms such as “metal ion binding”, “calcium ion binding”, and “small molecule binding”.
Remarkably, 31 out of the 37 protein-coding genes highlighted by McCusker et al. [11] as involved in PZQ response are assigned to eight of the modules found here. Importantly, the majority of the highlighted Smps (25 Smps, not including isoforms) belong to the four most prominent modules: 12 Smps are in the “floralwhite” module, 11 Smps are in the “brown2” module, 2 Smps are in the “plum” module, and 4 Smps are in the “plum2” module.

2.6. PZQ-Modulated lncRNAs Co-Localize with PZQ-Modulated Smps in Different Cell Types

To identify the cell types in which the RT-qPCR validated differentially expressed lncRNAs under PZQ treatment are mainly expressed, we searched their expression localization on the Schistosoma mansoni single-cell cluster atlas [33]. Then, we matched these results with protein-coding genes modulated by PZQ that had been highlighted by McCusker et al. [11], which are expressed in the same single-cell clusters (Figure 7).
The most prominent match we identified was SmLINC121232-IBu (Figure 7a), which co-localizes with Smp_086480 (SmTAL2), Smp_086530 (SmTAL3), and Smp_169200 (SmTAL11), at the tegument cell cluster (Figure 7a, b). All these transcripts were found to be downregulated in our analysis and are highlighted DEGs in McCusker et al. [11] due to the immunomodulatory activity of the tegument allergen-like proteins (SmTALs) they encode [34,35,36].
Additionally, two other matching genes were found: Smp_076320 (Myb/SANT-like DNA-binding domain-containing protein 3) (Figure 7d); which is significantly enriched in vitellocytes, matching the cell cluster expression profile of SmLINC110492-IBu (Figure 7c); and the Smp_105220 (Lymphocyte antigen 6B, or SmLy6B) (Figure 7f) which is significantly enriched in neurons and muscles cell clusters, matching the profile of SmLINC101519-IBu (Figure 7e).

3. Discussion

Here, we have shown that S. mansoni long non-coding RNA levels can be modulated in vivo by praziquantel (PZQ), the only drug currently used to treat schistosomiasis. Hundreds of lncRNAs were found to be differentially expressed (DE) in S. mansoni after PZQ exposure of mice infected for 4 or 7 weeks. This shows that there is an extensive reprogramming of the gene expression in S. mansoni worms that goes beyond the changes in protein-coding gene expression. Indeed, this observation is in line with previous studies showing that lncRNA levels can be modulated by drugs in S. mansoni [23] and in other eukaryotes [25,37,38,39,40].
Several studies have demonstrated and validated the modulation of lncRNA levels in other eukaryotes as a response to different drugs [25,37,38,39,40]. According to the recently published “ncRNADrug” database, a manually curated resource containing thousands of ncRNAs in current literature found to be associated with drug response and resistance [37], a set of 67 unique lncRNAs has been validated as DE under the regulation of 74 different molecules, whilst regulating the expression of 127 unique genes.
In this regard, we sought to investigate possible roles of the lncRNAs that we validated as DE, and we observed that genes found to be modulated in vivo by PZQ are grouped into 15 co-expression modules, four of which contain the selected DE lncRNAs and most of the DE protein-coding genes (Smps) found by McCusker et al. [11]. Being co-expressed in the same modules suggests that these lncRNAs are potential regulators of the protein-coding genes whose expression levels are also modulated by PZQ.
We also observed that some of the PZQ-modulated Smps co-localize at the same single-cell clusters with the RT-qPCR-validated DE lncRNAs. For instance, Smp_086480 (SmTAL2), a member of the Schistosoma mansoni tegument allergen-like protein family (SmTAL) co-localizes with SmLINC121232-IBu. It is described that the majority of known SmTALs function as calcium-ion binding proteins, an important step of calcium signaling in eukaryotic cells [34,35,36]. Previous studies have successfully demonstrated that, except for SmTAL3 and SmTAL5, the other 11 known SmTALs bind to calcium-ions [34,35]. Importantly, SmTAL1, 4, 5, and 8 were shown to interact with PZQ [35] and, therefore, these proteins may be involved in the worm’s response to PZQ along with SmLINC121232-IBu.
In addition, we observed that SmLINC110492-IBu co-localizes with the vitellaria-specific Smp_076320, a Myb/SANT-like DNA-binding domain-containing protein. These DNA-binding domains are associated with transcription factors, chromatin-remodeling proteins and other transcriptional regulation proteins [41,42,43]. Interestingly, an orthologue of this protein was recently identified in the trematode Fasciola hepatica as a candidate gene involved in drug resistance [44].
Finally, it is also worth mentioning that SmLINC101519-IBu co-localizes with the muscle- and neuron-enriched Smp_105220 (Lymphocyte antigen 6B, SmLy6B). It is believed that the Ly6 protein family is tegument surface specific and capable of inducing host immune responses [45,46], but it has also been shown that in adult worms their expression is up-regulated in intra-parasite tissues, including parenchyma and muscles [47]. In humans and other mammals these proteins are involved in cell proliferation, migration, cell–cell interaction, cell signaling, and in the immune response [47,48,49].
Interestingly, all of these genes are expressed in the four most prominent co-expression modules in terms of the proportion of differentially expressed genes, correlation to PZQ treatment, and the presence of the validated lncRNAs: “brown2”, “floralwhite”, “plum”, and “plum2”. On top of that, these four co-expression modules are enriched in the calcium-ion binding function, as well as other processes that may be related to PZQ response, such as transmembrane ion transport, response to stimulus and stress, small molecule metabolism, and DNA repair.
The fact that some of the validated differentially expressed lncRNAs co-localize with these protein-coding genes in the same cell cluster and co-expression modules, especially SmLINC101519-IBu and Smp_105220, which are both part of the “floralwhite” module, raises the hypothesis of a possible regulation of these protein-coding genes by the lncRNAs identified here, and regulation of these specific proteins may be coherent in response to PZQ treatment. The SmTAL products could be related to the classic well-studied PZQ response in S. mansoni [34,35]. Further, although there is uncertainty around the functions of Myb/SANT-like domain-containing protein and SmLy6B, one might speculate they could be associated to stress response due to PZQ treatment. This is in line with the vitellaria-specific Smp_076320 being down-regulated, and the Smp_105220 being up-regulated upon PZQ treatment.
The molecular mechanism of action of PZQ leading to worm death is still poorly understood, although the drug has been the main anthelmintic used to control schistosomiasis for over 40 years [5,6,8]. It has been shown that PZQ has deleterious effects on parasite musculature and tegument in vitro [50,51], however the effects of PZQ treatments in vivo were not clear until recently, given the short half-life of PZQ [52] and the importance of host immune system engagement for drug efficacy in animal models [53,54]. Recently, McCusker et al. [11] showed that worms harvested from mouse livers following sub-lethal PZQ treatment revealed drug-evoked changes in the expression of putative immunomodulatory and anticoagulant gene products, such as tegument-like allergens [35] and Kunitz-type protease inhibitors [55], which can reflect mechanisms of parasite immune-evasion in response to chemotherapy [11]. Since lncRNAs have been shown to act as regulators of the expression of protein-coding genes through different mechanisms [12,13], it is possible that some of the lncRNAs found here by re-analysis of these RNA-Seq datasets [11] as differentially expressed upon in vivo PZQ treatment may act by regulating the expression of these putative immunomodulatory and anticoagulant protein-coding genes [11].
Moreover, several lncRNAs have been implied in drug resistance mechanisms in other organisms, in particular human cancers, where their mechanisms of action in response to drug treatment are well studied. Several lncRNAs have been reported in association with increase of drug efflux and metabolism, inhibition of apoptosis, induced DNA repair, and protection from oxidative stress [16,56,57]. In this context, it is crucial to note that drug resistance in S. mansoni is a growing concern in public health, since the emergence of resistant strains has already been observed in laboratory and is particularly concerning in the field [5,6,8]. It is possible that similar lncRNA-associated drug resistance mechanisms described in other models are related to the S. mansoni lncRNAs found here to be regulated by PZQ.
There are several mechanisms by which lncRNA expression can be regulated, and that might underlie the drug induced changes in lncRNA expression. In general, similar to that of protein-coding genes, lncRNA expression regulation is determined by modulation of chromatin state and by post-transcriptional modifications [58,59]. The primary mechanism for the regulation of lncRNA expression is based on chromatin state. LncRNAs might be either down-regulated by DNA hypermethylation or up-regulated by histone deacetylation in promoter regions. They usually have low CpG promoters, which is indicative of DNA hypermethylation throughout evolutionary history, and are thus associated with fewer transcriptional activation histone marks (e.g., H3K4me3) resulting in lower expression levels. Additionally, many transcription factors acting upon lncRNAs are known to be shared with protein-coding genes, and finer regulation of promoter activity operates in association with microRNAs [58]. Post-transcriptional regulation impacting RNA stability is also important for their regulation. For example, there are more than 100 distinct modified nucleotides in lncRNAs [58]. It is also known that microRNA interaction with RNA-binding proteins (RBP) mediates decay of lncRNA [58,59], as well as enzymes that control decapping, deadenylation, and nucleolytic activity [59].
Unlike short RNAs [60,61,62,63,64,65,66], the mechanisms of regulation of lncRNAs are largely unknown in parasites. In addition, while short RNAs (especially microRNAs) have been more explored in various helminths [67,68,69,70], lncRNAs have received little attention, being identified by transcriptomic approaches only in a few helminths other than S. mansoni and in protozoans [19,71,72,73] or studied in a limited number of free-living nematodes [74,75]. Therefore, it is time to consider lncRNAs as regulated by or as possible drug targets also in parasitic diseases, especially because they show lower conservation in their primary sequences between species than protein-coding genes [12,13], which in principle would favor the development of therapeutical strategies with less adverse effects [20].
The choice of lncRNAs to be further validated as drug targets should be based on the appropriate selection of lncRNA candidates, as reviewed by Silveira et al. [19]. This selection should be guided by functional characterization of the lncRNA as well as by the demonstration of the lncRNA relevance to the parasite biology. We have recently shown that lncRNAs are key components intervening in S. mansoni adult worm homeostasis, being essential for adult worm pairing status and survival in the mammalian host [24]. Interestingly, two of the lncRNAs modulated by PZQ, found here, are also differentially expressed between paired and unpaired S. mansoni adult worms (SmLINC133371-IBu and SmLINC101519-IBu), highlighting their involvement in adult worm homeostasis.
In summary, this study adds a deeper layer on the understanding of the effects of PZQ in S. mansoni and sheds light on the relevance of looking at lncRNA regulation in response to drug treatment in parasites. These data are, to our knowledge, the first report of modulation of lncRNA levels by a drug in vivo in any helminth. The mechanisms involving lncRNAs and protein-coding gene expression networks, to be further explored, will provide insights into potential mechanisms for PZQ treatment failure or routes to anthelmintic drug resistance.

4. Materials and Methods

4.1. Re-Analyses of PZQ Treatment RNA-Seq Data

Public RNA-Seq data previously generated by McCusker et al. [11] from Schistosoma mansoni worms after in vivo PZQ treatment were downloaded from the SRA-NCBI database (PRJNA597909 and PRJNA602528). Libraries were previously generated using the TruSeq Stranded mRNA kit (Illumina, San Diego, CA, USA) and sequenced using the Illumina HiSeq 2500 system [11]. The datasets included both a time-course and a transcriptional response experiment involving adult and immature worms. For the time-course experiments, adult worms were harvested from mice (n = 4 per time point) at ten different time points (0, 0.25, 1, 3, 6, 9, 12, 24, 48, and 96 h) after the administration to mice of a single curative dose of PZQ (400 mg/kg) 7 weeks post-infection. The transcriptional response experiments were performed by harvesting worms 14 h after treating mice (n = 5 per group) with either vehicle control or a sub-lethal PZQ dose (100 mg/kg). These transcriptional response experiments were conducted on mice treated with PZQ at 4 weeks post-infection (yielding juvenile worms) or at 7 weeks post-infection (yielding adult worms).
As previously established by our group [23,24,31,76], the identification and differential expression analysis of lncRNAs and protein-coding genes were performed using the following pipeline: primarily, adapters and low-quality reads were filtered out using fastqc v 0.11.9 [77] and fastp v 0.20.0 [78]. The S. mansoni genome sequence v.7 (PRJEA36577) from the WormBase ParaSite resource (WBPS14), along with the lncRNA transcriptome identified and annotated by Maciel et al. [26], which comprises 16,583 lncRNA genes in addition to 12,693 protein-coding genes, were used as references. STAR v 2.7.3a [79] was used to align the reads, which were quantified using RSEM v 1.3.1 [80]. The counts were normalized with the TMM method using edgeR v 3.36.0 [81,82,83]. Two distinct approaches were considered for statistical analysis of differential expression: limma+voom [84] and svaseq+edgeR [81,85]. Only the genes with an FDR lower than 0.05 in both analyses were considered differentially expressed (DE). For the DE genes (DEGs) in the time-course experiment, fold-change was estimated by a generalized linear model.

4.2. Selection of LncRNAs for RT-qPCR Validation

We selected ten candidate lncRNAs for RT-qPCR differential expression validation based on the DEGs identified within all three experiments or in the time-course experiment exclusively. Six candidate protein-coding genes were also selected for RT-qPCR differential expression validation. Additionally, the five most stable non-DE protein-coding genes across all time points were chosen as candidate reference genes for the RT-qPCR assays, given the better suitability of condition-specific reference genes as normalizers over standard housekeeping genes [31,86]. NormFinder [87] and geNorm [88] software were employed, via the online tool RefFinder [89] (https://blooge.cn/RefFinder/, accessed on 9 January 2024), to determine the two most stable candidate reference genes (Smp_017650.1 and Smp_036130.1), which were used to normalize the Ct values of the targeted genes and calculate their relative expression (2−∆∆Ct).

4.3. RNA Extraction, Quantification, and Quality Assessment

To validate the findings from our RNA-Seq re-analyses, we used worms obtained in a time-course experiment, performed by McCusker et al. [11] to validate protein-coding genes DE by RT-qPCR assays [11]. This experiment employed the same treatment as previously described for the RNA-Seq experiments, and a similar set of time points after PZQ treatment (0, 6, 12, 24, and 48 h, with n = 6 biological replicates each).
Total RNA was extracted using a combination of TRIzol (Invitrogen, Carlsbad, CA, USA) and the RNeasy Micro Kit (Qiagen, Waltham, MA, USA) for purification, according to the manufacturer’s instructions. RNA samples were quantified using the Qubit RNA HS Assay Kit (Thermo Fisher Scientific) in a Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and their integrities were verified in a 2100 Bioanalyzer Instrument (Agilent Technologies, Santa Clara, CA, USA) with Agilent RNA 6000 Pico Kit (Agilent Technologies).

4.4. Primer Design, Reverse Transcription, and Quantitative PCR (RT-qPCR) Assays

All primer pairs were designed using Primer3 (https://bioinfo.ut.ee/primer3-0.4.0/, accessed on 9 January 2024), aiming for primer pairs designed between two exons whenever possible. Each primer was designed with the following parameters: maximum primer size of 25 nt, annealing temperature ranging from 57 °C to 63 °C, GC content between 30% and 60%, product size varying between 90 and 250 nt, and a maximum allowable length of a mononucleotide repeat set to three. The primer sequences and efficiencies are shown in Supplementary Table S6.
We used Superscript IV (Life Technologies, Carlsbad, CA, USA) for the reverse transcription, generating cDNA which was diluted eight times in water. The RT-qPCR assays were performed using 1X LightCycler 480 SYBR Green I Master Mix (Roche Diagnostics, Basel, Switzerland) with 800 nM of each oligonucleotide pair on a LightCycler 480 System instrument (Roche Diagnostics). All primers were tested in technical triplicates for all samples, in addition to a no template control (NTC).
To calculate differential gene expression, the delta-Ct method was applied [90]. Raw Ct values are shown in Supplementary Table S7. Differences in the expression of the assessed genes across all five time points were evaluated by a linear mixed-effects model [91]. This type of model accounts for the random effects observed in our samples, leading to a more accurate fit [91]. The model fitting and determination of significance (at a 0.05 significance level) were performed with the R package nlme v 3.1-164 [92].

4.5. Gene Expression Patterns of Selected Genes across S. mansoni Life-Cycle Stages

We generated plots showcasing the expression patterns for the selected lncRNAs and protein-coding genes across S. mansoni life-cycle stages (miracidium, sporocyst, cercaria, schistosomulum, adult male, and adult female). These plots were made using publicly available RNA-Seq normalized counts (TMM), from Silveira et al. [31], available at https://verjolab.shinyapps.io/Reference-genes/ (accessed on 9 January 2024).

4.6. Weighted Gene Co-Expression Network Analysis (WGCNA)

We performed a WGCNA on the RNA-Seq data from all the three in vivo PZQ experiments together, using the R package WGCNA v1.72.5 [32]. First, the transcripts count for the “time-course experiment”, “4-weeks experiment”, and “7-weeks experiment” were normalized applying a variance stabilizing transformation with DESeq2 v1.38.3 [93]. We considered the genes expressed only if their counts were greater than or equal to three in at least 75% of the samples. The sample SRR10947815 was removed from this analysis, because of its outlier behavior observed in the hierarchical clustering, which isolated it from any sample cluster.
To maximize scale-free topology properties, fit index was calculated as a function of a series of soft-threshold powers, in order to pick the best power (Figure S10a). Based on fit index (R2) and mean connectivity (k), power equal to 9 (R2 = 0.973, k = 86.36) was selected for network construction and module detection. Modules were constructed based on topological overlap matrix (TOM) dissimilarity (1 − TOM). The closest modules were then merged based on their eigengenes values (the first principal component, PC1, of the expression matrix of the corresponding module) dissimilarity distance, considering a 0.40 dissimilarity threshold (Supplementary Figure S10b).
Intramodular analysis was carried out on the merged modules to identify genes’ membership, i.e., the Pearson’s correlation between gene expression data and their corresponding module eigengene (Supplementary Figure S11). To assess whether the modules were enriched in DEGs, a one-tailed Fisher’s exact test was conducted, with the alternative hypothesis stating that the observed number of DEGs in a module is greater than expected. Furthermore, each module eigengene value was investigated with sample group (PZQ-treated or control) in a linear model (eigengene ~ sample group) in order to evaluate the association between modules and PZQ treatment.
Functional profiling of the protein-coding genes inside the modules was performed with gProfiler2 v0.2.3 [94] with g:SCS multiple testing correction method. Gene ontology (GO) terms annotation of S. mansoni protein-coding genes (Smp) were obtained using gProfiler’s “Smansoni_v7” database. REVIGO v.1.8.1 [95] was used to remove redundancy in enriched GO terms (Figure 6 and Supplementary Figure S13).

4.7. Single-Cell Clusters Search of LncRNAs and Protein-Coding Genes

The S. mansoni single-cell clusters atlas publicly available at Schistosoma mansoni LncRNAs Database (http://verjolab.usp.br:8081/ accessed on 6 March 2024) [33] was used to identify the cell type where each validated DE lncRNA, and protein-coding gene highlighted by McCusker et al. [11], is expressed and significantly enriched. Results of matching cell types between the lncRNA and the protein-coding genes were selected.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ncrna10020027/s1; Table S1: RNA-Seq samples analyzed and alignment statistics; Table S2: List of transcripts differentially expressed from S. mansoni adult worm couples obtained across 10 timepoints of 0, 0.25, 1, 3, 6, 9, 12, 24, 48 and 96 h from mice after treatment with a single, curative dose of PZQ (400 mg/kg) delivered by oral gavage at 7 weeks post-infection (PRJNA602528); Table S3: List of transcripts differentially expressed from S. mansoni adult worm couples obtained from mice 14 hours after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 4 weeks post-infection (PRJNA597909); Table S4: List of transcripts differentially expressed from S. mansoni adult worm couples obtained from mice 14 hours after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 7 weeks post-infection (PRJNA597909); Table S5: List of genes belonging to each module according to the WGCNA analysis; Table S6: List of primers used in RT-qPCR; Table S7: Cts (Cycle Threshold) of the lncRNAs and Smps obtained in the RT-qPCR assays; Figure S1: Clustering of RNA-Seq biological replicates assessed by principal 12 component analysis (PCA). Figure S2: Heatmap of differentially expressed protein-coding genes (Smps) detected by RNA-Seq of adult worm couples obtained from mice across 10 timepoints at 0, 0.25, 1, 3, 6, 9, 12, 24, 48 and 96 h after treatment (Treat) with a single, curative dose of PZQ (400 mg/kg) delivered by oral gavage at 7 weeks post-infection; Figure S3: Heatmaps of differentially expressed protein-coding genes (Smps) detected by RNA-Seq of adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 4 weeks post-infection; Figure S4: Heatmaps of differentially expressed protein-coding genes (Smps) detected by RNA-Seq of adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 7 weeks post-infection; Figure S5: Analysis of long non-coding RNAs (lncRNAs) differentially expressed upon distinct praziquantel (PZQ) treatment regimens and 5-azacytidine (5-AzaC) treatment; Figure S6: Validation by RT-qPCR of protein-coding genes (Smps) detected as differentially expressed upon praziquantel treatment in vivo; Figure S7: Correlation between RNA-Seq and RT-qPCR analyses; Figure S8: RT-qPCR results of long non-coding RNAs (lncRNAs) not detected as differentially expressed upon praziquantel treatment in vivo; Figure S9: RNA-Seq expression profiles of selected protein-coding genes (Smps) measured in control assays at different S. mansoni life-cycle stages; Figure S10: Weighted Gene Co-Expression Network Analysis (WGCNA) parameters and results for the RNA-Seq experiments data; Figure S11: Gene-Module correlation for the four modules containing the differentially expressed long non-coding RNAs (lncRNAs) tested by RT-qPCR; Figure S12: Clustering of RNA-Seq biological replicates assessed by modules eigengenes values; Figure S13: Functional profiling of the protein-coding genes inside “plum” and “plum2” modules.

Author Contributions

Conceptualization: M.S.A.; methodology: P.J.P. and M.S.A.; validation: P.J.P., A.F.-C. and M.S.A.; formal analysis: P.J.P., A.C.T. and M.S.A.; investigation: P.J.P., A.C.T., J.D.C., S.V.-A. and M.S.A.; resources: J.D.C., S.V.-A. and M.S.A.; data curation: P.J.P., A.C.T. and M.S.A.; writing—original draft preparation: P.J.P. and M.S.A.; writing—review and editing,: P.J.P., S.V.-A. and M.S.A.; visualization: P.J.P. and A.C.T.; supervision: M.S.A.; funding acquisition: S.V.-A. and M.S.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by a grant from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) Thematic grant number 2018/23693-5 to SV-A. PJP and AF-C received fellowships from Conselho Nacional de Desenvolvimento Científico e Tecnológico (Fundação Butantan/CNPq and 148745/2022-9, respectively). SV-A received an established investigator fellowship award from CNPq (306646/2019-6), Brazil.

Institutional Review Board Statement

The animal study protocol was approved by the Ethics Committee of the Medical College of Wisconsin (protocol codes AUA00006471 and AUA00006735). Animal work was carried out with the oversight and approval of the Laboratory Animal Resources facility at the Medical College of Wisconsin, adhering to the humane standards for the health and welfare of animals used for biomedical purposes defined by the Animal Welfare Act and the Health Research Extension Act.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data generated or analyzed during this study are included in this published article (and its Supplementary Information files).

Acknowledgments

We would like to thank Ana Paula Lopes Vidal for the administrative support.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in: the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

  1. Colley, D.G.; Bustinduy, A.L.; Secor, W.E.; King, C.H. Human schistosomiasis. Lancet 2014, 383, 2253–2264. [Google Scholar] [CrossRef] [PubMed]
  2. WHO. Fact Sheets: Schistosomiasis; 2023. Available online: https://www.who.int/news-room/fact-sheets/detail/schistosomiasis (accessed on 9 January 2024).
  3. Zacharia, A.; Mushi, V.; Makene, T. A systematic review and meta-analysis on the rate of human schistosomiasis reinfection. PLoS ONE 2020, 15, e0243224. [Google Scholar] [CrossRef] [PubMed]
  4. WHO. Global Health Estimates 2020: Disease Burden by Cause, Age, Sex, by Country and by Region, 2000–2019. Geneva, World Health Organization. 2020. Available online: https://www.who.int/data/gho/data/themes/mortality-and-global-health-estimates/global-health-estimates-leading-causes-of-dalys (accessed on 9 January 2024).
  5. McManus, D.P.; Dunne, D.W.; Sacko, M.; Utzinger, J.; Vennervald, B.J.; Zhou, X.N. Schistosomiasis. Nat. Rev. Dis. Primers 2018, 4, 13. [Google Scholar] [CrossRef] [PubMed]
  6. Vale, N.; Gouveia, M.J.; Rinaldi, G.; Brindley, P.J.; Gärtner, F.; Correia da Costa, J.M. Praziquantel for Schistosomiasis: Single-Drug Metabolism Revisited, Mode of Action, and Resistance. Antimicrob. Agents Chemother. 2017, 61, e02582-16. [Google Scholar] [CrossRef]
  7. Lewis, F.A.; Tucker, M.S. Schistosomiasis. Adv. Exp. Med. Biol. 2014, 766, 47–75. [Google Scholar] [CrossRef] [PubMed]
  8. Wilson, R.A. Schistosomiasis then and now: What has changed in the last 100 years? Parasitology 2020, 147, 507–515. [Google Scholar] [CrossRef]
  9. Park, S.K.; Friedrich, L.; Yahya, N.A.; Rohr, C.M.; Chulkov, E.G.; Maillard, D.; Rippmann, F.; Spangenberg, T.; Marchant, J.S. Mechanism of praziquantel action at a parasitic flatworm ion channel. Sci. Transl. Med. 2021, 13, eabj5832. [Google Scholar] [CrossRef] [PubMed]
  10. Le Clec’h, W.; Chevalier, F.D.; Mattos, A.C.A.; Strickland, A.; Diaz, R.; McDew-White, M.; Rohr, C.M.; Kinung’hi, S.; Allan, F.; Webster, B.L.; et al. Genetic analysis of praziquantel response in schistosome parasites implicates a transient receptor potential channel. Sci. Transl. Med. 2021, 13, eabj9114. [Google Scholar] [CrossRef]
  11. McCusker, P.; Rohr, C.M.; Chan, J.D. Schistosoma mansoni alter transcription of immunomodulatory gene products following in vivo praziquantel exposure. PLoS Negl. Trop. Dis. 2021, 15, e0009200. [Google Scholar] [CrossRef]
  12. Quinn, J.J.; Chang, H.Y. Unique features of long non-coding RNA biogenesis and function. Nat. Rev. Genet. 2016, 17, 47–62. [Google Scholar] [CrossRef]
  13. Ransohoff, J.D.; Wei, Y.; Khavari, P.A. The functions and unique features of long intergenic non-coding RNA. Nat. Rev. Mol. Cell Biol. 2018, 19, 143–157. [Google Scholar] [CrossRef]
  14. Statello, L.; Guo, C.J.; Chen, L.L.; Huarte, M. Gene regulation by long non-coding RNAs and its biological functions. Nat. Rev. Mol. Cell Biol. 2021, 22, 96–118. [Google Scholar] [CrossRef]
  15. Kopp, F.; Mendell, J.T. Functional Classification and Experimental Dissection of Long Noncoding RNAs. Cell 2018, 172, 393–407. [Google Scholar] [CrossRef] [PubMed]
  16. Liu, K.; Gao, L.; Ma, X.; Huang, J.J.; Chen, J.; Zeng, L.; Ashby, C.R., Jr.; Zou, C.; Chen, Z.S. Long non-coding RNAs regulate drug resistance in cancer. Mol. Cancer 2020, 19, 54. [Google Scholar] [CrossRef] [PubMed]
  17. Smallegan, M.J.; Rinn, J.L. Linking long noncoding RNA to drug resistance. Proc. Natl. Acad. Sci. USA 2019, 116, 21963–21965. [Google Scholar] [CrossRef]
  18. Melé, M.; Mattioli, K.; Mallard, W.; Shechner, D.M.; Gerhardinger, C.; Rinn, J.L. Chromatin environment, transcriptional regulation, and splicing distinguish lincRNAs and mRNAs. Genome Res. 2017, 27, 27–37. [Google Scholar] [CrossRef]
  19. Silveira, G.O.; Coelho, H.S.; Amaral, M.S.; Verjovski-Almeida, S. Long non-coding RNAs as possible therapeutic targets in protozoa, and in Schistosoma and other helminths. Parasitol Res. 2022, 121, 1091–1115. [Google Scholar] [CrossRef]
  20. Winkle, M.; El-Daly, S.M.; Fabbri, M.; Calin, G.A. Noncoding RNA therapeutics-challenges and potential solutions. Nat. Rev. Drug Discov. 2021, 20, 629–651. [Google Scholar] [CrossRef] [PubMed]
  21. Oliveira, K.C.; Carvalho, M.L.; Maracaja-Coutinho, V.; Kitajima, J.P.; Verjovski-Almeida, S. Non-coding RNAs in schistosomes: An unexplored world. An. Acad. Bras. Cienc. 2011, 83, 673–694. [Google Scholar] [CrossRef]
  22. Kim, H.C.; Khalil, A.M.; Jolly, E.R. LncRNAs in molluscan and mammalian stages of parasitic schistosomes are developmentally-regulated and coordinately expressed with protein-coding genes. RNA Biol. 2020, 17, 805–815. [Google Scholar] [CrossRef]
  23. Amaral, M.S.; Maciel, L.F.; Silveira, G.O.; Olberg, G.G.O.; Leite, J.V.P.; Imamura, L.K.; Pereira, A.S.A.; Miyasato, P.A.; Nakano, E.; Verjovski-Almeida, S. Long non-coding RNA levels can be modulated by 5-azacytidine in Schistosoma mansoni. Sci. Rep. 2020, 10, 21565. [Google Scholar] [CrossRef]
  24. Silveira, G.O.; Coelho, H.S.; Pereira, A.S.A.; Miyasato, P.A.; Santos, D.W.; Maciel, L.F.; Olberg, G.G.G.; Tahira, A.C.; Nakano, E.; Oliveira, M.L.S.; et al. Long non-coding RNAs are essential for Schistosoma mansoni pairing-dependent adult worm homeostasis and fertility. PLoS Pathog. 2023, 19, e1011369. [Google Scholar] [CrossRef] [PubMed]
  25. Jiang, W.; Qu, Y.; Yang, Q.; Ma, X.; Meng, Q.; Xu, J.; Liu, X.; Wang, S. D-lnc: A comprehensive database and analytical platform to dissect the modification of drugs on lncRNA expression. RNA Biol. 2019, 16, 1586–1591. [Google Scholar] [CrossRef]
  26. Maciel, L.F.; Morales-Vicente, D.A.; Silveira, G.O.; Ribeiro, R.O.; Olberg, G.G.O.; Pires, D.S.; Amaral, M.S.; Verjovski-Almeida, S. Weighted Gene Co-Expression Analyses Point to Long Non-Coding RNA Hub Genes at Different Schistosoma mansoni Life-Cycle Stages. Front. Genet. 2019, 10, 823. [Google Scholar] [CrossRef]
  27. Geyer, K.K.; Rodriguez Lopez, C.M.; Chalmers, I.W.; Munshi, S.E.; Truscott, M.; Heald, J.; Wilkinson, M.J.; Hoffmann, K.F. Cytosine methylation regulates oviposition in the pathogenic blood fluke Schistosoma mansoni. Nat. Commun. 2011, 2, 424. [Google Scholar] [CrossRef]
  28. Lu, L.J.; Randerath, K. Mechanism of 5-azacytidine-induced transfer RNA cytosine-5-methyltransferase deficiency. Cancer Res. 1980, 40, 2701–2705. [Google Scholar] [PubMed]
  29. Taylor, S.M.; Jones, P.A. Mechanism of action of eukaryotic DNA methyltransferase. Use of 5-azacytosine-containing DNA. J. Mol. Biol. 1982, 162, 679–692. [Google Scholar] [CrossRef] [PubMed]
  30. Geyer, K.K.; Munshi, S.E.; Vickers, M.; Squance, M.; Wilkinson, T.J.; Berrar, D.; Chaparro, C.; Swain, M.T.; Hoffmann, K.F. The anti-fecundity effect of 5-azacytidine (5-AzaC) on Schistosoma mansoni is linked to dis-regulated transcription, translation and stem cell activities. Int. J. Parasitol. Drugs Drug Resist. 2018, 8, 213–222. [Google Scholar] [CrossRef]
  31. Silveira, G.O.; Amaral, M.S.; Coelho, H.S.; Maciel, L.F.; Pereira, A.S.A.; Olberg, G.G.O.; Miyasato, P.A.; Nakano, E.; Verjovski-Almeida, S. Assessment of reference genes at six different developmental stages of Schistosoma mansoni for quantitative RT-PCR. Sci. Rep. 2021, 11, 16816. [Google Scholar] [CrossRef]
  32. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 2008, 9, 559. [Google Scholar] [CrossRef]
  33. Morales-Vicente, D.A.; Zhao, L.; Silveira, G.O.; Tahira, A.C.; Amaral, M.S.; Collins, J.J., 3rd; Verjovski-Almeida, S. Single-cell RNA-seq analyses show that long non-coding RNAs are conspicuously expressed in Schistosoma mansoni gamete and tegument progenitor cell populations. Front. Genet. 2022, 13, 924877. [Google Scholar] [CrossRef] [PubMed]
  34. Thomas, C.M.; Fitzsimmons, C.M.; Dunne, D.W.; Timson, D.J. Comparative biochemical analysis of three members of the Schistosoma mansoni TAL family: Differences in ion and drug binding properties. Biochimie 2015, 108, 40–47. [Google Scholar] [CrossRef] [PubMed]
  35. Carson, J.; Thomas, C.M.; McGinty, A.; Takata, G.; Timson, D.J. The tegumental allergen-like proteins of Schistosoma mansoni: A biochemical study of SmTAL4-TAL13. Mol. Biochem. Parasitol. 2018, 221, 14–22. [Google Scholar] [CrossRef] [PubMed]
  36. Thomas, C.M.; Timson, D.J. A mysterious family of calcium-binding proteins from parasitic worms. Biochem. Soc. Trans. 2016, 44, 1005–1010. [Google Scholar] [CrossRef] [PubMed]
  37. Cao, X.; Zhou, X.; Hou, F.; Huang, Y.E.; Yuan, M.; Long, M.; Chen, S.; Lei, W.; Zhu, J.; Chen, J.; et al. ncRNADrug: A database for validated and predicted ncRNAs associated with drug resistance and targeted by drugs. Nucleic Acids Res. 2024, 52, D1393–D1399. [Google Scholar] [CrossRef] [PubMed]
  38. Ghiam, S.; Eslahchi, C.; Shahpasand, K.; Habibi-Rezaei, M.; Gharaghani, S. Identification of repurposed drugs targeting significant long non-coding RNAs in the cross-talk between diabetes mellitus and Alzheimer’s disease. Sci. Rep. 2022, 12, 18332. [Google Scholar] [CrossRef] [PubMed]
  39. Chen, H.; Zhang, Z.; Zhang, J. In silico drug repositioning based on integrated drug targets and canonical correlation analysis. BMC Med. Genomics 2022, 15, 48. [Google Scholar] [CrossRef]
  40. Guo, H.; Liu, J.; Ben, Q.; Qu, Y.; Li, M.; Wang, Y.; Chen, W.; Zhang, J. The aspirin-induced long non-coding RNA OLA1P2 blocks phosphorylated STAT3 homodimer formation. Genome Biol. 2016, 17, 24. [Google Scholar] [CrossRef] [PubMed]
  41. Klempnauer, K.H.; Sippel, A.E. The highly conserved amino-terminal region of the protein encoded by the v-myb oncogene functions as a DNA-binding domain. EMBO J. 1987, 6, 2719–2725. [Google Scholar] [CrossRef]
  42. Biedenkapp, H.; Borgmeyer, U.; Sippel, A.E.; Klempnauer, K.H. Viral myb oncogene encodes a sequence-specific DNA-binding activity. Nature 1988, 335, 835–837. [Google Scholar] [CrossRef]
  43. Boyer, L.A.; Latek, R.R.; Peterson, C.L. The SANT domain: A unique histone-tail-binding module? Nat. Rev. Mol. Cell Biol. 2004, 5, 158–163. [Google Scholar] [CrossRef] [PubMed]
  44. Beesley, N.J.; Cwiklinski, K.; Allen, K.; Hoyle, R.C.; Spithill, T.W.; La Course, E.J.; Williams, D.J.L.; Paterson, S.; Hodgkinson, J.E. A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance. PLoS Pathog. 2023, 19, e1011081. [Google Scholar] [CrossRef] [PubMed]
  45. Chalmers, I.W.; Fitzsimmons, C.M.; Brown, M.; Pierrot, C.; Jones, F.M.; Wawrzyniak, J.M.; Fernandez-Fuentes, N.; Tukahebwa, E.M.; Dunne, D.W.; Khalife, J.; et al. Human IgG1 Responses to Surface Localised Schistosoma mansoni Ly6 Family Members Drop following Praziquantel Treatment. PLoS Negl. Trop. Dis. 2015, 9, e0003920. [Google Scholar] [CrossRef]
  46. Egesa, M.; Lubyayi, L.; Tukahebwa, E.M.; Bagaya, B.S.; Chalmers, I.W.; Wilson, S.; Hokke, C.H.; Hoffmann, K.F.; Dunne, D.W.; Yazdanbakhsh, M.; et al. Schistosoma mansoni schistosomula antigens induce Th1/Pro-inflammatory cytokine responses. Parasite Immunol. 2018, 40, e12592. [Google Scholar] [CrossRef]
  47. Farias, L.P.; Krautz-Peterson, G.; Tararam, C.A.; Araujo-Montoya, B.O.; Fraga, T.R.; Rofatto, H.K.; Silva, F.P., Jr.; Isaac, L.; Da’dara, A.A.; Wilson, R.A.; et al. On the three-finger protein domain fold and CD59-like proteins in Schistosoma mansoni. PLoS Negl. Trop. Dis. 2013, 7, e2482. [Google Scholar] [CrossRef]
  48. Kong, H.K.; Park, J.H. Characterization and function of human Ly-6/uPAR molecules. BMB Rep. 2012, 45, 595–603. [Google Scholar] [CrossRef]
  49. Loughner, C.L.; Bruford, E.A.; McAndrews, M.S.; Delp, E.E.; Swamynathan, S.; Swamynathan, S.K. Organization, evolution and functions of the human and mouse Ly6/uPAR family genes. Hum. Genomics 2016, 10, 10. [Google Scholar] [CrossRef]
  50. Fetterer, R.H.; Pax, R.A.; Strand, S.; Bennett, J.L. Schistosoma mansoni: Physical and chemical factors affecting the mechanical properties of the adult male musculature. Exp. Parasitol. 1978, 46, 59–71. [Google Scholar] [CrossRef] [PubMed]
  51. Xiao, S.H.; Catto, B.A.; Webster, L.T., Jr. Effects of praziquantel on different developmental stages of Schistosoma mansoni in vitro and in vivo. J. Infect. Dis. 1985, 151, 1130–1137. [Google Scholar] [CrossRef]
  52. Olliaro, P.; Delgado-Romero, P.; Keiser, J. The little we know about the pharmacokinetics and pharmacodynamics of praziquantel (racemate and R-enantiomer). J. Antimicrob. Chemother. 2014, 69, 863–870. [Google Scholar] [CrossRef]
  53. Panic, G.; Ruf, M.T.; Keiser, J. Immunohistochemical Investigations of Treatment with Ro 13-3978, Praziquantel, Oxamniquine, and Mefloquine in Schistosoma mansoni-Infected Mice. Antimicrob. Agents. Chemother. 2017, 61. [Google Scholar] [CrossRef] [PubMed]
  54. Reimers, N.; Homann, A.; Hoschler, B.; Langhans, K.; Wilson, R.A.; Pierrot, C.; Khalife, J.; Grevelding, C.G.; Chalmers, I.W.; Yazdanbakhsh, M.; et al. Drug-induced exposure of Schistosoma mansoni antigens SmCD59a and SmKK7. PLoS Negl. Trop. Dis. 2015, 9, e0003593. [Google Scholar] [CrossRef] [PubMed]
  55. Morais, S.B.; Figueiredo, B.C.; Assis, N.R.G.; Alvarenga, D.M.; de Magalhaes, M.T.Q.; Ferreira, R.S.; Vieira, A.T.; Menezes, G.B.; Oliveira, S.C. Schistosoma mansoni SmKI-1 serine protease inhibitor binds to elastase and impairs neutrophil function and inflammation. PLoS Pathog. 2018, 14, e1006870. [Google Scholar] [CrossRef]
  56. Wang, Y.; Qin, Z.; Cai, S.; Yu, L.; Hu, H.; Zeng, S. The role of non-coding RNAs in ABC transporters regulation and their clinical implications of multidrug resistance in cancer. Expert. Opin. Drug Metab. Toxicol. 2021, 17, 291–306. [Google Scholar] [CrossRef]
  57. He, J.; Zhu, S.; Liang, X.; Zhang, Q.; Luo, X.; Liu, C.; Song, L. LncRNA as a multifunctional regulator in cancer multi-drug resistance. Mol. Biol. Rep. 2021, 48, 1–15. [Google Scholar] [CrossRef] [PubMed]
  58. Wu, Z.; Liu, X.; Liu, L.; Deng, H.; Zhang, J.; Xu, Q.; Cen, B.; Ji, A. Regulation of lncRNA expression. Cell Mol. Biol. Lett. 2014, 19, 561–575. [Google Scholar] [CrossRef] [PubMed]
  59. Yoon, J.H.; Kim, J.; Gorospe, M. Long noncoding RNA turnover. Biochimie 2015, 117, 15–21. [Google Scholar] [CrossRef]
  60. Zheng, Y.; Cai, X.; Bradley, J.E. microRNAs in parasites and parasite infection. RNA Biol. 2013, 10, 371–379. [Google Scholar] [CrossRef]
  61. Britton, C.; Winter, A.D.; Gillan, V.; Devaney, E. microRNAs of parasitic helminths-Identification, characterization and potential as drug targets. Int. J. Parasitol. Drugs Drug. Resist. 2014, 4, 85–94. [Google Scholar] [CrossRef]
  62. Tritten, L.; Burkman, E.; Moorhead, A.; Satti, M.; Geary, J.; Mackenzie, C.; Geary, T. Detection of circulating parasite-derived microRNAs in filarial infections. PLoS Negl. Trop. Dis. 2014, 8, e2971. [Google Scholar] [CrossRef]
  63. Marks, N.D.; Winter, A.D.; Gu, H.Y.; Maitland, K.; Gillan, V.; Ambroz, M.; Martinelli, A.; Laing, R.; MacLellan, R.; Towne, J.; et al. Profiling microRNAs through development of the parasitic nematode Haemonchus identifies nematode-specific miRNAs that suppress larval development. Sci. Rep. 2019, 9, 17594. [Google Scholar] [CrossRef]
  64. Meningher, T.; Barsheshet, Y.; Ofir-Birin, Y.; Gold, D.; Brant, B.; Dekel, E.; Sidi, Y.; Schwartz, E.; Regev-Rudzki, N.; Avni, O.; et al. Schistosomal extracellular vesicle-enclosed miRNAs modulate host T helper cell differentiation. EMBO Rep. 2020, 21, e47882. [Google Scholar] [CrossRef] [PubMed]
  65. Liu, J.; Zhu, L.; Wang, J.; Qiu, L.; Chen, Y.; Davis, R.E.; Cheng, G. Schistosoma japonicum extracellular vesicle miRNA cargo regulates host macrophage functions facilitating parasitism. PLoS Pathog. 2019, 15, e1007817. [Google Scholar] [CrossRef] [PubMed]
  66. Leija-Montoya, A.G.; Gonzalez-Ramirez, J.; Martinez-Coronilla, G.; Mejia-Leon, M.E.; Isiordia-Espinoza, M.; Sanchez-Munoz, F.; Chavez-Cortez, E.G.; Pitones-Rubio, V.; Serafin-Higuera, N. Roles of microRNAs and Long Non-Coding RNAs Encoded by Parasitic Helminths in Human Carcinogenesis. Int J. Mol. Sci 2022, 23, 8173. [Google Scholar] [CrossRef] [PubMed]
  67. Shao, C.C.; Xu, M.J.; Alasaad, S.; Song, H.Q.; Peng, L.; Tao, J.P.; Zhu, X.Q. Comparative analysis of microRNA profiles between adult Ascaris lumbricoides and Ascaris suum. BMC Vet. Res. 2014, 10, 99. [Google Scholar] [CrossRef]
  68. Fontenla, S.; Rinaldi, G.; Smircich, P.; Tort, J.F. Conservation and diversification of small RNA pathways within flatworms. BMC Evol. Biol. 2017, 17, 215. [Google Scholar] [CrossRef]
  69. Macchiaroli, N.; Cucher, M.; Kamenetzky, L.; Yones, C.; Bugnon, L.; Berriman, M.; Olson, P.D.; Rosenzvit, M.C. Identification and expression profiling of microRNAs in Hymenolepis. Int. J. Parasitol. 2019, 49, 211–223. [Google Scholar] [CrossRef]
  70. Holz, A.; Streit, A. Gain and Loss of Small RNA Classes-Characterization of Small RNAs in the Parasitic Nematode Family Strongyloididae. Genome Biol. Evol. 2017, 9, 2826–2843. [Google Scholar] [CrossRef]
  71. Santos, L.N.; Silva, E.S.; Santos, A.S.; De Sa, P.H.; Ramos, R.T.; Silva, A.; Cooper, P.J.; Barreto, M.L.; Loureiro, S.; Pinheiro, C.S.; et al. De novo assembly and characterization of the Trichuris trichiura adult worm transcriptome using Ion Torrent sequencing. Acta Tropica. 2016, 159, 132–141. [Google Scholar] [CrossRef]
  72. Azlan, A.; Halim, M.A.; Azzam, G. Genome-wide identification and characterization of long intergenic noncoding RNAs in the regenerative flatworm Macrostomum lignano. Genomics 2020, 112, 1273–1281. [Google Scholar] [CrossRef]
  73. Ross, E.; Blair, D.; Guerrero-Hernandez, C.; Sanchez Alvarado, A. Comparative and Transcriptome Analyses Uncover Key Aspects of Coding- and Long Noncoding RNAs in Flatworm Mitochondrial Genomes. G3 2016, 6, 1191–1200. [Google Scholar] [CrossRef] [PubMed]
  74. Rodelsperger, C.; Menden, K.; Serobyan, V.; Witte, H.; Baskaran, P. First insights into the nature and evolution of antisense transcription in nematodes. BMC Evol. Biol. 2016, 16, 165. [Google Scholar] [CrossRef] [PubMed]
  75. Wei, S.; Chen, H.; Dzakah, E.E.; Yu, B.; Wang, X.; Fu, T.; Li, J.; Liu, L.; Fang, S.; Liu, W.; et al. Systematic evaluation of C. elegans lincRNAs with CRISPR knockout mutants. Genome Biol. 2019, 20, 7. [Google Scholar] [CrossRef] [PubMed]
  76. Amaral, M.S.; Santos, D.W.; Pereira, A.S.A.; Tahira, A.C.; Malvezzi, J.V.M.; Miyasato, P.A.; Freitas, R.P.; Kalil, J.; Tjon Kon Fat, E.M.; de Dood, C.J.; et al. Rhesus macaques self-curing from a schistosome infection can display complete immunity to challenge. Nat. Commun. 2021, 12, 6181. [Google Scholar] [CrossRef] [PubMed]
  77. Andrews, S. FastQC: A Quality Control Analysis Tool for High Throughput Sequencing Data, 0.11.9; 2010. Available online: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. (accessed on 1 December 2023).
  78. Chen, S.; Zhou, Y.; Chen, Y.; Gu, J. fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 2018, 34, i884–i890. [Google Scholar] [CrossRef]
  79. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29, 15–21. [Google Scholar] [CrossRef] [PubMed]
  80. Li, B.; Dewey, C.N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, 323. [Google Scholar] [CrossRef] [PubMed]
  81. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef] [PubMed]
  82. Chen, Y.; Lun, A.T.; Smyth, G.K. From reads to genes to pathways: Differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Res 2016, 5, 1438. [Google Scholar] [CrossRef]
  83. McCarthy, D.J.; Chen, Y.; Smyth, G.K. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012, 40, 4288–4297. [Google Scholar] [CrossRef]
  84. Liu, R.; Holik, A.Z.; Su, S.; Jansz, N.; Chen, K.; Leong, H.S.; Blewitt, M.E.; Asselin-Labat, M.L.; Smyth, G.K.; Ritchie, M.E. Why weight? Modelling sample and observational level variability improves power in RNA-seq analyses. Nucleic Acids Res. 2015, 43, e97. [Google Scholar] [CrossRef] [PubMed]
  85. Leek, J.T. svaseq: Removing batch effects and other unwanted noise from sequencing data. Nucleic Acids Res. 2014, 42, e161. [Google Scholar] [CrossRef] [PubMed]
  86. Dos Santos, K.C.G.; Desgagné-Penix, I.; Germain, H. Custom selected reference genes outperform pre-defined reference genes in transcriptomic analysis. BMC Genom. 2020, 21, 35. [Google Scholar] [CrossRef] [PubMed]
  87. Andersen, C.L.; Jensen, J.L.; Ørntoft, T.F. Normalization of real-time quantitative reverse transcription-PCR data: A model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004, 64, 5245–5250. [Google Scholar] [CrossRef] [PubMed]
  88. Vandesompele, J.; De Preter, K.; Pattyn, F.; Poppe, B.; Van Roy, N.; De Paepe, A.; Speleman, F. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol 2002, 3, RESEARCH0034. [Google Scholar] [CrossRef] [PubMed]
  89. Xie, F.; Wang, J.; Zhang, B. RefFinder: A web-based tool for comprehensively analyzing and identifying reference genes. Funct. Integr. Genomics 2023, 23, 125. [Google Scholar] [CrossRef]
  90. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
  91. Pinheiro, J.; Bates, D. Mixed-Effects Models in S and S-PLUS; Springer: New York, NY, USA, 2000. [Google Scholar]
  92. Pinheiro, J.; Bates, D.; DebRoy, S.; Sarkar, D.; EISPACK_authors; Heisterkamp, S.; Van_Willigen, B.; Ranke, J.; R_Core_Team. nlme: Linear and Nonlinear Mixed Effects Models. 2023. Available online: https://CRAN.R-project.org/package=nlme (accessed on 1 December 2023).
  93. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014, 15, 550. [Google Scholar] [CrossRef]
  94. Kolberg, L.; Raudvere, U.; Kuzmin, I.; Vilo, J.; Peterson, H. gprofiler2—An R package for gene list functional enrichment analysis and namespace conversion toolset g:Profiler. F1000Res 2020, 9, ELIXIR-709. [Google Scholar] [CrossRef]
  95. Supek, F.; Bosnjak, M.; Skunca, N.; Smuc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE 2011, 6, e21800. [Google Scholar] [CrossRef]
Ncrna 10 00027 g001
Figure 1. Heatmap of differentially expressed long non-coding RNAs (lncRNAs) detected by RNA-Seq of adult worm couples obtained from mice across ten time points at 0, 0.25, 1, 3, 6, 9, 12, 24, 48, and 96 h after treatment (Treat) with a single, curative dose of praziquantel (400 mg/kg) delivered by oral gavage at 7 weeks post-infection. Hierarchical clustering of differentially expressed lncRNAs (lines) from adult worm samples (columns) harvested at various time points following PZQ treatment of infected mice, as indicated by the color bar at the top and the blue color scale at right. These results were obtained by re-analyses of the RNA-Seq data from McCusker et al. [11] using the S. mansoni lncRNA reference transcriptome published by Maciel et al. [26]. Gene expression levels were measured by RNA-Seq and are shown as Z-scores, which are the number of standard deviations below (blue, downregulated) or above (red, upregulated) the mean expression value among treated and control samples for each gene, as indicated in the scale at right. 189 lncRNAs were considered significantly differentially expressed, being 127 long intergenic ncRNAs, 58 antisense lncRNAs and 4 sense lncRNAs (FDR < 0.05).
Figure 1. Heatmap of differentially expressed long non-coding RNAs (lncRNAs) detected by RNA-Seq of adult worm couples obtained from mice across ten time points at 0, 0.25, 1, 3, 6, 9, 12, 24, 48, and 96 h after treatment (Treat) with a single, curative dose of praziquantel (400 mg/kg) delivered by oral gavage at 7 weeks post-infection. Hierarchical clustering of differentially expressed lncRNAs (lines) from adult worm samples (columns) harvested at various time points following PZQ treatment of infected mice, as indicated by the color bar at the top and the blue color scale at right. These results were obtained by re-analyses of the RNA-Seq data from McCusker et al. [11] using the S. mansoni lncRNA reference transcriptome published by Maciel et al. [26]. Gene expression levels were measured by RNA-Seq and are shown as Z-scores, which are the number of standard deviations below (blue, downregulated) or above (red, upregulated) the mean expression value among treated and control samples for each gene, as indicated in the scale at right. 189 lncRNAs were considered significantly differentially expressed, being 127 long intergenic ncRNAs, 58 antisense lncRNAs and 4 sense lncRNAs (FDR < 0.05).
Ncrna 10 00027 g001
Ncrna 10 00027 g002
Figure 2. Heatmaps of differentially expressed long non-coding RNAs (lncRNAs) detected by RNA-Seq of adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of praziquantel (100 mg/kg) delivered by intraperitoneal injection at 4 (a) or 7 (b) weeks post-infection. Hierarchical clustering of differentially expressed lncRNAs (lines) from five adult worm sample replicates (columns) harvested 14 h following praziquantel (PZQ) treatment of infected mice or from five control (CTL) replicates, as indicated by the color bars at the top and the color legend at right. These results were obtained by re-analysis of the RNA-Seq data from McCusker et al. [11] using the S. mansoni lncRNA reference transcriptome published by Maciel et al. [26]. Gene expression levels were measured by RNA-Seq and are shown as Z-scores, which are the number of standard deviations below (downregulated, light blue) or above (upregulated, dark blue) the mean expression value among treated and control samples for each gene, as indicated in the scale at right. (a) In the 4-weeks experiment, 181 lncRNAs were considered significantly differentially expressed, being 108 long intergenic ncRNAs, 57 antisense lncRNAs and 16 sense lncRNAs (FDR < 0.05). (b) In the 7-weeks experiment, 348 lncRNAs were considered significantly differentially expressed, being 236 long intergenic ncRNAs, 95 antisense lncRNAs and 17 sense lncRNAs (FDR < 0.05).
Figure 2. Heatmaps of differentially expressed long non-coding RNAs (lncRNAs) detected by RNA-Seq of adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of praziquantel (100 mg/kg) delivered by intraperitoneal injection at 4 (a) or 7 (b) weeks post-infection. Hierarchical clustering of differentially expressed lncRNAs (lines) from five adult worm sample replicates (columns) harvested 14 h following praziquantel (PZQ) treatment of infected mice or from five control (CTL) replicates, as indicated by the color bars at the top and the color legend at right. These results were obtained by re-analysis of the RNA-Seq data from McCusker et al. [11] using the S. mansoni lncRNA reference transcriptome published by Maciel et al. [26]. Gene expression levels were measured by RNA-Seq and are shown as Z-scores, which are the number of standard deviations below (downregulated, light blue) or above (upregulated, dark blue) the mean expression value among treated and control samples for each gene, as indicated in the scale at right. (a) In the 4-weeks experiment, 181 lncRNAs were considered significantly differentially expressed, being 108 long intergenic ncRNAs, 57 antisense lncRNAs and 16 sense lncRNAs (FDR < 0.05). (b) In the 7-weeks experiment, 348 lncRNAs were considered significantly differentially expressed, being 236 long intergenic ncRNAs, 95 antisense lncRNAs and 17 sense lncRNAs (FDR < 0.05).
Ncrna 10 00027 g002
Ncrna 10 00027 g003
Figure 3. Analysis of long non-coding RNAs (lncRNAs) differentially expressed upon distinct praziquantel treatment regimens. The Venn diagram shows the number of lncRNAs that were detected as differentially expressed in one experiment alone, or at the same time in three different experiments: RNA-Seq from adult worm couples obtained from mice across ten time points at 0, 0.25, 1, 3, 6, 9, 12, 24, 48, and 96 h after treatment with a single, curative dose of PZQ (400 mg/kg) delivered by oral gavage at 7 weeks post-infection (“Time-course” experiment, pink circle); (ii) RNA-Seq from adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 4 weeks post-infection (“4-weeks” experiment, blue circle); (iii) RNA-Seq from adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 7 weeks post-infection (“7-weeks” experiment, green circle).
Figure 3. Analysis of long non-coding RNAs (lncRNAs) differentially expressed upon distinct praziquantel treatment regimens. The Venn diagram shows the number of lncRNAs that were detected as differentially expressed in one experiment alone, or at the same time in three different experiments: RNA-Seq from adult worm couples obtained from mice across ten time points at 0, 0.25, 1, 3, 6, 9, 12, 24, 48, and 96 h after treatment with a single, curative dose of PZQ (400 mg/kg) delivered by oral gavage at 7 weeks post-infection (“Time-course” experiment, pink circle); (ii) RNA-Seq from adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 4 weeks post-infection (“4-weeks” experiment, blue circle); (iii) RNA-Seq from adult worm couples obtained from mice 14 h after treatment with a single, sub-lethal dose of PZQ (100 mg/kg) delivered by intraperitoneal injection at 7 weeks post-infection (“7-weeks” experiment, green circle).
Ncrna 10 00027 g003
Ncrna 10 00027 g004
Figure 4. Validation by RT-qPCR of lincRNAs detected as differentially expressed upon praziquantel treatment in vivo. Ten lincRNAs detected by RNA-Seq as differentially expressed were selected for in vitro RT-qPCR assays, namely: (a) SmLINC101519-IBu, (b) SmLINC105115-IBu, (c) SmLINC110492-IBu, (d) SmLINC121232-IBu, (e) SmLINC133371-IBu, (f) SmLINC142502-IBu, (g) SmLINC159037-IBu, (h) SmLINC161393-IBu, (i) SmLINC163938-IBu, and (j) SmLINC172840-IBu. For each of the ten selected lincRNAs, the expression profiles obtained with RNA-Seq are shown on the right (blue line) as TMM (trimmed mean of M values), whereas the RT-qPCR results are shown on the left (black line). For the RT-qPCR data, mean ± SEM from three biological replicates is shown, and a linear mixed-effects statistical model was applied. * p  <  0.05, ** p  <  0.01, *** p  <  0.001, and **** p  <  0.0001.
Figure 4. Validation by RT-qPCR of lincRNAs detected as differentially expressed upon praziquantel treatment in vivo. Ten lincRNAs detected by RNA-Seq as differentially expressed were selected for in vitro RT-qPCR assays, namely: (a) SmLINC101519-IBu, (b) SmLINC105115-IBu, (c) SmLINC110492-IBu, (d) SmLINC121232-IBu, (e) SmLINC133371-IBu, (f) SmLINC142502-IBu, (g) SmLINC159037-IBu, (h) SmLINC161393-IBu, (i) SmLINC163938-IBu, and (j) SmLINC172840-IBu. For each of the ten selected lincRNAs, the expression profiles obtained with RNA-Seq are shown on the right (blue line) as TMM (trimmed mean of M values), whereas the RT-qPCR results are shown on the left (black line). For the RT-qPCR data, mean ± SEM from three biological replicates is shown, and a linear mixed-effects statistical model was applied. * p  <  0.05, ** p  <  0.01, *** p  <  0.001, and **** p  <  0.0001.
Ncrna 10 00027 g004
Ncrna 10 00027 g005
Figure 5. RNA-Seq expression profiles of selected lincRNAs measured in control assays at different S. mansoni life-cycle stages. The expression levels of the ten selected lincRNAs are shown, namely: (a) SmLINC101519-IBu, (b) SmLINC105115-IBu, (c) SmLINC110492-IBu, (d) SmLINC121232-IBu, (e) SmLINC133371-IBu, (f) SmLINC142502-IBu, (g) SmLINC159037-IBu, (h) SmLINC161393-IBu, (i) SmLINC163938-IBu, and (j) SmLINC172840-IBu. These S. mansoni lincRNAs were selected after re-analysis of RNA-Seq public datasets of parasites collected from mice treated with PZQ in the “time-course experiment” [11]. The y-axis shows the expression level (shown as TMM—trimmed mean of M values) for each lincRNA in control RNA-Seq assays, as compiled by Silveira et al. [31], at the stage indicated in the x-axis as follows: miracidia/sporocysts (M/S), cercariae (C), schistosomula (S), adult males (M), and adult females (F).
Figure 5. RNA-Seq expression profiles of selected lincRNAs measured in control assays at different S. mansoni life-cycle stages. The expression levels of the ten selected lincRNAs are shown, namely: (a) SmLINC101519-IBu, (b) SmLINC105115-IBu, (c) SmLINC110492-IBu, (d) SmLINC121232-IBu, (e) SmLINC133371-IBu, (f) SmLINC142502-IBu, (g) SmLINC159037-IBu, (h) SmLINC161393-IBu, (i) SmLINC163938-IBu, and (j) SmLINC172840-IBu. These S. mansoni lincRNAs were selected after re-analysis of RNA-Seq public datasets of parasites collected from mice treated with PZQ in the “time-course experiment” [11]. The y-axis shows the expression level (shown as TMM—trimmed mean of M values) for each lincRNA in control RNA-Seq assays, as compiled by Silveira et al. [31], at the stage indicated in the x-axis as follows: miracidia/sporocysts (M/S), cercariae (C), schistosomula (S), adult males (M), and adult females (F).
Ncrna 10 00027 g005
Ncrna 10 00027 g006
Figure 6. Analysis of the co-expression modules identified in WGCNA. (a) This plot shows the proportion of differentially expressed genes (DEGs) inside each of the 15 modules. The total number of genes inside each module is proportional to the point size; blue points (representing “brown2”, “floralwhite”, “plum” and “plum2”) indicate modules that contain the differentially expressed lncRNAs tested by RT-qPCR. The most significantly enriched and relevant gene ontology (GO) terms, including biological process and molecular function terms, are shown for “brown2” (b) and “floralwhite” (c) modules. The size of each point is proportional to the precision, also known as gene ratio (i.e., the proportion of genes in the input list that are annotated to the function). The colors show the statistical significance of the enrichment.
Figure 6. Analysis of the co-expression modules identified in WGCNA. (a) This plot shows the proportion of differentially expressed genes (DEGs) inside each of the 15 modules. The total number of genes inside each module is proportional to the point size; blue points (representing “brown2”, “floralwhite”, “plum” and “plum2”) indicate modules that contain the differentially expressed lncRNAs tested by RT-qPCR. The most significantly enriched and relevant gene ontology (GO) terms, including biological process and molecular function terms, are shown for “brown2” (b) and “floralwhite” (c) modules. The size of each point is proportional to the precision, also known as gene ratio (i.e., the proportion of genes in the input list that are annotated to the function). The colors show the statistical significance of the enrichment.
Ncrna 10 00027 g006
Ncrna 10 00027 g007
Figure 7. Single-cell clusters expression profiles of three RT-qPCR validated differentially expressed lncRNAs matching the expression profile of three differentially expressed protein-coding genes modulated under praziquantel exposure [11]. UMAP plots show the expression enrichment of genes and are colored by gene expression (blue = low, red = high). The scale represents log10(UMIs+1). The genes with expression enrichment in tegument are (a) SmLINC121232-IBu and (b) Smp_086480 (SmTAL2); in vitellocytes (c) SmLINC110492-IBu and (d) Smp_076320 (Myb/SANT-like DNA-binding domain-containing protein); and in neurons and muscles (e) SmLINC101519-IBu and (f) Smp_105220 (Lymphocyte antigen 6B).
Figure 7. Single-cell clusters expression profiles of three RT-qPCR validated differentially expressed lncRNAs matching the expression profile of three differentially expressed protein-coding genes modulated under praziquantel exposure [11]. UMAP plots show the expression enrichment of genes and are colored by gene expression (blue = low, red = high). The scale represents log10(UMIs+1). The genes with expression enrichment in tegument are (a) SmLINC121232-IBu and (b) Smp_086480 (SmTAL2); in vitellocytes (c) SmLINC110492-IBu and (d) Smp_076320 (Myb/SANT-like DNA-binding domain-containing protein); and in neurons and muscles (e) SmLINC101519-IBu and (f) Smp_105220 (Lymphocyte antigen 6B).
Ncrna 10 00027 g007
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Jardim Poli, P.; Fischer-Carvalho, A.; Tahira, A.C.; Chan, J.D.; Verjovski-Almeida, S.; Sena Amaral, M. Long Non-Coding RNA Levels Are Modulated in Schistosoma mansoni following In Vivo Praziquantel Exposure. Non-Coding RNA 2024, 10, 27. https://doi.org/10.3390/ncrna10020027

AMA Style

Jardim Poli P, Fischer-Carvalho A, Tahira AC, Chan JD, Verjovski-Almeida S, Sena Amaral M. Long Non-Coding RNA Levels Are Modulated in Schistosoma mansoni following In Vivo Praziquantel Exposure. Non-Coding RNA. 2024; 10(2):27. https://doi.org/10.3390/ncrna10020027

Chicago/Turabian Style

Jardim Poli, Pedro, Agatha Fischer-Carvalho, Ana Carolina Tahira, John D. Chan, Sergio Verjovski-Almeida, and Murilo Sena Amaral. 2024. "Long Non-Coding RNA Levels Are Modulated in Schistosoma mansoni following In Vivo Praziquantel Exposure" Non-Coding RNA 10, no. 2: 27. https://doi.org/10.3390/ncrna10020027

APA Style

Jardim Poli, P., Fischer-Carvalho, A., Tahira, A. C., Chan, J. D., Verjovski-Almeida, S., & Sena Amaral, M. (2024). Long Non-Coding RNA Levels Are Modulated in Schistosoma mansoni following In Vivo Praziquantel Exposure. Non-Coding RNA, 10(2), 27. https://doi.org/10.3390/ncrna10020027

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop