Epidemiological dynamics of an urban Dengue 4 outbreak in São Paulo, Brazil

Study site

Guarujá (23°59′37″S 46°15′23″W) is a coastal city in Santo Amaro Island, situated at the shore of the state of São Paulo, Brazil (Fig. 1A). The city is embedded in a tropical rain forest. It has a tropical humid climate that is characterized by having high average air temperature and rainfall. The average annual temperature is 24.7 °C (Min 18 °C | Max 31.3 °C) and the annual rainfall is 3,400 mm; February is the wettest month (average rainfall of 413 mm) and August is the driest one (average rainfall of 156 mm). The city main economic sources are seasonal tourism and port related activities. The estimated population in 2013 was 306,683 and the human population density was around 2,000 inhabitants per km². Official dengue figures by the Epidemiological Surveillance Center of the State (CVE) date back to 1997 and sum 24,000 total cases up to 2012. There is no detailed information concerning the previous exposure to distinct DENV serotypes and the municipality relies on the surveillance by the Adolfo Lutz Institute (The Central Public Health laboratory from The State of São Paulo).

Sample collection

In late 2012, our group at University of São Paulo joined efforts with the Guarujá Municipality’s office of epidemiological surveillance and with a local clinical laboratory analysis center (Itapema) to map the incidence of dengue in the city and obtain viral genetic data. Both institutions contributed with the collection of samples citywide and by performing preliminary immunochromatographic diagnostic tests for dengue (Kassim, 2011).

Patients of any age with symptoms and signs of dengue disease that were examined in Primary Health or Emergency Care Units and tested positive for IgM and/or NS1 were considered for the study. Symptoms included fever, frontal or retro-orbital headache, severe pain (muscles, bones, legs, joints, lower back or abdominal), nausea, vomiting, taste disturbance and anorexia. Signs included high fever (usually between 38.5 and 41 °C) persisting for 24 h, rash, hemorrhagic manifestations, hypotension and narrow pulse pressure.

Ethics statement

The Ethical Review Board of the Biomedical Science Institute at University of São Paulo approved this study (Statement 933/CEP). All adult subjects provided an informed written consent, and a parent or guardian of any child participant provided the written informed consent on their behalf.

Molecular testing

Viral RNA was extracted from serum samples with the QIAmp viral RNA mini kit (Qiagen, Venlo, Limburg, Netherlands) and the complementary DNA was synthesized using the SuperScript® VILO™ cDNA Synthesis Kit (Life Technologies, Carlsbad, California, United States). The GoTaq® Green Master Mix (Promega, Madison, Wisconsin, United States) was used for PCR amplifications of the envelope gene using the primers of Bennett et al. (2003). The ExoSAP-IT reagent was used for PCR Product Cleanup (Affymetrix, Santa Clara, California, United States) and the sequencing reaction was performed using the BigDye® Terminator v3.1 Cycle Sequencing Kit (Life Technologies). Sequencing reaction products were purified using the BigDye XTerminator Purification Kit (Life Technologies) and sequenced on an ABI PRISM® 3130 Genetic Analyzer (Life Technologies). Contigs were assembled using the program Codon Code aligner (www.codoncode.com/aligner).

Genetic analysis

Sequences were aligned using Muscle 3.8.31 (Edgar, 2004a; Edgar, 2004b) followed by visual inspection and manual editing with Mesquite 2.75 (Maddison & Maddison, 2014). Polymorphisms were analyzed with DNASP 5 (Librado & Rozas, 2009). The Tajima D statistic test was used to evaluate deviations from the neutral expectation of molecular evolution (Tajima, 1989). The package HyPhy v2.2 was used to screen for recombination (SBP-Single Breakpoint Recombination and GARD-Genetic Algorithms for Recombination Detection) and for positive selection; the dengue strain H780090 isolated in Boa Vista, RR-Brazil (29 November 2010) was used as a reference for selection analysis. Both genealogy-based, codon-site models Single Likelihood Ancestor Counting (SLAC) and the Fixed Effects Likelihood (FEL) methods were used to estimate the non-synonymous (dN) and synonymous (dS) rates of substitution (Delport et al., 2010; Kosakovsky Pond & Frost, 2005; Pond, Frost & Muse, 2005).

Phylogenetic analysis

The sequences obtained in this study were combined with a DENV-4 database from a previous study (Villabona-Arenas & Zanotto, 2011) to identify the genotype. This was achieved using high-throughput clustering with the UCLUST algorithm in the package USEARCH (Edgar, 2010).

The JModeltest software was used for the statistical selection of the best-fit model of nucleotide substitution under the Akaike information criterion (Darriba et al., 2012; Guindon & Gascuel, 2003). Sequences were dated according to the day of sampling and used for phylogenetic reconstruction and the estimation of the rate of evolutionary change (μ) (subs/site/year) using Bayesian Inference (IB) in Beast v2.3.1 (Bouckaert et al., 2014); the tree prior was a Birth-Death with Serial Skyline Sampling (BDSKY) (Stadler et al., 2013) (see Table S1 for parameterization). A Bayesian maximum clade credibility (MCC) tree was inferred from a set of plausible trees sampled at the stationary phase of four independent Markov Chain Monte Carlo (MCMC) runs with 200 million generations each using a relaxed (uncorrelated lognormal) molecular clock (Drummond et al., 2006). The convergence of parameters was assessed using Tracer v1.6 program (http://tree.bio.ed.ac.uk/software/tracer/) until all parameters estimates showed effective sample size (ESS) values over 200.

Time-varying reproduction numbers

We used the approach of Cori et al. (2013) to estimate the instantaneous reproduction number using the R-package EpiEstim (Cori et al., 2013; Salje et al., 2012). The method requires incidence data and a serial interval distribution (a gamma distribution with shift 1) that describes the time between the onset of symptoms in a primary case and the onset of symptoms of secondary cases. We used the weekly notifications from the municipality’s office of epidemiological surveillance for the epidemiological year of 2012–2013 (daily data was not available) and a serial interval distribution that reflected previous estimates of the dengue incubation period (Table S1); the outcome of a censored Bayesian time-to-event model estimated the dengue intrinsic incubation period (the time between a human being infected and the onset of symptoms due to the infection) around six days (95% CI 3–10 days) and the best-fitting temperature-dependent extrinsic incubation period around 6.5 days (95% CI 2–15 days) at 30 °C (Chan & Johansson, 2012). The analysis took into account the uncertainty in the serial interval distribution by integrating over a range of means and standard deviations of the serial interval (the mean and standard deviation were allowed to vary according to truncated normal distributions; Table S1) (Cori et al., 2013). We used a gamma prior distribution (Table S1) for the reproduction number that includes previous estimates (between 1.33 and 11.6) (Halstead, 2008).

Spatiotemporal dispersion pattern

In order to find the most parsimonious set of rates explaining the diffusion process along the sampled trees for a geolocated dataset, we used the Bayesian stochastic search variable selection (BSSVS) approach as implemented in Beast v2.3.1 (Bouckaert et al., 2014). The method assumes exchange rates in a continuous-time Markov chain (CTMC) to be zero with some prior probability (Lemey et al., 2009) and performs ancestral reconstruction on a single character, which represents the location of the taxa. Sampled taxa are associated with locations and the ancestral states of the internal nodes in the sampled trees can be reconstructed from the taxon locations; because the sampled trees are in units of time we can reconstruct the diffusion over time by following character transformation over the branches (Drummond & Bouckaert, 2015). Locations were represented by discrete groups of adjacent neighborhoods (discrete phylogeography) and a Bayes factor (BF) test was run to identify the rates contributing to the migration path with the software Spread v1.0.4 (Bielejec et al., 2011). The number of neighborhoods was reduced to a maximum of 10 localities (chosen by vicinity a by number of samples reported) in order to diminish sample-size bias (disproportionate sampling can strongly bias phylogeographic analyses because over sampled populations will more likely to be inferred as source populations).

Sampling

The year 2013 coincided with a steep rise in the confirmed cases of dengue fever in the State of São Paulo (Fig. 1B). Public Health authorities of Guarujá reported a total of 1,805 autochthonous dengue cases during this year.

We studied 505 PCR dengue-positive patients during the study In Guarujá. These samples were collected between December 2012 and July 2013. Serotyping determined 10 (1.9%) to be DENV-1, eight (1.5%) to be DENV-2, two to be DENV-3 (0.4%) and 505 to be DENV-4 (96.2%). DENV-4 was relatively new to the country and outbreaks had been reported throughout the country since it was first detected in Brazil in 2011. Preliminary results reporting the documentation of the co-circulation of the four serotypes was published elsewhere (Villabona-Arenas et al., 2014).

Complete envelope (E) gene sequences were obtained for 354 DENV-4 (1,485 bp-long). The remaining DENV-4 was not processed due to technical problems (e.g., did not yield sufficient viral RNA). These sequences were deposited in GenBank under the accessions KP703864–KP704217.

DENV-4 genetic diversity

A total of 109 sites (7.3%) of the envelope gene were polymorphic; 32 sites (2.1%) fell in the first and second codon position and 77 (5.2%) fell in the third codon position. There were a total of 95 haplotypes and 75 of them represented one unique sequence; the three most frequent haplotypes had 119 (34%), 38 (11%) and 30 (8.5%) sequences (Haplotype diversity, the probability that two haplotypes drawn uniformly at random from the population are not the same, was 0.86). The Tajima’s D value was of −2.56 (p < 0.001) (−2.48 for the combination of both first and second codon positions, −2.26 for the third codon position) evidenced an excess of low frequency polymorphisms relative to expectation, indicating population size expansions and/or purifying selection. In agreement with this, the overall rate of non-synonymous over synonymous changes dN/dS value of 0.15 (95% CI 0.10–0.21) for the entire gene suggested purifying selection. A few codons, which fell in the central and dimerization gene domains, showed statistically significant purifying selection (codons 92, 133, 184 and 225) at the significant level of 0.05. Although a few sites experienced an elevation on dN/dS there was no statistical evidence for adaptive evolution.

Evolutionary history and epidemiological dynamics

The 354 DENV-4 local sequences fell into the Latin-American cluster of viruses of Genotype II together with other Brazilian samples (see Data S1).

The best-fit model of nucleotide substitution was Tamura-Nei (TrN) with invariables sites. The mean evolutionary rate was 2.79 × 10⁻³ substitutions per site per year (95% HPDs: 2.06 × 10⁻³–3.77 × 10⁻³ substitutions per site per year). The estimates for the epidemiological parameters were: the sampling proportion was 1.4% (95% HPDs: 0.03–3.1%), the infectious period was 6.7 days (95% HPDs: 5–10 days) and the origin of the epidemic was the 21 December 2012 (95% HPDs: 16 December 2012–26 December 2012). Figure S1 show the extent to which prior information matched the posterior. Sampling from the prior analysis indicated that the posterior and prior traces were the same and that the overall constraints were not forcing the results.

Figure 2A shows the MCC tree with two clades early on in the epidemic. The mean time to the most recent common ancestor of these clades did not differ significantly, suggesting that both viral lineages diverged over similar time-scales, and then co-circulated.

The dynamic of R is presented in Fig. 2B. A value of the parameter R over 1.0 indicates that the disease will be able to spread in a population. For the time series analysis, the estimates are high during the first four months of the year and the curve decrease rapidly by the end of April with values below 1.0 in May and June. We did not included the birth-death skyline plots because this reconstruction may be misleading around the time we had severe sampling issues (due to the switch to clinical diagnosis) and not enough phylogenetic diversity was observed within the sampled genetic data (du Plessis & Stadler, 2015); additional BDSKY analyses where sampling proportions were estimated in a piece-wise manner over six different intervals did not solve this. Moreover, these methods not always reconstruct complex dynamics when other factors such as seasonality, spatial structure and vector dynamics are not incorporated (Rasmussen, Boni & Koelle, 2014).

Phylogeography of DENV-4 over the city

We were able to geolocate 286 patients (81%) (Fig. 1A) based on the addresses recorded by the Guarujá Municipal Health Department (Records were not available for the remaining patients); Fig. 2B compares the actual number of official cases reported and our DENV-4 geolocated sampling. Our first geolocated sample was collected in 2 January 2013 at the neighborhood Enseada. This location has a high number of residents (20,883 based on the 2010 census records) and is home to the largest beach concentrating mostly residents. Pae-Cará together with its neighbor Itapema, are the neighborhoods with more residents (26,054 and 26,070 respectively) followed by Morrinhos (24,387), Enseada and Jardim Boa Esperança (20,753). The digital map was provided by the Municipality’s office and represents the master plan for development and urban planning in the city. Figure 3 illustrates the overall discrete spatial diffusion over the urban area. These results were gauged from a full location-annotated MCC tree, available as Fig. 4, which evidenced an early widespread distribution of the virus in January. The introduction events into each discrete unit are depicted in Fig. 3A; these represent viral diffusion during the onset of the epidemic. These figures suggest that two localities, Enseada and Pae-Cará, were key virus sources. Later on, all regions become interconnected in terms of viral diffusion. The adjacent high-income, low-population density neighborhoods (Jardim Acapulco and Pernambuco) had no cases sampled. This is not explained by distance or lack of connection, because low-income areas nearby (the shantytowns in Mar e Céu, to the South, and Pereque Beach, to the North) had several cases during the epidemic. Bayes factor test of significant diffusion rates shows that another two localities (Morrinhos and Jardim Boa Esperança) played an important role during the course of the epidemic (Fig. 3B). The initial diffusion pattern reproduced to some degree the main access highways of the island: a north-south axis with Pae-Cará and a west-east axis over the littoral with Enseada (Fig. 3C).