All procedures performed on animals in this study were conducted following the animal husbandry guidelines of the Chinese Academy of Medical Sciences and with permission from the Experimental Animal Committee of the Chinese Academy of Medical Sciences with the Ethical Clearance Number IPB-2011-8.

All parasite materials were collected from the Xinjiang Uyghur Autonomous Region, China. For E. granulosus transcriptome sequencing, a unilocular cyst was isolated from a sheep liver infected with E. granulosus. After rinsing 10 times with normal saline, protoscoleces were aspirated from the hydatid cysts, washed three times with normal saline and then soaked in RNAlater solution (Ambion, CA, United States); samples were stored at -80°C until total RNA was isolated. Hydatid cysts of E. multilocularis were freshly isolated from male gerbils at 5–6 months after infection with E. multilocularis. After rinsing with normal saline, the cysts were opened and homogenized to a smooth consistency in 200 ml normal saline. The suspension was successively passed through 40, 80, 140, and 325 mesh metal sieves. The protoscoleces on the 325 mesh metal sieve were collected into a beaker and were washed five times with ice-cold normal saline. The storage method was as above.

For mRNA library construction and deep sequencing, total RNA was isolated using an RNeasy Mini kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Potential contaminating genomic DNA was removed from RNA samples using a Turbo DNA-free kit (Ambion, CA, United States). The quantity and integrity of the total RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, United States). Normalized starting quantities of total RNA were used to prepare two biological replicates of Illumina sequencing libraries for each species using a TruSeq^TM RNA sample preparation kit (Illumina, San Diego, CA, United States). Following the manufacturer’s instructions, libraries were constructed and sequenced using a paired-end read protocol with 100 bp of data collected per run on an Illumina HiSeq^TM 2000 instrument housed at MininGene Biotechnology Co. Ltd, Beijing, China. Data analyses and base calling were performed using the Illumina instrument software.

A Perl program was written to remove low-quality reads (more than half of the bases had a quality value less than 5). The retained high-quality reads were mapped to the published genomes of E. granulosus (Tsai et al., 2013; Zheng et al., 2013) and E. multilocularis (Tsai et al., 2013) by HISAT (Kim et al., 2015), and then assembled with StringTie (Pertea et al., 2015) to construct unique transcript sequences, using the parameters: -g -b -u -o. AS events in a gene were detected by SUPPA (Alamancos et al., 2015). Cuffquant (Trapnell et al., 2010) was used to quantify the expression of genes and transcripts, and the raw expression values were normalized by Cuffnorm. The expression values were calculated for each sample based on the number of fragments per kilobase of exon per million reads mapped. InterPro domain information (Mulder et al., 2003) was annotated by InterProScan (release 53.0) (Zdobnov and Apweiler, 2001) and functional assignments mapped with GO (Harris et al., 2004). WEGO (v3.3) (Ye et al., 2006) was used for GO classification. To obtain an overview of the gene pathway networks, KEGG analysis was performed by comparing the transcripts with the KEGG database (release 58) (Kanehisa et al., 2006) using BLASTX (Altschul et al., 1997; Ouyang et al., 2007) at E-values ≤ 1E^-10. A Perl script was used to retrieve KO information from the BLAST results, and then establish pathway associations between the transcripts and database. A Perl script was written to count the number of total genes and the number of AS genes in a certain functional pathway, then a significant score P-value was calculated using the R scripts under hypergeometric distribution. A function pathway with P-value < 0.05 was deemed as enrichment in AS genes.

To validate AS events experimentally, RT-PCR was performed for eight pairs of homologous genes that underwent different AS types in the two parasites. For each sample, 1 μg total RNA was reverse transcribed into first-strand cDNA using a SuperScript III Reverse Transcriptase kit (Invitrogen, Carlsbad, CA, United States) with oligo(dT)_12-18 primer. Primers were designed based on the consensus gene sequences of the two parasites, and the primer sequences are listed in Supplementary Table S1. PCR was performed in a 25 μl reaction system using High-Fidelity PCR Master Mix (NEB) and the procedure was as follows: initial denaturation at 98°C for 1 min; 98°C for 10 s, 60°C for 30 s, and 72°C for 30 s, for 35 cycles; and a final extension at 72°C for 5 min. The PCR products were visualized by 1.5% agarose gel electrophoresis analysis.

Using protein sequences of Schistosoma mansoni VALs (Chalmers et al., 2008) as query sequences, the initial protein sequences of Echinococcus VALs were retrieved from the non-redundant protein sequence database of the National Center for Biotechnology Information by BLASTP search (E-value cut-off: 10^-5). CD-HIT v4.5.4 software¹ was then used to eliminate redundant protein sequences. All remaining protein sequences were further examined for the presence of SCP/TAPS-representative protein domains by searching the Conserved Domain Database (Marchler-Bauer et al., 2011) of the National Center for Biotechnology Information. A phylogenetic relationship tree was built using the full-length amino acid sequences of Echinococcus VALs in the following steps. First, sequences were aligned using ClustalX (Larkin et al., 2007), then they were refined manually, and a phylogenetic tree was finally generated using MEGA 5.0 software (Tamura et al., 2011) by the neighbor-joining method (the bootstrap test was performed with 1000 replicates).

In this study, a total of 43,509,848 and 44,379,588 high-quality reads (100 bp/read) were obtained from worms of E. granulosus and E. multilocularis, respectively. Low-quality sequences were removed by a Perl program, and the retained high-quality reads were mapped to the published genomes of E. granulosus and E. multilocularis to construct unique transcript sequences. A total of 34,717,856 reads (79.79%) mapped to the E. granulosus genome reported by Tsai et al. (2013), 31,156,736 reads (71.61%) mapped to the E. granulosus genome reported by Zheng et al. (2013), and 38,882,179 reads (87.61%) mapped to the E. multilocularis genome (Tsai et al., 2013) (Table 1). Therefore, all the subsequent assemblies of reads and bioinformatics analyses, such as gene annotation and AS events, were based on the first mapping results for E. granulosus. Using the software StringTie, 24,550 (7,925 known and 16,625 novel transcripts) and 23,771 transcripts (8,432 known and 15,339 novel transcripts) were assembled for E. granulosus and E. multilocularis, respectively, and the assembly yielded 11,330 genes (6,815 known and 4,515 novel genes) for E. granulosus and 10,101 genes (7,051 known and 3,050 novel genes) for E. multilocularis (Table 1 and Supplementary Tables S2, S3), compared with the reference genome data. The mean length of transcripts for the two parasites was 2016 and 2086 bp, respectively (Table 1). All sequence data has been deposited in the GEO database² with an accession number of GSE59173.

The sequence reads could be classified into four types: exon, intron, intergenic, and spliced. The proportions of the four sequence types in the E. granulosus transcriptome were 56% (exons), 4% (introns), 15% (intergenic), and 25% (spliced) (Figure 1). The percentages of exons and introns in the E. multilocularis transcriptome were 55 and 3%, which is similar to the E. granulosus transcriptome; while the proportions of intergenic and spliced transcripts in the E. multilocularis transcriptome were 24 and 18%, which is significantly different from the proportions in the E. granulosus transcriptome (Figure 1). In our transcriptome database, the total exon lengths for E. granulosus and E. multilocularis were 22,249,494 bp (19.80% of the genome) and 21,805,251 bp (19.17% of the genome), which were longer than the total exon lengths in the reference genome databases for E. granulosus (15,213,358 bp, 13.54%) and E. multilocularis (16,112,293 bp, 14.16%) (Tsai et al., 2013).

A total of 6,826 AS events from 3,774 E. granulosus genes (33.31%) and 6,644 AS events in 3,611 E. multilocularis genes (35.75%) were bioinformatically predicted (Figure 2A and Supplementary Tables S4, S5). These AS events could be classified into seven types: intron retention, exon skipping, alternative donor site, alternative acceptor site, alternative first exon, alternative last exon, and mutually exclusive exon. The results showed that the proportions of the seven AS types were similar in E. granulosus and E. multilocularis. Intron retention (39 and 42%) was the most common AS type, followed by alternative donor site (21 and 17%), alternative acceptor site (17 and 18%), exon skipping (both 16%), and the least common splicing forms were alternative first exon (both 4%), alternative last exon (both 2%), and mutually exclusive exon (both 1%) in the two parasites, E. granulosus and E. multilocularis, respectively (Figure 2B and Supplementary Tables S4, S5).

To validate the AS events predicted by bioinformatic analysis, RT-PCR was carried out on eight pairs of homologous genes producing different types of AS events in the two parasites. Four pairs of homologous genes with intron-retention events were serine/arginine-rich splicing factor (EG.4660 and EM.3827), SNARE protein Sec22 (EG.5720 and EM.2557), hypothetical protein (EG.3709 and EM.5622), and hypothetical protein (EG.5720 and EM.2557). The other four pairs of homologous genes with exon skipping events were peroxidasin (EG.9768 and EM.6188), filamin (EG.8983 and EM.6048), adducin-related protein (EG.7075 and EM.1228), and heat shock protein 70 (EG.583 and EM.6781). Agarose gel electrophoresis showed that bands corresponding with the predicted AS events could be detected in the RT-PCR products (Figures 3, 4). Notably, some extra bands detected in the RT-PCR products were proved to be additional AS events from the corresponding genes by Sanger sequencing method (Figures 3, 4). The results were indicative that the bioinformatic prediction based on the primary sequencing data was reliable, and additional AS events for some genes, which could not be sequenced nor predicted in the present study, also occurred in the two parasites.

Gene Ontology analysis was performed to summarize and explore the functional categories of the genes that underwent AS events in the two parasites. A total of 2,364 AS genes of E. granulosus and 2,376 AS genes of E. multilocularis were annotated for different GO terms (Supplementary Tables S3, S4). The results showed that the proportions of enriched GO terms for AS genes were globally similar between E. granulosus and E. multilocularis based on biological process (Figure 5A), molecular function (Figure 5B), and cellular component (Figure 5C). For biological processes, genes with AS events were predominantly enriched in GO terms that were relevant to cellular process and metabolic process. For the molecular function category, the majority of AS genes were annotated with the term binding, followed by the term catalytic activity.

Moreover, 1,099 and 1,115 AS genes were annotated for KEGG pathways. The enrichment pathway analysis showed that these genes were significantly enriched in 26 pathways relating to cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems (P < 0.05, Table 2). For environmental information processing, the AS genes were mainly enriched in the phosphatidylinositol signaling system, cellular antigens, and Jak-STAT, VEGF, and Notch signaling pathways, suggesting that AS may play important roles in the response to host environmental stress. For metabolism, the most significantly enriched pathways were involved in purine, fatty acid, galactose, glycerolipid, or N-glycan metabolism or biosynthesis. For genetic information processing, the enriched AS genes were grouped mainly in basal transcription factors, RNA transport, and the mRNA surveillance pathway. For organismal systems, the significant enrichment of AS genes included those relating to development (osteoclast differentiation), the immune system (Fc epsilon RI, T cell receptor, and Toll-like receptor signaling pathway), and the endocrine system (GnRH signaling pathway) (Table 2). Further experimental verification is necessary to clarify the biological regulation mechanisms of these AS events in Echinococcus species.

Using AS gene sequences of E. multilocularis as the query database, we performed a BLASTX search to retrieve homologous genes in the following five organisms: E. granulosus, S. japonicum, S. mansoni, Caenorhabditis elegan, and Homo sapiens. A total of 3389 of E. multilocularis AS genes were identified as having homologous genes in the five organisms. Venn diagram analysis showed that there were three major groups of E. multilocularis AS genes: group I (conserved genes, 158 genes) (Figure 6A and Supplementary Table S6), for which homologous genes could be retrieved in all five organisms; group II (Platyhelminthes genes, 300 genes) (Figure 6A and Supplementary Table S7), for which homologous genes could only be retrieved in E. granulosus, S. japonicum, and S. mansoni; and group III (Cestoidea genes, 2360 genes) (Figure 6A and Supplementary Table S8), for which homologous genes could only be retrieved in E. granulosus. Moreover, the percentages of E. granulosus genes with AS events from the three groups were compared to investigate whether AS events were related to their gene sequence conservation. In group I, 63.3% of E. granulosus homologs with AS events were identified, which was lower than the percentages for group II (66.7%), and group III (67.8%) (Figure 6B). Meanwhile, we compared the proportions of S. japonicum genes with AS events from group I and group II using our published S. japonicum transcriptome data (Piao et al., 2014). The percentage of S. japonicum homologs with AS events in group I (67.1%) was lower than the percentage in group II (73.0%) (Figure 6C). The representative homologous gene pairs with AS events from the three groups are listed in Table 3.

The VAL genes are members of the sperm-coating protein/Tpx-1/Ag5/PR-1/Sc7 (SCP/TAPS) superfamily, which may be important for chronic host/parasite interactions (Chalmers et al., 2008; Farias et al., 2012). It has been found that members of the VAL family exhibit multiple expression patterns and AS in the life cycle of S. mansoni (Chalmers et al., 2008) and S. japonicum (Liu et al., 2015). Based on the genome databases of the two species, we identified 22 E. granulosus VALs and 19 E. multilocularis VALs with complete SCP/TAPS domains. Next, a phylogenetic tree was reconstructed to examine the phylogenetic relationships among VALs in the two Echinococcus species. The result showed that the Echinococcus VAL family was divided into three major groups, and the majority of genes had homologous counterparts between the two Echinococcus species (Figure 7). Using these nucleotide sequences of Echinococcus VALs as query sequences, local BLASTN searches were performed to identify the corresponding mRNA sequences in our RNA-Seq data. Nine E. granulosus VALs and nine E. multilocularis VALs were proved to be transcribed in protoscoleces (Figure 7), and the others may be expressed in other developmental stages or the amount of their transcription was too low to be detected by our study. Interestingly, AS transcripts of E. granulosus VALs (CDS24051.1 and CDS22593.1) and the E. multilocularis VAL (CDS35801.1) from group II were discovered in our RNA-Seq data (Figure 7). Therefore, AS events may be ubiquitous in the VAL family of parasitic platyhelminths, and is well worth further study.