Insights into the spatiotemporal dynamics of West Nile virus transmission in emerging scenarios

The incidence of West Nile fever (WNF) is highly variable in emerging areas, making it difficult to identify risk periods. Using clinical case records has important biases in understanding the transmission dynamics of West Nile virus (WNV) because asymptomatic infections are frequent. However, estimating virus exposure in sentinel species could help achieve this goal at varying spatiotemporal scales. To identify the determinants of inter-annual variation in WNV transmission rates, we designed a 15-year longitudinal seroepidemiological study (2005–2020) in five environmentally diverse areas of southwestern Spain. We modeled individual annual area-dependent exposure risk based on potential environmental and host predictors using generalized linear mixed models. Further, we analyzed the weight of predictors on exposure probability by variance partitioning of the model components. The analysis of 2418 wild ungulate sera (1168 red deer - Cervus elaphus - and 1250 Eurasian wild boar - Sus scrofa) with a highly sensitive commercial blocking ELISA identified an average seroprevalence of 24.9% (95% confidence interval (CI): 23.2–26.7%). Antibody prevalence was slightly higher in wild boar (27.5%; CI: 25.1–30.1%) than in deer (22.2%; CI: 19.8–24.7%). We observed a spatial trend in exposure, with higher frequency in the southernmost areas and a slight, although area-dependent, increasing temporal trend. Host-related predictors were important drivers of exposure risk. The environmental predictor with the highest weight was annual cumulative precipitation, while temperature variations were also relevant but with less weight. We observed a coincidence of spatiotemporal changes in exposure with the notification of WNF outbreaks in horses and humans. That indicates the usefulness of wild ungulates as sentinels for WNV transmission and as models to understand its spatiotemporal dynamics. These results will allow the development of more accurate predictive models of spatiotemporal variations in transmission risk that can inform health authorities to take appropriate action.


Introduction
The distribution of flaviviruses has increased considerably worldwide [1]. The vast majority of flaviviruses are vector-borne, and >50% of them are zoonotic and cause emerging diseases of variable animal and public health impact. Some flaviviruses of the Japanese encephalitis virus (JEV) serocomplex group, e.g., West Nile virus (WNV) or Usutu virus (USUV), have expanded their distribution range recently. This has favored the appearance of new West Nile fever (WNF) outbreaks and Usutu cases in livestock, humans, and wildlife [2]. Flavivirus infections often manifest clinically with signs of neuroinvasive disease that may progress to lethal disease [3]. In this regard, WNV and USUV are among the most impacting European flaviviruses [4]. However, JEV, Bagaza virus and other endogenous mosquito-borne flaviviruses have also been detected [5][6][7]. The flaviviruses of the JEV serocomplex are maintained by mosquitoes of the genus Culex and birds. However, their ecology is more complex and many other actors are involved in their local maintenance and transmission, e.g., abiotic environmental conditions and mammal presence/abundance [8,9].
In Spain, WNF emerged in horses, wild birds, and even humans at the beginning of the 21st century [2,10]. In 2007, WNV was confirmed as the cause of death of golden eagles (Aquila chrysaetos) in south-central Spain [11]. Neither human nor horse WNF case was reported until 2010 in the country [12]. Since then, the distribution area of WNV has expanded northwards. At present, lineages 1 and 2 are present in the southwest and northeast of Spain, respectively, with current incursions of both lineages in areas of the country distant from those in which they were initially detected. Two thousand and twenty was a remarkable epidemic year with 77 human WNF cases and eight casualties [2]. The short viremia and the high frequency of asymptomatic cases hinder early detection and favor the spread of WNV towards inland Iberia.
Culex spp. mosquitoes are widely distributed in Spain [13], and their remarkable adaptability to different environments, ranging from pond or reservoir waters to abandoned tires [14], makes them almost ubiquitous. In southern continental Spain, the most abundant Culex spp. are C. pipiens and C. theileri [15]. A female mosquito needs to feed on a viremic highly competent bird [4] or co-feed with an infected competent mosquito on a host [16] to become infected with WNV and transmit it while feeding again on a host. Culex pipiens feeds preferentially on birds but it also feeds on mammals, while C. theileri prefers mammals and can, eventually, feed on birds [14]. Likewise, the local abundance of Culex spp. and the transmission probability of WNV are shaped not only by weather conditions but also by host availability [17]. Wild ungulates experienced a demographic explosion in the Northern Hemisphere from mid 20th century [18] whose impact on mosquito dynamics is poorly understood. Wild ungulates are good WNV sentinels [8,9]. They would indeed be a better WNV monitoring group in Iberia when compared to birds because (1) they are widely distributed and abundant, (2) they are hunted in large numbers annually, and (3) their higher body size would potentially attract a higher number of mosquitoes as it occurs to ticks [19]. West Nile virus is the most detected Flavivirus in white-tailed deer [20] and feral pigs [21] in the USA, and also in Spanish wild ungulates [8,9]. Although the viremia they develop is insufficient for them to act as virus amplifiers, the serological reaction may be helpful to anticipate the occurrence of WNV outbreaks in areas of high risk of infection for susceptible hosts. No study has focused on unraveling WNV exposure dynamics determinants at small spatial scales. Thus, our goal was to use wild ungulates as WNV sentinels [8] and longitudinally study WNV exposure rate in the southwestern quarter of Spain to unravel virus dynamics in emerging areas of contrasting environmental favorability for its vectors. Wild ungulates could provide insights into the transmission ecology of WNV that would be helpful to progress in case prevention.

Study design
Five areas in southwestern Spain were selected as representative of the wild ungulate populations at risk of exposure to WNV in the country (Fig. 1). The selected areas were (1) "Doñana" National Park (A 1 ), (2) western "Sierra Morena" mountain chain (A 2 ), (3) central "Sierra Morena" mountain chain (A 3 ), (4) the "Guadiana River" valley (A 4 ), and (5) the "Montes de Toledo" mountain chain. Climatic and environmental conditions are remarkably diverse in the study areas, especially in A 1 if compared to the rest. In A 1 , a sub-humid thermo-Mediterranean climate with Atlantic influence predominates, where the average annual temperature is 17.6 • C and the average annual rainfall is 506.2 mm [22]. A 1 habitat includes large areas of marshland, sand dunes, and maritime pine (Pinus pinaster) woodland patches interspersed within Mediterranean scrubland. In A 2 , A 3 , A 4 , and A 5 , the predominant habitat consists of a combination of woodland areas of holm oak (Quercus ilex), cork oak (Q. suber), and gall oak (Q. faginea) combined with pine (Pinus spp.) forests and abundant Mediterranean scrub (Cystus spp., Rosmarinus spp., Erica spp. and Phillyrea spp.) interspersed with patches of grassland. Average annual temperatures are 17.4 • C, 15.3 • C, 15.0 • C and 14.5 • C in A 2 , A 3 , A 4 , and A 5 , respectively; the average annual rainfall for the 2010-2019 period [22] in these areas was 670.6 mm, 585.5 mm, 570.1 mm, and 520.9 mm, respectively. In A 4 and A 5 the predominant climate is continental Mediterranean.
Several studies have shown exposure of wild birds to flaviviruses in A 1 [23]. In both A 1 and A 2 , WNF cases have been consistently reported in horses and humans since 2010 [12], and Culex spp. vectors are abundant [24]. WNF was confirmed in 2014 in a horse farm in A 4 [25] and three positive birds were detected in 2007 in A 5 [11]. Further, WNF has been reported in horses and birds to the north and west of A 5 since 2014 [25]. Culex pipiens was recorded as the most abundant WNV vector by Durán-Martínez [15] in A 3 to A 5 followed in abundance by C. theileri. This evidence shows that the selected areas are appropriate for studying WNV emergence determinants in Spain.
We designed a longitudinal retrospective survey of WNV exposure in wild ungulates. We first calculated the sample size required to estimate annual seroprevalence per area. The estimation was based on the average antibody prevalence reported in wild ungulates in the region (4%) [8]. We adjusted the estimate for a large wild ungulate population (>2000 animals), with a precision of 7% and a confidence level of 95% [26]. Wild ungulate serum samples were collected in commercial hunting events annually in the study areas for 15 years. For ecological modeling purposes, the collection date was used to split samples into annual seasons (2005/2006 to 2019/2020). A season was defined as the period between April 1 in the year "t" and March 31 in the year "t + 1", so each season comprised the period of vector activity during which WNV can be transmitted [27]. In this way, we can ensure that a seropositive animal that was born in the spring of year "t" and sampled between April of that year and March of year "t + 1" was exposed to WNV during the mosquito activity season of its birth year.
Wild ungulates harvested by hunters, gamekeepers, or environment agents in commercial hunting events or population control trials were surveyed from 2005 to2020 (Table 1). Sampling was conducted in accordance with Spanish and EU regulations. Sera collected from 2418 animals (1168 red deer -Cervus elaphusand 1250 Eurasian wild boar -Sus scrofa) were selected. Data on spatial location, municipality, and survey date were recorded. Animals were sexed and their age was estimated through tooth eruption patterns [28] and split into (1) yearling (<1 year old), (2) juvenile (1-2 years old), and (3) adult (>2 years old) age classes.

Laboratory analyses
Multi-species blocking enzyme-linked immunosorbent assays (bELISA; INGEZIM West Nile COMPAC®, Ingenasa, Madrid, Spain) were performed to the selected sera for the detection of WNV antibodies following the manufacturer's recommendations. This commercial test is the most specific for WNV of all commercial WNV serological assays [29]. The bELISA displays high sensitivity (100%) and good specificity (79.5-96.5%) for anti-WNV anti-NPV antibodies [30].

Weather determinants
The population dynamics of WNV vectors, and thus of WNV transmission dynamics, is modulated by abiotic environmental conditions [31]. Unusual warm winters or rainy years could favor vector populations locally and maximize virus transmissibility. Therefore, stochastic variations in local weather conditions, rather than prevailing climatic conditions, were better suited to the objectives of our study. Meteorological parameters were collected from 2004 to 2020 from weather stations managed by the Spanish Meteorological Agency (AEMET). Different variables related to temperature and rainfall were estimated ( Table 2). Each individual animal was related to the estimated data from the spatially closest meteorological station. Annual mean values of the meteorological magnitudes were calculated from monthly records and taking into account the sampling season described above (April 1 to March 31).

Statistical analyses
To avoid common problems in data modeling, a comprehensive descriptive analysis was performed as recommended by Zuur et al. [32] and some highly correlated variables were excluded ( Table 2). All the continuous variables were rescaled prior to modeling to homogenize the scales of measures by applying natural logarithmic transformations. Once the exploratory analysis of the data was performed, the selected variables were included in all possible combinations of linear generalized mixed models considering survey area and year (season) as random factors. Models were ranked according to the corrected Akaike information criterion (AICc), using the 'dredge' function of the MuMin R package [33], and models with an AICc difference (ΔAICc) of less than two were included in the model averaging process to obtain the final best-fit model [see [34]]. We analyzed multicollinearity effects by estimating the variance inflation factor for model parameters [35]. A final step was to estimate the relative contribution of each model parameter to the variation of the probability of exposure. For this, we partitioned the coefficient of determination R 2 of the predictors included in the best model by estimating the semi-partial (part) R 2 and the  structural coefficients with the R package 'partR2' [36].

Descriptive findings
The sample size required per season and study area was 31, so we adjusted available samples in our sera bank to this number whenever possible. Sample selection was performed to balance sex and age classes in the study spatiotemporal scale unit (see Table 3). Six hundred and twenty-six of the 2418 animals (266 deer and 360 wild boar) were less than one year old, 735 (329 deer and 406 wild boar) were between one and two years old, and 957 (551 deer and 406 wild boar) were more than two years old. Sex was recorded in 2312 animals (1129 males and 1183 females), of which 1140 were red deer (574 males and 566 females) and 1172 were wild boar (555 males and 617 females).
In total, 603 of the 2418 animals were antibody positive ( Fig. 1). In all areas, seroprevalence values were higher in wild boar than in red deer (see Table 3).

Spatiotemporal trends in seroprevalence
Exposure to WNV was observed during the study period in all study areas (Fig. 2). A WNV transmission peak was observed in 2013/2014 with antibody prevalence reaching 43.2%. However, the temporal exposure dynamics was locally variable. In A 1 , two prevalence peaks were observed in 2013/2014 and in 2017/2018, both with rates above 60%. The prevalence peaks in A 2 were of similar levels to those in A 1 , with two peaks observed in 2009/2010 and 2019/2020 and a period of high prevalence between 2013/2014 and 2015/2016 with respect to the average values of the whole study period in this area. These peaks in A 1 and A 2 occurred mainly in years when numerous outbreaks in equids and human cases of WNF were declared in southwestern Spain [25] (Fig. 2). In several years we observed a temporal coincidence between incidence peaks in animals less than one year of age in A 1 and A 2 and the notification of WNF cases. In A 3 , several minor prevalence peaks (2010/ 2011, 2012/2013 and 2017/2018 to 2018/2019) were observed in the population with values around 30-40%. In A 4 , the temporal evolution of WNV exposure was more stable with a peak in prevalence and incidence between 2013/2014 and 2014/2015 that coincided with the notification of the only documented case of WNF in equids in this area. Finally, the exposure evolution scenario in A 5 was stable, with a peak prevalence around 40% in 2012/2013 to 2013/2014 in which we also observed a notable increase in estimated incidence with respect to the previous season, and a peak of lower prevalence (around 20%) between 2017/ 2018 and 2018/2019.
Increasing seroprevalence patterns were observed in the core study region and in areas A 1 , A 2 , and A 4 in contrast to the detrended exposure pattern observed in A 3 and A 5 . Disregarded local variations in exposure trends, seroprevalence local patterns displayed high disparity between contiguous years. Population seroprevalence patterns evolved over time in a manner very similar to the variation in the estimated incidence of new cases in animals <1 year of age (Fig. 2). Seasonal prevalence and incidence peaks did indeed temporarily match in most of the seasons and areas.

Exposure risk determinants
Ten models were best fitted to explain variation in WNV exposure probability (ΔAICc<2; Supplementary Table 1). The resulting average model included: 1) individual host predictors such as age and host species; and 2) weather predictors. Statistically significant host-realted predictors included age class and species. This revealed a higher probability of WNV exposure in >1 year-old animals relative to yearlings, and a higher risk in wild boar than in red deer (Fig. 3). Annual and summer accumulated rainfall were both also selected as relevant modulating predictors. An overall positive association was observed with increasing annual rainfall whereas exposure probability was higher in the driest summers ( Table 4). The largest risk weights for individual predictors were age and annual rainfall (Fig. 4). Although less relevant, temperatures were also selected as risk-modulating predictors in the best-fitted model.

Discussion
Our study reconfirms the presence of WNV antibodies among Iberian wild ungulate populations using a highly sensitive blocking ELISA test as reported by Boadella et al. [8], García Bocanegra et al. [9], and Caballero-Gómez et al. [37] in Spain. No study has previously dealt with unraveling the determinants of the small spatial scale transmission dynamics of WNV in a longitudinal approach, so our findings are novel and relevant to understanding and better predicting local transmission risks for improved prevention.
The observed seroprevalence is above the one previously reported in wild ungulates in Spain [8,9] and in other countries [20,38,39]. Considering that the first cases of WNF occurred in 2010, it is noteworthy that in all our 15 years of study, the seroprevalence was above 12%. The higher mean antibody prevalence observed in our study compared to those reported previously is probably related to the increasing environmental favorability for WNV vectors with time as the spread of WNV in the region indicates. In fact, Boadella et al. [8] observed a 10-fold increase in the seroprevalence in this region after the notification of the first WNF outbreaks in 2010, indicating the clear expansion of the virus that we confirm in this study and that the remarkable increase in WNF notifications throughout Spain corroborates [25]. The early detection in time of WNV antibodies in our study is consistent with the detection of antibodies after 2004 in wild birds [23]. It also agrees with reports from raptors in 2007 [11] and with WNVpositive mosquitoes in 2008 and 2009 [40]. Altogether, these findings confirm the active circulation of WNV in southern Spain well in advance Table 3 Summary of population (p) and yearling (y) specific bELISA results (No. of bELISA positive vs. the no. of sera analyzed, and the proportion of seropositive animals -in % -within brackets) throughout study area along the 15-year study period. Data are split for wild ungulates (WU: red deer + Eurasian wild boar), red deer (RD) and Eurasian wild boar (WB).  Fig. 2). This indicates that WNV may have been prevalent before outbreak reporting. Most probably, sporadic WNF cases went unnoticed by veterinary and medical practitioners as it occurred with Crimean-Congo hemorrhagic fever virus in the country [41]. Our results do not show a steep increasing trend in exposure as was observed with other emerging infections, e.g., bluetongue virus [42]. However, the time trend was positive in those areas where WNF outbreaks were reported (A 1 , A 2 , and A 4 ), which may indicate that changes in environmental conditions have favored the emergence of WNF. In fact, winter temperatures show an increasing trend over the study period in these areas (see Supplementary Fig. 1), which may have favored vector survival rates in winter (perhaps also of WNV avian reservoirs) and resulted in an increase in transmission. Boadella et al. [8] reported higher seroprevalence values in wild boar when compared to red deer. Several hypotheses could aid in explaining this finding that was also observed in our study. Wild boar prey on small animals and scavenge on the carcasses of various species [43], including birds. Several birds, e.g., waterfowl, corvids, or raptors, are highly sensitive to WNV infection [44] and their infected carcasses could be serving as wild boar food. Alternative hypotheses include potential differences in host selection preferences of WNV vectors, differences in space use between red deer and wild boar, or even differences in their daily activity patterns in relation to the host-seeking activity patterns of WNV vectors. The increasing pattern of exposure risk with age is in line with expectations because the risk of interacting with an infected mosquito increases with lifetime and WNV antibodies may persist for >5 years in mammals [45]. However, even though an expected long persistence of WNV antibodies, the evolution of seroprevalence in yearlings (indicators of exposure in the sampling year) and in the population (Fig. 2) overlapped, indicating that annual variations in exposure frequency occur in the overall population and that annual exposure rates can be considered as good indicators of infection incidence. Therefore, modeling the annual risk on the population in a particular year would be indicative of the real risk of exposure in this year. This could be, within an emerging scenario, the consequence of large interannual variations in environmental conditions that generate yearly variations in vector abundance together with the low average WNV antibody prevalence (when compared to other vector-pathogen models such as Hyalomma spp. ticks and Crimean-Congo hemorrhagic fever virus [46]).
The observed spatial seroprevalence pattern agrees with the highest number of outbreaks being reported in A 1 and A 2 [25]. This could be a consequence of differences in environmental conditions triggering variations to the burdens of mosquito vectors as well. In this regard, the diversity and abundance of Culex spp. are much lower in zones A 3 , A 4 and A 5 [15,47] where a dry continental climate predominates, compared to those reported in the two sub-humid thermo-Mediterranean climate zones with Atlantic influence, A 1 and A 2 [24,48]. This fact can be clearly linked to the positive influence of the accumulated rainfall on the risk of exposure to WNV, especially relevant in inland dry continental areas of Spain. The rainfall regime, unlike temperatures, did not show a clear trend during the study period in any of the areas (Supplementary Fig. 1), but a very changeable pattern between years as is usual in a Mediterranean climate. The rainfall filling semi-permanent or seasonal ponds with water is crucial for egg laying and development of mosquito larval stages, and thus a relevant parameter on their abundance [14], consequently resulting in a higher risk of WNV transmission. However, when water availability does not depend directly on the local rainfall regime but on that of its hydrographic basin, as is the case of A 1 , the dependence of these will be less in comparison with non-aquatic ecosystems such as those of the other four study areas. Although being a relevant modulating factor for vector populations, the model did not indicate a significant effect of temperature on risk even though different temperature-related parameters were selected within the best-fitted model. Differences between average winter and summer temperatures  Table 4 Output of average selected model including the selected predictors (Acronyms are shown in Table 2), the estimate and its associated standard error (SE), the statistic (z) and the p-value. The temporal patterns of variation in WNV exposure in our study neatly aligned with the record of equine and human outbreaks since 2010 in or near these areas. In 2010, 36 outbreaks in horses and two human cases were recorded in southwestern Spain [2,25], coinciding with the seroprevalence peaks observed in 2009/2010 and 2010/2011 in areas A 1 and A 2 . In 2013, 35 outbreaks in equines were recorded in the same area, and again we observed a marked increase in exposure in A 1 and A 2 . In area A 4 , the only equine outbreak was recorded in 2014, when we observed the highest seroprevalence level in the area. Further, between 2018 and 2020, we observed high seroprevalence levels in A 1 and A 2 . In 2020, 136 equine outbreaks and 77 human cases were reported in southern Spain [2,25]. This finding suggests that monitoring WNV exposure evolution in wild ungulates (commercially hunted and the subject of the Spanish national wildlife health surveillance program) could alert to the potential occurrence of WNF outbreaks. The combination of wild ungulate monitoring with wild bird active and passive surveillance, livestock monitoring, and mosquito surveillance would significantly improve our chances to predict the risk of WNF emergence in specific areas of Spain and thus reduce its impact. Our goal was not to predict spatiotemporal risks, but our findings highlight the usefulness of building predictive models based on sentinels to improving preventive capabilities against this disease.

Conclusions
West Nile virus was circulating in Spain long before the emergence of WNF. The main determinants of its transmission, excluding host-related predictors, are environmental and probably mostly related to vector abundance. We reconfirm the usefulness of wild ungulates as sentinels to monitor WNV transmission dynamics, but recommend that an integrated WNV surveillance program including other wildlife species, domestic animals and vectors that are relevant determinants of WNV Fig. 4. Forest plots for comparison of part R 2 (coefficient of determination) for model predictors (A), inclusive R 2 (B), structure coefficients (C) and beta weights (C) including confidence intervals (CI) for the general population risk model. The acronyms for each predictor are included as described in Table 2. ecology is designed and implemented in the framework of an expected increase in WNF incidence in Spain in the near future.

Funding
This study was funded by the Spanish Ministry for the Science and Innovation (MCI) through the Spanish Research Agency (AEI) and EU FEDER (E-RTA2015-0002-C02 & CGL2017-89866-R grants), and by the Regional Government of Castilla-La Mancha (JCCM) and the European Social

Ethics statement
The study was conducted on red deer and wild boar that were shot by hunters during commercial hunting events or by environment agents as part of population control measures. Therefore, they were not shot deliberately for this study and no ethical permission was required to gather samples.

Declaration of Competing Interest
The authors declare no conflict of interest.

Data availability
Data will be made available on request.