Introduction

After the introduction of the Asian lineage of chikungunya virus (CHIKV) into the Americas in 2013, and subsequent detection of the East/Central/South African (ECSA) lineage in 2014 in Bahia state Northeast Brazil, more than 2.9 million infections have been reported in Brazil up to 2019 [1]. Common clinical manifestations of the disease include fever, muscle pain, rash and severe joint pain, which may last for months to years [2]. It is argued that 82% of infections caused by the ECSA lineage are symptomatic, and there are suggestions of lineage-specific infection outcomes [3]. The ECSA lineage seems to be the main genotype currently circulating in Brazil since its introduction in the northeastern region in 2014, despite detection of the Asian lineage in that same year in northern Brazil [4]. From previous studies addressing interregional dispersion of the virus [4], CHIKV outbreaks have been registered in northern, northeastern, southeastern and midwestern states of Brazil between 2016 and 2019 [5–8]. The number of probable cases reached 87,687 in 2018 and 132,205 in 2019, with the midwestern region accounting for the second largest number of cases (15.8%) registered in 2018 in Brazil [9]. Despite the large number of cases, much is unknown about the genomic diversity and evolution of CHIKV lineages currently circulating in the midwestern region of Brazil, as well as their dispersion dynamics in South American countries that border this Brazilian region. Thus, to investigate the genomic diversity and evolution of CHIKV, we analyzed 24 new CHIKV genomes generated by next generation sequencing, providing additional information on the introduction and spread of the ECSA lineage in Midwest Brazil as well as in the Americas.

Materials and methods

Viral RNA isolation and quantitative real-time RT-PCR

Serum samples were submitted to nucleic acid purification using the QIAmp Viral RNA Mini kit (Qiagen), following manufacturer’s recommendations. The CHIKV RNA detection by RT-qPCR was performed using a Real Time RT-qPCR protocol adapted from [10], using Promega GoTaq Probe 1-Step RT-PCR System Kit in Bioer LineGene 9660 equipment. Samples were selected for sequencing based on Ct-value <32 to maximize genome coverage of clinical samples by nanopore sequencing [11] (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Epidemiological data for the sequenced samples.

https://doi.org/10.1371/journal.pntd.0009290.t001

Results

We obtained CHIKV RT-qPCR positive clinical samples as part of the genomic surveillance project called ZIBRA 2 (https://www.zibra2project.org/), which aimed to perform, from a lab on wheels, genome sequencing of arboviruses circulating in Midwest Brazil, in the Mato Grosso state. Supported by the Brazilian Ministry of Health and PAHO/WHO, ZIBRA 2 has carried out genomic surveillance of arboviruses in the northern, northeastern and southeastern regions of Brazil [5, 6, 24]. We used the portable MinION sequencer and an amplicon approach [12] to generate 24 partial and near complete CHIKV genomes from serum samples provided by the Central Public Health Laboratory and collected during the 2018 outbreak in the midwestern state of Mato Grosso (MT). These samples, most of which were collected in March 2018, had an average Ct value of 24.48 (ranging from 18.59 to 31.21) and were from patients with average 29 years of age, and in their majority were female (87.5%) living in the city of Cuiaba, MT (Table 1 and Fig 1A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. CHIKV transmission dynamics in Brazil (2014–2019).

(A) Map of Brazil and Americas showing sampling location of the CHIKV genomes from this study (white circles). (B) Weekly notified chikungunya cases normalized per 100K individuals per Brazilian region in 2015–2020 (until EW06). Incidence (cases per 100K population) is presented in log10 for visual purposes. Epidemic curves are coloured according to geographical macro region: SE = Southeast, NE = Northeast, MW = Midwest, N = North, S = South. (C) Time-scaled phylogeographic tree of 120 complete and near-complete CHIKV genome sequences from the ECSA genotype sampled in Brazil and Americas. Colours represent different sampling locations according to the legend on the left of the tree. Tip circles (white) represent the genome sequences generated in this study. Skyline plot (black and gray lines) is superimposed. Relative genetic diversity is represented here as a surrogate for the product of effective population size and generation time. The solid black line represents the mean relative genetic diversity and the light gray area around the line represents the 95% HPD interval of the estimates.

https://doi.org/10.1371/journal.pntd.0009290.g001

The 24 new CHIKV genome sequences from MT had a mean genome coverage of 81.79% (coverage range 43.7%-93.6%). This genome coverage obtained is considered sufficient to perform phylogenetic inferences, according to a study that showed the occurrence of a decrease in phylogenetic accuracy when genome coverage is reduced from 40% to 20% [25].

Genomic data obtained in this study belonged to the ECSA lineage, as confirmed by the chikungunya virus typing tool (https://www.genomedetective.com/) [14], and clustered with other Brazilian sequences from previous outbreaks reported in other geographic regions (Fig 1C). These new genome sequences were submitted to GenBank under the accession numbers MN428504-MN428527.

Fig 1B shows the CHIKV weekly cases normalized per 100K individuals notified between 2015 and 2020 (until epidemiological week–EW 06) in five Brazilian regions: Southeast (SE), Northeast (NE), Midwest (MW), North (N), and South (S). Five CHIKV epidemic waves were found (2016–2020), characterized by a reduction in total cases per year from 2018 to 2020. Although without further information it is difficult to assert the drivers of this reduction in incidence, it is likely that significant herd-immunity has accumulated since CHIKV’s introduction, as suggested in other studies [26, 27]. Weekly reported incidence reveals that chikungunya was mostly reported in the Northeast region between 2015 and 2018, the South region had the lowest incidence in the entire time period, and that contrary to all other regions, the Southeast has presented an increasing incidence trend over the years (Fig 1B).

To explore in more detail the evolutionary relationship of these new MT sequences in the Brazilian context and to infer the ancestral location of CHIKV strains circulating in MT, we used a Bayesian discrete phylogeographic approach, employing the uncorrelated relaxed molecular clock and a Bayesian Skyline coalescent model (and a linear regression of root-to-tip genetic distance against sampling date which revealed sufficient temporal signal, r2 = 0.70; S1 Fig) in a dataset comprising the new MT sequences described in this study plus 96 published CHIKV genomes (Fig 1C). The Bayesian demographic reconstructions (Fig 1C, superimposed to the tree) provided evidence of seasonal oscillations, although with wide credible intervals, highlighting a slight but gradual decline of the median Effective population size (Ne) estimated over 2018–2020 which appears to be in agreement with decreasing incidence in time (Fig 1B).

Our Maximum clade credibility (MCC) tree showed that the new 2018 MT isolates formed a single well-supported clade, with posterior support of 1.0 (Fig 1C). Interestingly, these new isolates did not group with the other six previously published sequences sampled in 2017 also from Cuiabá-MT [28]. As these 2017 MT sequences also formed a well-supported clade (posterior support = 1.0), our MCC tree topology suggests that at least two independent introductions of the ECSA lineage occurred in MT. A previous study reported the circulation of CHIKV ECSA lineage in MT and indicated, by RT-qPCR, that the earliest human case of the disease recorded in that state dates to July 2015 [8]. From our time-measured tree, we estimated the time of the most recent common ancestor (TMRCA) of the two independent introduction events in MT to be between late September 2016 (May 2016 to December 2016, 95% HPD) for the first introduction event, and early February 2017 (September 2016 to June 2017, 95% HPD) for the second event. This interval does not include the earliest CHIKV positive case previously reported in MT and this divergence might display an absence of sufficient data because of non-sampling of earliest isolates [8].

The phylogeographic analysis showed that CHIKV was most likely introduced in MT from a dispersion event originating in the northeastern region of Brazil (location probability 0.98). Clades comprising other isolates from 2017–2018 and sampled in the northeastern region illustrate the persistence and re-emergence of the ECSA lineage in the northeastern region of Brazil since its introduction in 2014. In addition, our tree shows the 2018 MT outbreak clade is closely related to the clade containing isolates from Paraguay sampled in 2018, when 1,237 CHIKV cases were reported in that country [1]. We estimate that the TMRCA of the isolates from Paraguay dates back to December 2017 (June 2017 to May 2018, 95% HPD) and probably originated in the midwestern region of Brazil (location probability 0.90). This would be the second event of cross border transmission from Brazil, as our tree also shows isolates from Haiti sampled in 2016 clustered with an isolate from Northeast Brazil also sampled in 2016 (posterior support = 0.97, location probability 0.97). These results indicate a possible role of Brazil as a source for cross-border dispersion of the CHIKV ECSA lineage to other countries in the Americas since its introduction into the country.

Discussion

More than 930,000 cases have been notified since CHIKV was first detected in Brazil in 2014 (4). Despite this large burden of disease, much is unknown about the origins of the virus responsible for the Brazilian outbreaks. To get more insight regarding CHIKV dispersion through different Brazilian regions and South American countries we generated 24 partial and near complete genome sequences from the 2018-CHIKV-ECSA epidemic registered in the state of Mato Grosso (MT), Midwest of Brazil, using a combined strategy of a mobile sequencing mission through this region, genomic, and epidemiological analysis.

Epidemiological data revealed yearly patterns of CHIKV transmission with a reduction in the number of reported cases over 2018–2020, likely a consequence of an expected, gradual accumulation of herd-immunity over the 7 years since its introduction in 2014. Phylogeographic reconstructions suggest that at least two independent introductions of the ECSA lineage occurred in MT from a dispersion event originating in the northeastern region of Brazil and estimated the time of the most recent common ancestor (TMRCA) to be between late September 2016 for the first introduction event, and early February 2017 for the second event. Furthermore, our analysis suggests that the midwestern Brazilian region appears to have acted as a source of virus transmission towards Paraguay, a bordering South American country to the Brazilian Midwest.

In summary, our data reveals a complex pattern of CHIKV transmission between epidemic seasons and sampled locations and suggests that Brazil has played a role as source for international dispersion (enhanced by cross-border transmission to other Americas countries such as Paraguay and Haiti) of the CHIKV-ECSA genotype to other American countries. Those results highlight the utility of combining genomic, epidemiological and evolutionary methods to understand ongoing mosquito-borne epidemics. Our analyses further indicate that additional data is required to better identify routes of CHIKV-ECSA genotype transmission into Brazil, and to understand its transmission dynamics through other American countries.

Supporting information

Acknowledgments

We thank all personnel from Health Surveillance System from the Mato Grosso and Brazilian Ministry of Health that helped with samples, sources, epidemiological data collection and transportation of reagents.