Full-Length Genome Sequencing of SARS-CoV-2 Directly from Clinical and Environmental Samples Based on the Multiplex Polymerase Chain Reaction Method

The new coronavirus, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is spreading worldwide with the number of confirmed cases increasing dramatically^{[1, 2]}. Coronaviruses have the largest genome among RNA viruses and have caused two major outbreaks^{[3, 4]}. The new coronavirus (SARS-CoV-2) is the third coronavirus spread globally in the past 17 years. According to the World Health Organization, SARS-CoV-2 has affected more than 155.6 millon patients in 222 countries as of 6 May 2021 and has become a major global health concern (https://www.who.int/emergencies/diseases/novel-coronavirus-2019).

Rapid genome sequencing of the virus directly from clinical and environmental specimens is important for real-time genomic surveillance in managing virus outbreaks^[5]. Genomic sequencing directly from clinical and environmental samples is challenging owing to low sensitivity. Therefore, we designed a multiplex PCR method for virus genome sequencing of SARS-CoV-2 directly from clinical and environmental samples with reference to the whole-genome sequencing method for Zika virus^[6]. SARS-CoV-2 can be detected in the bronchoalveolar lavage fluid, sputum, nasal swab, throat swab, feces, blood, and saliva^{[7, 8]}. Hence, in this study, we validated the method with different types of these samples.

The complete protocol is schematically represented in Figure 1. The first step of the protocol was to design overlapping primer pairs to cover the entire genome (Supplementary Table S1 available in www.besjournal.com). The multiplex PCR comprised a set of 102 oligonucleotide primer pairs and the amplicons generated by the primer pairs spanned the target genome (Figure 1A). The subsequent step was the amplification and sequencing of fragments that were performed with Nanopore GridION X5 to obtain reads of 400 bp or Illumina paired-end library protocol, allowing reads of approximately 300 bp. During genome sequencing, some amplicons were over-sequenced, while others were under-sequenced (Figure 2). This occurs owing to differences in the amplification process and genome content.

Figure 1.  Sequencing schemes employed in the study. (A) Schematic showing expected amplicon products for each pool of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) genome. We designed and optimized polymerase chain reaction (PCR) primers to generate amplicons that would span the SARS-CoV-2 genome. We designed 102 primer pairs, which covered 100% of the SARS-CoV-2 genome. The predicted forward primers (blue arrows) and reverse primers (red arrows) scaled according to the SARS-CoV-2 virus coordinates. (B) Workflow of the multiplex PCR method for sequencing SARS-CoV-2 directly from clinical specimens.

Figure 2.  Study of the six diluted strains using the multiplex polymerase chain reaction sequencing method on a Nanopore GridION X5. (A) Analysis of the severe acute respiratory syndrome coronavirus 2 genome assembly at various read depths. The longest contiguous sequence produced at each read depth as a fraction of the full genome length of six diluted samples is shown. (B) Summary of statistical analysis for the sequencing results of six strains.

We evaluated the sensitivity of the multiplex PCR sequencing method using a 10-fold gradient dilution of the reference strain (EPI_ISL_402119). The sequencing results for Nanopore GridION X5 showed that the mapping reads of six gradient dilution samples against the SARS-CoV-2 reference (EPI_ISL_402119) were 3,733.734 (99.71%), 3,639.845 (99.73%), 3,819.822 (99.65%), 2,934.107 (99.27%), 682.764 (91.12%), and 55.617 (44.7%), respectively (Figure 2B). We identified that as the cycle threshold (Ct) value increased, the percentage of SARS-CoV-2 reads and the percentage of SARS-CoV-2 genome coverage showed a downward trend (Figure 2A). Furthermore, when the Ct value was < 37, the coverage of the SARS-CoV-2 genome was > 95%. Interestingly, when the Ct value was > 40, the coverage of the SARS-CoV-2 genome was < 50%. The median depth and interquartile range (IQR) of coverage of all samples are listed in Figure 2B, whereas the results on the Illumina MiSeq platform are not presented.

In this study, we used 14 different types of clinical and environmental specimens, including bronchoalveolar lavage fluid, sputum, nasal swabs, throat swabs, and feces samples, from patients infected with the new coronavirus and environmental samples. An in-depth summary, including detailed information on sequence reads, depth distributions, and genome coverage per sample, of the outputs from all samples is presented in Table 1. In comparison with previous studies and despite the long genome of the virus, sequencing of the SARS-CoV-2 genome was very specific. The number of reads that mapped against the SARS-CoV-2 sequence was very high (> 96% for 10 of the samples), revealing that the primer design and the biological protocol led to some noise in the sequencing data. A total of 11 full-length or near-full-length SARS-CoV-2 genomes (> 99% genome coverage) were obtained from 14 libraries.

Sequencing depth was not uniform among the amplicons along the genome; some amplicons were over-sequenced (median depth: 8,832×), while others were only sequenced a few times or not at all. Quantitative results of samples showed that all samples with relatively low Ct values (high viral load) led to a complete or nearly complete assembly. Three of the samples, hCoV-19_Sample3, hCoV-19_Sample10, and hCoV-19_Sample14, with higher Ct values (lower viral load), led to a final assembly covering only 77.15%, 96.57%, and 98.99% of the genome sequence, respectively. The significant benefit of utilizing a multiplex PCR method for virus genome sequencing over a Sanger sequencing or an unbiased metagenomic approach was the substantial increase in the number of reads specific to the viral genome^{[6, 9]}. The number of merged reads was between 150754 and 2628756, providing sufficient sequencing data for reliable analysis of the 14 samples. We constructed a phylogenetic tree based on the full-length genome sequences of SARS-CoV-2 derived from sequencing (Supplementary Figure S1 available in www.besjournal.com).

Figure S1.  Phylogenetic tree based on full-length genome sequences of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) derived by sequencing

A limitation of this study is that our method is not suitable for identifying novel viruses as primers are SARS-CoV-2 specific. Amplicon sequencing may result in incomplete genome coverage, especially when lower abundance viral genomes are present, and the loss of both 5’ and 3’ regions. To obtain complete genomes, it may be necessary to replace the problematic primers or adjust their concentration according to other primers.

In summary, our method showed advantageous prospects for generating new coronavirus sequences directly from clinical and environmental samples; however, further studies are required to confirm this finding. On a long-term basis, this method can potentially be used as a routine laboratory test to aid in treatment, vaccine design and deployment, infection control strategies, and surveillance.