Risk factors and characteristics of H5N1 Highly Pathogenic Avian Influenza (HPAI) post-vaccination outbreaks.

Highly pathogenic avian influenza (HPAI) virus H5N1 is now endemic in South-East Asia but HPAI control methods differ between countries. A widespread HPAI vaccination campaign that started at the end of 2005 in Viet Nam resulted in the cessation of poultry and human cases, but in 2006/2007 severe HPAI outbreaks re-emerged. In this study we investigated the pattern of this first post-vaccination epidemic in southern Viet Nam identifying a spatio-temporal cluster of outbreak occurrence and estimating spatially smoothed incidence rates of HPAI. Spatial risk factors associated with HPAI occurrence were identified. Medium-level poultry density resulted in an increased outbreak risk (Odds ratio (OR) = 5.4, 95% confidence interval (CI): 1.6-18.9) but also climate-vegetation factors played an important role: medium-level normalised difference vegetation indices during the rainy season from May to October were associated with higher risk of HPAI outbreaks (OR = 3.7, 95% CI: 1.7-8.1), probably because temporal flooding might have provided suitable conditions for the re-emergence of HPAI by expanding the virus distribution in the environment and by enlarging areas of possible contacts between domestic waterfowl and wild birds. On the other hand, several agricultural production factors, such as sweet potatoes yield, increased buffalo density, as well as increased electricity supply were associated with decreased risk of HPAI outbreaks. This illustrates that preventive control measures for HPAI should include a promotion of low-risk agricultural management practices as well as improvement of the infrastructure in village households. Improved HPAI vaccination efforts and coverage should focus on medium poultry density areas and on the pre-monsoon time period.

INTRODUCTION

Highly pathogenic avian influenza (HPAI) H5N1 has affected many millions of poultry producers and their families in South-East Asia for more than a decade. HPAI had been successfully controlled in Hong Kong and Thailand by culling of birds on infected premises and intense surveillance [13, 15]. However, countries with a large concentration of backyard poultry and small commercial poultry enterprises have major difficulties controlling the disease. In Viet Nam HPAI has been established since 2003 [2] and it is believed that the virus spread from China to Viet Nam, with the North and the South having different virus sublineages [8]. An intensive vaccination program initiated in 2005 provided lessons to the international community on how the disease can be successfully controlled. Following vaccination, no outbreaks occurred almost over an entire year. Re-emergence of H5N1 began in December 2006 in the Ca Mau province of Southern Viet Nam spreading quickly to other provinces in the region. No human cases were reported as a result of these outbreaks [8].

Our objectives were to (1) describe the spatio-temporal patterns of this epidemic, (2) calculate the incidence rate of outbreak occurrence and (3) identify risk factors for HPAI outbreak occurrence in this post-vaccination period.

MATERIALS AND METHODS

Study area

The study population consisted of all communes south of the Mekong River, Tien River and the Co Chien River system, representing 869 communes in total. A digital map of the administrative boundaries at province, district and commune level of Viet Nam for the study area was obtained from the Viet Nam government’s Department of Animal Health (DAH). The point locations of all commune centroids in the study area were calculated using the geographical information system software ArcView 9.2 (ESRI Inc., 1999–2006) with the ArcView extension software Spatial Analyst 9.2 (ESRI Inc., 1999–2006).

HPAI outbreak information

We used three data sources to summarize all reported HPAI outbreaks from the end of 2006 to the beginning of 2007: (1) outbreaks officially reported to the World Organisation of Animal Health (OIE)1, (2) data published on the DAH website (in Vietnamese)2 and (3) a small set of data from farm-level outbreak investigations conducted by an Australian Centre of International Agricultural Research (ACIAR) project in the Mekong Delta of Viet Nam3. The two websites (OIE and DAH) contained mainly outbreak information at commune level. Village outbreaks reported in the first two data sources and all data from the third data source were aggregated at the commune level. A total of 152 outbreaks were reported between 6th December 2006 and the 29th January 2007.

Space-time cluster analysis

The retrospective space-time scan statistic [6] for point data (commune centroids) was used to identify space-time clusters of outbreaks. The analysis was conducted using the Bernoulli model in the software SatScan4. The study period was the interval from three weeks prior to the first until three weeks after the last outbreak. This was done in order to take account of the incubation period for HPAI which ranges between 2–14 days with up to 21 days in rare cases and to account for the excretion period of HPAI virus following infection [5, 12]5. All communes with at least one outbreak were specified as case communes with all outbreak dates used in the analysis, in some cases accounting for multiple outbreaks within a commune occurring at the same time. All communes without HPAI outbreaks were defined as control communes and were allocated a randomly selected date within the time span from the first to the last outbreak date, since the software does not calculate this automatically [9]. The temporal dimension investigated in this analysis was aggregated into analysis units each representing 7 days. The maximum spatial cluster size used was 20% of the population at risk, while 50% of the total study period represented the maximum temporal cluster size. Only the most-likely cluster identified by the analysis was considered. The analysis was run for 999 iterations, meaning that the smallest possible P-value was 0.001.

Spatially-adjusted incidence rate

Incidence rate calculations were conducted at district level (77 districts) using data for all 869 communes. The risk period started three weeks before the first case occurred (15th November 2006) and ended on 19th February 2007. If no outbreak occurred in a particular commune, the commune was at risk for the whole period, i.e. the maximum risk period was therefore 96 days. If several outbreaks occurred in a commune, only the first outbreak was used in the calculation. The sum of days at risk and the count of occurrences per district were calculated and used for the incidence rate calculation [3].

The spatial analysis software GeoDa 0.9.5-i.5 (Luc Anselin and The Regents of the University of Illinois, 1998–2004) was used to calculate the incidence rate per district. The incidence rate per district was spatially weighted using the spatial empirical Bayes method depending on local sample sizes. A Rook contiguity weight matrix with one order of contiguity was used to represent the neighbouring communes in the analysis.

Risk factor analysis

Spatial random effect logistic regression was used to investigate the association between HPAI outbreak occurrence at commune level and several risk factors. A commune was considered to be an outbreak commune if HPAI occurred at least once during the study period. All communes without HPAI outbreaks were used as controls.

Risk factors investigated included socio-economic and agricultural variables that were considered to have some importance for poultry management. In addition, population and climate-vegetation data were used. These risk factor data were obtained from a variety of sources: the Food and Agriculture Organisation (FAO) Artemis and GeoNetwork, the LandScan Global Population 1998 database and the 2001 Agricultural Census database of the General Statistics Office of Viet Nam. Spatial information was also considered, including the road length and river length within a commune as indicators of density of poultry transport networks. These spatial risk factors were calculated using the ArcView software extension Spatial Analyst 9.2 (ESRI Inc., 1999–2006).

Of the 869 communes in the study region, 259 communes (including 250 controls and 9 cases) had some missing data for at least one risk factor. A missing data imputation using linear regression was conducted to replace these missing values. Non-missing data were used to predict the values of missing data and these were then varied by adding a random value, which in turn was the regression residual from a randomly selected case without missing values [1]. The data file with the imputed values was then used for further analysis.

For ease of interpretation and assessment of non-linearity of effects, all 36 continuous-scale risk factor variables were converted into three categories representing terciles of their distribution. Odds ratios and 95% confidence intervals were calculated for each risk factor using logistic regression in the univariable analysis. In the multivariable analysis, a logistic regression analysis with forward stepwise variable selection was conducted to identify risk factors statistically significant at P < 0.05 using the likelihood ratio test. First-order interaction terms were tested for statistical significance. The addition of a random effect with a spherical spatial covariance structure provided the best model fit. The effect estimates for the final model were converted into odds ratios and their 95% confidence intervals.

The missing data imputation was conducted in SPSS 15.0 (SPSS Inc., 1989–2006), the logistic regression models in STATA 10 (StataCorpLP, 2007) and the final logistic regression analysis with a spatial random effect in SAS for Windows Version 9.1 using the PROC GLIMMIX procedure (SAS Institute, Cary, NC, USA).

RESULTS

A total of 152 outbreaks were recorded over the risk period of 96 days, amongst a total of 66 affected communes. The commune that reported the first occurrence on the 6th December 2006 was also the commune that reported the highest number of outbreaks (13 outbreaks). The location of the significant (P ≤ 0.001) most likely spatio-temporal outbreak cluster comprising of 102 out of 869 communes in the study area is displayed in Figure 1. A total of 90 HPAI outbreaks occurred within the cluster area (radius = 42.5 km). The risk of an HPAI outbreak was 10.2 times higher within the cluster compared to the outside area.

Figure 1.

Communes in the Mekong Delta of Viet Nam with reported HPAI outbreaks between December 2006 and February 2007. Outbreak locations are displayed as black dots (some communes had several outbreaks). The most likely spatio-temporal cluster is encircled with its cluster centre is shown by the star. (A color version of this figure is available at www.vetres.org.)

The neighbourhood-adjusted district level Bayesian incidence rates for HPAI outbreaks are shown in Figure 2. The highest incidence rate for HPAI occurrence was estimated for the district Tran Van Thoi in the province Ca Mau.

Figure 2.

Incidence rate of HPAI outbreak occurrence at district level for the 2006/2007 HPAI outbreak period in the Mekong Delta of Viet Nam. Incidence rates were calculated using spatial empirical Bayesian neighbourhood adjustment. The legend shows categorized percentiles of the incidence rates with the number of districts in each percentile category listed in brackets. (A color version of this figure is available at www.vetres.org.)

The effect estimates for all 35 variables investigated in the univariable analysis are displayed in Table I. The variables included in the final multivariable model are shown in Table II. HPAI outbreak occurrence was associated with sweet potato yield, buffalo density, the rescaled normalised difference vegetation index (NDVI) median from May–October, the proportion of households with electricity and poultry flock density.

Table I.

Results of the univariate analysis of potential risk factors for the 2006/2007 HPAI outbreak period in the Mekong Delta of Viet Nam.

Variable	Categories	Odds Ratio	95% CI	P
River length within a commune (in km)				< 0.001
	≤ 8.18	1
	8.18–20.6	2.19	0.97–4.92
	≥ 20.7	4.73	2.24–9.97
Major road length within a commune (in km)				0.019
	≤ 0.00	1
	0.01–3.51	0.49	0.22–1.10
	≥ 3.52	1.42	0.83–2.44
Re-scaled NDVI annual mean				0.027
	≤ 145.40	1
	145.41–163.14	2.36	1.23–4.54
	≥ 163.15	1.55	0.77–3.10
Re-scaled NDVI annual median				0.016
	≤ 147.15	1
	147.16–164.00	2.37	1.24–4.54
	≥ 164.01	1.48	0.73–3.02
Re-scaled NDVI mean May–October				0.002
	≤ 141.14	1
	141.15–166.15	1.8	1.01–3.21
	≥ 166.16	0.58	0.28–1.22
Re-scaled NDVI median May–October				0.004
	≤ 144.00	1
	144.01–168.00	1.79	1.00–3.18
	≥ 168.01	0.61	0.29–1.27
Re-scaled NDVI mean November–April				0.005
	≤ 148.22	1
	148.23–170.79	0.5	0.27–0.91
	≥ 170.80	0.38	0.20–0.73
Re-scaled NDVI median November–April				0.007
	≤ 148.00	1
	148.01–171.03	0.42	0.22–0.78
	≥ 171.04	0.47	0.26–0.86
Number of households				0.002
	≤ 1995.00	1
	1995.01–3009.52	2.76	1.31–5.83
	≥ 3009.53	3.24	1.55–6.77
Population				< 0.001
	≤ 9768.00	1
	9768.01–14683.00	2.07	0.95–4.51
	≥ 14683.01	3.98	1.94–8.19
Population density (per km2)				< 0.001
	≤ 308.00	1
	308.01–652.35	0.43	0.24–0.75
	≥ 652.36	0.13	0.05–0.31
Proportion of females				0.037
	≤ 0.51	1
	0.51–0.52	0.88	0.50–1.55
	≥ 0.52	0.44	0.22–0.87
Proportion of young (≤ 14 years)				< 0.001
	≤ 0.32	1
	0.33–0.35	2.88	1.26–6.59
	≥ 0.36	4.98	2.27–10.91
Proportion of working age (15–65 years)				0.001
	≤ 0.61	1
	0.62–0.63	0.73	0.42–1.27
	≥ 0.64	0.29	0.14–0.60
Proportion of old (> 65 years)				< 0.001
	≤ 0.04	1
	0.05–0.059	0.45	0.26–0.79
	≥ 0.06	0.11	0.04–0.27
Average No. of persons per household				< 0.001
	≤ 4.77	1
	4.78–5.02	2.57	0.98–6.72
	≥ 5.03	8.66	3.63–20.65
Proportion households with electricity				< 0.001
	≤ 0.26	1
	0.27–0.53	0.22	0.11–0.43
	≥ 0.54	0.22	0.11–0.44
Proportion households with TV				0.004
	≤ 0.43	1
	0.44–0.56	3.05	1.50–6.18
	≥ 0.57	2.29	1.10–4.76
Proportion households with a radio				0.007
	≤ 0.45	1
	0.46–0.53	2.69	1.31–5.51
	≥ 0.54	2.6	1.27–5.36
Proportion of people below poverty line (as defined by Minot and Goletti 2000, [7])				< 0.001
	≤ 0.37	1
	0.38–0.46	1.6	0.92–2.78
	≥ 0.47	0.33	0.15–0.75
Percentage of agricultural land				< 0.001
	≤ 76.78	1
	76.79–80.89	0.11	0.33–0.36
	≥ 80.90	1.31	0.78–2.21
Area used for aquaculture (per km2 of total land)				< 0.001
	≤ 2.42	1
	2.43–56.90	3	0.62–14.55
	≥ 59.91	45.85	11.06–190.14
Paddy yield (Quintal per ha)				< 0.001
	≤ 44.30	1
	44.31–46.95	0.04	0.02–0.11
	≥ 46.96	0.06	0.02–0.15
Sweet potato yield (Quintal per ha)				< 0.001
	≤ 124	1
	125–202	0.08	0.03–0.18
	≥ 203	0.06	0.2–0.26
Maize yield (Quintal per ha)				< 0.001
	≤ 36	1
	37–43	15.58	5.49–44.20
	≥ 44	8.6	2.94–25.19
Income (average per household) (in Vietnamese Dong)				< 0.001
	≤ 737 000	1
	738 000–802 000	1.5	0.88–2.57
	≥ 803 000	0.19	0.07–0.50
Pig density (No. per km2)				< 0.001
	≤ 0.53	1
	0.54–0.96	0.97	0.57–1.63
	≥ 0.97	0.1	0.03–0.34
Cattle density (No. per km2)				< 0.001
	≤ 0.02	1
	0.03–0.13	0.95	0.56–1.61
	≥ 0.14	0.55	0.13–0.23
Buffalo density (No. per km2)				< 0.001
	≤ 0.00	1
	0.01–0.02	0.1	0.04–0.25
	> 0.02	1.16	0.66–2.01
Percentage of urban land use per total area of commune				0.009
	≤ 0.00	1
	≥ 0.01	0.38	0.02–0.85
Proportion of chicken farms				0.001
	≤ 0.47	1
	0.48–0.66	0.67	0.39–1.14
	≥ 0.67	0.23	0.09–0.56
Poultry farm all density (No. farms per km2)				< 0.001
	≤ 0.93	1
	0.94–1.50	0.13	0.02–0.09
	≥ 1.51	0.89	0.04–0.19
Poultry flock density (No. of flocks > 50 heads per km2)				< 0.001
	≤ 0.57	1
	0.58–1.14	1.64	0.97–2.76
	≥ 1.15	0.13	0.04–0.43
Slaughter liveweight density (t per km2)				< 0.001
	≤ 0.92	1
	0.93–3.45	0.5	0.30–0.85
	≥ 3.46	0.07	0.02–0.15
Egg production density (No. per km2)				< 0.001
	≤ 36 956.29	1
	36 956.30–75 451.56	0.12	0.06–0.23
	≥ 75 451.57	0.04	0.01–0.15

Table II.

Results of the spatial random effect logistic regression analysis for risk factors associated with HPAI outbreak occurrence in 2006/2007 in the Mekong Delta of Viet Nam.

Variable	Categories	Odds Ratio	95% CI	P
Sweet potato yield (Quintal per ha)				< 0.001
	≤ 124	1
	125–202	0.03	0.01–0.11
	≥ 203	0.04	0.01–0.36
Buffalo density (No. per km2)				< 0.001
	≤ 0.00	1
	0.01	0.01	0.00–0.05
	≥ 0.02	0.26	0.07–0.94
Re-scaled NDVI median May–October				0.005
	≤ 144.00	1
	144.01–168.00	3.69	1.68–8.07
	≥ 168.01	2.37	0.84–6.67
Proportion households with electricity				0.004
	≤ 0.26	1
	0.27–0.53	0.27	0.12–0.64
	≥ 0.54	0.38	0.16–0.89
Poultry flock density (No. of flocks > 50 heads per km2)				0.007
	≤ 0.57	1
	0.58–1.14	5.44	1.56–18.98
	≥ 1.15	0.61	0.09–4.21

DISCUSSION

Viet Nam was one of the first countries to introduce large-scale poultry vaccination against HPAI and it is very likely that this strategy contributed to the apparent cessation of cases in poultry and humans for a period of about 12 months starting from November 2005. However, in December 2006, HPAI H5N1 outbreaks were reported on poultry farms in three southern Vietnamese provinces followed by four more provinces in January 2007. The outbreaks were believed to be triggered by illegal hatching and raising of unvaccinated ducklings6 ,7. Uncontrolled hatchings were banned by the veterinary authorities in Viet Nam since the 19th February and only registered breeders were allowed to produce ducklings (Directive 3065/2005 by the Ministry of Agriculture of Viet Nam on the 4th February 2005). Also poor vaccination techniques and vaccine quality as well as low vaccination coverage were considered as possible reasons for the re-emergence of HPAI in Viet Nam7.

All HPAI outbreaks were clustered in the southern part of the Mekong region, in an area with a dense network of channels and islands, where duck production represents a significant farming income. The highest incidence rate of HPAI was found in the district where the first outbreak and the highest number of outbreaks during this epidemic occurred, with all of these outbreaks occurring in the first month of the epidemic. This indicates that HPAI spread from the initially infected areas, where the virus presence accumulated in the face of unsuccessful control, to the neighbouring areas in a somewhat radial pattern. Spread of HPAI is believed to occur through mechanisms such as scavenging of domestic poultry, via transport of birds and movement by people [14]. The lower Mekong Delta is also a bread basket for many agricultural products with rice being the most important crop in Viet Nam and in the whole region8. Rice cropping intensity has previously been associated with increased risk of avian influenza outbreaks [4]. We were able to show that the risk of HPAI outbreaks is actually reduced in areas where agriculture activities other than rice production such as production of sweet potatoes are widespread. Sweet potato is an important food for people and is also commonly used as pig feed [10]. Although the Mekong Delta accounts for a relatively small area of sweet potato cultivation compared to the Red River Delta, it produces by far the highest average yields, due in part to hotter weather and a longer growing period9. Field experiments to improve low-yield summer-autumn (July–October) rice production in the Mekong Delta indicated that rotation of rice with sweet potato production significantly improved the rice yield. Although speculative at this point, diversifying rice production to multiple cropping might be beneficial in two ways: by increasing rice yield and perhaps by indirectly reducing the risk of HPAI outbreaks during high risk periods when sweet potatoes are planted instead of rice. The latter may be due to sweet potato leaf from the harvest being unsuitable as sole source of duck feed and primarily used for feeding of pigs, buffaloes and cattle. Therefore harvested sweet potato fields will not satisfy the dietary requirements of ducks. This is in contrast to rice paddies, where relatively significant amounts of grain remain in the paddy field following the rice harvest. As a consequence, grazing of ducks may be less common on sweet potato fields, thereby resulting in a reduced HPAI risk. However, as duck grazing and crop cultivation is often conducted by different producers with different resources targeting a range of markets for each product, a potential change in agricultural production methods may not be economically acceptable for some of them.

Buffalo density was also associated with a reduced risk of HPAI outbreaks. Buffaloes are used for ploughing, for milk and meat production and for the production of manure as fertilizer. The number of buffaloes has declined dramatically over the past years in the Mekong Delta due to increasing use of tractors and reduction of grazing land. In 1995, the Mekong Delta’s buffalo population represented only 4.95% of Viet Nam’s total buffalo population10. In our results the effect of higher buffalo density may represent a higher density of non-cultivating farm management systems where buffaloes are raised for milk and meat production rather than ploughing, and thereby reflect from an HPAI perspective less risky agricultural activities compared with rice production.

Supply of electricity in the household might be a surrogate for household individual wealth or of general infrastructure in the region and the isolation of particular households. Village-based households will have a higher chance of being supplied with electricity, while isolated households in the channel system of Mekong might be more difficult to provide with a permanent power supply. These isolated households are likely to be less easily accessible for veterinary authorities, and therefore may be more likely to be missed during a vaccination campaign.

Most interesting was that medium poultry density, which probably represents density of small commercial poultry enterprises, was associated with an increased outbreak risk. Areas of higher density of poultry had a lower risk (although not significant) – suggesting that larger poultry producers have higher levels of biosecurity or perhaps a better HPAI vaccination coverage.

The re-scaled medium NDVI index for the wet season was associated with an increased risk of HPAI outbreaks. Climatic vegetation related factors in the wet season influencing HPAI occurrences have been mentioned by others [11]. Re-scaled NDVI values range from 0 to 200, with values less than 100 representing water and other non-vegetation surfaces and values greater than 100 representing vegetation surfaces11. Therefore, the re-scaled NDVI categories identified in our study represent waterways with widespread vegetation – which probably refers to areas where duck may scavenge. These areas might also be potential virus reservoirs. Towards the end of the rainy season in September and October, the combination of floodwaters from the rivers, local rainfall and tidal inundation can result in the flooding of large areas of the delta. This might also influence the abundance of a variety of wild birds associated with these wetlands and might pose a further risk for HPAI occurrence.

The epidemic wave investigated here could perhaps be also related to the increased transport frequency and production intensity in preparation for the Tet festival which was celebrated on the 18th February 2007, with the holiday period in Viet Nam beginning one to two weeks before this date. Previous work also suggested that the HPAI epidemics in 2004 and 2005 in Viet Nam were associated with higher production and widespread poultry trade prior to the Tet festival [11]. This repeated HPAI temporal outbreak pattern signifies that intensified HPAI awareness, alertness, surveillance combined with good vaccination coverage during the pre-Tet months should be a focus of HPAI control strategies in Viet Nam.

In our dataset some of the risk factor variables had missing values. In such a situation, statistical analysis procedures will usually exclude the observations which have missing values in any of the risk factor variables, which results in a loss of statistical power and may introduce bias (if the values are not missing at random). An alternative approach would be to replace the missing values with imputed values, so that statistical power levels are not compromised. However, missing data imputation using regression techniques may itself introduce bias, because all observations with the same values for the risk factor variables will be imputed using the same value for the missing variable. To address this problem we applied a stochastic substitution, which uses the regression technique but adds the regression residual from a randomly selected case from the set of cases without missing values to the predicted value [1].

Overall our study has shown that in the post-vaccination period HPAI outbreaks occurred as a tight cluster. Previous work described mainly pre-vaccination risk factors for HPAI. Pfeiffer et al. [11] indicated that bird density and vegetation indices were associated with increased risk of HPAI outbreaks – these were also the most important risk factors in our post-vaccination analysis. However, we were able to show that in particular the rainy season with its temporal flooding might have provided suitable conditions for the re-emergence of HPAI by influencing the abundance of domestic ducks and wild birds, by expanding the areas of possible contacts between waterbirds and poultry and perhaps by spatially expanding the virus distribution in the environment. Two HPAI vaccination rounds are usually conducted by veterinary authorities in Viet Nam (March–April and August–September). The temporal flooding might have also hindered access to some poultry flocks thereby during the second, rainy-season vaccination campaign leaving some flocks unvaccinated.

HPAI is a global animal health problem, a zoonosis and a threat to local economies. Our study provided insights with respect to how the disease re-emerged despite of a range of control measures including vaccination being in place and it highlighted important HPAI risk factors for the South and South-East Asian subcontinent. Furthermore we illustrated that certain agricultural management practices together with an improved infrastructural and socio-economic status of households have the potential to decrease the risk of HPAI outbreaks – information that will be essential for policy makers for defining improved HPAI control strategies in the region and globally.