Introduction

Since 2020, an increase in both H5Nx and H5N1 cases has been observed (Fig. 1). A variant of the Gs/Gd H5N1 viruses belonging to the H5 clade 2.3.4.4b has led to an unprecedented number of deaths in wild birds and poultry in many countries in Africa, Asia, and Europe (Huang et al., 2023). In 2021, the virus spread to North America, and in 2022, to Central and South America (FAO, WHO and WOAH situation analysis, July 12, 2023) spilling over into poultry farms and infecting an alarming number of wild terrestrial, marine and domestic mammals (Elsmo et al., 2023; Leguia et al., 2023; Neumann and Kawaoka, 2024; Peacock et al., 2024; Plaza et al., 2024); marking an unprecedented expansion in the geography and impact of HPAI. In 2022, 67 countries across five continents reported H5N1 HPAI outbreaks in poultry and wild birds to WOAH, resulting in over 131 million domestic poultry deaths or cullings (WOAH, 2023). In 2023, another 14 countries, mostly in the Americas, reported outbreaks. Several mass death events in wild birds were caused by influenza A(H5N1) clade 2.3.4.4b viruses (Klaassen and Wille, 2023; WOAH, 2023).

Figure 1:

In that epidemiological context, understanding the environmental factors associated with the risk of HPAI circulation remains crucial. Extensive exploration of risk factors for HPAI presence, spread, and persistence has been conducted in various countries and regions. These studies have identified domestic waterfowl, several anthropogenic variables (human population density, distance to roads) and indicators of water presence as important factors for risk of H5N1 local circulation (Gilbert and Pfeiffer, 2012). More recent work on global suitability for HPAI have pinpointed that host-related variables such as poultry density are the strongest contributors to HPAI persistence (Dhingra et al., 2016; Gierak and Śmietanka, 2021; Martin et al., 2011).

In addition, poultry intensification, international poultry trade, live bird markets, and wild bird migratory routes have been recognised as playing a key role in the transmission and spread of HPAI (Vandegrift et al., 2010). Lastly, evidence of climate change impacting the dynamics of HPAI has also been discussed (Prosser et al., 2023). For example, changes in rainfall alter the distribution, abundance, and quality of wetlands and can impact waterfowl populations (Forcey et al., 2007; Gilbert et al., 2008; Vandegrift et al., 2010). However, recent outbreaks raise questions on the ability of previous models to assess HPAI H5 ecological suitability.

Understanding the interplay between environmental variables and avian influenza becomes increasingly crucial to have an in-depth understanding of the risk factors that govern AIVs circulation and spread. Following a previous study dedicated to the ecological niche of H5NX and H5N1 (Dhingra et al., 2016), we aimed to address the following questions: (i) what factors can explain the observed increase in H5Nx and H5N1 cases since 2020? (ii) Do we observe a change or an extension of the ecological niche beyond its traditional spatial coverage? (iii) are we able to predict the areas at risk based on ecological niche models trained on occurrence data before 2020? To this end, we used HPAI H5N1 and H5Nx occurrence data collected before 2020 to train ecological niche models and tested whether these models could accurately predict areas ecologically suitable for HPAI circulation post-2020. Specifically, we employed a boosted regression trees (BRT) approach to estimate the HPAI ecological suitability given local environmental conditions. We also computed the diversity indices of wild bird species involved in HPAI occurrences, helping to understand how shifts in species diversity might correspond with the environmental factors considered in our ecological niche modelling analyses. Our results provide a better understanding of HPAI dynamics by identifying key environmental factors driving the increase in H5Nx and H5N1 cases in poultry and wild birds, investigating potential shifts in their ecological niches, and improving the prediction of at-risk areas. This knowledge might aid policymakers, epidemiologists, and wildlife conservationists in better targeting surveillance and control measures.

Results

In this study, we apply an ecological niche modelling approach using the boosted regression trees (BRT) method to evaluate the suitability of environments for H5Nx and H5N1 HPAI viruses in wild and domestic birds across two periods: 2015-2020 and 2020-2022. Specifically, we trained distinct ecological niche models for each combination of species (wild or domestic birds), virus strain (H5Nx or H5N1), and study period. In addition to evaluating each model separately, models trained on data from the 2015-2020 period were tested to predict occurrences in the 2020-2022 period, allowing us to evaluate the potential of pre-2020 models to predict a more recent distribution of H5N1 and H5Nx cases.

We assess the influence of a broad range of spatial predictors — including anthropogenic, topographical, land cover, eco-climatic, and poultry density — on virus occurrences (Fig. S1). We use occurrence records from the EMPRES-i database, and sample pseudo-absence points to account for the lack of virus detection data in some areas. These pseudo-absence points were distributed based on human population density, with more pseudo-absence points sampled in areas of higher population to reflect the greater surveillance efforts in those regions (Fig. S2).

The accuracy of the models was validated through three types of cross-validation: (i) a standard cross-validation with a random and stratified divide between training and validation sets, (ii) a spatial cross-validation based on the approach used by Dhingra and colleagues (2016) and that consists in clustering presence and pseudo-absence points into folds using reference presence points, and (iii) a second spatial cross-validation based on a blocks generation technique (Valavi et al., 2019). Models using spatial cross-validation techniques are less affected by spatial autocorrelation, as shown by higher spatial sorting bias (SSB) values (Fig. S3). Of the seven datasets evaluated, four show better performance with the reference points-based spatial cross-validation method, suggesting it is less impacted by spatial autocorrelation compared to the block generation technique. Consequently, we have proceeded with the reference points approach for consistency with Dhingra et al. (2016).

Using pre-2020 data, our models demonstrate a relatively robust capability to predict the distribution of post-2020 occurrence data: the Area Under the Curve (AUC) values for ecological niche models trained occurrence data in domestic birds before 2020 range between 0.74 and 0.77 when evaluated with occurrence data in domestic birds after 2020, indicating good predictability. In comparison, ecological niche models trained on occurrence data in domestic birds after 2020 display slightly higher predictive performances when evaluated with the same data, with AUC values ranging between 0.78 and 0.83 (Table S1). These findings underscore the pre-2020 models’ effectiveness in forecasting the recent geographic distribution of ecological suitability for H5Nx and H5N1 occurrences.

Previous literature reviews (Gilbert and Pfeiffer, 2012; Dhingra et al., 2016) have summarised the predictor variables commonly linked to occurrences of HPAI. As detailed in the Materials and Methods section and in the table in Supplementary Information Resources S1, we explored four sets of environmental variables considered by Dhingra and colleagues (Dhingra et al. 2016): host variables (set 1), land cover variables (set 2), eco-climatic variables (set 3), and a risk-based selection of variables performed by Dhingra and colleagues (set 4). Set 3 shows the lowest predictive performance for both domestic and wild bird cases during the two periods and, in contrast, sets 2 and 4 demonstrate the highest predictive accuracy for wild and domestic birds, respectively (Fig. S3). Set 4 includes a combination of host variables with cultivated and managed vegetation, open water areas, distance to water, and the annual mean of land surface temperature. For domestic birds, set 1, which includes host variables, shows high predictive performance before 2020, but its performance decreases after 2020 in favour of set 4. The ecological niche models trained on wild and domestic bird cases, respectively with the sets 2 and 4, are associated with AUC values exceeding 0.77, indicating relatively good predictions of H5 occurrences with these sets of environmental variables.

Figure 2 displays the response curves of the most important variables with their respective relative influence (RI) (Table S2). For the H5N1 and H5NX ecological niche models for domestic birds, densities of intensive chicken population, duck population, and human population emerge as significant predictors, each with RI values exceeding 5%. Notably, there is a sharp increase in the RI for H5N1 in areas with high densities of intensive chicken populations — from 8.5% to 30.4% since 2020. Similarly, the RI of cultivated and managed vegetation has doubled for both strains post-2020. Moreover, the response curves associated with these predictors show a positive correlation, indicating that higher values of these predictors are linked with an increased likelihood of HPAI occurrences given local environmental conditions. This result indicates a trend towards increased HPAI susceptibility in environments characterised by intensive agricultural and vegetation management practices.

Figure 2:

In contrast, eco-climatic factors such as land surface temperature and precipitation (set 3) showed only moderate influence. We also observe a decrease in the importance of duck population density for both subtypes after 2020, possibly due to the increasing diversity of bird species involved in the transmission dynamics of H5N1 and H5NX (Fig. S4). These findings highlight the crucial role of anthropogenic and host-related variables in accurately predicting HPAI occurrences.

In the ecological niche models trained for wild bird cases, urban and built-up areas are associated with H5N1 and H5Nx outbreaks (Fig. 2), and have the highest relative importance (RI) prior to 2020. Specifically, the RI values were 54.5% for the period before 2020 and decreased to 39.3% for predictions extending beyond 2020. This decline in RI values may indicate a reduced bias in observation data; typically, dead wild birds are more frequently found near human-populated areas, which could influence earlier data sets. After 2020, we observe an increase in the importance of habitats, such as deciduous broadleaf trees (6.1%), mixed and other tree regions (11.4%), and herbaceous vegetation for H5Nx (9.5%). However, as illustrated by the response curves (Fig. 2), the presence of H5N1 and H5Nx in wild birds declines in areas with greater tree cover, while it increases in regions dominated by herbaceous vegetation. Open water areas consistently show high RI values across all time periods and virus strains, particularly for H5Nx before 2020 (25.5%) and H5N1 after 2020 (22.0%), highlighting their significant role in the ecological models of HPAI outbreaks.

The risk maps generated show the predicted ecological suitability for H5N1 and H5NX avian influenza viruses across two time periods, 2015-2020 and 2020-2022 (Fig. 3), and are also accessible on the following link for a dynamic visualisation of the results: https://mood-platform.avia-gis.com/core. These maps reveal both broad regional risks and pinpointing specific locales of heightened susceptibility, especially in Europe and Asia (India and China). Our models highlight extensive spread of H5N1 and H5Nx HPAI among domestic birds in South Korea, Japan, Singapore, Malaysia, Vietnam, Cambodia, Thailand, and the Philippines, as well as in European countries such as the United Kingdom, France, the Netherlands, Germany, Italy, Ukraine, and Poland. Additionally, several areas in African nations, including Nigeria and South Africa, are identified as suitable environments for both H5N1 and H5Nx occurrence in scenarios before and after 2020. Regions in North America and South America — including Bolivia, Brazil, Colombia, Ecuador, Chile, Peru, and Venezuela — demonstrate ecological suitability for H5Nx and H5N1, with the period after 2020 showing a marked increase in the risk of local H5N1 circulation.

Figure 3:

Our findings remain overall consistent with the ones previously reported by Dhingra and colleagues (Dhingra et al., 2016), who used data from January 2004 to March 2015 for domestic poultry, even if we can note some differences: North America, West Africa, eastern Europe and Bangladesh have a higher risk of H5 infection before 2016, while our maps mainly highlight regions in China, South-East Asia, and Europe (Fig. S5). In India, both analyses consistently pinpoint strong risk areas for H5N1 and H5Nx, though our data suggests these areas have higher ecological suitability. Similar to results reported by Dhingra and colleagues, we observe an increase in the ecological suitability estimated for local H5N1 circulation in South America’s domestic bird populations post-2020. Finally, in North America, while Dhingra and colleagues predicted higher risk areas for H5Nx, our findings highlight similar regions but with lower ecological suitability.

For H5Nx and H5N1 cases in wild birds, the estimated ecological niches are mostly related to environmental factors such as open water and distance to water, which result in ecological suitability hotspots estimated near coasts and main rivers. For both H5Nx and H5N1, however, isolated areas on the risk map should be interpreted with caution. These isolated areas may result from sparse data, model limitations, or local environmental conditions that may not accurately reflect true ecological suitability. The risk maps reveal that high-risk areas have expanded after 2020. North America appears to be more susceptible to H5Nx in particular near the Great Lakes region. This expansion is evident in Russia, South America, and North America aligning with major bird migration routes. Notably, while H5Nx has never been reported in Australia, the models indicate potential ecological suitability there as well. Additionally, regions in China, India, and Europe also display expanded areas of ecological suitability.

In Table 1, we report the estimation of avian species diversity indices for species involved in HPAI transmission for the pre-2020 and post-2020 periods. We observe variations between these two periods both in the overall bird population and within specific wild bird species groups, including marine species. For all birds, the Shannon index increased from 4.23 before 2020 to 5.03 after 2020, indicating a more diverse infected bird population in the latter period. The Simpson index also rose from 0.97 to 0.98, suggesting a slightly higher concentration of certain species post-2020. Sea birds exhibited a similar upward trend in diversity post-2020, with the Shannon index increasing from 2.19 to 2.26. This rise in Shannon index for sea birds suggests a broadening of species diversity within this category from 2020 onwards.

Table 1.

When considering birds affected by the H5N1 strain of HPAI, there is a notable increase in the Shannon index from 2.00 before 2020 to 2.95 after 2020, accompanied by a rise in the Simpson index from 0.83 to 0.91. These increases could be indicative of greater diversity and a more even distribution of species affected by H5N1 after 2020. In contrast, birds not affected by the H5N1 strain showed a slight decrease in diversity post-2020. The Shannon index decreases from 3.42 to 3.26, and the Simpson index decreases from 0.95 to 0.92. Given that HPAI particularly affected sea birds >2020, we further explore the yearly distribution of sea bird families. As depicted in Fig. S4, sea birds belonging to the families Laridae, Sulidae, Ardeidae, Pelecanidae, Spheniscidae, Alcidae, Phalacrocoracidae were more prevalent >2020 in comparison to <2020. Especially, the families Laridae and Sulidae have a special increased prevalence from 2021 onwards.

Discussion

Using an ecological niche modelling approach, our study presents an in-depth analysis of the ecological suitability for the risk of local circulation of H5Nx and H5N1 HPAI viruses across different temporal and ecological contexts. By conducting boosted regression tree (BRT) analyses and incorporating a diverse set of spatial predictors, our study provides insights into the environmental and anthropogenic factors influencing the distribution of these viruses in both wild and domestic bird populations. For domestic birds, our results concur with the those previously reported by Dhingra and colleagues (Dhingra et al., 2016), indicating that ecological niche models incorporating host variables (chicken and duck population densities) outperform those based on land-use or eco-climatic factors in terms of predictive performance. Contrariwise, in our analysis of H5N1 and H5Nx HPAI in wild birds, the most significant predictors were anthropogenic factors, notably urban areas and distance to open water. However, observation biases might affect these findings since outbreaks in wild birds are commonly reported near human populations.

In our analyses, intensive chicken population density emerges as a significant predictor, underscoring the influence of poultry farming practices on HPAI spread (Dhingra et al., 2016; Vincenti-Gonzalez et al., 2024) and revealing a complex interplay between ecological niches and the anthropogenic landscape. Similarly to the results of Dhingra and colleagues (Dhingra et al., 2016), land-use and climatic predictors do not play an important role in the niche ecological models, even for wild birds. Additionally, we found that while duck population density is a significant predictor for HPAI occurrence in domestic birds before 2020, its importance has decreased post-2020, particularly within the H5Nx ecological niche model. This observation is particularly noteworthy, as ducks have historically been closely linked to areas of persistence and evolution of H5N1 HPAI (Gilbert and Pfeiffer, 2012). However, our findings, in line with earlier studies such as the one of Dhingra and colleagues (Dhingra et al., 2016), highlight a significant shift in the epidemiological dynamics of HPAI. Specifically, while the association between the H5Nx clade 2.3.4.4 and ducks has for instance been documented in South Korea (Hill et al., 2015), it is becoming increasingly evident that HPAI is now more strongly associated with high-density chicken farming and anthropogenic factors. This shift suggests that the virus could now be circulating more among farms (farm to farm transmission) (FAO, 2023), with fewer introductions by wild birds — particularly wild ducks, which have previously acted as “Trojan horses” for H5N1 HPAI (Kim et al., 2014). Consequently, we could be witnessing a shift from the traditional duck/rice ecosystems towards a scenario where intensive chicken farming and a wider diversity of wild birds play a more prominent role in the transmission dynamics of HPAI. Therefore, monitoring areas with high intensive chicken densities and regular surveillance of wild birds remains crucial for the early detection and management of HPAI outbreaks.

A difference between the ecological niche models previously trained by Dhingra and colleagues (Dhingra et al., 2016) and our models is the pronounced increase in ecological suitability for H5N1 in domestic bird populations estimated in India and China. More generally, our models highlight the extensive spread of H5N1 and H5Nx HPAI among domestic birds in Asia previously identified as “reassortment sink areas” (Dhingra et al., 2018). Reassortment refers to the process where influenza viruses exchange genetic material, potentially leading to new, more virulent strains. The extent of virus suitability in Asia highlights a potential change in the dynamics of virus spread, possibly driven by intensified poultry farming practices and increased human-wildlife interactions in these regions (Bahl et al., 2016; Gilbert et al., 2017).

The recent outbreaks of infection in cattle within the United States have been partially delineated by our HPAI risk mapping approach, highlighting major hotspots predominantly in the Great Lakes region, including Michigan. From March to July 2024, approximately 26 herds in Michigan were affected (U.S. Department of Agriculture, 2024). Additionally, our risk maps identify localised hotspots around urban centres in Colorado and Texas, where herds have also been infected. These outbreaks appear to stem from a single introduction event in Texas by wild birds (Worobey et al., 2024a, 2024b). This trend may be indicative of an increasing risk of future introductions from wild birds, as we observed a significant rise in the diversity of infected bird species after 2020, particularly among sea birds.

The latest strain of H5N1 (H5 clade 2.3.4.4b), believed to have emerged in 2014, shows adaptations for infecting a wide range of wild birds and mammals (Caliendo et al., 2022; Graziosi et al., 2024). By 2020, its spread among wild birds was three times faster than in farmed poultry due to mutations enabling the virus to infect diverse species (CDC, 2024; Xie et al., 2023). Unlike previous epidemic waves in Europe, outbreaks in the last years have also continued through summer, confirming that the virus’s current circulation is no longer seasonal but nesting in the wild bird population (Sciensano, 2024). Therefore, the virus is able to infect a wider range of species than previous forms, as it circulates all seasons (Sacristán et al., 2024). The densely packed colonies of sea birds have also facilitated the rapid virus spread, resulting in high mortality in these sea bird colonies (Boulinier, 2023; Bregnballe et al., 2024). For instance, the Netherlands lost up to 80% of its Sandwich Terns in weeks, and the UK saw significant outbreaks, especially on Scottish islands (Knief et al., 2024). In South America, similar trends have been observed, with notable die-offs in seabird populations along coastal regions, particularly impacting species such as penguins, pelicans and cormorants (Leguia et al., 2023). Since then, in early 2023, the Peruvian outbreak had spread to marine mammals, particularly affecting the South American sea lion, which also began to experience a mass die-off (Campagna et al., 2024; Leguia et al., 2023).

The increase in the number of wild bird species being infected by HPAI, particularly H5N1, could be attributed to several factors such as viral evolution, increased interaction between domestic and wild birds, climate and land-use changes, among others. Regarding land-use changes, large parts of Europe have become increasingly fragmented due to the expansion of urban and transport infrastructure (European Environment Agency (EEA), 2022). This habitat loss and fragmentation can significantly impact water availability and trophic resources (Santos-Tovar et al., 2024), which could in turn aid the spread and transmission of avian influenza viruses due to relocation of birds and outbreaks (Yin et al., 2022). Climate change has also been reported as a factor that could affect wild bird distribution in different ways (Gilbert et al., 2008). As global climate conditions change, avian migratory patterns and routes are also changing (Hitch and Leberg, 2007; Van Doren, 2022). Furthermore, higher temperatures and extreme weather have resulted in large-scale population shifts in a range of temperate species (Dezfuli et al., 2022). Overall, these elements therefore raise the question of the role of climate and land-use changes in the recent shifting epidemiology of HPAI in wild birds.

The increasing diversity of infected bird species is part of a broader resurgence of H5N1 globally. Phylogenetic analysis from Europe revealed significant genetic diversity in HPAI-H5N1 viruses since October 2020, likely due to multiple reassortment events with LPAI viruses. Notable introductions from Russia to Europe were observed through autumn migrations of wild birds (ECDC, 2021). Additionally, the transatlantic transport of HPAI H5N1, documented in Canada in December 2021, suggests the virus can spread across continents, a route previously seen only with LPAI viruses (Caliendo et al., 2022; Dusek et al., 2014; Huang et al., 2014). In November 2022, H5N1 (clade 2.3.4.4b) was introduced into South America, presumably through the movements of migratory wild birds travelling south during the boreal winter. The virus spread rapidly and caused significant mortality in wild birds and marine mammals across multiple countries (Leguia et al., 2023; (Banyard et al., 2023; Campagna et al., 2024; Gamarra-Toledo et al., 2023; WOAH, 2024).

Despite the robustness of our ecological niche models, some limitations must be acknowledged. Firstly, the accuracy of our models is inherently dependent on the quality and completeness of the occurrence data used, which is influenced by the different reporting systems in place in different countries. Secondly, the definitions of key explanatory variables, such as intensive versus extensive production, may vary across regions, adding complexity to global predictions. To address these limitations, future work could focus on developing country-level or regional models that incorporate localised data and context-specific definitions. Thirdly, underreporting and biases in outbreak reports, particularly in wild bird populations, can skew model predictions and underestimate areas at risk. Finally, another limitation of our model for wild bird populations is the lack of incorporation of migratory patterns and behaviours. Migratory birds have complex movement patterns influenced by seasonal changes, food availability, and habitat conditions. Our models, using static spatial predictors, do not capture these dynamic routes and stopover sites, leading to potential inaccuracies in predicting areas at risk for HPAI spread. In future work, the integration of dynamic migratory data and environmental changes into our ecological niche models could provide a more comprehensive and accurate risk assessment for HPAI spread among wild bird populations.

Despite these limitations, our study provides significant contributions to the understanding of H5N1 and H5Nx HPAI virus dynamics. By integrating a wide range of spatial predictors and using ecological niche modelling techniques, we provide a quantitative analysis on the factors driving virus distribution. Overall, our study underscores the importance of continuous surveillance and the need for adaptive management strategies to address the evolving threat of HPAI viruses globally.

Materials and Methods

HPAI occurrence data acquisition

The dataset containing highly pathogenic avian influenza (HPAI) occurrences, was constructed by extracting relevant information from the EMPRES-i database spanning the period from January 5, 2015, to March 7, 2023. This dataset encompasses comprehensive details pertaining to H5N1 and H5Nx subtypes, including the family of birds involved in transmission, observation and reporting dates, precise location coordinates, and the respective countries where the reports originated. The dataset comprises a total of 17,410 H5Nx and 8,383 H5N1 reported cases.

Environmental data acquisition

Previous literature reviews (Gilbert and Pfeiffer, 2012; Dhingra et al., 2016) have summarised the predictor variables commonly linked to occurrences of HPAI. Here, we explored four sets of environmental variables considered by Dhingra and colleagues (Dhingra et al. 2016), and which are detailed in the table in Supplementary Information Resources S1: host variables (set 1), land cover variables (set 2), eco-climatic variables (set 3), and a risk-based selection of variables performed by Dhingra and colleagues (set 4).

Set 1 includes log10-transformed extensive and intensive chicken density (Gilbert et al., 2015), duck density (Robinson et al., 2014) and human population density from the Gridded Population of the World (Center For International Earth Science Information Network-CIESIN-Columbia University, 2017) database. The chicken density data includes GIS layers representing the global distribution of chickens, categorised into extensive and intensive systems. It uses the model published by Gilbert et al. (2018) applied to the Gridded Livestock of the World (GLW) version 3 chicken layers (Gilbert et al., 2022a). The global duck distribution data were computed using the GLW version 4 (Gilbert et al., 2022b). Set 2 includes land cover information, which was obtained from the global 1-km Consensus Land Cover database by (Tuanmu and Jetz, 2014) Tuanmu and Jetz (2014), and identifies different land cover categories by providing a percentage of representation of each land cover category within raster cells of approximately one square kilometre (http://www.earthenv.org/landcover.html). Additionally, this dataset was supplemented by a layer about the distance of each spatial point to the open water. Set 3 encompasses seasonal and large-scale patterns of eco-climatic indices like daytime land surface temperature (day LST) and the normalised difference vegetation index (NDVI), as documented by Scharlemann and colleagues (Scharlemann et al., 2008). Finally, set 4 is a selection from the previous sets based on prior epidemiological knowledge (Gilbert and Pfeiffer, 2012; Dhingra et al. 2016).

Ecological niche modelling

The different ecological niche modelling analyses were performed with the boosted regression trees (BRT), available in the R package “dismo” (Hijmans et al., 2017). BRT is a machine learning approach that can be used to generate a collection of sequentially fitted regression trees optimising the predictive probability of occurrence given local environmental conditions (Elith et al. 2006) The predictive probability, interpreted as ecological suitability for local virus circulation, ranges from “0” (unsuitable environmental conditions) and “1” (highly suitable environmental conditions). BRT allows us to model complex non-linear relationships between the response and various predictor variables (Elith et al. 2008). It has also been demonstrated that the BRT approach had a superior predictive performance compared to alternative modelling approaches (Elith et al., 2006).

Overall, we performed ten replicated BRT analyses based on each of the seven datasets of occurrence records considered in our study: (i) H5Nx cases in wild birds < 2020 (01/01/2015 - 31/12/2019), (ii) H5N1 cases in wild birds > 2020 (01/01/2020 - 06/03/2023), (iii) H5Nx cases in wild birds > 2020 (01/01/2020 - 06/03/2023), (iv) H5Nx cases in domestic birds < 2020 (01/01/2015 - 31/12/2019), (v) H5N1 cases in domestic birds < 2020 (01/01/2015 - 31/12/2019), (vi) H5N1 cases in domestic birds > 2020 (01/01/2020 - 06/03/2023), and (vii) H5Nx cases in domestic birds > 2020 (01/01/2020 - 06/03/2023). Of note, there was not enough occurrence data to train ecological niche models on H5N1 cases in wild birds < 2020 (only 33 occurrence data from eight different countries — Bangladesh, Bulgaria, China, India, Kazakhstan, Nepal, Romania, and Russia — between the 01/01/2015 and 31/12/2019).

We ran BRT analyses using both occurrence and “pseudo-absence” data, the latter corresponding to randomly sampled absence points in the study area (Elith et al., 2006). To address the heterogeneity of AIV surveillance efforts and to avoid misclassifying low-surveillance areas as unsuitable for virus circulation, we trained the ecological niche models only considering countries in which five or more cases have been confirmed. For China and Russia, given the size of these countries, we however considered the administrative polygons of level 1 (admins-1; provinces and federal districts, respectively) instead of the country polygons for the selection of administrative entities with at least five confirmed cases to be included in the study area. In practice, this implied that the study area used to train our models were restricted to those countries/admins-1 and that we only sampled pseudo-absence points in such restricted study areas. In addition, the sampling pseudo-absence points was based on the human population distribution to further account for potential within country differences in surveillance intensity (Dhingra et al. 2016): the probability to sample a pseudo-absence from a given grid cell was proportional to the log-transformed human population density in that cell. The human population density raster was retrieved from the Gridded Population of the World (GPW) project, which was downscaled to a resolution of 5 arcminutes. As pointed out by Dhingra and colleagues, the absence of the target pathogen in a region may simply be due to the fact that the pathogen has not (yet) been introduced in that region (Dhingra et al. 2016). To limit the impact of pseudo-absences that would have been sampled in those regions, a minimum distance of 10 km and a maximum distance of 1,000 km to the nearest presence points were therefore introduced to discard all simulated pseudo-absence points falling outside that range (Phillips et al. 2009, Dhingra et al. 2016). While we initially simulated ten times more pseudo-absence than presence points, we eventually randomly retained three times more pseudo-absence points than presence points in each country/admin-1 included in the study area.

Given that spatial data typically exhibits autocorrelation, employing a standard cross-validation procedure may result in an overestimation of the BRT model’s accuracy (Randin et al. 2006). To circumvent this issue, we employed a spatial cross-validation procedure to select the optimal number of trees in each BRT model. Specifically, we tested and compared two distinct spatial cross-validation procedures: (i) the procedure introduced by Dhingra and colleagues and based on (five) reference presence points used to cluster the presence and pseudo-absence points in (five) folds according to their nearest distance to a reference point, and (ii) the procedure introduced by Valavi and colleagues and based on a block generation method (Valavi et al. 2019). In the first spatial cross-validation procedure, the five reference presence points were randomly selected among all presence points. If each of the five generated spatial folds did not gather at least 10% of the presence points, the random selection of reference points were re-initiated. To compare the performance of the three cross-validation procedures applied in this study (the standard procedure and the spatial procedures based on the reference points or blocks generation), we estimated the spatial sorting bias (SSB) metric (Hijmans, 2012) ranging from 0 and 1, an SSB near 0 and 1 indicating an extreme impact and the absence of an impact of spatial autocorrelation on the model training, respectively.

All BRT analyses were run and averaged over 10 cross-validated replicates, with a tree complexity set at 4, an initial number of trees set at 10, a learning rate of 0.01, a tolerance parameter set to 0.001, a step size of 5, and considering 5 spatial folds. We evaluated the inferences using the area under the receiver operating characteristic (ROC) curve, also simply referred to as “area under the curve” (AUC). We further used the AUC metric to assess the capacity of our models trained on occurrence data <2020 to predict >2020 distribution of occurrence data. To assess how each environmental factor contributed to the different BRT models, we computed their relative influence (RI) in those models and their response curves showing how ecological suitability varies with one specific factor. RI values were obtained by evaluating the number of times a specific environmental factor was selected for splitting a tree, weighted by the squared improvement to the model as a result of each split, averaged over all trees (Elith et al. 2008).

We explored the four sets of environmental variables assembled by Dhingra and colleagues (Dhingra et al. 2016) and described above in the “Environmental data acquisition” subsection. The BRT analyses presented here were eventually performed based on set 2 for wild bird cases and set 4 for domestic bird cases. This selection is aligned with the methodology of Dhingra and colleagues (2016), prioritising variables which demonstrated the best predictive accuracy and clear epidemiological significance to minimise the potential for coincidental correlations. The RI values obtained when analysing each set of environmental variables separately are all reported in a table in Supplementary Information Resources S2.

Bird diversity indices calculation

Diversity index calculations were performed using the HPAI dataset. Cases lacking information on the date of collection, bird family, and subtype were excluded from the analysis. The diversity of HPAI cases across bird families was quantified using two indices: the Shannon diversity index and the Simpson diversity index. The Shannon index accounts for evenness of bird families (Shannon, 1948), whereas the Simpson index provides a measure of dominance (Simpson, 1949), emphasising the most abundant families in the dataset. Simpson index goes between 0 and 1, where 1 represents infinite diversity and 0, no diversity. Calculations were done using two time points: before and after 2020, and for various data classifications: overall and annually, across all bird families and specifically marine families, as well as for all subtypes (HxNx) and those identified as H5N1 subtype, wild birds, and poultry. These indices were calculated using the R package “vegan” (Oksanen et al., 2023): the dataset was first aggregated by year, HPAI subtype and bird family, and we then used the “diversity” function to compute the Shannon and Simpson indices for each year.

Supplementary figures

Figure S1 (part 1):

Figure S1 (part 2):

Figure S2 (part 1):

Figure S2 (part 2):

Figure S3:

Figure S4:

Figure S5:

Acknowledgements

MFVG and SD acknowledge funding from the European Union Horizon 2020 project MOOD (grant agreement n°874850). MCD is supported by the UKRI GCRF One Health Poultry Hub (Grant No. B/S011269/1), one of twelve interdisciplinary research hubs funded under the UK government’s Global Challenge Research Fund Interdisciplinary Research Hub initiative. MCD and SD are supported by the BELSPO project BE-PIN. SD also acknowledges support from the Fonds National de la Recherche Scientifique (F.R.S.-FNRS, Belgium; grant n°F.4515.22), from the Research Foundation - Flanders (Fonds voor Wetenschappelijk Onderzoek - Vlaanderen, FWO, Belgium; grant n°G098321N) and from the European Union Horizon 2020 project LEAPS (grant agreement n°101094685).

Addtional information

Data and code accessibility

R scripts and related files needed to run the analyses, as well as Supplementary Information Resources S1 (description and source of each environmental variable included in the original sets of variables) and S2 (mean relative influence computed for the environmental variables within each considered set of variables), are all available at https://github.com/sdellicour/h5nx_risk_mapping.