Retrospective Spatio-Temporal Dynamics of Dengue Virus 1, 2 and 4 in Paraguay

Dengue virus (DENV) has been a major public health concern in Paraguay, with frequent outbreaks occurring since early 1988. Although control measures have been implemented, dengue remains a significant health threat in the country, and continued efforts are required for prevention and control. In response to that, in collaboration with the Central Public Health Laboratory in Asunción, we conducted a portable whole-genome sequencing and phylodynamic analysis to investigate DENV viral strains circulating in Paraguay over the past epidemics. Our genomic surveillance activities revealed the co-circulation of multiple DENV serotypes: DENV-1 genotype V, the emerging DENV-2 genotype III, BR4-L2 clade, and DENV-4 genotype II. Results additionally highlight the possible role of Brazil as a source for the international dispersion of different viral strains to other countries in the Americas emphasizing the need for increased surveillance across the borders, for the early detection and response to outbreaks. This, in turn, emphasizes the critical role of genomic surveillance in monitoring and understanding arbovirus transmission and persistence locally and over long distances.


Introduction
Emerging and re-emerging arboviral diseases such as the dengue (DENV), Zika (ZIKV), and chikungunya (CHIKV) viruses pose a significant public health concern in Latin America, including Paraguay [1]. The country's sub-tropical climate with high temperatures and humidity provides favorable conditions for the transmission of arboviruses, which are viruses transmitted to humans and other animals by blood-feeding arthropods such as mosquitoes [2]. Paraguay, a landlocked country in South America that shares borders with Brazil, Argentina, and Bolivia, has been severely affected by recent epidemics of ZIKV, CHIKV, and dengue fever [3].
DENV is the most common arboviral disease in Paraguay, with regular outbreaks occurring since 1988. All four DENV serotypes (DENV-1-4) have been detected in the country, with multiple serotypes co-circulating in some seasons. DENV-1 was first detected in 1988 and has been responsible for several epidemics in subsequent years [1]. DENV-2 was first isolated in Paraguay in 2001 and circulated at low transmission rates until 2005. Since 2010, the circulation of DENV-2 has been confirmed again and has caused significant epidemics over the years, resulting in an increase in severe cases and deaths [1]. DENV-3 was first detected during small dengue outbreaks in Paraguay in 2002 and 2003, and it has continued to circulate at a low level since then. As for DENV-4, it was first isolated in Paraguay in 2012, and it has continued to circulate seasonally without causing any major epidemics [1].
Despite the implementation of several measures to control the spread of arboviral diseases, little is still known regarding the viral landscape in the country [2]. Genomic surveillance is an extra layer towards control and outbreak response that provides key information about the diversity of circulating viral strains [4]. In this context, using a combination of portable whole-genome sequencing and genomic epidemiology, we found several cases of patients presenting febrile illness in different districts to be infected with different dengue virus serotypes, including DENV-1, DENV-2, and DENV-4. This study highlights the importance of active viral monitoring to prevent future epidemics, given the co-circulation of multiple DENV serotypes, and the crucial role of genomic surveillance in monitoring and understanding the transmission and persistence of arboviruses, both locally and over long distances.

Materials and Methods
Serum samples (n = 99) retrieved from patients presenting with symptoms compatible with an arboviral infection were collected by the Laboratorio Central de Salud Publica of Paraguay in Asunción for molecular diagnosis. Samples were submitted first to nucleic acid extraction using the QIAamp Viral RNA Mini Kit (QIAGEN) and then subjected to real-time reverse transcription PCR to detect ZIKV, CHIKV, and DENV serotypes 1-4 as described previously [5][6][7]. Positive samples were selected for sequencing based on the cycle threshold value ≤ 35 and the availability of demographic metadata such as sex, age, and municipality of residency. Extracted RNA was first converted to cDNA using the SuperScript IV Reverse Transcriptase kit (Invitrogen) and subjected to sequencing multiplex PCR (35 cycles) as previously described [8]. DNA library preparation was conducted using the Ligation Sequencing kit (Oxford Nanopore Technologies) and Native Barcoding Expansion 1-96 kit (Oxford Nanopore Technologies), following the reaction conditions as previously described in [8]. Sequencing was performed for up to 24 h on a MinION device and consensus sequences were obtained using the Genome Detective software [9].
The arbovirus genotyping tool was used to investigate sequence genotypes [10]. Subsequently, to analyze the origins and spatial dynamics of different dengue virus serotypes detected in Paraguay, we downloaded all sequences assigned as DENV-1, DENV-2, and DENV-4 from GenBank collected up to 15 March 2023. We excluded sequences without a sampling date and location as well as sequences that covered less than 50% of the virus genome.
All sequences were aligned using MAFFT [11] and manually edited to remove artifacts using Aliview [12]. The GTR nucleotide substitution model, which was inferred as the best-fit model by the ModelFinder application implemented in IQ-TREE2 [13], was used to estimate maximum Likelihood (ML) phylogenetic trees. The tree topology's robustness was determined using 1000 bootstrap replicates. TempEst [14] was used to assess the presence of a temporal signal, and the BEAST package [15] was used to infer time-scaled phylogenetic trees. To estimate the most appropriate molecular clock model for the Bayesian phylogenetic analysis, we used a stringent model selection analysis that included both pathsampling (PS) and steppingstone (SS) procedures [16]. For all datasets, the uncorrelated relaxed molecular clock model was chosen by estimating marginal likelihoods using the codon-based SRD06 model of nucleotide substitution and the nonparametric Bayesian Skyline coalescent model. To model the phylogenetic diffusion of detected community transmission clades we used a flexible relaxed random walk diffusion model [17,18] that accommodates branch-specific variation in rates of dispersal with a Cauchy distribution and a jitter window site of 0.01 [19,20]. Latitude and longitude coordinates were assigned to each sequence. BEAST v1.10.4 was used for the MCMC analyses, which were run in duplicate for 100 million interactions and sampled every 10,000 steps in the chain. Tracer was used to assess convergence for each run (effective sample size for all relevant model parameters > 200). After discarding the initial 10% as burn-in, MCC trees for each run were summarized using TreeAnnotator. Finally, we extracted and mapped spatiotemporal information embedded in the posterior trees using the R package 'seraphim' version 1.0 [20].
Epidemiological time series data for confirmed, suspected, and probable DENV infections were downloaded from the PAHO website [21]. When presenting the series, we show "all cases" (all categories summed) and "confirmed" (confirmed category) separately. Index P estimations for Paraguay were obtained from Taishi et al. (2023) for Paraguay that were available per month and up to the year 2019 [22].
Seven of the seventeen departments reported the co-circulation of multiple viral strains ( Figure 1A). The samples had a mean RT-qPCR cycle threshold value of 24, and their associated average age was 34 for females and 38 for males (years), with 57% of patients identified as female ( Table 1). All patients were classified as having dengue with no warning signs, and no one was notified as severe. Seven of the seventeen departments reported the co-circulation of multiple viral strains ( Figure 1A). The samples had a mean RT-qPCR cycle threshold value of 24, and their associated average age was 34 for females and 38 for males (years), with 57% of patients identified as female ( Table 1). All patients were classified as having dengue with no warning signs, and no one was notified as severe.
Overall, the 99 novel DENV genome sequences had a mean genome coverage of 90.0% (coverage range 70-94%), with a mean depth of >1000× for all samples. All sequenced samples were collected between February 2014 and December 2022 and were made available on the GenBank database under the accession numbers OQ567774-OQ567872 (details in Table 1). The DENV typing tool classified all DENV-1 genomes generated here as genotype V, all DENV-2 genomes as genotype III, and all DENV-4 genomes as genotype II (Table 1).
Notified DENV cases revealed yearly epidemic waves recurrently taking place in the late summer months ( Figure 1B). The estimated transmission potential, as measured by the mosquito-viral suitability index P, was found to capture the seasonal timing of notified cases over the years. However, these estimations were only available per month up to the year 2019. We compared the dynamics of three climatics (humidity, temperature, and precipitation) and index P versus notified DENV cases by calculating the Spearman's correlation coefficient ( Table 2). The highest coefficient was found between the index P and DENV notifications at 0.55, which increased to 0.71 when applying a lag of +1 month to the index P series ( Figure 1B). Overall, the 99 novel DENV genome sequences had a mean genome coverage of 90.0% (coverage range 70-94%), with a mean depth of >1000× for all samples. All sequenced samples were collected between February 2014 and December 2022 and were made available on the GenBank database under the accession numbers OQ567774-OQ567872 (details in Table 1). The DENV typing tool classified all DENV-1 genomes generated here as genotype V, all DENV-2 genomes as genotype III, and all DENV-4 genomes as genotype II (Table 1).
Notified DENV cases revealed yearly epidemic waves recurrently taking place in the late summer months ( Figure 1B). The estimated transmission potential, as measured by the mosquito-viral suitability index P, was found to capture the seasonal timing of notified cases over the years. However, these estimations were only available per month up to the year 2019. We compared the dynamics of three climatics (humidity, temperature, and precipitation) and index P versus notified DENV cases by calculating the Spearman's correlation coefficient ( Table 2). The highest coefficient was found between the index P and DENV notifications at 0.55, which increased to 0.71 when applying a lag of +1 month to the index P series ( Figure 1B).
To explore the phylodynamic history of DENV-1, we then combined our 47 newly generated sequences with other DENV-1 genotype V genomes available on GenBank (n = 501). Phylogenetic analysis revealed that the novel isolates were organized into four distinct clades, named hereafter as clades I, II, III, and IV ( Figure 2A).  To explore the phylodynamic history of DENV-1, we then combined our 47 newly generated sequences with other DENV-1 genotype V genomes available on GenBank (n = 501). Phylogenetic analysis revealed that the novel isolates were organized into four distinct clades, named hereafter as clades I, II, III, and IV (Figure 2A).  Our analysis additionally suggested that the Brazilian southeast region was the leading hub responsible for multiple viral introductions in Paraguay. Among those clades, clade III resulted in a robust cluster, suggesting a likely persistence of that viral strain in the country over a period of 2 years (Figure 2A). To investigate the evolution of clade III in more detail, we used a smaller data set (n = 21) derived from this clade individually. An analysis of substitution rate constancy revealed a strong correlation between the sampling date and the root-to-tip divergence in this data set (r 2 = 0.70, coefficient correlation= 0.84). Phylogeographic analyses of clade III allowed the reconstruction of viral movements across different districts in Paraguay ( Figure 1B) and suggested a mean time of origin as mid-May 2015 (95% highest posterior density (HPD): 1 May 2015 to 6 June 2015). Viruses from this clade spread from the midwestern district (Distrito Capital) towards the southeast and later to the midwest, as indicated by isolates from the Presidente Hayes region (Figure 2A).
To explore the phylodynamic of DENV-2, we also performed a phylogenetic analysis of the 27 newly generated sequences plus 620 complete genome sequences of DENV-2 III available on GenBank ( Figure 2B). Interestingly, our analysis allowed the detection of the emerging BR-4 L2 clade of DENV-2 genotype III in Paraguay. This clade was first detected at the end of 2019 in Brazil [21], and the introduction into Paraguay was likely mediated by the Brazilian Midwest region, reflecting a cross-border transmission that occurred in early July 2021 (95% highest posterior density (HPD): 21 December 2019 to 16 May 2021).
To explore the phylodynamic history of DENV-4, we performed a phylogenetic analysis of the 25 newly generated sequences plus 132 complete genome sequences of the DENV-4 genotype II available on GenBank ( Figure 2C). We found that 96% (24 of 25) of the novel isolates from Paraguay fell within a single, large, well-supported monophyletic clade (bootstrap score BS = 100%). We also identified a single isolate outside the main clade. This isolate, sampled in March 2016, falls in a clade containing sequences from southeast Brazil, reinforcing the possible role of Brazil as a source for the international dispersion of different viral strains to other countries from the Americas. To gain more insight into the transmission dynamics of the large DENV-4 II clade from Paraguay, we built a dataset containing only isolates from Paraguay belonging to this clade (n= 33) and submitted it to phylodynamic inference. A regression of genetic divergence from root to tip against sampling dates confirmed that this subset had sufficient temporal signal (coefficient correlation = 0.82, r 2 = 0.70). Phylogeographic analyses ( Figure 2C) suggested a mean time of origin as late-April 2018 (95% highest posterior density (HPD): 16 April 2018 to 6 May 2018). Viruses from this clade spread from the midwestern district (Distrito Capital) towards the southeast and later to the northeast, as indicated by isolates from the Alto Parana region ( Figure 2C).

Discussion
Since 1988, Paraguay has experienced significant dengue epidemics with multiple serotypes circulating concurrently during epidemic seasons. This increases the risk of secondary infections, which have a higher likelihood of severe illness and fatal outcomes [22].
DENV-1 was first detected in 1988 and has since been associated with ongoing epidemics [22]. DENV-2 was detected in 2001 and has caused significant epidemics linked to an increase in fatalities and severe cases since 2010 [22][23][24]. DENV-3 was first confirmed during small outbreaks in 2002 and 2003 and has since circulated at a low proportion, while DENV-4 was first reported in 2012 and has since circulated seasonally without causing significant epidemics [1]. The circulation of these serotypes has been characterized by multiple viral introduction events in the country, some of which were connected to viral strains circulating in other South American countries, such as Brazil [3].
Despite the co-circulation of different viral strains and the high disease burden, much remains unknown about the origins and routes of transmission of such viruses in the country. To obtain a better understanding of the DENV evolution in Paraguay, using whole genome sequencing and epidemiological analysis, we generated 99 nearly complete genome sequences collected between February 2014 and December 2022.
The newly generated DENV-1 sequences were classified as genotype V and formed four distinct clades (I-IV), one of which, clade III, resulted in a robust cluster, indicating the strain's likely persistence in the country over a two-year period. Our findings also indicated that the Brazilian southeast region was the primary source of multiple viral introductions in Paraguay for this serotype. All new DENV-2 complete genome sequences obtained in this study were identified as belonging to genotype III, which has previously been found in South American epidemics. Interestingly, our analysis enabled the detection in Paraguay of the emerging DENV-2 genotype III, the BR-4 L2 clade, which was first detected in Brazil by the end of 2019. In addition, we found evidence of the co-circulation of the DENV-4 genotype II during the study period. Our findings indicated that most of the novel isolates from Paraguay from this serotype belonged to a single, large, well-supported monophyletic clade, implying the serotype's persistence in the region over a three-year period. Surprisingly, this serotype is not the most common in other South American countries, including Brazil, where this viral strain has not been seen since 2020.
Notably DENV-3 was not detected in the study. However, it is unclear why DENV-3 is not currently circulating in Paraguay. One possibility is that the virus has been eradicated or reduced to such low levels that it is not detected by routine surveillance methods. Alternatively, DENV-3 may still be present in the country, but at a very low prevalence or in localized areas that were not sampled in the study. It is also possible that the epidemiological factors or transmission dynamics of DENV-3 in Paraguay are different from those of other serotypes, making it less likely to circulate or cause large outbreaks. Further research is necessary to gain a better understanding of the circulation and persistence of DENV-3 in Paraguay.
Overall, our study reveals a complex pattern of DENV transmission between epidemic seasons and sampling locations, with the southeastern and midwestern Brazilian regions acting as sources for international dispersion. Our findings highlight the need to increase genomic surveillance across borders to enable early detection and response to viral outbreaks. However, the scarcity of complete DENV genome sequences in South America limits our ability to characterize the molecular epidemiology of viral strains at a regional level, emphasizing the importance of increasing sequencing efforts to improve real-time data generation, sharing, and representativeness.