Nanopore Sequencing of SARS-CoV-2: Comparison of Short and Long PCR-tiling Amplicon Protocols

Comparison of PCR-tiling amplification protocols

The PCR-tiling amplification combined with nanopore sequencing was employed for genome sequence analysis of 152 SARS-CoV-2 isolates from Slovakia (S1 Table). The genome sequences were obtained using primer schemes generating either ∼400-bp (Artic Network version V3, https://github.com/artic-network/artic-ncov2019), ∼2-kb (35), or ∼2.5-kb long amplicons (27). In some experiments, we loaded the same sequencing libraries to both the standard and Flongle flow cells. This allowed us to compare and evaluate the MinION runs and to analyze the problems affecting the performance of nanopore sequencing.

To compare primer sets for short and long amplicons and sequencing devices, three different batches (UKBA-2, UKBA-3 and UKBA-4 in Table 1) consisting of 10-12 multiplexed samples were sequenced using multiple strategies. Fig 1 shows the comparison of the fraction of samples in a sequencing run successfully sequenced at various cut-offs measuring the total amount of sequencing normalized by the number of samples in the run. In both batches UKBA-2 and UKBA-3, 2-kb amplicons clearly outperform 400-bp amplicons. Sequencing of a mixture of longer and shorter amplicon pools provided comparable results to sequencing longer amplicons alone, perhaps because the mixture was enriched in the long amplicons. Finally, the Flongle and standard flow cells are similarly successful at comparable sequencing volumes. However, there are two disadvantages to using the Flongle flow cells. First, the Flongle cannot be washed and reused, its entire capacity is used for a single experiment. Second, since there is a large variance in the amount of data produced by a single Flongle flow cell (in our experiments, the number of active pores in Flongles ranged between 18 and 67 pores and produced between 110 and 830 Mbp - see Table 1), the capacity may be insufficient to completely recover sequences of 10 or more multiplexed samples. We consider as an important advantage that the runs using the standard flow cells can be terminated when sufficient data is collected, and thus these flow cells can be reused in further experiments after washing with nuclease containing buffer (i.e., EXP-WSH003 or EXP-WSH004). Moreover, the standard flow cells allow simultaneous sequencing of a greater number of barcoded samples with a longer run.

• Download figure
• Open in new tab

Fig 1. The fraction of successfully sequenced multiplexed samples over time.

A sample is considered as successfully sequenced if the resulting sequence produced by the Artic pipeline has fewer than 500-bp (a) or 3-kb (b) marked as missing bases. Each run is represented by several time points, each point showing the percentage of successfully sequenced barcodes (y-axis) upon reaching a specified amount of sequenced data per barcode (x-axis).

Note that batch UKBA-2 included samples with low product concentrations after PCR amplification. As a result, three samples (barcodes 02, 06 and 11) could not be completed reliably even after combining data from all six sequencing runs. Batches UKBA-3 and UKBA-4 contained only samples with Cq values from RT-qPCR below 26. S1 Figure shows the amount of missing sequence in individual samples plotted against possible explanatory variables, namely the Cq values, amplicon concentration, and RNA sample storage time prior to amplification. Although the expected trends are in some cases observable, they are not followed universally.

Using the Artic pipeline for further analysis, sequencing reads must first pass a series of filters to ensure no barcode bleeding and to remove possible contamination. The number of reads passing these filters and used for the identification of variants in the final step of the pipeline varied between runs. In our experiments their fraction comprise between 14 and 55% (Fig 2A). Majority of failed reads are due to the low quality or incompleteness (groups A-C) and comprise between 41-78% of all reads. While there are no clear differences between short and long amplicon protocols, with 2-kb amplicons these low-quality reads are apparently more prevalent on the Flongle runs compared to the standard flow cells.

• Download figure
• Open in new tab

Fig 2. Reasons for discarding reads in the Artic pipeline.

The sequencing reads must pass through a series of filters to ensure correct sample assignment and the read quality. A: reads without barcode identification. B: reads with only one barcode (Artic pipeline requires barcodes on both ends to ensure that the whole read was sequenced and to decrease the probability of barcode bleeding). C: low-quality reads (base caller quality less than 7). D: reads do not align to the SARS-CoV-2 reference. E: reads are too short (likely due to fragmentation) or F: too long (i.e. chimeric reads); lengths between 1500 and 3000 for 2-kb amplicons, between 350 and 619 for 400-bp amplicons. The reads passing all filters are included in group G. Panel (A): Summary per run. Panel (B): Detailed per-barcode analysis for UKBA-2 samples, 2-kb amplicons, standard flow cell.

Interestingly, in some runs, up to 6% of reads that pass the base quality filters do not map to the target reference genome. In particular, four samples in batch UKBA-2 of 2-kb amplicon run (barcodes 02, 07, 08 and 11) have a very high fraction of non-target reads (Fig 2B). The majority (82-96%) of these reads map to the human genome, and a smaller fraction (0.3-9%) map to bacterial genomes, including the species colonizing human oral cavity and respiratory tract (e.g., Actinomyces graevenitzii, Haemophilus parainfluenzae, Leptotrichia spp., Prevotella spp., Pseudomonas aeruginosa, Rothia mucilaginosa, Streptococcus pneumoniae, S. mitis, S. parasanguinis, S. salivarius, Tannerella forsythia, Veillonella parvula). All four samples showed a lower viral load (i.e., Cq value > 30) in RT-qPCR assays, and the amplification in the PCR-tiling protocol resulted in lower product yield. Human and bacterial reads represent artefacts apparently resulting from a non-specific amplification of contaminating nucleic acids present in clinical samples.

We have also observed that some amplicons originate from sub-genomic RNAs that co-purify with the SARS-CoV-2 genomic RNA. It has been demonstrated that the amount of sub-genomic RNAs correlates with the disease severity. As these molecules are strongly repressed in asymptomatic patients (38), their proportion in the sequencing data can serve as a molecular marker. The most abundant reads are derived from the N mRNA (39). The sub-genomic RNAs are generated in the process of the virus replication/transcription (5) and start with a leader sequence originating from the untranslated 5’ end of the viral genome, followed by a downstream sequence containing a particular open reading frame. The leftmost primer in both 400-bp and 2-kb primer sets investigated in this study is contained within the leader sequence. This facilitates amplification of sub-genomic RNAs with appropriate right primers (Fig 3).

• Download figure
• Open in new tab

Fig 3. Sub-genomic RNAs

(black) and amplicons of primer pool 1 from the 2-kb primer set (blue), visualized by the UCSC genome browser (40).

Table 2 lists the fraction of selected sub-genomic RNAs among reads that could be aligned to the SARS-CoV-2 genome. These fractions are relatively low, with the remaining sub-genomic RNAs being even more rare. However, the fractions vary among the samples. In UKBA-2 run with 2-kb amplicons, the highest fraction of 14.3% was observed for the gene N mRNA in barcode 07 and the fraction of 7.5% was observed for the ORF3a mRNA in barcode 11. Some of these sub-genomic amplicons are discarded from the analysis as too short, while others lead to lower coverage in areas not covered by the sub-genomic RNA (Fig 4).

• Download figure
• Open in new tab

Fig 4. Coverage along the genome in two MinION runs for batch UKBA-2.

In both runs, an initial portion of the run containing on average 40-Mbp of sequencing data per barcode was used. Coverage values higher than 1000 were clipped at this value and are shown in blue. Coverage below 20 (default Artic cutoff) is shown in red. Medians of 10-bp windows are shown for smoothing. The very starts and ends of the genome are not covered by amplicons, and are thus displayed in red. Shaded area in the left column corresponds to amplicon 13. Some barcodes have a visible dip in the coverage at the left end of this amplicon; this difference in coverage is caused by reads originating from sub-genomic RNAs corresponding to the gene S. Similar plots for additional runs are shown in S2 Fig.

From these pilot experiments, we conclude that even though 400-bp amplicons have a lower percentage of discarded reads (Fig 2), they produce fewer finished sequences at a comparable overall amount of sequence data (Fig 1). The reason is a very uneven coverage of individual amplicons (Fig 4). This is observed in both sets of primers, but for the 400-bp amplicons we see a much lower coverage in the worst covered regions (Fig 5). Additional sequencing runs (UKBA-6, UKBA-10, UKBA-11, and UKBA-12) were performed with long 2-kb amplicons on standard MinION flow cells with similar results (Fig 2A; S2 Fig).

• Download figure
• Open in new tab

Fig 5. Coverage distribution in different sequencing runs.

For each barcode, coverage by reads passing the Artic filter was computed along the genome and the distribution of the coverage values was summarized as a violin plot (blue), cropped at coverage 1000. Orange dots represent median coverage and green dots 10% percentile (corresponding to a approx. 3-kb portion of the genome with the lowest coverage). In all runs, an initial portion containing on average 40-Mbp of sequencing data per barcode was used.

To investigate if a different primer scheme for generating long amplicons can solve the problem with uneven coverage (in particular, see amplicon 13 in Fig 4 which partially covers the S gene region important for identification of the SARS-CoV-2 Variants of Concerns), we also tested the 2.5-kb primer panel (27). Except for the leftmost primer, the primer positions in this panel differ from those of the 2-kb scheme. We have performed two sequencing runs with 2.5-kb primer set (UKBA-19, UKBA-21). In the first experiment, we have noticed an almost complete drop of coverage in the last amplicon derived from the 3’ end of the genome; for the second experiment, we have replaced the primers for the right-most amplicon with the right-most primers from the 2-kb panel, which mitigated the issue. Comparing the coverage of individual amplicons between the 2-kb and 2.5-kb schemes (Fig 6), the coverage in the 2.5-kb scheme indeed appears to be more even. Fig 5 illustrates that our modification of the 2.5-kb scheme leads to a particularly small difference between the median coverage and coverage of the lowest 10% of the genome, which may result in fewer regions with insufficient coverage. However, we have also noticed a higher percentage of failed reads, with only 24% (UKBA-19) and 16% (UKBA-21) reads passing all filters and being usable for variant identification. Further analysis revealed a notable increase in single-barcode reads (group B) and shorter than expected reads (group E), pointing to difficulties in amplifying and sequencing longer fragments. More experiments are required to determine whether the 2.5-kb scheme results in more fully-assembled genomes over the 2-kb scheme.

• Download figure
• Open in new tab

Fig 6. Average coverage along the genome

for seven runs with 2-kb amplicons (batches UKBA-2,3,4,6,10,11,12) and two runs with 2.5-kb amplicons (UKBA-19 with the original primer set and UKBA-21 with the last primer pair replaced by its counterpart from the 2-kb scheme). Each line depicts the average coverage over all samples in a run at the time point when 40-Mbp per sample was sequenced on average. Medians of 50-bp windows are shown for smoothing. Note a drop-out in the amplicon 13 (2-kb scheme) which covers a 3’ end of orf1b and about a third of the S gene including the region associated with mutations in Variants of Concern such as B.1.1.7.

Collection of samples and RNA preparation

Oropharyngeal swabs of patients with suspected COVID-19, collected between March 30, 2020 and March 19, 2021, were preserved in 2-3 ml of viral transport medium and delivered to the laboratory of the Biomedical Research Centre of the Slovak Academy of Sciences in Bratislava, Slovakia in frame of the routine RT-qPCR diagnostics for SARS-CoV-2. Initially (UKBA-2 samples), 100 µl of the swab medium was used for the RNA extraction using the Zymo Research Quick-RNA™ Viral 96 Kit (Zymo Research, Irvin, California, USA). Resulting RNA was eluted to 20 µl of nuclease free water. For all other specimens, the Biomek i5 Automated Workstation (Beckman Coulter, Indianapolis, Indiana, USA) was employed using the RNAdvanced Viral kit (Beckman Coulter, Indianapolis, Indiana, USA). In this case, RNA was extracted from 200 µl of swab medium and eluted to 40 µl of nuclease free water.

Real-time quantitative PCR (RT-qPCR)

In frame of the routine RT-qPCR diagnostics, presence of SARS-CoV-2 RNA was detected by vDETECT COVID-19 RT-qPCR kit, rTEST COVID-19 RT-qPCR kit or rTEST COVID-19 RT-qPCR ALLPLEX kit (MultiplexDX, Bratislava Slovakia) targeting RNA-dependent RNA polymerase (RdRp) and Envelope (E) genes. The RT-qPCR assays were carried on QuantStudio(tm) 5 Real-Time PCR System (Applied Biosystem, Foster City, California, USA).

Library preparation and DNA sequencing

The sequencing libraries were constructed using a ligation kit (SQK-LSK109) essentially as described in a PCR-tiling of COVID-19 virus protocol (PTC_9096_v109_revF_06Feb2020; Oxford Nanopore Technologies) with minor modifications. Briefly, RNA samples extracted from swabs positive for the presence of SARS-CoV-2 in RT-qPCR assay (quantification cycle (Cq) values 13.46 - 32.03; S1 Table) were converted into cDNA using a SuperScript IV reverse transcriptase (Thermo Fisher Scientific) or LunaScript® RT SuperMix Kit (New England Biolabs). For each sample, the overlapping amplicons were generated using a Q5® Hot Start High-Fidelity DNA polymerase (New England Biolabs) and the primer pools spanning the SARS-CoV-2 genome sequence (i.e., 400-bp Artic nCoV-2019 V3 panel (https://github.com/artic-network/artic-ncov2019) purchased from Integrated DNA Technologies (IDT, cat.no. 10006788) and the 2-kb (35) and 2.5-kb schemes (27), custom synthesized by Microsynth). The same cycling program was used for all amplicon types (i.e., 30 sec initial denaturation at 98°C, followed by 25 to 35 cycles of 15 sec at 98°C (denaturation) and 5 min at 65°C (combined annealing and polymerization), and cooling to 4°C). The amplifications were performed in two separate reactions and the overlapping amplicons were pooled, purified using an equal volume of AMPure XP magnetic beads (Beckman Coulter) and quantified using a Qubit 3.0 spectrophotometer and dsDNA Broad Range Assay Kit (Thermo Fisher Scientific). About 50-75 ng (400-bp amplicons) and 250-300 ng (2 and 2.5-kb amplicons) of each SARS-CoV-2 isolate were treated with NEBNext Ultra II End repair / dA-tailing Module (New England Biolabs). The samples were then barcoded using EXP-NBD104 (barcodes 1-12) or EXP-NBD114 (barcodes 13-24) kits (Oxford Nanopore Technologies) and NEBNext Ultra II Ligation Master Mix (New England Biolabs). Barcoded samples were pooled and purified using 0.6 volume of AMPure XP magnetic beads. The AMII sequencing adapter (Oxford Nanopore Technologies) was ligated to about 75 ng (400-bp amplicons) or 300 ng (2 and 2.5-kb amplicons) of barcoded pools using Quick T4 DNA ligase (New England Biolabs) and the sequencing libraries were purified using 0.6 volume of AMPure XP magnetic beads. About 20 ng (400-bp amplicons) and 90 ng (2 and 2.5-kb amplicons) of the libraries were loaded on an R9.4.1 flow cell (FLO-MIN106). The sequencing was performed using a MinION Mk-1b device (Oxford Nanopore Technologies). For sequencing on the Flongle flow cells (FLO-FLG001), the library preparation was the same, except that one third to one half of the library was loaded compared to the amount used for the standard flow cell.

Data processing

Nanopore sequencing data were base called and demultiplexed using Guppy v.3.4.4. Variant analysis was performed using Artic analysis pipeline v.1.1.3. (https://artic.network/ncov-2019/ncov2019-bioinformatics-sop.html) using recommended settings. Minimum and maximum read lengths in the Artic guppyplex filter were set to 350 and 619 for the 400-bp amplicons and 1500 and 3000 for both the 2 and 2.5-kb amplicons, respectively. For batch UKBA-2, the final sequences were produced by first combining sequencing reads from both standard and Flongle runs with the same primer set and running the Artic pipeline. Subsequently the results for the two primer sets were combined so that regions sufficiently covered by at least one amplicon set were considered as finished. The same process was used in batch UKBA-3, but there was only data from standard flow cells available. Subsequent batches were based on 2 or 2.5-kb amplicons sequenced on a standard flow cell.

To compare different primer sets and sequencing devices, reads were also demultiplexed at the less strict default Guppy settings and aligned to various reference genomes by minimap2 v. 2.13-r852-dirty (41). Reference genomes include the SARS-CoV-2 genome MN908947.3 (1), the human genome version hg19 downloaded from the UCSC genome browser (40), and the database for bacterial species typing included in the Japsa software (42). To detect sub-genomic RNAs, reads were aligned to transcripts downloaded from the UCSC genome browser by minimap2, and classified as sub-genomic, if the alignment to a sub-genomic RNA has at least 5 matches more than the best alignment to the reference genome. Read coverage was computed using genomecov tool from BEDTools (43) with options -bga -split.

To compare the results for various sequencing data volumes, reads were ordered by the sequencing finish time and the initial portion with the desired total length was selected and used for the analysis in the Artic pipeline. To compare batches with a different number of samples, the cutoffs were expressed as the average amount per barcode.

Data availability

The SARS-CoV-2 genome sequences determined in this study were deposited in the GISAID database (http://www.gisaid.org). The accession IDs are shown in S1 Table. The sequencing reads mapping to the reference genome were deposited to ENA under study PRJEB44303.