Epidemiological and genomic investigation of chikungunya virus in Rio de Janeiro state, Brazil, between 2015 and 2018

Abstract

Since 2014, Brazil has experienced an unprecedented epidemic caused by chikungunya virus (CHIKV), with several waves of East-Central-South-African (ECSA) lineage transmission reported across the country. In 2018, Rio de Janeiro state, the third most populous state in Brazil, reported 41% of all chikungunya cases in the country. Here we use evolutionary and epidemiological analysis to estimate the timescale of CHIKV-ECSA-American lineage and its epidemiological patterns in Rio de Janeiro. We show that the CHIKV-ECSA outbreak in Rio de Janeiro derived from two distinct clades introduced from the Northeast region in mid-2015 (clade RJ1, n = 63/67 genomes from Rio de Janeiro) and mid-2017 (clade RJ2, n = 4/67). We detected evidence for positive selection in non-structural proteins linked with viral replication in the RJ1 clade (clade-defining: nsP4-A481D) and the RJ2 clade (nsP1-D531G). Finally, we estimate the CHIKV-ECSA's basic reproduction number (R0) to be between 1.2 to 1.6 and show that its instantaneous reproduction number (Rt) displays a strong seasonal pattern with peaks in transmission coinciding with periods of high Aedes aegypti transmission potential. Our results highlight the need for continued genomic and epidemiological surveillance of CHIKV in Brazil, particularly during periods of high ecological suitability, and show that selective pressures underline the emergence and evolution of the large urban CHIKV-ECSA outbreak in Rio de Janeiro.

Fig 1. Map of Brazil and Rio de Janeiro state, state-level weekly CHIKV incidence in the 2016–2018 period and R_t estimate.

(A) Map of Brazil with all states colored in light gray, except for Rio de Janeiro, in green. (B) Rio de Janeiro map with municipality borders delimitation. The municipality of Duque de Caxias, where samples for genome sequencing were collected, is colored yellow. (C) State level weekly incidence data. Bars correspond to the number of CHIKV cases. (D) Estimate for Rt between 2016 and 2018 (left y-axis). The solid line indicates mean values, while the ribbon indicates the 95% confidence interval. The dashed line represents index P values, a measure of transmissibility potential for the vector Aedes aegypti (right y-axis). The dotted line marks the critical epidemic threshold (R_t = 1). Maps shapefiles were downloaded from Instituto Brasileiro de Geografia e Estatística (IBGE) at https://www.ibge.gov.br/en/geosciences/territorial-organization/territorial-meshes/18890-municipal-mesh.html (last accessed 26 June 2023).

Fig 2

Symptoms distribution, RT-qPCR Ct values dynamics and their correlation with genome sequencing coverage. (A) Distribution of symptoms exhibited by all CHIKV positive patients in the Duque de Caxias cohort. (B) The time lag between symptoms onset and sample collection dates exhibits correlation with RT-qPCR Ct values. As infection proceeds, viral loads decrease (Cts increase) likely due to immunological response. (C) Negative correlation between RT-qPCR Cts and genome sequencing coverage. Sequences characterized from samples with higher viral load (lower Cts) tend to exhibit higher coverage, although no strong statistical correlation was inferred on a linear model (p = 0.08).

Fig 3. Maximum likelihood phylogenetic analysis of global and ECSA-American datasets.

(A) Phylogenetic tree inferred from the global dataset. All new characterized genomes clustered within the ECSA-Br clade. Names of lineages and relevant clades are indicated. SH-aLRT statistical support values for these clades are shown close to their defining nodes (colored in red). (B) Phylogeny inferred from the filtered ECSA-American dataset. Tip shapes are colored according to sampling location (Centre-West region: purple, North region: yellow, Northeast region: red, Rio de Janeiro state: green). Sequences generated in this study are highlighted with red circles around the tip shapes. Clades composed mostly by RJ sequences are indicated along their SH-aLRT support values. (C) The root-to-tip regression plot, which indicates a strong temporal signal (R² = 0.72, slope = 5.19 x 10⁻⁴). Scale bars represent substitutions per site (s/s).

Fig 4. Bayesian time scaled phylogeographic reconstruction for the clade ECSA-Br.

(A) The molecular clock phylogeny annotated with discrete trait reconstructions. Colors indicate estimated ancestral locations (Centre-West region: purple, North region: yellow, Northeast region: red, Rio de Janeiro state: green). Tip shapes mark sequences generated in this study. The x-axis depicts the timescale, while the density plots indicate the posterior distributions estimated for the age of clades ECSA-American, RJ1 and RJ2. Posterior probabilities for both RJ clades are shown. Positively selected mutations detected with MEME and FEL models (NS1: D351G and NS4: A481D) are exhibited on the branches where they occurred according to the ancestral states reconstruction performed (TreeTime). The inset marks the date of arrival of the index case in the Northeast region. (B) Skygrid reconstruction plot for the clade RJ1. A separate analysis was performed with only RJ sequences from the clade RJ1, allowing the reconstruction of the dynamics of variation of viral effective population size in the state (left y-axis). The analysis reveals CHIKV genetic diversity varied over time, with periods of high diversity matching peaks observed in incidence data (right y-axis).