key: cord-0882256-hp6kt4ud
authors: Kurpas, Monika; Jaksik, Roman; Kimmel, Marek
title: Evolutionary analysis of genomes of SARS-CoV-2-related bat viruses suggests old roots, constant effective population size, and possible increase of fitness
date: 2022-03-01
journal: bioRxiv
DOI: 10.1101/2022.02.28.482287
sha: 0b0393716204fc42e19063d36afdcb77bb599378
doc_id: 882256
cord_uid: hp6kt4ud

It is of vital practical interest to understand the co-evolution of bat β-coronaviruses with their hosts, since a number of these most likely crossed the species boundaries and infected humans. Complete sequences of 47 consensus genomes are available for bat β-coronaviruses related to the SARS-CoV-2 human virus. We carried out several types of evolutionary analyses using these data. First, using the publicly available BEAST 2 software, we generated phylogenetic trees and skyline plots. The roots of the trees, both for the entire sequences and subsequences coding for the E and S proteins as well as the 5’ and 3’ UTR regions, are estimated to be located from several decades to more than a thousand years ago, while the effective population sizes remained largely constant. Motivated by this, we developed a simple estimator of the effective population size in a Moran model with constant population, which, under the model is equal to the expected age of the MRCA measured in generations. Comparisons of these estimates to those produced by BEAST 2 shows qualitative agreement. We also compared the site frequency spectra (SFS) of the bat genomes to those provided by the Moran Tug-of-War model. Comparison does not exclude the possibility that overall fitness of the bat β-coronaviruses was increasing over time as a result of directional selection. Stability of interactions of bats and their viruses was considered likely on the basis of specific manner in which bat immunity is tuned, and it seems consistent with our analysis.

It seems to be a predominant belief that the recent β-coronavirus (Severe Acute Respiratory Syndrome Coronavirus 2, SARS-CoV-2), which causes acute respiratory disease COVID-19, is of zoonotic origin [23] . Other hypotheses, including genetic engineering of the virus, are considered controversial in the research community [2] . For previous epidemics such as the Severe Acute Respiratory Syndrome (SARS) and Middle East Respiratory Syndrome (MERS) caused by β-coronaviruses, it has been proven that the reservoir of viruses capable of infecting humans is, among other, the population of bats [23] . It is an interesting question, what is the nature of the interaction of the population of the viruses with bat host.

In this study, we used consensus genomic sequences of bat β-coronaviruses to infer the pattern of their evolution and to estimate the relevant parameters. The simulation analysis was supported by a phylogenetic analysis carried out using the Bayesian Evolutionary Analysis by Sampling Trees 2 (BEAST 2) [4] software, which allows the creation of a phylogenetic tree, estimation of the age of the most recent common ancestor (MRCA) of the studied sequences and estimation of size of the studied population using the Bayesian Skyline algorithm. We cross-checked the sophisticated and comprehensive BEAST 2 estimates by using the maximum parsimony analysis [11, 12] . We also developed estimates of mutation rate and effective population size under the hypothesis of neutrality and under sampling that provided genomes captured at different dates.

We also proposed to use a version of the Tug-of-War model introduced by McFarland [18] , to model the site frequency spectra (SFS) of the population of viruses and, by comparison with the data, to further refine conclusions from phylogenetic analysis.

In addition, we carried out the analyses using the 3' and 5' UTR regions which may include regulatory elements and the genes coding for the envelope (E) and spike (S) proteins. This latter is a glycoprotein that binds to the ACE2 receptor present in the cytoplasmic membrane of the host cells and thus allows the virus to penetrate them [22] .

The analysis was carried out using 47 nucleotide sequences of SARS genomes isolated from bat hosts, downloaded from the National Center for Biotechnology Information (NCBI) Virus database [1] (full list of sequences can be found in Supplementary Information). Samples were dated from the period 2004-2017.

The sequences were aligned using MUSCLE multiple sequence aligner [9, 8] based on which we assessed the variability at each genomic position relative to the consensus sequence. For each variant we calculated the fraction of sequences in which an alternative allele was observed, and summarized the data for each gene. Gene positions inside the created consensus sequence were determined using local sequence alignment, based on coordinates of SARS genes obtained from the Reference Sequence (RefSeq) database under accession number NC 004718.3. RefSeq database, created by NCBI, is an open access, annotated and curated collection providing single records for each natural biological molecule (DNA, RNA or protein) for major organisms including viruses, bacteria and eukaryotes. Based on aligned samples, we deduced the ancestral sequence; see further on. Mutation frequency study was carried out in the R computer language, using seqinr, msa and Biostrings libraries and visualized using Gviz and ggplot2.

Individual parts of the genome (3'UTR, 5'UTR regions; envelope and spike gene coding fragments) were cut from the alignment based on Multiple Sequence Alignment (MSA) annotations (see Supplementary Information) . The length of the complete sequence in the alignment is 30821 nucleotides.

BEAST 2 (Bayesian Evolutionary Analysis Sampling Trees 2) is a cross-platform program for Bayesian phylogenetic analysis of molecular sequences. It uses Markov chain Monte Carlo (MCMC) method to estimate rooted, time-measured phylogenies. BEAST 2 uses MCMC algorithm to average over tree space, so that each tree is weighted proportional to its posterior probability.

BEAST's graphical user interface, BEAUti (Bayesian Evolutionary Analysis Utility), allows to specify the inference parameters and create .xml files used for calculation. We performed analyses using the following parameters in corresponding BEAUti booktabs:

• Partitions -allowing browsing uploaded nucleotide/amino acid sequences and division into partitions corresponding to reading frames.

In our case nucleotide sequences were either not partitioned (3'UTR, 5'UTR regions) or partitioned (gene E and S; whole genome) in accordance with the assumption of higher mutation rates of third nucleotide in the codon (split {1,2} + 3) with an appropriate reading frame.

We linked trees and molecular clocks of both partitions, to generate one set of phylogenetic trees. Site models remained independent.

• Tip Dates -allowing appropriate collection date assignment for each sequence.

• Site Model -allowing choosing substitution model for nucleotide/amino acid sequence. We used HKY (Hasegawa, Kishino, and Yano [16] ) substitution model with estimated substitution rate and empirical nucleotide frequencies.

• Clock Model -allowing choosing clock model type. In our inference we used the strict clock model. We analyzed the output from BEAST 2 using Tracer (version 1.7.1) [20] , which graphically and quantitively summarizes the distributions of continuous parameters and provides diagnostic information. We used Tracer also to obtain data for Bayesian Skyline plots. [11] .

This program carries out the unrooted maximum parsimony tree-making algorithm, analogous to the Wagner trees [17] , on DNA sequences. The method of Fitch [13] was used to count the number of base changes needed for a given tree. In the maximum parsimony method the phylogenetic tree that minimizes the total number of character-state changes is sought.

The assumptions of Dnapars method [10] are briefly listed. Each site evolves independently.

Different lineages evolve independently. The probability of a base substitution at a given site is small over the lengths of time involved in a branch of the phylogeny. The expected amounts of change in different branches of the phylogeny do not vary by so much that two changes in a highrate branch are more probable than one change in a low-rate branch. The expected amounts of change do not vary enough among sites that two changes in one site are more probable than one change in another.

' The input for this method was the MSA described earlier on. All parameters were left unchanged except the parameter 5: "Print sequences at all nodes of tree". As a result we obtained the sequences at tree nodes and the ancestral sequence.

The most-parsimonious tree often underestimates the evolutionary change that has occurred. In order to make tree estimation more reliable we used the bootstrap method seqboot from PHYLIP package. As an input MSA was used. As an output the set of 100 bootstrapped samples was obtained. In the next step it was used in the DNApars method in order to obtain tree estimates with the parameter "Analyze multiple data sets" enabled.

In the final step we checked the support for specific forks in consensus tree using Consense algorithm from PHYLIP package. The analysis was performed using "Majority Rule extended" parameter.

We used for simulations a version of the Moran Model, inspired by the Tug-of-War model of

McFarland [18] . Consider a population of a fixed number N of individuals, each of them characterized by a pair of integers γ i = (α i , β i ), corresponding to the numbers of drivers and passengers in its genotype, respectively. This pair determines the fitness x i of the i-th individual by the formula

where s > 0 and d ∈ (0, 1) are parameters describing selective advantage of driver mutations over passenger mutations.

This model is a version of the "Tug-of-War" model introduced by McFarland [18] . We provide hypotheses for the complete version. In this paper we will also use the neutral version, with

The entire population may be identified with the vector

The accompanying vector

of fitness values determines the future evolution of the population, seen as being under drift and selection pressure, as follows:

The time to the first selection/drift event for the entire Γ is exponentially distributed with

After this time, one individual dies and is replaced by an exact copy of one of the remaining individuals, the probability that the i-th individual dies and is replaced by the j-th (j = i) being x j (N −1)Σx . The process then continues with Γ modified by replacing its i-th coordinate by its j-th coordinate.

Moreover, each individual may, after an independent exponentially distributed time with parameter λ, and independently of other individuals, undergo a mutation event, changing its state to either (α + 1, β) or (α, β + 1) with probabilities p ∈ (0, 1) and q = 1 − p, respectively.

Mathematical properties of this model are subject of separate publications, in preparation. Here, we used stochastic simulations to explore predictions of this model. The following list (see also Fig. 1 ) summarizes the notation used. Expectation of the pairwise difference under assumption that l = k equals to

since the expected count of mutations in time T lk equals µT lk , while under Moran model and the ISM the expected time to the MRCA equals N , so that the expected count of pairwise differences between k and l, equals 2N µ.

Given sample values of T k , k = 1, 2, . . . , n and d lk , 1 ≤ k < l ≤ n, we obtain estimating equations for µ and 2N , by minimizing the sum of squares

Differentiating with respect to µ and 2N and equating the derivatives to 0, we further obtain

where Σd = 1≤k<l≤n d lk , ΣT = 1≤k<l≤n T lk , ΣdT = 1≤k<l≤n d lk T lk , and ΣT 2 = 1≤k<l≤n T 2 lk Statistical properties of the resulting estimators were investigated by simulation and proved satisfactory, although there exists a trade-off between estimates of the two parameters (see Appendix).

When applying the method to the selected subunits of the genome we performed sequence quality control to avoid large intervals of deletions or missing nucleotides. We detected such structures in 3'UTR and 5'UTR regions, possibly due to genetic material degeneration or sequencing errors. For this reason we excluded 4 sequences from the 3'UTR and 5 sequences from the 5'UTR dataset. The remaining sequences were trimmed so that they did not contain excess nucleotide deletions.

Inference from evolutionary models of DNA often exploits summary statistics of the sequence data, a common one being the so-called Site Frequency Spectrum. In a sequencing experiment with a known number of sequences, we can estimate for each site at which a novel somatic mutation has arisen, the number of genomes that carry that mutation. These numbers are then grouped into sites that have the same number of copies of a mutant. be present in all genomes in the sample; these are not segregating and are called truncal mutations.

Using MATLAB software we created observed site frequency spectrum of bat β-coronaviruses relative to the ancestral sequence obtained with use of PHYLIP package. Except standard nucleotide symbols, the ancestral sequence contained letters that (according to IUPAC rules [5] ) represent ambiguity and are used when more than one kind of nucleotide could occur at that position. It also contains "?" symbols, which represents that in addition to one or more bases, a deletion may or may not be present. While evaluating the SFS, we counted mutations in accordance with these rules -if the difference in the sequence at a given site fell within the area of ambiguity we did not count this case as a mutation. 

We prepared 100 bootstrapped samples using the seqboot method from PHYLIP package run on 47 aligned bat β-coronavirus sequences. Subsequently we used the Dnapars method from the PHYLIP package to create most parsimonious tree for each sample (in some cases more than one tree was found). In order to explore the created set of trees we used the consense method from PHYLIP package. The resulting consensus tree (Fig. 3) shows the most probable tree for this dataset. At each node, presented is the number of trees supporting the arrangement of branches displayed, which was determined using the extended majority rule (see Methods). 

We used BEAST 2 to obtain phylogenetic trees based on the whole genomes (Fig. 4) and four genome fragments ( Supplementary Information, Fig. S2 , S3, S4, and S5). Based on tip dates applied together with sequences the algorithm is able to estimate the age of most recent common ancestor of aligned sequences. Below we compare estimates for complete sequences (Fig. 5) with sequences trimmed in a way described in section 2.5 (Fig. 6 ). Whiskers mark the value range. Whiskers mark the value range.

Despite truncation of the sequences, overall MRCA age distribution is similar in both cases and the differences are not statistically significant -boxplots show overlapping interquartile range.

Based on genomic sequences generated as described in Section 2.5 we built the distance matrices and including collection dates for each one of them we calculated estimators using the least square (LS)based method. Table 1 includes the comparison between MRCA age estimated by BEAST 2 and by our method, which was used also to calculate mutation rate. Let us notice that the estimates of effective population sizes based on the two methods are quite different as are the inferred time scales (see Discussion). However, both methods suggest large age of the most recent common ancestors of the whole genomes as well as the 4 genomic subsets we considered. Also, estimated mutation rates from the least-square method, suggest high conservation of some parts of the genome (the UTR regions and the envelope gene) as opposed to fast-mutating spike gene. Interpretation of these estimates is suggested in the Discussion section. 

We compared our calculated effective population size with population size estimated by Bayesian Skyline algorithm. In line with our results, BEAST 2 predicted higher effective population size for the whole genome and for the spike gene coding fragment and than for envelope gene or 3' and 5'

UTR regions (Table 1 and Fig. S6 ).

In order to check whether the trimming of deletions or missing nucleotides did not affect qualitatively the estimates, we compared Bayesian Skyline plots for complete (Fig. S6 in the Supplement) and trimmed (Fig. S7 in the Supplement) regions of genome. Our results do not show significant changes of shape or estimated values in both cases.

Bayesian Skyline plots in most cases indicated no significant changes in population size over time. We do not observe sharp changes or bottlenecks in the population size.

The Moran Tug-of-War model [18, 7] described in Methods, which includes the stochastic processes of mutation, genetic drift and selection, allows mapping of the evolutionary dynamics occurring in real populations. As explained, some mutations are "driver" mutations that increase fitness of individuals, and other are "passenger" mutations that are neutral or decrease fitness.

In the simulations we first used the neutral version of the model, in which there was no selective advantage of driver mutations over passenger mutations, with parameters s=d=0 and p=0.5, so 

In this study, we analyzed the genomic sequences of bat coronaviruses, using general use phylogenetic software packages PHYLIP and BEAST 2 as well as the least-square (LS) method and compared these results to predictions of different theoretical models.

As seen in Table 1 , the estimated effective population size is smaller for the 3' UTR and 5' UTR regions and the envelope gene than for the spike gene, what implicates that the former are more conserved than the latter. This is consistent with the hypothesis that UTRs play a distinctive role in gene regulation, but also with higher mutability with spike, which is seen in the human SARS- In bats, the virus-host interaction is regulated by a distinctive action of the innate immune system, different from that in other mammals [3] . Interestingly, the fragments of viral genomes are often permanently inserted into the genomes of their carriers, including the regions undergoing expression and recombination. It is hypothesized that stretches of viral genetic material found in tolerant to pathogens, which may explain why they are considered one of the major reservoirs of viral diseases. As in the organisms of bats, also in humans, the genomes of the virus undergo constant mutations. Research carried out on SARS-CoV-2 sequences has shown that despite the short time that has elapsed since the outbreak of the pandemic, the genetic material contains mutations specific for different regions of the world [14] . Therefore, it is crucial to understand the mechanisms that govern both virus evolution (increasing adaptation to new hosts) and virus-host co-evolution.

Accession numbers of 47 genomic sequences of bat SARS-CoV-related β-coronaviruses analyzed. 

We carried out a number of simulations to illustrate the statistical properties of the least-square method. We repeated simulations involving exclusion of random 20% of sequences from the test sample and calculation of the least square estimators based on the remaining ones. These studies show a negative correlation between the estimates of mutation rate and effective population size, this latter equal to the time to MRCA. This effect is expected if the time to MRCA is much higher than the time span of the observations.

The quantity which is best estimable seems to be θ = 2N µ. The following graphs (Fig. S8) illustrate results obtained for subsequences of the viral genomes. 

The proximal origin of SARS-CoV-2

Novel insights into immune systems of bats

BEAST 2: a software platform for Bayesian evolutionary analysis

Nomenclature for incompletely specified bases in nucleic acid sequences: recommendations 1984

Statistical inference for the evolutionary history of cancer genomes

Probability models for DNA sequence evolution

MUSCLE: a multiple sequence alignment method with reduced time and space complexity

MUSCLE: multiple sequence alignment with high accuracy and high throughput

Maximum likelihood and minimum-steps methods for estimating evolutionary trees from data on discrete characters

PHYLIP (phylogeny inference package), version 3.5 c

Toward defining the course of evolution: minimum change for a specific tree topology

Phylogenetic network analysis of SARS-CoV-2 genomes

The age of a mutation in a general coalescent tree

Dating of the human-ape splitting by a molecular clock of mitochondrial DNA

Quantitative phyletics and the evolution of anurans

Tug-of-war between driver and passenger mutations in cancer and other adaptive processes

Figtree-version 1.4.3, a graphical viewer of phylogenetic trees

Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Systematic Biology

The potential role of endogenous viral elements in the evolution of bats as reservoirs for zoonotic viruses

Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2

A pneumonia outbreak associated with a new coronavirus of probable bat origin

BKM-726/RAU1/2020) granted by Polish Ministry of Science and Higher Education. Roman Jaksik was supported by 2016/23/D/ST7/03665 grant from Polish National Science Centre. Marek Kimmel was supported by the NSF/DMS Rapid Collaborative grant to

All authors reviewed and approved the final version.

The authors declare that they have no competing interests.

All relevant data are within the manuscript and its Supporting Information files.