SARS-CoV-2 Sequence Characteristics of COVID-19 Persistence and Reinfection

Data search and selection criteria

We conducted a systematic literature review in PubMed through March 8, 2021 for cases of persistent COVID-19 using the search term “((covid or sars-CoV-2) AND (persistent or persistence or prolonged)) AND (sequence or evolution)”. A search for COVID-19 reinfection reports was made using the terms “(covid or sars-CoV-2) AND (reinfection)”. Both peer-reviewed and preprint results were evaluated. We used the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) for reviewing literature and for reporting search results. Additional preprints that appeared through Google search and that met our criteria were also included. For cases of reinfection, papers were included if the authors described it as a case of reinfection diagnosed >30 days after the initial infection and if whole genome SARS-CoV-2 sequences or sites of mutations relative to a reference sequence (e.g., Wuhan-Hu-1) from both infection time-points were available. Of the 291 results from the search, 14 articles met the inclusion criteria and were included in the present report along with 2 additional preprints that were identified (Supplemental Figure 1A).

• Download figure
• Open in new tab

Figure 1.

Maximum-likelihood phylogenetic tree of sequences from persistent COVID-19 cases (Pe1-Pe7), COVID-19 reinfection cases (Re1-Re16), the variants of concern B.1.1.7 and B.1.351, and globally sampled sequences from GISAID.

Persistent cases were included if the authors described it as a case of persistent COVID-19 infection and if longitudinal whole genome SARS-CoV-2 sequences were available. The search returned 129 results, 7 of which met the inclusion criteria and were included in the present report along with one other preprint (Supplemental Figure 1B). Only sequences obtained directly from patient respiratory tract samples were included in our analysis to exclude the possibility of sequence changes during the ex vivo culture process. Three cases were excluded due to uncertainty in their classification as either reinfection or persistent infection cases (Supplemental Methods, Supplemental Table 1, Supplemental Figure 2).

• Download figure
• Open in new tab

Figure 2.

Comparison of viral sequences from reinfection cases. (A) Circos plot showing location of nucleotide changes in the reinfecting sequence relative to the initial infection sequence for each of the 20 cases. Inner ring indicates nucleotide position in kilobases. Synonymous changes are in green, nonsynonymous changes in orange, deletions in black. (B) Nucleotide (NT) substitution frequency pooled across all reinfection cases for each SARS-CoV-2 gene. Dashed line indicates global substitution frequency across the whole genome. Substitution frequency for each gene was compared to the substitution frequency in the rest of the genome using a Fisher’s exact test. P-values were corrected for multiple comparisons using the Bonferroni correction. * <0.05, ** <0.01 and ***<0.001 (C) Nucleotide changes in the second infection relative to the first infection by clinical disease severity. Mutations shown for the whole genome and S gene. P=0.67, Mann Whitney test. (D) Circos plot showing location of nucleotide mutations from the second infection relative to other viruses circulating at the same time in the same geographic region. Only rare mutations present in <1% of contemporaneous community sequences are shown. (E) Number of rare nucleotide polymorphisms at each time point relative to circulating sequences in the community. P=0.26, Wilcoxon matched-pairs signed rank test. ORF: open reading frame, S: Spike, E: Envelope, M: Membrane, N: Nucleocapsid.

Sequences were analyzed for mutations using NextClade (https://clades.nextstrain.org/) and snp-sites (https://github.com/sanger-pathogens/snp-sites). The degree of reinfection severity, either more or less severe compared to the first infection, was classified based on an explicit determination by the authors of each article or by comparing symptoms, duration of illness, and hospitalization status between both episodes.

Sequencing dataset compilation and phylogenetic tree construction

The sequencing dataset contained a total of 262 globally representative SARS-CoV-2 genomes selected from GISAID and sequences from the reinfection and persistence cases (Supplemental Methods; Supplemental Data 1). The sampled sequences were chosen to be representative of global sequence diversity throughout the time course of the pandemic. Sequences of variants of concern B.1.1.7 and B.1.351 were also included. Nucleotide sequence alignment was performed using MAFFT (Multiple Alignment using Fast Fourier Transform) [12]. Best-fit nucleotide substitution was calculated using model selection followed by maximum likelihood (ML) phylogenetic tree construction using IQ-Tree with 1000-bootstrap replicates [12].

Mutation analysis

For reinfection cases, mutations were determined in two ways. First, nucleotide and amino acid changes were identified for the reinfection sequences relative to the first infection sequence. The frequency of nucleotide or amino acid changes within each gene was compared to the frequency of changes in the remainder of the genome by Fisher’s exact tests with a Bonferroni correction (for multiple comparisons). The relationship between disease severity and number of nucleotide or amino acid changes in the genome was assessed using a Mann-Whitney test. Second, to identify unique characteristics of reinfecting viruses, each of the first and reinfection sequences were compared to circulating sequences in the community as defined by the same NextStrain clade sampled within one month obtained from the same geographic location uploaded to GISAID (https://www.gisaid.org/; Supplemental Table 2, Supplemental Methods, Supplemental Data 2). Rare mutations were determined as polymorphisms that were present only in the reinfecting sequence (not the initial variant) and found in less than 1% of contemporaneous community sequences. Mutation locations are graphically represented in Circos plots [13].

For persistent infections, sequence changes were assessed at two time intervals: before or after convalescent plasma or monoclonal antibody treatment. Sequences sampled before convalescent plasma or antibody treatment were compared to the first sequence sampled. For sequences sampled after convalescent plasma or antibody treatment, sequence changes (both nucleotide and amino acid) were determined relative to the last pre-treatment sequence. Linear regression was used to estimate the rate of viral changes between two intervals. The slope of the trendline was compared to the latest global clock rate (March 29, 2021) as estimated by NextStrain (https://nextstrain.org/ncov/global/).

Time-measured phylogenetic analysis

The temporal signal of the ML tree was examined in TempEst [14] regressing on root-to-tip divergence, and outliers were inspected in the distribution of residuals. A high degree of clock-like behavior in the whole dataset was observed (R² = 0.721) in root-to-tip regression analysis with the slope rate as 8.26E-4 and the rough ancestral time of the sample was calculated as 2019.84. This suggests that the whole dataset has a realistic temporal signal and it is appropriate for an estimation of temporal parameters. No outliers were found in this sample. To further examine the temporal signal in the sequences from persistent patients (especially these with > 2 sequences), separate root-to-tip regression analysis also supported temporal signal for a time-measured phylogeny. To compare the evolutionary rates between the reported persistent infections and the general population infections, time-measured phylogenetic reconstruction was conducted in Bayesian Evolutionary Analysis Sampling Trees (BEAST) v1.10.4 [15]. Nine partitions, including eight persistent patients and the global sequences, were used as separate groups of taxa, to estimate separate evolutionary rates. Due to large uncertainties with small samples, persistent patients with only two viral sequences were excluded from this analysis. A general time reversible (GTR) model was applied with gamma-distributed rate variations among sites. A lognormal relaxed molecular clock was used with an initial mean of 0.0008 and a uniform prior ranging from 0.0 to 1.0. A logistic growth tree prior was applied. Four independent Bayesian Markov Chain Monte Carlo (MCMC) chains of 100 million generations were performed with a sampling step every 10,000 generations to yield 10,000 trees per run. To ensure a sufficient effective sample size ESS > 200, the convergence of three runs was diagnosed in Tracer v 1.7.1 (http://tree.bio.ed.ac.uk/software/tracer/) for all parameters. LogCombiner v1.10.4 as part of the BEAST software package was used to combine the multiple runs to generate log and tree files after appropriate removal of the burn-in from each MCMC chain. The comparison of the evolutionary rates from the combined log file is analyzed and visualized in R v4.0.2 (https://www.r-project.org/).

Statistical analysis

Nonparametric Wilcoxon rank sum or matched pairs signed rank tests were used to compare the number of amino acid changes between sequences. Statistical analyses were performed using GraphPad Prism 9 (GraphPad Software, San Diego, CA).

Sequence analysis of reinfection cases

A total of twenty cases from sixteen reports were included in this analysis (Table 1) [2-7, 16-25]. A broad range of age groups were represented and 90% were under the age of 70 years. Most (80%) of the cases had no reported comorbidities and while one patient had diabetes and end-stage renal disease, none had high-level immunosuppression. The interval between diagnosis of the first infection and the second infection ranged from 44 days to 282 days with a median of 113.5 days. Five patients had more severe illness during the second infection, while six had less severe symptoms on reinfection, including two who were asymptomatic on reinfection. Two cases were asymptomatic in both infections, five cases reported the same severity for both infections and no information on infection severity was available for two cases (Table 1). Six cases reported reinfection with a virus from the same clade.

Phylogenetic analysis demonstrated distinct branching for the two sequences in each of the reinfection cases, corroborating results discussed in the original reports (Figure 1). We compared nucleotide and amino acid changes in the reinfecting viral sequence compared to the initial sequence and found a median of 17.5 nucleotide changes (range 9-37) and 9 amino-acid changes (range 6-24) compared to the original sequence (Figure 2A). The nucleotide changes between the initial and reinfecting sequences were distributed across the SARS-CoV-2 genome, with significantly higher frequencies of changes in ORF8 (P<0.001) and N (P=0.001) (Figure 2B). A similar pattern was observed with amino acid changes (Supplemental Figure 3A). All but two reinfection cases had at least one substitution or deletion in the S gene (Supplemental Table 3). Next, we assessed whether reinfection with a more divergent second virus resulted in more severe disease. We found no significant differences in the number of nucleotide or amino acid changes in the reinfecting virus compared to the original viral variant when categorized by the severity of the reinfection (Figure 2C; Supplemental Figure 3B). Both the initial and reinfecting SARS-CoV-2 variants were similar to the sequences circulating in the community at the time of reinfection. The initial infecting variant harbored a median of only 2 rare nucleotide mutations compared to contemporaneous circulating variants in the community and the reinfecting variant contained a median of only 1 rare nucleotide mutation (Figure 2D-E; Supplemental Figure 3C).

• Download figure
• Open in new tab

Figure 3.

SARS-CoV-2 mutation location and evolutionary rate in the persistent COVID-19 cases. (A) Nucleotide (NT) substitution frequency pooled across all persistent cases for each SARS-CoV-2 gene. Dashed line indicates global substitution frequency across the whole genome. Substitution frequency for each gene was compared to the substitution frequency in the rest of the genome using a Fisher’s exact test. P-values were corrected for multiple comparisons using the Bonferroni correction. * <0.05, ** <0.01 and ***<0.001 (B) Nucleotide changes in samples taken prior to convalescent plasma or monoclonal antibody treatment relative to first sampled sequence in each persistently infected patient. Regression line and 95% confidence bands are shown. Purple dash-dotted line is global rate estimate obtained from NextStrain. (C) Nucleotide changes in samples taken after convalescent plasma or monoclonal antibody treatment relative to last sample taken prior to treatment in each persistently infected patient. Regression line and 95% confidence bands are shown. (D) Substitution rate (nucleotide substitutions per site per year) of sampled global SARS-CoV-2 sequences relative to persistent patients based on Markov chain Monte Carlo time-measured phylogenetic reconstruction. Box plots show median and interquartile ranges of estimated substitution rates. The mean, median, and 95% highest posterior density (HPD) interval can be found in Supplemental Table 4.

Sequence analysis of persistent COVID-19 cases

A total of nine cases from seven reports describing persistent infection were retrieved from our literature search. Of these nine cases, all but one had B cell immunodeficiency [8-10, 26-29]. Four were treated with B cell-depleting therapy for lymphoma or autoimmune disorders, while four had B cell lymphomas treated with chemotherapy (Table 2). One patient had advanced HIV infection with a CD4+ count of 0 cells/mm³ and diminished CD19+ cell counts. The median length of infection was 154 days and 33% of the cases ended in death. One patient had asymptomatic disease throughout [9]. Four patients were treated with convalescent plasma at least once during their illness [9, 10, 26, 28], and one patient was treated with the monoclonal antibodies casirivimab and imdevimab [8].

Phylogenetic analysis revealed that, for each of the nine patients, sequences formed a distinct cluster, confirming what was found in the original reports (Figure 1). New mutations emerging over time were detected in all of the persistent COVID-19 patients with further changes identified after treatment with convalescent plasma or monoclonal antibodies (Supplemental Figure 4). Mutations occurred with significantly higher frequency in S (P<0.001) and ORF7a (P=0.02) and lower frequency in ORF1a (P=0.02) (Figure 3A; Supplemental Figure 5A). The rate of viral evolution was plotted for each patient both for the interval before and after convalescent plasma/antibody treatment. Before antiviral treatment, the rate of sequence changes over time appeared faster than the NextStrain estimate for the global rate of SARS-CoV-2 evolution (dotted purple line) (Figure 3B; Supplemental Figure 5B). Treatment with convalescent plasma or antibody cocktail treatment was insufficient to halt intra-host viral evolution (Figure 3C; Supplemental Figure 5C).

We also performed time-measured phylogenetic reconstruction with the pre-treatment persistent sequences to compare the rate of intra-host viral evolution in persistent COVID-19 to the rate of community-driven evolution. This analysis provided further evidence that SARS-CoV-2 evolution appeared faster in these persistent infection individuals compared to the rate in the general public population, though substantial uncertainties are shown in these estimates given the limited sequence sampling in each patient (Figure 3D; Supplemental Table 4).