Clade-level Spatial Modelling of HPAI H5N1 Dynamics in the Mekong Region Reveals New Patterns and Associations with Agro-Ecological Factors

Abstract

The highly pathogenic avian influenza (HPAI) H5N1 virus has been circulating in Asia since 2003 and diversified into several genetic lineages, or clades. Although the spatial distribution of its outbreaks was extensively studied, differences in clades were never previously taken into account. We developed models to quantify associations over time and space between different HPAI H5N1 viruses from clade 1, 2.3.4 and 2.3.2 and agro-ecological factors. We found that the distribution of clades in the Mekong region from 2004 to 2013 was strongly regionalised, defining specific epidemiological zones, or epizones. Clade 1 became entrenched in the Mekong Delta and was not supplanted by newer clades, in association with a relatively higher presence of domestic ducks. In contrast, two new clades were introduced (2.3.4 and 2.3.2) in northern Viet Nam and were associated with higher chicken density and more intensive chicken production systems. We suggest that differences in poultry production systems in these different epizones may explain these associations, along with differences in introduction pressure from neighbouring countries. The different distribution patterns found at the clade level would not be otherwise apparent through analysis treating all outbreaks equally, which requires improved linking of disease outbreak records and genetic sequence data.

Figure 1. Descriptive information.

(a) Study area. Three localities were marked on this map: A: The central plain of Thailand (the green hatched area); B: The Mekong Delta; C: The Red river Delta. The dark grey areas show the mask used in this study and corresponds to International Union for Conservation of Nature (IUCN) conservation areas, permanent water bodies and city centres where poultry is presumable absent. The graticule is composed of a 5-degree increments (http://www.naturalearthdata.com/) and the coordinate system is ‘EPSG:3148’. The data used to produce these maps were all from public sources, and the country limit data are from the FAO Global Administrative Unit Layers (GAUL) database. This figure was built with the R-3.2.3 software (https://cran.r-project.org/). (b) Summary diagram of the models and outbreaks clade allocation steps. b. 1. General model (red box): Step 1: Extraction of oubreak dataset (EMPRES-i (http://empres-i.fao.org/); Department of Animal Health (Hanoi, Viet Nam)); Step 2: Random sampling of pseudo-absences from the population at risk and computation of ‘P(H5N1)’, the HPAI H5N1 presence probability. The random sampling procedure of pseudo-absences is repeated for each bootstrap of the analysis. b. 2. Clade models (blue box): procedure for linking outbreaks to HPAI H5N1 sequences and the clade models. Step 3: A buffer of a given distance is drawn around each outbreak location (outbreak data) and the sequence data were extracted (Bioinformatics OpenFlu database (http://openflu.vital-it.ch/); Influenza Research Database (IRD, http://www.fludb.org/)); Step 4: Sequences located within that buffer are identified. Step 5: One of these sequences is randomly selected and used for analysis. The entire allocation procedure is repeated for all outbreaks and bootstraps of the analysis.

Figure 2. Distribution maps of outbreaks.

Distribution of outbreaks (yellow) used for the HPAI H5N1 general model (top left), and maps of one realization of the outbreaks with clade assignation for the three time periods. Note that the colours of points individual points is transparent such as to show areas with overlapping points with more saturated colours. The mask corresponds to International Union for Conservation of Nature (IUCN) conservation areas, permanent water bodies and city centres where poultry is presumable absent. The graticule is composed of a 5-degree increments (http://www.naturalearthdata.com/) and the coordinate system is ‘EPSG:3148’. The data used to produce these maps were all from public sources, and the country limit data are from the FAO Global Administrative Unit Layers (GAUL) database. This figure was built with the R-3.2.3 software (https://cran.r-project.org/).

Figure 3. Suitability model predictions.

General HPAI H5N1 suitability model (P(HPAI)) for the entire time period (a); and for each clade P(Clade) in their respective time periods (b), with colour intensity proportional to the predicted suitability for a particular clade. Note that in Red, Green, Blue (RGB) visualisation, pixels in black correspond to pixels where the predicted suitability for all three clades was close to zero, as the intensity of the colour varies with suitability. The mask corresponds to conservation areas of International Union for Conservation of Nature (IUCN), permanent water bodies and city centres where poultry is presumable absent. The graticule is composed of a 5-degree increments (http://www.naturalearthdata.com/) and the coordinate system is ‘EPSG:3148’. The data used to produce these maps were all from public sources, and the country limit data are from the FAO Global Administrative Unit Layers (GAUL) database. This figure was built with the R-3.2.3 software (https://cran.r-project.org/).

Figure 4. Mean relative contribution of each predictor variable to the prediction for the general and clade models for individual time periods.

The relative contribution plots show the statistic importance of each predictor variable in the model. It is measured as a percentage on the Y-axis for each predictor variable: human population density (HpopDn, human/km² - log10), chicken density (ChDn, head/km² - log10), duck density (DuDn, head/km² - log10), cropping intensity (CropLd, intensity scale from 0 to 4), accessibility (travel-time) to major cities (Access, min), and the proportion of land occupied by water (Wapc). Log10 notes a logarithm transformation of base 10 of predictor variables.

Figure 5. Partial dependence plots or BRT profiles for the four top predictor variables for the general suitability model and for individual clade models.

The partial dependence plots show the predicted dependence between the dependant variable of BRT models on the X-axis (the probability of virus presence) and each predictor variable on the Y-axis. Four predictor variables were selected for this figure: human population density (HpopDn, human/km² - log10), chicken density (ChDn, head/km² - log10), duck density (DuDn, head/km² - Log10) and accessibility (travel-time in minute) to major cities (Access, min). Log10 notes a logarithm transformation of base 10 of predictor variables. The dashed line represents the mean profile, whilst transparent lines represent each bootstrap. The black transparent ticks at the bottom of each plot represent the observed distribution of predictor variables for one bootstrap and the corresponding dataset.