Impact of SARS-CoV-2 Mutations on PCR Assay Sequence Alignment

Real-time reverse transcription polymerase chain reaction (RT-PCR) assays are the most widely used molecular tests for the detection of SARS-CoV-2 and diagnosis of COVID-19 in clinical samples. PCR assays target unique genomic RNA regions to identify SARS-CoV-2 with high sensitivity and specificity. In general, assay development incorporates the whole genome sequences available at design time to be inclusive of all target species and exclusive of near neighbors. However, rapid accumulation of mutations in viral genomes during sustained growth in the population can result in signature erosion and assay failures, creating situational blind spots during a pandemic. In this study, we analyzed the signatures of 43 PCR assays distributed across the genome against over 1.6 million SARS-CoV-2 sequences. We present evidence of significant signature erosion emerging in just two assays due to mutations, while adequate sequence identity was preserved in the other 41 assays. Failure of more than one assay against a given variant sequence was rare and mostly occurred in the two assays noted to have signature erosion. Assays tended to be designed in regions with statistically higher mutations rates. in silico analyses over time can provide insights into mutation trends and alert users to the emergence of novel variants that are present in the population at low proportions before they become dominant. Such routine assessment can also potentially highlight false negatives in test samples that may be indicative of mutations having functional consequences in the form of vaccine and therapeutic failures. This study highlights the importance of whole genome sequencing and expanded real-time monitoring of diagnostic PCR assays during a pandemic.

Real-time reverse transcription polymerase chain reaction (RT-PCR) assays are the most widely used molecular tests for the detection of SARS-CoV-2 and diagnosis of COVID-19 in clinical samples. PCR assays target unique genomic RNA regions to identify SARS-CoV-2 with high sensitivity and specificity. In general, assay development incorporates the whole genome sequences available at design time to be inclusive of all target species and exclusive of near neighbors. However, rapid accumulation of mutations in viral genomes during sustained growth in the population can result in signature erosion and assay failures, creating situational blind spots during a pandemic. In this study, we analyzed the signatures of 43 PCR assays distributed across the genome against over 1.6 million SARS-CoV-2 sequences. We present evidence of significant signature erosion emerging in just two assays due to mutations, while adequate sequence identity was preserved in the other 41 assays. Failure of more than one assay against a given variant sequence was rare and mostly occurred in the two assays noted to have signature erosion. Assays tended to be designed in regions with statistically higher mutations rates. in silico analyses over time can provide insights into mutation trends and alert users to the emergence of novel variants that are present in the population at low proportions before they become dominant. Such routine assessment can also potentially highlight false negatives in test samples that may be indicative of mutations having functional consequences in the form of vaccine and therapeutic failures. This study highlights the importance of whole genome sequencing and expanded real-time monitoring of diagnostic PCR assays during a pandemic.
Whole genome sequencing (WGS) based genetic surveillance is critical for tracking and forecasting pathogen evolution. Phylogenetic analyses can estimate the spatiotemporal mobility of pathogens between communities around the world over the course of a pandemic. This has been essential for identifying a circulating Variant Being Monitored (VBM), Variant of Concern (VoC), or Variant of Interest (VoI) of SARS-CoV-2, the causative agent of coronavirus disease 2019 . Since the emergence of SARS-CoV-2 in late 2019, sustained transmission among the human host has generated numerous variants with specific phenotypic attributes. Mutations that potentially confer increased transmission, pathogenicity (4), immune escape, or resistance to therapeutics can be identified when sequenced and interpreted by combining in silico inferences with in vitro functional, animal, and epidemiological studies to link genotypes to phenotypes. Without genomic surveillance, some mutations may go unnoticed if they do not impact diagnostic assays.
Polymerase chain reaction (PCR) assays are susceptible to false negatives that result from mutations that weaken primer annealing. For example, a mutation or mutations causing one primer set to fail within a multiplex assay-called a partial assay failure-is best exemplified by the Alpha (B.1.1.7) variant. This variant was discovered first in the UK in December, 2020 by the so-called S gene target failure (SGTF) (aka the S gene dropout) in a Thermofisher TaqPath PCR test kit that targets 3 different regions (ORF1ab, E gene, and S gene) of the viral genome (5). A 6-nt deletion ( 69-70 aa) resulted in the SGTF while the other two targets of the kit were positive for a given sample (5). After analyses determined the cause of the target failure was a new mutation, SGTF was used as a proxy for the presence of Alpha variant (B.1.1.7) in test samples. There are other examples of less severe assay failures in some commercial SARS-CoV-2 PCR tests (6,7). Such failures affect biosurveillance network responses and public health policy decisions for disease control, containment, and prevention as they form the basis for our understanding of how an outbreak is progressing. Thus, systematic and periodic assessment of real-time PCR assay performance is critical for maintaining assay specificity and decreasing false negatives (6-10).
The Global Initiative on Sharing Avian Influenza Data (GISAID) organization and Nextstrain project provide genomic data, metadata, and phylodynamic analysis for monitoring ongoing outbreaks (11,12). Continuous rapid sampling and sequencing of SARS-CoV-2 and sharing of data throughout the course of this pandemic provided a unique opportunity for visualizing signature erosion over time. Specifically, since January 20, 2020, we have used the PCR Signature Erosion Tool (PSET) to periodically evaluate assays in silico against all SARS-CoV-2 genome sequences from the GISAID EpiCoV TM database and share results on virological.org (13). PSET has been used previously to evaluate the impact of genomic drift on current and proposed assays for Ebolavirus and Mammarenavirus (14,15).
Here, we present our assessment of diagnostic PCR assays in silico by applying PSET to regularly collected SARS-CoV-2 sequences. Recent alignment-based studies have been conducted with comparable results (8,9). Our approach significantly increases the number of assays and subject sequences tested, breaks-down alignment rates over time by lineage, compares the mutation rate of assay target regions to the rest of the genome, and calculates the number of sequences producing false negatives for multiple assays.

PSET Analysis
In this study we tested the in silico performance of 43 PCR assays (Supplementary Table 1 Table 2) and sequences downloaded from the GISAID EpiCoV TM database on July 7, 2021 (11,13,(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26) (27,28). Additional follow-up analysis was performed on the B.1.1.529 VoC as it emerged later. The PSET definition for each assay target was based on a reference amplicon sequence with delimited primer and probe regions (Figure 1). Twenty nucleotides of additional sequence context outside the amplicon were also included to inform alignment at the 5 ′ and 3 ′ ends. Context and inter-primer regions were obtained via global-local alignment of the assay primers to the SARS-CoV-2 reference genome by running the glsearch36 program of the FASTA suite (29).
The first phase of PSET analysis queries the assay target definition against a BLAST+ database to search for matching subjects based on local alignment with the blastn program (30). Expansion of ambiguous DNA codes of the query is required for compatibility. For example, "GAWTAYA" has two, representing four expansions: "GAATACA, " "GATTACA, " "GAATATA, " and "GATTATA." Only the first permutation is queried. An additional step re-evaluates the BLAST+ identity statistics by replacing the expanded query with the original one to account for ambiguous base similarity. Subjects with ≥85% identity to the amplicon region are then extended to cover the query range and extracted.
The second phase queries the corresponding primers separately against the library of extracted sequences to search for matching subjects based on global-local alignment with the glsearch36 program, which is compatible with ambiguous DNA codes. Subjects with ≥90% identity are kept and aggregated by subject accession. A true positive (TP) is called if all primers aligned to the subject with the required identity and strand arrangement such that the primer would hypothetically amplify the target. The special case of a perfect true positive (PT) is called when there is 100% identity. Otherwise, a false negative (FN) is called. If the subject contains an N in any of the primer alignment regions, an N is appended to the category name (TPN and FNN). Doing so avoids removing low-quality, yet potentially informative sequences. An additional unknown (UNK) category is included to account for subjects that failed alignment during the first phase. In the case of near-neighbor analysis, a false positive (FP) or true negative (TN) is called when the taxonomy FIGURE 1 | PSET assay definition. The assay definition includes square brackets and parentheses to delimit the primer sequences. Note that a probe is optional and that any sequence outside of the amplicon region is considered context for alignment purposes. This example corresponds to the cdc_n1 assay.
of the subject differs from the assay target. FPs with 100% identity are perfect false positives (PF).

Variation Analysis
For variant analysis, the local GISAID EpiCoV TM database used for the PSET analysis was filtered to include only unambiguous DNA codes, resulting in a subset of 961,051 sequences. Single nucleotide variations (SNVs) and insertions or deletions (indels) were calculated by running the nucmer and show-snps programs of the MUMmer4 suite (46). Parameters for nucmer included NCBI GenBank accession no. NC_045512.2 (47) as the reference and flags to search the forward strand of each query (-f) with 28 threads (-t). The show-snps program transformed the resulting delta file into a tabular file (-T) of SNP/indel calls. An R script loaded the variants identified via MUMmer4 and split the data into two groups according to overlap with assay primer and probe target regions. The percentage of positions with ≥n mutations was calculated, where n was in the range [0, 500]. A two-sample Kolmogorov-Smirnov procedure tested the null hypothesis that the observed mutation percentages arose from the same distribution for the assay target and non-target genomic regions.

Omicron Analysis
The Omicron (B.1.1.529) wave emerged after the initial analysis was completed. Accordingly, an updated analysis targeted this VoC. On February 16, 2022, genomes were downloaded manually from the EpiCoV TM search page. Filters were enabled to include complete sequences with high coverage and collection dates. The PSET workflow analyzed the final set of 48,358 sequences (Supplementary Table 3) as previously described. Figure 2 shows the target region of each assay with respect to the reference SARS-CoV-2 genome (NCBI GenBank accession no. NC_045512.2), SNVs, and indels (47). Mutations with <1% prevalence (∼99.62% of all positions with observed mutations) were excluded from the figure but remain listed in Supplementary Table 4. The figure also plots the assay target regions on the genome. Plotting all the regions together with respect to mutations helps visualize signature erosion due to genetic drift.

Assay Regions and Genomic Variation
On the other hand, when the assay target regions are compared specifically against the non-target regions, assay target regions with mutations are present in higher percentage of genomes (Figure 3). The two-sample Kolmogorov-Smirnov test obtained a p value of 2.2E−16, rejecting the null hypothesis at a significance level of 0.05 in favor of the two-sided alternative.

Assay Alignment Confusion Matrix
The assays were able to correctly detect the vast majority of the 1,690,689 genome sequences according to the in silico PSET analysis. Of the 43 assays, 34 aligned with 100% identity to over 1.6 million subjects. Otherwise, TP rates exceeded 98%, except for the Young-S and China_N assays, which exhibited high FN rates at ∼47% and ∼59% respectively. The TP percentage was calculated as the sum of PT, TP, and TPN divided by the total number of subjects. Alignments exceeded identity threshold for most corresponding subject amplicon sequences, only failing on average ∼0.74% (UNK%) of the time. Table 1 and Supplementary Table 7 include a full breakdown by assay. In the separate Omicron study, most of the assays (40 / 43) were able to detect sequences with a TP rate of ≥96%, with Young_S, Thailand_WH-NIC_N, and China_N assays having reached TP rates of 75.64%, 2.98%, and 0.34% respectively (Supplementary Table 8). We also looked at the specificity of the assays against near-neighbor sequences (Supplementary Table 9), where 8 showed significant false positives (FP) and perfect FPs (PFs) (those with 100% identity).

Assay Alignment Identity Over Time
Alignment identity of all assays against sequences of different lineages of SARS-CoV-2 was assessed from 2020-03-15 to 2021-07-05. Figure 4 depicts heat maps of the TP rate over time for each VoC/VoI with a corresponding line graph of the cumulative log-total number of sequences. Only GISAID sequences with collection date metadata specifying year, month, and day were included in the heatmap. The graph reveals sudden changes in assay alignment identity with respect to variant abundance. Alignment identity remained constant for most assays, with perfect or near-perfect TP rates. However, some interesting patterns were observed. Some assays targeting the S and N genes performed poorly. For example, Young-S exhibited low TP rates within the B.  In both cases, the WHO called last. We found no evidence of a call from the ECDC in these cases. The B.1.429 case shows the initial call occurring during a period of low alignment identity for the Young-S assay and the beginning of a decline for the Japan_NIID_WF-1_Seq_F519 assay. Despite the varying times at which these lineages were identified, when looking at the cumulative increase of the GISAID sequences for each, it is readily apparent that they initially occurred at very low proportions. This increase, and the resulting VoC/VoI designations, underscore the importance of continual monitoring from multiple organizations offering different perspectives and criteria with respect to spatiotemporal trends.

In silico False Negative Distribution by Assay
As assay targets were distributed throughout the genome, we were also interested in seeing whether emerging variants possessed mutations that resulted in more than one assay  Table SC shows that the majority (97%) of sequences producing a single FN were observed to originate from the China_N (87%) or Young-S (10%) assays. The majority (765,813) of sequences produced FNs in two different assays, China_N and Young-S, while 643,876 sequences did not produce an FN in any of the 43 assays. We observed a substantial drop-off in the number of sequences which caused FNs in 3 or more assays. This suggests that it is unlikely for a sequence to produce FNs in more than 2 assays and that a multiplex panel with any of the other 41 assays evaluated here would likely perform successfully.

Signature Erosion Analysis
The rapid increase of novel SARS-CoV-2 variants raises concerns regarding the efficacy of PCR-based diagnostic assays. This study assesses the potential signature erosion of the current COVID-19 real-time diagnostic assays using the GISAID sequence database. We find that 41 out of the 43 assays continue to perform very well throughout the observed timeframe even as new variants arose. Additional analysis completed after the emergence of Omicron revealed that 40 of 43 assays continued to perform well.
Variants containing mutations within the target regions for primer and probe hybridization were observed to have a substantial effect on predicted assay alignment identity. The characteristic S: 69/70 mutation of the Alpha variant occurs within the forward primer of the Young-S assay target, while mutations N:R203K and N:G204R occurred within the borders of the China_N assay target, resulting in PSET calling significant FN rates of those assays. The China_N assay failed completely for the lineages B.1.1.7, P1, and P2. This may be due to a "GGG" to "AAC" mutation in the forward primer, as observed in previous studies (8,9). The 3-nt substitution mutation could potentially reduce the China_N forward primer binding affinity, which could significantly increase the FN rate, especially since it occurs at the 5 ′ -end. The B.1.1.529 lineage contains an N: 31/33 deletion (48) which falls inside the borders of the Thailand_WH-NIC_N assay and results in PSET predicting significant FN rates. This may be a similar situation as described above where primer binding affinity is affected by this type of mutation. Various point mutations at the primer binding motifs are observed in all the lineages, which did not affect overall performance according to evaluation of in silico alignments. The exception was for lineage B.1.258, which bears a complementary sequence to the China_N primer resulting in perfect alignment. The China_N and Young-S assays experience 56× and 44× the FN rate of all other assays combined, respectively. In almost all assays, there is a percentage of mutations in one of the targets of the assay.

Recent Studies and Future Direction
The COVID-19 pandemic has demonstrated the need for the in silico monitoring of PCR-based diagnostic assay performance. Furthermore, given that surveillance network

ID
Assay From  To  Gene  PT  TP  TPN  FN  FNN  UNK   1  Japan_NIID_WH-1_F501  483  896  ORF1ab  1,657,346  20,369  178  89  116  12,591   2  Japan_NIID_WH-1_F509  491  837  ORF1ab  1,640,106  37,112  429  201  268  12,573   3  Japan_NIID_WH-1_Seq_F519  501  823  ORF1ab  1,634,109  28,526  433  6, data are widely available to aid in quick-turnaround research and development of replacement assays, in silico monitoring has become an important early warning indicator (10). Our approach is consistent with recent alignment-based studies by confirming the overall high assay target sequence identity and detection of the China_N mutation (8,9). We also observed potential early signs of signature erosion prior to official variant calls. Future research can further elucidate validation of in vitro and in silico predictions while considering the clinical and public health relevance. For example, the PSET algorithm could incorporate mismatch position, since mismatches near the 3 ′ end of the probe sequence potentially have a greater impact on assay performance (49,50). Additional wet lab experiments can help systematically select alignment thresholds and motivate algorithmic refinements.

CONCLUSION
This study highlights the consequences of mutations in SARS-CoV-2 genomes on PCR-based diagnostic assays. The importance of real-time monitoring of molecular assay alignment in silico is highlighted by the discovery of the Alpha variant. The failure of one primer/probe combination within the multiplex assay targeting a region of a deletion in this variant underscored the need for extensive sequencebased surveillance. Sustained transmission and proliferation of viruses such as SARS-CoV-2 in a global pandemic leads to rapid evolution and accumulation of mutations that confer advantageous phenotypes, such as potentially evading diagnostics, therapeutics, and vaccines. Real-time in silico monitoring of assay signature erosion allows for the redesign and refinement of diagnostic assays to address assay and medical countermeasure failures to avoid dire failures of medical countermeasures.

ACKNOWLEDGMENTS
We gratefully acknowledge the work of all authors who generated and submitted the numerous sequences to GISAID's EpiCoV TM Database on which this study is based on. We also acknowledge the contributions and data provided by the GISAID organization. Special thanks to Dr. Yan Yi for reviewing and advice on the statistical analysis.