Dynamics and ecology of a multistage expansion of Oropouche virus in Brazil

Oropouche virus (OROV)¹ is an arthropod-borne virus first identified in 1955 in Oropouche, a village in Trinidad and Tobago¹. OROV causes a febrile illness² that can progress to severe neurological disease^1,3 and has been associated with maternal vertical transmission⁴. This re-emerging virus circulates primarily among wildlife, such as non-human primates and sloths, with antibody detections in rodents and birds^5,6. During the last 60 years, OROV has caused around 30 documented human outbreaks in the Amazon region^7,8,9,10. The midge Culicoides paraensis serves as the primary vector for human transmission^9,11,12, but other suspected secondary vectors include Culex quinquefasciatus^13,14,15, Coquillettidia venezuelensis, Aedes (Ochlerotatus) serratus^1,6, Aedes aegypti, Aedes albopictus and C. quinquefasciatus^15,16. However, there is no evidence to date of vector competence in any mosquito species¹⁷.

In March 2024, the Pan American Health Organization issued an alert in response to a rapid increase in Oropouche fever cases across several countries^18,19. Brazil has been particularly affected, reporting not only the most cases, but also severe complications linked to OROV infection, including fatal cases^4,18,20, maternal and fetal impacts^21,22, with ~8,318 confirmed cases in Brazil alone as of epidemiological week 40 of 2024. Recent epidemiological data and genomic investigations in Brazil^{23,24,25,26,27} have described the recent expansion of OROV into previously non-endemic regions and other states outside the Amazon Basin, where OROV had not been reported or had only sporadic detections in past decades. Some of these studies have identified reassortment events in the virus genome that may have contributed to its changing epidemiology^{23,24,25,26,27}. While the exact role of reassortment in the adaptation of OROV to new environments remains to be fully understood, it has possibly impacted its spread into new ecological niches (that is, the ranges of environmental conditions suitable for the local circulation of the virus)^{23,24,25,26,27}.

As with other arboviruses²⁸, recent changes in ecological context, such as deforestation, urbanization, human mobility and climate change, have possibly contributed to the emergence of OROV in new regions²⁹. In particular, environmental disruption pushes non-human mammal reservoirs and vectors into closer contact with human populations, facilitating viral spread³⁰. Additionally, human activities^31,32 such as urban expansion and altered land use increase the risk of transmission of OROV to humans in sylvatic and peri-urban settings, where vectors such as C. paraensis thrive^33,34,35. Recent studies have used phylogenetic methods to reveal the timing and spatial paths of the emergence and spread of OROV_BR-2015-2024 lineage in Amazonian subregions²³ or have investigated environmental factors linked with OROV transmission using the association of epidemiological data with land cover types^26,27. Nonetheless, epidemiological, ecological and genomic data have not yet been analysed in an integrated way and much remains unknown about the broader disease ecology of OROV, particularly concerning environmental correlates of local circulation and of its recent geographical range expansion. This lack of comprehensive understanding hampers effective risk assessment and preparedness efforts, both in Brazil and across the Americas.

Our study aimed to test key epidemiological hypotheses regarding OROV disease ecology and its range expansion, both spatially and temporally, specifically: (1) that the dispersal history and dynamics of OROV in Brazil followed a multistep pattern over different time periods; (2) that recent changes in OROV distribution reflect a spatial redefinition of its ecological niche, with environmental suitability expanding from historically restricted Amazonian areas to newly permissive regions across Brazil; and (3) that subtle variations in C. paraensis suitability may have facilitated the recent expansion of OROV. We integrated georeferenced pathogen genomes, epidemiological and geospatial data through complementary analytical frameworks.

Dispersal history and dynamics of OROV in Brazil

The epidemiological dynamics of OROV expansion in Brazil in late 2023 and 2024 show a two-stage process: a rapid rise in cases in Amazonian states, namely Acre, Amapá, Amazonas, Pará, Rondônia and Roraima, as well as parts of Mato Grosso and Maranhão, which together encompass most of the Amazon Basin, particularly centred around Manaus (capital of the state of Amazonas), followed by widespread circulation in other parts of the country (Supplementary Fig. 1). To investigate this epidemiological pattern, we applied a continuous phylogeographic approach using over ~500 genomes generated by previous studies^23,24 and sampled between 2022 and 2024. Building on prior knowledge of reassortment events across genome segments^23,24, we conducted separate phylogeographic reconstructions for segments L, M and S, acting as experimental replicates for downstream analyses. Our analysis extracted spatiotemporal dispersal information from the branches of 100 spatially explicit phylogenies subsampled from each post-burn-in posterior distribution, revealing new insights about OROV dispersal events between sampled regions, as detailed below (Fig. 1a).

The earliest lineage dispersal events were restricted to the Amazon basin and inferred to be before 2020, with the virus gradually spreading to other Brazilian states in a southeasterly direction in 2023 and 2024. The apparent discrepancies between the time inferred for the most recent common ancestor (tMRCA) of each segment (Fig. 1a) probably reflect slightly different evolutionary histories among those segments. Indeed, reassortment between genomic segments means that they have distinct genealogies. However, these tMRCA estimates are not substantially different as they have large overlapping 95% highest posterior density (HPD) intervals, which could also be related to the weak yet substantial temporal signal evaluated for those alignments (Methods; March 2017 to September 2020 for segment L, August 2016 to December 2019 for segment M and July 2015 to April 2018 for segment S). By further examining the spatial dissemination of OROV lineages across Brazil, we found that the virus reached a maximal wavefront distance (furthest extent of the wavefront from the location of the inferred tree root position) of over 3,000 km from its epidemic origin through the entire dissemination period (Fig. 1b). This reconstruction further captured two distinct expansion phases during which the wavefront distance increased rapidly, indicating the introduction into new areas in 2024 (Fig. 1b). The rapid expansion of the wavefront in 2024 corresponded to the increase in cases in areas outside the Amazon (Supplementary Fig. 1). Furthermore, diffusion coefficient estimates changed through time, with increases in the first 6 months of 2023 and 2024, corresponding to the expansion phases and consistently detected in all three OROV segments (Fig. 1c), reflecting a substantial dispersal capacity. In 2024, a non-negligible proportion of inferred lineage dispersal events occurred within the states of Mato Grosso, Minas Gerais and Bahia (~46%, ~45% and ~43% for segments L, M and S, respectively), that is some Brazilian states where there was no previous sustained OROV circulation.

When inspecting the association between the geographic distance (km) travelled by each phylogenetic branch and the associated diffusion velocity (km² per day) of the inferred lineage dispersal event, we observed a large group of phylogenetic branches associated with relatively short dispersal distance (<20 km) and slow diffusion (<4 km² per day; Fig. 1d), as well as a subsequent group of phylogenetic branches that correspond to faster long-distance dispersal events (that is, lineage dispersal events associated with notably high diffusion velocity; Fig. 1d). This supports a multiscale expansion process with a combination of short-distance dispersal events and fast long-distance jumps, with some of the latter probably reflecting human-mediated virus movements through air traffic. We also identified signals of substantial isolation-by-distance (IBD), with Pearson correlations between patristic and log-transformed geographic distances between samples close or greater than 0.5 for all three segments (0.470 (95% HPD = [0.360, 0.599]) for segment L, 0.575 (95% HPD = [0.286, 0.607]) for segment M and 0.684 (95% HPD = [0.326, 0.710]) for segment S); the patristic distance being defined as the sum of the branch lengths that link two tip nodes in a phylogenetic tree. Overall, these dispersal metrics provide a clear indication of rapid long-distance dispersal events during the recent OROV expansion in 2024 beyond the Amazon Basin, probably human-mediated through air travel, followed by more localized viral circulation.

Ecological factors associated with dispersal of OROV in Brazil

To investigate the ecological factors associated with the spatial expansion of OROV in Brazil, we analysed virus dispersal history in relation to 28 environmental factors (including land-use, climatic and demographic variables; Supplementary Fig. 2). This includes a species distribution model for C. paraensis, the recognized vector for OROV^1,9,11,12,36. Using high-resolution environmental rasters for each covariate (Supplementary Fig. 2), we extracted environmental values at the geographic location of each phylogenetic node. These correspond to inferred positions through continuous phylogeographic analysis for internal nodes and to sampling locations (down to municipality level) for tip nodes (Supplementary Fig. 3). Specifically, we assessed shifts in dispersal environments relative to four specific periods: (1) before and (2) after the re-emergence and epidemic expansion in the Amazon (< and > mid-2023), (3) during the Amazon-restricted transmission phase (<2024) and (4) during the national expansion phase (>2024).

Before formally testing these associations, we first visually explored the temporal evolution of the environmental conditions associated with the dispersal locations of inferred lineages. For this purpose, we extracted and averaged the environmental values crossed by lineage dispersal during successive time slices of ~2 weeks and while considering a sliding window of ~1 month. These visual explorations indicate that, when comparing dispersal events before and after 2023–2024, different environmental conditions for OROV circulation are already noticeable (Fig. 2): certain dispersal environments seem to shift at the mid-2023 time point, while others at the 2024 time point (Fig. 2 and Supplementary Fig. 3). For instance, before 2024, dispersal events tended to occur in areas with relatively higher evergreen broadleaf forest cover and higher mean annual temperature (Fig. 2 and Supplementary Fig. 3). Interestingly, most lineage dispersal locations before 2024 were associated with a mean temperature of ~27 °C, which probably corresponds to the thermal optimum for OROV transmission suitability within our model, rather than the ecological optimum for C. paraensis abundance³⁷. However, dispersal environments were already shifting towards areas with comparatively higher human population density and urbanization around mid-2023 (Fig. 2 and Supplementary Fig. 3). With strong Bayes factor (BF) support (>20)³⁸, our analyses demonstrated that these trends in dispersal environments were consistent across the posterior distribution of trees obtained through continuous phylogeographic inference.

To statistically test the association of environmental conditions with the spread of OROV lineages, we used a landscape phylogeographic approach to analyse the values extracted at the tree node positions for all 28 environmental factors considered in this study³⁹. Specifically, for the four distinct time periods mentioned above (<mid-2023, >mid-2023, <2024 and >2024), we tested whether inferred OROV lineages tended to preferentially circulate in or avoid circulating in certain environmental conditions. Statistical support (BF) was obtained by comparing the results from phylogeographic reconstructions with a null dispersal model, in which a new continuous diffusion process was randomly simulated along the same tree topologies. During both transmission phases divided by the 2024 cut-off, our results reveal strong support (that is, BF > 20 for at least two out of three segments) for preferential circulation of inferred OROV lineages in areas with comparatively higher population density and urbanization, lower evergreen broadleaf forest cover, areas associated with cocoa cultivation and higher ecological suitability for the C. paraensis vector (Fig. 2 and Supplementary Table 1), which is also found in areas without banana and cacao plantations. We also found strong support for a preferential circulation of inferred OROV lineages in areas associated with banana cultivation in the expansion phase (>2024), as well as in the pre-expansion phase (<2024) for segment L (Supplementary Table 1). These results are consistent with existing knowledge about C. paraensis larvae developing in microhabitats of decaying debris from banana and cacao plantations^40,41,42. Specific to the Amazon-only transmission phase (<2024), our results also reveal a preferential circulation in areas with comparatively higher mean annual temperatures (BF > 20; Fig. 2c). Of note, we also tested the association of forest loss with the OROV spread but it did not lead to strong support in any of the time periods considered. However, an examination of the environmental conditions associated with lineage dispersal locations between transmission phases before and after mid-2023 revealed compelling differences: for instance, for viral lineages inferred before mid-2023, there was no strong support for a preferential circulation of inferred OROV lineages in areas more ecologically suitable for the C. paraensis vector, areas associated with a comparatively higher mean annual temperature, lower evergreen broadleaf forest coverage and denser human populations, although there was still a preferential circulation in more urbanized areas (Supplementary Table 1). This suggests ecological differences underlying endemic OROV circulation within the Amazon before the re-emergence and rapid epidemic expansion at the end of 2023. Overall, our findings reveal similar trends across the three segments analysed (Supplementary Table 1).

Mapping ecological niches for OROV transmission and range expansion in Brazil

The rising number of human OROV detections provides an opportunity to apply modelling approaches to identify areas ecologically suitable for local virus circulation and potential human infections. We used an ensemble modelling approach to reveal OROV transmission suitability across Brazil. The suitability index ranges from 0 (unsuitable conditions) to 1 (highly suitable conditions) and illustrates the potential geographic areas where environmental conditions are most favourable for OROV transmission (Fig. 3a). The disease presence points used as model input represent OROV circulation leading to human cases from 450 geocoded sampling locations in Brazil (based on molecular testing and sequencing records) from 1957 to 2024 (~85% corresponding to 2023–2024). Given this large time span, climatic variables were matched to the corresponding decade of the occurrence point. In the context of presence-only ecological niche modelling, where true disease absence data are unavailable, we sampled pseudo-absence points. Pseudo-absence points sampling was informed by a kernel density estimate of human population density to reflect surveillance efforts, which we assumed is proportional to human population density, with an exclusion radius around presence points (Supplementary Fig. 4a). The model incorporated the same environmental covariates used in the landscape phylogeography analyses (Supplementary Fig. 1), while testing for overfitting along with model performance (Supplementary Table 3). A principal component analysis was performed to assess multicollinearity between the environmental covariates (Supplementary Fig. 4b), resulting in a final selection of eight variables: annual mean temperature, evergreen broadleaf forest cover, croplands cover, elevation, human population density, annual mean monthly precipitation, as well as cocoa- and banana-harvested area coverages. We also report variability among our model predictions as a measure of uncertainty (Supplementary Fig. 4c). Uncertainty is higher in the central regions of Brazil (Supplementary Fig. 4c), where limited surveillance and sparse sampling reduce the accuracy of our predictions. Overall, our findings indicate that the highest ecological suitability for local OROV circulation (>0.7) is concentrated in the northern regions (Fig. 3a), particularly within the Amazon Basin, which aligns with previous studies identifying this area as a significant epicentre for the virus^23,24. Additionally, moderate- to high-suitability areas extend towards the northeast and central-western parts of Brazil, particularly in states such as Pará, Maranhão, Bahia and Mato Grosso, suggesting an expansion into previously non-endemic regions. Our model also highlights the potential for OROV transmission to move beyond traditional transmission range, affecting urban and peri-urban areas, particularly in the northeast (Fig. 3a), where human interaction with vectors may amplify transmission risks. This trend is especially notable for highly suitable areas around the coast of Brazil, where ~111.28 million people—corresponding to 54.8% of the population—reside (Fig. 2b). For cities situated within or near extensive forested areas, proximity to natural environments complicates the identification of precise exposure risks, as human–wildlife interactions and potential transmission events may occur in peri-urban areas, within the urban fringe or deeper within the surrounding forests.

To determine the individual contributions of each environmental factor to our ecological suitability prediction, we further calculated their relative importance (RI) in the resulting ecological niche models. We found that evergreen broadleaf forest cover and human population density contributed the most to the model predictions, followed by precipitation and banana and cocoa agricultural lands (Fig. 3d). This aligns with factors identified in the independent landscape phylogeographic analyses and provides important insights into the disease ecology of OROV. We plotted response curves to assess the relationship between the environmental factors and OROV transmission suitability. These curves illustrate how ecological suitability varies with changes in one factor while all others are kept constant at their mean. We observe a clear tipping point in ecological suitability with mean monthly precipitation, where environments with <2.5 mm mean monthly precipitation appear unsuitable for OROV transmission (ecological suitability ~0) whereas suitability increases considerably above this threshold (ecological suitability >0.4; Fig. 3d). For temperature, there is a clear increase in OROV suitability in environments with annual mean temperatures above 25–27 °C (Fig. 3d), reflecting ecological suitability for viral circulation rather than vector abundance. OROV transmission suitability decreases considerably at high elevations, while relatively higher human population density, evergreen broadleaf forest cover and higher areas of banana and cocoa agriculture appear to be associated with an enhanced risk of local OROV circulation (Fig. 3d).

To investigate whether the expansion of OROV in Brazil in late 2023 and 2024 was associated with an expansion of its ecological niche, we compared ecological niche models trained on pre-mid-2023, pre-2024 and all available occurrence records. Specifically, we compared both the resulting maps of OROV ecological suitability and the capacity of each category of models to predict the distribution of most recent occurrence data. The suitability ranges obtained with the pre-2024 models (Fig. 3b) highlight a similar spatial distribution of highly suitable regions compared with that obtained with the full models (Fig. 3a). While the full suitability estimates show a slight expansion of suitable areas, the pre-2024 model was able to predict the 2024 occurrence points with a relatively high performance (true skills statistic (TSS) = 0.60, area under curve (AUC) = 0.85; Supplementary Table 2). However, the pre-mid-2023 suitability range (Fig. 3c) shows clear differences from the full model, particularly in the landscape of ecological suitability within the Amazon; and the models had a lower predictive ability for occurrence points sampled after mid-2023 (TSS = 0.47, AUC = 0.76; Supplementary Table 2). While some of this effect could be due to the limited sample size of the pre-mid-2023 data (n = 89), this also suggests a shift in ecological suitability after mid-2023, associated with OROV amplification within the Amazon, before circulation in other parts of the country. This shift could have facilitated viral circulation outside the Amazon, which was probably associated with OROV lineages recently reaching areas that were newly ecologically suitable for local OROV transmission. This is supported by our earlier findings of long-distance and rapid-dispersal events from continuous phylogeography.

Integrating results from the ecological niche modelling with the inferred dispersal histories, we further assessed the estimated viral and vector ecological suitability at each dispersal location associated with OROV circulation across space and time (Fig. 3e,f). Examining viral lineage dispersal locations inferred before mid-2023, after mid-2023 and in 2024, it appears that for most of its dispersal history, sampled OROV lineages circulated in regions of moderate viral and vector ecological suitability (~0.4–0.5), consistent across the three temporal models we estimated and the viral and vector ecological niche models (Fig. 3e,f). This was then followed by a peak in both viral and vector ecological suitability values, reaching ~0.6–0.8 associated with OROV dispersal from mid- to late-2023, coinciding with unprecedented epidemic expansion within the Amazon (Fig. 3e,f and Supplementary Fig. 1). This peak in ecological suitability was probably a result of the virus being introduced to areas of the Amazon that had a combination of favourable environments and relatively higher population densities, such as Manaus, leading to an amplification of transmission. This expansion in a highly suitable and more densely populated environment, in turn, probably facilitated the spread of the pathogen beyond its usual transmission range. This further emphasizes the need for improved surveillance in blind-spot regions where transmission suitability is high, but where cases may be under-reported or genomic data are scarce. In such regions, introductions could rapidly lead to amplified outbreaks. After 2024, ecological suitability associated with dispersal locations of sampled viral lineages became more variable, with values either remaining high (in the full OROV ecological niche model) or declining again (in the past OROV and C. paraensis ecological niche models). This variability is probably caused by the virus expanding into new areas, including those with less favourable environmental conditions (for example, lower mean annual temperatures—Fig. 2d, higher elevation and lower precipitation—Supplementary Fig. 3), despite relatively high population densities. In fact, the city of Manaus is situated within the Amazon rainforest, a unique characteristic that differs significantly from the primarily rural areas where OROV expansion has been observed.

The recent emergence and expansion of OROV into previously non-endemic regions in the Americas underscores the critical need for a deeper understanding of the factors associated with its spread. In response to the 2023–2024 outbreaks, recent studies have described the emergence of a reassortant lineage^23,24, investigated spatiotemporal movement dynamics at a smaller scale^18,19,23,24, reported severe clinical outcomes⁴ and documented biological differences in the new lineage⁴³. Some studies have investigated the ecological factors associated with the local circulation of OROV in parts of the country using epidemiological data^26,27, but not while also integrating insights from its recent dispersal history outside the Amazon Basin. Overall, large gaps remain in our understanding of the sudden large dispersal expansion beyond the Amazon⁴⁴.

Compared with previous approaches, our study provides an integrative analytical framework to examine how ecological factors may be associated with viral dispersal characteristics derived from phylogeographic reconstructions. This design enables us to explore these associations in a time-varying manner across the growth of the epidemic. We complement this with ecological niche modelling, a well-established method for assessing the spatiotemporal distribution of suitable environments for target organisms, including viruses. However, ecological niche models are observational tools with inherent limitations, especially for populations expanding into new environments that may not yet fully occupy their ecological niche. Their outputs should therefore be interpreted cautiously and ideally in combination with additional lines of evidence. Integrating phylogeographic and ecological niche modelling analyses provided complementary perspectives that allow us to address two related but distinct questions: the potential influence of environmental covariates on viral lineage dispersal dynamics and the ecological suitability associated with the risk of local viral circulation. Accordingly, for the first time we provide estimates of the strength of association of covariates using two methodological approaches, with hypothesis testing performed with genomic data, including previous comprehensive OROV genomic datasets beyond the Amazon region^18,19,23,24 and in a time-varying manner to test the timing of the association of different environmental factors with the OROV range expansion. Unlike previous investigations that focused either on spatiotemporal lineage movements within the Amazon or limited phylogenetic placements outside endemic regions^18,19,23,24,and others exploring environmental factors impacting OROV through associations between epidemiological data with land cover^26,27, our integrative approach explicitly combines phylodynamic reconstructions with ecological niche modelling to formally test associations between environmental factors and OROV dispersal. This combined framework enables a national-scale perspective and provides complementary insights into both the timing of viral expansion and the ecological conditions shaping it.

Environmental conditions play a pivotal role in shaping OROV transmission dynamics⁴⁵. Zones undergoing land cover transitions, particularly those involving deforestation and agricultural activities, have emerged as hotspots for virus spread^46,47. As these areas transition from sylvatic (forest) environments to more urbanized or agricultural landscapes, the resulting habitat changes bring vectors and reservoir hosts into closer contact with human populations, creating new opportunities for OROV transmission and these must be studied using an integrated approach⁴⁸. Overall, the results of our analyses of environmental factors were consistent across both the phylogeographic and case-based ecological niche modelling frameworks. A notable exception was broadleaf forest cover, which exhibited a negative association in the phylogeographic approach but a positive association in the ecological niche modelling analyses. This contrast could potentially reflect a dynamic in which deforestation in areas of high forest cover create ecological conditions that facilitate OROV transmission. In this context, forest cover would then exhibit opposing associations depending on whether the signal captures dispersal dynamics—where forest loss may be more influential—or local circulation, where remaining forest cover may provide suitable habitats for vectors and increase opportunities for human exposure.

Through the landscape phylogeographic and ecological niche modelling analyses conducted in this study, we propose a mechanism for the expansion of OROV, first within and then beyond the Amazon. We detect shifts in ecological suitability and differences in the environmental space associated with dispersal events preceding epidemic amplification within the Amazon. Our findings support a two-step expansion process, in which circulation within a highly suitable environment led to a rapid outbreak in the Amazon, subsequently facilitating spread beyond the usual transmission range. This pattern is supported by independent reconstructions of all three genomic segments of the virus. While previous studies have identified the Amazon region as the primary source of OROV emergence in Brazil^23,24, re-emergence or introduction of the reassortant lineage in the vicinity of Manaus probably exposed a large human population and contributed to outbreak amplification. Higher suitability in this zone may reflect a combination of favourable environmental conditions and relatively denser human populations within the Amazon. Our results further suggest that circulation of OROV in non-endemic regions of Brazil may be self-limiting, particularly in less favourable environments.

More broadly, our study highlights the value of integrating phylodynamics and ecological niche modelling to investigate eco-epidemiological dynamics of re-emerging arboviruses and urban amplification. Using complementary analytical approaches and a comprehensive genomic dataset, we identify key environmental factors associated with OROV transmission. Both approaches implicate higher human population density and temperature, as well as banana and cocoa cultivation, as factors associated with OROV dispersal and circulation, consistent with previous findings²⁶. These insights can inform public health planning and mitigation measures, including targeted vector control around agricultural areas near urban centres. While not confirmatory, our findings indicate that OROV expansion was associated with increasingly favourable and well-connected environments. Genomic reassortment may have contributed to enhanced transmission efficiency or immune escape, as suggested by prior hypotheses. The concentration of high-suitability areas in coastal regions of Brazil is of particular concern, given dense populations and the cocirculation of several arboviruses. Repeated introductions into these areas could lead to further amplification, especially in immunologically naive populations or in the presence of variants capable of immune evasion⁴³. Transport infrastructure probably also shapes OROV spread, particularly in the Amazon region where roads and waterways act as major mobility corridors. Several suitable regions, especially in the northeastern and central-west areas, remain surveillance blind spots, underscoring the need to prioritize active surveillance for early detection and outbreak prevention.

This study needs to be interpreted in light of several limitations. Our landscape phylogeographic approach is affected by the pattern of the sampling effort, as approximately half the node locations are sequenced tips, making them prone to sampling biases³⁹. This means that this approach cannot currently test the true drivers of transmission but rather tests the strength of associations of dispersal environments. Additionally, we lack a species distribution model for vertebrate host species (for example, sloths) that typically compose the natural reservoir of OROV, which limits our ability to fully understand host–vector-specific environmental suitability. To mitigate this bias, we test associations against a null dispersal model generated through simulations and only consider factors with the strongest BF support (>20). Also of note was that this study was designed to only examine the environmental factors associated with dispersal expansion in Brazil, such that the potentially critical role of reassortants in viral adaptation and fitness was not addressed here. Our phylogeographic analyses focused solely on the reassortant lineage that emerged in Brazil. To elucidate the role of these evolutionary processes in the ability of the virus to adapt to new niches, hosts or vectors, our approach could be extended to analyse evolutionarily distinct lineages of OROV that previously circulated in Brazil or neighbouring regions. Vectors could also have played a central role in the recent range expansion of OROV, with the primary vector, C. paraensis, potentially adapting to new landscapes. Yet, the limited data on this species and the lack of direct evidence of such adaptation highlight the need for a more systematic surveillance to better characterize the distribution of C. paraensis. Moreover, ecological niche models are constrained by the availability and representativeness of georeferenced data and should therefore be interpreted cautiously. The interaction between the virus, its vectors and its reservoir hosts is crucial for both the colonization of new transmission zones and the maintenance of transmission in established areas. Finally, while deforestation has been discussed as a potential driver of the 2023–2024 OROV outbreak, in our analyses we did not find strong support for an association between forest loss and viral lineage dispersal locations. However, this absence of supported association should not be interpreted as evidence that forest loss is unimportant for OROV transmission. Rather, its effect may already be captured through correlated land-use variables in our models, particularly agricultural expansion. This may also reflect time-lags between forest cover loss and downstream impacts or limitations in the spatiotemporal resolution of the available data to fully capture these dynamics.

The geographic expansion of OROV into new regions highlights the pressing need for integrating environmental monitoring into public health frameworks. To effectively predict and mitigate the risks posed by OROV and other arboviruses, surveillance efforts must account for the complex interaction between environmental changes, vector ecology and human behaviour. As OROV continues to adapt to new ecological niches, potentially driven by combinations of genetic evolution and both natural and anthropogenic factors, a deep understanding of these dynamics will be essential for developing targeted intervention strategies to control its spread and minimize its public health impact.

OROV genomic data

Complete genome sequences of the S, M and L segments of the OROV, obtained from the first extra-Amazon OROV cases reported in the states of Bahia (northeast Brazil), Minas Gerais (southeast Brazil), Mato Grosso (midwest Brazil) and Paraná (south Brazil), were analysed together with corresponding segments from recently published full-length OROV genomes belonging to the Brazilian 2022–2024 sublineage²⁴. The OROV sequences used correspond to GenBank accession IDs PQ168520–PQ247806 and PP153945–PQ065491. For each genomic segment (n = 545), several sequence alignments were generated using MAFFT^49,50 and subsequently manually inspected and curated in AliView⁵¹ to remove sequencing artefacts and alignment errors. To ensure appropriate quality control for downstream evolutionary analyses, genomic regions identified by RDP5⁵² as having probably arisen through recombination or reassortment were systematically excluded from the analyses. Specifically, these regions were masked by replacing the corresponding alignment positions with gap characters (‘–’), resulting in genome alignments that explicitly exclude recombinant or reassortant signals, following previously established and validated approaches⁵². Sequences exhibiting extensive recombination or reassortment across most of the genome were excluded from further analyses (n = 43 for segment S and n = 1 for segment L). The resulting masked alignments were used for all downstream evolutionary and phylodynamic analyses.

Disease occurrence data

Disease occurrence data were compiled from several sources. Epidemiological data on OROV cases were retrieved from the Brazilian Ministry of Health, accessible at https://www.gov.br/saude/pt-br/assuntos/saude-de-a-a-z/o/oropouche. The dataset includes confirmed case reports from all Brazilian states where OROV cases have been notified, including but not limited to Acre, Alagoas, Amazonas, Bahia, Ceará, Minas Gerais, Pará, Rio de Janeiro and São Paulo. The data cover the years 2023 and 2024, organized by epidemiological week and include key variables such as the municipality, state, year of occurrence and the corresponding epidemiological week of reporting. These data were geocoded at the municipality level and occurrence deduplicated by month. Additional OROV occurrence data were gathered and geocoded from all records on the Global Biodiversity Information Facility (years 1957 to present) and GenBank (years 2015 to present) databases. After deduplicating all occurrence records, we obtained a total of 450 unique sampling location points covering the years 1957 to 2024, with most (~85%) sampled in 2023–2024.

Geospatial data

We tested several environmental factors both as associations with viral dispersal locations and as covariates in our ecological niche model. These factors included human population density, main land cover and climatic variables within the study area (Brazil). Each environmental factor was described by a raster that defines its spatial heterogeneity. Supplementary Table 3 details the source and resolution of each original raster file. In particular, the C. paraensis distribution map was obtained from previous work which used 78 cleaned occurrence points from Brazil to fit a random forest (RF) model (TSS = 0.687). Key predictors of suitability included banana- and sugarcane-harvested area, humidity, soil moisture and cassava-harvested area^35,53. Each raster was cropped to match our study area (Brazil) by using Brazil shapefiles from the “rnaturalearth” package in R.

Phylogeographic reconstruction and lineage dispersal statistics

To model the spatiotemporal spread of OROV using spatially explicit phylogeographic reconstruction using the continuous diffusion model implemented in the software package BEAST 1.10⁵⁴, we used three different data subsets containing the 2022–2024 OROV sequences from various regions of Brazil of the S (n = 501), M (n = 545) and L (n = 544) segments. Before conducting the phylogeographic analyses, we assessed the strength of the molecular clock signal in each data subset using the root-to-tip regression method available in TempEst v1.5.3⁵⁵. Preliminary BEAST reconstructions revealed an outgroup of sequences which diverged from the main clade ~60 years ago for each segment; these were subsequently discarded from our analyses. Temporal signal was evaluated for the resulting datasets yielding to correlation coefficients close to or above 0.5 (S, 0.4972; M, 0.635; L, 0.5123), which corresponds to a weak yet substantial temporal signal. We reconstructed the spread of OROV lineages within Brazil by using a flexible relaxed random walk (RRW) diffusion model⁵⁶, which accommodates branch-specific variation in dispersal rates, with a Cauchy distribution and a jitter window size of 0.01 (ref. ⁵⁷). The latitude and longitude coordinates of each sample were used in this analysis. Markov Chain Monte Carlo analyses were run in BEAST v1.10.4, with chains of up to one billion iterations each, sampling every 100,000 steps in the chain. The chains were stopped when convergence was reached following the removal of burn-in states. Convergence of each run was assessed using Tracer v.1.7.1, ensuring that the effective sample size for all relevant model parameters was >200 (ref. ⁵⁸). Maximum clade credibility trees were retrieved and annotated using TreeAnnotator after discarding burn-in samples, the number of which was also determined in Tracer. Finally, the R package “seraphim”⁵⁹ was used to extract and map the spatiotemporal information embedded in the posterior trees. We further used “seraphim” to estimate three dispersal statistics from these movement vectors for each segment: maximal wavefront distances, weighted diffusion coefficients⁶⁰, measuring the dispersal capacity of viral lineages and an IBD signal measured as the Pearson correlation between the patristic and log-transformed geographic distances computed for each pair of tip nodes⁶¹.

Landscape phylogeographic analyses

To investigate the association between environmental conditions (Supplementary Fig. 2) and dispersal locations of inferred OROV lineages, we first conducted a visual exploration and then formally tested these associations using a landscape phylogeographic approach³⁹. For the visual exploration, we reported the environmental values explored by phylogenetic branches using the “spreadValues” function implemented in the R package “seraphim”⁵⁷. Specifically, this function extracts and averages the environmental values crossed by the phylogenetic branch segments occurring during each successive time slice; and this individually for several trees sampled from the post-burn-in posterior distribution of trees inferred by each segment-specific continuous phylogeographic analysis. We here considered 100 posterior trees from each phylogeographic analysis, 200 time slices each spanning ~2 weeks from early July 2015 to early May 2024, as well as a sliding window of ~1 month.

To assess the tendency of inferred viral lineages to preferentially circulate within or avoid circulating in specific environmental conditions, we then compared the distribution of mean environmental values extracted at node positions in inferred trees (E_estimated) with those extracted at node positions in trees whose dispersal history had been resimulated under a null dispersal model (E_simulated). To generate such a null dispersal model, an RRW diffusion process was simulated along each tree topology used for the phylogeographic analyses. These RRW simulations were performed using the “simulatorRRW1” function of the R package “seraphim” from the sampled precision matrix parameters estimated by the phylogeographic analyses. Simulations conducted under this null dispersal model⁶² aim to generate a set of lineage dispersal histories (phylogeographic reconstructions) (1) based on the same sampling effort (the same number of samples), (2) maintaining the tree topology and branch lengths (the evolutionary relationships among samples), (3) starting from the same epidemic origin (the location inferred for the most ancestral node of each posterior tree), (4) avoiding non-accessible areas (for example, marine areas) and (5) constrained within the study area (here delimited by the minimum convex hull polygon drawn around all internal and tip node locations retrieved from the set of posterior trees inferred by the continuous phylogeographic analysis). Specifically, the aim of those simulations is therefore to generate a null dispersal model imitating the actual dispersal history but in which the environmental factors have no determined impact on the dispersal location of viral lineages.

From these simulations, values at node positions (E_simulated) constitute the distribution of mean environmental values explored under a dispersal scenario that is not impacted by any underlying environmental condition. For each environmental factor and segment-specific phylogeographic reconstruction, we then compared the distribution of E_estimated values computed from posterior trees with the distribution of E_simulated values retrieved from the same tree topologies along which an RRW diffusion process had been resimulated. Specifically, we approximated a BF support equal to (p_e/(1 − p_e))/(0.5/(1 − 0.5)). To test if viral lineages tended to avoid circulating within a particular environmental factor e, p_e was defined as the frequency at which E_estimated < E_simulated; and to test if viral lineages tended to preferentially circulate within a particular environmental factor e, p_e was defined as the frequency at which E_simulated < E_estimated. Following the scale of interpretation in ref. ³⁸, we here highlight BF values >20 considered as strong supports.

Ecological niche modelling

Ecological niche models are built using a variety of statistical methods, each varying in complexity and underlying assumptions about the interaction between species occurrences and environmental factors⁶³. Recent studies have shown that disparities among different model structures can be very large, making model selection difficult⁶⁴. It is important to note that ecological niche models may achieve comparable performance metrics yet produce very different suitability maps, as also highlighted in ref. ⁶⁵. An alternative is to use an ensemble of models to avoid selecting one single best model but instead to use a group of methods for inference. In other words, the presence of a species might be well classified by some models and misclassified by others, such that making use of an ensemble model can reduce the predictive uncertainty of a single model by combining predictions⁶⁶. In this study, we applied this approach to investigate the distribution of OROV transmission by creating an ecological suitability map based on the occurrence of OROV disease in Brazil and relevant environmental variables (section on Geospatial Data above). We defined ecological suitability as the output of ecological niche modelling, estimating the likelihood that a given location provides environmental conditions favourable for the potential local circulation of OROV. For this analysis, we aggregate all raster maps to the lowest resolution available (~27 km²). We used an ensemble of seven statistical, machine learning and envelope models: generalized linear model, generalized additive model, boosted regression trees, RF, classification tree analysis, surface range envelope and maximum entropy model. For RF, we built a conservative RF (500 trees with a maximum depth of 30) to avoid overfitting; detailed parameters for all models are given in Supplementary Table 6. To assess the potential expansion of OROV ecological niche in Brazil during late 2023 and 2024, we computed and compared ecological niche models using occurrence data from different time periods. Specifically, we created three ensemble models using pre-mid-2023 (89 occurrences), pre-2024 (133 occurrences) and the full model with all available occurrence records (450 occurrences).

To assess the performance of the models we use block cross validation⁶⁷. This is a spatially explicit method used to assess model performance by dividing the study area into geographic blocks. Instead of randomly splitting data, this technique ensures that training and testing data are spatially independent, reducing the risk of spatial autocorrelation. TSS and receiver-operating characteristic curve (AUC) are then used to evaluate the predictive performance of the models based on the test (validation) dataset. TSS is equivalent to sensitivity + specificity − 1 and ranges from −1 to 1; value of 1 indicates perfect classification, 0 means the model is no better than random guessing and negative TSS indicates the model performs worse than random guessing. AUC ranges from 0 to 1; AUC of 1 indicates perfect model performance, 0.5 indicates no discrimination (that is, the model is no better than random guessing) and <0.5 indicates the model performs worse than random guessing. We only retained models with a TSS score of >0.7 to build the ensemble model (Supplementary Fig. 4d). The mean probabilities from each model were then computed and we weighted the predictions of each model according to its performance during training, giving more weight to better-performing models. The weights ensure that higher-quality models contribute more to the final ensemble prediction. The different resulting ‘suitability indexes’ are then combined to get a single value per site. Once the ensemble predictions are generated, the ensemble model itself is evaluated using the same metrics applied to the individual models; TSS and AUC.

Disease presence points used as input represent OROV circulation occurrence from 450 unique sampling locations in Brazil (molecular testing and sequencing records) from the years 1957 to 2024. Occurrence points with available collection dates were matched to corresponding climatic variables by month. Specifically, temperature and precipitation data for each point were extracted from the monthly climate layers matching the collection month. The sampling of pseudo-absences was done at a 1:1 ratio with presence points⁶⁸ and informed by a distribution of presence points and a human population density kernel density estimate. Our aim was to sample absences in proportion to the rate of presence points while giving higher priority to areas with greater population density, ensuring more focused sampling in regions where human populations are denser. This approach helps us to account for disease testing biases in more urbanized areas in the pseudo-absence distribution. Pseudo-absences were also selected within a perimeter of 50–300 km around presence points, ensuring that absence points are neither too close to presence points (to avoid the same niche) or too far (to promote localized sampling strategy). This approach was tested against a random sampling strategy (Supplementary Table 4) for pseudo-absences and after a sensitivity analysis testing the impact of maximum radius distance around presence points (varying from 150 km to 450 km) on model performance (Supplementary Table 5). The minimum distance was kept at 50 km to ensure that pseudo-absence points do not coincide spatially with presence point locations (pixels)^68,69. This was tested on the full model only.

To determine the independent contributions of each variable to our suitability prediction, we further calculated their RI). Because the RF model demonstrated the highest predictive performance, RI and response curves were derived from this model, whereas the ensemble model was used to build the suitability maps. In the case of RF models, RI is computed by assessing how frequently a variable is selected for splitting at tree nodes, weighted by the squared improvement in model performance resulting from each split and averaged across all trees⁷⁰. Higher RI values indicate greater relative contribution of that variable to the predictive performance of the model. During variable selection, predictors with low contribution or high collinearity with other variables were excluded to avoid redundancy. Forest loss was among the variables removed during this process, as it was strongly correlated with forest cover and did not improve model performance. We also produced response curves to visualize the effect of each variable on suitability predictions in the RF models. These response curves allow us to observe how changes in a single variable influence the predicted outcome, while other variables are held constant (at their mean). By examining these relationships, we gain insights into how each variable individually contributes to the overall predictions of the model.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.