Using multiple sampling strategies to estimate SARS-CoV-2 epidemiological parameters from genomic sequencing data
=================================================================================================================

* Rhys P. D. Inward
* Kris V. Parag
* Nuno R. Faria

## ABSTRACT

SARS-CoV-2 virus genomes are currently being sequenced at an unprecedented pace. The choice of viral sequences used in genetic and epidemiological analysis is important as it can induce biases that detract from the value of these rich datasets. This raises questions about how a set of sequences should be chosen for analysis, and which epidemiological parameters derived from genomic data are sensitive or robust to changes in sampling. We provide initial insights on these largely understudied problems using SARS-CoV-2 genomic sequences from Hong Kong, China, and the Amazonas State, Brazil. We consider sampling schemes that select sequences uniformly, in proportion or reciprocally with case incidence and which simply use all available sequences (unsampled). We apply *Birth-Death Skyline* and *Skygrowth* methods to estimate the time-varying reproduction number (*R**t*) and growth rate (*r**t*) under these strategies as well as related *R*** and date of origin parameters. We compare these to estimates from case data derived from *EpiFilter*, which we use as a reference for assessing bias. We find that both *R**t* and *r**t* are sensitive to changes in sampling whilst *R*** and the date of origin are relatively robust. Moreover, we find that analysis using unsampled datasets, which reflect an opportunistic sampling scheme, result in the most biased *R**t* and *r**t* estimates for both our Hong Kong and Amazonas case studies. We highlight that sampling strategy choices may be an influential yet neglected component of sequencing analysis pipelines. More targeted attempts at genomic surveillance and epidemic analyses, particularly in settings with limited sequencing capabilities, are necessary to maximise the informativeness of virus genomic datasets.

## INTRODUCTION

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an enveloped single-stranded zoonotic RNA virus belonging to the *Betacoronavirus* genus and *Coronaviridae* family1. It was first identified in late 2019 in a live food market in Wuhan City, Hubei Province, China2. Within a month, SARS-CoV-2 had disseminated globally through sustained human-to-human transmission. It was declared a public health emergency of international concern on the 30th of January 2020 by the World Health Organisation3. Those infected with SARS-CoV-2 have phenotypically diverse symptoms ranging from mild fever to multiple organ dysfunction syndromes and death4.

Despite the implementation of non-pharmaceutical interventions (NPIs) and rollout of vaccination programmes in many countries to control their epidemics, as of the 16th July 2022, over 557 million SARS-CoV-2 cases and 6.3 million deaths have been reported worldwide5. These NPIs can vary within and between countries and include restrictions on international and local travel, school closures, social distancing measures and the isolation of infected individuals and their contacts6. The key aim of NPIs is to reduce epidemic transmission, often measured by epidemiological parameters such as the time-varying effective reproduction number (*R**t* at time *t*) and growth rate (*r**t*), which both provide updating measures of the rate of spread of a pathogen (see Table 1 for detailed definitions)7,8.

View this table:
[Table 1:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/T1)

Table 1: 
Key parameters and definitions for SARS-CoV-2

However, there is currently great difficulty in estimating and comparing epidemiological parameters derived from case and death data globally due to disparities in molecular diagnostic surveillance and notification systems between countries. Further, even if data are directly comparable, the choice of epidemiological parameter can implicitly shape insights into how NPIs influence transmission potential9,10. As such, there is a need to supplement traditional estimates with information derived from alternative data sources, such as genomic data11, to gain improved and more robust insights into viral transmission dynamics12,13.

Phylodynamic analysis of virus genome sequences have increasingly been used for studying emerging infectious diseases, as seen during the current SARS-CoV-2 pandemic14–17, recent Ebola virus epidemics in Western Africa18 and the Zika epidemic in Brazil and the Americas19,20. Transmissibility parameters such as the basic reproduction number (*R***), *R**t* and *r**t* can be directly inferred from genomic sequencing data or from epidemiological data, while other epidemiological parameters such as the time of the most recent common ancestor (TMRCA) of a given viral variant or lineage can only be estimated from genomic data. This is of particular importance for variants of concern (VOC), genetic variants with evidence of increased transmissibility, more severe disease, and/or immune evasion. VOC are typically detected through virus genome sequencing and only limited inferences can be made using epidemiological data alone21.

Currently, SARS-CoV-2 virus genomes from COVID-19 cases are being sequenced at an unprecedented pace providing a wealth of virus genomic datasets22. There are currently over 11.9 million genomic sequences available on GISAID, an open-source repository for influenza and SARS-CoV-2 genomic sequences23. These rich datasets can be used to provide an independent perspective on pathogen dynamics and can help validate or challenge parameters derived from epidemiological data. Specifically, the genomic data can potentially overcome some of the limitations and biases that can result from using epidemiological data alone. For example, genomic data are less susceptible to changes at the government level such as alterations to the definition of a confirmed case and changes to notification systems24,25. Inferences from virus genomic data improve our understanding of underlying epidemic spread and can facilitate better-informed infection control decisions26. However, these advantages are not straightforward to realise. The added value of genomic data depends on two related variables: sampling strategy and computational complexity.

The most popular approaches used to investigate changes in virus population dynamics include the Bayesian Skyline Plot27 and Skygrid28 models and the Birth-Death Skyline (BDSKY)29. These integrate Markov Chain Monte Carlo (MCMC) procedures and often converge slowly on large datasets30. As such, currently available SARS-CoV-2 datasets containing thousands of sequences become computationally impractical to analyse and sub-sampling is necessary. Although previous studies have examined how sampling choices might influence phylodynamic inferences30–34, this remains a neglected area of study35, particularly concerning SARS-CoV-2 for which sequencing efforts have been unprecedented 36. To our knowledge, there are no published studies concerning SARS-CoV-2 which explore the effect that sampling strategies can have on the phylodynamic reconstruction of key transmission parameters. Incorrectly implementing a sampling scheme or ignoring its importance can mislead inferences and introduce biases30,37.

This raises the important question that motivates our analysis: how should sequences be selected for phylodynamic analysis and which parameters are sensitive or robust to changes in different sampling schemes. Here we explore how diverse sampling strategies in genomic sequencing may affect the estimation of key epidemiological parameters. We estimate *R***, *R**t*, and *r**t* from genomic sequencing data under different sampling strategies from a location with higher genomic coverage represented by Hong Kong, and a location with lower genomic coverage represented by the Amazonas state, Brazil. We then compare our estimates against those derived from reference case data. By benchmarking genomic inferences against those from case data we can better understand the impact that sampling strategies may have on phylodynamic inference, bolster confidence in estimates of genomic-specific parameters such as the origin time (or TMRCA) and improve the interpretation of estimates from areas with heterogeneous genomic coverage.

## METHODS

### Empirical Estimation of the Reproduction Number, Time-varying Effective Reproduction Number, and Growth Rate

#### Epidemiological Datasets

Two sources of data from the Amazonas state, Brazil, and one source of data from Hong Kong were used to calculate empirical epidemiological parameters. For the Amazonas state, case data from the SIVEP-Gripe (Sistema de Informação de Vigilância Epidemiológica da Gripe) SARI (severe acute respiratory infections) database from the 30th of November 2020 up to 7th of February 2021 were used. Here we were interested in cases caused by the Gamma VOC first detected in Manaus14. The number of Gamma cases was calculated by using the proportion of Gamma viral sequences uploaded to GISAID within each week (Supplementary Figure 1). For Hong Kong, all case data were extracted from the Centre of Health Protection, Department of Health, the Government of the Hong Kong Special Administrative region up to the 7th of May 2020. Due to lags in the development of detectable viral loads, symptom onset and subsequent testing38; the date on which each case was recorded was left shifted by 5 days within our models39 to account for these delays in both datasets.

#### Basic Reproduction Number

The *R*** parameter was estimated using a time series of confirmed SARS-CoV-2 cases from both Hong Kong and the Amazonas state. To avoid the impact of NPIs, only data up to the banning of mass gathering in Hong Kong (27th March 2020) and until the imposition of strict restrictions in the Amazonas state (12th January 2021) were used. We estimated *R*** from weekly counts of confirmed cases using maximum likelihood methods. The weekly case counts were assumed to be Poisson distributed and were fitted to a closed Susceptible-Exposed-Infectious-Recovered (SEIR) model (equation (1)) by maximising the likelihood of observing the data given the model parameters (Table 2). Subsequently, the log-likelihood was used to calculate the R by fitting β, the effective contact rate.

View this table:
[Table 2:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/T2)

Table 2: 
This shows the parameter estimates used within the deterministic SEIR model.

![Formula][1]</img>  To generate approximate 95% confidence intervals (CIs) for *R***, non-parametric bootstrapping was used with 1000 iterations.

#### Time-varying Effective Reproduction Number

To estimate *R**t* from case time series data the *EpiFilter* method44 was used. *EpiFilter* describes transmission using a renewal model; a general and popular framework that can be applied to infer the dynamics of numerous infectious diseases from case incidence45. This model describes how the number of new cases (incidence) at time *t* depends on *R**t* at that specified time point and the past incidence, which is summarised by the cumulative number of cases up to each time point weighted by the generation time distribution. *EpiFilter* integrates both Bayesian forward and backward recursive smoothing. This improves *R**t* estimates by leveraging the benefits of two of the most popular *R**t* estimation approaches: *EpiEstim* 46 and the Wallinga-Teunis method47. *EpiFilter* minimises the mean squared error in estimation and reduces dependence on prior assumptions, making it a suitable candidate for deriving reference estimates. We use these to benchmark estimates independently obtained from genomic data. We assume the generation time distribution is well approximated by the serial interval (SI) distribution46, which describes the times between symptom onsets between an infector–infectee pair. We describe the specific SI distributions that we used in the next subsection.

#### Growth Rate

After *R**t* has been inferred, the Wallinga-Lipsitch equation for a gamma distributed generation time distribution (equation (2)) was used to estimate the exponential epidemic *r**t* 48. The SI for Hong Kong was derived from a systematic review and meta-analysis49 and a study exploring SI in Brazil was used for the Amazonas datasets50. The SI was assumed to be gamma distributed. The gamma distribution is represented by gamma (ε, γ) with ε *and* γ being the shape and scale parameters respectively. ![Formula][2]</img>  

### SARS-CoV-2 Brazilian Gamma VOC and Hong Kong datasets

All high-quality (<1% N, or non-identified nucleotide), complete (>29 kb) SARS-CoV-2 genomes were downloaded from GISAID23 for Hong Kong (up to 7th May 2020) and the Amazonas state, Brazil (from 30th November 2020 up to 7th February 2021). Using the Accession ID of each sequence, all sequences were screened and only sequences previously analysed and published in PubMed, MedRxiv, BioRxiv, virological or Preprint repositories were selected for subsequent analysis. For both datasets, sequence alignment was conducted using MAFFTV.751. The first 130 base pairs (bp) and last 50 bps of the aligned sequences were trimmed to remove potential sequencing artefacts in line with the Nextstrain protocol52. Both datasets were then processed using the Nextclade pipeline for quality control ([https://clades.nextstrain.org/](https://clades.nextstrain.org/)). Briefly, the Nextclade pipeline examines the completeness, divergence, and ambiguity of bases in each genetic sequence. Only sequences deemed ‘good’ by the Nextclade pipeline were selected. Subsequently, all sequences were screened for identity and in the case of identical sequences, for those with the same location, collection date, only one such isolate was used. Moreover, PANGO lineage classification was conducted using the Pangolin22 (v2.3.3) software tool ([http://pangolin.cog-uk.io](http://pangolin.cog-uk.io)) on sequences from the Amazonas state and only those with the designated P.1/Gamma lineage were selected (Supplementary Figure 1).

### Maximum Likelihood tree reconstruction

Maximum likelihood phylogenetic trees were reconstructed using IQTREE253 for both datasets. A TIM2 model of nucleotide substitution with empirical base frequencies and a proportion of invariant sites was used as selected for by the ModelFinder application54 for the Hong Kong dataset. For the Brazilian dataset, a TN model of nucleotide substitution55 with empirical base frequencies was selected for. To assess branch support, the approximate likelihood-ratio test based on the Shimodaira–Hasegawa-like procedure with 1,000 replicates56, was used.

### Root-to-tip regression

To explore the temporal structure of both the Brazilian and Hong Kong dataset, TempEst v.1.5.357 was used to regress the root-to-tip genetic distances against sampling dates (yyyy-mm-dd). The ‘best-fitting’ root for the phylogeny was found by maximising the R2 value of the root-to-tip regression (Supplementary Figure 2). Several sequences showed incongruent genetic diversity and were discarded from subsequent analyses. This resulted in a final dataset of N = 117 Hong Kong sequences and N = 196 Brazilian sequences. The gradient of the slopes (clock rates) provided by TempEst were used to inform the clock prior in the phylodynamic analysis.

### Subsampling for analysis

Four retrospective sampling schemes were used to select a subsample of Amazonas and Hong Kong sequences. Each sampling period was broken up into weeks with each week being used as an interval according to a temporal sampling scheme (without replacement). This temporal sampling scheme was based on the number of reported cases of SARS-CoV-2.

The temporal sampling schemes that we explored were:

*   **No sampling strategy applied:** All sequences were included without a sampling strategy applied (equivalent to opportunistic sampling).

*   **Proportional sampling**: Weeks are chosen with a probability proportional to the value of the number of incident cases in each epi-week.

*   **Uniform sampling:** All weeks have equal probability.

*   **Reciprocal-proportional sampling:** Weeks are chosen with a probability proportional to the reciprocal of the incident number of cases in each epi-week.

These sampling schemes were inspired by those recommended by the WHO for practical use in different settings and scenarios58. Proportional sampling is equivalent to representative sampling, uniform sampling is equivalent to fixed sampling whilst the unsampled data includes all sampled sequences. Reciprocal-proportional sampling is not commonly applied in practice and was used as a control within this study.

### Bayesian Evolutionary Analysis

Date molecular clock phylogenies were inferred for all sampling strategies applied to the Amazonas and Hong Kong dataset using BEAST v1.10.459 with BEAGLE library v3.1.060 for accelerated likelihood evaluation. For both the Amazonas and Hong Kong datasets, a HKY substitution model with gamma-distributed rate variation among sites and four rate categories was used to account for among-site rate variation61. A strict clock molecular clock model was chosen. Both the Amazonas and Hong Kong dataset were analysed under a flexible non-parametric skygrid tree prior62. Four independent MCMC chains were run for both datasets.

For the Amazonas dataset, each MCMC chain consisted of 250,000,000 steps with sampling every 50,000 steps. Meanwhile, for the Hong Kong dataset, each MCMC chain consisted of 200,000,000 steps with sampling every 40,000 steps. For both datasets, the four independent MCMC runs were combined using LogCombiner v1.10.459. Subsequently, 10% of all trees were discarded as burn in, and the effective sample size of parameter estimates were evaluated using TRACER v1.7.263. An effective sample size of over 200 was obtained for all parameters. Maximum clade credibility (MCC) trees were summarised using Tree Annotator59.

### Phylodynamic Reconstruction

#### Estimation of the Basic and Time-varying Effective Reproduction Numbers

The Bayesian birth-death skyline (BDSKY) model29 implemented within BEAST 2 v2.6.564 was applied to estimate the time-varying transmissibility parameter *R**t* (Table 3). An HKY substitution model with a gamma-distributed rate variation among sites and four rate categories61 was used alongside a strict molecular clock model. The selected number of intervals for both datasets was 5, representing *R**t* changing every 2.5 weeks for the Hong Kong datasets and every 2 weeks for the Brazilian datasets, with equidistant intervals per step. An exponential distribution was used with a mean of 36.5y-1 for the rate of becoming infectious, assuming a mean duration of infection of 10 days15. A uniform distribution prior was used for the sampling proportion, which models changes in case ascertainment. Four independent MCMC chains were run for 50 million MCMC steps with sampling every 5000 steps for each dataset. These MCMC runs were combined using LogCombiner v2.6.5.64 and the effective sample size of parameter estimates evaluated using TRACER v1.7.263. We obtained an effective sample size above 200 for all parameters (indicating convergence) and plotted all results using the bdskytools R package ([https://github.com/laduplessis/bdskytools](https://github.com/laduplessis/bdskytools)).

View this table:
[Table 3:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/T3)

Table 3: 
Values and priors for the parameters of the BDSKY model. s/s/y=substitutions per site per year.

#### Estimation of Growth Rates

For each dataset, a scaled proxy for *r**t* was obtained from the *Skygrowth* method65 within R. *Skygrowth* uses a non-parametric Bayesian approach to apply a first-order autoregressive stochastic process on the growth rate of the effective population size. The MCMC chains were run for one million iterations for each dataset on their MCC tree with an Exponential (10−5) prior on the smoothing parameter. The *Skygrowth* model was parameterised assuming that the effective population size of SARS-COV-2 could change every two weeks. To facilitate a comparison of the scaled proxy for *r**t* estimated by *Skygrowth* and exponential *r**t* estimated by *EpiFilter*, the *r**t* estimated by the *Skygrowth* method was rescaled to the exponential growth rate. This was achieved by adding a gamma rate variable to the scaled proxy for *r**t*, which assumed a mean duration of infection of 10 days15, to calculate *R**t*. Subsequently, the Wallinga-Lipsitch equation (Equation 2) was used to convert *R**t* into the exponential growth rate48.

### Comparing Parameter Estimates from Genetic and Epidemiological Data

To compare estimates derived from epidemiological and genetic data the Jensen-Shannon divergence (DJS)66, which measures the similarity between two probability mass functions (PMFs), was applied. The DJS offers a formal information theoretic evaluation of distributions and is more robust than comparing Bayesian credible intervals (BCIs) since it considers both the shape and spread of a given distribution. The DJS is a symmetric and smoothed version of the Kullback-Leibler divergence (DKL) and is commonly used in the fields of machine learning and bioinformatics. The DKL between two PMFs, P and Q, is defined in equation (3) below67, with *x* spanning the common domain of those PMFs. ![Formula][3]</img>  To calculate the PMF for each epidemiological parameter, the cumulative probability density function was extracted for each model, converted to a probability density function and a discretisation procedure was applied to generate the associated PMF.

The Jensen-Shannon distance (JSD) is a metric which takes the square-root of the total DJS and is the metric that we used to compare parameter estimations from differing sampling strategies. The JSD can be calculated using Equation 4 with P and Q representing the two probability distributions and DKL as the KL divergence. A smaller JSD metric indicates that two probability distributions (P and Q) are more similar with a Jensen-Shannon distance of 0 uniquely indicating that both distributions are equivalent. The mean JSD was taken over all intervals for the BDSKY and *Skygrowth* models to obtain an overall measure of the level of estimated similarity across the epidemic trajectory. We do not expect the JSD to perfectly align with the 95% highest posterior density interval if the shapes of distributions from different schemes are very different. ![Formula][4]</img>  

### Data availability

Code and data for reproducing the analyses presented in this study are freely available at [https://github.com/rhysinward/Phylodynamic-Subsampling](https://github.com/rhysinward/Phylodynamic-Subsampling).

## RESULTS

### Sampling Schemes

#### Hong Kong

Hong Kong reacted rapidly upon learning of the emergence of SARS-CoV-2 in Wuhan, Hubei province, China, by declaring a state of emergency on the 25th of January 2020 and by mobilising intensive surveillance schemes in response to initial cases68. This appeared to be successful in controlling the first wave of cases. However, due to imported cases from Europe and North America, a second wave of SARS-CoV-2 infections emerged prompting stricter NPIs such as the closure of borders and restrictions on gatherings 68. Following these measures, the incidence of SARS-CoV-2 rapidly decreased (Figure 1). Hong Kong has a high sampling intensity with 11.6% of confirmed cases sequenced during our study period. Further, Hong Kong has high quality case data with a high testing rate through effective tracing of close contacts, testing of all asymptomatic arriving travellers and all patients with pneumonia69.

![Figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F1.medium.gif)

[Figure 1.](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F1)

Figure 1. 
Confirmed incident SARS-CoV-2 cases from Hong Kong until 7th of May 2020. The arrows represent policy change-times68.

The number of cases within Hong Kong for each week was used to inform the sampling schemes used within this study. This resulted in the unsampled scheme having N = 117 sequences, the proportional sampling scheme having N = 54 sequences, the uniform sampling scheme having N = 79 and the reciprocal-proportional sampling scheme having N = 84 sequences (Supplementary Figure 3).

#### Amazonas

The Amazonas state of Brazil had its first laboratory confirmed case of SARS-CoV-2 in March 2020 in a traveller returning from Europe70. After a first large wave of SARS-CoV-2 infections within the state that peaked in early May 2020 (Figure 2), the epidemic waned, cases dropped, remaining stable until mid-December 2020. The number of cases then started growing exponentially, ushering in a second epidemic wave. This second wave peaked in January 2021 (Figure 2) and coincided with the emergence of a new SARS-CoV-2 VOC, designated P.1/Gamma14.

![Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F2.medium.gif)

[Figure 2.](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F2)

Figure 2. 
Confirmed incident SARS-CoV-2 cases from Amazonas state, north Brazil until 7th of February 2021. The arrows represent key policy change-times72.

To combat this second wave, the Government of the Amazonas state suspended all non-essential commercial activities on the 23rd of December 2020 ([http://www.pge.am.gov.br/legislacao-covid-19/](http://www.pge.am.gov.br/legislacao-covid-19/)). However, in response to protests, these restrictions were reversed, and cases continued to climb. On the 12th of January, when local transmission of P.1/Gamma was confirmed in Manaus, capital of Amazonas state71, NPIs were re-introduced ([http://www.pge.am.gov.br/legislacao-covid-19/](http://www.pge.am.gov.br/legislacao-covid-19/)) which seemed to be successful in reducing the case incidence in the state. However, cases remained comparatively high (Figure 4). Amazonas has a sampling intensity with 2.4% of suspected P.1/gamma cases sequenced during our study period.

The number of cases within the Amazonas state informed the sampling schemes used within this study. This resulted in the unsampled scheme having N = 196 sequences, the proportional sampling scheme having N = 168 sequences, the uniform sampling scheme having N = 150 and the reciprocal-proportional sampling scheme having N = 67 sequences (Supplementary Figure 4).

### Root-to-tip Regression

The correlation (R2) between genetic divergence and sampling dates for the Hong Kong datasets ranged between 0.36 and 0.52 and between 0.13 and 0.20 for the Amazonas datasets (Supplementary Figure 2). This implies that the Hong Kong datasets have a stronger temporal signal. This is likely due to the Hong Kong datasets having a wider sampling interval (106 days) compared to the Amazonas datasets (69 days). A wider sampling interval can lead to a stronger temporal signal73. The gradient (rate) of the regression ranged from 1.16×10−3 to 2.09×10−3 substitutions per site per year (s/s/y) for the Hong Kong datasets and 4.41×10−4 to 5.30×10−4 s/s/y for the Amazonas datasets.

### Estimation of Evolutionary Parameters

The mean substitution rate (measured in units of number of s/s/y) and the TMRCA was estimated in BEAST, for both datasets, and the estimation from all sampling schemes was compared.

#### Hong Kong

For Hong Kong, the mean substitution rate per site per year ranged from 9.16×10−4 to 2.09×10−3 with sampling schemes all having overlapping Bayesian credible intervals (BCIs) (Supplementary table 2; Supplementary Figure 5A). This indicates that the sampling scheme did not have a significant impact on the estimation of the clock rate. Moreover, the clock rate is comparable to estimations from the root-to-tip regression and to early estimations of the mean substitution rate per site per year of SARS-CoV-2 (Duchene et al., 2020).

Molecular clock dating of the Hong Kong dataset indicates that the estimated time of the most common recent ancestor was around December 2020 (Figure 3B; Supplementary Table 2). This is a few weeks before the first confirmed case which was reported on the 18th of January 2021. Once again, all sampling strategies have overlapped BCIs and with the range in means differing by around three weeks, a relatively short time scale, suggesting that the sampling scheme does not significantly impact the estimation of the TMRCA.

![Figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F3.medium.gif)

[Figure 3.](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F3)

Figure 3. 
*R*** estimated from BDSKY (using sequence data) and TMRCA for Hong Kong and Brazil. Figure 1A and B represent Hong Kong and Figure 1C and D represent the Amazonas. The central line represents the posterior mean and with intervals representing 95% highest posterior density interval.

#### Brazil

For the Gamma VOC in the Amazonas state, the mean substitution rate ranged from 4.00×10−4 to 5.56×10−4 s/s/y with all sampling schemes having overlapped BCIs (Figure 3D, Supplementary Table 2; Supplementary Figure 5B). This indicates that sampling strategy does not impact the estimation of the clock rate, supporting findings from the Hong Kong dataset. This also supports estimations from the root-to-tip analysis (Supplementary Figure 2).

Molecular clock dating estimated a TMRCA mean around late October to early November (Figure 3D; Supplementary Table 2). This is around five weeks before the date of the first P.1 case identified in Manaus used in our study. All sampling schemes have overlapping BCI consistent with the conclusion from the Hong Kong data that TMRCA is relatively robust to sampling.

### Estimation of Basic Reproduction Number

We found from using genomic data, Hong Kong had a posterior mean *R*** estimate of 2.07 (Figure 3A) across all sampling strategies. Using a proportional sampling strategy gave the highest posterior mean *R*** estimate of 2.38 with the unsampled sampling strategy giving the lowest posterior mean *R*** estimate of 1.87. Overall, Brazil had a higher posterior mean *R*** estimate with a value of 2.24 (Figure 3B) across all sampling strategies. The uniform sampling strategy yielded the highest posterior mean *R*** estimate of 2.50 while the unsampled sampling strategy gave the lowest one of 1.82. Using case data, we found similarly found that Hong Kong had a lower *R*** of 2.17 (95% credible interval (CI) = 1.43 - 2.83) when compared to Amazonas which had a *R*** of 3.67 (95% CI = 2.83 – 4.48). All sampling schemes for both datasets were characterised by similar *R*** values (Figure 3) indicating that the estimation of *R*** is robust to changes in sampling scheme.

### Time-varying Reproduction number and Growth rate

We estimate *R**t* and *r**t* for local SARS-CoV-2 epidemics in Hong Kong and Amazonas, Brazil. Our main results showing these two parameters and JSD metrics are shown in Figures 4-8.

![Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F4.medium.gif)

[Figure 4:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F4)

Figure 4: 
*R**t* estimated from both the BDSKY and *EpiFilter* methods for Hong Kong. Titles indicate the sampling scheme used in panels A-D. The light-shaded area represents the 95% highest posterior density interval. The solid line represents the mean *R**t* estimate with *EpiFilter* in green and BDSKY in blue. The black line plots the number of cases. We refer to Figure 1 for a brief description of key events 1–3. The Jensen Shannon Distance (JSD) is given in panel E and ranks the sampling strategies based on how similar the BDSKY estimates under those strategies are to those derived from *EpiFilter* (smaller values are better). Panel F provides the pairwise JSD between the BDSKY estimates under different sampling strategies, showing often appreciable difference among strategies.

#### Hong Kong

We applied the BDSKY model to estimate the *R**t* for each dataset subsampled according to the different sampling strategies (Figure 4). We compared these against the *R**t* from incidence data, derived from *EpiFilter*. Based on the proportional sampling scheme, which had the lowest JSD (Figure 4E), we initially infer a super-critical *R**t* value, with a mean around 2, that appears to fall swiftly in response to the state of emergency and the rapid implementation of NPIs. A steady transmission rate subsequently persisted throughout the following weeks around the critical threshold (*R**t* = 1). This period is followed by a sharp increase in *R**t*, peaking at a mean value of 2.6. This is likely due to imported cases from North America and Europe68. This led to a ban on international travel resulting in a sharp decline in *R**t* (Figure 2). However, this decline lasted around a week with the mean *R**t* briefly increasing until more stringent NPIs such as the banning of major gatherings were implemented. Following this, the *R**t* continued its sharp decline falling below the critical threshold, with transmission becoming sub-critical (Figure 4). The proportional sampling scheme showed the most divergence from all other sampling schemes whilst the uniform and reciprocal-proportional sampling schemes were almost identical (Figure 4F).

These results were mirrored in the estimation of *r**t*. (Figure 5), where estimates derived from the proportional sampling scheme showed the least divergence (Figure 5E). There was an initial decline in the *r**t*, which steadied at a value of ∼ 0, indicating that epidemic stabilisation had occurred. This stable period is followed by an increase in *r**t* peaking at around a 0.050 d-1 (Figure 5B). In response to NPIs, the *r**t* starts to decrease, falling below 0, indicating a receding epidemic. The rate of this decline peaks at around -0.075 d-1 (Figure 5B). Unlike the estimation of *R**t* (Figure 4), the unsampled sampling scheme showed the most divergence from all other sampling schemes (Figure 5F). It also has a high divergence from estimates derived from *EpiFilter* when compared the proportional sampling scheme which was the most closely related to *EpiFilter* (Figure 5E). Once again, the uniform and reciprocal-proportional schemes are the most closely related (Figure 5E).

![Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F5.medium.gif)

[Figure 5:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F5)

Figure 5: 
*r**t* estimated from both the *Skygrowth* and *EpiFilter* methods for Hong Kong. Titles indicate the sampling scheme used in panels A-D. The light-shaded area represents the 95% highest posterior density interval. The solid line represents the mean *r**t* estimate with *EpiFilter* in orange and *Skygrowth* in blue. The black line refers to the number of cases. We refer to Figure 1 for a brief description of key events 1–3. The Jensen Shannon Distance (JSD) is given in panel E and ranks the sampling strategies based on how similar the *Skygrowth* estimates under those strategies are to those derived from *EpiFilter* (smaller values are better). Panel F provides the pairwise JSD between the BDSKY estimates under different sampling strategies, showing often appreciable difference among strategies.

#### Brazil

The uniform, reciprocal-proportional, and proportional sampling schemes all showed a similarly low JSD (Figure 6E). Based on these sampling schemes, we initially infer super-critical transmission (*R**t* > 1) with a mean value of 3 (Figure 6). From this point, the *R**t* declines, although it remains above the critical threshold (*R**t* = 1) for much of the study period. Sub-critical transmission (*R**t* < 1) was only reached after the re-imposition of NPIs. This implies that initial restrictions, such as the suspension of commercial activities, were likely insufficient for suppressing spread. Only after more stringent restrictions were imposed did *R**t* become sub-critical. However, there is no evidence of a sharp decrease in *R**t* once restrictions were re-imposed, which may suggest limited effectiveness. The unsampled sampling scheme again showed the most divergence from all other sampling schemes (Figure 6F) and the highest divergence from the case data estimate (Figure 6E) with the uniform and proportional sampling schemes showing the most similarity. As such, applying no sampling strategy/opportunistic sampling leads to, from the perspective of comparing to *EpiFilter*, the most biased estimates.

![Figure 6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F6.medium.gif)

[Figure 6:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F6)

Figure 6: 
*R**t* estimated from both the BDSKY and *EpiFilter* methods forAmazonas, Brazil. Titles indicate the sampling scheme used in panels A-D. The light-shaded area represents the 95% highest posterior density interval. The solid line represents the mean *R**t* estimate with *EpiFilter* in green and BDSKY in blue. We refer to Figure 2 for a brief description of key events, including 5 which corresponds to the second lockdown. Event “a’’ corresponds to the suspension of commercial activities in Manaus; event “b” corresponds to the resumption of commercial activities in Manaus72. The Jensen Shannon Distance (JSD) is given in panel E and ranks the sampling strategies based on how similar the BDSKY estimates under those strategies are to those derived from *EpiFilter* (smaller values are better). Panel F provides the pairwise JSD between the BDSKY estimates under different sampling strategies, showing often appreciable difference among strategies.

Based on the proportional sampling scheme, which had the lowest JSD (Figure 7E) we infer a steady decline in *r**t* which matches the pattern seen with the *R**t* value (Figure 7). The initial *r**t* implied a 0.250 d-1. Subsequently, the *r**t* falls over the study period. *r**t* falls below 0 after the re-imposition of NPIs declining at -0.030 d-1 by the end of the study period. There is no evidence of any noticeable declines in *r**t* when interventions were introduced indicating that they might not have significantly impacted the growth rate of P.1/gamma. The unsampled sampling scheme was again most divergent from other sampling schemes as well as from estimates derived from *EpiFilter* with the uniform and reciprocal-proportional being most similar.

![Figure 7:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F7.medium.gif)

[Figure 7:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F7)

Figure 7: 
*r**t* estimated from both the *Skygrowth* and *EpiFilter* methods for Amazonas, Brazil. Titles indicate the sampling scheme used in panels A-D. The light-shaded area represents the 95% highest posterior density interval. The solid line represents the mean *r**t* estimate with *EpiFilter* in orange and *Skygrowth* in blue. We refer to Figure 2 for a brief description of key events, including 5 which corresponds to the second lockdown. Event “a’’ corresponds to the suspension of commercial activities in Manaus; event “b” corresponds to the resumption of commercial activities in Manaus72. The Jensen Shannon Distance (JSD) is given in panel E and ranks the sampling strategies based on how similar the *Skygrowth* estimates under those strategies are to those derived from *EpiFilter* (smaller values are better). Panel F provides the pairwise JSD between the BDSKY estimates under different sampling strategies, showing often appreciable difference among strategies.

## Discussion

In this study, we applied phylodynamic methods to available SARS-CoV-2 sequences from Hong Kong and the Amazonas state of Brazil to infer their key epidemiological parameters and to compare the impact that various sampling strategies have on the phylodynamic reconstruction of these parameters.

We estimated the basic reproductive number of SARS-CoV-2 in Hong Kong to be 2.17 (95% CI = 1.43-2.83). This supports previous estimates of the initial *R*** in Hong Kong68,74 which estimated *R*** to be 2.23 (95% CI = 1.47-3.42). For the Amazonas state in Brazil, we estimated the *R*** to be 3.67 (95% CI = 2.83 – 4.48). Whilst the population of Amazonas State may not be fully susceptible to P.1/Gamma14,82, this should not affect the comparison among sampling schemes. We found that *R*** is robust to changes in sampling schemes (Figure 3A and C).

For the Hong Kong dataset, the proportional sampling scheme was superior to all other sampling schemes in estimating *R**t*. It successfully predicted the initial super-critical *R**t*, its decline in response to rapid NPIs, and subsequent increase and decline during the second wave of infections (Figure 4B). This was in comparison to the uniform sampling scheme, which provided the worst (largest) JSD (Figure 4D) and an *R**t* estimate that was largely insensitive to NPIs. The proportional sampling scheme, alongside the uniform sampling scheme, best estimated *r**t* (Figure 5B and C). In contrast, for the Amazonas dataset, the uniform sampling scheme best estimated the *R**t* and *r**t* (Figure 6C) whilst the proportional sampling scheme best captured *r**t* (Figure 7C). It captured both its initial super-critical *R**t* and high *r**t* alongside their subsequent decline.

We found that estimates from all sampling schemes were distinct from those obtained using the unsampled data and that on some instances the sampling schemes were also appreciably different from one another (see panel F in Figures 4-7) with the uniform and reciprocal-proportional sampling strategies being most similar. This highlights how different sampling schemes can produce significantly differing estimates of epidemiological parameters and underscores the need for considering sampling and its potential impact on estimations.

Our *R**t* estimates are consistent with previous estimates of Gamma VOC’s transmissibility in Amazonas state14. This contrasted with the unsampled data in which the *r**t* increased at the end of the period (Figure 7A). This highlights that unlike *R***, both *R**t* and *r**t* are sensitive to changes in sampling and that even related epidemiological parameters like *R**t* and *r**t* may require different sampling strategies to optimise inferences.

Molecular clock dating of the Hong Kong and Amazonas dataset has revealed that the date of origin is relatively robust to changes in sampling schemes. For Hong Kong, SARS-CoV-2 likely emerged in mid-December 2019 around 5 weeks before the first reported case on the 22nd of January 202068. The Amazonas dataset revealed that the date of the common ancestor of the P.1 lineage emerged around late October 2020 to early November, around 5 weeks before the first reported case on the 6th of December14, with all BCI’s overlapping for each sampling strategy. Like the molecular clock dating, we found that the molecular clock rate was robust to changes in sampling strategies in both datasets with all sampling strategies having overlapped BCI’s (Supplementary Table 2 and Supplementary Figure 5). For the Hong Kong dataset, its clock rate is comparable to early estimations of the mean substitution rate per site per year of SARS-CoV-213. However, the clock rate estimated for the Brazilian dataset is lower than the initial 8.00×10−4 s/s/y which is used in investigating SARS-CoV-275 and that has been used in previous analyses of Gamma VOC76. This initial estimation of evolutionary rate was estimated from genomic data taken over a short time span at the beginning of the pandemic introducing a time dependency bias. By using a more appropriate clock rate it can improve tree height and rooting resulting in more robust parameter estimations77.

Treating sampling times as uninformative has been shown to be inferior to including them as dependent on effective population size and other parameters by several previous studies30,31,34,78. Whilst these studies did not consider the estimation of epidemiological parameters, they highlight the potential of systematic biases being introduced into the phylodynamic reconstruction by not using a sampling scheme or by assuming an incorrect model for how sampling schemes introduce information. This was supported by our results as phylodynamic inferences with no sampling strategy applied had the poorest overall performance for both Hong Kong and the Amazonas state. This implies that sampling design choices can significantly impact phylodynamic reconstruction, and that exploration of sampling strategies is needed to obtain the most reliable estimates of key epidemiological parameters.

While our results provide rigorous insight into the dynamics of SARS-CoV-2 and the impact of sampling strategies in the Amazonas state and Hong Kong, there are limitations. The *Skygrowth* and BDSKY models do not explicitly consider imports into their respective regions. This is particularly relevant for Hong Kong as most initial sequences from the region were sequenced from importation events79 which can introduce error into parameter estimation80. However, as the epidemic expanded, more infections were attributable to autochthonous transmission79, and the risk of error introduced by importation events decreased. Moreover, while sampling strategies can account for temporal variations in genomic sampling fractions there is currently no way to account for non-random sampling approaches in either the BDSKY or *Skygrowth* models81. It is unclear how network-based sampling may affect parameter estimates obtained through these models82 presenting a key challenge in molecular and genetic epidemiology. Spatial heterogeneities were also not explored within this work. This represents the next key step in understanding the impact of sampling as spatial sampling schemes would allow the reconstruction of the dispersal dynamics and estimation of epidemic overdispersion (*k*), a key epidemiological parameter.

Finally, we compared our phylodynamic estimates against epidemiological inferences derived from incident case data from Hong Kong and Amazonas state, two settings with very different diagnostic capacity. While Hong Kong has high quality case data with a high testing rate69, there is a large underreporting of SARS-CoV-2 cases in the Amazonas state72,83. Future epidemiological modelling work is needed to compare parameter estimates obtained from case data, death data and excess death data across different settings. This will improve the benchmarks we use to compare sequence-based estimates against.

This work has highlighted the impact and importance that applying temporal sampling strategies can have on phylodynamic reconstruction. Whilst more genomic datasets from a variety of countries and regions with different sampling intensities and proportions are needed to create a more generalisable sampling framework and to dissect any potential cofounders, this study has demonstrated that genomic datasets that commonly feature opportunistic sampling (i.e., there is no deliberate strategy design) can introduce significant uncertainty and biases in the estimation of epidemiological parameters. This finding signifies the need for more targeted attempts at performing genomic surveillance and epidemic analyses particularly in resource-poor settings with limited genomic capability.

## Data Availability

All data produced in the present study are available upon reasonable request to the authors 

## Role of the Funding Sources

N.R.F. acknowledges support from Wellcome Trust and Royal Society Sir Henry Dale Fellowship (204311/Z/16/Z), Bill and Melinda Gates Foundation (INV-034540) and Medical Research Council-Sao Paulo Research Foundation (FAPESP) CADDE partnership award (MR/S0195/1 and FAPESP 18/14389-0) ([https://caddecentre.org](https://caddecentre.org)). K.V.P. acknowledges support from grant reference MR/R015600/1, jointly funded by the UK Medical Research Council (MRC) and the UK Department for International Development (DFID) and from the NIHR Health Protection Research Unit in Behavioural Science and Evaluation at University of Bristol. R.P.D.I acknowledges support from European Union Horizon 2020 project MOOD (#874850).

## CRediT authorship contribution statement

R.P.D.I, K.V.P and N.R.F conceived and designed the study, R.P.D.I wrote and performed the analyses. R.P.D.I wrote the manuscript which was edited and supervised by K.V.P and N.R.F. All authors have contributed to and approved the manuscript for submission.

## Supplementary Figures and Tables

![Supplementary Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F8.medium.gif)

[Supplementary Figure 1:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F8)

Supplementary Figure 1: 
The proportion of P.1 sequences compared to non-P.1 sequences from Amazonas State, Brazil found on GISaid (Shu and McCauley, 2017).

![Supplementary Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F9.medium.gif)

[Supplementary Figure 2:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F9)

Supplementary Figure 2: 
Root-to-tip genetic distances to sample collection dates for the SARS-CoV-2 genome datasets used in this study: A-D represents Hong Kong and E-H represent Amazonas State. Plots are based on the maximum likelihood trees rooted by maximising R2. The linear regression trend lines are shown to data points, corresponding to the genome sequences (represented with black dots).

![Supplementary Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F10.medium.gif)

[Supplementary Figure 3:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F10)

Supplementary Figure 3: 
Number of sequences for each week and sampling scheme for Hong Kong dataset.

![Supplementary Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F11.medium.gif)

[Supplementary Figure 4:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F11)

Supplementary Figure 4: 
Number of sequences for each week and sampling scheme for Amazonas dataset.

View this table:
[Supplementary Table 1:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/T4)

Supplementary Table 1: 
TMRCA and mean substitution rate both with 95% BCI for each sampling strategy for Hong Kong and Amazonas datasets alongside the Jensen-Shannon distance. Full posterior distribution of the TMRCA and substitution rates obtained under the different sampling strategies can be found in Figure 3B and D and Supplementary Figure 5.

![Supplementary Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/07/22/2022.02.04.22270165/F12.medium.gif)

[Supplementary Figure 5:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/F12)

Supplementary Figure 5: 
Mean substitution rate (s/s/y) for Hong Kong and Brazil. Figure 1A represents Hong Kong with Figure 1B representing the Amazonas. The central line represents the posterior mean and with intervals representing 95% Highest Posterior Density Interval

View this table:
[Supplementary Table 2:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/T5)

Supplementary Table 2: 
Accession ID of each Hong Kong sequence for each sampling strategy used within this study

View this table:
[Supplementary Table 3:](http://medrxiv.org/content/early/2022/07/22/2022.02.04.22270165/T6)

Supplementary Table 3: 
Accession ID of each Amazonas State, Brazil sequence for each sampling strategy used within this study

## Footnotes

*   5 Jointly supervised this work

*   We have performed additional analyses comparing how similar distributions are between sampling strategies for both the time-varying effective reproduction numbers (Rt) and rates of epidemic growth (rt). We have also moved supplementary figure 1 to the main text, fixed typing errors and clarified salient features of the manuscript to improve understanding.

*   Received February 4, 2022.
*   Revision received July 21, 2022.
*   Accepted July 22, 2022.


*   © 2022, Posted by Cold Spring Harbor Laboratory

This pre-print is available under a Creative Commons License (Attribution 4.0 International), CC BY 4.0, as described at [http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)

## References

1.  1.Gorbalenya, A. E. et al. The species Severe acute respiratory syndrome-related coronavirus: classifying 2019-nCoV and naming it SARS-CoV-2. Nature Microbiology 5, 536–544 (2020).
    
    

2.  2.Zhu, N. et al. A Novel Coronavirus from Patients with Pneumonia in China, 2019. New England Journal of Medicine 382, 727–733 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMoa2001017&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

3.  3.World Health Organisation. Public Health Emergency of International Concern (PHEIC). (2020).
    
    

4.  4.Verity, R. et al. Estimates of the severity of coronavirus disease 2019: a model-based analysis. The Lancet. Infectious diseases 20, 669–677 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1473-3099(20)30243-7&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32240634&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

5.  5.World Health Organisation. Coronavirus disease (COVID-19) Weekly Epidemiological Update and Weekly Operational Update. [https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports](https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports) (2022).
    
    

6.  6.European Centre for Disease Prevention and Control. Guidelines for the implementation of non-pharmaceutical interventions against COVID-19 Key messages General considerations on NPI to control COVID-19. (2020).
    
    

7.  7.Anderson Vegari, C., Baggaley, R., Hollingsworth, T. D. D. & Maddren, R. The Royal Society SET-C Reports. Reproduction number (R) and growth rate (r) of the COVID-19 epidemic in the UK: methods of estimation, data sources, causes of heterogeneity, and use as a guide in policy formulation [report unpublished]. The Royal Society 1–86 (2020).
    
    

8.  8.UK Health Security Agency. The R value and growth rate. [https://www.gov.uk/guidance/the-r-value-and-growth-rate](https://www.gov.uk/guidance/the-r-value-and-growth-rate) (2022).
    
    

9.  9.Parag, K.V., Thompson, R.N. & Donnelly, C.A. Are epidemic growth rates more informative than reproduction numbers?. J R Stat Soc Series A, 1–11 (2022).
    
    

10. 10.Dushoff, J. & Park, S. W. Speed and strength of an epidemic intervention. Proceedings of the Royal Society B: Biological Sciences 288, 20201556 (2021).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1098/rspb.2020.1556&link_type=DOI) 

11. 11.World Health Organisation. Genomic sequencing of SARS-CoV-2 A guide to implementation for maximum impact on public health. (2021).
    
    

12. 12.Jombart, T. et al. Bayesian Reconstruction of Disease Outbreaks by Combining Epidemiologic and Genomic Data. PLOS Computational Biology 10, e1003457. (2014).
    
    

13. 13.Duchene, S. et al. Temporal signal and the phylodynamic threshold of SARS-CoV-2. Virus Evolution 6, (2020).
    
    

14. 14.Faria, N. R. et al. Genomics and epidemiology of the P.1 SARS-CoV-2 lineage in Manaus, Brazil. Science 372, 815 LP –821 (2021).
    
    

15. 15.Nadeau, S. A., Vaughan, T. G., Scire, J., Huisman, J. S. & Stadler, T. The origin and early spread of SARS-CoV-2 in Europe. Proceedings of the National Academy of Sciences 118, e2012008118 (2021).
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxNzoiMTE4LzkvZTIwMTIwMDgxMTgiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMi8wNy8yMi8yMDIyLjAyLjA0LjIyMjcwMTY1LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 

16. 16.Romano, C. M. & Melo, F. L. Genomic surveillance of SARS-CoV-2: A race against time. The Lancet Regional Health - Americas , 100029 (2021).
    
    

17. 17.Volz, E. et al. Evaluating the Effects of SARS-CoV-2 Spike Mutation D614G on Transmissibility and Pathogenicity. Cell 184, 64–75.e11 (2021).
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

18. 18.Dudas, G. et al. Virus genomes reveal factors that spread and sustained the Ebola epidemic. Nature (2017) doi:10.1038/nature22040.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature22040&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28405027&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

19. 19.Faria, N. R. et al. Establishment and cryptic transmission of Zika virus in Brazil and the Americas. Nature 546, 406–410 (2017).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature22401&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28538727&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

20. 20.Grubaugh, N. D. et al. Genomic epidemiology reveals multiple introductions of Zika virus into the United States. Nature 546, 401–405 (2017).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature22400&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28538723&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

21. 21.Harvey, W. T. et al. SARS-CoV-2 variants, spike mutations and immune escape. Nature Reviews Microbiology 19, 409–424 (2021).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/S41579-021-00573-0&link_type=DOI) 

22. 22.Rambaut, A. et al. A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology. Nature Microbiology 5, 1403–1407 (2020).
    
    

23. 23.Shu, Y. & McCauley, J. GISAID: Global initiative on sharing all influenza data - from vision to reality. Euro surveillance : bulletin Europeen sur les maladies transmissibles = European communicable disease bulletin 22, 30494 (2017).
    
    

24. 24.Tsang, T. K. et al. Effect of changing case definitions for COVID-19 on the epidemic curve and transmission parameters in mainland China: a modelling study. The Lancet. Public health 5, e289–e296 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S2468-2667(20)30089-X&link_type=DOI) 

25. 25.de Souza, W. M. et al. Epidemiological and clinical characteristics of the COVID-19 epidemic in Brazil. Nature Human Behaviour 4, 856–865 (2020).
    
    

26. 26.Dolan, P. T., Whitfield, Z. J. & Andino, R. Mapping the Evolutionary Potential of RNA Viruses. Cell Host and Microbe 23, 435–446 (2018).
    
    

27. 27.Drummond, A. J., Rambaut, A., Shapiro, B. & Pybus, O. G. Bayesian Coalescent Inference of Past Population Dynamics from Molecular Sequences. Molecular Biology and Evolution 22, 1185–1192 (2005).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/molbev/msi103&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15703244&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000228700200004&link_type=ISI) 

28. 28.Gill, M. S. et al. Improving Bayesian population dynamics inference: a coalescent-based model for multiple loci. Molecular biology and evolution 30, 713–724 (2013).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/molbev/mss265&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23180580&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000314891800018&link_type=ISI) 

29. 29.Stadler, T., Kühnert, D., Bonhoeffer, S. & Drummond, A. J. Birth–death skyline plot reveals temporal changes of epidemic spread in HIV and hepatitis C virus (HCV). Proceedings of the National Academy of Sciences 110, 228LP–233 (2013).
    
    

30. 30.Hall, M. D., Woolhouse, M. E. J. & Rambaut, A. The effects of sampling strategy on the quality of reconstruction of viral population dynamics using Bayesian skyline family coalescent methods: A simulation study. Virus evolution 2, vew003–vew003 (2016).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/ve/vew003&link_type=DOI) 

31. 31.Parag, K. V, du Plessis, L. & Pybus, O. G. Jointly Inferring the Dynamics of Population Size and Sampling Intensity from Molecular Sequences. Molecular Biology and Evolution 37, 2414–2429 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/molbev/msaa016&link_type=DOI) 

32. 32.Stack, J. C., Welch, J. D., Ferrari, M. J., Shapiro, B. U. & Grenfell, B. T. Protocols for sampling viral sequences to study epidemic dynamics. Journal of the Royal Society, Interface 7, 1119–1127 (2010).
    
    

33. 33.de Silva, E., Ferguson, N. M. & Fraser, C. Inferring pandemic growth rates from sequence data. Journal of The Royal Society Interface 9, 1797–1808 (2012).
    
    

34. 34.Karcher, M. D., Palacios, J. A., Bedford, T., Suchard, M. A. & Minin, V. N. Quantifying and Mitigating the Effect of Preferential Sampling on Phylodynamic Inference. PLoS computational biology 12, e1004789–e1004789 (2016).
    
    

35. 35.Frost, S. D. W. et al. Eight challenges in phylodynamic inference. Epidemics 10, 88–92 (2015).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.epidem.2014.09.001&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25843391&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

36. 36.du Plessis, L. et al. Establishment and lineage dynamics of the SARS-CoV-2 epidemic in the UK. Science 371, 708–712 (2021).
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzNzEvNjUzMC83MDgiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMi8wNy8yMi8yMDIyLjAyLjA0LjIyMjcwMTY1LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 

37. 37.Hidano, A. & Gates, M. C. Assessing biases in phylodynamic inferences in the presence of super-spreaders. Veterinary Research 50, 74 (2019).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s13567-019-0692-5&link_type=DOI) 

38. 38.Gostic, K. M. et al. Practical considerations for measuring the effective reproductive number, Rt. PLOS Computational Biology 16, e1008409 (2020).
    
    

39. 39.Pullano, G. et al. Underdetection of cases of COVID-19 in France threatens epidemic control. Nature 590, 134–139 (2021).
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

40. 40.The World Bank. Population, total - Hong Kong SAR, China. [https://data.worldbank.org/indicator/SP.POP.TOTL?locations=HK](https://data.worldbank.org/indicator/SP.POP.TOTL?locations=HK) (2021).
    
    

41. 41.IBGE. Population Projections. [https://www.ibge.gov.br/en/statistics/social/population.html](https://www.ibge.gov.br/en/statistics/social/population.html) (2020).
    
    

42. 42.Byrne, A. W. et al. Inferred duration of infectious period of SARS-CoV-2: Rapid scoping review and analysis of available evidence for asymptomatic and symptomatic COVID-19 cases. BMJ Open 10, 1–16 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1136/bmjopen-2020-038832&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.ncbi.nlm.nih.gov/pubmed/32978187&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

43. 43.McAloon, C. et al. Incubation period of COVID-19: a rapid systematic review and meta-analysis of observational research. BMJ Open 10, e039652 (2020).
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoiYm1qb3BlbiI7czo1OiJyZXNpZCI7czoxMjoiMTAvOC9lMDM5NjUyIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjIvMDcvMjIvMjAyMi4wMi4wNC4yMjI3MDE2NS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 

44. 44.Parag, K. V. Improved estimation of time-varying reproduction numbers at low case incidence and between epidemic waves. PLOS Computational Biology 17, e1009347 (2021).
    
    

45. 45.Fraser, C. Estimating Individual and Household Reproduction Numbers in an Emerging Epidemic. PLoS ONE 2, e758 (2007).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pone.0000758&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17712406&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

46. 46.Cori, A., Ferguson, N. M., Fraser, C. & Cauchemez, S. A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics. American Journal of Epidemiology 178, 1505–1512 (2013).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/aje/kwt133&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24043437&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

47. 47.Wallinga, J. & Teunis, P. Different Epidemic Curves for Severe Acute Respiratory Syndrome Reveal Similar Impacts of Control Measures. American Journal of Epidemiology 160, 509–516 (2004).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/aje/kwh255&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15353409&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000223938000001&link_type=ISI) 

48. 48.Wallinga & Lipsitch. How generation intervals shape the relationship between growth rates and reproductive numbers. Proceedings of the Royal Society B: Biological Sciences 274, 599–604 (2007).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1098/rspb.2006.3754&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17476782&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000243354200019&link_type=ISI) 

49. 49.Rai, B., Shukla, A. & Dwivedi, L. K. Estimates of serial interval for COVID-19: A systematic review and meta-analysis. Clinical epidemiology and global health 9, 157–161 (2021).
    
    

50. 50.Prete, C. A. et al. Serial interval distribution of SARS-CoV-2 infection in Brazil. Journal of travel medicine 28, 1–3 (2021).
    
    

51. 51.Katoh, K., Misawa, K., Kuma, K. & Miyata, T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Research 30, 3059–3066 (2002).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/nar/gkf436&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12136088&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000177154300016&link_type=ISI) 

52. 52.Hadfield, J. et al. Nextstrain: real-time tracking of pathogen evolution. Bioinformatics 34, 4121–4123 (2018).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10/gdkbqx&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

53. 53.Minh, B. Q. et al. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. Molecular Biology and Evolution 37, 1530–1534 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/molbev/msaa015&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32556291&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

54. 54.Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. ModelFinder: fast model selection for accurate phylogenetic estimates. Nature Methods 14, 587–589 (2017).
    
    

55. 55.Tamura, K. & Nei, M. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Molecular Biology and Evolution 10, 512–526 (1993).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/oxfordjournals.molbev.a040023&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=8336541&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1993LD11400002&link_type=ISI) 

56. 56.Anisimova, M., Gil, M., Dufayard, J.-F., Dessimoz, C. & Gascuel, O. Survey of branch support methods demonstrates accuracy, power, and robustness of fast likelihood-based approximation schemes. Systematic biology 60, 685–699 (2011).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/sysbio/syr041&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21540409&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

57. 57.Rambaut, A., Lam, T. T., Max Carvalho, L. & Pybus, O. G. Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen). Virus evolution 2, vew007–vew007 (2016).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/ve/vew007&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27774300&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

58. 58.World Health Organisation. Guidance for surveillance of SARS-CoV-2 variants Interim guidance. (2021).
    
    

59. 59.Suchard, M. A. et al. Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus evolution 4, vey016–vey016 (2018).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/ve/vey016&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29942656&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

60. 60.Ayres, D. L. et al. BEAGLE 3: Improved Performance, Scaling, and Usability for a High-Performance Computing Library for Statistical Phylogenetics. Systematic Biology 68, 1052–1061 (2019).
    
    

61. 61.Hasegawa, M., Kishino, H. & Yano, T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution 22, 160–174 (1985).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/BF02101694&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=3934395&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1985AWB1600007&link_type=ISI) 

62. 62.Hill, V. & Baele, G. Bayesian Estimation of Past Population Dynamics in BEAST 1.10 Using the Skygrid Coalescent Model. Molecular Biology and Evolution 36, 2620–2628 (2019).
    
    

63. 63.Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Systematic biology 67, 901–904 (2018).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/sysbio/syy032&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29718447&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

64. 64.Bouckaert, R. et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLOS Computational Biology 15, e1006650 (2019).
    
    

65. 65.Volz, E. M. & Didelot, X. Modeling the Growth and Decline of Pathogen Effective Population Size Provides Insight into Epidemic Dynamics and Drivers of Antimicrobial Resistance. Systematic Biology 67, 719–728 (2018).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/sysbio/syy007&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29432602&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

66. 66.Lin, J. Divergence measures based on the Shannon entropy. IEEE Transactions on Information Theory 37, 145–151 (1991).
    
    

67. 67.Kullback, S. & Leibler, R. A. On Information and Sufficiency. The Annals of Mathematical Statistics 22, 79–86 (1951).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.transproceed.2017.12.043&link_type=DOI) 

68. 68.Cowling, B. J. et al. Impact assessment of non-pharmaceutical interventions against coronavirus disease 2019 and influenza in Hong Kong: an observational study. The Lancet Public Health 5, e279–e288 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/s2468-2667(20)30090-6&link_type=DOI) 

69. 69.Wu, P. et al. Suppressing COVID-19 Transmission in Hong Kong: An Observational Study of the First Four Months. SSRN (2020) doi:10.21203/rs.3.rs-34047/v1.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.21203/rs.3.rs-34047/v1&link_type=DOI) 

70. 70.Nascimento, V. A. do et al. Genomic and phylogenetic characterisation of an imported case of SARS-CoV-2 in Amazonas State, Brazil. Memórias do Instituto Oswaldo Cruz 115, (2020).
    
    

71. 71.Faria, N. R. et al. Genomic characterisation of an emergent SARS-CoV-2 lineage in Manaus: preliminary findings. Virological [https://virological.org/t/genomic-characterisation-of-an-emergent-sars-cov-2-lineage-in-manaus-preliminary-findings/586](https://virological.org/t/genomic-characterisation-of-an-emergent-sars-cov-2-lineage-in-manaus-preliminary-findings/586) (2021).
    
    

72. 72.Sabino, E. C. et al. Resurgence of COVID-19 in Manaus, Brazil, despite high seroprevalence. Lancet (London, England) 397, 452–455 (2021).
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

73. 73.Drummond, A. J., Pybus, O. G., Rambaut, A., Forsberg, R. & Rodrigo, A. G. Measurably evolving populations. Trends in Ecology & Evolution 18, 481–488 (2003).
    
    

74. 74.Zhao, S. et al. Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak. International journal of infectious diseases : IJID : official publication of the International Society for Infectious Diseases 92, 214–217 (2020).
    
    

75. 75.Andersen, K. G., Rambaut, A., Lipkin, W. I., Holmes, E. C. & Garry, R. F. The proximal origin of SARS-CoV-2. Nature Medicine (2020) doi:10.1038/s41591-020-0820-9.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-020-0820-9&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32284615&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

76. 76.Naveca, F. G. et al. COVID-19 in Amazonas, Brazil, was driven by the persistence of endemic lineages and P.1 emergence. Nature Medicine 27, 1230–1238 (2021).
    
    

77. 77.Boskova, V., Stadler, T. & Magnus, C. The influence of phylodynamic model specifications on parameter estimates of the Zika virus epidemic. Virus Evolution 4, (2018).
    
    

78. 78.Liu, Q. et al. Population Genetics of SARS-CoV-2: Disentangling Effects of Sampling Bias and Infection Clusters. Genomics, Proteomics & Bioinformatics 18, 640–647 (2020).
    
    

79. 79.Adam, D. C. et al. Clustering and superspreading potential of SARS-CoV-2 infections in Hong Kong. Nature Medicine 26, 1714–1719 (2020).
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

80. 80.Parag, K. v, Cowling, B. J. & Donnelly, C. A. Deciphering early-warning signals of SARS-CoV-2 elimination and resurgence from limited data at multiple scales. Journal of The Royal Society Interface 18, 20210569 (2022).
    
    

81. 81.Vasylyeva, T. I. et al. Phylodynamics helps to evaluate the impact of an HIV prevention intervention. Viruses 12, 1–15 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3390/v12040404&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F07%2F22%2F2022.02.04.22270165.atom) 

82. 82.Volz, E. M., Koelle, K. & Bedford, T. Viral phylodynamics. PLoS computational biology 9, e1002947–e1002947 (2013).
    
    

83. 83.Buss, L. et al. Three-quarters attack rate of SARS-CoV-2 in the Brazilian Amazon during a largely unmitigated epidemic. Science 371, 288–292 (2021).
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzNzEvNjUyNi8yODgiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMi8wNy8yMi8yMDIyLjAyLjA0LjIyMjcwMTY1LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==)

 [1]: /embed/graphic-3.gif
 [2]: /embed/graphic-4.gif
 [3]: /embed/graphic-7.gif
 [4]: /embed/graphic-8.gif