A SARS-CoV-2 Delta variant containing mutation in the probe binding region used for RT-qPCR test in Japan exhibited atypical PCR amplification and might induce false negative result.

INTRODUCTION: A recent pandemic of SARS-CoV-2 infection has caused severe health problems and substantially restricted social and economic activities. RT-qPCR plays a vital role in the diagnosis of SARS-CoV-2 infection. The N protein-coding region is widely analyzed in RT-qPCR to diagnose SARS-CoV-2 infection in Japan. We recently encountered two cases of SARS-CoV-2-positive specimens showing atypical amplification curves in the RT-qPCR.
METHODS: We performed whole-genome sequencing of 63 samples (2 showing aberrant RT-qPCR curve and 61 samples infected with SARS-CoV-2 simultaneously in the same area) followed by Phylogenetic tree analysis.
RESULTS: We found that the viruses showing abnormal RT-qPCR curves were Delta-type variants of SARS-CoV-2 with a single nucleotide mutation in the probe-binding site. There were no other cases with the same mutation, indicating that the variant had not spread in the area. After searching the database, hundreds of variants were reported globally, and one in Japan contained the same mutation. Phylogenetic analysis showed that the variant was very close to other Delta variants endemic in Japan but quite far from the variants containing the same mutation reported from outside Japan, suggesting sporadic generation of mutant in some domestic areas.
CONCLUSIONS: These findings propose two key points: i) mutations in the region used for SARS-CoV-2 RT-qPCR can cause abnormal amplification curves, and ii) various mutations can be generated sporadically and unpredictably; therefore, efficient and robust screening systems are needed to promptly monitor the emergence of de novo variants.

Introduction

Since its first report from a seafood market in Wuhan, China, in 2019 [1,2], coronavirus disease 2019 (COVID-19) has spread over six continents within a short period; the World Health Organization (WHO) declared COVID-19 a global pandemic in March 2020 [3]. As of November 4, 2021, there have been more than 247.4 million confirmed cases of COVID-19 worldwide, with over 5 million deaths attributed to it [4]. A metagenomic RNA sequencing study by Wu et al., in 2020 revealed that the virus causing COVID-19 is a novel strain of coronavirus, closely associated with severe acute respiratory syndrome coronavirus (SARS-CoV) and the middle east respiratory syndrome coronavirus (MERS-CoV). Therefore, it was named “severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)" [2].

The whole viral genome sequence of SARS-CoV-2 was first reported in January 2020 [2]. Since then, several quick and reliable diagnostic tests have been developed, including reverse transcription-quantitative polymerase chain reaction (RT-qPCR)-based diagnostic methods [5,6], CRISPR/Cas12 based methods [7], specific high-sensitivity enzymatic reporter unlocking (SHERLOCK) assay [8], and loop-mediated isothermal amplification (LAMP) [9]. To date, the United States Food and Drug Administration (USFDA) has issued emergency use authorization (EUA) for 266 molecular diagnostic tests that were developed based on the sequence information provided by the SARS-CoV-2 genome [10]. Although each method has advantages and disadvantages, RT-qPCR based methods have been the most common diagnostic methods applied to detect SARS-CoV-2 in patient samples. Several protein-coding sequences of the SARS-CoV-2 genome have been used for qRT-PCR, such as genes encoding RNA-dependent RNA polymerase (RdRP) [11], Spike (S), Nucleocapsid (N) [12] and Envelope (E) proteins [13].

In Japan, the multiplex Quantitative reverse-transcription PCR (RT-qPCR) based diagnostic method was first developed by the National Institute of Infectious Disease (NIID) in 2020 [14] and applied for the mass diagnosis of SARS-CoV-2 in 2020 [15]. NIID suggested the N1 and N2 regions of the N protein-coding sequences as the target for RT-qPCR amplification (Fig. 1 A). Since then, 26.6 million tests have been performed in Japan, among which 6% were positive for COVID-19 [16]. When the study was conducted, 124,358 diagnostic tests were performed in Kumamoto Prefecture in Japan, of which 14,392 tests were positive [17].

Fig. 1

(A) Location of SARS-CoV-2 genome where the RT-qPCR primer-probe sets target for detecting SARS-CoV-2 in patient samples. Both primer-probe sets (N1 and N2) target the N protein-coding sequences. N1_F, N1_R, N1_P are the primers and probe of the N1 set and N2_F, N2_R. N2_P are the primers and probe of the N2 set (Figure not drawn to scale). (B) RT-qPCR curve to detect SARS-CoV-2 using NIID recommended N2 primer-probe set produced aberrant curve in two samples (A8167 and A8168) compared to other positive samples in the same run.

RT-qPCR based testing plays a vital role in identifying infected patients and confining them where uninfected individuals can perform social activities [18]. Therefore, accuracy in detecting SARS-CoV-2 in patients plays a crucial role in controlling the transmission of this highly contagious pathogen. False-negative results put much stress on the containment efforts for COVID-19. Successful RT-qPCR based detection is dependent on the efficient binding of primers and probes to the target areas. Any changes in the target nucleotide sequence can significantly lower the binding affinity of the primers and probes to their target region, resulting in a pseudo-negative diagnosis [19,20]. Consequently, the resulting undetected mutants can also result in future waves of COVID-19 [21]. Although SARS-CoV-2 has a lower mutation frequency compared to other RNA viruses because of its RNA proofreading mechanism [22], its higher transmission rate makes SARS-CoV-2 a good candidate to observe for any unforeseen mutation within the RT-qPCR primer and probe binding regions so that the quality of the detection methods could be ensured.

This article reports a point mutation in the N gene of SARS-CoV-2 that resulted in an atypical RT-qPCR curve during the diagnostic test, leading to a dubious diagnostic interpretation. Sequencing results revealed the presence of a single point mutation at the probe-target site in the N gene. This analysis shows that sequencing can play a problem-solving role in RT-qPCR based diagnostic complications and raise caution to institutes performing RT-qPCR tests to detect SARS-CoV-2 infection.

Materials and methods

Regulatory approvals

In this study, samples were collected by the Kumamoto City Medical Association Inspection Center (hereinafter referred to as the Kumamoto PCR center). The downstream analyses were performed at the Division of Genomics and Transcriptomics in the Joint Research Center for Human Retrovirus Infection, Kumamoto University [Approved by the Ethics Review Committee of the Faculty of Life Sciences Kumamoto University (Approval no. 2223)].

Specimen collection and storage

For this study, sputum and/or nasopharyngeal swab samples were collected from the Kumamoto city area from August 16, 2021, to September 8, 2021, as part of the regular COVID-19 diagnostic service. The collected samples and extracted RNAs were preserved at −80 °C until further experiments were performed according to the guidelines of the NIID, Japan [23].

RNA extraction and qRT-PCR

For RT-qPCR analysis and sequencing, RNA was extracted from the sputum and/or nasopharyngeal swab samples (140 μL) using the QIAamp Viral RNA Mini Kit (QIAgen) and EZ1 Virus Mini Kit v2.0 (QIAgen) according to the manufacturer's instructions. The final RNA product was eluted in 60 μL of buffer AVE (RNase free water with 0.04% sodium azide), and 5 μL of the eluted RNA was used for qRT-PCR, and the remaining RNA was stored at −80 °C for further analysis. RT-qPCR was performed on a LightCycler96 System (Roche, Basel, Switzerland) according to the protocol suggested by Shirato et al. [14] to detect the N gene using the One Step PrimeScript™ RT-PCR Kit (Perfect Real Time) (Takara). The N2 primer-probe set was used as described by Shirato et al. [14].

SARS-CoV-2 whole-genome sequencing (WGS)

cDNA synthesis, viral sequence enrichment, library amplification, and indexing were performed using the QIAseq DIRECT SARC-CoV-2 kit (QIAgen), following the manufacturer's recommendations. QIAseq DIRECT UDI set-A (QIAgen) was used for multiplexing the samples. SARS-CoV-2 libraries of 25 μL volume for each sample were prepared at the end of the process. The quality of the enriched SARS-CoV-2 libraries was evaluated by electrophoresis with a TapeStation 4150 system (Agilent Technologies). Finally, the prepared libraries were denatured and subjected to sequencing using MiSeq reagent Micro and Nano Kits (Version 2, 300 cycles) in the MiSeq desktop sequencing system (Illumina).

Data analysis

Upon sequencing in the Illumina MiSeq sequencer, one FASTQ file (Read1) was generated for each sample. Adapter sequences were trimmed from the Read1 using the ‘cutadapt’ tool [24]. After the adapter trimming step, a cleaning step was performed using the PRINSEQ tool [25] and Read1 with Phred score >20 were used for downstream analysis. Adapter trimmed and cleaned reads were then aligned to the SARS-CoV-2 reference genome NC_045512 (isolate Wuhan-Hu-1) using the BWA-MEM algorithm [26]. Subsequently, the Samtools program [27] was used to remove multiple aligned reads, and Freebayes (Version 1.2.0) command-line tools [28] were used to call the variants from the aligned reads and create a variant call format (VCF) file. Finally, the aligned files and VCF files were visualized using the Integrative Genomics Browser (IGV) [29]. The resulting consensus SARS-CoV-2 genomes were deposited at the Global Initiative on Sharing All Influenza Data (GISAID) [30] (Supplementary data 1).

For secondary data analysis, the Pangolin COVID-19 Lineage Assigner tool was used to assign the nomenclature for the viral genomes proposed by Rambaut et al. [31] and the CoV-GLUE web application tool was used to check for any novel mutations within the viral sequences [32].

Phylogenetic tree analysis

For local phylogenetic tree construction, 61 SARS-CoV-2 strains were collected from Kumamoto city between September 1, 2021, and September 8, 2021. At first, those 61 viral sequences along with 2 sequences containing G29234A mutation were aligned to the SARS-CoV-2 reference genome NC_045512 (isolate Wuhan-Hu-1) using Geneious Prime (version 2020.2.4) (https://www.geneious.com/). Then, the aligned sequences were used to construct a maximum-likelihood phylogenetic tree using phyML 3.0 [33] with Smart Model Selection (SMS) [34].

SARS-CoV-2 genomes and associated metadata were downloaded from the GISAID database [30] (accessed on October 18, 2021) to construct a global phylogenetic tree. Global genome collections were downloaded from the “Region-specific Auspice source file” of the GISAID database, resulting in 3580 viral genomes (globally distributed random viral strains collected from December 2019 to October 2021) for global phylogenetic tree construction. Globally distributed viruses with the G29234A mutation were identified using the ‘substitution’ tool from the GISAID database (accessed on October 18, 2021) [30]. Then, the resulting 250 viral genomes containing the G29234A mutation were added to the Nextstrain build [35] as ‘focal’ sequences. The “global genome collections” were kept as ‘contextual’ genomes to build the phylogenetic tree. No non-human viral hosts were considered for tree construction.

Results

Real-time PCR showed an atypical curve

NIID initially recommended a multiplex RT-qPCR system for detecting SARS-CoV-2 using two primer and probe sets (N1 and N2) [14]. Subsequently, according to the 3rd edition of “Guidelines for the operation of the new coronavirus (SARS-CoV-2) test methods, NIID recommends using N2 primer-probe set with one enforcement to detect viral RNA in the samples during the epidemic period in which outbreaks continue in multiple municipalities [36]. Following the instructions from NIID, the Kumamoto PCR center used the N2 primer-probe set for diagnosis purposes. In this RT-qPCR experiment, 20 samples collected from the Kumamoto city area in the same time frame (Supplementary data1) were run in duplicate. From the RT-qPCR curve, it was seen that two of the samples had lower C values (Fig. 1, Table 1 ); however, the fluorescence intensity of samples A8167 and A8168 was not amplified in the same way as the other samples during the same run (Fig. 1). This unusual amplification of the RT-qPCR curve gave rise to the possibility of some variations in the primer-probe target sequences of the N gene of those samples.

Table 1

Cycle threshold (C) values of 20 samples collected from Kumamoto city area from August 16, 2021, to September 7, 2021.

	Ct values
Sample	Read 1	Read 2
A5631	33.78	33.83
A5715	33.19	33.15
A5806	35.14	34.65
A5993	32.92	32.92
A6131	30.11	30.28
A6203	28.35	28.32
A6338	33.39	33.46
A6450	29.05	29.02
A6563	32.23	32.5
A6844	30.58	30.57
A6881	27.12	27.03
A6954	30.13	30.06
A6979	28.86	28.81
A7122	27.27	27.25
A7123	31.15	31.11
A7410	26.26	26.26
A7456	32.04	31.92
A7885	30.38	30.4
A8167	20.19	20.02
A8168	24.3	24.36

SARS-CoV-2 whole-genome sequencing revealed G29234A mutation

Based on the hypothesis mentioned earlier, whole-genome sequencing was performed initially on seven samples (five samples showing characteristic RT-qPCR curve and two samples showing aberrant RT-qPCR curve from the same run) following the protocol mentioned in the Materials and Methods section. IGV visualization of the VCF files revealed the presence of a G→A mutation at 29234 loci of the SARS-CoV-2 genome of A8167 and A8168 samples, whereas no mutation was found in the same locus for the rest of the samples (Fig. 2 A). Coverage depth analysis revealed that the reads covering 29,234 loci were enriched with mutated base A (n = 2296, 99.9%), whereas the wild type base G was supported by only three reads (0.1%) for sample A8167 (Fig. 2B). A similar pattern was observed in sample A8168 as well, where the mutated A base was supported by 1,066 reads (100%) compared with the wild type G base with zero (0%) reads in the mentioned position of the SARS-CoV-2 genome (Fig. 2B). This mutation overlaps with the N gene probe used in the recommended N2 primer-probe set (Fig. 1A). It can be assumed that the G29234A mutation in samples might be a cause of the aberrance in the RT-qPCR curve by considering the sequencing data obtained from the samples and their corresponding RT-qPCR curve.

Fig. 2

(A) Integrative Genomics Viewer (IGV; https://software.broadinstitute.org/software/igv/) screenshot depicts the Variant Call Format (VCF) file showing point mutation for the seven sequenced samples. Single nucleotide mutations were observed at the 29234 loci of the SARS-CoV-2 genome for both A8167 and A8168 samples. The orange bar labelled “NIID_2019-nCOV_N2_P″ corresponds to the region covered by the N2 probe suggested by Shirato et al. [14]. (B) The BAM files show the reads and coverage for A8167 (left) and A8168 (right). Numbers in the box correspond to the raw read number and the percentage of the reads at the 29234 loci of the SARS-CoV-2 genome.

Phylogenetic analysis showed the similarities among the local strains in Kumamoto

Whole-genome sequencing of another 61 samples collected from Kumamoto city within the same time frame was performed to investigate whether the variant containing the G29234A mutation was endemic or sporadic in the Kumamoto area. The resulting sequences showed no variation in the 29234 loci of the SARS-CoV-2 genome (Supplementary Fig. 1). For further confirmation, a local phylogenetic tree was constructed using the variants containing the G29234A mutation and the variants without the mutation collected from Kumamoto city within the same time frame (Fig. 3 ). The variants containing the G29234A mutation clustered within the same clade (AY29; Delta variant Japanese sub-lineage [37]) in the local phylogenetic tree. The horizontal branch distances in the tree did not deviate much for the mutated samples when compared with the rest of the samples from the same clade. These results confirmed that the mutant samples (A8167 and A8168) are genetically similar to the local strains but differ in mutation landscape and that the G29234A mutation was sporadic at the time of reporting in the local area.

Fig. 3

Maximum likelihood phylogenetic tree constructed using Geneious prime (version 2020.2.4) (https://www.geneious.com/prime/) with 63 locally collected samples from Kumamoto city. The samples represent the same locality within the same time frame. The scale bar indicates the number of substitutions per site. The Wuhan-hu-1 (NC_045512; highlighted in blue) is placed as the root sequence of the tree, and the red-highlighted samples contain the G→A mutation at 29234 loci of the SARS-CoV-2 genome. AY is the alias used for B.1.617.2 lineage (Delta variant).

Next, SARS-CoV-2 sequences deposited in the GISAID database were analyzed to determine whether the same mutated variant was distributed across different geographical locations. When highly covered genomes (>29,000 bp), excluding the low-coverage genomes (>5% Ns), were searched for the mutations, 250 samples were found to have the same G→A mutation at 29234 loci from 3,341,006 submitted genomes to the GISAID database (mutation frequency 0.007%; accessed on October 18, 2021). After creating the global phylogenetic tree with randomly sampled sequences from different locations across the world along with the mutated viral samples, several endemic clusters were identified throughout the course of the SARS-CoV-2 pandemic (Fig. 4 ). Initially, a sporadic variant group with G29234A mutation was observed in Canada in April 2020, followed by two isolated clusters in England and Germany between December 2020 and June 2021. Nevertheless, the same mutation was observed in groups mainly concentrated in the USA and Europe in recent times (Fig. 4). Another phylogenetic tree was constructed with the globally distributed G29234A variants and the variants sampled in the Kumamoto city containing the G29234A mutation (Fig. 5 ) to understand a possible source of the mutated variant found in Kumamoto city. The resultant tree clearly showed cluster separation between the G29234A variants sampled from the rest of the world and the samples collected from the Kumamoto city (Fig. 5; zoomed in-right panel). Considering the factors such as relative cluster position of G29234A mutants collected from the Kumamoto area (Figs. 3 and 5) to the other mutant strains in the tree and the sample collection dates (Supplementary data2), it could be inferred that the mutated variants might not have migrated from outside. Instead, they acquired these mutations sporadically in Japan.

Fig. 4

Nextstrain [35] generated phylogenetic tree with viral strains containing A allele at 29234 loci of the SARS-CoV-2 genome (highlighted in orange).

Fig. 5

Maximum likelihood phylogenetic tree built in Geneious prime (version 2020.2.4) (https://www.geneious.com/prime/) showed the relationships of viral strains collected from Kumamoto area from September 1, 2021 to September 8, 2021 with the viral strains having G29234A mutation from around the world. The Wuhan-hu-1 (NC_045512; highlighted in blue) is placed as the root sequence of the tree, and the green-highlighted samples represent the viral strains collected from Kumamoto city. Red-bold faced viral strains (Zoomed in-Right panel) represent viral strains collected from the Kumamoto area, which contain G29234A mutation.

Discussion

From the position of primers and probes suggested by Shirato et al. [14], it could be said that the presence of G→A mutation at 29234 bp [position relative to reference sequence NC_045512.2] played a role in obtaining the abnormal RT-qPCR curves while diagnosing the patients using the NIID recommended N2 primer-probe set in Kumamoto PCR center. According to the “Guidelines for the operation of the new coronavirus (SARS-CoV-2) test methods (3rd Edition)", the sensitivity of the N2 primer-probe set was higher than that of the N1 set with all test samples during the protocol development. Similar results were found even in samples with low viral titers [36]. Therefore, the NIID advised all test centers to follow their recommended operation instructions during different endemic conditions. NIID recommended using the N2 primer-probe set with one enforcement during the epidemic period in which outbreaks continue in multiple municipalities [36]. Therefore, the similar use of the N2 primer-probe set might have been adopted in other prefectural test centers during the ongoing epidemic waves in Japan. Given this situation, the possibility of the presence of the G29234A mutation in the circulating strains might jeopardize some test results by returning a dubious RT-qPCR curve [19,20].

Initially, the G29234A mutation showed a homoplasy trait. It appeared spontaneously in two clades at two time points: in 20C (clustered in Canada) during the early pandemic period in 2020 and 20I (also known as Alpha strain) from December 2020 to June 2021, mainly concentrated in England and Germany (Fig. 4). However, from the recent trait, it was seen that the mutation had been spreading among the Delta variants of SARS-CoV-2. Although the frequency of this mutation was 0.007% compared to the total submitted genomes in the GISAID database at the time of reporting, the existence of this mutation in multiple clusters among the 21J (Delta) clade (Fig. 4, Supplementary Fig. 2) in different geological locations at the same time raised the possibility of local spreading. It also increased the probability that the mutations observed in the Kumamoto area might have migrated from external sources. However, the phylogenetic analysis revealed that the strains that acquired the G29234A mutation in Japan were distantly located in a separate cluster with other strains collected from the Kumamoto area compared to those with the same mutation in different regions of the world (Fig. 5), which suggested sporadic acquisition of the mutations in Japan instead of being migrated from other external sources.

While searching for a possible causal mechanism for the G292324A mutation, we found a 5′-UGG (edited G is underlined) trinucleotide motif in the genomic strand of SARS-CoV-2. This complementary strand is 5′-CCA (edited C is underlined), suggesting the sequence preference of RNA for the APOBEC3G cytosine deaminase [38,39]. However, at least three studies with SARS-CoV-2 genomic sequence analyses have shown that changes in C→U are biased in the viral genomic stand [[40], [41], [42]]. Therefore, the G292324A mutation may not be due to deamination by APOBEC3 proteins. Further investigation is required to confirm the involvement of APOBEC3 enzymes and/or other mechanisms in the G292324A mutation. In addition, the G29234A mutation caused a change in the amino sequence from glycine (G) to serine (S) at 321-position of the N protein. To determine why this mutation in the N protein did not sustain the other clades during the early pandemic and why this mutation had been spreading within the Delta variants requires further investigation.

The above findings indicate two key points. First, mutations in the primer-probe target regions could cause an atypical RT-qPCR curve, leading to false-negative results. We would like to encourage examiners to check the amplification curve if the data is available. Second, various mutations can occur sporadically and unpredictably. Therefore, efficient and robust screening systems are deemed necessary to promptly monitor the emergence of new variants of interest. Based on the experience we faced this time, we would like to suggest ways to monitor de novo mutations' emergence. Whole viral genome sequencing should be performed with 5–10% of total positive samples to monitor any unexpected mutation depending on human resources and facility availability. Alternatively, known mutations could be analyzed using specific probes that recognize each known mutation. This approach would be useful, especially in areas where WGS is unavailable or when we need to increase the examination coverage for some particular mutations. An improved diagnostic test that can evaluate multiple mutations simultaneously is urgently required to establish an efficient monitoring system for variant viruses.

Authors' contribution

All authors meet the IGMJE authorship criteria. YS and TK conceived and designed the study; YO collected the samples; SAR, YO and MBH performed the experiments; SAR and YS analyzed the data; SAR, YS, TK and TI wrote and edited the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by research grants from the (AMED) (JP21fk0410023, JP21wm0325015, JP21jm0210074) to Y. Satou, the , (JST) , Shin-Nihon Foundation of Advanced Medical Research and an intramural grant from the to T. Ikeda. The funders had no role in the study design, data collection, data interpretation, or discussion regarding submission for publication.

Data availability statement

The SARS-CoV-2 genome sequences were downloaded from GISAID, subject to GISAID's terms and conditions (https://www.gisaid.org/registration/terms-of-use/).

Declaration of competing interest

The authors declare no conflict of interest.