Characterization of viral RNA splicing using whole-transcriptome datasets from host species

RNA alternative splicing (AS) is an important post-transcriptional mechanism enabling single genes to produce multiple proteins. It has been well demonstrated that viruses deploy host AS machinery for viral protein productions. However, knowledge on viral AS is limited to a few disease-causing viruses in model species. Here we report a novel approach to characterizing viral AS using whole transcriptome dataset from host species. Two insect transcriptomes (Acheta domesticus and Planococcus citri) generated in the 1,000 Insect Transcriptome Evolution (1KITE) project were used as a proof of concept using the new pipeline. Two closely related densoviruses (Acheta domesticus densovirus, AdDNV, and Planococcus citri densovirus, PcDNV, Ambidensovirus, Densovirinae, Parvoviridae) were detected and analyzed for AS patterns. The results suggested that although the two viruses shared major AS features, dramatic AS divergences were observed. Detailed analysis of the splicing junctions showed clusters of AS events occurred in two regions of the virus genome, demonstrating that transcriptome analysis could gain valuable insights into viral splicing. When applied to large-scale transcriptomics projects with diverse taxonomic sampling, our new method is expected to rapidly expand our knowledge on RNA splicing mechanisms for a wide range of viruses.

As increasing number of next-generation sequencing (NGS) datasets are being produced from various-omics initiatives, transcriptome sequencing of flora and fauna for a specific developmental stage/condition gains its popularity in biological research. Transcriptomics is implemented in discoveries of novel transcripts, SNPs, gene splicing and fusion, in determination of gene structure, function and regulation, and in quantification of expression levels 1 . It has already contributed a great deal of understanding to the mechanisms of functional elements, genes and transcripts 2,3 .
RNA splicing plays a vital role in genetics by increasing mRNA and protein diversities and by regulating gene expressions, providing an important link between genetic variation and disease [4][5][6][7] . Alternative splicing (AS) is one of the major mechanisms in increasing the diversity of proteins translated from a limited number of genes in metazoans 8,9 . The spliceosome complex, composed of at least 170 proteins and several small nuclear RNAs (snRNAs), is the key structure responsible for splicing in eukaryotes 10 . The complex defines exons/introns in transcribed RNAs by three major sequence elements: the 5′ splice site (donor site), the 3′ splice site (acceptor site), and the branch point 4,11 . When compared with annotated genome sequences, transcriptome sequencing could identify gene splicing isoforms and expression patterns associated with biological functions 12,13 .  After virus detection and false positive (sequence matched to non-viral subjects) removal, the assembled transcriptome sequences had best matches with AdDNV in A. domesticus and with PcDNV in P. citri (Fig. 1B, Supplementary Table S1 and Supplementary Text S1). Near full-length consensus genomes of AdDNV (NCBI Accession No: KX145610) and PcDNV (NCBI Accession No: KX145609) were called based on the templates of viral reference genome, with 5,259 bases (96.94% of the AdDNV reference KF015278.1 with 6,084 mapped reads) and 5,220 bases (97.03% of the PcDNV reference NC004289 with 33,604 mapped reads), respectively. The missing regions were located at 5′ and 3′ ends for both viruses and were replaced by Ns for the following analysis. Eleven single nucleotide variants (SNVs) were detected for AdDNV and eight SNVs for PcDNV ( Fig. 2A and B, Supplementary Table S2 and Supplementary Data file S1). None of these SNVs were located in splicing sites while some SNVs resulted in nonsynonymous mutations in ORF translations, i.e., six out of eleven SNVs in AdDNV and five out of eight in PcDNV (Supplementary Data file S1 Column N). Phylograms based on genome sequences and deduced proteins confirmed that AdDNV and PcDNV were closely related species within Ambidensovirus (Supplementary Fig. S1 and Supplementary Text S1) 43,49-57 . Expression of viral genes in transcriptome. Compared to the unspliced aligner BWA 58 (Supplementary   Table S2), TopHat2 59,60 obtained much greater depth coverages (Supplementary Table S3 and Supplementary Data file S2) due to successful alignments of fragmented (gapped) transcripts onto reference genomes. AdDNV_1KITE obtained 6,090 mapped reads (998,487 bases) with an average depth coverage of 189×, and a highest coverage of 659× on a single site in the NS region (Fig. 2C). The PcDNV_1KITE consensus sequence had 40,101 reads mapped against the genome at an average coverage of 1,325× and a highest coverage of 2,785× on a single site in the VP region (Fig. 2D). To examine whether the difference in sequencing depth may affect the viral gene expression patterns, we carried out downsampling analyses using proportions (1/10 and 1/20) of the PcDNV dataset (details provided in the Validation Section).
Several conserved protein domains play important roles in DNA replication, gene expression, infection and transfection in DNVs 46,61,62 . Conserved domains, including Parvo_NS1 located in NS proteins, Pavo_coat_N, Annotations of (C) AdDNV and (D) PcDNV. Coverage (Y-axis) of each nucleotide position (X-axis) was plotted for AdDNV_1KITE and PcDNV_1KITE. Six readingframes and previously described genes were represented using information provided by NCBI, including: start/stop codons (short blue/red vertical bars), transcription directions (black arrows, from top to bottom: forward reading frames +1, +2, +3 and reverse reading frames −1, −2, −3), ORFs (solid gray boxes). Virus genes (blue boxes), proteins (red boxes) and conserved motifs (black boxes) were represented according to the NCBI annotations. BWA mapping profiles (green lines), TopHat2 mapping profiles (purple lines), TopHat2 gap mapping profiles (yellow lines, the number of both splicing and non-splicing reads, correspond to splicing junctions) were represented according to the mapping results. AdDNV introns reported in existing studies include: In (nt 223 to  Denso_VP4 and phospholipase A2 (PLA2) motif located in VP proteins appeared to be highly expressed in both transcriptomes 63,64 (Fig. 2 and Supplementary Text S1). In the meantime, the unevenness of coverage also suggested possibilities of novel transcripts or other special characters in particular regions that were prone or hard to be enriched during the library construction.
Detected splicing junction regions matched to previously reported AdDNV introns In (nt 223 to 855), Ib (nt 4403 to 4544) and II (nt 4260 to 4434), which were determined by Sanger sequencing of RT-PCR products 48 48 . In short, two novel introns (AdDNV_ I2, _I3) and three known introns (AdDNV_I1, _I4, _I5) were detected while a previously reported intron (Ia, nt 4403 to 4758) was not detected in this research. Reads supported splicing junction models were showed in Fig. 3A by Integrative Genomics Viewer (IGV) 68 . Nine long open reading frames (ORFs) ranging from 207 to 2,451 nt were detected in the AdDNV_1KITE genome ( Table 2, Fig. 4), six of which had been previously validated by experiments 48 . Five of the ORFs did not require splicing whereas the other four were splicing products.
Introns of AdDNV_I1, _I4 and _I5 were supported by 88, 135 and 18 reads, respectively (Table 1,  Supplementary Table S4). These detected introns were congruent with experimental results from a previous study 48 , demonstrating that the NGS approach and bioinformatics applied in our study were reliable. The missing of a previously reported AdDNV (Intron Ia, nt 4403 to 4758) might be caused by inter-individual differences. In addition, two novel introns (AdDNV_I2 and AdDNV_I3) were detected at low read numbers (3 and 7, respectively, Table 1). Three out of four introns were validated using RT-PCR (details provided in the Validation Section). AdDNV_I2 was confirmed demonstrating that NGS was sensitive in detecting introns at very low expression levels. The splicing junctions and AS pattern of PcDNV are reported in the present study for the first time (Fig. 3B). The A5SS and A3SS introns occurred in PcDNV with high frequencies. Seven out of eight NS introns belonged to A5SS and all shared the same 3′ site at nt 879 (Fig. 3E). In the VP encoding region, four splicing islands were detected and all of them contained either A5SS or A3SS or both modes ( Fig. 3E and Table 1). Five in 23 introns had only a single read support suggesting limited function if there were any. All introns belonged to canonical GT-AG introns 69 . Twenty-eight ORFs ranging from 531 to 2,517 nt were detected in the PcDNV genome ( Table 2, Fig. 4). Five of the ORFs, as in AdDNV, did not require RNA splicing whereas the others (23 out of 28) were splicing products.
Different from AdDNV, PcDNV lacked the major AdDNV_I1 intron in the NS region. Instead, PcDNV displayed a set of seven A5SS introns with a shared receptor junction at nt 879. With different donor positions, the PcDNV NS introns produced a set of novel proteins composing NS ORF1 and NS ORF2 (Fig. 4C). The dominant NS splicing was PcDNV_NS_ORF7_I2. Other seven splicing events in NS region produced more additional isoforms of the ORF1-ORF2 protein and ORF1-ORF3 protein (Fig. 4C). 17 out of the 28 ORFs were located in VP region. The most common VP splicing in PcDNV was PcDNV_VP_ORF8_I16 (Fig. 4D). Like the previously described AdDNV_VP_ORF6_I4, PcDNV_VP_ORF8_I16 eliminated the stop codon of PcDNV_VP_ORF5 and joined PcDNV_VP_ORF4 (homologue of AdDNV-VP2) reading frame. Downsampling tests of PcDNV showed similar results (details in Validation Section), demonstrating that the differences observed between the two viruses were not solely caused by sequencing depth. More detailed descriptions of gene products of the two viruses could be found in Supplementary Text S1. Detected splicing junction models of PcDNV_1KITE: Solid gray areas represented the TopHat2 mapping profiles and each color-coded block represented a splicing junction. Red and purple blocks were forward and reverse junctions, respectively. The edge of each block represented the coverage of supporting reads and the length of a block represented the location of a splicing event. The number near each block was the coverage of supporting reads. The middle bridge showed the intron region from the splicing event. The block thickness represented frequency (the number of supporting reads) of the intron. Splice site compositions for donor sites, branch sites and acceptor sites of all GT-AG type introns in AdDNV_1KITE (Panel C) and PcDNV_1KITE (Panel D) were displayed using WebLogo. The overall height of each stack indicated the sequence conservation at that position, measured in bits. Proteins mediating the GT-AG splicing were labelled as snRNP (small nuclear ribonucleoproteins) and SR (splicing regulatory proteins). (E) Log2 scale of reads density of introns in the genome alignment: The Y-axis showed the expression levels (the number of reads) of intron related splicing events. Introns with forward junctions (red labels, at NS region) and reverse junctions (blue labels, at VP region) of PcDNV_1KITE (top half) and AdDNV_1KITE (bottom half) were shown in the genome alignment. Multiple splicing events (orange labels) were also displayed.  non-splicing ratios also produced transcripts with high FPKM/TPM values. In AdDNV, the ratio of AdDNV_I4 was the highest in AdDNV (Supplementary Table S4) and AdDNV_NS_ORF2 (Fig. 4A) had the highest FPKM indicating that the encoded AdDNV_NS1 was the most abundantly expressed protein ( Table 2) while AdDNV_ VP_ORF6_I4 was the mostly expressed VP isoform (VP1) (Fig. 4B).
Gene expressions of PcDNV_1KITE. PcDNV_I16 had the highest splicing/non-splicing ratio in PcDNV, followed by PcDNV_I23 (Supplementary Table S4). Based on the FPKM values, PcDNV_NS_ORF7_I2, PcDNV_ NS_ORF6_I1 and PcDNV_NS_ORF1 were the most abundantly expressed PcDNV_NS proteins (Fig. 4C, Table 2). Although PcDNV_VP_ORF8_I16 encoded a VP protein similar to AdDNV_VP_ORF6_I4 (AdDNV_ VP1), its FPKM value was second to that of PcDNV_VP_ORF10_I23 which encoded a novel protein without any conserved Ambiensovirus VP motifs 48,57,72 ( Fig. 4B and D, Table 2). The PcDNV_NS proteins had smaller FPKM values than the PcDNV_VP proteins, suggesting that PcDNV_NS proteins were expressed less abundantly than the PcDNV_VP proteins (Table 2). This was different to the situation found in AdDNV.

Expression patterns of viral genes. Expression patterns of the viral transcripts were significantly different between
AdDNV and PcDNV on two levels, i.e., transcript isoforms and expression abundance ( Fig. 4 and Table 2). Most of the spliced transcripts had low FPKM values or even effectively zero count (for those with small effective lengths), suggesting that these rare splicing products were unlikely to be responsible for any fundamental viral function 48,57,72 . On the other hand, as splicing products dominated in both NS and VP of AdDNV and PcDNV, RNA splicing played essential roles in these two densoviruses. Differential splicing resulted in remarkable divergence of the viral transcriptome (Figs 3 and 4). The two viruses were phylogenetically closely related ( Supplementary Fig. S1) therefore were expected to adopt similar gene expression strategies including AS 48,52 . Indeed, the splicing junctions were located in both NS and VP regions in AdDNV and PcDNV (Fig. 3A and B) and IR, A5SS and A3SS modes of alternative splicing were also detected for PcDNV in this study. The majority of the PcDNV GT-AG splicing occurred at a region similar to the AdDNV splicing hotspot (Fig. 3E), i.e., PcDNV_I16 (covered by 1,994 reads) and AdDNV_I4 (covered by 135 reads) were both positioned at the genome alignment region from nt 4,596 to 4,418. This structural consistency indicated that this splicing event was conserved across species and likely played an important role in densoviruses (Fig. 3E, Table 1, Supplementary Table S4).
On the other hand, novel splicing events, including the canonical GT-AG introns with SE modes (Table 1), were discovered for PcDNV. There were more splicing events in PcDNV than in AdDNV, even if the same level of sequencing depth was tested in each sample. The question whether or not such a difference may be caused by the number of host individuals sampled remains open. In general, the PcDNV consensus was less restricted than that of AdDNV, particularly at both donor (+4G and +5T) and acceptor positions (−3C and −5T) ( Fig. 3C and D, Supplementary Fig. S3, Supplementary Text S1). All five AdDNV introns belonged to the canonical intron GT-AG type 67 . It is worth noting that the consensus sequences of these GT-AG introns (nG|GTAnGTnG for donor and TnTTGCAG|An for acceptor, Fig. 3C) were different to the corresponding consensus sequences of GT-AG introns in PcDNV (AG|GTAAnnnn for donor and AnTTACAG|AT for acceptor with all junctions, AG|GTAAnnnn for donor and AATTACAG|AT for acceptor when junctions of very low frequencies (<3 reads) were excluded, Fig. 3D, Supplementary Fig. S3). Both viral consensus sequences had differences with the conserved consensus of U2 type GT-AG introns (AG|GTRAGT for donor and YYTTYYYYYYNCAG|G for acceptor) or U12 type introns (|RTATCCTTT for donor and TTCCTTRAY for branch sites) ( Fig. 3C and D, Supplementary Fig. S3, Supplementary Text S1) 69,73 . Validating splicing products using bioinformatics and experimental methods. Downsampling of the PcDNV dataset. To examine potential effects of sequencing depth on resultant splicing patterns, we randomly subsampled 1/10 and 1/20 of the PcDNV dataset and performed identical analyses on these subsamples. The 1/10 subsample had 2,067,041 reads in total, with 8,428 reads mapped to PcDNV after filtering (Supplementary Table S3). 12 out of 23 splicing junctions recovered from the full dataset were detected, including all major junctions (read depth > 50) and some (4 out of 14) minor junctions with low read depths (read depth < 10, Supplementary Tables S3 and S4). Similarly, the 1/20 subsample had 1,034,316 reads in total with 4,806 reads mapped to PcDNV after filtering. Although the mapped viral reads in the 1/20 subsample of PcDNV was fewer than that from AdDNV (6,090), more splicing junctions were detected (9 as oppose to 5 in AdDNV) ( Supplementary Tables S3 and S4). These results showed that increases of sequencing depth enhanced the number of rare splicing junctions but had no impact on the detection of major splicing junctions.
Splicing results from RNA-seq aligner STAR. We also analyzed splicing patterns using STAR 74 . In AdDNV, the five introns detected by TopHat2 were also supported by STAR (Supplementary Tables S3 and S4). In PcDNV, STAR detected 28 introns (23 by TopHat2), including 23 GT-AG introns that were also detected by TopHat2 and five additional rare splicing junctions supported by single read (Supplementary Table S4).
Junctions validation using RT-PCR method. RNA extracts remaining from the production of the two transcriptomes were used for RT-PCR validations (Supplementary Table S7). Primers were designed based on ORF sequences revealed by viral genome assemblies and annotation from this study to amplify regions of the splicing junctions. Several ORFs could share the same primers (Table 3, Supplementary Table S8, Supplementary Fig. 4). As expected, RT-PCR amplified ORFs with the most abundant junctions or highest expression levels sharing the same primers. AdDNV_I2, _I4 and _I5 (read depths of 3, 135 and 18, respectively) and PcDNV_I1, _I2, _I16 and _I23 (read depths of 43,275,1994 and 1699, respectively) were validated by Sanger sequencing of the RT-PCR products 75 . However, some rare junctions with low coverage depths were not confirmed by this method. The results showed that deep sequencing was more sensitive in detecting rare junctions than the RT-PCR based approach. On the other hand, it is also possible that splicing junctions supported by a single read might be caused by sequencing and/or mapping errors.

Discussion
Virus infections are common in most eukaryotic organisms 30 . Deep sequencing of transcriptomes coupled with bioinformatics pipeline developed in the present study can readily detect transcripts of the target organism (host) as well as those of the pathogens. In our pipeline, we firstly screened the assembled sequences (scaffolds) for potential viral sequences against a customized viral database instead of the full Nt database. The resulting viral candidate sequences were then screened against the Nt database to remove false positives (sequences matched to non-viral subjects). The outcome of this two-steps virus screening is the same as directly using the assembled sequences to screen against the Nt database. However, virus screening using two large datasets (assembled sequences and Nt database) consumes much greater computer resource. Therefore, our pipeline for virus detection improved computational efficiency without compromising on accuracy. Given that a wide range of transcriptomes for non-model organisms have been de novo assembled by a series of large-scale transcriptome projects, including the 1KITE 42,76,77 , this practice provides an effective pathway to characterize virus diversity across various lineages of life.
In addition, our study demonstrated that transcriptome sequencing was an effective and accurate approach to improve the understandings of gene catalogues, expression levels, and RNA splicing patterns for pathogens. Viral gene expression profiles can be identified from the transcriptomes (Figs 3 and 4). Transcriptomic analysis in a phylogenetic context can help to elucidate functional conservativeness and novelty of expressed genes. For instance, phylogenetically closely related viruses, AdDNV and PcDNV, exhibited both conserved transcriptions and lineage-specific profiles. Common VP1 isoforms were abundantly transcribed using highly conserved splicing events (AdDNV_VP_ORF6_I4 and PcDNV_VP_ORF6_I16, Fig. 4), suggesting that these proteins were essential for the virus survival. On the other hand, differences in structural and non-structural proteins inferred from viral transcripts were found between the two viruses, including novel sequence divergences and distinctive variations in expression ratios and levels.
The diverse splicing mechanisms in PcDNV seemed to have directly led to more varieties of gene splicings than that in AdDNV (Fig. 3E, Table 1, Supplementary Tables S4 and S5) and these novel introns inturn produced new protein isoforms. The number of pooled individual specimens in the 1KITE project was justified according to the body weight of insects to provide sufficient RNA for transcriptome sequencing, and it could cause limitations in sample variation and inadequacies in detecting of inter-individual differences in the current study. Variations introduced by sequencing depth (1,325× versus 189×), number of pooled individual specimens (150 versus 1), temporal variation in gene expressions and inter-individual difference may have also contributed to some of the unique patterns observed in PcDNV, especially some of the rare introns. In future research, analyses on more virus-host pairs with controlled individual numbers could help to clarify the effects of different factors. RT-PCR based validations verified the majority of major novel splicing junctions, but most junctions with low read support could not be confirmed. Therefore, the possibility of artificial errors cannot be completely ruled out for expressed junctions at very low read coverage. Empirical validation of rare splicing junctions is a challenging task. Although the Sanger sequencing of RT-PCR products has successfully verified a number of major and minor junctions, the conventional method is far not as sensitive as NGS. Therefore, interpretation of rare slicing junctions should be cautious until additional technologies are developed for validation and functional testing.
Nevertheless, transcriptome sequencing of host and co-expressed pathogen creates a unique opportunity to examine host-virus association. To our knowledge, our report on AS consensus sequences were the first description for PcDNV and AdDNV. Viruses rely on the host machinery for RNA biology and can co-evolve with the host splicing 78 . It is likely that the observed divergence of viral RNA splicing patterns (Fig. 3) were influenced by both host and viral factors. Additional pipelines for analyzing the host splicing patterns may help us to understand the virus-host interaction and co-evolution in the future.

Methods
Virus detection. Raw reads were assembled using SOAPdenovo-Trans 79 with the following settings: "-K 31 -i 20 -e 3 -M 3 -L 100" 42 . We then searched for matches against a customized virus database (described below) using the Basic Local Alignment Search Tool (BLASTn, version 2.2.26) 80 , including 6 steps (Fig. 1B).
Customized virus database. All known viruses in non-redundant nucleotide (Nt) databases and their corresponding taxonomic identity were downloaded from the GenBank (November 2014) and served as the virus reference database, which was much smaller than the complete Nt database but more dedicated to viruses. This modification improved screening efficiency and alleviated computational demand.
Virus search. Assembled transcriptome scaffolds were searched for sequence homology against the customized virus reference database using a local BLASTN algorithm (e value < 1e-5). Only query sequences with a match-length ≥200 nt and identity ≥90% were retained as candidate viral sequences for further analyses.
Overlapping candidate viral sequences were merged by combining the BLASTN results using an in house Perl script (available at https://github.com/linzhi2013/Virusfishing) with improved procedures in query selection (merged candidate viral sequences versus whole scaffolds or sequencing reads) 81 .
Removal of false positives. Candidate viral sequences were compared to the complete Nt database using BLASTN aiming to identify false positives. Sequences with higher scores to non-virus subjects were deemed as false positives and subsequently removed from downstream analyses. align all raw reads of the transcriptome onto virus reference genomes with default parameters. After removal of PCR redundancy using SAMtools (If multiple read pairs have identical external coordinates, only the pair with highest mapping quality was retained) and reads with more than two mismatches, consensus viral sequences were obtained using SAMtools (version 0. 1.19) 82 , with ambiguous sites substituted by the base of highest allele frequency and with missing sites replaced by Ns.
Genome coverage and annotation. ReSeqTools (version 0.23) 83 and self-developed Perl scripts were used to calculate the genome coverage and read depth. Viral gene information was downloaded from the GenBank. All translation start site and translation termination site at both directions were annotated based on prediction of amino acid sequences.
Phylogram construction. The virus phylogram was inferred using both genome sequences and protein sequences of the assembled viruses, reference viruses and other nine related viruses from the same family downloaded from GenBank. MUSCLE 84 was applied to conduct multiple sequence alignments (MSA) and MEGA5 was used for tree construction using the maximum likelihood method 85 and Neighbour-Joining method 86 with 1000 bootstraps (MEGA5 software) 87 .

Gene expression analyses. Gene expressions were analysed in five steps:
Detection of alternative splicing. TopHat2 (version tophat-2.0.7) 59,88 was applied to indentify RNA splicing patterns based on junction signals. Alternative splicing events were identified accroding to gene splicing patterns. All transcriptome reads were mapped to reference viral genome assemblies and host database using paramters "-r 10 -i 50 -I 2000-library-type fr-unstranded -G" according to the TopHat2 manual, respectively. TopHat2 aligns reads that are spanning across gaps onto a reference more efficiently than the unspliced aligners, such as BWA and Bowtie 60,89 . RNA mapping characteristics were obtained from high-quality BAM files with unique mapping reads after removing of PCR redundancy and reads with more than two mismatches.

Determination of AS junctions, introns and gene expression patterns.
Splicing junction locations and counts were obtained from spliced sequences that might contain one or two junctions per read. Non-canonical splicing sites with low number of supporting (n < 3) reads were filtered out from further analysis. And the putative branch site regions were identified in introns by searching regions relative to the 3′ splicing site 67 , using an in house Perl script. Sequence logos of all detected junctions and junctions with coverage depth ≥ 3 were generated by WebLogo 90 and consensus sequences were generated by the tool cons from EMBOSS 91 , respectively. Introns and splicing patterns were identified according to the letters on donor and acceptor sites of the splice junctions 69,73,92,93 . The splicing level of an intron was recorded as read counts mapped to a junction, similar to the method by MATS 94 in calculating the exon inclusion level. Integrative Genomics Viewer (IGV) 68,95 was used to show the putative viral AS and intron models.
Annotation of viral open reading frames. Firstly, the putative coding regions of the viral genome were identified using an in house Perl script. Putative expression products with amino acid length > 30 aa were defined according to the open reading frames (ORFs) from the unspliced/spliced sequences. The resulting amino acid sequences were annotated using the online searching tool BLASTp at NCBI 80 . ORFs and spliced isoforms were named following the names described in the NCBI annotation with splicing sites noted as: I for intron and MI for multiple introns of a splicing event.
Comparison between related viral species. Multiple aligments were conducted for genome sequences using ClustalW 96 to compare splicing junction locations between viral species. The graphical mapping details of final alignments were drawn using R.
Calculation of expression levels. A local estimation of each junctions were used to show their expression levels. The number of splicing and non-splicing reads located in junction regions were caculated based on the alignment Bam files. The number of non-splicing reads with donor or acceptor sites (reads spanning the exon-intron) were also caculated. We considered one ORF as one gene to calculate the expression levels. The expression level related values: i.e., the counts of positions of valid fragment (effective length), the sum of the posterior probability of each read coming from this part over all reads (expected count), the transcripts per million (TPM), and fragments per kilobase of transcript per million mapped reads (FPKM), were calculated using RSEM 70 and Kallisto 71 . A relative expression level of each viral gene was measured using the FPKM value of one gene divided by the FPKM value of the highest expressed gene of the same virus.
Validation of splicing products. Downsampling. 1/10 and 1/20 subsamples of PcDNV related whole transcriptome reads were randomly selected from total reads using the toolkit Seqtk 97 and analyzed using the same analytical pipeline conducted for the full transcriptomes to examine potential impact of sequencing depth.
Junction detection using STAR. STAR (version 2.5.3a) was applied to check the detected junctions using the parameters "-runMode genomeGenerate-sjdbOverhang 149-genomeChrBinNbits 12-genomeSAindexNbases 7 -sjdbGTFfile …" for STAR index construction and "-alignSJoverhangMin 5-alignIntronMin 50-alignIntron-Max 2000-outSJfilterOverhangMin 5 5 5 5" for alignment. Then the alignment BAM file was used to define junctions using the script "Splicing_Search.pl" in pipeline after filtering with the same parameters as described in the TopHat2 step. DNaseI treatment and reverse transcription. Total mRNAs were obtained from the 1KITE project directly and the quality was evaluated using a Qubit Fluorometer. To eliminate DNA contamination, 1 μg of total RNA was treated with 1 U DNaseI (Promega, US) and incubated at 37 °C for 30 min. Then samples (1 μg) of total RNA were reversely transcribed in a 25 μL reaction mixtures with oligo(dT) 15 primers and M-MLV reverse transcriptase (Promega, US).
PCR. Splicing was validated by PCR of the cDNAs using primers designed from viral genome sequences from this study (Supplementary Table S8). The PCR mixture (50 μL) consisted of 4 μL cDNA template, 0.2 μM primers, 200 μM dNTP mix and 2.5 U TaKaRa LA Taq polymerase. The amplification conditions were denaturing at 94 °C for 5 min, then 40 cycles of denaturing at 94 °C for 30 s, annealing at 50 °C for 30 s, and extension at 72 °C (the extension time depends on the length of PCR amplicon, 1 kb/min), and final extension at 72 °C for 10 min.