Analyses of Early ZIKV Genomes are Consistent with Viral Spread from Northeast Brazil to The Americas

1. Introduction

Zika virus (ZIKV) is an emerging arthropod-borne virus that causes severe neurotropic diseases, which include Guillain-Barré and Congenital Zika Syndrome [1]. The virus belongs to the Flavivirus genus of the Flaviviridae family. It has a positive-sense single-stranded RNA genome of approximately 11 kb [2,3] and a substitution rate of 7.55×10^-4 to 1.66×10^-3 substitutions per site per year [4]. ZIKV was first isolated in 1947 from a sentinel monkey in the Zika Forest, Uganda. ZIKV was described as causing sporadic infections for half a century. However, in 2007 and 2013, large outbreaks were reported in Micronesia [5] and French Polynesia [6]. In March 2015, ZIKV was identified in Brazil [7,8] and further spread to more than 40 countries in the Americas [9,10]. In February 2016, the World Health Organization (WHO) declared the ZIKV outbreak a Public Health Emergency of International Concern (PHEIC) due to its association with microcephaly and neurological diseases as Guillain-Barré syndrome [11,12]. This PHEIC declaration was lifted in November 2016 [13].

The revolution in molecular epidemiological surveillance with the advent of next-generation sequencing and advanced phylogenetic analyses allows to uncover unobserved transmission dynamics and helps to explain the introduction, spread, and evolution of infectious diseases [11,14]. Molecular dating analyses based on ZIKV sequences from patients identified during the Brazilian outbreak revealed that this virus possibly entered Brazil in the second half of 2013 [14]. Additionally, Lednicky et al. (2016) [15] and Campos et al. (2018) [16] described that ZIKV in Brazil possibly originated from Haiti. Both discoveries revealed, at the time of the event, a gap in our national surveillance system that hampered early detection and interventions to curb the spread of the virus in the country.

However, these findings need validation with new strategies of phylogenetic analyses and new sequences. This is because factors such as the limited number of ZIKV sequences compared to the proportion of cases reported during the outbreak in Brazil and the low sequence quality from the beginning of the epidemic (breadth of coverage < 70%) can directly affect the credibility of pathogen transmission dynamics reconstruction [17].

The low number of high-quality sequences available for ZIKV can be explained by several factors, including low viral load in ZIKV infection [18], low detection of ZIKV-infected patients due to a large number of asymptomatic cases (75% to 80%) or mild/unspecific symptoms [19], abrupt reduction in ZIKV cases [20] and the challenges faced by the Brazilian health system to perform real-time genomic surveillance. Retrospective sequencing of samples from the early stages of the outbreak in Brazil has the potential to increase the accuracy and precision of transmission dynamics reconstruction, which can help inform public health strategies for future outbreaks.

2. Materials and Methods

During the ZIKV outbreak, our group actively recruited patients with arboviral-like symptoms in Salvador (n=948) [21] and Campo Formoso (n=230) [22], Bahia. A total of 21 cases of viremic ZIKV infection were identified between May 2015 to June 2016 in Bahia, Brazil. An amplicon tiling multiplex approach described by Quick et al. (2017) [23] was applied in conjunction with the MinION platform (Oxford Nanopore Technologies) for sequencing. For further details on diagnosis, sequencing and phylogenetic analysis, we refer to the Supplementary Appendix.

3. Results and Discussion

The sequencing efforts yielded near full-length sequences (>80% reference coverage, Genbank No. KJ776791.1) for 14 ZIKV-positive samples, with an average sequencing depth ranging from 377x to 1,438x (Table 1). The clinical and social-demographic characteristics are listed in Table S1.

A comprehensive dataset was compiled with 14 new genomes and 206 curated near full-length ZIKV genomes (Table S2). We performed a time-calibrated discrete phylogeographic reconstruction (details and references in Supplementary Appendix) to explore the geographic spread and transmission of ZIKV in Bahia and Brazil, as shown in Figure 1.

The newly sequenced ZIKV genomes analyzed in this study are distributed over four distinct clades, containing sequences from diverse regions of Brazil (North, Northeast, and Southeast) as well as other countries in the world, thus providing valuable insights into the migration patterns of ZIKV.

Ten of the newly sequenced isolates (TRDP238, TRDP256, TRDP257, TRDP274, TRDP282, TRDP300, TRDP309, TRDP317, TRDP333, and ZK0110) clustered in a clade together with genomes from Brazil (Southeast and Northeast), Singapore and Haiti (Clade I, Figure 1B). The estimated time to the most recent common ancestor (tMRCA) was August 2013, with 95% Bayesian high posterior density (HPD) between January 2013 and February 2014. The Singapore isolate was sequenced from a traveler from São Paulo (Brazil) to Singapore in May 2016 and was described as the first confirmed case of ZIKV infection in this country [24]. The isolates from Haiti have been previously associated to isolates from Brazil by Lednicky et al. (2016) [15] and Campos et al. (2018) [16], which also suggested that the 2014 Haitian ZIKV strain led to the spread of ZIKV in Brazil. However, our phylogeographic reconstruction with the new isolates, reveals a paraphyletic clustering of viruses from Haiti with respect to Brazilian lineages, suggesting that ZIKV spread from Brazil to Haiti. In addition, this reconstruction suggests that an outbreak in Salvador was responsible for seeding other locations in Northeast and Southeast Brazil and later spreading to Haiti. Noteworthy, this had already been suggested by Massad et al. (2017) [25], but based only on considering the limited numbers of ZIKV infections in Haiti compared to those obtained from the Brazilian epidemic.

Two of the newly sequenced isolates (TRDP173 and ZK0152) clustered together in a clade with genomes isolated from Brazil (Southeast and North) and USA (Clade II, Figure 1C), with a tMRCA estimated around February 2014 (95% HPD: July 2013, September 2014). The USA isolate was reported in a study by Grubaugh et al. (2019) [10], which clustered exclusively with Brazilian sequences from the Southeast [10]. However, in our analyses, the USA sequence formed a monophyletic clade with the TRDP173 isolated from Bahia.

The TRDP252 isolate clustered together in a clade with sequences from the Brazil (Southeast, North and Northeast), and Dominican Republic (Clade III, Figure 1D), with a tMRCA estimated around October 2013 (95% HPD: April 2013, March 2014). This Dominican Republic isolate was reported in a study by Metsky et al. (2017) [26], which clustered with Brazilian sequences, including from the Southeast, similar to our results.

The last isolate (TRDP433), sampled in July 2015, clustered together with sequences from Northeast Brazil, with a tMRCA estimated around March 2013 (95% HPD: September 2013, August 2014). All the isolates in this last clade were from Salvador with sampling dates between May 2015 and January 2016 (Clade IV, Figure 1E).

By analyzing sequences from other countries with the addition of new sequences from the Northeast region, we provide further evidence supporting the central role of ZIKV migration from the Northeast to the Southeast of Brazil and its subsequent spread outside the country. Unfortunately, the long time between the MRCA and the sequence dates does not allow more confident conclusions about the migration routes of ZIKV.

Although ZIKV was detected almost simultaneously in Brazil in the Salvador and Natal metropolitan regions during the period of 2015-2016, before the present study, only 5 sequences with >80% reference coverage (Genbank No. KJ776791.1) originating from this region were publicly available. The poor genomic sampling from this crucial moment and place of the Brazilian ZIKV outbreak limits the power of phylogenetic reconstructions and could bias the interpretation of disease transmission dynamics and may result in incorrect interpretations regarding the origin of ZIKV introduction and its routes of spread in Brazil.

Of note, since 2016, several independent efforts have been made in genomic surveillance of ZIKV to cover most of the Brazilian territories over different time periods, such as the prominent ZiBRA project that produced more than 140 sequences (Zika in Brazil Real Time Analysis) [23,27,28,29,30,31]. However, these great efforts were unable to rescue enough samples from this specific period due to logistic limitations in storing samples for long periods in Central Public Health Laboratories (LACEN).

5. Conclusions

Here we provide more credibility in the patterns of ZIKV spread from Northeast Brazil to the Americas based on the generation of 14 extra ZIKV sequences with good coverage (>80%), representing ~75% of the new publicly available dataset of this period, increasing our understanding of the dynamics of arbovirus infections worldwide and lending support to future public health responses against this virus.