Abstract
The COVID-19 pandemic has initiated an unprecedented worldwide effort to characterize its evolution through the mapping of mutations of the coronavirus SARS-CoV-2. The early identification of mutations that could confer adaptive advantages to the virus, such as higher infectivity or immune evasion, is of paramount importance. However, the large number of currently available genomes precludes the efficient use of phylogeny-based methods. Here we present CoVtRec, a fast and scalable Topological Data Analysis approach for the surveillance of emerging adaptive mutations in large genomic datasets. Our method overcomes limitations of state-of-the-art phylogeny-based approaches by quantifying the potential adaptiveness of mutations merely by their topological footprint in the genome alignment, without resorting to the reconstruction of a single optimal phylogenetic tree. Analyzing millions of SARS-CoV-2 genomes from GISAID, we find a correlation between topological signals and adaptation to the human host. By leveraging the stratification by time in sequence data, our method enables the high-resolution longitudinal analysis of topological signals of adaptation. We characterize the convergent evolution of the coronavirus throughout the whole pandemic to date, report on emerging potentially adaptive mutations, and pinpoint mutations in Variants of Concern that are likely associated with positive selection. Our approach can improve the surveillance of mutations of concern and guide experimental studies.
Introduction
The COVID-19 pandemic, caused by the coronavirus SARS-CoV-2, has led to millions of lost human lives and devastating economic impact worldwide. As the virus continues to spread through the world, it is acquiring new mutations in its genome, and although most mutations will be deleterious or neutral, a few of them could be advantageous for the virus, for instance, by increasing its infectivity or by helping it to avoid the immune system. As more people develop immune protection by previous viral infections or through vaccination, it is important to rapidly and effectively identify mutations that could confer the virus some fitness advantage [1, 2].
One approach to identify adaptive mutations consists of experimentally mutating many positions and testing the effect on certain phenotypes, like binding to the human receptor or immune evasion [3–6]. However, experimental approaches are limited by the vast number of possible variations, especially when different genetic backgrounds are considered [7].
A more data-driven approach, which solely relies on the genomic information of the virus, is to look for mutations in a particular genomic locus that occur multiple times. In this approach, one usually reconstructs an estimated phylogenetic tree and identifies mutations that occur in independent branches [8–16]. If a mutation gives some sort of advantage to the virus, we expect it to occur in several places in the phylogeny independently and its prevalence to increase with time. For instance, the D614G mutation in the Spike gene was identified early in the pandemic and is now found in virtually all virus isolates [17].
In the COVID-19 pandemic, an unprecedented worldwide effort to sequence viruses resulted in a growing number of millions of genomes available to the scientific community [18]. Ideally, one would like to leverage all this genomic information at real-time to rapidly report the emergence of potential mutations of concern [1, 2]. Several phylogeny-based tools and methods are available [7, 19–31]. Phylogenetic approaches, however, become daunting as the number of sequences increases, and are computationally challenging when the number of genomes exceeds the tens of thousands [13, 19, 32]. In addition, homoplasies caused by the independent reemergence of mutations, as well as recombination events, confuse the generation of phylogenetic trees [33, 34]. It has moreover been observed that certain methods for constructing a single optimal phylogeny for SARS-CoV-2 can be problematic as the number of isolates is large while the genetic diversity is low [32].
Here we establish a novel method, based on Topological Data Analysis (TDA) and implemented in the CoVtRec pipeline (see https://tdalife.github.io/covtrec), that can efficiently identify emerging potentially adaptive mutations without the need to choose a single optimal phylogenetic tree in a vast range of equally plausible tree reconstructions. It systematically detects convergent events in viral evolution merely by their topological footprint in the genome alignment, thus overcoming limitations of current phylogenetic inference techniques. Thanks to our use of highly optimized algorithms, it easily scales to hundreds of thousands of distinct genomes. Our approach has the ability to map the dynamics of topological signals of adaptation at daily resolution over any period of time, which enables not only the real-time monitoring but also the retrospective assessment of adaptive processes in the evolution of the coronavirus. We evaluate how our new method compares, both in performance and in results, with state-of-the-art phylogeny-based methods for the construction and analysis of phylogenetic trees, such as UShER [29], IQTree [35] and HomoplasyFinder [9], and for the assessment of fitness effects of mutations as developed by Bloom & Neher [7].
We perform CoVtRec analyses of large alignments with millions of genomic sequences shared via GISAID, the global data science initiative [18, 36], in three different phases of the pandemic: (i) in February 2021, retrospectively monitoring the convergent evolution of the virus from the beginning of the pandemic in a largely immunologically naive population; (ii) in January 2022, in the presence of increasing host immunity and dominant variants, with the Omicron variant on the rise; and (iii) in March 2023, in order to obtain a comprehensive picture of the dynamics of adaptation over the whole course of the pandemic. Invoking published experimental and epidemiological data into our analysis, we investigate the correlation between topological signals and host adaptation. As we show, our method can identify adaptive mutations at an early stage, often several months before they become recognizable by their prevalence in the population. We characterize potentially adaptive mutations in different phases of the pandemic, and identify mutations in Variants of Concern that are likely due to convergent evolution.
Results
Quantification of topological recurrence
Chan et al. [37] initiated the use of persistent homology, a method from Topological Data Analysis, to extract global evolutionary features from large genomic datasets. This method detects topological cycles in the dataset, which correspond to reticulate events in the phylogeny and may cause phylogenetic inference methods to produce ambiguous tree topologies. All information is compiled into a stable and unbiased descriptor known as a persistence barcode (see Figure 1, Figure 2 and Methods).
(A) and (B) show possible evolutionary histories on the level of individuals in the example of a genome with only two nucleotides. The coloring of the edges corresponds to the acquisition of a specific mutation, while the coloring of the nodes represents individuals carrying this mutation. Convergent evolution (A) or recombination (B) leads to the presence of four alleles for which there is no single consistent phylogeny (four-gamete test). On the genomic level, genetically identical individuals cannot be distinguished and incompatible phylogenies are represented by a topological cycle in the corresponding phylogenetic network (C).
(A) Persistent homology detects reticulate events in viral evolution by means of a persistence barcode. Each bar in the barcode (A2) corresponds to a topological cycle in the reticulate phylogeny (A1). Bars at small genetic distance scales are expected to correspond mainly to homoplasies, while recombination events typically produce topological features at larger scales. (B) SNV cycles are topological cycles in the reticulate phylogeny for which adjacent sequences differ by single nucleotide variations (SNV) only. Under the assumption of single substitutions per site, any SNV in each such cycle appears twice, independently of each other, and distributed across both possible lineages. The topological recurrence index (tRI) of a specific mutation is the total number of SNV cycles in the reticulate phylogeny in which this mutation is acquired. In the displayed example phylogeny, the red edges indicate the acquisition of a specific mutation, while the red nodes represent viruses carrying this mutation. The mutation is acquired in two SNV cycles (shaded in blue) and therefore has a tRI of 2.
There are several scenarios that can lead to reticulate events. For instance, if a genome of an organism imports genetic material from a different genome, in lateral gene transfer for instance, we will observe that parts of the newly generated genome resemble the parent, while others resemble the genome of the organism that exported that material. Recombination and reassortments are common phenomena observed in viruses where two parental strains co-infect the same host cell generating a new virus containing genetic material from both parental strains. But similarity between genomes can also be generated at smaller scales, when the same mutation occurs independently twice, making the two strains more similar than expected. Persistent homology captures all these events, and also the scale of the events. Although in general it requires care to infer the biological origin of a given topological cycle, in viral evolution one expects bars at small scales to correspond mostly to homoplasies, while well-supported recombination events typically produce topological features at larger scales, as entire blocks of genetic material are exchanged in the process [37–39].
Here we define a novel index of recurrence that is based on persistent homology and does not rely on a possibly ambiguous tree reconstruction. We use a specifically designed algorithm, implemented in Ripser [40], that associates to relevant bars in the persistence barcode explicit SNV cycles, given by a series of isolates that approximates all evolutionary steps as faithfully as possible in terms of single nucleotide variations (SNV). We define the topological recurrence index (tRI) of a specific mutation as the total number of SNV cycles in the reticulate phylogeny that contain an edge corresponding to this mutation. This index provides a lower bound for the number of independent occurrences of the mutation in the phylogeny and is therefore a measure for convergent evolution (see Figure 2 and Methods).
Our method enables the temporal mapping of topological signals of recurrence at daily or monthly resolution by leveraging the stratification by time in genomic datasets given by the collection dates of genomic isolates (see Supplementary Information Figure 13, Figure 5, Figure 8 and Methods).
Comparison with phylogeny-based methods
Standard phylogenetic methods for the detection of convergent evolution are based on the reconstruction of a phylogenetic tree and therefore have an unfavorable scaling, due to the rapid growth of the number of trees representing evolutionary histories that are compatible with the observations [32]. Persistent homology provides a new approach by measuring convergent evolution purely in terms of topological cycles, without the need to construct any phylogenetic trees. It enables a rapid and scalable analysis of large alignments with hundreds of thousands of distinct sequences. A performance analysis showed that when applied to sample alignments of up to 5,000 SARS-CoV-2 sequences, the topological method was faster than state-of-the-art phylogeny-based methods such as UShER and HomoplasyFinder [9, 29, 35] (see Figure 3). Moreover, the topological recurrence analysis of current alignments with up to 13 million SARS-CoV-2 sequences was accomplished in less than a day (see Methods).
(A) Basic runtime comparison between Topological Data Analysis (TDA)-based methods and phylogeny-based methods for random samples of up to 5,000 SARS-CoV-2 genomes. We used IQTree [35] and UShER [29] to reconstruct phylogenetic trees. The subsequent homoplasy analysis was performed with HomoplasyFinder [9]. For the TDA-approach we used Hammingdist [41] to generate genetic distance matrices and Ripser [40] for the subsequent computation of persistence barcodes (see Methods). (B) Runtime and memory usage for the computation of persistence bar-codes with Ripser [40] for monthly sub-alignments of the GISAID alignment with 161,024 genetically distinct whole genome sequences covering the first year of the pandemic (see Methods).
The UShER phylogenetic tree [42] has recently been used to estimate fitness effects of mutations in the evolution of SARS-CoV-2 by means of a “fitness index” [7]. We compared results from this fitness analysis with results from our topological recurrence analysis based on GI-SAID data as of March 2023. For Spike gene amino acid changes, we found a strong correlation between positive topological recurrence index (tRI) and positive “fitness index” (Fisher’s exact test, p < 10e 10) with a Spearman correlation coefficient of 0.63 (p < 10e 99) (see Supplementary Information Figure 12).
Topological analysis of the first year of the pandemic
We analyzed topological signals for convergent evolution of the coronavirus SARS-CoV-2, both across the whole genome and on the Spike gene, during the first year of the pandemic from its beginning in December 2019 until February 2021. To that end, we performed a topological recurrence analysis for a curated alignment of 303,651 high-quality SARS-CoV-2 whole genome sequences from GISAID [18, 36] (see Methods). The resulting persistence barcode features 2,899 bars, 58% of which (corresponding to 1,684 SNV cycles) concentrate at small genetic distance scales 2 and are therefore expected to be associated mainly with homoplasic events (see Supplementary Information Figure 1). We did not observe bars at genetic distances that are associated with well-supported recombination events, which is compatible with previous reports that recombination did not play a substantial role in the evolution during the first year of the pandemic [30, 43].
The topological recurrence index (tRI) is designed in such a way that it naturally excludes recombination events. It is defined in terms of SNV cycles built exclusively from single nucleotide variations (see Figure 2 and Methods), while recombination events in SARS-CoV-2 typically involve the exchange of larger parts of the genome. This has also been documented for later phases of the pandemic, and in particular for circulating recombinant forms during the Omicron waves [44, 45]. It is therefore unlikely that recombination events significantly contribute to tRI signals in SARS-CoV-2 evolution.
In large genomic datasets, topological cycles in persistent homology may arise randomly, causing a certain amount of noise in the topological recurrence analysis. However, from simulations we inferred that at least 60% of the topological cycles found in the GISAID dataset must be real features due to increased mutation probabilities, selection, recombination, or sequencing errors (see Supplementary Information Figure 9 and Methods). We implemented several filtration rules to make sure that sequencing errors have no significant effect on results of the topological recurrence analysis (see Methods).
We detected 401 non-synonymous and 299 synonymous mutations across the whole genome showing significant tRI signal (see Figure 4 and Supplementary Information Table 1). Here signals with tRI ≥ 2 are statistically significant (p < 0.05; see Methods).
Topological recurrence analysis for the whole genome and the Spike gene from December 2019 until February 2021. The tables in (A1) and (B1) list the topological recurrence index (tRI) and the prevalence for mutations with strongest topological signal. In each table the last column reports filtration rules recommended in Turakhia et al. [33, 46] and De Maio et al. [47] for highly homoplasic sites, problematic sites and potential artifacts that are normally removed or masked in SARS-CoV-2 phylogenetic analyses. For a complete list see Supplementary Information Table 1 and Supplementary Information Table 2. (A2) Histogram showing the distribution of topological signals across the whole genome. In every region of the genome, reticulate events play a crucial role in the evolution of the virus. (B2) Histogram showing the distribution of topological signals across the Spike gene. There is an accumulation of topologically recurrent mutations in the S1 subunit and the S1-S2 cleavage region, as well as in the signal peptide at the beginning of the Spike gene.
Several of the mutations with strong tRI signal, such as G11083T, C21575T, C15324T and C3037T, are known to be highly homoplasic for various reasons and cause stability issues in the construction of phylogenetic trees [11, 32, 33, 47]. While the corresponding sites are typically masked in phylogenetic analyses, this is not necessary in our Topological Data Analysis approach as the computation of the topological recurrence index (tRI) does not rely on the construction of phylogenetic trees. We prefer to take an unbiased approach and include additional information on highly homoplasic and potentially problematic sites only in the downstream analysis of mutations with significant topological signal. Furthermore, we augment our tables of tRI scores by common masking schemes suggested in Turakhia et al. [33, 46] and De Maio et al. [47] (see Figure 4 and Supplementary Information Table 1).
Host RNA editing processes have been reported to be a factor in SARS-CoV-2 mutagenesis [48–50]. We investigated the role of innate immune editing in the topological analysis of the first year of the pandemic (see Supplementary Information Figure 7). Among mutations on the whole genome flagged by a significant tRI signal, we observed a preponderance of C>T transitions with a strong bias towards the positive-sense strand: 58% of the non-synonymous mutations, and 73% of the synonymous mutations were C-to-T changes. A cosine similarity analysis between the trinucleotide signature of mutations flagged with significant tRI and COSMIC signatures [51] revealed a high resemblance with signature SBS 2 that is attributed to the action of the APOBEC family of cytidine deaminases (cosine similarity of 0.39 for non-synonymous mutations and of 0.44 for synonymous mutations; see Methods). These findings suggest that topological signals are also driven by C>T hypermutation caused by APOBEC-mediated human host responses [52–54].
Notably, we observed that 43% (299 out of 700) of the topologically recurrent single nucleotide variations on the whole genome were synonymous. The excess and imbalance in frequencies of C-to-T transitions (73% synonymous vs. 58% non-synonymous) among topologically recurrent mutations (see Supplementary Information Figure 7) suggests that the action of APOBEC enzymes is one explanation for this observation. Moreover, it is known that synonymous mutations are not always neutral or almost neutral [55], and it is likely that the topological recurrence index also captures signals of the adaptation of the SARS-CoV-2 codon usage to the human host via codon optimization and deoptimization, processes which can affect SARS-CoV-2 protein translation efficiency [56–59]. A substantial number of fitness increasing synonymous mutations was also reported in the phylogeny-based assessment of fitness effects by Bloom & Neher [7, 60]. These observations suggest that the assessment of the adaptiveness of topologically recurrent synonymous mutations requires care, for example by invoking other indicators of adaptation. Our downstream analysis of the potential adaptiveness of topologically recurrent mutations in the present work focuses on protein changing substitutions as they are more interpretable in terms of protein selection.
The topological recurrence analysis can also be carried out for single genes, based on alignments of appropriately truncated genome sequences (see Methods). In this case, evolutionary processes outside the specific gene are ignored, which typically leads to the creation of more topological cycles and hence to a more detailed picture of the ongoing convergent evolution in the respective gene. For the Spike gene alignment, we found a total of 322 non-synonymous and 196 synonymous Spike mutations with significant tRI signal (see Figure 4, Supplementary Information Table 2 and Methods). Here signals with tRI 8 are statistically significant (p < 0.05; see Methods). We observed a distinct accumulation of topological signals in the S1 subunit, which is associated with host receptor recognition and contains epitopes for antibody binding [61–65], and in the Spike protein signal peptide as well as in the S1-S2 cleavage region [66] (see Figure 4).
We found a strong topological signal for the mutation S:D614G. A comparative time series analysis (tRI vs. prevalence) covering the first year of the pandemic revealed a steady increase in both tRI and prevalence until the new variant eventually superseded the wild type in the population (see Figure 5 and Methods). This mutation is known to increase transmissibility [17, 67] and in vitro infectiousness [68–70]. The mutations S:A222V and S:S477N are among those with strongest topological signal for recurrence (see Figure 5). Both are associated with lineage B.1.177 / 20E (EU1) which emerged in Europe in mid-2020, and S:S477N is now also seen in the Omicron variant B.1.1.529 [71]. While S:S477N is known to affect the binding affinity to the ACE2 receptor [3, 72] leading to a slight increase in fitness, there is no conclusive evidence yet whether or not S:A222V also results in higher transmissibility [73]. Our time series analysis for the latter shows that the particularly strong tRI signal is still rising after the initial surge in prevalence in the European summer of 2020 (see Figure 5). This suggests that S:A222V is notably recurrent, independently of its impact on viral fitness.
Each panel shows a comparative time series analysis (topological recurrence index (tRI) vs. prevalence) with monthly timesteps ranging from December 2019 until February 2021, in comparison with an ancestral state reconstruction obtained with standard phylogenetic methods. The tRI time series maps the dynamics of topological signals for specific mutations and therefore allows for the retrospective assessment of adaptive processes over a longer period of time. (A) Topological footprint of a highly recurrent mutation in the example of the mutation ORF1ab:G11083T. A time series analysis shows a monthly increase of the tRI, while the prevalence stays low. This indicates that the mutation has been reemerging frequently and steadily since the beginning of the pandemic. The consistently low prevalence suggests that the mutation is neutral or deleterious, as a beneficial substitution would be expected to establish itself in larger subpopulations. (B, C) Recurrence of the mutations S:A222V and S:S477N persists after its initial surge in prevalence in mid-2020. (D) The mutation S:D614G shows a pattern typical for an adaptive mutation—after a rapid increase in tRI and prevalence the tRI reaches a plateau once the mutation has become dominant, superseding the wild type in the early phase of the pandemic. All ancestral state reconstructions are based on the Nextstrain tree [19] of a curated subsample of 3,507 sequences from the GISAID dataset as of 3 March 2021 [18, 36] (see Methods). The respective amino acid substitutions are highlighted in purple.
We noticed that in particular on the receptor-binding domain (RBD), several of the topologically significant amino acid changes (tRI 8) are found in Variants of Interest (VOI) or Variants of Concern (VOC) [71], with a distinct accumulation in the receptor-binding motif (see Table 1, Supplementary Information Figure 2 and Supplementary Information Figure 3). Specifically, the substitutions S:N501Y, S:E484K, S:L452R, S:K417N, S:F490S and S:S494P all result in reduced binding of polyclonal convalescent plasma [4] and exhibit a distinct increase in tRI and prevalence starting in late 2020 (see Supplementary Information Figure 3). This pattern is likely due to selective pressures induced by immune evasion in a host population with rising immunity. While both S:N501Y and S:N501T produce comparable tRI signals and induce a slight antibody escape, the fact that S:N501Y has a comparatively high prevalence of 19%, and is seen in several VOCs, is probably due to the additional increase in ACE2-binding affinity caused by the asparagine-to-tyrosine substitution. Topological signals for the mutations S:Y453F and S:F486L exclusively originated from a small subpopulation in minks [74–76]. The fact that tRI signals remain low (tRI = 8) despite becoming significant already in June/May 2020 suggests that both mutations have an adaptive effect in minks but do not confer a significant fitness advantage in humans.
The table displays the top ten amino acid substitutions on the receptor-binding domain with significant topological recurrence index (tRI 8), plus a few more selected mutations, together with the tendency of the tRI and the date of initial acquisition of a significant tRI signal. Several of the mutations with strong tRI signal occur in lineages designated as VOI/VOC [71]. For the full table see Supplementary Information Figure 2.
Correlation with host adaptation
Our analysis of the first year of the pandemic revealed that on the receptor-binding domain, there is a strong correlation between significant topological signals (tRI 8) and an increase in ACE2-binding affinity [3] compared to the wild type (Fisher’s exact test, p < 0.01; Figure 6). We did not find a similar correlation between significant tRI and an increase in plasma antibody escape [4], which is plausible as immune evasion has become a relevant factor for the evolution of the virus only towards the end of 2020. However, among those mutations with strong topological signal of convergence (tRI > 12) we found a correlation between increased mean plasma antibody escape > 0.025 and the initial acquisition of a significant tRI signal after October 2020 (Fisher’s exact test, p < 0.05; Figure 6). This indicates that beginning in late 2020, immune escape became an increasingly important evolutionary factor.
Both diagrams visualize contingency tables underlying Fisher’s exact test. (A) There is a strong correlation between a significant topological recurrence index (tRI 8) and an increase in ACE2-binding affinity in the year from December 2019 until February 2021. (B) Among highly topologically reccurent amino acid changes (tRI > 12), more recent ones show increased mean plasma antibody escape. There is a significant correlation between increased mean plasma antibody escape > 0.025 and the initial acquisition of a significant tRI signal after October 2020 (Fisher’s exact test, p < 0.05). Experimental data taken from Starr et al. [3, Table S2] and Greaney et al. [4, Table S3].
The latter observation led us to investigate to what extent topological signals of recurrence for amino acid changes on the Spike gene are modulated by the dynamic SARS-CoV-2 fitness landscape [77] shaped by increasing host immunity and changes in viral genetic background through the emergence of new dominant virus variants. To that end, we performed a topological recurrence analysis of GISAID data exclusively from the United Kingdom covering the whole pandemic (December 2019 until March 2023; see Figure 7 and Methods). We recorded a significant increase in amino acid variations newly flagged by significant tRI during the second half of the year 2020, with a distinct accumulation in the S1-S2 cleavage region, which is a determinant of SARS-CoV-2 infectiousness and transmissibility [61, 66, 78], and to a lesser extent in the N-terminal domain and the receptor binding domain, which contain epitopes for antibody binding [61–65]. This is compatible with the observation that host adaptation during the first year of the pandemic was mainly enhancing infectiousness and transmissibility in a largely immunologically naive host population, as manifested in the Alpha variant, while immune escape became an evolutionary factor only towards the end of the year 2020 when infection-induced population immunity was on the rise [77]. After a longer period with only few tRI signals we observed another significant increase in amino acid variations newly flagged by significant tRI beginning in the second half of the year 2021, with a distinct accumulation in the N-terminal domain and the receptor binding domain in the S1 subunit as well as in the fusion peptide and the heptapeptide repeat sequence 2 in the S2 subunit, which are targets for the action of neutralizing antibodies [62–65]. This confirms that during this second phase of the pandemic, immune evasion became the dominant factor of adaptation in the presence of widespread population immunity, as manifested in the Delta and Omicron variants, after the great majority of the UK population had received a second dose of vaccination by the third quarter of 2021 [77].
Dynamics of the interaction between topological signals of recurrence and host immunity/waves of dominant variants in the United Kingdom. For each month from January 2020 until March 2023, the heatmap shows the density of Spike gene amino acid variations for which the topological signal (tRI) turned significant during that month (density is centered moving average with window size of 101 sites). Topological recurrence analysis based on GISAID data from the UK covering the whole pandemic from December 2019 until March 2023 (see Methods); additional epidemiological information taken from Carabelli et al. [77, Figure 4].
(A) Specificity and sensitivity of tRI-based vs. prevalence-based identification of variable Spike amino acid sites likely to become highly prevalent (prevalence contained in 10% percentile of all mutations) three months thereafter. (A1) Reliability of identifications for each of the months October 2020 through February 2021. While false-positive rates of tRI-based and prevalence-based analyses are similar and remain small, the false-negative rate of tRI-based analyses is better than that of prevalence-based analyses (29% vs. 50%). (A2) Contingency table for tRI-based analyses for February 2021. Variable Spike amino acid sites were flagged if they exhibited a significant and increasing tRI signal from August until November 2020, and designated prevalent if their prevalence was contained in the 10% percentile of all mutations in February 2021. (B) Typical comparative time series pattern (tRI vs. prevalence) of an emerging adaptive mutation as observed in the escape mutation S:E484K. The tRI rises to a significant level (tRI = 5) in September 2020, while the prevalence stays very low < 0.02% and shows a visible increase to 0.06% by a factor of 3 for the first time several weeks later in October 2020. Both tRI and prevalence show a rising tendency in February 2021. (C) Cycle localization in the genome alignment yields geographic and temporal information about the independent acquisition of S:E484K. The diagram shows SNV cycles created by reticulate events involving the Zeta variant (P.2) in Brazil in October 2020, and within the Alpha variant (B.1.1.7) in England in January 2021. (D) Typical comparative time series pattern observed in emerging adaptive Spike mutations, with the onset of significant tRI signals preceding the rise in prevalence by several weeks.
We observed that many, though not all, of the mutations that are of high prevalence also tend to show a strong topological signal (see Figure 4, Supplementary Information Table 1 and Supplementary Information Table 2). While from a conceptual point of view tRI and prevalence are independent quantifiers of viral evolution, this phenomenon can largely be explained by the observation that adaptiveness is a confounding factor for both of them: evolutionary advantages lead to an increase in both tRI and prevalence. One of the advantages of the topological recurrence index, however, lies in its great sensitivity and specificity towards signals associated with convergent evolution. While the prevalence depends heavily on the dynamics of the pandemic, including founder effects, the topological signal does not.
Surveillance of emerging potentially adaptive mutations
The topological recurrence analysis has the capability to identify specific mutations at an early stage that are likely to confer a fitness advantage to the virus and may become prevalent in the future. We observed that the sensitivity and specificity of the topological recurrence index (tRI) is sufficient to detect signals of convergent evolution already at very low mutation frequencies, often several months before the respective mutation becomes recognizable by its prevalence in the population. As an initial assessment of the reliability of this method, we used topological signals to identify mutations at a given point in time that are likely to become highly prevalent three months later. A variable amino acid site was flagged if it exhibited a significant and increasing tRI signal that persisted from August until November 2020. We found that flagged mutations were strongly correlated with a high prevalence contained in the 10% quantile of all mutations in February 2021 (Fisher’s exact test, p < 10e 10; Figure 8). Results of similar quality were obtained for mutations flagged in any of the months October 2020 through January 2021. In contrast, a similar analysis based on prevalence (replacing tRI) had about the same false-positive rate, but a significantly higher false-negative rate, especially at low mutation frequencies (see Figure 8). Further validation comes from the observation that the onset of a significant tRI signal precedes the actual rise in prevalence by several weeks for the fitness increasing mutations S:L18F [79], K417N [4], L452R [4, 80], S:E484K [4, 80], S:N501Y [3, 4] and S:P681H [81, 82] seen in the VOCs Alpha, Beta, Gamma, Delta and Omicron [71] (see Figure 8). Moreover, looking at the four VOCs Alpha, Beta, Gamma and Delta which occurred during the first year of the pandemic, among defining Spike variable amino acid sites we found a substantial number of sites expressing a significant topological signal of recurrence already months before the variant was first reported (see Supplementary Information Figure 5). In practice, one may wish to improve the tRI-based identification of emerging potentially adaptive mutations by invoking other indicators of positive selection, such as the prevalence and its tendency, or experimental data. Focusing on the RBD, we found that the amino acid changes S:A348S, S:N354D, S:P384L, S:N439K, S:G446V, S:A475V, S:E484G, S:F490S, S:N501T and S:Y508H developed a significant tRI signal over the course of the first year of the pandemic, showed a rising prevalence in many cases, and were all associated with an increased mean plasma antibody escape > 0.01 [4, 83], but had low prevalence < 5% and had not been seen in any VOI/VOC as of February 2021 (see Table 1, Supplementary Information Figure 2 and Supplementary Information Figure 4). Mutations at these RBD residues are likely to confer a fitness advantage to the virus and might therefore appear in future variants. In fact, several months after our analysis was completed, the immune escape mutations S:F490S [4, 80] and S:G446S [4, 84], which had not been observed in notable variants before, occurred in the Lambda [85, 86] and Omicron [87] variants, which were designated as VOI/VOC in June/November 2021 [71]. By investigating explicit representatives of topological cycles in the genomic dataset, we were able to extract geographic and temporal information about the independent acquisition of topologically recurrent mutations during the pandemic (see Figure 8).
Ongoing convergent evolution and the Omicron variant
We analyzed topological signals of adaptation also at later phases of the pandemic, by applying our CoVtRec surveillance pipeline to two curated GISAID Spike gene alignments [18, 36] covering the period from December 2019 until January 2022 and March 2023, respectively (see Methods). As of January 2022, a total of 312 non-synonymous and 192 synonymous Spike mutations with significant tRI signal were found (see Figure 9, Supplementary Information Table 4 and Methods). Here signals with tRI 76 are statistically significant (p < 0.05; see Methods). We observed that on the receptor-binding domain, amino acid changes with significant topological signal accumulate in the receptor-binding motif (see Figure 9).
(A) The table displays selected amino acid substitutions on the Spike gene with significant topological recurrence index (tRI 76). For the full table see Supplementary Information Table 4. (B) Heatmap of all amino acid variations across the RBD showing topological signals of convergence. Significant signals (tRI 76; shown in red) accumulate in the receptor-binding motif. (C) Topological signals of convergence for the defining Spike mutations in the Omicron variant (B.1.1.529; BA.1, BA.2 and BA.3 [88, 90, 91]). 38% of these mutations are significantly topologically recurrent with tRI 76 (highlighted in boldface) and are likely due to convergent evolution. Corresponding tRI time series analysis charts are provided in Supplementary Information Figure 13.
Of particular interest are the defining amino acid substitutions in the Omicron variant B.1.1.529 [87, 88], which was designated as VOC in November 2021 [71] and was rapidly spreading all over the world in January 2022 [89]. Our analysis revealed that 38% (18 out of 48) of the Spike amino acid substitutions seen in B.1.1.529 (BA.1, BA.2 and BA.3 [90, 91]) —with 50% (15 out of 30) in the sublineage BA.1 and 50% (14 out of 28) in the sublineage BA.2— show significant topological signal of recurrence (see Figure 9). This includes, among others, the following Spike amino acid substitutions: S:G339D, S:N440K, S:S477N, S:T478K and S:N501Y, which are known to enhance binding to the human ACE2 receptor; S:K417N, S:G446S and S:E484A, which may reduce polyclonal antibody binding; S:P681H at the S1-S2 cleavage region [66, 78, 92]. Our findings provide further insight into the possible origins of the Spike mutations observed in the Omicron variant [93]. In fact, all 18 substitutions with significant topological recurrence index (tRI 76) likely arose from convergent evolution, while the emergence of substitutions showing weak topological signal may be due to other reasons (see Figure 9).
Finally, a topological recurrence analysis covering the whole pandemic from December 2019 through March 2023 revealed a total of 328 non-synonymous and 212 synonymous Spike mutations with significant tRI signal (see Supplementary Information Table 5 and Methods). Here signals with tRI 75 are statistically significant (p < 0.05; see Methods). These mutations are presently topologically recurrent. They potentially confer a fitness advantage to the virus and could appear in future variants.
Discussion
In the current COVID-19 pandemic, the early identification of emerging adaptive mutations in large SARS-CoV-2 genomic datasets is of paramount importance. Such mutations could be associated with vaccine resistance or higher transmissibility, among other concerning attributes [94]. We present here an effective and unbiased method, based on a technique from Topological Data Analysis known as persistent homology, that can rapidly identify the presence of these mutations and is able to efficiently deal with the ever-increasing wealth of genomic sequencing data created by global public health surveillance. While phylogeny-based approaches rely on the construction of a phylogenetic tree, often in combination with phylogenetic model assumptions, our method quantifies the potential adaptiveness of mutations merely by their topological footprint in the genome alignment and is agnostic about biological reasons. It enables the longitudinal analysis of the dynamics of topological signals of adaptation at daily resolution over any given period of time. This makes it a useful tool both for the real-time monitoring and for the retrospective characterization of the convergent evolution of the coronavirus starting from the very beginning of the pandemic.
In benchmark tests, our method performed faster than state-of-the-art phylogeny-based tools such as IQTree, UShER and HomoplasyFinder [9, 29, 35]. In particular, it accomplishes the analysis of current alignments with millions of SARS-CoV-2 genomes within just one day. Orthogonal biological validation for our method comes from a strong correlation between mutations flagged as adaptive by our Topological Data Analysis approach and mutations flagged as fitness increasing in a recent phylogeny-based assessment of fitness effects of mutations by Bloom & Neher [7].
Drawing from data from GISAID [18, 36], we applied our CoVtRec analysis pipeline to characterize the evolution of the coronavirus in three different phases of the pandemic. The analysis of the first year ranging from December 2019 until February 2021 revealed a total of 700 (518) topologically recurrent mutations distributed across the whole genome (Spike gene), which shows that convergence plays a significant role in the evolution of SARS-CoV-2. This is compatible with results from phylogeny-based analyses of the adaptive evolution during the first year of the pandemic [13, 95]. We found that our method reliably identifies highly homoplasic sites across the whole genome, making it a useful tool in the design of unbiased masking schemes in phylogenetic inference.
We studied the correlation between topological signals of recurrence and host adaptation, including experimental [3, 4] and epidemiological data [77] into our topological recurrence analysis. For amino acid changes on the receptor-binding domain we found a correlation between significant topological signals and an increase in ACE2-binding affinity compared to the wild type of the virus. Moreover, our analysis suggests that topological signals of recurrence are modulated by the SARS-CoV-2 fitness landscape. We inferred that during the first year of the pandemic, adaptation was mainly enhancing infectiousness and transmissibility in a largely immunologically naive host population, while in the second phase of the pandemic beginning in the year 2021, immune evasion became an important factor of adaptation in the presence of increasing population immunity.
As we observed, our method likely also captures topological footprints of APOBEC-mediated host RNA editing processes which can give rise to recurrent hypermutable C>T transitions. These mutations may show significant topological recurrence index (tRI) in our analysis but could nevertheless not be fitness increasing [52, 53]. We found that this in particular applies to synonymous mutations, which make up a substantial proportion of single nucleotide variations with significant topological recurrence index. However, it is likely that other processes, such as adaptation to the human codon usage, could also drive topological signals of synonymous mutations. Caution is therefore required in the biological interpretation of significant topological signals of recurrence.
We demonstrated that our method is sufficiently sensitive to detect emerging adaptive mutations already at an early stage when mutation frequencies are still very low, and provides geographic and temporal information about virus isolates involved in concrete reticulate events. We identified mutations on the Spike gene that showed significant topological signals of recurrence with a rising tendency in February 2021. These mutations had the potential to confer a fitness advantage to the virus in a largely immunologically naive host population. In fact, several of these candidate mutations have later received special attention as features of the Lambda and Omicron Variants of Interest/Concern [71], not until several months after our analysis was completed.
We explained how our results can aid the understanding of possible reasons for the occurrence of certain mutations in Variants of Concern such as the Omicron variant. We found that as of January 2022, 38% of the Spike gene mutations observed in the Omicron variants are likely due to convergent evolution.
The topological analysis of the current phase of the pandemic in March 2023 revealed a total of 540 topologically recurrent mutations on the Spike gene. Updated CoVtRec analyses are available at https://tdalife.github.io/covtrec.
Epistasis has been reported to play a key role in SARS-CoV-2 evolution [5, 96–101], but the phylogeny-based computational characterization of this phenomenon across all sites along the whole genome remains challenging due to the rapidly growing combinatorial complexity in possible interactions [102]. While in the present study we confined ourselves to the assessment of the adaptive potential of single mutations in the presence of specific genetic backgrounds, our method opens new ways for the computational mapping of the epistatic interaction in the evolution of the coronavirus. For example, high-resolution time series analyses can help to identify networks of interacting mutations by comparing the onsets of significant topological signals of adaptation, possibly after restricting to fixed genetic backgrounds. In addition, topological footprints of networks of interacting mutations should be accessible by means of evaluating topological cycles in which adjacent genomes are allowed to differ by more than just one nucleotide substitution.
Our method has several limitations. First, in contrast to phylogeny-based approaches, our approach based on topological footprints will resolve homoplasies only on small genetic distance scales and will not distinguish between forward and backward mutations. Also, in its present form the topological analysis does not include insertions and deletions. Second, sequencing errors [103, 104] may tamper topological signals to some extent. We implemented several filtration rules in order to control this effect. In this regard it is also important to note that persistent homology is robust with respect to noise in the genomic dataset [105]. Third, the topological recurrence analysis is agnostic about biological reasons, and mutations flagged by significant topological recurrence index (tRI) need not always be adaptive. The downstream analysis of mutations with significant tRI therefore requires care and should invoke additional biological information, such as masking schemes for problematic sites, epidemiological data or experimental evidence of adaptation available from the literature.
Based on these insights, we propose persistent homology as a surveillance system for emerging potentially adaptive mutations in the ongoing evolution of the coronavirus SARS-COV-2 as well as for the retrospective characterization of the dynamics of adaptation for specific mutations. We foresee this capability also in future pandemics of other pathogens. Our method allows for a targeted mapping of topologically recurrent mutations in any region of the genome, compiling a shortlist of flagged mutations which can guide and motivate the experimental study of adaptive effects of specific mutations, also in genes other than the Spike gene. We envision a combined effort between public health organizations with extensive sequencing of viral genomes, the computational characterization of potentially adaptive variants, and the experimental phenotypic characterization of these variants.
Methods
Data acquisition and data preparation
We used three separate datasets of SARS-CoV-2 genome sequences provided by the GISAID EpiCoV Database [18, 36]. The first dataset was obtained by downloading all available SARS-CoV-2 whole genome sequences as of 28 February 2021, isolated from human and animal hosts, that carried the following attributes: “complete”, “high coverage”, “low coverage excluded”, “collection date complete”. This dataset comprised 450,473 sequences. We aligned all sequences with MUSCLE (v3.8.31) [106], using as reference genome the sequence Wuhan/Hu-1 with accession number EPI ISL 402125, truncated at the start codon of ORF1ab (reference site position 266) and the stop codon of ORF10 (reference site position 29,674). Subsequently all sequences containing at least one ambiguous site (letter “N”) were removed. This resulted in an alignment comprising 303,651 complete SARS-CoV-2 genomes of length 30,130nt. The second dataset is the GISAID alignment msa 0126.fasta from January 2022. It comprises 6,703,606 SARS-CoV-2 whole genome sequences that have been aligned to the reference sequence Wuhan/WIV04 with GISAID accession number EPI ISL 402124 using MAFFT (v7.497) [107]. Sequences were truncated to the Spike gene (reference site positions 21,563 to 25,384), and finally sequences containing any letters other than A, C, T or G were removed. This resulted in an alignment comprising 318,851 distinct complete SARS-CoV-2 Spike genes of length 4,796nt. The third dataset is the GISAID alignment msa 0324.fasta from March 2023. It comprises 13,766,674 SARS-CoV-2 whole genome sequences that have been aligned to the reference sequence Wuhan/WIV04 with GISAID accession number EPI ISL 402124 using MAFFT (v7.497) [107]. Sequences were truncated to the Spike gene (reference site positions 21,563 to 25,384), and finally sequences containing any letters other than A, C, T or G were removed. Subsequently two sub-alignments were selected: one with sequences exclusively from the United Kingdom (England, Wales, Scotland and Northern Ireland), resulting in an alignment comprising 91,063 distinct complete SARS-CoV-2 Spike genes of length 7,942nt; and a second one with sequences from all over the world that occur in at least two identical copies, resulting in an alignment comprising 220,978 distinct complete SARS-CoV-2 Spike genes of length 7,942nt. A list of accession numbers of all GISAID sequences used in this work, along with an acknowledgement of the contributions of both the submitting and the originating laboratories, can be retrieved from https://doi.org/10.55876/gis8.230731by.
The experimental data on viral phenotypes by Starr et al. and Greaney et al. was downloaded from [3, Table S2] and [4, Table S3]. The values for the mean plasma antibody escape used in Supplementary Information Figure 2 were computed by averaging mutation escape values for every substitution in [4, Table S3] over all subjects.
The data on fitness effects of Spike gene amino acid changes as of March 2023 by Bloom & Neher was downloaded from [108, Spreadsheet aa fitness.csv]. Subsequently the columns “gene”, “aa site”, “aa” and “fitness” were extracted.
The data on filtration rules for highly homoplasic sites, problematic sites and potential artifacts in SARS-CoV-2 sequence alignments by Turakhia et al. and De Maio et al. was taken from [47] and downloaded from [46, Spreadsheet problematic sites sarsCov2.vcf]. Subsequently the columns “POS”, “FILTER” and “INFO” were extracted.
The data on COSMIC mutational signatures for Genome GRCh37 by COSMIC was downloaded from [109]. Subsequently the column “SBS2” was extracted.
Distance matrices
We used Hammingdist (v0.19.0) [41] to compute the genetic distance matrix of a given alignment of genome sequences. For any pair of sequences in the alignment, this matrix gives the Hamming distance between the two sequences, which is the number of site positions at which the nucleotides in the two aligned sequences differ. Noteworthy, our convention in this work is that insertions and deletions (dashes “-” in aligned sequences) do not contribute to the genetic distance.
From the whole genome alignment covering the first year of the pandemic we created 15 time buckets, each ranging from December 2019 to one of the months between December 2019 and February 2021. For each time bucket, a time bucket sub-alignment of all genetically distinct sequences whose collection dates belong to the given time bucket was created by selecting isolates by their date stamp and removing genetically identical sequences (Hamming distance = 0). The largest time bucket sub-alignment ranging from December 2019 to February 2021 contained 161,024 genetically distinct sequences. Then for each such sub-alignment the corresponding genetic distance matrix, which is a sub-matrix of the distance matrix of the whole alignment, was derived. We obtained 15 distance matrices of whole genome time bucket sub-alignments. This process was repeated for all sub-alignments after truncating sequences to the Spike gene (reference site positions 21,563 to 25,384). We obtained 15 distance matrices of Spike gene time bucket sub-alignments.
The genetic distance matrices for the other two Spike gene alignments (as of January 2022 and March 2023) were computed in a similar fashion. Subsequently a Vietoris-Rips transformation as implemented in the MuRiT algorithm was applied [110].
Topological Data Analysis and viral evolution
Topological Data Analysis (TDA) is a field of data science that aims to study the shape of large datasets, by extracting topological structures and patterns. Such topological structures have associated dimensions: structures of dimension zero can be thought of as the connected components, and structures of dimension one are essentially the loops, or topological cycles, of the dataset. Structures of higher dimensions can also be defined, but are also more difficult to interpret. Here we are interested in reticulate evolutionary processes, thus we choose to focus on topological structures in dimension one, since topological cycles can be interpreted as signals of divergence from phylogenetic trees (see Figure 2).
Datasets often come as point clouds: in our setting, each point corresponds to a virus genome sample, and lies in a high-dimensional space where each nucleotide of the genome is a dimension. A common way to extract the phylogenetic network from this point cloud simply amounts to connecting samples as soon as their genetic distance is less than a given threshold r > 0. This results in a (neighborhood) graph, whose set of cycles provides candidates for the topological structures in dimension one of the true underlying network. However, a main limitation of this approach comes from the fact that relevant topological structure typically appears at multiple scales (see Supplementary Information Figure 8).
The most common way to handle this issue in Topological Data Analysis is to actually compute and track the cycles for all possible values of r ranging from 0 to +. As r increases, some cycles can appear, and some already existing cycles can disappear, or get filled in. The whole point of Topological Data Analysis is to record, for each cycle, its radius of appearance, or birth time, and radius of disappearance, or death time. This is called the persistent homology of the point cloud. The construction, based on a scale parameter r, can be summarized as follows. The input is a distance matrix describing the dataset, considered as a finite metric space. First, consider the geometric graph at scale r, whose vertices are the data points, with any two points connected by an edge whenever their distance is at most r. Generalizing this construction, the Vietoris–Rips complex at scale r connects any subset of the data points with a simplex (an edge, a triangle, a tetrahedron, or a higher-dimensional generalization thereof) whenever all pairwise distances in the subset are at most r. A Vietoris–Rips complex is a particular type of simplicial complex, a higher-dimensional generalization of graphs which is of crucial interest in algebraic topology, in particular in homology theory. The family of Vietoris–Rips complexes for all parameters r is called the Vietoris–Rips filtration. It provides a multiscale method to extract cycles of various sizes, and to encode them in a so-called persistence barcode: each bar, or interval, in this barcode corresponds to a cycle representing a topological feature (a reticulate evolutionary process in our case), and the bar endpoints correspond to its radii of birth and death (the maximum genetic distance between consecutive samples forming the cycle, and, roughly, the maximum pairwise genetic distance between samples forming the cycle).
Each bar indicates the presence of a reticulate event, implying that the evolutionary history cannot be fully explained in terms of a single phylogenetic tree. The mathematical background of this phenomenon is a classical theorem due to Rips, which asserts that trees have trivial persistent Vietoris–Rips homology [111]. The corresponding cycle in the associated reticulate phylogeny can then be localized in the sequence alignment by tracing it back to the isolates that constitute the reticulate event. Moreover, the length of the bar represents the cycle size. In our case, this corresponds to the length of the reticulate evolutionary process, which allows to distinguish, for instance, between homoplasies and recombinations (see Figure 2).
Computation of persistent homology
Ripser is a state-of-the-art software for the computation of persistent homology based on the topological construction of Vietoris–Rips complexes, developed by one of the authors [40]. For the computation of the persistence barcode, Ripser resorts to various optimizations, which are crucial when handling datasets of the size considered here. Notably, Ripser computes persistent cohomology, which is not based on cycles but instead on cocycles, often described intuitively as cuts that tear open a hole. In order to obtain the requisite cycles representing the features in persistent homology, we used a custom version of Ripser that subsequently carries out a second computation, this time based on cycles instead of cocycles. While a naive computation based on homology would be prohibitively expensive, the previous computation of the persistence barcode based on cocycles makes the subsequent computation of representative cycles feasible.
For our computations, we used a customized version of Ripser to compute the representative cycles for the persistent homology of the Vietoris–Rips filtration associated to the genetic distance matrix for each time bucket sub-alignment (whole genome and Spike gene). As we are only interested in SNV cycles, the computation of persistence barcodes for the time bucket sub-alignments was restricted to small genetic distance scales (Ripser scale parameter threshold set to 2), which greatly increases the speed of the computation.
The homological features identified by persistent homology admit different representative cycles. In order to obtain cycles that fit tightly to the data points, our customized version of Ripser uses a method called exhaustive reduction [112, 113], which can be roughly summarized as follows. Whenever a representative cycle contains an edge that also appears in another cycle as the longest edge, a tighter representative can be obtained by replacing the edge with the remaining edges from the other cycles, which all have shorter length.
Topological features are statistically significant
We estimated the expected number of topological cycles in persistent homology that are created by random homoplasic events in the GISAID dataset covering the first year of the pandemic with 161,024 genetically distinct whole genome sequences collected from December 2019 until February 2021. To this end, we simulated several evolutionary scenarios under the following assumptions: uniform probability distribution for substitutions across the genome, no variations in fitness, and zero recombination rate.
We generated forward simulations of viral evolution based on a Wright-Fisher model using SANTA-SIM (v1.0) [114] with fixed parameters: number of generations (N = 10, 000), number of sequences sampled from the population per time step (n = 15), recombination rate (ρ = 0), and variable parameters: mutation rate per site per generation, effective initial population, carrying capacity, population growth rate per generation. We considered five scenarios: In scenarios I-III we varied the mutation rate under the assumption of fixed population size, while in scenarios IV and V we investigated the effects of logistic growth of the viral population (see Supplementary Information Figure 9). The range of mutation rates in scenarios I-III were chosen such that the diversity in the simulated phylogenies are in close correspondence to the observed diversity in the GISAID dataset. While a mutation rate of µ = 0.75E 7 substitutions per generation per site systematically underestimates the maximal distances to the root, the highest value of µ = 1.25E 7 produces slightly larger maximal values. In fact, scenario II with µ = 1.00E 7 reproduced the observed maximal distance accurately and provides a good approximation of the GISAID dataset (see Supplementary Information Figure 9). Major differences between simulations and the GISAID dataset are likely due to epidemiological phenomena in the ongoing pandemic such as variable population size and sequencing rate, and the spread of certain variants.
For each of the simulated datasets, we computed its persistent homology in dimension one using Hammingdist [41] and Ripser [40]. In order to keep overall computational expenses at a reasonable level, we resorted to extrapolations from smaller simulated datasets to the size of the GISAID dataset with 161,024 sequences. For all scenarios we produced 100 simulations for each of the following values of the effective population p (resp. carrying capacity c): 100, 500, 1000, 2500, 5000, 7500, 104, 105. Additionally, we included five simulations for p = 106 to achieve a better support of the extrapolation fit. For each value of p we randomly chose 60% of the simulations as training data, used to determine the parameters of different models in a non-linear least squares fit, while the remaining 40% were reserved for later validation and comparison of the models.
For each scenario we considered a quadratic, cubic, powerlaw and exponential model for the observed points (xi, yi), and linear and powerlaw fits for the squared residuals (yi yfit(xi))2 in the training data (see Supplementary Information Figure 11). In each model, we then used the resulting fits mean(x) and var(x) as estimators for the mean and variance of an underlying Panjer (a, b, 0)-class distribution [115, 116]. The quantiles of the observed number of cycles in the training data fit the quantiles of the Panjer distribution with corresponding mean and variance remarkably well (see Supplementary Information Figure 10). We then determined the likelihood L = Πi PPanjer(y = yi | mean(xi), var(xi)) to observe the validation data {(xi, yi)}. For each model, the corresponding log-likelihoods are listed alongside the corresponding fits in Supplementary Information Figure 11. According to the log-likelihoods, the variance of the Panjer distribution is generally best described by a powerlaw behaviour. An exception is scenario II, for which the small sample of 5 simulations at p = 106 has an uncharacteristically small variance that skews the fits and corresponding likelihoods. Among the models that assume a powerlaw dependence of the variance, again with exception of scenario II, the cubic-powerlaw model yields maximum likelihoods. Finally, we determined the 95% prediction intervals for the expected numbers of random cycles by use of the cubic-powerlaw extrapolation of mean and variance of a Panjer distribution (see Supplementary Information Figure 9). The validation data of scenarios I, IV and V, which were all based on the same mutation rate, are well described by the prediction intervals of scenario I. The prediction intervals of scenario V differ significantly from the other two scenarios only at high numbers of distinct sequences. This difference arises because simulations in scenario V generally produce fewer distinct sequences than scenario I and IV, such that a steeper extrapolation is not sufficiently penalized. Hence, the prediction intervals of scenario V illustrate the error margins of the extrapolations, but are not likely to faithfully represent the expected number of one-dimensional cycles. We also observe that higher mutation rates in scenarios II and III lead to smaller numbers of one-dimensional cycles in the dataset.
In conclusion, the 95% prediction interval of scenario V yields an upper bound between 1023 and 1171 expected random cycles in a dataset comparable to the GISAID dataset with 161,024 distinct sequences. Moreover, since the diversity of the GISAID dataset is better approximated by scenario II than by scenarios I, IV or V, it is reasonable to rely on the prediction interval of scenario II, which predicts that in 95% of the cases we expect between 362 and 408 random cycles (see Supplementary Information Figure 9).
Topological recurrence analysis
We performed a topological recurrence analysis both for the whole genome and for the Spike gene. In principle, this analysis can be done for any region of the genome, by working with intervals in the sequence alignment consisting of appropriately truncated genomes [38, 39]. As it turns out, the number of topological cycles, and hence the resulting tRI scores, will heavily depend on the choice of these intervals. Generally speaking, as evolutionary processes outside the specific genome region are ignored, genetic distances between isolates will decrease. This typically leads to the formation of more topological cycles at small genetic distance scales, and hence to stronger tRI signals. There is a trade-off between the strength of tRI signals and the genomic scope of the analysis. This explains why in our analysis of the first year of the pandemic, tRI signals in the Spike gene analysis are normally stronger than corresponding signals in the whole genome analysis (see Figure 4). We remark that there is in general no relation between tRI scores for different genomic regions.
Regarding the alignments covering the first year of the pandemic, we proceeded as follows: For each time bucket sub-alignment (whole genome and Spike gene) a complete list of SNV cycles (topological cycles all of whose edges correspond to single nucleotide variations) in this alignment was generated from the corresponding Ripser output. For each edge in an SNV cycle the endpoints of the edge correspond to a pair of uniquely determined sequences in the alignment that differ in exactly one nucleotide site position and hence determine an SNV. Then for each such SNV, its topological recurrence index (tRI) is by definition the total number of all SNV cycles containing an edge that gives rise to the given SNV. We restricted our analysis to SNVs with the following two properties: (i) one of the two nucleotides involved in the SNV agrees with the nucleotide in the reference sequence EPI ISL 402125 at that site position, and (ii) the SNV is isolated in the sense that at the two preceding and following site positions the nucleotides are the same as in the reference sequence. These two conditions ensure that the corresponding SAAV is uniquely determined by the SNV. We used custom code implemented in Python to compute the tRI of each such SNV for every time bucket sub-alignment (whole genome and Spike gene). Moreover, for every whole genome time bucket sub-alignment the prevalence of every SNV was computed as the quotient of the number of all sequences carrying that SNV by the number of all sequences in that sub-alignment. Note that the sub-alignments entering into this computation consisted of genetically distinct sequences. Finally, for every SNV the measurements of both tRI (whole genome and Spike gene) and prevalence for all time buckets were combined into a time series analysis chart.
Even if all SNV cycles arose through random processes, it is expected that the resulting tRIs are distributed uniformly among all observed mutations. So the probability for a given mutation to have tRI k is given by a binomial distribution where the number of trials corresponds to the number of mutations in SNV cycles, and the probability for success is the inverse of the number of mutations that are realized in the dataset. From this we deduce that in the whole genome analysis a tRI 2 is highly significant (p < 0.01), while for the Spike gene analysis any signal with tRI 8 is significant (p < 0.05).
For the topological recurrence analysis of the Spike gene alignments as of January 2022 and March 2023, we used Vietoris-Rips transformations in multipersistent homology as implemented in the MuRiT algorithm [110] to compute tRI time series analysis charts at daily resolution from the natural stratification by time induced by collection dates of viral isolates. For each mutation, we computed the tRI growth rate as the 14 days moving average of tRI data. We remark that the tRI growth rate is tampered by computational artifacts at the time when the tRI first attains a positive value. Moreover, towards the end of the period covered by the sequence data the sampling rate is normally still low, which will result in reduced tRI growth rates.
The effect of sequencing errors on topological results
Persistent homology is robust with respect to noise in datasets [105]. Nevertheless, sequencing errors in the SARS-CoV-2 alignments [103, 104] might tamper results of the topological recurrence analysis. We implemented the following filtration rules in order to control this effect. During data preparation, sequences not labeled “high coverage” and “complete” by GISAID, or containing any letters other than A, C, T or G, were removed from the alignment. Moreover, in subsequent analyses we only used significant tRI signals, while non-significant signals were discarded. Here the tRI significance level was computed from a test statistics with the null hypothesis assuming that tRI scores are distributed randomly across all sites on the genome. The combination of these filtration rules keeps the effect of sequencing errors on tRI signals to a minimum.
Our tRI significance levels are corroborated by the following estimate of the expected amount of erroneous tRI signal per mutation. Here we assume that the probability of a sequencing error at a given nucleotide site is of the order of 10−3 for typical GISAID consensus genomes. The creation of a tRI signal for a given SNV requires the formation of a topological SNV cycle, a process that involves at least two edges associated with some SNV (see Figure 1). The probability for this to be due to sequencing errors is therefore of the order of 10−3 10−3 = 10−6. The GISAID alignments we used contain about 300,000 resp. 4,000,000 distinct genomes. Hence for a given mutation the expected errors in tRI are 10−6 300, 000 = 0.3 resp. 10−6 4, 000, 000 = 4, which is below the tRI significance levels we observed in our analyses.
In view of the large number of available sequences, a common method to contain sparse sequencing errors is to confine the analysis to subalignments of genomes that appear in multiple identical copies in the alignment. We observed that the qualitative structure of the tRI landscape along the Spike gene did not change under this additional filter and characteristic and significant signals persisted (see Supplementary Information Figure 6). On the other hand, the topological analysis is sensitive enough to detect cycles that arise from unique genomes in the alignment which may nevertheless be correctly sequenced and carry important information about the evolution of the virus. As these valuable signals might otherwise have been lost, we preferred not to use this additional filter in our analyses.
Performance analysis
We performed a basic runtime comparison between Topological Data Analysis (TDA)-based methods and standard phylogeny-based methods for a random sample of 5,000 SARS-CoV-2 genomes, and ten nested sub-alignments thereof, drawn from the GISAID alignment covering the first year of the pandemic. We used IQTree [35] (v2.1.3, with default settings and fast search option) and UShER [29] (v0.5.6, with default settings) to reconstruct phylogenetic trees. Here for each sub-alignment, UShER placed the 500 newly added samples using a mutation-annotated tree obtained from the previous sub-alignment. The subsequent homoplasy analysis was performed with HomoplasyFinder [9] (with default settings). For the TDA-analysis we used Hammingdist [41] (v0.15.0) to generate genetic distance matrices, and a custom version of Ripser [40] (with scale parameter threshold set to 2) for the subsequent computation of persistence barcodes and representative cycles. All computations were carried out on a machine with Intel Xeon E7-4850 v4 processors and 128 kernels. The resulting runtimes for each sample are shown in Figure 3.
The computation of the genetic distance matrix for the whole genome alignment covering the first year of the pandemic with 303,651 sequences was carried out with Hammingdist [41] (v0.13.0) on a machine with Intel Xeon Gold 6230R processors and 52 kernels. The runtime was 57 minutes and the memory usage was 36 gigabytes for the whole genome analysis (49 seconds and 2 gigabytes for the Spike gene analysis).
The computation of the persistence barcodes for all monthly sub-alignments of the corresponding alignment with 161,024 genetically distinct sequences was carried out with Ripser [40] (with scale parameter threshold set to 2) on an Intel Xeon Gold 6230R processor. The runtime and memory usage for each sub-alignment are shown in Figure 3. For the largest time bucket sub-alignment ranging from December 2019 to February 2021 with 161,024 genetically distinct sequences, the runtime was 45 minutes and the memory usage was 49 gigabytes for the whole genome analysis (59 seconds and 2.1 gigabytes for the Spike gene analysis).
For the GISAID alignment as of March 2023 with 13,766,674 SARS-CoV-2 Spike gene sequences the total computing time for the topological recurrence analysis was 12.7 hours.
Ancestral state reconstruction analysis
For the study of the evolutionary histories of topologically highly recurrent mutations (see Figure 5) we performed ancestral state reconstruction analyses using Mesquite (v3.61) [117]. As input we used a curated alignment of 3,507 genome sequences and its corresponding Maximum-Likelihood tree, downloaded from Nextstrain [19] on 3 March 2021. The tree was rooted using the oldest sequence available (EPI ISL 406798, collected on 26 December 2019). We inferred the evolution of each amino acid of interest along this SARS-CoV-2 tree using a parsimony approach.
Mutational patterns and signatures in topological signals
We performed a frequency analysis for nucleotide transitions with significant tRI signal across the whole genome during the first year of the pandemic from December 2019 until February 2021 (see Supplementary Information Figure 7). The analysis was carried out separately for non-synonymous mutations and synonymous mutations, based on tRI scores taken from Supplementary Information Table 1. We analyzed the local sequence context of these transitions by determining the signature of relative allele frequencies of trinucleotide motifs centered at a given single nucleotide variation. Here the relative allele frequency was computed as the quotient of the total number of trinucleotide motifs N[X>Y]N centered at a mutation X>Y with significant tRI by the total number of triplets NXN in the reference genome Wuhan/Hu-1 EPI ISL 402125. In this way, we obtained mutational signatures both for non-synonymous and synonymous mutations. We compared the resulting signatures with the COSMIC signatures [51] for Genome GRCh37 taken from [109] by restricting to trinucleotide motifs centered at any of the six substitutions C>A, C>G, C>T, T>A, T>C and T>G. Subsequently we computed the cosine similarity between our signatures and COSMIC signature SBS 2 for non-synonymous (cosine similarity of 0.39) and synonymous mutations (cosine similarity of 0.44).
Data availability
The SARS-CoV-2 genome data used in this work are available from the GISAID EpiCov Database [18, 36] at https://www.gisaid.org. To view the contributors of each individual sequence with details such as accession number, Virus name, Collection date, Originating Lab and Submitting Lab and the list of Authors, visit https://doi.org/10.55876/gis8.230731by. Experimental data on viral phenotypes by Starr et al. and Greaney et al. is available from [3, Table S2] and [4, Table S3]. Data on fitness effects of Spike gene amino acid changes by Bloom & Neher is available from [108]. Data on filtration rules for highly homoplasic sites, problematic sites and potential artifacts in SARS-CoV-2 sequence alignments by Turakhia et al. and De Maio et al. is available from [46, 47]. Data on COSMIC mutational signatures for Genome GRCh37 by COSMIC is available from [109].
Code availability
Code used for the analyses is available at https://github.com/ssciwr/hammingdist and https://github.com/Ripser/ripser/tree/tight-representative-cycles. All other code is available from the corresponding authors upon request.
Author contributions
M.B., L.H., M.N., J.P.G., R.R., A.O. designed the study; M.B., L.H., A.O. curated data; M.B., M.C., L.H., A.O., J.P.G. performed computational analyses; U.B., M.B., L.H., M.N., A.O. developed and implemented software; M.B., L.H., A.O. acquired computing resources; M.B., L.H., A.O. drafted the manuscript; all authors contributed to the final version of the paper.
Competing interests
R.R. is a founder of Genotwin, he is member of the Scientific Advisory Board of AimedBio and consults for Arquimea Research. The other authors declare no competing interests.
Additional information
Supplementary Information is available for this paper.
Correspondence and requests for materials should be addressed to M.B., L.H. or A.O.
Supplementary Information
Supplementary Information Table 1. External spreadsheet containing full results of the topological recurrence analysis for the whole genome as of February 2021. The table lists mutations together with their topological recurrence index (tRI) and prevalence. All mutations with statistically significant tRI 2 are included. The table is augmented by filtration rules (columns “FILTER” and “INFO”) taken from Turakhia et al. [33, 46] and De Maio et al. [47] for highly homoplasic sites that are normally removed or masked in SARS-CoV-2 phylogenetic analyses.
Supplementary Information Table 2. External spreadsheet containing full results of the topological recurrence analysis for the Spike gene as of February 2021. The table lists mutations together with their topological recurrence index (tRI) and prevalence. All mutations with tRI 2 are included, but only a tRI 8 is statistically significant. The table is augmented by filtration rules (columns “FILTER” and “INFO”) taken from Turakhia et al. [33, 46] and De Maio et al. [47] for highly homoplasic sites that are normally removed or masked in SARS-CoV-2 phylogenetic analyses.
Supplementary Information Table 3. External spreadsheet containing a sublist of the list in Supplementary Information Table 2 featuring all mutations on the receptor-binding domain together with their topological recurrence index (tRI) and prevalence.
Supplementary Information Table 4. External spreadsheet containing full results of the topological recurrence analysis for the Spike gene as of January 2022. The table lists mutations together with their topological recurrence index (tRI). All mutations with tRI 2 are included, but only a tRI 76 is statistically significant. The table is augmented by filtration rules (columns “FILTER” and “INFO”) taken from Turakhia et al. [33, 46] and De Maio et al. [47] for highly homoplasic sites that are normally removed or masked in SARS-CoV-2 phylogenetic analyses.
Supplementary Information Table 5. External spreadsheet containing full results of the topological recurrence analysis for the Spike gene as of March 2023. The table lists mutations together with their topological recurrence index (tRI). All mutations with tRI 2 are included, but only a tRI 75 is statistically significant. The table is augmented by filtration rules (columns “FILTER” and “INFO”) taken from Turakhia et al. [33, 46] and De Maio et al. [47] for highly homoplasic sites that are normally removed or masked in SARS-CoV-2 phylogenetic analyses.
Persistence barcode representing the persistent homology in dimension one of the GISAID dataset [18, 36] covering the first year of the pandemic from December 2019 until February 2021, comprising 161,024 genetically distinct high-quality SARS-CoV-2 whole genomes (see Methods). Each of the 2,899 bars in the barcode corresponds to a topological cycle in the reticulate phylogeny. The rich topology of the dataset indicates a multitude of reticulate events that shaped the evolution of the virus in the course of the pandemic. A total of 1,684 bars, which is 58% of all bars, concentrate at small genetic distance scales 2 and are therefore expected to be associated mainly with homoplasic events.
(A) Table of all amino acid changes with statistically significant topological recurrence index (tRI 8), plus two more selected mutations. The table provides the tRI together with its tendency, the initial acquisition date of a significant tRI signal (p < 0.05), the prevalence together with its tendency, notable Pango lineages containing the mutation [88], ACE2-binding affinity as in Starr et al. [3, Table S2], and mean plasma antibody escape as in Greaney et al. [4, Table S3] (see Methods). Mutations with rising tRI signal and mean plasma antibody escape > 0.01 (shaded) potentially confer a fitness advantage to the virus and are therefore candidates that might appear in future new variants. (B) Heatmap of all amino acid variations across the RBD showing any topological signal of convergence. There is a distinct accumulation of signals in the receptor-binding motif, while other regions on the RBD, notably residues 390-435, show only few signals of convergence.
Comparative time series analysis charts (tRI vs. prevalence) as of February 2021 for amino acid changes on the Spike gene seen in notable lineages that have been designated as VOIs/VOCs [71]. Significant tRI values are highlighted in boldface.
Comparative time series analysis charts (tRI vs. prevalence) for mutations on the receptor-binding domain that showed a significant tRI signal with rising tendency in February 2021 and are associated with an increased mean plasma antibody escape > 0.01 [4]. These mutations had low prevalence < 5% and were not seen in any VOI/VOC as of February 2021 (see Supplementary Information Figure 2), but mutations at corresponding residues are likely to confer a fitness advantage to the virus and might therefore appear in future variants. In fact, the immune escape mutations S:G446S [4, 84] and S:F490S [4, 80] later appeared in the Lambda [85, 86] and Omicron [87] variants, which were designated as VOI/VOC in June/November 2021 [71].
For each time bucket ranging from December 2019 until February 2021, the tables display the topological recurrence index (tRI) for all defining variable Spike gene amino acid sites in the given VOC. Significant tRI values are highlighted in boldface. For all four VOCs, a substantial number of defining variable amino acid sites showed a significant topological signal already months before the variant was first reported.
The panels display topological signals along the Spike gene as of February 2021 for different filtration parameters n. Sequence alignments were filtered by collecting only those sequences that appear in at least n identical copies in the alignment. (A) The tRI landscape for n = 1, based on the whole alignment with no additional filtration as it is used in the analyses in this study. (B) and (C) display the tRI landscape after additional filtration with parameters n = 5 and n = 10. The qualitative structure of the tRI landscape does not change: while the overall strength of tRI signals diminishes, characteristic and significant signals persist.
Frequency analysis of nucleotide substitutions across the whole genome showing significant topological signal during the first year of the pandemic (December 2019 until February 2021). The panels display total frequencies of single nucleotide variations for non-synonymous mutations (A1) and synonymous mutations (B1). There is an excess in C-to-T transitions which suggests that the topological recurrence index also captures the action of APOBEC-mediated host RNA editing processes. Relative frequencies of trinucleotide motifs centered at single nucleotide variations for non-synonymous variations (A2) and synonymous variations (B2) (see Methods). Statistical analyses based on Supplementary Information Table 1.
Each point represents a sample, and we display the geometric graphs (A), the resulting Vietoris–Rips complexes at different scales (B) and the persistence barcode in dimension one (C). If one only chooses one scale, one might either see nothing, or detect the small red cycle but miss the large yellow one, or vice versa. A solution to handle this issue is to characterize each cycle with its scale of appearance and disappearance: the red cycle induces a red bar in the barcode, and similarly for the yellow cycle.
(A) Simulations were generated with SANTA-SIM [114] for five distinct scenarios with varying growth rate g and mutation rate µ. The 95% prediction intervals for the number of one-dimensional cycles in each scenario are based on the extrapolation of a Panjer distribution for an increasing number of distinct sequences in the simulated phylogenies (see Methods). For each scenario, the validation dataset shown in the plot is well-described by the corresponding prediction intervals. The simulations demonstrate that for each time bucket, the number of topological cycles observed in GISAID data is significantly larger than the expected number of topological cycles that are due to noise (i.e. random topological cycles in simulated data). As of February 2021, scenario V suggests that less than 40% (1,171) of all cycles (2,899) in the GISAID alignment are due to noise. (B) Comparison of genetic distances to the root in simulated data vs. distances to the Wuhan/Hu-1 reference sequence EPI ISL 402125 in the GISAID dataset. Scenarios I, IV and V with low mutation rate systematically underestimate the maximal distance, while the highest mutation rate in scenario III yields larger distances. The mutation rate of scenario II describes the maximal distance and overall diversity well. Differences to the GISAID data are expected to be due to real-world effects like variation of population growth, belated up-take in sequencing efforts, and enhanced spread of certain variants. (C) Parameters and prediction intervals in scenarios I-V. Scenarios I-III vary over a range of mutation rates that roughly capture the diversity of the GISAID dataset. Scenarios IV and V probe the influence of logistic population growth. For all scenarios we produced 100 simulations for each of the values of the carrying population c 105, and five simulations for c = 106.
For each value of carrying population c we determined the mean and variation of the observations and used these as parameters for the Panjer distribution.
For each scenario we fit quadratic, cubic, powerlaw, and exponential models to the observed number of one-dimensional cycles in simulations. Then we fit a linear and powerlaw model to the corresponding residuals as an estimate for the variance of the data. The quality of each model is evaluated through the log-likelihood to observe the validation dataset given a certain model.
Topological recurrence index (tRI) vs. fitness index [7] for Spike gene amino acid changes as of March 2023. (A) We found a strong correlation between positive tRI and positive fitness index (Fisher’s exact test, p < 10e 10) with a Spearman correlation coefficient of 0.63 (p < 10e 99). Significant outliers, for which tRI and fitness index show divergent tendencies, are marked with red circles and are discussed in panel (B). The topological recurrence analysis is based on GISAID data [18, 36] covering the whole pandemic from December 2019 until March 2023, comprising 13,766,674 genetically distinct high-quality SARS-CoV-2 Spike genes (see Supplementary Information Table 5 and Methods). The data on fitness effects (fitness index) is taken from Bloom & Neher [108]. The statistical correlation analysis was performed for Spike gene amino acid changes listed in the file aa fitness.csv (see Methods). (B) Assessment of the potential adaptiveness for mutations marked with a red circle in (A) for which tRI and fitness index show significantly divergent tendencies. Among these, all mutations flagged with significant tRI 75 show signs of adaptation such as frequent occurrence in independent lineages, occurrence in Variants of Concern or experimental evidence. The assessment is based on phylogenetic data taken from Cov2Tree [119] and epidemiological data taken from covSPECTRUM [91], as well as experimental results from the literature.
Time series analysis charts for Omicron Spike gene amino acid changes with positive topological signal listed in panel (C) of Figure 9. Each chart shows the topological recurrence index (tRI, in red) and the tRI growth rate (14 days moving average, in blue) from January 2020 until January 2022 at daily resolution, where in the upper diagram the shaded region marks the level of significance (see Methods). The tRI growth rate is tampered by computational artifacts at the time when the tRI first attains a positive value, and towards the end of the covered period in January 2022 when the sampling rate in the underlying sequence alignment may still be low.
Acknowledgements
The authors gratefully acknowledge all data contributors, i.e. the Authors and their Originating laboratories responsible for obtaining the specimens, and their Submitting laboratories for generating the genetic sequence and metadata and sharing via the GISAID Initiative [18, 36], on which this research is based. An acknowledgement table can be retrieved from https://doi.org/10.55876/gis8.230731by. The authors acknowledge the use of de.NBI Cloud and the support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen and the German Federal Ministry of Education and Research (BMBF) through grant no 031 A535A. They thank M. Hanussek for IT support and early access to VALET [118]. The authors further acknowledge support from the Interdisciplinary Center for Scientific Computing at Heidelberg University and the development work of the Scientific Software Center of Heidelberg University carried out by L. Keegan and D. Kempf. This research was supported by the DFG Collaborative Research Center SFB/TRR 109 “Discretization in Geometry and Dynamics”. M.B. was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 - 390900948 (the Heidelberg STRUCTURES Excellence Cluster). L.H. thanks the Evangelisches Studienwerk Villigst for their support. M.N. was supported by the Vector Stiftung (“Topological Genomics”). A.O. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 281869850 (RTG 2229). This research was funded by the Federal Ministry of Education and Research (BMBF) and the Baden-Württemberg Ministry of Science as part of the Excellence Strategy of the German Federal and State Governments (KIT Centers, “Topological Genomics”), and by the Vector Stiftung (“Topological Genomics”).
Footnotes
Major revisions; new analyses added.
References
- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.
- 11.↵
- 12.
- 13.↵
- 14.
- 15.
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.
- 21.
- 22.
- 23.
- 24.
- 25.
- 26.
- 27.
- 28.
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.
- 58.
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.
- 64.
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.
- 98.
- 99.
- 100.
- 101.↵
- 102.↵
- 103.↵
- 104.↵
- 105.↵
- 106.↵
- 107.↵
- 108.↵
- 109.↵
- 110.↵
- 111.↵
- 112.↵
- 113.↵
- 114.↵
- 115.↵
- 116.↵
- 117.↵
- 118.↵
- 119.↵
