INTRODUCTION

The COVID19 pandemic is caused by SARS- CoV-2 which was first reported in December 2019 in Wuhan, China. As of February 2021, it has become a global public health problem, affecting more than 111 million people, and resulting in more than 2.4 million deaths ¹. Peru reported its first case on March 6, 2020. Ten days later, the government decreed a state of national emergency that included quarantine measures and a curfew ². However, transmission was widespread within the country. In addition, the number of positive cases by molecular tests increased to more than 364,000 and by rapid tests to more than 840,000 ³.

Epidemiological tools are of great importance during a pandemic, especially during the decision-making process, because they can provide a picture of the initial and current transmission dynamics. The use of these tools is supported by the early isolation of positive cases, contact tracing, and selection of appropriate community and individual mitigation measures ⁴. One of the most important estimates during the early stages of the pandemic is the basic reproductive number (R0) ^{5 , 6}, which represents the number of secondary cases expected from a primary case in a population in which all individuals are susceptible to the disease; it measures the potential speed of transmission during the initial phase of the pandemic. R0 was estimated at 2.3 (95% confidence interval [95% CI]: 2.0-2.5) for Lima ^{2 )} and between 2.36 and 5.24 for Peru ^{6 )} during March 2020.

Another important parameter is the effective reproductive number (Rt), which measures the expected number of new infections caused by an infectious individual during the course of the pandemic, when there are already non-susceptible individuals in the population because they were already infected or are immune to the disease ⁷. If Rt>1, the cases will increase exponentially; if Rt<1, the number of infected individuals will decrease ⁸. The Rt for Peru was estimated at 2.36 (95% CI: 2.11- 2.63) during March 2020 ⁷.

The availability of free data in public repositories encourages open science and rapid public health decision-making ⁹. For example, the first SARS-CoV-2 genome, published in January 2020, enabled the early development of RT-PCR molecular tests ¹⁰, which are used worldwide to confirm the presence of the virus. Genomic information is also used in phylodynamic analyses to obtain estimates of the transmission of emerging diseases, and in real-time molecular evolution processes, thus collaborating with epidemiological surveillance ¹¹.

Rt was calculated with models that are similar to the epidemiological SIR (susceptible, infected, recovered), using phylogenetic information to obtain transmission and recovery rates ¹². The first genome of a Peruvian case was sequenced by the Instituto Nacional de Salud (INS) in March 2020 ¹³. As of February 2021, there were 113 Peruvian case genomes from March and April 2020 deposited in the public database Global Initiative on Sharing All Influenza Data (GISAID).

The aim of our study was to use the genomic information available from Peru to estimate the Rt during the first two months of the epidemic and to determine the dynamics of initial transmission by integrating molecular analyses with epidemiological estimates.

THE STUDY

This was a descriptive study. The information was obtained from the GISAID database in February 2021 and corresponds to March and April 2020. The unit of analysis was each of the genomes available for that period. We used 113 genomes that met the following criteria: complete genome (>29,000 base pairs) and a percentage of ambiguous nucleotides of less than 1%. Accession numbers can be found in Supplementary Material 1. Of the samples analyzed, 59 were from Lima and 54 from other regions: 8 from Áncash, 7 from Arequipa, 2 from Cajamarca, 2 from Callao, 2 from Cusco, 1 from Huánuco, 5 from Ica, 7 from Junín, 3 from La Libertad, 12 from Lambayeque and 5 from Loreto. The unit of time for the analysis was a calendar day. Sequences were aligned with the MAFFT program ¹⁴ and the SARS-CoV-2 reference genome (GenBank identification code: NC_045512.2). The multiple sequence alignment was manually edited in the Geneious v.2020.2.3 program to remove the flanks corresponding to the untranslated regions (UTRs). The final alignment size resulted in 29,409 sites. This alignment can be found in Supplementary Material 2.

Estimation of Rt from genetic data is possible using the Birth-Death Skyline (BDSKY) model ¹⁵. The model is based on the birth-death process, in which individuals generate new disease transmission events at a rate λ (birth) and become noninfectious at a rate δ (death), due to host recovery, behavioral changes, or successful treatment. The effective reproductive number estimated from the genomes is calculated by dividing both rates.

The parameters needed to model these rates were estimated from phylogenies generated by Bayesian analysis in the BEAST2 program ¹⁶. This program uses Markov chain Monte Carlo (MCMC) and estimates the parameters of the BDSKY model at each generation of the MCMC. Ten MCMCs were run, with 50 million generations each. Each chain was sampled every 50,000 generations, the first 25% of the results from each chain were discarded (burnin). The results of the 10 MCMCs were combined into a single file (Supplementary Material 3).

From an initial (a priori) distribution, and a significant number of generations, the values tend to appear in a stationary (posterior) distribution. When this happens, the parameters are said to have achieved “convergence”. To accept the convergence of MCMC, 200 or more is the needed estimated sample size (ESS) in the posterior distribution.

A nucleotide substitution rate of 22,905 substitutions per year was used to calibrate the phylogenies. This rate was obtained from the NextStrain genomic surveillance portal ¹⁷. The result of the division between the rate and the size of the reference genome (29,903 nucleotides) was 0.00077 substitutions/site/year. Similar mutation rates have been used in other studies on introductions and early SARS-CoV-2 transmission ^{18 , 19}.

The a priori parameters are shown in Table 1. We carried out the following: 1) Calculation of the rate of becoming non-infectious (becomeUninfectiousRate), expressed in units per year, and resulting from dividing the days of the year by the average time of days during which an individual is infectious (365/15 ~25 units). We chose 25 as the initial mean of a normal distribution. 2) The common ancestor of the Peruvian samples was estimated between late December 2019 and January 2020; therefore, we chose a normal distribution with mean 0.3 (~110 days ago) as the initial distribution. 3) The Rt was estimated using a lognormal distribution with mean 1.65 and standard deviation 1 as the parameter, which implies 0.1<Rt<7.1. The Rt was calculated for each week, approximately.

The information used in this study comes from the public database GISAID. None of the data analyzed had any personal identification. This research is aligned with the need for information to address a health emergency.

FINDINGS

The stationary distribution of MCMCs, with convergence (ESS > 200), was achieved by estimating the Rt (Table 2, Supplementary Material 3). We calculated the rate of becoming non-infectious at 39.7, HPD 95%: 31.76-47.62 (HPD: highest probability density interval). This value corresponds to an estimated 365/39.7 ~9 days, HPD 95%: 7.66 to 11.5 days, in which the individual would remain infectious. The Rt was greater than 1 day before the first report of SARS-CoV-2 in Peru (Figure 1, Supplementary Material 4). Rt decreased during the second week of March and reached an Rt<1 value a few days before the start of quarantine. The lowest Rt values were found between March 19 and 25 (Supplementary Material 4), after the quarantine was decreed. Rt fluctuated close to 1 during the following weeks until rising again in the second half of April 2020.

Figure 1 Estimated values of Rt by the Birth-Death Skyline model. Uncertainty is represented by a distribution with higher probability density intervals at 95% (HPD 95%) (shaded area). The Rt value decreases during the first two weeks of March and increases moderately during the last week of April 2020. 

DISCUSSION

During the first ten days of pandemic in Peru, the Rt reported was 2.36-5.24 ⁷. Another report calculated an Rt of 4.4 for March 2020 and a reduction to 3.2 during the first two weeks of quarantine ⁴. Our research shows concordance with these estimates and with non-peer-reviewed reports ⁸. The higher initial values could be due to initial uncertainty in the calculation of Rt during the first few days ⁸.

Control of the epidemic can be achieved by reaching a value of Rt<1. According to our results, the quarantine allowed Peru to maintain this value for at least one week after the quarantine was decreed. However, the decrease in Rt before the beginning of the quarantine could be explained by the measures taken by the government the previous week, such as suspending classes and flights from Europe and Asia, as well as prohibiting meetings of more than 300 people ².

Rt<1 was maintained during the first week of quarantine probably due to the strict conditions applied during the first days. The subsequent increase in Rt after that week may have been caused by recurrent crowding in public places such as markets or banks. For example, the failure of the sex-specific mobilization restriction (April 3-10, 2020), which caused crowding on days intended for females.

Although early strict mitigation measures such as confinement were taken, it did not prevent people from moving out of their homes. Confinement has worked in countries with formal labor markets, but in Peru, the suppression of informal economic activities — which reaches more than 72% — is an unsustainable measure in the context of our social dynamics ²⁰, that may have contributed to people leaving their homes and enabled the transmission of the virus.

Regarding the limitations of our study, the number of genomes that were sampled is small compared to the total cases reported in Peru during March and April 2020. However, other studies have used similar methodologies to estimate the Rt, assuming a sampling frequency of 10^-3 and 10^-5 (1/1,000 to 1/100,000) genomes per reported cases ^{18 , 19}. At the end of April 2020, there were more than 36,000 cases, so the sampling frequency in Peru would be 113/36,000 ~10-³. Even if this official figure was lower than the actual existing cases, the value would still be within the range used in other studies.

On the other hand, geographic representation should also be improved to determine the value of the Rt at the regional level with a greater detail. This is another limitation of our study and is related to the public data available in the GISAID database. Therefore, we suggest increasing the geographic diversity of the reported samples.

In conclusion, by using genomic information, it is evident that the Rt decreased considerably during the first half of March, reached its lowest value the week after the beginning of quarantine, and then increased moderately since the second half of April.

It is clear that more sequencing initiatives are needed to strengthen genomic surveillance in our country. The power of this analysis facilitates follow-up of new mutations, foreign variants, and local variants and provides data on the genetic diversity of the virus, its molecular evolution and its phylogeographic transmission patterns in our country. In this study, the genomic information obtained was analyzed with an epidemiological approach to estimate the speed of virus transmission during the early stage of the pandemic.

Since genomic data complement epidemiological studies, their integration is necessary for a future national epidemiological surveillance system that will allow public health decisions to be made based on evidence provided by open and multidisciplinary science. It is also necessary to promote the training of skilled human resources who can understand and integrate epidemiological methods with bioinformatics tools.