To perform an extensive comparison of the methods on a wide array of plant viruses with different genome types, five plant samples were selected, containing different viruses, most of which (samples I, II, and IV) have already been very well characterized for viruses using HTS and targeted detection methods (Pecman et al., 2017). Either infected plant leaf samples (samples I, II, III, and V) or an infected seed sample (sample IV) were used (Table 1; Supplementary Tables S1, S2). Additionally, leaves of healthy tobacco plants were used as a negative control. RNA was isolated (Figure 1) from all leaf samples using the RNeasy Plant Mini Kit (Qiagen) including a DNase step (RNase-Free DNase Set, Qiagen) according to the manufacturer's instructions. From the seed sample (sample IV), RNA was isolated using a combination of CTAB buffer and RNeasy Plant Mini Kit (Qiagen) as described in the study of Adams et al. (2009) with minor modification: incubation with 4 M LiCl was done at 4°C overnight.

A part of the RNA extracted from sample I was further processed: ribosomal RNA was depleted from the extract using the RiboMinusTM Plant Kit for RNA-seq (Invitrogen# A1083808), obtaining two versions of sample I: with and without ribosomal RNA depletion (Figure 1, for more details refer to Supplementary Tables S1, S2).

The direct RNA sequencing kit (SQK-RNA002) and cDNA-PCR sequencing kit (SQK-PCS108) required polyA+ RNA as input RNA. Therefore, for each sample, polyA tailing of RNA was performed using E. coli poly(A) polymerase (NEB# M0276). The mixture was incubated at 37°C for 15 min. The reaction was stopped by directly proceeding to the clean-up step with Agencourt® AMPure® XP beads (Beckman Coulter) using 1.8 (AMPure® XP beads): 1 (poly(A) tailing mixture) ratio.

The direct RNA sequencing kit (SQK-RNA002, version DRS_9080_v2_revB_22Nov2018, Oxford Nanopore Technologies) was used to prepare sequencing libraries for all the samples included in the study. The cDNA-PCR sequencing kit (SQK-PCS108, version PCS_9035_v108_revH_26Jun2017, Oxford Nanopore Technologies) was additionally used to prepare libraries from sample I and sample I with ribosomal-depleted RNA.

For all but one library, the recommended amount of RNA input was used (Oxford Nanopore Technologies protocols); however, when preparing the library for direct RNA sequencing of rRNA-depleted totRNA (sample I), 278 ng of RNA, instead of 500 ng of RNA, was used (due to the lower extraction yield). Each library was then sequenced on a separate flow cell (R9.5.1) for 46–48 h using MinION device and MinKNOW software (v18.12.6). The reads were base-called using Guppy v3.1.5 and command: rna_r9.4.1_70bps_hac.cfg / dna_r9.4.1_bps450_hac.cfg –device auto –u_substitution false.

Sequencing libraries for each sample were prepared using total RNA and the ScriptSeq™ Complete Plant Leaf Kit (production discontinued, Illumina, USA) which included a ribosomal RNA depletion step. The libraries were sequenced on a MiSeq (Illumina, USA) using a 2x300-bp (V3) cartridge.

In the first part of the data analysis, the aim was to investigate how well different sequencing approaches compare in terms of virus detection from complete datasets (Figure 1). To achieve this, established in-house virus detection workflows were used and virus presence was reported as follows.

Illumina reads were analyzed using a pipeline in CLC Genomic Workbench (v12, v21) and additional Diamond BLASTX analysis, as described below. Quality control was performed, then, adapters were removed from all reads, and additionally, reads were trimmed by quality (quality limit = 0.05) and length (all reads shorter than 30 nts were discarded). Trimmed reads were mapped to viral RefSeq (NCBI database, updated 19.05. 2019) and de novo assembled. Contigs (longer than 100 nts) were mapped to viral RefSeq (NCBI database, updated 19.05.2019) and unmapped contigs were further analyzed by searching for conserved protein domains using Pfam analysis (v32) (refer to Supplementary Tables S3–S6 for details about the used parameters). Additionally, de novo assembled contigs were analyzed with Diamond BLASTX (v0.9.22) (Buchfink et al., 2015) against the NCBI nr database (June 2018). Diamond outputs were taxonomically classified and visualized using Megan 6.19.2 (Huson et al., 2007).

To analyse nanopore sequencing data, a similar pipeline was constructed using tools, which enable analysis of long reads. The statistics and quality of MinION reads were checked using the programs NanoQC v0.8.1, NanoStat v1.1.2, and NanoPlot v1.20.1 (De Coster et al., 2018). The read plots generated using the NanoQC v0.8.1 were inspected for each sample individually and were used to determine how to trim them (length of reads, head of reads, and tail of reads) (Supplementary Table S1) using program NanoFilt 2.5.0 (De Coster et al., 2018). The trimmed reads were again quality checked with NanoStat v1.1.2 and then mapped to the viral RefSeq (NCBI database, updated 19.05.2019) using minimap2 (v2.16-r922) and the commands: minimap2 -ax splice -uf -k14 for direct RNA reads and minimap2 -ax map-ont for cDNA-PCR reads. Reads were also analyzed using Diamond BLASTX (v0.9.22) (Buchfink et al., 2015) with the –frameshift 15–range-culling–sensitive command option. All reads were de novo assembled by combining different programs using a Pomoxis (https://github.com/nanoporetech/pomoxis) inspired approach: after fast mapping (minimap2) (Li, 2018) and de novo assembly (Miniasm) (Li, 2016), two rounds of the contig correction using racon (Vaser et al., 2017) were run. The script together with corresponding parameters is shown in Supplementary Data 1. The assembled contigs were analyzed using BLASTn against the NCBI nt database and visualized with Megan 6.19.2 (Huson et al., 2007). TASVd (sample V) was not detected by mapping direct RNA MinION reads to the viral RefSeq database, so in the next step, more closely related sequence from NCBI GenBank (KY810784) was used as the reference for reads and contig mapping.

To be able to perform comparisons of different sequencing approaches for different viruses or viroids, complete or near complete consensus genomic sequences of viruses in the samples were obtained. Since some of the samples were identical to the ones from a previous study (Pecman et al., 2017), these reference sequences were already available. For phaseolus vulgaris alphaendornavirus 2 (PVeV2, sample V) and tomato spotted wilt orthotospovirus TSWV (sample III), a reference consensus was generated based on the mapping of Illumina reads and contigs, as previously described (Pecman et al., 2017). For pepino mosaic virus (PepMV), two divergent strains (80% nucleotide identity) were detected in the sample I (PepMV-EU and PepMV-Ch2); thus, in this case, the complete genome sequences of KF718832.1 and JX866666.1 were used in subsequent comparisons as described previously (Pecman et al., 2017). For tomato yellow leaf curl virus (TYLCV) (KY810789), tomato chlorosis virus (ToCV) (KY810786), (KY810787), cauliflower mosaic virus (CaMV) (KY810770), and cabbage cytorhabdovirus 1 (CCyV1) (KY810772), the reference sequences described in the study of Pecman et al., (2017) were used in first step, but due to some mismatches after mapping Illumina reads and contigs to those reference sequences, few nucleotides were changed and “new” consensus genome sequences were used for the purpose of the following analysis only (Supplementary Data 2). In the case of sample V, only PVeV2 complete genome was covered by reads by both approaches, and thus, only this endornavirus was included in further analyses.

To calculate the viral sequence content in the complete datasets, trimmed reads for each dataset were mapped to the corresponding reconstructed or selected viral or viroid reference sequences (Section Construction or Selection of Consensus Genomes Used in Further Analyses). For each sample–virus–sequencing type combination, three parameters were reported: the percentage of mapped reads, average depth (the number of times a nucleotide is covered by a sequencing read averaged across the complete reference genome sequence), and fraction of reference covered by reads.

The next step enabled the investigation of (I) average nanopore reads and contig identities (compared to the corresponding reference sequences) and (II) identities of consensus genome sequences generated after mapping the reads to the reference sequence, from now on named “consensus sequence identity”. (I) To determine the average identities of nanopore sequencing reads (proxy of error rate), and the identities of de novo assembled contigs (generated from those reads), compared to reference sequences, reads and contigs were mapped to corresponding viral genome reference sequences. Identities were calculated using Minimap2 PAF output (Pairwise mApping Format) (Li, 2018). In this way, the average BLAST-like alignment identity was calculated for each mapping of nanopore sequencing data, either with reads or with contigs.

(II) If sequencing errors in nanopore reads are relatively random and if there is a substantial number of reads covering the reference genome sequence after mapping, the resulting consensus sequence should be “error corrected”. To test this, extracted consensus sequences (or their fragments–if whole genomes were not obtained) derived from mapping nanopore or Illumina reads to reference sequences were aligned with original reference sequences and pairwise identities were calculated using CLC Workbench Genomics v12, v21. For this analysis, consensus sequences (or their fragments if whole genomes were not obtained) were used if they had at least 1x average coverage in the read mapping step.

To be able to compare the datasets obtained by different sequencing approaches, complete datasets were rarefied to obtain subsamples with comparable numbers of nucleotides. Depending on the original size (Supplementary Tables S1, S2), datasets from different sequencing approaches generated a different number of subsamples. The largest datasets, contained 1,500 million nucleotides, followed by 1,300, 1,100, 900, 700, 500, 300, 200, 100, 50, 30, and 10 millions of nucleotides (Supplementary Table S7). Each set of subsampled reads was randomly generated using CLC Genomic Workbench (v12, v21) for Illumina reads and Seqtk Sample for MinION reads. Subsampling was repeated five times for each of the subsample sizes for each sample.

The rarefied subsets were used for further comparative analysis of different sequencing approaches. The read subsamples were mapped to reconstructed or selected reference sequences. The reads were also de novo assembled and the resulting contigs were mapped to the selected reference sequences (for nanopore approaches as described above, for Illumina approach refer to Supplementary Table S3–S5 for details about the used parameters). Finally, the fractions of reference sequences covered by reads or contigs for different subsamples were calculated and visualized as line or bar plots.

Longer reads or contigs mapping to small, circular genomes, could influence mapping efficiency (Visser et al., 2016). All mappings were performed to an artificially constructed viroid sequence, which was made by joining 10 repeated viroid genome sequences. Parameters for mapping contigs from Illumina rRNA-depleted totRNA dataset to viroids were adjusted (for samples I and IV). The length fraction parameter [CLC Genomic Workbench (v12, v21)] was set to 0.5 (50 %) instead of the 0.90 (90 %) used for viruses. Every mapping was individually inspected, and in cases where contigs were longer than the reference sequence, the fraction of reference covered by contig was reported as 100%.

In the case of sample I, in some subsamples, uneven coverage of the two different PepMV strain genome sequences was observed; thus, we performed additional analyses to test whether those observations are the consequence of the presence of two strains of the same virus present in the dataset. For each sequencing approach, one subsample (in which unequal contig coverage was observed for the two viral strains) was chosen (Figure 4C, see ^*). For each of those chosen subsamples, further analyses were implemented: (i) mapping original subsampled reads to the strain better covered by contigs in the original analysis, (ii) de novo assembly of the unmapped reads, and (iii) mapping newly assembled contigs to the corresponding reference genome sequence of the other present PepMV strain.

To obtain the proxy of the speed of virus or viroid detection and compare different employed MinION sequencing approaches, the sequencing time needed to obtain 50% of viral or viroid genome covered by generated reads was calculated. For this, a script (get_cumulative_yield_table.py, Supplementary Data 3) was used to calculate the cumulative yields of reads in gigabases for every 10 min of sequencing from MinION basecalling output files (summary.txt, refer to Supplementary Table S8). According to the rarefaction analysis from section Comparison of the Performance of Nanopore and Illumina Sequencing For Detection of Viruses on Rarefied Datasets (the number of nucleotides estimated to cover more than 50% of a reference viral genome) for each virus or sequencing type combination, the time point at which this was reached was calculated according to the cumulative yields of the reads during sequencing.

Using the pipeline described above for detection of viral sequences, all viruses except one were detected in the samples employing the MinION direct RNA sequencing approach (Table 1). Using this approach, ToMV, which was present in sample I in an extremely low titer (Figure 2, Figure 3A), was not detected. TASVd (sample IV) was at first not detected by mapping reads to the viral RefSeq database; however, when using a more closely related sequence from NCBI GenBank (KY810784) as the reference, a few reads of this viroid were detected.

When using the three additional MinION sequencing approaches: direct RNA sequencing of rRNA-depleted totRNA, cDNA-PCR sequencing of totRNA, and cDNA-PCR sequencing of rRNA-depleted totRNA, all of the viruses or viroids present in the samples were detected (Figure 2).

The analysis of complete datasets showed that Illumina sequencing of rRNA-depleted totRNA resulted in a markedly higher relative fractions of viral reads compared to MinION direct RNA sequencing of totRNA (Figure 3A), which was expected due to the inclusion of the ribosomal RNA depletion step in the Illumina approach. These differences were reflected also when performing additional analyses on rarefied datasets (Figure 3); however, they were dependent on the amount of virus reads present in the original dataset.

As noted in section Nanopore Sequencing Using MinION Device Resulted in Comparable Detection of Plant Viruses/Viroids as Illumina Sequencing, no reads of ToMV were detected by MinION direct RNA sequencing of totRNA even in the complete dataset. Comparisons of rarefied subsamples further showed that MinION direct RNA sequencing performed comparably well to Illumina sequencing in cases in which a high fraction of specific viral reads was present in the samples (Figure 3). For MinION direct RNA datasets, in which virus sequences were present at more than 0.5% (PepMV in sample I and CaMV in sample II), relatively high fractions of corresponding viral genomes were covered by reads and contigs even in the smallest subsamples (Figure 3). Rarefaction analysis also showed very similar performance between the two approaches for PVeV2 (sample V)–in both cases, the fraction of genome covered by reads or contigs dropped correspondingly with the decreased dataset sizes. For TSWV, with 0.2% of the reads mapped to the viral genome, sharp drops in fractions of reference covered by reads or contigs were observed at smaller subsample data sizes for MinION direct RNA sequencing approach. For TYLCV, none of the two approaches enabled reconstruction of the complete genome, and Illumina sequencing approach performed only slightly better considering the two investigated parameters. In this analysis, for the remaining viruses and viroids (ToMV, ToCV, CCyV1, STV, CLVd, and TASVd), the MinION direct RNA sequencing method resulted in lower fraction of the genome covered by reads and contigs when compared to Illumina sequencing. In several cases, de novo assembly did not produce any contigs for the corresponding viruses, even for the complete datasets (data not shown).

Direct RNA sequencing using the MinION enabled detection of most of the plant viruses or viroids that were previously confirmed in the same samples with rRNA-depleted totRNA Illumina sequencing. However, it showed somewhat inferior performance in the systematic comparisons described in a previous section. Thus, a single sample, containing five viruses and one viroid (sample I), was selected and used to further explore the performance of three additional sample preparation or nanopore sequencing library preparation approaches, which included either rRNA depletion, sequencing of PCR-amplified cDNA, or both (Figure 1), to investigate how much of the performance deficit is due to the sample preparation or library preparation method and how much is due to the platform.

Of these approaches, the one which is most comparable to the method used for Illumina sequencing (cDNA-PCR sequencing of rRNA-depleted totRNA) resulted in very similar fractions of viral reads in the sequenced datasets (Figure 4). This approach resulted in lower fractions of specific virus or viroid reads for almost all viruses than observed for both nanopore cDNA-based sequencing approaches or for the Illumina sequencing approach. Compared to MinION direct RNA sequencing without ribosomal RNA depletion, it resulted in an increased fraction of specific viral reads for four out of 7 viruses or viroids, including the detection of one virus (ToMV), which was not detected using direct RNA sequencing alone.

MinION sequencing of cDNA-PCR without ribosomal RNA depletion also resulted in relatively high fractions of specific viral reads. For four out of seven viruses or viroids, the numbers of corresponding reads were even higher than in the rRNA-depleted dataset sequenced by the same method (Figure 4).

Moreover, even though for most viruses, the fractions of viral reads did not increase, when including rRNA depletion to direct RNA sequencing, the rarefaction analysis showed improved performance also in this case, for all but one virus, TYLCV (Figure 4).

For viroids, the Illumina rRNA-depleted totRNA sequencing approach resulted in higher fractions of viroid reads than any of the MinION sequencing approaches (Figures 3A, 4A).

De novo assembly of reads from sample I, which contained two strains of PepMV in some subsamples, resulted in an assembly of contigs corresponding only to one strain (Figures 3C, 4C). The effect was observed, when using either Illumina or nanopore sequencing approaches. After the removal of the reads of one or the other PepMV strains and performing de novo assembly again (as described in Section Comparison of the Performance of Nanopore and Illumina Sequencing for Detection of Viruses on Rarefied Datasets), the artifact was no longer observed. Additionally, for some subsamples of MinION data, the assembled contigs were longer than the reference genome sequence (Figure 4C). Further investigation (visual inspection of nanopore reads mapping to the corresponding contigs) revealed mistakes or artifacts in some mapped reads, which likely led to artifactual de novo assembly of the corresponding contigs.

Closer investigation and comparison of average read identities (proxy of sequencing error rate) for different MinION sequencing approaches for sample I revealed the lowest read identities when using direct RNA sequencing of rRNA-depleted totRNA (minimum 84.67%, average 87.43%). The highest average read identities were observed when using cDNA-PCR sequencing of the rRNA-depleted totRNA approach (minimum 90.36%, average 91.5%). For the same approach, also the highest average contig identities were observed (minimum 95.08%, average 99.73%). Even though the inclusion of the ribodepletion step resulted in a slightly decreased mean quality score (Supplementary Table S2), as was already observed in other studies (Wongsurawat et al., 2019), the calculation of average MinION read identity (Figure 2) did not show any marked differences.

As expected, the pairwise identities of consensus viral sequences compared to generated or selected references were higher in the case of Illumina sequencing approach for all viruses or viroids for which a calculation was possible (>99.5% in all cases). For MinION sequencing approaches, pairwise identities of consensus viral sequences compared to generated or selected references were higher than 98%, when average sequencing depth values were 5x or more, except in the case of ToMV (sample I), where, despite the sequencing depth of 13.27, this was 93.27%. Upon visual inspection of the mapping files, we observed that this was a consequence of uneven coverage by reads (the 5‘ of the viral genome was covered by very few reads), which contributed to the final lower average pairwise identity of consensus viral genome sequence.

A relatively short time (>30 min) was needed for retrieving sufficient data (covering at least 50% of viral genomic sequences), in cases, in which fractions of specific viral reads in samples were higher than 0.2% (PepMV, CaMV, TSWV; Figures 3, 4; Table 2) using any of the sample preparation approaches. Further, in the case of PVeV2, for which 0.0043% of reads in the sample were mapped to the virus reference sequence, ~2.5 h was needed to retrieve sufficient amount of data. For cases, in which virus or viroid reads were present in samples in very small fractions, a relatively long time was needed to retrieve sufficient amount of data; e.g., for CLVd (0.0005% of reads were mapped to viroid reference sequence), the sequencing should last at least around 12 h to retrieve enough data. Equivalent observations were made for several other viruses or viroid, sequenced within the sample I (Table 2). In some cases, especially when very low fractions of specific viral reads were observed in the samples, retrieving enough data to cover >50% of genome sequence was not achieved even if the sequencing lasted for 46–48 h (Table 2).