Surveillance efforts aimed at early detection of avian influenza (AI) in wild birds has markedly increased in the European Union (EU) since the emergence of the Asiatic H5N1 highly pathogenic (HP) AI virus (AIV) strain in 2005. In addition to the severe economic losses that HPAI infection inflicts on commercial and non-commercial poultry industries in affected regions, its zoonotic potential, plus fears that the disease may lead to a pandemic in humans triggered a social alarm that has probably influenced the large number of samples collected throughout the EU in wild birds [Reference Martinez1].

A number of routes of entry have been proposed or described for AIV, including the trade of asymptomatic domestic ducks [Reference Kim2, Reference Quan3], importation of infected meat [Reference Irvine4, Reference Harder5], and illegal introduction of captive birds [Reference Dudley6]. However, during 2005–2008 almost half of the H5N1 HPAI-infected European countries reported initial outbreaks in wildlife with no further outbreaks recorded in poultry. Moreover, the risk of introduction of HPAI into Spain through the legal importation of poultry has been estimated as negligible [Reference Sánchez-Vizcaíno7]. For those reasons and although other routes of entry should be considered when designing a broad surveillance programme for AI, it is likely that the risk of HPAI epidemics in Spain is mostly influenced by the probability of introduction of AI-infected wild birds.

Sampling of wild birds is a labour-intensive, costly, and time-consuming task that has not been exempt from discussion at the decision-making level in the EU and in other regions affected by the disease. The Asiatic H5N1 strain was the first HP AIV isolated from multiple species of wild birds, and epidemiological factors expected to influence the sampling design, such as pathogenesis, incubation period, excretion titres, and immune status of the wild population, were largely unknown at the time the sampling schemes were designed. Historically, birds from the Anseriformes and Charadriiformes orders have been considered reservoirs of AIVs [Reference Brown and Stallknecht8]. Consequently, EU surveillance schemes were designed to target periods of time and geographical areas in which it was most likely to find a large number or density of birds from those orders that could have migrated from infected areas and that were likely to have contact with local poultry. Aspects of animal behaviour, such as gregariousness of bird species, were also used to target the sampling scheme to areas at high risk for HPAIV infection [9]. Moreover, field and experimental studies were conducted in order to elucidate epidemiological aspects of the disease that could help to identify, in particular, which species were at highest risk for the disease. However, those studies suggested that parameters related to the transmission and spread of the disease were highly variable among species [Reference Brown and Stallknecht8, 10]. For those reasons, and because the design of this approach does not consider factors such as risk of infection at the origin of migration or the distance of migration, wild-bird sampling schemes that are based only on the expected distribution of the susceptible wild-bird population may not be the most appropriate approach to maximize the probability of early detection of HPAI viruses in the EU.

An alternative approach to overcome these limitations may be implementation of a risk-based surveillance system that takes into account the historical information on places and times of the year at which HPAI outbreaks were most likely to occur, as suggested by the identification of significant space–time clustering [Reference Iglesias11]. Unfortunately, such an approach cannot be implemented in countries in which HPAIV incursions have been incidental, such as Spain, or in countries where the disease has never been reported. Alternatively, quantification of the risk of H5N1 HPAI incursion may be approximated using information on the probability of wild birds' direct contact with areas in which the infection is likely to occur.

Spanish wetlands are used as wintering sites by about 1·5 million aquatic birds [Reference Marti and Del Moral12]. In line with EU standards and following the general approach recommended by EFSA [9, 13], H5N1 HPAI surveillance efforts in wild birds in Spain have been based both on an active and passive collection of samples. Collection of samples for AI surveillance in Spain is administered by each of the 17 autonomous communities, which are divided into 50 provinces and 497 veterinary units (VUs). For purposes of active surveillance, the minimum number of samples to be collected by each of the 17 autonomous communities of Spain has been specified by the national government [14]. This minimum number of samples was estimated based on the abundance of wild birds in each of the four areas into which Spain has been divided, and that have been generically referred to as North, Mediterranean, South, and Central.

In October 2009, a HPAI outbreak was reported for the first time in a Spanish layer farm. The outbreak was located in Almoguera, province of Guadalajara, autonomous community of Castile-La Mancha, and was caused by a H7 HPAIV strain. The outbreak affected a population of 308 640 birds, with morbidity, mortality, and fatality rates at the time of detection of 9·72%, 9·72%, and 100%, respectively. The entire farm was depopulated [15].

The objective of this report was to quantify the spatial variation in the relative risk for H5N1 HPAI virus introduction into Spain considering the probability of migratory movements from areas of Europe in which outbreaks are most likely to occur. Furthermore, the association between risk of introduction and probability of sampling under the current surveillance scheme as applied in the country was also estimated.

Relative risk for H5N1 HPAI introduction into Spain was quantified at the VU level, which is the smallest administrative level at which actions are taken on the prevention and contingency planning for wild bird species in Spain. The VUs of Ceuta and Melilla, and of the Canary Islands, which are not located in continental Spain, were not assessed in this study. Thus, 492 (99%) of the 497 VUs were included in the analysis.

The probability that an H5N1 HPAI-infected water bird migrates from any location in Europe into Spain during winter (P _e) was estimated as the product of the probability of H5N1 HPAI infection estimated for each of 135 500 20×20 km² cells c in Europe (P _i) and the probability that a wild bird that originated from cell c reaches Spain (P _r). A detailed description of the procedure used to estimate P _i is available elsewhere [Reference Martinez1]. Briefly, a co-kriging model was used to quantify the value of P _i throughout Europe, using data on the incidence of HPAI outbreaks reported from 2005 to 2007 and on the distribution of the susceptible avian population, both wild and domestic. The value of P _r was computed by fitting an inverse beta function to the observed distribution of the distance-dependent probability of bird movements based on ringing recovery information gathered by the Spanish Office of Migratory Species (SPMS, Ministry of Environment, and Rural and Marine Affairs, MARM). In this way, it was assumed that birds departing from cells located closer to Spain had a greater likelihood of reaching the country compared to those from distant locations.

Information in the SPMS database includes location (x, y coordinates) and dates of ringing and recovery, bird species, and ring ID for all birds identified from 1960 (8875 records). A broad spectrum regarding the wintering period in Spain was considered to take into account variation in migration patterns within species. The distribution of expected patterns of autumn migration was fitted using a large period of time (50 years) in order to characterize, with greater confidence, the true variation in the data. Only records that matched the following requirements were used for the analysis:

1. (1) The bird belonged to one of the 25 wintering aquatic species that, as a consequence of their abundance, distribution, and gregarious habits in Spain, are believed to present most of the risk for H5N1 HPAIV introduction into the country [Reference Martinez16].

2. (2) The bird was recovered in Spain but ringed in another European country.

3. (3) The movement, as recorded in the database, was on a southern (i.e. from north to south), southeastern, southwestern, or western direction.

4. (4) The bird was recovered in Spain between July and March (>92% of total data), which corresponds to the time from the early autumn migration until the end of the wintering period.

5. (5) Time between ringing and recovery was <3 weeks, which approximates the duration of virus excretion [10].

After data cleaning was performed, 4331 of 8875 bird movement records met the specified criteria and were available to fit the inverse beta function used to model the value of P _r.The geographical extension of Spain was divided into 43 200 5×5 km² cells in a matrix of 240 columns and 183 rows, and the value of P _e was computed as the product of P _i and P _r for each cell i that received at least one of the 4331 bird movements recorded here.

Inverse distance weighting (IDW) was used to produce an isopleth map of the value of P _e throughout Spain (P _w) (Fig. 1). The IDW technique interpolated the values of P _e estimated for each single cell i that received at least one of the 4331 bird movements (n=4306), to compute the value of P _w for every single cell c in Spain. Values of P _w were computed assuming that the influence of a given cell i on any other cell c decreased with distance d, so that:

where P represents a scale factor. The value of P that resulted in the smallest error of the estimates, which was computed by comparing P _w with P _e in every cell i in which values of P _e were available, was 1. Consequently, a value of P=1, which is equivalent to assume a linear function for the influence of distance, was assumed for the computation of IDW. The distribution of P _e was bimodal, with most (87·2%) of the 4306 cells that received at least one movement having values of P _e<0·0025, some cells (12·8%) showing values of P _e>0·01, and no cell showing values of 0·0025<P _e<0·01. Accuracy of the predictions was evaluated by comparing the values of P _w with P _e for those locations for which values of P _e were available. Using a cross-validation process, it was estimated that values of P _w were significantly higher for cells with P _e>0·01, compared to cells with P _e<0·0025 (Mann–Whitney test, P<0·01) and that the error of the predictions (P _e – P _w) was spatially independent (Moran's I=−0·007). These results suggest that IDW, which is a relatively simple interpolation method, resulted in accurate predictions of P _w.

Fig. 1. Quantile distribution of the inverse distance weighting (IDW) value of the probability that an H5N1 HPAI-infected water bird migrates from any location in Europe into Spain during winter (P _e).

For each VU_i, the relative risk (RR_i) of HPAI introduction was estimated as the values of P _w predicted for cells within the VU (P _wi) divided by the lowest value of P _wi computed for any VU_i in Spain. Therefore, a value of RR=1 was assumed for the VU for which the minimum value of P _wi was estimated and, for any other VU, the value of RR increases proportionally to the increment in the risk.

The value of RR_i was used to compute the number of birds that would have been sampled per VU (m _e) if the number of birds sampled in 2007 throughout Spain (∑m _i) had been distributed proportionally to the risk as

The difference (d _i) between m _ei and the number of birds sampled per VU in 2007 (m _i) was computed for each VU_i. The value of d _i gives an estimate of the deviation between the number of samples collected per VU in 2007 (m _i) and the number of samples that would have been collected (m _ei) if a risk-based sampling design had been applied. The values of m _i were provided by MARM.

Abundance of aquatic wild birds per VU_i (a _i) was computed for the species selected for the analysis using the Spanish Society of Ornithology (SEO/BirdLife 2007) database, in which bird count per wetland and per species is gathered for the month of maximum abundance of wintering birds, i.e. January, 2007.

Correlations between the probability that a bird was sampled from a certain VU_i (P _mi=m _i/∑m _i), the probability that a bird was present (abundance) in a certain VU (P _ai=a _i/∑a _i), and the probability that a sample would have been obtained from the VU if a risk-based sampling scheme had been used (Pm _ei=m _ei/∑m _ei) was assessed by computing Spearman's correlation coefficients. A Spearman correlation coefficient >0·6 (P<0·05) was assumed to indicate a strong association between the parameters.

Geographical areas in which the value of RR_i was significantly higher or lower (P<0·05) than that expected under the null hypothesis of even distribution of those parameters throughout Spain was identified using the two-tailed normal model of the scan statistic [Reference Kulldorff17]. A large number of geographical circles were alternatively placed on the centroid of each VU; circles had a maximum radius that included, at the most, 50% of the population at risk. The mean value of RR_i was computed within each circle and compared to that expected under the null-hypothesis of even distribution of RR_i, which was computed using 999 replications of a Monte Carlo process. Circles in which the mean value of RR_i was significantly (P<0·05) lower or higher than that expected were assumed as clusters of low or high values of RR_i, respectively. The procedure was repeated using d _i, rather than RR_i, as the variable of interest, to identify geographical clusters of VUs in which the value of d _i was significantly (P<0·05) lower or higher than that expected.

The @Risk 4.2 (Palisade Corp., USA) for Excel® software was used to fit the data to theoretical distributions, ArcGIS 9.2 (ESRI™) and its Geostatistical Analyst (ESRI) extension were used to create maps and to compute the IDW, SPSS 15.0 for Windows (SPSS Inc., USA) was used to compute the Spearman correlation tests, and SatScan 8.0 (www.satscan.org) was used to run the spatial scan statistics tests.

The risk of a H5N1 HPAI incursion, arising from wild, migratory birds, was unevenly distributed throughout Spain. Interestingly, areas estimated to be at the highest risk for H5N1 HPAI introduction followed a northeast–southwest direction through continental Spain. Although bird movements within Spain can be in any direction, the northeast–southwest direction is one of the main migration pathways in Spain. Coincidently with previous studies that identified the highest risk for AI introduction in northwestern and eastern regions of Spain [Reference Martinez16], the scan statistic identified two significant clusters (P<0·05) of areas at high RR for H5N1 HPAI introduction into Spain (Fig. 2). The RR in the primary cluster was 2·17 times higher than the mean risk of the country and had a radius of 108 km that included 43 VUs mostly located within the provinces of Zaragoza, Soria, La Rioja, Álava, and Navarra. The secondary cluster had 61-km radius and included eight VUs, most of which were located in the provinces of Sevilla and Huelva, and in which the RR was 3·74 times higher than the mean RR estimated in Spain.

Fig. 2. Quantile distribution by veterinary unit (VU) of the relative risk (RR) of HPAI introduction and the two significant clusters identified by the scan statistic as high RR areas of HPAI H5N1 introduction into Spain.

Data on aquatic birds' abundance data and on AI surveillance sampling were available for 71% and 48% VUs, respectively. Correlations between observed sampling (P _mi) and wild water-bird abundance (P _ai) (n=178, Spearman's ρ=0·36, P<0·05), between observed (P _mi) and estimated (P _me) sampling (n=236, Spearman's ρ=−0·10, P=0·12), and between wild water-bird abundance (P _ai) and estimated sampling (P _me) (n=350, Spearman's ρ=−0·07, P=0·22) were low (Spearman's ρ <0·6), suggesting low association between the number of samples collected per VU, abundance of susceptible species in the VU, and risk for HPAI introduction into the VU.

It is likely that the spatial distribution of sampling for HPAIV practised in Spain in recent years has been influenced by social, political, and economic factors and by the urgency to design and implement a sampling protocol with limited available information. These drawbacks may have resulted in a deviation of the distribution of collected samples compared to those expected if abundance- or risk-based designs were used. Moreover, factors influencing the risk-based sampling scheme include the probability of HPAIV infection at the origin, frequency and distribution of migration, and the probability that an infected bird reaches Spain. None of those factors are considered in an abundance-based sampling scheme, which results in major differences in the spatial distribution of the probability of sampling.

Although most of the VUs (83%) collected fewer samples than the expected number if samples had been collected using a risk-based design, as indicated by values of d _i>0 (Fig. 3), such deviation does not necessarily mean that the efforts made by Spain for early detection of HPAIV infection have been insufficient. Certainly, previous studies have suggested that wild bird surveillance for HPAIV was more intense in Spain, compared to other EU countries and considering the expected risk of infection of the country [Reference Martinez1]. However, the results presented here suggest that such an increase in the number of samples collected would have been most efficient if the samples had been collected from regions at the highest risk for HPAIV infection.

Fig. 3. Quantile distribution by veterinary unit (VU) of the absolute difference (d _i) between expected sampling (m _e) and observed sampling (m _i).

There were no geographical clusters of VUs in which the mean value of d _i was significantly lower or higher than that expected under the null hypothesis of even distribution of d _i, as indicated by the results of the scan statistic.

The results presented here may be biased as a consequence of the non-probabilistic scheme used for the recovery of bird rings, which may have biased the risk towards those regions in which recovery was most likely to occur or towards regions that receive birds that are more easily identified than others. Although many of the rings were recovered during hunting activities, they represent <45% of the total recoveries in Spain and most of the remaining 55% recoveries were conducted by qualified ornithologists and included non-hunting species such as larids or waders. Although it cannot be considered as a validation of the results presented here, it is noteworthy that the HPAIV outbreak reported in October 2009 was located in a VU at high risk for HPAIV infection which, for that reason, would have been intensively sampled if a risk-based surveillance approach had been in place. However, no samples were collected from this VU in 2007, suggesting that although bias in our predictions cannot be ruled out, the scheme proposed in the current study is likely to result in more accurate results than the sampling schemes currently used in Spain and other European countries.

The risk estimates presented here are based on predictions of the probability of infection at origin and will remain valid and accurate as long as the conditions observed at the time they were computed remain stable and constant. For that reason, risk estimates and risk-based sampling schemes such as the ones estimated and proposed here must be revised and updated on a regular basis in order to adjust them to reflect changes in ecological and epidemiological conditions.

Many of the recent HPAI epidemics in European poultry were associated with contact with infected wild birds. Quantification of the province-specific risk imposed by wild-bird migration is important because, along with information on poultry or poultry-farm density and biosecurity, it provides the basic information required to formulate a risk-based surveillance plan for either or both wild- and domestic-bird populations.

In conclusion, the current study quantified the risk for introduction of HPAIV into Spain associated with migration of wild birds. Results were used to propose a risk-based design for HPAIV surveillance that considers probability of migration of wild birds and probability of infection at origin would be more effective in terms of probability of detection of HPAIV in Spain than the sampling schemes recently used or proposed in the country, mitigating the consequences of such incursions. Ultimately, the results presented here could contribute to increase the efficacy and efficiency of surveillance systems for HPAIV and to mitigate the impact of HPAI in Spain, other European countries and throughout the world where a similar approach is designed and implemented.