Introduction
Tetracapsuloides bryosalmonae is an obligate endoparasite that causes proliferative kidney disease (PKD) in salmonid fish. PKD is considered as one of the most serious emerging temperature-dependent parasitic freshwater diseases of salmonids in the Northern Hemisphere, causing significant losses on farms and in natural populations (Okamura et al., Reference Okamura, Anderson, Longshaw, Feist and Canning2001, Reference Okamura, Hartikainen, Schmidt-Posthaus and Wahli2011). Tetracapsuloides bryosalmonae belongs to the Class Myxozoa within the Phylum Cnidaria and belongs to the relative species poor Subclass Malacosporea (Giribet and Edgecombe, Reference Giribet and Edgecombe2020). The parasite has a two-host life-cycle, alternating between an invertebrate host (several species of freshwater bryozoans) and a fish host (salmonids) (Okamura et al., Reference Okamura, Hartikainen, Schmidt-Posthaus and Wahli2011). Fish host invasion occurs via epidermal mucous cells (e.g. gills, skin; Longshaw et al., Reference Longshaw, Le Deuff, Harris and Feist2002), followed by dissemination to internal organs, primarily the kidney and spleen via the vascular system. Sporogonic development takes place in the kidney of permissive fish hosts, and spores infective to bryozoans are released in fish urine (Hedrick et al., Reference Hedrick, Baxa, De Kinkelin and Okamura2004). The development of PKD symptoms is temperature-dependent. Exposure to water temperatures ≥15°C elicits a strong immune response, kidney hypertrophy and anaemia, which may lead to death (Hedrick et al., Reference Hedrick, MacConnell and de Kinkelin1993; Bettge et al., Reference Bettge, Wahli, Segner and Schmidt-Posthaus2009). As a consequence, the decline of native salmonid fish populations in several countries has been linked with a combination of increasing temperatures and emergence of PKD (Wahli et al., Reference Wahli, Knuesel, Bernet, Segner, Pugovkin, Burkhardt-Holm, Escher and Schmidt-Posthaus2002; Okamura et al., Reference Okamura, Hartikainen, Schmidt-Posthaus and Wahli2011; Sudhagar et al., Reference Sudhagar, Kumar and El-Matbouli2020).
Despite the high economic and ecological significance of PKD, the genomic resources for T. bryosalmonae are limited to a few hundred expressed sequence tags and some other targeted loci (e.g. Saulnier et al., Reference Saulnier, Brémont and De Kinkelin1996; Holland et al., Reference Holland, Okamura, Hartikainen and Secombes2011; Hartikainen et al., Reference Hartikainen, Gruhl and Okamura2014, Reference Hartikainen, Filippenko, Okamura and Vasemägi2015; Carraro et al., Reference Carraro, Hartikainen, Jokela, Bertuzzo and Rinaldo2018). Only very recently, an intersect transcriptome of T. bryosalmonae consisting of de novo assembled transcripts present in both rainbow trout (Oncorhynchus mykiss) and bryozoan host (Fredericella sultana) has been reported (Faber et al., Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021). However, the T. bryosalmonae transcriptome recovered from the dead-end host O. mykiss might lack genes essential for the parasite to reach sporogony in natural fish host. Knowledge of the parasite genomic composition and gene expression is therefore needed to better understand the underlying mechanisms of parasite life stages, host–parasite interactions, resistance, tolerance and virulence. Moreover, genomic information may help to predict potential new drug targets and drugs/drug-like compounds against the parasite (Cowell and Winzeler, Reference Cowell and Winzeler2019). The climate-change driven emergence of PKD underlies the urgency for the identification of effective chemotherapeutic targets against T. bryosalmonae that can be used to protect stocks and improve the welfare of farmed salmonids.
In this study, we report the de novo assembly and annotation of the partial T. bryosalmonae transcriptome based on dual RNA sequencing (Westermann et al., Reference Westermann, Gorski and Vogel2012) of both parasite and host utilizing infected brown (Salmo trutta) trout kidney tissue. We also identify T. bryosalmonae genes that are differentially expressed in relation to parasite load (PL) to obtain a better understanding of host–parasite interactions. Based on searches against known drug target proteins in the ChEMBL database, we identify four promising protein drug targets (consisting of six genes) whose suppression might obstruct the parasite growth and mitigate infections. Thus, these findings prime future experimental approaches to develop and test new control strategies against PKD. Furthermore, the assembled T. bryosalmonae transcriptome will be useful in understanding the parasite biology and serve as a valuable reference for characterizing the parasite development during infection.
Materials and methods
Sampling and measurement of kidney hyperplasia
Sampling and phenotyping methods of juvenile (0+) brown trout used in this study were reported in Debes et al. (Reference Debes, Gross and Vasemägi2017). DNA was extracted from kidney tissue as in Aljanabi and Martinez (Reference Aljanabi and Martinez1997). Fourteen fish kidneys corresponding to eight infected and six uninfected fish from three wild populations were selected. All the infected fish originated from River Mustoja (Estonia) and were collected in September, 2014, whereas the uninfected fish originated from two parasite-free rivers (Preedi and Vodja, Estonia; three per river) and were collected in September, 2015. Sampling from multiple populations was necessary as T. bryosalmonae prevalence in the study area is often 100% (Dash and Vasemägi, Reference Dash and Vasemägi2014; Debes et al., Reference Debes, Gross and Vasemägi2017), prohibiting the sampling of both infected and uninfected individuals from the same location. We also quantified a typical disease trait associated with PKD, kidney swelling (renal hyperplasia), estimated as a body-length-standardised kidney depth, as described by Debes et al. (Reference Debes, Gross and Vasemägi2017).
PL measurement
PL (a proxy for the inverse of resistance) was determined from infected fish kidney tissue by qPCR using previously described protocols (Bruneaux et al. Reference Bruneaux, Visse, Gross, Pukk, Saks and Vasemägi2016; Debes et al. Reference Debes, Gross and Vasemägi2017). For each sample, we quantified the occurrence of two target DNA sequences: a 166 bp fragment of T. bryosalmonae 18S rDNA sequence (GenBank Accession U70623) and 74 bp fragment of the S. salar nuclear DNA sequence (GenBank Accession BT049744.1, Prefoldin subunit 6). Simultaneous quantification of both DNA fragments enabled us to standardize the amount of parasite DNA relative to host DNA. Three technical replicates were run for both targets in each sample. Based on PL, the infected samples were divided into categories of moderate (n = 4) and severe (n = 4) infection. On average the moderately infected fish exhibited lower kidney swelling measurements, compared to severely infected fish (~twice the kidney depth; Supplementary Fig. 1).
RNA isolation, library preparation and RNA sequencing
We preserved half of the kidney cross-section for RNA isolation by freezing them in liquid N2 and storing them at −80°C. The other half of the kidney cross-section was stored in the same way but used for estimation of PL using qPCR. Total RNA was extracted from the kidney cross-section of the 14 individuals using NucleoSpin® RNA kit. RNA quality was tested using the Advanced Analytical Fragment Analyzer. Quantity of each total RNA sample was measured with a Nanodrop ND-2000 spectrophotometer. Library preparation was done according to Illumina TruSeq® Stranded mRNA Sample Preparation Guide (part # 15031047) starting from 100 ng of total RNA. The steps involved in cDNA library constructions were as follows: purification of poly-A containing mRNA, fragmentation of mRNA using divalent cations, first-strand (using reverse transcriptase and random primers) and second strand (using DNA Polymerase I and RNase H) synthesis of cDNA, purification and enrichment with PCR. A unique Illumina TruSeq indexing adapter was ligated to each sample during the adapter ligation step to be able to pool several samples for sequencing. The quality of the prepared libraries was evaluated with Advanced Analytical Fragment Analyzer. The concentration of each library was quantified with Qubit® Fluorometric Quantitation, Life Technologies. Due to adapter peaks, the libraries were purified twice. The samples were normalized and pooled for the automated cluster preparation that was carried out using Illumina cBot station. The 14 libraries were pooled into a single pool. The samples were sequenced with Illumina HiSeq 2500 instrument using TruSeq v3 sequencing chemistry in three lanes. Paired-end sequencing with 100 bp read length was used, followed by 6 bp index run. The base calling was performed using Illumina's standard bcl2fastq2 software. RNA sequencing was performed at the Finnish Functional Genomics Centre, Turku Centre for Biotechnology.
Data processing and removal of host RNA sequence
The quality assessment of sequenced reads was performed with FastQC 0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) both before and after pre-processing. Trimmomatic (v. 0.35) was used for adapter trimming and read pre-processing (parameters: ILLUMINACLIP = TruSeq3-s.e..fa; HEADCROP = 10; SLIDINGWINDOW = 4:20; LEADING = 5; TRAILING = 5; MINLEN = 40). In order to assemble a host-free parasite transcriptome, the filtered reads were first mapped to the Atlantic salmon genome (GCF_000233375.1) and the unmapped reads were subsequently mapped to the northern pike (Esox lucius, Esociformes) genome (GCF_000721915.2) using Hisat2 (v. 2.0.1; score-min = L,0,−0.4) (Kim et al., Reference Kim, Langmead and Salzberg2015). We used chromosome level assembly of Atlantic salmon genome (Lien et al., Reference Lien, Koop, Sandve, Miller, Kent, Nome, Hvidsten, Leong, Minkley, Zimin, Grammes, Grove, Gjuvsland, Walenz, Hermansen, Von Schalburg, Rondeau, Di Genova, Samy, Olav Vik, Vigeland, Caler, Grimholt, Jentoft, Inge Våge, De Jong, Moen, Baranski, Palti, Smith, Yorke, Nederbragt, Tooming-Klunderud, Jakobsen, Jiang, Fan, Hu, Liberles, Vidal, Iturra, Jones, Jonassen, Maass, Omholt and Davidson2016), as it is evolutionarily close to brown trout (Ahmad et al., Reference Ahmad, Debes, Palomar and Vasemägi2018; Lemopoulos et al., Reference Lemopoulos, Uusi-Heikkilä, Huusko, Vasemägi and Vainikka2018), and at the time of analysis, genome assembly of S. trutta was not available. Esociformes represents the closest order to Salmoniformes (Rondeau et al., Reference Rondeau, Minkley, Leong, Messmer, Jantzen, von Schalburg, Lemon, Bird and Koop2014). The paired-end reads from the infected fish samples that were not aligned to either the salmon or pike genome were further taxonomically classified by using Kraken2 (Wood et al., Reference Wood, Lu and Langmead2019) with a custom database. The custom Kraken2 database consisted of two fish genomes (Atlantic salmon and Northern pike), four available myxozoan genomes (Enteromyxum leei, Sphaeromyxa zaharoni, Kudoa iwatai and Thelohanellus kitauei), both host (Salmo trutta, 3177 sequences, mean GC content = 43.2%) and parasite (T. bryosalmonae, 741 sequences, mean GC content = 33.6%) nucleotide sequences available in NCBI and ENA (Supplementary Table 1), as well as, unpublished contigs for T. bryosalmonae, derived from its bryozoan host F. sultana (mean GC content = 27.1%; H. Hartikainen unpublished), along with Kraken2-provided archaea, bacteria, plasmid, viral, human, fungi, protozoa (protists) and UniVec_Core sequences. The reads that were classified within the NCBI bony fish division were excluded before the de novo assembly of the parasite transcriptome.
De novo transcriptome assembly of T. bryosalmonae
The putative parasite transcriptome was assembled using Trinity (v. 2.5.1) assembler (Grabherr et al., Reference Grabherr, Haas, Yassour, Levin, Thompson, Amit, Adiconis, Fan, Raychowdhury, Zeng, Chen, Mauceli, Hacohen, Gnirke, Rhind, Di Palma, Birren, Nusbaum, Lindblad-Toh, Friedman and Regev2011) with the default parameters. CD-HIT-EST (Fu et al., Reference Fu, Niu, Zhu, Wu and Li2012) was used to cluster 95% identical transcripts and isoforms. The paired-end reads from both the uninfected samples that were not mapped to any fish genome and the non-bony fish reads from the infected samples were aligned on the clustered transcripts using Bowtie2 (Langmead and Salzberg, Reference Langmead and Salzberg2012). The read count matrix for all samples was generated using the summarize Overlaps function implemented in R for further evaluation and differential expression analyses. In order to filter out transcripts originating from the host or from potential contamination, the longest isoforms of the assembled transcripts were subjected to Kraken2 classification and similarity searching using BLASTN and BLASTX (Altschul et al., Reference Altschul, Gish, Miller, Myers and Lipman1990) against NR databases (only the top hits, dust = no, soft masking = false, E-value cut-off = 1 × 10−4).
Filtering and completeness assessment of the parasite transcriptome
The putative parasite transcripts were filtered using the following criteria: (i) at least 40 reads mapped, (ii) 90% or more of the mapped reads originate from infected fish and (iii) retained only those transcripts that matched to GenBank's invertebrate division or no hit was found using Kraken2, BLASTN and BLASTX. In order to assess the completeness of the transcriptome, these transcripts were subjected to Benchmarking Universal Single-Copy Orthologs (BUSCO,v. 3.0) (Simão et al., Reference Simão, Waterhouse, Ioannidis, Kriventseva and Zdobnov2015) analysis against 978 Metazoan single-copy orthologs.
Recently, Faber et al. (Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021) reported an intersect transcriptome of T. bryosalmonae, common to both rainbow trout (O. mykiss) and bryozoan host (F. sultana) which consisted of 7362 contigs showing 50.1% completeness (excluding fragmented hits) against the BUSCO eukaryotic dataset (303 genes). Underestimation of completeness for the myxozoan genomes/transcriptomes using BUSCO metazoan set has already been reported (Hartigan et al., Reference Hartigan, Kosakyan, Pecková, Eszterbauer and Holzer2020; Yahalomi et al., Reference Yahalomi, Atkinson, Neuhof, Chang, Philippe, Cartwright, Bartholomew and Huchon2020). Therefore, in order to identify BUSCO metazoan genes that are either entirely missing or present in all myxozoans, we analysed genome annotations-based coding sequences (CDS) and transcriptome assemblies for seven species: Thelohanellus kitauei, Kudoa iwatai, Henneguya salminicola, Myxobolus squamalis, Myxobolus cerebralis, Sphaerospora molnari and Ceratonova shasta (Supplementary Table 2). The T. kitauei transcriptome assembly was not available and therefore only CDS extracted from the genome using available genome feature file (GFF) were used. For K. iwatai, the genome annotation was not available and therefore, only the transcriptome reads were aligned to the genome using splicing-aware HISAT2 followed by Stringtie (Pertea et al., Reference Pertea, Pertea, Antonescu, Chang, Mendell and Salzberg2015) to assemble alignment into potential transcripts (reference guided transcriptome assembly). The splicing events in the K. iwatai genome were summarized using ASTALAVISTA web server (Foissac and Sammeth, Reference Foissac and Sammeth2007). The assembled transcripts were then extracted from the genome using gffread (Pertea and Pertea, Reference Pertea and Pertea2020). BUSCO genes were considered as present in myxozoans even if it was a fragmented match.
Annotation of putative parasite transcripts
The functional annotation of the parasite transcripts was performed using two sources, namely; KEGG automatic annotation server (KAAS; Moriya et al., Reference Moriya, Itoh, Okuda, Yoshizawa and Kanehisa2007) and Trinotate (Bryant et al., Reference Bryant, Johnson, DiTommaso, Tickle, Couger, Payzin-Dogru, Lee, Leigh, Kuo, Davis, Bateman, Bryant, Guzikowski, Tsai, Coyne, Ye, Freeman, Peshkin, Tabin, Regev, Haas and Whited2017). For the KEGG Orthology search, the nucleotide sequences of these transcripts were searched against a manually curated set of 39 species (including three cnidarians; Nematostella vectensis, Exaiptasia pallida and Hydra vulgaris). The assigned KEGG ortholog (KO) identifiers were searched for pathways in the reference, as well as, the organism specific pathways (for the three cnidarians) using KEGG mapper (https://www.genome.jp/kegg/tool/map_pathway1.html). The assigned KOs were also tested for pathway enrichment against all the reference KOs using the enrichKEGG (minGSSize = 1, maxGSSize = 5000, other default parameters) function implemented in clusterprofiler. Organism specific pathway enrichment using KO identifiers was not possible in clusterprofiler.
Open reading frames (ORFs; minimum of 50 amino acid length) for the parasite transcripts were identified using Transdecoder (Haas et al., Reference Haas, Papanicolaou, Yassour, Grabherr, Blood, Bowden, Couger, Eccles, Li, Lieber, Macmanes, Ott, Orvis, Pochet, Strozzi, Weeks, Westerman, William, Dewey, Henschel, Leduc, Friedman and Regev2013) using homology searches (BLASTP and hmmscan (Eddy, Reference Eddy2011) results of longest ORFs against NCBI refseq database and pfam, respectively). The reduced value of minimum ORF length was used to consider the shorter transcripts (200–300 bp, n = 526) and fragmented genes for the subsequent annotation. Because myxozoans are underrepresented in the public databases and difficult to annotate, ORFs > 150 bp length (50 aa) were retained even when no homology was found. The resulting protein sequences were then used for BLASTP searches and for signal peptide prediction in Trinotate annotation suite. For Trinotate annotation (report E-value <1 × 10−3), we used the results from BLASTP and BLASTX homology searches against Swissprot, transmembrane helices predictions using TMHMM Server (v. 2.0), pfam domain search using hmmscan, signal peptides prediction using SignalP 5.0, and BLASTP homology search against Refseq proteins. For the genes with multiple isoforms and protein translations, only the hits with the maximum number of GO terms were reported. Finally, the ribosomal RNAs were predicted using RNAmmer 1.2 Server (Lagesen et al., Reference Lagesen, Hallin, Rødland, Stærfeldt, Rognes and Ussery2007).
Overlap between two T. bryosalmonae transcriptomes
In order to find the overlap between our T. bryosalmonae transcriptome and the recently reported T. bryosalmonae intersect transcriptome (from the dead-end vertebrate host, rainbow trout [O. mykiss], and the bryozoan F. sultana; Faber et al., Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021), we carried out reciprocal BLASTN. A transcript (longest isoform in case of multiple isoforms) was considered as present in Faber et al. (Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021) transcriptome if its top hit (95% minimum identity) matched the best hit in the reciprocal BLASTN. As the scaffolds in Faber et al. (Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021) transcriptome were relatively long (N50 = 2479) compared to ours (N50 = 657), for the reciprocal BLASTN, we allowed up to 3 hits per scaffold when Faber et al. (Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021) transcriptome was queried against our transcripts. All isoforms of the same gene were considered as present if the longest isoform was present in Faber et al. (Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021) transcriptome and vice versa. The GO terms between the genes found in Faber et al. (Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021) transcriptome vs the genes unique to our study were compared using WEGO 2 (Ye et al., Reference Ye, Zhang, Cui, Liu, Wu, Cheng, Xu, Huang, Li, Zhou, Zhang, Bolund, Chen, Wang, Yang, Fang and Shi2018).
Differential gene expression analysis of the parasite transcripts
In order to identify potential links between differences in gene expression and parasite development, gene-level differential expression (DE) analysis of T. bryosalmonae between the moderately and severely infected fish was carried out using the DEseq function implemented in the R package DESeq2 (Love et al., Reference Love, Huber and Anders2014). For this purpose, the read counts of the different isoforms of the same genes were summed up to obtain gene-level read counts. An explorative default significance level (FDR <0.1) was used to identify differentially expressed genes since we were mostly interested in detecting candidates that need further validation using RT-qPCR.
Identification of drug targets
In order to screen for potential drug targets, T. bryosalmonae transcriptome was searched against ChEMBL protein sequences using TBLASTN (soft masking = false, E-value cut-off = 1 × 10−3). To identify parasite-specific targets and avoid potential drug toxicity to the host, only those transcripts lacking BLASTX hits against the Salmo salar proteome (soft masking = false, E-value cut-off = 1 × 10−3) were retained. To increase the sensitivity of the homology search, we focused on targets with an alignment length >100 amino acids based on TBLASTN search. The protein sequences of the putative parasite transcripts corresponding the top matched drug targets (based on E-value) were subjected to homology-based structure modelling in Phyre2 (http://www.sbg.bio.ic.ac.uk/phyre2) using the intensive mode (Kelley et al., Reference Kelley, Mezulis, Yates, Wass and Sternberg2015). Only those sequences were subject to homology modelling which had a conserved alignment at active or binding sites of the drug targets. The best-supported structures were subsequently superimposed with the template structures using Chimera (v. 1.14; Pettersen et al., Reference Pettersen, Goddard, Huang, Couch, Greenblatt, Meng and Ferrin2004) for visualization.
Results
After the quality control and filtering, 493.1 Million (401.6 M paired-end and 91.5 M single-end) reads from the host and parasite were retained (range 22.9–47.2 M reads per sample). Among these, 266 M (220.7 M paired-end and 45.3 M single-end) reads originated from infected fish. After mapping on fish reference genomes, a total of 7.6 M out of 220.7 M (3.4% of the total) paired-end reads from infected fish, that were not assigned to bony fishes by Kraken2, were assembled into 16499 transcripts (12676 genes) using Trinity. Individual transcripts with 95% or more similarity were clustered into 14741 merged transcripts (12621 genes) using CD-HIT-EST with default parameters (Table 1).
Table 1. Summary statistics of the de novo transcriptome assembly of T. bryosalmonae
The overall average alignment rates of infected and uninfected samples against merged transcripts were 83.73 and 33.24%, respectively. Because parasite reads should originate only from infected fish, we used this criterion as an additional step to filter out host transcripts. More specifically, if >10% of the aligned reads of the particular transcript originated from the uninfected fish, the transcript was not assigned as a putative parasite sequence.
A large proportion of merged transcripts were not assigned to any taxonomic group (Fig. 1A, Kraken2 = 6508, BLASTN = 10283 and BLASTX = 6440) reflective of the highly diverged nature of myxozoans from other cnidarians (Chang et al., Reference Chang, Neuhof, Rubinstein, Diamant, Philippe, Huchon and Cartwright2015) and a small number of available myxozoan genomes. Among the sequences that were classified using BLASTN, the majority of hits were from vertebrate and salmonid fishes (Fig. 1B). However, relatively more hits to cnidarians (Hydra vulgaris and Nematostella vectensis and Exaiptasia pallida) were identified using BLASTX (Fig. 1C). Only 5993 (40%), 130 (0.9%) and 22 (0.2%) transcripts were assigned to Myxozoa by Kraken2 (Fig. 1D), BLASTX and BLASTN, respectively. Among these, 11 transcripts (Fig. 1A) were designated as myxozoan by all three searches and consisted of following T. bryosalmonae genes: antigen PKX101 (seven transcripts including four isoforms of the same gene), minicollagen (Ncol-1, DN1177_c0_g1_i1 and DN14479_c0_g1_i1), DNA repair protein RAD51-like protein (DN15328_c0_g1_i1, excluded because of low coverage) and 60S ribosomal protein L18 (RP-L18, DN3431_c0_g1_i1). BLASTN also assigned some additional transcripts to known genes of T. bryosalmonae, including 28S large subunit ribosomal RNA gene (DN11757_c1_g2_i1 and DN11757_c1_g4_i1), S-adenosylmethionine synthetase (DN7631_c0_g1_i1) and ribosomal protein rpl23a-like gene (DN1241_c0_g1_i1), respectively.
Fig. 1. Summary of the taxonomic classification of assembled transcripts. Bar plots in (A) indicate the number of transcripts assigned (or unclassified) to the most frequent GenBank's taxonomic divisions. Venn diagrams illustrate the classification agreement between the BLASTN, BLASTX and Kraken based assignments. Within the invertebrates and vertebrates, only myxozoans and bony fish are shown in the Venn diagrams, respectively. The bar plots in (B–D) show the species to which most of the transcripts were assigned by BLASTN, BLASTX and Kraken, respectively.
Transcriptome characterization and completeness
The overall distribution of GC content of 14741 merged transcripts was bimodal (Fig. 2A), suggesting that the vertebrate host and myxozoan parasite genomes have different GC content. Previous studies have identified myxozoan transcriptomes and genomes to be predominantly AT-rich (Yang et al., Reference Yang, Xiong, Zhou, Huo, Miao, Ran, Liu, Zhang, Feng, Wang, Wang, Wang and Yao2014; Foox et al., Reference Foox, Ringuette, Desser and Siddall2015; Yahalomi et al., Reference Yahalomi, Atkinson, Neuhof, Chang, Philippe, Cartwright, Bartholomew and Huchon2020). After filtering, the putative parasite transcripts (Fig. 2A; n = 3427; mean GC = 31.5%), from here on referred to as T. bryosalmonae transcriptome, and transcripts belonging to bony fishes were well separated (Fig. 2A; n = 1027; mean GC = 48.56%). The GC content of the T. bryosalmonae transcriptome was similar to other myxozoan genomes (GC content 23.6–37.5%) (Yahalomi et al., Reference Yahalomi, Atkinson, Neuhof, Chang, Philippe, Cartwright, Bartholomew and Huchon2020). The GC content of brown trout (43.9%) (Malachowicz et al., Reference Malachowicz, Wenne and Burzynski2017) also broadly matched with the GC content of the transcripts belonging to bony fishes. The T. bryosalmonae transcriptome consisted of more than 2 million bases with an average transcript length of 587 bp (range: 201–3388 bp). Among these transcripts, 929 represented isoform reconstructs of 407 genes, while the remaining 2498 transcripts corresponded to genes without isoforms. The read counts for the large majority of the putative parasite transcripts were positively correlated with the independent estimate of relative PL obtained via qPCR assay (Fig. 2B; 92.3% of transcript r ≥ 0.5) suggesting that our bioinformatic pipeline (Supplementary Fig. 2) successfully retained T. bryosalmonae reads. Moreover, the positive correlation estimates suggested that the coverage of the majority of assembled transcripts was sufficient to accurately reflect the relative PL.
Fig. 2. Characterization of T. bryosalmonae transcripts: (A) The distributions of GC content of all assembled transcripts (grey; n = 14741; mean GC = 35.51%), putative parasitic sequences (light grey; n = 3427; mean GC = 31.67%) and transcripts belonging to bony fishes in BLAST results (black; n = 1027; mean GC = 48.56%). (B) The distribution of Pearson's correlation coefficients between the raw read counts of putative parasite transcripts and the qPCR-inferred relative parasite load.
The assessment of T. bryosalmonae transcriptome completeness using BUSCO analysis indicated that only 26% of the core single copy metazoan genes were completely (14.5%) or partially (11.6%) present in our assembly (Supplementary Table 2 and 3). The recently published T. bryosalmonae transcriptome (Faber et al., Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021) on the other hand exhibited two times higher completeness as 48.3% of the metazoan single copy genes were completely (38.1%) or partially (10.1%) present. A total of 58 BUSCO genes (32 complete single; 2 complete duplicated and 24 fragmented) were identified in our T. bryosalmonae transcriptome assembly which were missing in the published transcriptome Faber et al. (Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021). A similar pattern of relatively low BUSCO completeness can be observed in other published myxozoan transcriptome and genome annotations (Hartigan et al., Reference Hartigan, Kosakyan, Pecková, Eszterbauer and Holzer2020; Yahalomi et al., Reference Yahalomi, Atkinson, Neuhof, Chang, Philippe, Cartwright, Bartholomew and Huchon2020). We subsequently evaluated the presence and absence of BUSCO core single copy metazoan genes in the published myxozoan parasite datasets (Fig. 3A–B, Supplementary Table 3). Despite variation in coverage and completeness of individual transcriptomes and genomes, we identified a set of BUSCO metazoan genes found in all myxozoans (n = 231, Fig. 3A). Furthermore, we observed a set of BUSCO metazoan genes (n = 232, 23.7%) that were missing in all seven myxozoan transcriptome and genome annotations. This suggests that these genes are either missing or highly divergent in myxozoans (Fig. 3B).
Fig. 3. BUSCO single-copy metazoan genes and proteases in Myxozoa (A) Venn diagram showing the number of BUSCO single-copy genes identified (complete + fragmented) in each species. The T and G after the name of species represent de novo assembled transcriptome and CDS extracted from the genomes, respectively. (B) Venn diagram showing the number of BUSCO single-copy genes missing in each species. (C) Distribution of proteases in current T. bryosalmonae transcriptome in relation to earlier studies. The numbers in each category represent protease classes found in each study.
KEGG ortholog assignment and ribosomal RNA
Out of 3427 putative parasite transcripts, only 1742 transcripts (1509 genes) were assigned to 1172 unique KEGG orthologs (KO) genes, most of which belong to enzymes (657 transcripts; 38%), membrane trafficking (243 transcripts; 14%), spliceosome (206 transcripts; 12%), exosome (173 transcripts; 10%) and chromosome and associated proteins (165 transcripts; 9%) (Supplementary Table 4). Among the enzymes, which were the most abundant proteins, 14% and 11% were protein kinases and peptidases (proteases) and inhibitors, respectively (Supplementary Table 4). The proteases were further distributed into Cysteine (32%), Metallo (27%), Threonine (16%), Serine (16%), and Aspartic proteases (8%) (Supplementary Table 5). Even though the Aspartic proteases were the least numerous proteases (5 genes), their expression was higher than for other protease classes, suggesting potential functional importance. Some of the proteases with the highest expression (average FPKM >1000) were: cathepsin D (CTSD, 2 genes), puromycin-sensitive aminopeptidase (NPEPPS), cytosol aminopeptidase (LAP3) and mitochondrial-processing peptidase subunit beta (PMPCB).
In both reference and cnidarian-specific pathway searches, a large proportion of these transcripts were assigned to metabolic pathways (ko01100; 262 transcripts), spliceosome (ko03040; 141 transcripts), ribosome (ko03010; 74 transcripts), RNA transport (ko03013; 72 transcripts), endocytosis (ko04144; 59 transcripts), oxidative phosphorylation (ko00190; 54 transcripts) and phagosome (ko04145; 54 transcripts) (Supplementary Table 6, Supplementary Fig. 3A). Similar to what has been found for T. kitauei (Yang et al., Reference Yang, Xiong, Zhou, Huo, Miao, Ran, Liu, Zhang, Feng, Wang, Wang, Wang and Yao2014), many T. bryosalmonae transcripts were assigned to genes involved in multicellular metazoan signalling pathways e.g. mTOR (40 transcripts) and Wnt (22 transcripts) signalling pathways (Supplementary Table 6). The reference-based pathway search also assigned genes to different human disease pathways, which, however, were not observed during organism-specific searches (Supplementary Fig 3B). The enrichment analysis against reference pathways revealed that disease-related genes are enriched among our transcripts and there is an overlap of genes in different disease-related pathways (Supplementary Fig. 4). The excess of human disease-related genes have also been reported in the cnidarian starlet sea anemone (N. vectensis) (Sullivan and Finnerty, Reference Sullivan and Finnerty2007), which most likely indicates that disease-associated genes in humans are biased towards ancient genes (Domazet-Loso and Tautz, Reference Domazet-Loso and Tautz2008).
Among the transcripts with KO gene assignment, 430 transcripts corresponded to 197 gene isoforms. Tetracapsuloides bryosalmonae genes with multiple isoforms included cyclin-dependent kinase 5 (CDK5; 5 isoforms), annexin A7/11 (ANXA7_11; 4 isoforms), ATPase family AAA domain-containing protein 3A/B (ATAD3A_B; 4 isoforms), histone-binding protein RBBP4 (RBBP4; 4 isoforms), SWI/SNF-related matrix-associated actin-dependent regulator of chromatin subfamily A member 2/4 (SMARCA2_4; 4 isoforms), cyclin T (CCNT; 4 isoforms), electron-transferring-flavoprotein dehydrogenase (ETFDH; 4 isoforms), V-type H + -transporting ATPase subunit H (ATPeV1H; 4 isoforms). Among the five CDK5 isoforms, three have premature stop codons while the other two were partial sequences. Similarly, two of the ANXA7_11 isoforms included premature stop codons. The genome-guided transcriptome assembly of K. iwatai (assembled for BUSCO, 23896 transcripts, 10556 genes) enabled us to assess the genome-wide pattern of splicing events in the myxozoans. In K. iwatai, a total of 23130 splicing events occurred in 5533 genes and intron retention was more frequent (6826, 29.5%) followed by alternative acceptor sites (3437, 14.8%), alternative donor sites (2861, 12.4%) and exon skipping (952, 4.1%). Multiple isoforms of the above-mentioned T. bryosalmonae genes (except ETFDH) were also observed in K. iwatai reference guided transcriptome assembly (Supplementary Table 7) further confirming that these isoforms are not the artefacts of de novo assembly. Thus, the presence of different isoforms may also contribute in increasing the functional diversity of the encoded proteins in Myxozoa. However, further experimental studies are needed to assign and validate the functions of different isoforms of T. bryosalmonae genes.
Only a single ribosomal RNA (DN11757_c1_g4_i1) was predicted from our parasite transcripts using RNAmmer showing >2000 bp alignment overlap (99% BLASTN identity) with the partial 28S large subunit ribosomal RNA gene of T. bryosalmonae (FJ981822.1).
Protein coding transcripts and GO annotation
A total of 4046 protein sequences were initially produced from the translation of the parasite transcriptome, which is higher than the number of input transcripts (3324) because multiple protein sequences per transcripts were produced due to short minimum protein length (50 aa). These included; 467 complete (11.5%, mean length 153.8 aa), 12095’- partial (missing start codon) (29.9%; mean length 157.6 aa), 1847 internal (45.6%; mean length 171.6 aa) (missing both start and stop codons) and 5233’- partial (missing stop codon) proteins (12.9%; mean length 143.7 aa).
By employing these sequences in the Trinotate annotation pipeline, a total of 1944 genes were assigned to 25195 GO terms comprising; 11325 (45%) biological processes, 7955 (31.6%) cellular component and 5915 (23.5%) molecular functions. The most abundant GO terms for cellular component, biological process and molecular function included integral component of the membrane, protein phosphorylation and ATP binding, respectively (Supplementary Fig. 5). The complete Trinotate annotation report along-with KEGG ortholog gene identifiers and transcript and protein sequences are provided in Supplementary Table 8.
Comparison between the two T. bryosalmonae transcriptomes
Using an adjusted reciprocal BLASTN, 81.5% (2793 transcripts, 2377 genes) of the transcripts were assigned to 1615 contigs of the recently reported intersect transcriptome (from bryozoan and rainbow trout, Faber et al., Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021). A total of 634 transcripts (528 genes), did not have a hit within the intersect transcriptome and are hence unique to our assembly (Supplementary Fig. 6A). In comparison of the transcripts that matched to the Faber et al., Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021 transcriptome, the unique transcripts had higher relative proportion of genes related to Ribosome (5.13%), Ribosome biogenesis (3%), Ubiquitin system (2%) and Peptidases and inhibitors (1.6%) based on BRITE hierarchies (Supplementary Fig. 6B). GO terms comparison using WEGO 2 showed that the unique T. bryosalmonae transcripts have higher (P < 0.05) proportion of genes related to the cellular component terms: Endoplasmic reticulum membrane (GO:0005789) and Mitochondrial protein complex (GO:0098798), and following molecular function terms: Structural molecule activity (GO:0005198) and Structural component of ribosome (GO:0003735) (Supplementary Fig. 6C and D).
Differentially expressed T. bryosalmonae genes linked to PL
A total of 55 T. bryosalmonae genes were differentially expressed at FDR <0.1 between severely and moderately infected fish (Fig. 4; Supplementary Table 9). Similar to the whole transcriptome annotation, only a proportion of differentially expressed genes (n = 31) retrieved functional annotation using KAAS and Trinotate. In the severely infected samples, the polar capsule related minicollagens (TBNCOL-1, TBNCOL-2 and TBNCOL-4) were downregulated, indicating that severely and moderately infected hosts may differ in the extent of polar capsule formation of T. bryosalmonae. Furthermore, in severely infected fish several lysosome and lipid metabolism-related genes were downregulated. These include cathepsin L (CTSL; 2 genes), Niemann-Pick C2 protein (NPC2), hydroxymethylglutaryl-CoA synthase (HMGCS), hydroxymethylglutaryl-CoA reductase (NADPH) (HMGCR) and PI-PLC X domain-containing protein 2 (PLCX2). Interestingly, differential expression of CTSLs has recently been found between the proliferative blood- and sporogonic gill stages of the myxozoan S. molnari (Hartigan et al., Reference Hartigan, Kosakyan, Pecková, Eszterbauer and Holzer2020). Therefore, it is possible that cathepsin L is also linked to sporogony in T. bryosalmonae. In severely infected fish downregulation was also inferred for the expression of glutathione peroxidase (gpx; the highest mean expression level among differentially expressed genes; AvgFPKM in the Supplementary Table 9), which is an antioxidant that protects the parasite from the damage caused by free radicals released from the host immune cells (Changklungmoa et al., Reference Changklungmoa, Chaithirayanon, Cheukamud, Chaiwichien, Osotprasit, Samrit, Sobhon and Kueakhai2018).
Fig. 4. Differential transcript expression and GO enrichment: (A) Heatmap showing variance stabilized expression of the differentially expressed transcripts between the moderately and severely infected hosts. (B) Volcano plot showing the fold-change differences between the moderately and severely infected hosts. Light blue and red points represent the down- and up-regulated genes in severely infected fish, respectively.
Among the upregulated genes in the severely infected fish, four were retrotransposon genes (annotated as Retrovirus-related Pol polyprotein from type-1 retrotransposable element R2; pol). The upregulated gene with the highest fold change corresponded to low-density lipoprotein receptor class A-like protein 3 (MT070967.1, 100% BLASTN identity) gene. Other upregulated genes included Talin (TLN), TRAF2 and NCK interacting kinase (TNIK), Tetraspanin-9 (TSN9), Vacuolar-sorting protein (snf7), importin-5 (IPO5), calcium-dependent secretion activator (CADPS), Rho GTPase-activating protein 39 (ARHGAP39), synaptosomal-associated protein 25 (SNAP25), AGPAT8 lysocardiolipin and lysophospholipid acyltransferase (LCLAT1) and large subunit ribosomal protein L3e (RP-L3e).
Drug targets and homology modelling
After stringent screening of 3427 putative parasite transcripts against known drug target proteins in the ChEMBL database, four potential chemotherapeutic targets, Endoglycoceramidase (two genes; EGCase1 and EGCase2), Legumain-like protease, Carbonic anhydrase 2 and Pancreatic lipase-related protein 2 (two genes) were identified (Table 2). The two EGCase transcripts encode 283 and 382 amino acid long partial protein sequences with 43.8% similarity. We subsequently individually modelled the 3D structures of the two EGCase and CA2 protein sequences by using Phyre2. For both EGCase1 and EGCase2, >90% of residues were modelled at >90% confidence during homology modelling. Both structures showed a high overlap with the template structure (Fig. 5, Supplementary Fig. 7 and 8). The structure of the EGCase2 (longer alignment length with the template) is shown in Fig. 5 in which the low confidence residues (92–114) were unreliably modelled into a loop because they were missing in the template sequence (Supplementary Fig. 8). The template EGCase chosen by Phyre2 (B chain of endoglycoceramidase ii from Rhodococcus sp. (PDB ID 2OYL) for the homology modelling was bound to an inhibitor molecule (cellobiose-like imidazole) using polar interactions with Lys66, His135, Asp137, Asn232, Glu233 and Glu351 (Caines et al., Reference Caines, Hancock, Tarling, Wrodnigg, Stick, Stütz, Vasella, Withers and Strynadka2007). Apart from Lys66, which was missing from the partial template sequence, the corresponding His65, Asp67, Asn172, Glu173 and Glu277 in the T. bryosalmonae EGCase2 (DN1345_c0_g1_i1, Fig. 5c) acquired the same positions in the binding site (Fig. 5C). The T. bryosalmonae EGCase1 (DN11176_c0_g2_i1) exhibits three amino acids occupying the similar 3D confirmation as the template structure's binding site amino acids (Supplementary Fig 7C). These conserved binding site residues of T. bryosalmonae Endoglycoceramidase II represent a promising target for the chemotherapeutic drugs with cellobiose-like imidazole inhibitor-like structures.
Table 2. Potential drug targets of the de novo assembled T. bryosalmonae transcripts
Fig. 5. Prediction of the EGCase2 (DN1345_c0_g1_i1) and CA2 (DN419_c0_g2_i1) three-dimensional structures in T. bryosalmonae. (pink). (A) The parasite EGCase2 (green) structure is superimposed on the chain B of the template endoglycoceramidase ii from rhodococcus sp. In EGCase2, the inhibitor (yellow) is shown in ball and stick representation. Similar binding cavities between template and target can be seen around the inhibitor molecule. (B) The same structure shown in ribbons to illustrate the consistency between the template and target structures. (C) The interactions between the EGCase II amino acids and inhibitor at the binding site are shown with black lines as reported in Caines et al. (Reference Caines, Hancock, Tarling, Wrodnigg, Stick, Stütz, Vasella, Withers and Strynadka2007). (D) The parasite CA2 (green) structure is superimposed on the chain A of the Schistosoma mansoni carbonic anhydrase (SmCA, pink). The surface is made transparent to show the active site. (E) The same structure shown in ribbons to illustrate the consistency between SmCA and the predicted T. bryosalmonae structures. (F) The active site (His76, Gln107, Glu122, and Thr missing) and Zinc binding Histidines (His109, His111 and His1136) residues of T. bryosalmonae CA2 also acquired similar orientation as the SmCA (Da'dara et al., Reference Da'dara, Angeli, Ferraroni, Supuran and Skelly2019).
More than 90% of residues of the T. bryosalmonae CA2 were modelled at >90% confidence using Phyre2. The predicted structure (Fig. 5D and E, green) was superimposed on chain A of the parasite Schistosoma mansoni carbonic anhydrase (SmCA, PDB ID 6QQM, pink). In SmCA, the active site residues are His88, Gln115, Glu129, and Thr231 and the Zinc binding Histidines are His117, His119 and His142 (Da'dara et al., Reference Da'dara, Angeli, Ferraroni, Supuran and Skelly2019). These important residues are also conserved in T. bryosalmonae CA2 (Active site: His76, Gln107, Glu122 and Zinc binding Histidines: His109, His111 and His136) while only Threonine is missing. These residues also acquired similar orientation as in the SmCA (Fig. 5F) thus the inhibitors used for the S. mansoni carbonic anhydrase might also be effective for T. bryosalmonae CA2. Finally, we did not predict 3D structures of LGMN and PLRP2s as they are partial sequences and lacked important active site or binding site residues based on alignments with ChEMBL hits (Supplementary Fig. 9 and 10).
Discussion
This study reports the first annotated transcriptome of T. bryosalmonae derived from wild, infected brown trout using dual RNA sequencing. By employing multi-step filtering of parasite reads followed by de novo assembly, our work provides new insights into the differentially expressed T. bryosalmonae genes in the fish and identifies six proteins that could act as potential targets for the development of antiparasitic drugs against PKD. As such, this work demonstrates the utility of the assembled T. bryosalmonae transcriptome and promotes experimental approaches to develop and test new control strategies against PKD.
A major challenge for the de novo transcriptome assembly using dual-RNA sequencing data is to correctly separate host and parasite RNA reads (Schulze et al., Reference Schulze, Schleicher, Guthke and Linde2016; Videvall et al., Reference Videvall, Cornwallis, Ahrén, Palinauskas, Valkiūnas and Hellgren2017). In earlier studies targeting myxozoan transcriptomes within fish hosts, RNA was isolated from the appendages visible by eye, e.g. plasmodia from the K. iwatai infected gilthead seabream (Chang et al., Reference Chang, Neuhof, Rubinstein, Diamant, Philippe, Huchon and Cartwright2015) or cysts from the T. kitauei infected common carp (Yang et al., Reference Yang, Xiong, Zhou, Huo, Miao, Ran, Liu, Zhang, Feng, Wang, Wang, Wang and Yao2014), and after careful washing steps that were followed by in silico filtering of reads to avoid inclusion of host sequences during the parasite transcriptome or genome assembly. Yet, more recent works have shown that host tissue contamination is still an important issue in next generation sequencing of myxozoans (Foox et al., Reference Foox, Ringuette, Desser and Siddall2015; Alama-Bermejo et al., Reference Alama-Bermejo, Meyer, Atkinson, Holzer, Wiśniewska, Kolísko and Bartholomew2020; Hartigan et al., Reference Hartigan, Kosakyan, Pecková, Eszterbauer and Holzer2020). Furthermore, isolation of microscopic T. bryosalmonae specimens from salmonid kidney tissue is challenging and our initial RNAseq data was dominated (more than 97% of total sequences) by brown trout RNA reads. Therefore, in order to produce a host-free parasite transcriptome, we employed a rigorous multi-step strategy, comprising host genome mapping, pre-assembly metagenomic classification and post-assembly filtering. Recently, Alama-Bermejo et al. (Reference Alama-Bermejo, Meyer, Atkinson, Holzer, Wiśniewska, Kolísko and Bartholomew2020) tackled the problem of host contamination removal from the myxozoan C. shasta transcriptome by employing a two-step filtering approach. In the first step, host reads were removed by aligning reads to host and parasite reference genomes. During the second step (transcript-level filtering), Alama-Bermejo et al. (Reference Alama-Bermejo, Meyer, Atkinson, Holzer, Wiśniewska, Kolísko and Bartholomew2020) used BLASTN against both host and parasite genomes to further exclude the host transcripts, while contaminant sequences (of bacterial, fungal or viral origins) were excluded based on BLASTX (Alama-Bermejo et al., Reference Alama-Bermejo, Meyer, Atkinson, Holzer, Wiśniewska, Kolísko and Bartholomew2020). In contrast, T. bryosalmonae transcriptome was not published at the time of the analysis. In addition, only a good quality genome of Salmo salar, representing a genetically close sister species to brown trout, was available for us at the time of the analysis. Therefore, to overcome these limitations, we included multiple myxozoan genomes in the Kraken2 database to classify the unmapped reads and assembled parasite transcriptome using reads that were not classified to ‘bony fishes’. Furthermore, reads from the uninfected fish that mapped on our assembled transcripts provided a useful control step which enabled to effectively remove both known and unknown transcripts originating from the fish host. Some homology between host and parasite short reads is expected and therefore we allowed up to 10% of the total reads per transcript to map to the uninfected fish. Similar to the second step transcript-level filtering in C. shasta (Alama-Bermejo et al., Reference Alama-Bermejo, Meyer, Atkinson, Holzer, Wiśniewska, Kolísko and Bartholomew2020), we also excluded transcripts that showed significant (e.g. bacterial, viral and vertebrate) BLAST/Kraken2 hits during the post-assembly filtering, thus any potential bacterial and human contamination during the sample collection and library preparation was also filtered out. Finally, the majority of transcripts with positive correlation between raw read counts and PL measured by an independent approach (qPCR; Bruneaux et al., Reference Bruneaux, Visse, Gross, Pukk, Saks and Vasemägi2016; Debes et al., Reference Debes, Gross and Vasemägi2017) indicated that our pipeline successfully removed the host sequences and captured true RNA reads from T. bryosalmonae.
Myxozoans are highly reduced metazoans with extremely small diverged genome sizes (22.5–188.5 Mb), which makes it challenging to evaluate the completeness of assembled transcriptomes and genomes (Chang et al., Reference Chang, Neuhof, Rubinstein, Diamant, Philippe, Huchon and Cartwright2015; Hartigan et al., Reference Hartigan, Kosakyan, Pecková, Eszterbauer and Holzer2020; Yahalomi et al., Reference Yahalomi, Atkinson, Neuhof, Chang, Philippe, Cartwright, Bartholomew and Huchon2020). By using a large set of published data, we identified a set of 231 BUSCO genes that were detected in all myxozoans. Based on this smaller set of BUSCO genes, the estimated T. bryosalmonae transcriptome assembly completeness is close to 56.3% (130 of 231). Thus, we suggest that this core set of 231 BUSCO genes can be used as an alternative to estimate the completeness of the assemblies in future myxozoan studies. Conversely, we also observed that 23.7% of metazoan single-copy gene genes were not identified in any published myxozoan genome and transcriptome assemblies. As the myxozoans have gone through an extreme reduction in body plan and genome size (Chang et al., Reference Chang, Neuhof, Rubinstein, Diamant, Philippe, Huchon and Cartwright2015), it is possible that at least some of these genes have been truly lost during myxozoan evolution. However, we cannot exclude, at this point, the alternative explanation that unidentified BUSCO genes are still present in myxozoans, but because of highly divergent nature, were not detected using BUSCO.
Currently, there are no licenced drugs available for efficient PKD treatment. Earlier studies demonstrated an efficiency of malachite green (triarylmethane dye; Clifton-Hadley and Alderman, Reference Clifton-Hadley and Alderman1987) and fumagillin (antibiotic of the fungus Aspergillus fumigatus; Le Gouvello et al., Reference Le Gouvello, Pobel, Richards and Gould1999) for treatment against PKD. Malachite green has been used for many fish species to control fungal attacks, protist infections and other diseases, such as caused by helminths (Srivastava et al., Reference Srivastava, Sinha and Roy2004). Fumagillin has been shown to prevent clinical outbreaks of whirling disease caused by the myxosporean parasite Myxobolus cerebralis in rainbow trout (El-Matbouli and Hoffmann, Reference El-Matbouli and Hoffmann1991). In addition, fumagillin is an effective antibiotic against the microsporidian pathogen Nosema apis in honey bees (Burnham, Reference Burnham2019). However, because of the accumulation of malachite green in fish and its potential teratogenic and carcinogenic properties, its use is, unsurprisingly, restricted for food fish. Similarly, it has been shown that fumagillin treatment leads to a dramatic loss in haematopoietic tissue in rainbow trout (Wishkovsky et al., Reference Wishkovsky, Groff, Lauren, Toth and Hedrick1990). Therefore, given the high economic losses incurred by PKD on the freshwater salmonid aquaculture industry of Europe and North America (Feist, Reference Feist2004; Burkhardt-Holm et al., Reference Burkhardt-Holm, Giger, GUttinger, Ochsenbein, Peter, Scheurer, Segner, Staub and Suter2005; Borsuk et al., Reference Borsuk, Reichert, Peter, Schager and Burkhardt-Holm2006; Okamura et al., Reference Okamura, Hartikainen, Schmidt-Posthaus and Wahli2011), novel treatment targets and compounds are urgently needed.
The assembly of the T. bryosalmonae transcriptome enabled us to identify potential drug targets while avoiding potential toxicity to the salmonid host due to drug target homology. The top drug candidate (based on the lowest BLAST E-value) endoglycoceramidase (EGCase) is an enzyme that hydrolyses acidic and neutral glycosphingolipids to oligosaccharides and ceramides. Our homology modelling analysis predicted that the T. bryosalmonae EGCase is expected to have the similar binding site as the original template and hence the same inhibitors, such as iminocyclitols, may affect this enzyme activity in T. bryosalmonae. Among the free-living Cnidaria, EGCase has been previously reported to occur in Cyanea nozakii (Jellyfish) (Horibata et al., Reference Horibata, Okino, Ichinose, Omori and Ito2000) and Hydra vulgaris (Hydrozoa) (Horibata et al., Reference Horibata, Sakaguchi, Okino, Iida, Inagaki, Fujisawa, Hama and Ito2004). Interestingly, EGCase is involved in a unique dietary glycosphingolipids catabolism pathway in Hydra (Horibata et al., Reference Horibata, Sakaguchi, Okino, Iida, Inagaki, Fujisawa, Hama and Ito2004). However, until this study, endoglycoceramidases have not been reported in Myxoza. The presence of two copies of this gene in a reduced cnidarian genome makes it possible that the glycosphingolipids catabolism pathway also plays a role in the parasite lipid metabolism.
The second identified drug target legumain-like protease, an asparaginyl endopeptidase, occurs in the gut of the hard tick (Ixodes ricinus) where it participates in the digestion of blood meal haemoglobin (Sojka et al., Reference Sojka, Hajdušek, Dvořák, Sajid, Franta, Schneider, Craik, Vancová, Burešová, Bogyo, Sexton, McKerrow, Caffrey and Kopáček2007). A similar legumain enzyme is also present in the parasitic blood fluke Schistosoma mansoni (Davis et al., Reference Davis, Nanduri and Watson1987). The proteolytic activity of legumain can be reduced by Aza-peptidyl Michael acceptors, which are inhibitors for the legumains (Götz et al., Reference Götz, James, Hansell, Dvořák, Seshaadri, Sojka, Kopáček, McKerrow, Caffrey and Powers2008; Ovat et al., Reference Ovat, Muindi, Fagan, Brouner, Hansell, Dvořák, Sojka, Kopáček, McKerrow, Caffrey and Powers2009; Rendón-Gandarilla et al., Reference Rendón-Gandarilla, Ramón-Luing, de los, Ortega-López, Rosa de Andrade, Benchimol and Arroyo2013). For example, Aza-peptidyl Michael acceptors inhibit legumain activity in the living cells of protist Trichomonas vaginalis and reduce the parasitic cytoadherence to the host cells (Rendón-Gandarilla et al., Reference Rendón-Gandarilla, Ramón-Luing, de los, Ortega-López, Rosa de Andrade, Benchimol and Arroyo2013). Similarly, Aza-peptidyl Michael acceptors (Götz et al., Reference Götz, James, Hansell, Dvořák, Seshaadri, Sojka, Kopáček, McKerrow, Caffrey and Powers2008) and epoxide inhibitors (Ovat et al., Reference Ovat, Muindi, Fagan, Brouner, Hansell, Dvořák, Sojka, Kopáček, McKerrow, Caffrey and Powers2009) have been shown to inhibit legumain in other evolutionary distant parasite species, such as S. mansoni and I. ricinus.
The third drug target, Carbonic anhydrase 2 (CA2), belongs to the plasma membrane anchored alpha carbonic anhydrases, which catalyse the reversible hydration of carbon dioxide and aid in intracellular pH regulation. Carbonic anhydrases are considered as potent anti-parasitic drug targets against protists, such as Trypanosoma, Leishmania (Vermelho et al., Reference Vermelho, Capaci, Rodrigues, Cardoso, Mazotto and Supuran2017) and, recently, against blood flukes, schistosomes (Da'dara et al., Reference Da'dara, Angeli, Ferraroni, Supuran and Skelly2019; Angeli et al., Reference Angeli, Pinteala, Maier, Simionescu, Da'dara, Skelly and Supuran2020). In the predicted three-dimensional model, the active site- and Zn binding residues of T. bryosalmonae CA2 retained similar positions as in SmCA, suggesting that similar therapeutic inhibitors may be used to target this protein in T. bryosalmonae. The schistosome alpha carbonic anhydrase plays a role in parasite virulence as its suppression using RNAi significantly reduced the establishment of parasite within host (Da'dara et al., Reference Da'dara, Angeli, Ferraroni, Supuran and Skelly2019). Several sulphonamide, thiol and hydroxamate substances have been shown to inhibit CA2 in Trypanosoma and Leishmania schistosomes (Vermelho et al., Reference Vermelho, Capaci, Rodrigues, Cardoso, Mazotto and Supuran2017; Angeli et al., Reference Angeli, Pinteala, Maier, Simionescu, Da'dara, Skelly and Supuran2020).
The fourth drug target, pancreatic lipase-related protein 2 (PLRP2), is a secretory lipase enzyme that in vertebrates hydrolyses both phospholipids and galactolipids and helps in milk fat digestion (Berton et al., Reference Berton, Sebban-Kreuzer, Rouvellac, Lopez and Crenon2009). This gene is also present in genomes of several cnidarians, including Exaiptasia pallida (NCBI Gene IDs: 110254178 & 110231612), Stylophora pistillata (NCBI Gene IDs: 111321164 & 110253264) and Nematostella vectensis, (NCBI Gene IDs: 5496642). Presence of PLRP2 in invertebrates and cnidarians suggest a potential role in nutrient acquisition and lipid metabolism. In humans, this enzyme has been employed as a target for the lipase inhibitors, such as Orlistat, for obesity treatment (Ballinger and Peikin, Reference Ballinger and Peikin2002). Interestingly, Orlistat has also been shown to inhibit in vitro growth of a protist parasite Giardia duodenalis (Hahn et al., Reference Hahn, Seeber, Kolodziej, Ignatius, Laue, Aebischer and Klotz2013). Altogether, earlier work on endoglycoceramidases, legumain-like proteases, carbonic anhydrases and lipase inhibitors suggests that the identified potential anti-parasitic drug targets represent promising chemotherapeutic candidates also against myxozoan T. bryosalmonae.
The pathways and enzymes assigned by the KEGG orthology provided important functional information about the genetic make-up of T. bryosalmonae. Among identified pathways, the endocytosis pathway is expected to have an important role in T. bryosalmonae. Specifically, endocytosis is involved in acquiring nutrients in other Myxozoa (Yang et al., Reference Yang, Xiong, Zhou, Huo, Miao, Ran, Liu, Zhang, Feng, Wang, Wang, Wang and Yao2014). In addition, endocytosis of the cellular content of the host phagocyte (that has already engulfed the parasite cell) has been observed in T. bryosalmonae (Morris and Adams, Reference Morris and Adams2008). Endocytosis has also been observed during the transfer of material from the primary cell to the enclosed secondary and tertiary cells in T. bryosalmonae (Morris and Adams, Reference Morris and Adams2008), suggesting that endocytosis may have wider functions than acquiring nutrients. Furthermore, it has been shown that the internal cells (both secondary and tertiary) of T. bryosalmonae are formed via the engulfment of other cells by primary cell (Morris, Reference Morris2010) and thus, the same molecular machinery i.e. endocytosis pathway genes are likely to be involved in these processes.
Peptidases (proteases) were the second most abundant enzymes found T. bryosalmonae. These are important catabolic enzymes that hydrolyze proteins into smaller polypeptides or single amino acids. A key role of the proteases in the nutrient digestion has already been proposed for the myxozoan T. kitauei (Yang et al., Reference Yang, Xiong, Zhou, Huo, Miao, Ran, Liu, Zhang, Feng, Wang, Wang, Wang and Yao2014). Furthermore, proteases are known virulence factors in many parasites and considered as important targets for the antiphrastic drug development (McKerrow et al., Reference McKerrow, Caffrey, Kelly, Loke and Sajid2006; Pina-Vazquez et al., Reference Pina-Vazquez, Reyes-Lopez, Ortiz-Estrada, de la Garza and Serrano-Luna2012; Siqueira-Neto et al., Reference Siqueira-Neto, Debnath, McCall, Bernatchez, Ndao, Reed and Rosenthal2018). More recently, proteases have also been suggested as potential drug targets in myxozoan S. molnari (Hartigan et al., Reference Hartigan, Kosakyan, Pecková, Eszterbauer and Holzer2020). In the current T. bryosalmonae assembly, the distribution of KEGG assigned proteases showed higher abundance of cysteine (32%) and metallo (27%) proteases, while a smaller proportion of serine (16%), threonine (16%) and aspartic proteases (8%) were identified (Fig. 3C). In contrast, recently published T. bryosalmonae assembly (Faber et al., Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021) reported a slightly higher proportion of metallo proteases in relation to cysteine and serine proteases. This difference may reflect the host-specificity of proteases as the former transcriptome was common to both rainbow trout and bryozoan host, while the latter was based on a native salmonid host, brown trout. Furthermore, despite low completeness of current transcriptome assembly, the protease composition of T. bryosalmonae was very similar to S. molnari, C. shasta, K. iwatai and M. cerebralis. Noteably, T. kitauei exhibited a clearly different protease composition compared to other myxozoans (Fig. 3C).
Between severely and moderately infected fish, altogether, 55 T. bryosalmonae genes were differentially expressed. Most interestingly, a higher expression of minicollagens, which are expected to be only expressed during myxozoan sporogenesis, in moderately infected fish suggests that more parasites had reached the sporogonic stage in moderately infected relative to the severely infected fish. Similarly, Niemann-Pick C2 protein and ammonium transporter Rh, whose homologs are found in the polar capsule of the myxozoan C. shasta (ecdysteroid-regulated 16 kDa-like and ammonium transporter Rh type C 2-like respectively) (Piriatinskiy et al., Reference Piriatinskiy, Atkinson, Park, Morgenstern, Brekhman, Yossifon, Bartholomew and Lotan2017) were upregulated in moderately infected fish. Thus, it is possible that fish with low PL shows increased abundance of sporogonic stages of T. bryosalmonae. However, further experimental time-series analysis is needed to confirm this hypothesis. Among the genes upregulated in severely infected fish was Talin, which is related to the cytoskeletal adhesion to the extracellular matrix (Klapholz and Brown, Reference Klapholz and Brown2017). Recently, the upregulation of the Talin gene in C. shasta has been shown to associate with the transformation from a localized (intestine) to a proliferated infection of other organs in rainbow trout (Alama-Bermejo et al., Reference Alama-Bermejo, Holzer and Bartholomew2019). Thus, Talin, as a putative virulence factor, may also play role in proliferation and adhesion in T. bryosalmonae.
Based on the 231 metazoan genes detected in all published myxozoan transcriptomes/genomes using BUSCO, we estimated that the assembled transcriptome recovered approximately ~56% of T. bryosalmonae genes (~92% for intersect transcriptome of Faber et al., Reference Faber, Yoon, Shaw, de Paiva Alves, Wang, Qi, Okamura, Hartikainen, Secombes and Holland2021). The removal of transcripts with low expression levels (less than 40 reads), is one likely reason behind the low BUSCO score and genes with low expression levels are most likely underrepresented in the current transcriptome assembly. Secondly, mapping of the reads against multiple fish genomes is expected to efficiently exclude teleost sequences but it also removes parasite reads homologous to the host (Schulze et al., Reference Schulze, Schleicher, Guthke and Linde2016). This, however, reduces potential toxicity due to nucleotide level homology between the host and parasite and facilitates drug target discovery process. Thirdly, RNA sequencing has known bias against extremes of GC content (Swindell et al., Reference Swindell, Xing, Voorhees, Elder, Johnston and Gudjonsson2014; Sheng et al., Reference Sheng, Vickers, Zhao, Wang, Samuels, Koues, Shyr and Guo2016) and the AT-richness of the T. bryosalmonae transcripts may have resulted in an underrepresentation of the parasite reads in the raw data. This might be one potential reason why only ~10% protein sequences in the current T. bryosalmonae transcriptome assembly are complete (i.e. having start and stop codons at the opposite ends). In addition, our filtering pipeline most likely removed a proportion of the true parasite reads from the further assembly, reflecting the trade-off between false positive and negative identification of parasite reads in dual-RNA sequencing experiments. Finally, despite applying two different strategies for the transcriptome annotation, one fourth (24.6%) of the transcripts still remained uncharacterized. This demonstrates the urgent need for expanding myxozoan genomic and proteomic resources. Despite the limitations, this study illustrates the value and utility of our T. bryosalmonae transcriptome assembly by enabling the identification of novel drug targets providing new biological insights into both parasite biology and disease progression.
In conclusion, we developed an in silico pipeline to efficiently separate transcripts of T. bryosalmonae from the kidney tissue of a native fish host, followed by de novo transcriptome assembly. As such, this study represents an important milestone towards developing new effective control strategies against PKD. More generally, we expect that the transcriptome assembly of T. bryosalmonae represents a valuable resource for future studies on genomic architecture and evolution of the myxozoans and will facilitate research on molecular mechanisms underpinning PKD.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S003118202100010X.
Acknowledgements
The authors wish to acknowledge CSC-IT Center for Science, Finland, for generous computational resources. The authors are also thankful to Dr Jason W. Holland, Dr Ashlie Hartigan and Dr Gema Alama-Bermejo for providing transcriptome sequences for T. bryosalmonae, S. molnari and C. shasta, respectively.
Financial support
The work was supported by the University of Turku Foundation (F.A. grant for PhD studies), the Ella and Georg Ehrnrooth Foundation (F.A. grant for PhD studies), Estonian Ministry of Education and Research (R.G., institutional research funding project IUT8-2), Estonian Research Council (R.G., grant number PRG852) the Academy of Finland (A.V., grant number 266321), German Research Foundation research fellowship (P.V.D., grant number DE 2405/1-1) and the Swiss National Foundation (H.H., Sinergia Program grant CRSII3_147649).
Conflicts of interest
The authors declare there are no conflicts of interest.
Ethical standards
Collection of the studied fish was carried based on permits 46/2014 and 63/2015 issued by Estonian Ministry of Environment. The fish were euthanized in accordance to the principles described in ‘DIRECTIVE 2010/63/EU OF THE EUROPEAN PARLIAMENT AND OF THE COUNCIL of 22 September 2010 on the protection of animals used for scientific purpose’. The authors have followed the principles of the 3Rs (Replacement, Reduction and Refinement) and have involved the minimum number of animals to produce statistically reproducible results.
Open access
