Red deer reveal spatial risks of Crimean-Congo haemorrhagic fever virus infection

Received: 15 July 2021 Revised: 24 September 2021 Accepted: 28 September 2021 DOI: 10.1111/tbed.14385 ORIGINAL ARTICLE Red deer reveal spatial risks of Crimean-Congo haemorrhagic fever virus infection Raúl Cuadrado-Matías1Beatriz Cardoso1,2Miriam A. Sas3 Ignacio García-Bocanegra4Isolde Schuster3David González-Barrio1,5 Sven Reiche6Marc Mertens3David Cano-Terriza4Laia Casades-Martí1 Saúl Jiménez-Ruiz1,4Jordi Martínez-Guijosa1Yolanda Fierro7 Félix Gómez-Guillamón8Christian Gortázar1Pelayo Acevedo1 Martin H. Groschup3Francisco Ruiz-Fons1 1Health & Biotechnology (SaBio) Group, Instituto de Investigación en Recursos Cinegéticos (CSIC-UCLM-JCCM), Ciudad Real, Spain 2CIBIO/InBio, Centro de Investigação em Biodiversidade e Recursos Genéticos, Universidade do Porto, Porto, Portugal 3Institute of Novel and Emerging Infectious Diseases, Friedrich-Loeffler-Institut, Insel Riems, Germany 4Grupo de Investigación en Sanidad Animal y Zoonosis (GISAZ), Dpto. de Sanidad Animal, Universidad de Córdoba, Córdoba, Spain 5Parasitology Reference and Research Laboratory, Spanish National Centre for Microbiology, Madrid, Spain 6Department of Experimental Animal Facilities and Biorisk Management, Friedrich-Loeffler-Institut, Insel Riems, Germany 7Yolfie Properties, Ciudad Real, Spain 8Consejería de Agricultura, Ganadería, Pesca y Desarrollo Sostenible, Junta de Andalucía, Málaga, Spain Correspondence FranciscoRuiz-Fons,Health &Biotechnology (SaBio)Group,Instituto deInvestigaciónen RecursosCinegéticos (IREC),Rondade Toledo 12,13071 CiudadReal,Spain. Email:[email protected] Funding information FCT-Fundaçãoparaa Ciênciae Tecnologia,Grant/Award Number: 2020.04872.BD;Junta deComunidades deCastilla-LaMancha, Grant/Award Numbers:SBPLY/19/180501/000321, 2018/12504,BES-2015-072206; EuropeanSocial Fund,Grant/AwardNumber:SBPLY/19/180501/000321;Universidad deCastilla-La Mancha,Grant/Award Numbers:PRE2018-083801, PEJ2018-003155-A, 2018/12504;Instituto deSaludCarlos III, Grant/Award Number:CD19CIII/00011; SpanishMinistry fortheScience andInnovation,Grant/AwardNumbers:CGL2017-89866R,E-RTA2015-0002-C02-02, PRE2018083801,PEJ2018-003155-A,BES-2015072206 Abstract Crimean-Congo haemorrhagic fever virus (CCHFV) continues to cause new human cases in Iberia while its spatial distribution and ecological determinants remain unknown. The virus remains active in a silent tick-animal cycle to which animals contribute maintaining the tick populations and the virus itself. Wild ungulates, in particular red deer, are essential hosts for Hyalomma ticks in Iberia, which are the principal competent vector of CCHFV. Red deer could be an excellent model to understand the ecological determinants of CCHFV as well as to predict infection risks for humans because it is large, gregarious, abundant and the principal host for Hyalomma lusitanicum. We designed a cross-sectional study, analysed the presence of CCHFV antibodies in 1444 deer from 82 populations, and statistically modelled exposure risk with host and environmental predictors. The best-fitted statistical model was projected for peninsular Spain to map infection risks. Fifty out of 82 deer populations were seropositive, with individual population prevalence as high as 88%. The highest prevalence of exposure to CCHFV occurred in the southwest of the Iberian Peninsula. Climate and ungulate abundance were the most influential predictors of the risk of exposure to the virus. The highest risk regions were those where H. lusitanicum is most abundant. Eight This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Transboundary and Emerging Diseases published by Wiley-VCH GmbH Transbound Emerg Dis. 2021;1–16. wileyonlinelibrary.com/journal/tbed 1 2CUADRADO-MATÍAS ET AL. ofthenineprimaryhumancasesoccurredinorborderingtheseregions,demonstrating that the model predicts human infection risk accurately. A recent human case of CCHF occurred in northwestern Spain, a region that the model predicted as low risk, pointing out that it needs improvement to capture all determinants of the CCHFV infection risk. In this study, we have been able to identify the main ecological determinants of CCHFV, and we have also managed to create an accurate model to assess the risk of CCHFV infection. KEYWORDS epidemiology, orthonairovirus, prevention, risk mapping, tick, zoonosis 1INTRODUCTION The causes underlying the spatial spread of vector-borne pathogens (VBP) have been investigated in some highly relevant model VBP of human, agriculture or conservation impact, for example Borrelia burgdorferi s.l., tick-borne encephalitis virus, West Nile virus or bluetongue virus (Jacquot et al., 2017; Jaenson et al., 2012; Kilpatrick, 2011; Mannelli et al., 2012). These examples show that VBP spread on their (flying) vectors and hosts in search of new suitable ecosystems.Sometimes, VBP find favourableconditions tocolonizenew areas by adapting to local autochthonous vertebrates and invertebrates (Kilpatrick, 2011). VBP can also spread on their competent vectors and hosts to the neighbourhood (Leighton et al., 2012). If pathogens are plastic enough to adapt to new (host and vector) environments, their chances to colonize and spread to new vast unexplored rich niches are high. Ticks cannot fly, but some manage to take a free and comfortable flight that transports them to distant places. Ticks exploit migratory birds that fly across large distances in a short period between reproduction and wintering headquarters (Estrada-Peña et al., 2021). Resident hosts also spread ticks at shorter distances (Buczek et al., 2020; Ruiz-Fons & Gilbert, 2010). The colonization of new ecosystems by VBP may be silent if pathogens spread without attracting the attention of health/conservation authorities or researchers. Crimean-Congo haemorrhagic fever virus (CCHFV), a tick-borne Orthonairovirus of African origin and currently endemic to the African, Asian and (southeastern) European continents, may have been flying in infected ticks on migratory birds for thousands of years across the western Mediterranean Basin (Palomar et al., 2013). The virus might also have passively exploited the movements of their vectors in internationallytradedlivestock(Muhanguzietal.,2020) as well asinanthropogenic wildlife translocations (Tsao et al., 2021) to reach southwestern Europe. It is not known how long CCHFV has been circulating in southwestern Europe. There was no evidence of the virus in Iberia until 1985 when CCHFV antibodies appeared in two humans during a serological survey for tick-borne pathogens in Portugal (Filipe et al., 1985). In 2010, the virus was detected in Hyalomma lusitanicum ticks collected on red deer (Cervus elaphus) in west-central Spain (Estrada-Peña et al., 2012), and later human clinical cases emerged in 2016, 2018, 2020 and 2021 in western Spain. Recently, the first known human clinical case of CCHF in Spain has been traced back to 2013 (Negredo et al., 2021). Ten human cases have been reported to date in the country, three of whom resulted in death (30% fatality rate): one case in 2013, one case in 2016, two cases in 2018, three cases in 2020 and two cases in 2021. No human clinical case has been reported in Portugal to date. The primary cases reported had clinical symptoms of haemorrhagic fever, but these may represent only a small proportion of the real infection cases that occur in Spain and Portugal annually. This perception relies on the high proportion of infections that may go unnoticed due to mild, nonspecific clinical signs (Bente et al., 2013). Additionally, a proportion of non-fatal cases of haemorrhagic fever may go undiagnosed, for example one of the primary cases reported in 2018 was diagnosed during a retrospective survey on fever cases of unknown origin in west-central Spain, and the 2013 case was also a patient with high fever of undiagnosed aetiology. The virus has been detected in several areas of Iberia, mostly in H. lusitanicum ticks collected (mainly) on red deer (Cajimat et al., 2017; Negredo et al., 2019). A study undertaken by the Spanish Ministry for Agriculture (MSCBS, 2019) revealed higher antibody prevalence in wildlife when compared to livestock in specific areas of Spain, suggesting that the virus was widespread in a sylvatic cycle on the mainland. Recently, CCHFV-positive ticks were found in new areas in southcentral and southern Spain. This observation provides more evidence of an enzootic but not homogeneous circulation of the virus in Iberia (Moraga-Fernández et al., 2020). PrimordialtickvectorsofCCHFVarespeciesinthegenus Hyalomma (Spengler & Bente, 2017). Two Hyalomma ticks are well established in Iberia, H. lusitanicum and H. marginatum. The former predominates in abundance over H. marginatum in southwestern Iberia, where nine of the ten human cases occurred. This region of Iberia has abundant populations of red deer and other wild ungulates (Acevedo et al., 2007, 2008) that maintain high burdens of H. lusitanicum ticks (Ruiz-Fons et al., 2006, 2013; Valcárcel et al., 2016). The red deer is, notably more thanother wild ungulates,the primaryhost of H. lusitanicum (andother) ticks in this region, and it is also widely distributed in Iberia (Bencatel et al., 2019; Palomo et al., 2007;Figure1). Domestic ruminants are also abundant in this region, and they are managed frequently under extensive production systems in large extensions of land. Domestic ruminants interact with coexisting wildlife, so H. lusitanicum and CUADRADO-MATÍAS ET AL.3 FIGURE 1 Spatial location of the surveyed red deer populations (red dots) and local sample size in relation to red deer (blue shaded squares) distribution in the Iberian Peninsula (Bencatel et al., 2019; Palomo et al., 2007) H. marginatum ticks are frequent in extensively raised domestic ruminants even though acaricides are spread over animals or administered in feed at different times of the year within the main tick activity period (spring to autumn). CCHFV infection causes no disease in ungulates and other mammals. However, infected animals may either replicate and transmit the virus to uninfected ticks, allow co-feeding transmission among infected and uninfected ticks, ease venereal transmission between adult ticks at mating or favour virus maintenance in the tick population by feedingticks (Estrada-Peña et al., 2013a;Ruiz-Fonset al., 2012).The virusmaybe maintained in the tickpopulation bytransovarial transmission. Direct exposure to ungulate carcass fluids may also be a relevant transmission pathway at the animal–human interface (Shahhosseini et al., 2018). However, in the EU, there is no evidence of transmission from domestic or wild ungulates to abattoir or game meat professionals during carcass dressing (ECDC, 2020). No vaccine is available to protect humans against CCHFV, so preventionisthe only measure for avoidingnew infections.Preventingtick bites include a series of recommendations to people in contact with domestic or wild animals or with the environment where they live, such as wearing appropriate clothes, carrying out a thorough inspection of the skin to remove ticks after any field activity, or using repellents against ticks. However, informing about high-risk areas would significantly increase awareness where it is more needed and result in better prevention strategies. In Iberia, the spatial distribution of CCHFV is currently unknown, even though the virus has a huge potential to spread as it is occurring with West Nile virus in this region (ECDC, 2020). Attempts to map the distribution of CCHFV in the human population have been unfruitful because humans do not frequently become exposed to tick bites (Monsalve Arteaga et al., 2020). Farmers often treat extensively produced livestock with acaricides that reduce exposure to CCHFV. In contrast, wild ungulates host high amounts of (Hyalomma) ticks, are infrequently treated against them and are widely distributed in Iberia (Bencatel et al., 2019; Palomo et al., 2007). Red deermeetmost of therequisitesthat would makethemgood indicators of the risk of CCHFV infection in Iberia: (i) they are abundant; (ii) they are gregarious; (iii) they are widely distributed and (iv) they host high numbersofHyalomma spp.(andother)ticks.Wethus hypothesizedthat theintimaterelationship ofreddeerand Hyalomma ticksin Iberiamight aid in mapping the risk of CCHFV exposure to inform public health authorities and the public as a preventative measure and thus reduce human CCHF cases. Understanding the fundamentals of the ecological background of the enzootic cycle of CCHFV would additionally result in insights for the future control of this emerging zoonosis in Iberia. 2MATERIALS AND METHODS 2.1 Survey design The study focused on the Iberian Peninsula, a 596,740 km2land heterogeneous in climatic, orographic, ecosystemic, and socioeconomic terms. We designed a cross-sectional survey based upon the hypothesis that red deer are exposed to bites from CCHFV infected ticks at a higher rate than other Iberian ungulates (Ruiz-Fons et al., 2013) and being abundant and widely distributed in the region as well, they would provide a realistic map of the risk of infection by CCHFV for humans. That would render red deer as exceptional indicators for the spatial distribution of CCHFV. Therefore, this would be very useful 4CUADRADO-MATÍAS ET AL. for informing the health authorities about the areas where the risk of exposure to CCHFV is higher so that preventative measures may be taken. To achieve this, we chose to estimate the rate of exposure of red deer to CCHFV by detecting the presence of specific antibodies in blood serum. The study design needed to be based on a representative sample of red deer populations in Iberia, so we checked the distribution range of the red deer (Bencatel et al., 2019;Palomoetal., 2007). The unit of study was the epidemiological population. For this study, we defined an epidemiological population as the group of red deer individuals inhabiting a specific territory under the management ofasingle authority (local/regional administration,hunters’ association or landowner). With this classification, the effect of the set of specific management measures to which deer were exposed locally could be homogenized, which depends exclusively on the goals of the local manager. These may vary from those of neighbouring populations subject to the decisions of another manager. We chose this classification system because many of the study populations are artificially restricted to a particular territory (range around 300 to 12,000 has) by large game fences (Acevedo et al., 2007, 2008). We calculated the minimum number of samples required to estimate antibody prevalence in the studyunits atthe previously known circulationratesinwestern Europe (Spengler et al., 2016a) with a 95% confidence level and an accepted 10% error using the proportion calculator of Epitools (Sergeant, 2018). Whendesigning the cross-sectionalsurvey,few serological surveyshad been conducted on CCHFV in domestic ungulates in western Europe (Spengler et al., 2016a). Reported antibody prevalence did not exceed 2%, so we estimated the required sample size for a 2% expected proportion. We selected serum samples of red deer collected between 2008 and 2016. When gathering samples, we also collected data on the surveyed location and individuals. We recorded the geographical coordinates of the surveyed sites with portable GPS devices. All the samples were collected from hunter-harvested red deer shot during commercial/social hunting events, or after official population control events carried out by environment agents in protected areas. We performed the sampling according to Spanish and EU regulations. We did not require any ethical approval from authorities because we did not shoot animals deliberately for the survey. 2.2 Serological analyses The presence of specific CCHFV antibodies in serum samples was estimated using a species-independent in-house competitive ELISA (cELISA) developed at the Friedrich-Loeffler Institute in Germany (Schuster et al., 2016). This cELISA uses recombinant nucleocapsid (N) protein from the CCHFV strain Kosovo Hoti (GenBank accession no. DQ133507). This protein was expressed in Escherichia coli and purified by nickelchelate affinity chromatography (Qiagen, Hilden, Germany). A mixture of 1:1 N protein and glycerin was incubated at 4◦C for 1 h. Two micrograms of this mixture were diluted in a solution of PBS and 10% glycerin and used to coat 96-well plates (PolySorp immunoplates, Nunc, Roskilde, Denmark). The sealed plates were incubated overnight at 4◦C. Two hundred microlitres of blocking buffer (3% BSA in PBS-Tween20 0.05%) were added per well and incubated at 37◦C for 1 h. Finally, the plates were washed 4 times with PBS-Tween20 1% (250 µl/well) before use. Test serum samples were heat-inactivated at 56◦C for 30 s. Each serum was diluted 1:4 in a solution of 1% BSA in PBS-Tween20 0.1%. Ninety-five microlitres of the diluted samples and controls were added in duplicate to the corresponding wells of 96-well plates and later incubated for 1 h at 37◦C. Five microlitres of monoclonal antibodies (mAbs) (obtained after immunization of BALB/c mice with the CCHFV N protein; see Schuster et al., 2016), diluted 1:50 in a 1% BSA in PBS-Tween20 0.5% solution, were added to the wells containing the samples and controls. After incubation for 10 minutes at RT, 100 µl/well of a H2SO4 solution were added to stop the reaction. The plates were incubated for 2 h at 37◦C, washed and to each well was added 100 µl of TMB solution (Bio-Rad, Munich, Germany). After incubation for 10 s at RT, 100 µl/well of a H2SO4 solution were added to stop the reaction. The optical density (OD) of the wells was read at a 450 nm/620 nm wavelength in an automatic Multiskan SkyHigh Microplate Spectrophotometer (Thermo Fisher Scientific, Waltham, USA). The results of the spectrophotometry were used to estimate the percent inhibition (PI) of each sample compared to the negative controls (NC): PI =100 −(sample-OD/NC-OD) ×100. Confirmed CCHFV seropositive cattle serum and seronegative confirmed sheep serum were employed as positive and negative controls, respectively. Test performance was estimated according to variations in plate type, coating conditions, amount of the antigen, incubation temperatures and times, buffers, washing times, serum/mAB/conjugate dilutions, substrate incubation time and application of mABs. The cELISA was validated by comparative analysis using samples from several animal species (domestic and wild) and humans, as described in Schuster et al. (2016). Briefly, the validation was carried out using animal sera (cattle, sheep, goat) from enzootic countries in southeastern Europe (previously classified as antibody-positive using adapted commercial serological screening tests for CCHFV) and sera of experimentally infected rhesus macaques as positive controls. Animal and huma sera from a non-enzootic country (Germany) were employed as negative controls. The cut-off point for maximizing the diagnostic sensitivity and specificity of the cELISA was determined by the receiving operator characteristics (ROC) analysis. PI values above 49 were considered positive, and those below 37 were considered negative. Those between 37 and 49 were considered inconclusive. The estimated diagnostic sensitivity for the competitive ELISA was 95%, and its diagnostic specificity reached 99% (Schuster et al., 2016) at the established cutoff point. The sera with homogeneous results in the duplicate analysis (both positive or both negative) were classified as positive or negative, respectively. Sera with contrasting results or with two inconclusive results were retested to classify them as negative or positive. The prevalence of antibodies was estimated as the ratio between the number of positive and analysed sera in percentage terms. We controlled uncertainty in prevalence values by associating the exact Clopper– Pearson exact 95% confidence interval (CI) to each prevalence value. CUADRADO-MATÍAS ET AL.5 2.3 Spatial cluster analysis We initially explored the spatial dependence of antibody prevalence in the study populations to identify the highest risk areas for red deer in Iberia. This was carried out by implementing a spatial cluster analysis with SaTScanTM v9.6 software (Kulldorf, 2018) to identify both high and low relative risk areas (RR), RR >1andRR<1, respectively. The Bernoulli model was employed for the cluster analysis (Coleman et al., 2009; Kulldorff, 1999) without including temporal parameters due to the cross-sectional nature of the survey. We set a circular spatial windowfortheclusterswith a maximum populationsizeatriskof50% with no overlapping between neighbouring clusters. Analyses were run with 9999 replications. The pvalue was estimated with the Gumbel approximation to infer the significance of cluster RR. Clusters were considered significant at p<.05. The clusters were hierarchically numbered and organized according to the pvalue. 2.4 Risk factor analysis To understand the factors determining variations in the probability of exposure to CCHFV, we undertook statistical modelling with a series of explanatory variables selected from the host individual and host population factors (Table 1). We gathered environmental variables as well because these may modulate exposure to CCHFV by influencing host and tick population dynamics (Ruiz-Fons et al., 2012, 2013). This approach was aimed at providing an overview of the factors determining exposure at the population level (as defined in this study) that could help design future strategies for reducing CCHFV transmission at the host-tick interface. Two purely spatial, nine climatic, two topographical, seven habitat, six host population and two individual host variables were initially selected based on their potential to modulate exposure to CCHFV tick vectors (see Appendix 1) and on availability at the study scale (see Acevedo et al., 2010). A large portion of the host population predictors and all environmental predictors were estimated at UTM 10 ×10 km spatial scale to cover the whole range of study populations (300–12,000 has). Environmental and host population predictors were only available in Spain. Therefore, any of the four surveyed Portuguese populations close to the Spanish border (Figure 1) were linked with data from the closest Spanish UTM 10 ×10 km square and that in central Portugal was discarded for statistical analyses. We thoroughly checked data to rule out any potential interference in statistical modelling (Zuur et al., 2010).Thisinitialstep enabled potentialoutliervaluesinthe predictors to be evaluated and controlled, that, when identified, were ruled out by logarithmic transformation. Logarithmic transformation was applied to all the continuous climatic predictors selected for modelling to homogenize the range of scales in the measures. A Pearson correlation plot wasbuilt withcontinuous predictorsusing the‘corrplot’package ofR in RStudio (Wei & Simko, 2017). Any collinearity (r≥|0.7|) was removed whenselectingthepredictorsformodelling. Wecheckedforanypotential dependence in the predictors and that sample size was balanced among classes in categorical variables. We finally checked for potential meaningful interactions between predictors. We selected two that could modulate host–tick interaction patterns: (i) deer sex and ageclass interaction because previous findings show the combined effect of sex and age driving the amount of H. lusitanicum ticks on red deer (Ruiz-Fons et al., 2013) and (ii) slope and soil permeability interaction because steep slope terrains may counteract the water retention potential of low permeable soils by runoff effects whereas highly permeable soils on flat terrains may lose water through drainage, which affects tick survival and abundance, and hence CCHFV transmission to deer. Thereafter, the individual risk of exposure to CCHFV (‘ecchfv’; positive/negative; N=1247) was modelled with the selected covariates in a logistic regression modelling approach [ecchfv ∼mg +sex +age +sex:age +lg(mxtMO) +lg(ts) +lg(ndvis) +lg(ndvi) +lg(sl):lg(sp) + dfav +bfav +lg(ld)] (McCulloch et al., 2008). All possible models were built and ranked by increasing corrected Akaike information criterion (AICc) with the ‘dredge’ function of the R ‘MuMIn’ statistical package (Barton, 2009). Those models displaying a difference of under 2 units in the AICc, focusing on the model with the lowest AICc (ΔAICc <2; Supplementary Table S1), were selected for model averaging with the ‘model.avg’ function of the ‘MuMIn’ R package. A subset of deer populations located within the main areas of distribution of H. lusitanicum in Iberia (southwestern Iberian quarter; Figure 1; Estrada-Peña et al., 2013b) was selected to seek for predictors modulating CCHFV exposure in the presence of one of its main vectors, thus replicating the risk modelling approach detailed above. We checked potential autocorrelational effects on model coefficients for the two model sets by estimating the generalized variance-inflation factors (GVIF) for model parameters (‘vif’ function of the ‘car’ R package; Fox & Weisberg, 2019), and we controlled this in the selected models. Finally, the relative proportions of the variation in the CCHFV exposure risk of selected model predictors – grouped into environmental (spatial, climatic, topographic and habitat predictors), host population and host individual factors (Table 1) – were calculated by variation partitioning (‘varPart’ function) using the ‘modEvA’ R package (Barbosa et al., 2013; see also Real et al., 2003). 2.5 Risk mapping To map the risk of exposure to CCHFV using red deer as an indicator of the transmission risk from infected ticks, a new model for mainland Spain was built, but only with (host) population and environmental predictors. That approach assumes that red deer sex and age population structures and the diversity of management measures are homogeneous across Spain. Therefore, there are similar effects on vector and virus dynamics across its distribution range. The model-building approach was as described previously. Before mapping the predicted exposure risk, we analysed the ability of the logistic regression model to make accurate predictions and thereby prevent (control) biased predictions (see Gude et al., 2009). The internal validation of the (average) logistic regression model was performed to estimate prediction error (accuracy) by internal bootstrap validation with the ‘boot’ R package 6CUADRADO-MATÍAS ET AL. TABLE 1 Set of explanatory predictors gathered for risk factor analysis Factor Predictor Description (unit of measure) Variable type Unit scale Average (range) Spatial lat Latitude (decimal degrees) Continuous Population 38.70 (36.25–43.26) long Longitude (decimal degrees) Continuous Population −4.57 (−8.37 to −5.34) Climatic mtMO Mean temperature May–October (◦C) Continuous UTM 10 ×10 km 20.4 (12.0–23.3) mntMO Mean of minimum temperatures May–October (◦C) Continuous UTM 10 ×10 km 14.9 (5.7–9.2) mxtMO Mean of maximum temperatures May–October (◦C) Continuous UTM 10 ×10 km 27.7 (18.1–30.6) ts Temperature seasonality (◦C) Continuous UTM 10 ×10 km 4.4 (2.7–8.5) apr Mean annual precipitation (mm) Continuous UTM 10 ×10 km 667.8 (367.6–1318.9) prMO Mean precipitation May–October (mm) Continuous UTM 10 ×10 km 164.2 (102.3–439.5) prs Precipitation seasonality (mm) Continuous UTM 10 ×10 km 64.2 (27.7–89.1) arad Mean annual radiation (Kwh/m2/day) Continuous UTM 10 ×10 km 2027.9 (1896.8–2063.8) acrad Accumulated annual radiation (Kwh/m2/day) Continuous UTM 10 ×10 km 24334.1 (22761.1–24765.0) Topographic alt Mean altitude (masl) Continuous UTM 10 ×10 km 582.7 (21.0–1783.8) sl Slope (degrees) Continuous UTM 10 ×10 km 4.3 (0.3–12.6) Habitat ndvi Mean normalized difference vegetation index Continuous UTM 10 ×10 km 145.5 (100.3–196.8) ndvis NDVI seasonality Continuous UTM 10 ×10 km 10.3 (1.9–33.7) sp Soil permeability Continuous UTM 10 ×10 km 1.6 (1.0–3.0) wp Woodland (%) Continuous UTM 10 ×10 km 35.2(1.5–78.0) scrp Scrublands (%) Continuous UTM 10 ×10 km 31.1 (4.7–71.9) svp Scarcely vegetated areas (%) Continuous UTM 10 ×10 km 0.7 (0.0–28.6) irrp Irrigated land (%) Continuous UTM 10 ×10 km 2.0 (0.0–68.4) Host population ld Livestock density (ind/Ha) Continuous UTM 10 ×10 km 0.55 (0.00–1.61) cd Cattle density (ind/Ha) Continuous UTM 10 ×10 km 0.11 (0.00–0.48) dfav Environmental favourability for red deer Continuous UTM 10 ×10 km 0.62 (0.14–0.92) wbfav Environmental favourability for wild boar Continuous UTM 10 ×10 km 0.48 (0.07–0.95) rfav Environmental favourability for roe deer Continuous UTM 10 ×10 km 0.37 (0.05–0.95) mg Estate management practices Categorical Population 1−3 Host individual sex Red deer sex class Categorical Individual Male/female age Age class of individual red deer Categorical Individual Yearling/juvenile/sub-adult/ adult Note: Variables selected for modelling are marked in bold italics. (Efron & Tibshirani, 1994) and by split-sample validation with the ‘caret’ R package (Kuhn, 2008). The logit equation from the output of the average model enabled the predicted probability of exposure to be estimated for every UTM 10 ×10 km square on the Spanish mainland. Those probabilities were represented at the UTM 10 ×10 km square level and mapped using ArcMap 10.5 (ESRI, Redlands, CA, USA) software. 3RESULTS The minimum largest required sample size per population was eight individuals for a 10% precision estimate in a large population (>400 individuals). The study was based on the retrospective analysis of red deer serum samples collected within the framework of other research projects. Therefore, we adapted our survey to the serum banks CUADRADO-MATÍAS ET AL.7 FIGURE 2 Results of the serological survey of CCHFV in terms of the epidemiological population in a map displaying the distribution of the red deer at the UTM10 ×10 km square unit level in the Iberian Peninsula (blue shaded squares; Bencatel et al., 2019; Palomo et al., 2007) available. That resulted in 1444 serum samples collected in 82 populations, 77 in Spain and 5 in Portugal, between 2008 and 2016 (Figure 1). The serological analyses identified that 367 sera had antibodies against CCHFV (367/1444; 25.4%, CI: 23.2–27.7). Positive samples were spatially widespread in Iberia (Figure 2).Fiftyofthe82surveyed populations (61%) had at least one positive sample. Prevalence values ranged from 0% to 88%. Populations with ≥10% antibody prevalence were all in the southwestern quarter of Iberia, except for a single population in eastcentral Spain (Figure 2). The rest, either displayed low antibody prevalence (<10%; n=9, 11%) or no positive sample (n=20, 24.4%; 0%), were distributed across Iberia. The spatial cluster analysis identified three statistically significant clusters with an increased risk of exposure to CCHFV and four additional low RR clusters (Figure 3). High RR clusters included (i) cluster 2 in the ‘Toledo Mounts’ system (radius: 17.4 km; RR: 3.32) in southcentral Spain; (ii) cluster 3 on the southernmost tip of mainland Spain (radius: 144.1 km; RR: 2.34) and (iii) cluster 5 in inland east-central Spain (radius: 0 km; RR: 3.03). Meanwhile, clusters 1 (radius: 247.4 km), 4 (104.7 km), 6 and 7 (single populations) represented areas of low RR of exposure to CCHFV in northeastern, southeastern and specific isolated red deer populations in southwestern Iberia, respectively. These results reflect that southwestern Iberia is a hotspot for the transmission of CCHFV. However, there are significant local variations within the large-scale risk area. The fundamental average environmental and hostpopulation parameters onlyslightly differed (not statistically; data not shown) between clusters of high or low relative risk except for the seasonality in the NDVI (Table 2). Populations within clusters of low RR were more seasonal in vegetation productivity than high RR areas. We built the model for the whole set of red deer populations surveyed in Iberia with 1247 samples; we excluded 197 samples due to missing data. Five models were within the ΔAICc <2 range established for model averaging (Supplementary Table S1). Average model calibration was appropriate according to the Hosmer–Lemeshow test (p=.105). The model also had good discriminatory power (AUC =0.725). Host individual/population parameters and environmental/topographic covariates were significant predictors for the risk of exposure of individual red deer to CCHFV in Iberia. We observed no gender effects, but a decreasing risk was observed with the increasing age class. Yearlings were at a higher risk of exposure to CCHFV than older individuals. The host individual factor explained 4.6% of the variation in virus exposure risk (Supplementary Figure S1). As expected, host population parameters played an important role in modulating the risk of exposure to CCHFV, including wild ungulate and livestock demographic traits and deer management. This factor explained the 8.6% risk variation. The increasing environmental favourability for the red deer and the high livestock density positively affected the risk of exposure to the virus. However, areas more favourable in terms of the environment for the Eurasian wild boar had a negative (not statistically significant) effect (Table 3). Red deer from populations being managed for hunting or farming were more prone to exposure to CCHFV than those in open, loosely managed ones. Environmental covariates also had a significant role in modulating the risk of exposure to the virus, either by favouring (NDVI, soil permeability or temperature) or by preventing (NDVI seasonality, terrain slope or temperature seasonality) exposure. The environmental factor explained the highest proportion of variation in the risk of exposure of deer to CCHFV itself (39.7%). The aggregate effects of host individual 8CUADRADO-MATÍAS ET AL. FIGURE 3 Spatial range of CCHFV high (red shadowed) and low (blue shadowed) relative risk (RR) prevalence clusters in relation to surveyed red deer populations. Populations included in any of the low or high RR identified clusters are shown distinctly to red deer populations not included in any cluster (white dots). The numbers on the right upper side of each cluster correspond to cluster numbers as shown in Table 2 and population factors explained the 17.7% variation in risk. However, the highest explained variation in proportion was achieved by the sum effects of (host) population and environmental predictors (64.7%; Supplementary Figure S2). Risk modelling within the main spatial range of H. lusitanicum in Iberia resulted in six models displaying AICc differences below the two established units (Supplementary Table S1). We built the model with 1025 samples from 63 deer populations in southwestern Iberia. Calibration of the average model was poor according to the Hosmer– Lemeshow test (p<.05), albeit good discriminatory power was maintained (AUC =0.717). In this core area, (host) individual traits were less significant predictors of virus exposure risk (Table 4). We included age and sex for model averaging, but their effect was not statistically significant. The (host) individual factor could indeed explain only 2.2% of the variation in exposure risk (Supplementary Figure S2). Host population and environmental variables were good predictors on this spatial scale. The host population factor explained 8.5% of variation in the model. The environment factor accounted for 52.5% of it. Jointly, these two factors could explain 66.3% of variation in the model. The index of environmental favourability for the red deer lost power with the regional model (probably because it was on a national basis and calculated with presence/absence – not abundance – data). However, livestock density (a more reliable estimate of local densities calculated with census data collected at the local veterinary unit level) still displayed a significant positive relationship to exposure risk. Deer farms in this region had a higher risk of exposure when compared to deer farming on the scale of mainland Spain. On this regional scale, the NDVI had the effect of boosting virus exposure to deer. In contrast, on the premise of an expected higher vegetation production with increasing annual rainfall, annual precipitation on this regional scale had a slight albeit significant negative effect on exposure risk. In contrast, seasonality in precipitation displayed positive influences. Finally, we observed a negative joint effect of the slope/soil permeability interaction on infection risk. The model built to predict exposure risk in mainland Spain corroborated the observed effects of (host) population and environmental factors with the peninsular scale model (Table 5). That is, the positive effects of deer favourability index, livestock density and soil permeability combined with the negative influence of wild boar favourability index, NDVI and temperature seasonality. The model enabled predicted risk at the UTM 10 ×10 km spatial scale on mainland Spain (Figure 4) to be represented. The predictive accuracy of this model, estimated by internal validation, showed a bootstrapped AUC of 0.714 that matched the apparent AUC (0.704) obtained from the bootstrap sample (R=200). Predictive accuracy was estimated with a split-sample procedure using 50% of the data as a training set and the remaining 50% as a test dataset. Estimated accuracy was 74%. Therefore, the model could not classify 26% of the samples correctly. CUADRADO-MATÍAS ET AL.9 TABLE 2 Descriptive parameters for significant predictors of the risk of exposure to CCHFV (identified by the selected overall risk model) of the high and low relative risk (RR) spatial clusters identified as statistically distinct High RR clusters Low RR clusters Parameter Cluster 2 (RR: 3.32) n=2 Cluster 3 (RR: 2.34) n=14 Cluster 5 (RR: 3.03) n=1 Cluster 1 (RR: 0.08) n=6 Cluster 4 (RR: 0.00) n=8 Cluster 6 (RR: 0.00) n=1 Cluster 7 (RR: 0.00) n=1 mxtMO (◦C) 27.4 (27.2–27.7) 28.2 (26.8–30.6) 25.6 (NA) 24.6 (21.6–27.4) 24.6 (23.0–26.2) 29.6 (NA) 28.9 (NA) ts (◦C) 5.1 (5.0–5.3) 3.1 (2.7–3.9) 5.4 (NA) 5.3 (3.8–8.4) 5.8 (5.1–6.6) 3.8 (NA) 4.5 (NA) apr (mm) 606.9 (584.7–629.1) 798.7 (595.7–952.8) 526.9 (NA) 590.4 (367.6–1086.1) 649.2 (413.8–1035.9) 612.8 (NA) 563.1 (NA) prs (mm) 50.4 (48.6–52.2) 85.1 (80.7–89.1) 26.2 (NA) 40.0 (27.7–49.2) 57.1 (47.1–70.5) 80.0 (NA) 60.2 (NA) sl (degrees) 3.6 (3.1–4.0) 5.1 (0.3–9.7) 3.4 (NA) 3.6 (0.8–8.7) 8.6 (4.8–11.4) 5.3 (NA) 4.5 (NA) sp (index) 2.4 (2.3–2.6) 1.3 (1.0–3.0) 1.8 (NA) 1.7 (1.0–2.7) 2.4 (1.1–3.0) 1.7 (NA) 1.0 (NA) ndvi (index) 150.5 (144.6–156.4) 163.6 (115.1–196.8) 128.1 (NA) 135.3 (100.3–177.5) 126.2 (105.5–148.2) 153.1 (NA) 141.1 (NA) ndvis (index) 5.1 (4.0–6.2) 8.7 (2.0–20.3) 3.6 (NA) 10.1 (4.1–17.9) 3.9 (1.9–7.4) 14.6 (NA) 17.2 (NA) dfav 0.802 (0.724–0.879) 0.509 (0.185–0.847) 0.631 (NA) 0.447 (0.144–0.830) 0.715 (0.541–0.923) 0.690 (NA) 0.785 (NA) wbfav 0.606 (0.514–0.698) 0.276 (0.131–0.532) 0.850 (NA) 0.745 (0.313–0.955) 0.737 (0.464–0.834) 0.488 (NA) 0.570 (NA) ld (N◦/hectare) 0.530 (0.529–0.531) 0.526 (0.112–1.182) 0.223 (NA) 0.566 (0.405–0.780) 0.413 (0.193–0.923) 1.182 (NA) 0.937 (NA) Note: The values displayed for any predictor include the average value of the set of populations included in the cluster and the range (except for single population clusters). The RR and the number of populations included in each cluster are shown. Cluster numbers correspond to those in Figure 3. Predictors and abbreviations are fully described in Table 1. NA: not estimated. 4DISCUSSION A set of preliminary results suggested that CCHFV could be enzootic in Iberia, or at least in part of the Peninsula (Moraga-Fernández et al., 2020; MSCBS, 2019). Our study confirms circulation of CCHFV in Iberia, where specific antibodies were present in red deer in northern, central and southern Spain as early as 2009 – in keeping with previous findings (Filipe et al., 1985), and where antibody prevalence values were very high locally. These results suggest that CCHFV is enzootic in Iberia. The spatial risk of CCHFV in Iberian red deer populations shapes the distribution range for H. lusitanicum ticks (EstradaPeña et al., 2013b). The seropositive cases detected on the northern shores of Iberia are most probably related to infected Hyalomma ticks transported on hosts (birds, livestock) because the local abundance of other potential vector species is low (Ruiz-Fons et al., 2006). Hyalomma marginatum appears at a low abundance in cattle in the south of the Basque Country (north-central Spain; García-Pérez A, personal communication). This concurs with observations in southern France and Germany (Chitimia-Dolber et al., 2019; Vial et al., 2016). That proves there is a tick trade-off among Iberian regions that is significant in the spread of tick-borne pathogens, and it illustrates the potential establishment of CCHFV in northern Iberia where there are abundant wild and domestic ungulate populations. A recent time-series survey in southern Spain shows that exposure to CCHFV within a specific red deer population is a dynamic process (Unpublished data). Our cross-sectional approach may thus have some temporal biases affecting the findings. Even a pure cross-sectional design, with samples gathered at a single point in time, would have biases because each population may follow its independent dynamics (e.g. Casades-Martí et al., 2020). Even though temporal variations may exist locally, on large spatial scales, these variations should be blurred in heterogeneous spatial distribution scenarios. The modelling output confirms the influential role of the main Hyalomma tick hosts (domestic ruminants and red deer) on the risk of exposure to CCHFV in Iberia. Livestock had a higher-than-expected role in the risk of exposure to CCHFV in the core distribution area of H. lusitanicum ticks, thus suggesting that acaricide treatments may be insufficient, inappropriate or ineffective. The effect of hosts might be associated with the role that ruminants play in maintaining Hyalomma ticks. Ruminants also may help transmit CCHFV, thereby collaborating in maintaining high local virus prevalence in the tick population (e.g. Moraga-Fernández et al., 2020). The Eurasian wild boar also hosts 16 CUADRADO-MATÍAS ET AL. of year monthly average measures) – were included as potential climatic predictors related to water. We selected solar radiation registers (mean and accumulated values) as they may affect soil moisture where CCHFV vectors live. Two variables were selected as potential topographical predictors of the risk of exposure to CCHFV because of their direct influence on vector abundance: altitude and slope. The altitude was obtained at a 100 ×100 m resolution scale from the Land Processes Distributed Active Archive Center (http://LPDAAC.usgs.gov). We estimated the average altitude at the UTM 10 ×10 km square. The slope was calculated, based on altitude, using the Idrisi SLOPE command (Eastman, 2004) and averaged with the spatial unit of the study. We selected habitat variables (vegetation productivity, land cover and soil traits) as potential predictors of CCHFV exposure risk because of their effects on the vectors and the hosts (e.g. Acevedo et al., 2010). We included two predictors related to the Normalized Difference Vegetation Index (NDVI) – average annual NDVI and NDVI seasonality. We considered NDVI registers as proxies for the contribution of the moisture of the soil (Nicholson & Farrar, 1994)to vegetation productivity. Thus, we included these as indices of the hydric stress that ticks suffer in soil. We derived NDVI values from a monthly dataset obtained from the NASA Goddard DAAC website (http//daac.gsfc.nasa.gov/data/avhrr/) at a resolution scale of 1000 m in 18 years (1982–2000). We estimated NDVI seasonality as the coefficient of variation of average monthly measures within a specific UTM 10 ×10 km square. Land cover variables indicating the proportion of coverage by woodlands, scrubland, scarce vegetation or irrigated crops were obtained at a 250 ×250 m resolution from the CORINE Land Use/Land Cover database (EEA, 2000) and transferred to the UTM 10 ×10 km square level. We characterized soil permeability in three classes from a map of the Spanish Geological and Mining Institute (IGME, 1979) of ground-water aquifers. Average soil permeability values were estimated on the UTM 10 ×10 km square level. Hosts are principal determinants of tick population dynamics. Host population traits may also be relevant predictors of the risk of exposure to CCHFV. Wild and domestic ungulates host high burdens of ticks (Castellà et al., 2001; Ruiz-Fons et al., 2013), contributing significantly to tick abundance (Ruiz-Fons & Gilbert, 2010; Ruiz-Fons et al., 2012). Mammals, including ungulates, additionally replicate CCHFV for a short time after infection and transmit it to feeding ticks (Spengler et al., 2016b). Mammals allow CCHFV co-feeding transmissionbetweenticks,itshorizontaltick-to-tickvenerealtransmissionand its vertical transmission which enable adult females to produce abundant infected offspring. Our aim was not to unravel the many traits of the community of hosts (for the ticks and the virus) that modulate the risk of exposure but rather to find significant associations that help us gain an insight into this unknown system. For this purpose, we selected variablesindicating therelativeabundance ofthe major potential hosts for CCHFV ticks in Iberia, that is, domestic ruminants, red deer, roe deer (Capreolus capreolus) and Eurasian wild boar (Sus scrofa). Livestockcensus datawere obtainedfrom the Spanish Ministry forEcological Transition (https://www.miteco.gob.es/es/) on a regional veterinary unit level for 2008. Data from the veterinary units were downscaled to UTM 10 ×10 km squares. Data on wild ungulate abundance on a large scale are scarce in Iberia. Therefore, we chose to estimate species abundance using the values for environmental favourability for red deer, roe deer and Eurasian wild boar calculated in a previously published study to model Culicoides imicola abundance in Spain (Acevedo et al., 2010). Iberian red deer populations are managed for different purposes and management schemes, which may influence exposure to ticks and thus to CCHFV. Public territories under protection and social hunting grounds promote natural management of wild ungulates with minimum intervention; these may include occasional deliveries of animal feed to attract deer to a specific area before hunting or specific habitat management practices, among other minor interventions. On lands devoted to hunting (hunting estates, game reserves), management practices are highly variable; these may range from occasional feeding before hunting to high-wire fencing with (occasional) artificial year-round feeding. Deer farming in Iberia is an extensive activity. Animals live in batches in large, fenced plots with natural vegetation, are provided with food year-round and are managed for health and productivepurposesseveraltimesperyear(González-Barrioetal., 2015a). Specifically, farmed deer are frequently treated with acaricides. The intensity of managing any deer population in the study was classified as low, medium or high. Several host traits may modulate the exposure individuals have to ticks that ultimately determine CCHFV transmission. In the red deer– H. lusitanicum system, we previouslyfoundsound differencesrelated to deer sex and age (Ruiz-Fons et al., 2013), so these two host individual variables were considered as potentially significant predictors of the risk of exposure to CCHFV.