Genomic characterization of Dengue virus circulation in Ethiopia

Patient characteristics and study settings

We conducted a cross-sectional hospital-based study in Afar, Gambela, and Dire Dawa city administration (Supplementary Figure 1). Between December 2022 and November 2023, we enrolled outpatients and inpatients presenting with fever (>37.5°C) at three major public hospitals in Dubti (Afar region), Dilchora (Dire Dawa), and Gambela. These hospitals serve as sentinel sites for acute febrile illnesses in the country and were selected due to their previous history of arbovirus outbreaks and serological evidence for arboviruses (5, 18).

Recruitment of participants

We obtained informed consent from all participants aged 18 years and older. For individuals aged 1 to 17 years, we obtained informed consent from a parent or guardian and, where appropriate, assent from the child.

Eligibility/exclusion criteria

The inclusion criteria for the study were: a fever of more than 37.5 °C, following the 2009 WHO Dengue classification (19), joint pain, rash, headache, myalgia, retro-ocular discomfort, stomach pain, and hemorrhagic presentation.

Patients with severe, and established chronic clinical illnesses such as people living with HIV, malignancies, and known metabolic disorders were excluded from the study.

Sample collection and processing

Each participant provided five milliliters of whole blood, which was collected into serum separator tubes. At the Afar site, blood sample collection took place between April 25, 2023, and May 15, 2023. The serum was isolated and collected in tubes containing 1 ml DNA/RNA shield and stored at −20 °C during field collections. It was then transported in a cold chain to the laboratory of the Ethiopian Public Health Institute (EPHI) and stored at −20 °C until analysis for PCR and sequencing. Blood samples from Dire Dawa and Gambela sites were collected between June 2023 and October 2023, stored at −20 °C, and then shipped in a cold chain to the laboratory of the EPHI, where they were stored at −80 °C until analysis.

RNA extraction and PCR amplification

We performed RNA extraction within five days of receiving the serum samples at the EPHI laboratory from each site. The RNA extraction was carried out using a Bioer NPA-32P instrument (Bior Technology, Zhejiang, China), which is a high-throughput automated workstation. We used the MagaBio Plus RNA Purification (Bior Technology, Zhejiang, China) for the extraction process, following the manufacturer’s instructions. Briefly, we prepared 96-well plates and added 300 µl of thawed and briefly vortexed serum samples to columns #1 and #7. The plate was then placed into the instrument, which agitated the plates every 10 seconds for 9 minutes until the extraction was completed, resulting in 70 µl RNA extracts.

The remaining serum samples from the Afar region were stored at −20°C for later whole genome sequencing. Meanwhile, samples from the Dire Dawa and Gambela sites were stored at −80°C. In February 2024, some of the qualifying serum samples were sent to the CERI and KRISP laboratories in South Africa for whole genome sequencing.

On the same day that we extracted the RNA, we used the CDC Trioplex real-time RT-PCR assay to screen for arboviral infections such as DENV, CHIKV, and ZIKV. This assay includes a set of published oligonucleotide primers and dual-labeled hydrolysis Taqman probes. To run the test, we combined 10 µl of sample RNA with 12.5 µl of PCR master mix reaction buffer (SuperScript™ III One-Step RT-PCR System with Platinum™ Taq DNA Polymerase, ThermoFischer Scientific), 1 µM virus-specific primers, 0.3 µM dengue-specific probe, 0.15 µM chikungunya-specific probe, 0.15 µM Zika virus-specific probe, and nuclease-free water in a 96-well optical PCR plate to make a final reaction volume of 25 µl. We then carried out reverse transcription (RT) for 30 minutes at 50°C, followed by RT inactivation for 2 minutes at 95°C, fluorescence detection for 15 seconds at 95°C, and annealing for 1 minute at 60°C. Additionally, the assay includes an internal control reaction targeting human endogenous ribonuclease P (RP) to ensure the extracted test specimen contains amplifiable RNA. The fluorescent signal intensity was captured using a real-time PCR instrument called QuantStudio5 Real-Time PCR system (ThermoFisher Scientific).

Dengue serotyping

Following PCR screening, samples positive for DENV were serotyped from the same RNA extract using CDC DENV 1-4 rRT-PCR Multiplex assay on the QuantStudio5 Real-Time PCR system. The assay employed DENV 1-4 specific primers and probes provided by the CDC, USA. To perform the assay, we mixed 5 µl of extracted RNA templates with 12.5 µl of the master mix (SuperScript™ III One-Step RT-PCR System with Platinum™ Taq DNA Polymerase, ThermoFischer Scientific). The assay was run in a multiplex format, with each serotype-specific probe labeled with different probes as per the manufacturer’s protocol.

Dengue sequencing

After serotyping, we selected specimens with PCR screening Ct values below 26 for sequencing. Sequencing was carried out with technical support from CERI and KRISP. The primer scheme for DENV sequencing and other arboviruses is available on the CLIMADE GitHub (https://github.com/CERI-KRISP/CLIMADE/tree/master/Protocols/Arboviruses). The COVIDSeq protocol was adapted to perform the library preparation, followed by sequencing on the Illumina Miseq platform. This protocol is available on the CLIMADE GitHub website, and it is described in detail by Elaine Vieira Santos, Debora Glenda Lima de La Roque in protocols.io (https://protocols.io/view/pathogen-whole-genome-sequencing-multiplexed-ampli-cgwbtxan).

The libraries were prepared using the Illumina COVIDSeq protocol (Illumina Inc, USA). Briefly, RNA was reverse transcribed to cDNA using random hexamers. The Dengue genome was amplified using the two pools of primers specific for the Dengue serotype 1-4. PCR amplicons were tagmented using the EBLTS (Enrichment BLT), which is a process that fragments and tags the PCR amplicons with adapter sequences. Adaptor ligated amplicons were further amplified using the unique 10 base pair Index 1 (i7) adapters and Index 2 (i5) adapters (IDT for Illumina-PCR Indexes Set 1) for each sample. The pooled amplicon library was quantified using a Qubit 2.0 fluorometer (Invitrogen, Inc.) and diluted to a 4 nM. The 4nM library was denatured and diluted to a final loading concentration of 12pM. Dual indexed paired-end sequencing was performed on the Illumina Miseq using a v3 600 cycle flow cell.

Bioinformatics analysis

Following base calling and demultiplexing of the sequence run, fastq files were processed using the Genome Detective analysis pipeline version 2.13.3 (https://www.genomedetective.com/). Upon completion, consensus fasta and Binary Alignment Map (BAM) files were retrieved from Genome Detective for each individual sample. Sequence genotyping and lineage classification was performed using the Genome Detective Dengue Typing Tool (https://www.genomedetective.com/app/typingtool/dengue/), based on a newly developed nomenclature for DENV classification (20). The resulting genomic sequences have been made available on the GISAID Epi-Arbo database (EPI_ISL_19229161 - EPI_ISL_19229193).

Sequence Alignment and Phylogenetic Analysis

DENV reference datasets were retrieved from two primary databases: the National Centre for Biotechnology Information (NCBI) GenBank and the GISAID Epi-Arbo web interface (https://gisaid.org/). Any duplicate entries were removed. Sequence genotyping and lineage classification was performed using the Genome Detective Dengue Typing Tool (https://www.genomedetective.com/app/typingtool/dengue/), based on a newly developed nomenclature for DENV classification (20). The retrieved sequences corresponding to the lineages detected from newly sequenced Ethiopian genomes were subjected to initial quality control to remove any unverified sequences and incomplete records (i.e. geographical location and sampling dates). The final dataset used for the phylogenetic analyses in this study consisted of 33 sequences generated in this study from Ethiopia, as well as 2348 publicly available sequences of DENV3 genotype III major lineage B (corresponding to one of the transmission lineages detected in Ethiopia) and 990 sequences of DENV1 genotype III major lineage A (corresponding to the other transmission lineage detected in Ethiopia). Our DENV1 and DENV3 datasets were aligned against the appropriate DENV serotype reference genome (DENV1: NC_001477.1; DENV3: NC_001475.2) using Nextalign v1.3.0 alignment tool (21).

Maximum likelihood trees for each of the two genotypes were generated using IQ-TREE v 2.2.2.2 with 1,000 bootstraps (22). The nucleotide substitution model TIM2+F+I+G4 was selected for DENV1 and GTR+F+I+G4 for DENV3, respectively, according to Bayesian Information Criterion using Model Finder (23). The molecular clock signal was evaluated using TempEst v1.5.3 (24). Potential outliers that violated the molecular clock assumption were removed prior to inferring a time-scaled tree (TreeTime v 0.10.0) using a molecular clock rate of 5.015 x10^-4 and 1.225×10^-4 for DENV1 and DENV3 respectively nucleotide substitutions per site per year, as determined by the clock function in TreeTime.

Time-calibrated Bayesian phylogenies

Time-calibrated Bayesian phylogenetic trees were constructed to estimate the time of emergence (time to the most common recent ancestor, TMRCA) of the transmission lineages circulating during the 2023 outbreak, and the likely introduction routes. The latest release of BEAST (v.1.10.4) was used along with the BEAGLE library v3.2.0 to improve computational speed (25, 26). A relaxed clock model was applied for all analyses along with the HKY+G+I nucleotide substitution model. A constant population coalescent model assumption was used. Subsequently, all analyses were performed in two independent runs of 100 million iterations each. Convergence of MCMC chains was checked using Tracer v.1.7.1 (27) and 10% of initial chains was discarded as burn-in. Post-burn-in samples were pooled to summarize parameter estimates using LogCombiner and TreeAnnotator, including posterior probability of each parameter and maximum clade credibility (MCC) trees.

Visualisations

The phylogenetic trees were visualized using Figtree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/) and other figures using R and ggplot (28).

Air Travel Data

We evaluated travel data generated from the International Air Transport Association (IATA) to quantify passenger volumes originating from international airports into. IATA data accounts for ∼90% of passenger travel itineraries on commercial flights, excluding transportation via unscheduled charter flights (the remainder is modeled using market intelligence).

Ethical considerations

The study obtained ethical approval from Addis Ababa University (Ref. no. CNCSDO/175/15/2023), and Institutional Review Boards (IRBs) of the EPHI (EPHI-IRB-536-2023).

Sample collection and dengue epidemiology in Ethiopia

Out of 891 febrile patients screened for DENV, CHIKV, and ZIKV infections, only dengue virus was detected in this study. The overall proportion with a positive PCR for DENV was 10.4% (93 out of 891 patients). The test positivity was slightly higher in men than women (Table 1). DENV was isolated from two of the sites; the proportion with a positive DENV PCR was 13.7% (41/300) in Afar and 17.9% (52/291) in Dire Dawa. None of the 300 samples from Gambela were positive for DENV. Virus genotyping shows evidence of co-circulation of DENV1 and DENV3 throughout the outbreak period (Figure 1B). Of the 93 isolates, 88 (94.6%) were identified as belonging to the DENV3 serotype, while the remaining 5 (5.4%) were identified as DENV1. At the Afar site, only DENV3 was detected (n=36), whereas at the Dire Dawa site, both DENV1 (n=5) and DENV3 (n=47) were detected (Table 1).

Virus genome sequencing

A total of 55 viral isolates were selected based on their PCR screening Ct values from Afar region and Dire Dawa City administration. The Ct values of the isolates ranged from 15.5 to 25.6. A total of 20 DENV3 were selected from Afar and 35 (5 DENV1 and 30 DENV3) were selected from Dire Dawa sites. Out of these, 54 isolates were successfully sequenced, and only one sample failed to be sequenced due to low sample volume. Among the successfully sequenced isolates, 33 gave near whole genome sequences (>90% genome coverage). Of these, all of the 5 DENV1 identified as DENV1 genotype III major lineage A (DENV1_IIIA), while the remaining 28 were identified as DENV3 genotype III major lineage B (DENV3_IIIB), with 1 isolate from Afar and 27 from Dire Dawa (Figure 2).

• Download figure
• Open in new tab

Figure 2. Sequencing process and coverage results.

A) Sample selection, genotyping and sequencing workflow. B) CT values versus sequence coverage for all sequenced specimens.

Samples collected from Dire Dawa demonstrated an expected association between a high viral load (lower Ct values) and high sequence coverage whereas most samples from Afar resulted in low sequence coverage genomes despite low Ct values. The Ct values were not repeated after storing the serum samples; we only used the first screening PCR Ct values. The sample collection and storage freezers used were different for both sites. The samples from Afar site were collected in DNA/RNA shield during the outbreak in May 2023 and stored for a long time at −20°C, with intermittent temperature fluctuations of the freezer. In contrast, samples collected from Dire Dawa and Gambela sites were relatively recent and stored at −80°C until they were sequenced.

Evolutionary history of sequences belonging to Ethiopia’s 2023 outbreak

A phylogenetic tree was reconstructed with 990 global DENV1_IIIA aligned sequences, including five from Ethiopia sequenced in this study (Figure 3A). For this lineage, all sequences from Ethiopia belonged to a single transmission cluster, which clustered monophyletically, with other African sequences (Democratic Republic of Congo and Djibouti), suggesting cryptic transmission in the region since 2019 (Figure 3B). TMRCA of the Ethiopian clade specifically to mid 2021 (Figure 3B). Basal to this African clade are sequences mostly of Asian origin suggested an introduction into Africa from Asia, which we can date back to the year 2018 (Figure 3B). For the DENV3_IIIB lineage, we constructed a phylogenetic tree with 2,348 global sequences, including the ones from Ethiopia (Figure 3C). Sequences from Ethiopia belong to two distinct clades on the tree (Figure 3D, E). While the first clade consisted of sequences generally of Asian and African origin, the Ethiopian sequences, with a TMRCA of early 2021 (95% HPD: 2020 – 2022) cluster monophyletically with a sequence from Italy’s 2023 outbreak (29). We infer that the common ancestor between the Ethiopian isolates and the Italy one existed around mid-2019 (95% HPD: mid-2017 - 2021) with long branches separating the two locations. The Ethiopia-Italy clade is supported with >70% posterior support at the relevant internal nodes on the MCC tree. This suggests some level of cryptic transmission in unsampled areas from where both Ethiopia and Italy could have received viral introductions. However, given historical ties between the two countries and the remaining high connectivity, Italy being the European country with the second most air travel passengers into Ethiopia (Supplementary Figure 2), an introduction from there is plausible. The second relevant clade of this lineage shows a clear introduction from Asia to Ethiopia late 2021 (95% HPD: 2021 - 2023), with potential persistent transmission since then (Figure 3E). While the closest relatives to this Ethiopian clade are sequences from India, it remains possible, again given long branch lengths on the tree, that intermediate unsampled regions are involved, including country such as Saudi Arabia (30), which is the top countries with largest passenger volume to Ethiopia (Supplementary Figure 2).

• Download figure
• Open in new tab

Figure 3. Phylogenetic analysis of Ethiopian genomes.

A) Timed maximum likelihood phylogeny of the DENV1 genotype III major lineage A clade, containing Ethiopian genomes sequenced in this study. The tree is coloured by sequence locations. B) Zoom-in of a sub-clade of the tree to show close evolutionary relationships related to the Ethiopian DENV1 genomes. C) Timed maximum likelihood phylogeny of the DENV3 genotype III major lineage B clade, containing Ethiopian genomes sequenced in this study. The tree is coloured by sequence locations. D) Maximum clade credibility (MCC) tree of sub-tree containing cluster 1 of Ethiopian DENV3 genomes. E) Maximum clade credibility (MCC) tree of sub-tree containing cluster 2 of Ethiopian DENV3 genomes.