Worldwide Niche and Future Potential Distribution of Culicoides imicola, a Major Vector of Bluetongue and African Horse Sickness Viruses

We modelled the ecoclimatic niche of Culicoides imicola, a major arthropod vector of midge-borne viral pathogens affecting ruminants and equids, at fine scale and on a global extent, so as to provide insight into current and future risks of disease epizootics, and increase current knowledge of the species' ecology. Based on the known distribution and ecology of C. imicola, the species' response to monthly climatic conditions was characterised using CLIMEX with 10′ spatial resolution climatic datasets. The species' climatic niche was projected worldwide and under future climatic scenarios. The validated model highlights the role of irrigation in supporting the occurrence of C. imicola in arid regions. In Europe, the modelled potential distribution of C. imicola extended further West than its reported distribution, raising questions regarding ongoing process of colonization and non-climatic habitat factors. The CLIMEX model highlighted similar ecological niches for C. imicola and the Australasian C. brevitarsis raising questions on biogeography and biosecurity. Under the climate change scenarios considered, its' modelled potential distribution could expand northward in the Northern hemisphere, whereas in Africa its range may contract in the future. The biosecurity risks from bluetongue and African horse sickness viruses need to be re-evaluated in regions where the vector's niche is suitable. Under a warmer climate, the risk of vector-borne epizootic pathogens such as bluetongue and African horse sickness viruses are likely to increase as the climate suitability for C. imicola shifts poleward, especially in Western Europe.

In the present study, occurrences published under one of these synonymic names were subsequently named C. imicola.
Geographical coordinates were derived from published coordinates, locality coordinates (654 data), or maps (727 data). In the latter case, maps were imported as digital scans in ArcGIS software and georeferenced. We subsequently up-scaled the data at 10' granularity using CliMond's spatial grid [7], resulting in a geo-database of 879 distinct cells where the insect has been recorded. When several occurrences were recorded in a single 10' cell, one record only was kept ensuring consistent resolution with climate data and removing duplicate records published several times.

Appendix S2: Model development
Parameters of temperature, soil moisture and cold stress indexes ( Figure SI1) were adjusted from initial values to match modelled survival (ie. positive EI) with reported occurrences.
The model was firstly developed only on training dataset, and then tested on verification dataset and finally tested using the validation dataset (Figure 1). At each step the goodness of fit and sensitivity were evaluated following Anderson et al. [81] and Webber et al. [82] (Table SI2).
Training dataset Culicoides imicola is widespread in the training zone (covering the native range) and the model sensitivity was significant given a high prevalence value (Table SI2). To test it, the immature model was then projected on the verification area.

Verification dataset
The projection of the immature model is shown in Figure SI3a with 5 unfitted locations highlighted with red crosses versus 126 locations correctly modelled as suitable (black symbols). All the 19 grid cells from the verification dataset located in equatorial climate (points located in China, Thailand, India and Vietnam; Figure 1) were modelled as suitable. For the 3 unfitted points located in Greece, the survival was limited by cold stress and subsequently we adjusted it by slightly reducing the impact of cold stress (parameter change in the insert of Figure SI1, completed model in Figure SI3b).
The other 2 unfitted occurrences were located in the Arabian Peninsula, areas experiencing an arid climate ( Figure SI3a) but where there was no evidence of irrigation to explain the presence of the insect at these locations whereas the irrigation scenario correctly reclassified 24 outliers from Egypt to the Sultanate of Oman (black crosses in Figure SI3).
The goodness of fit of the completed model is high given a small prevalence value of the verification dataset (Table SI2).

Validation dataset
The strength of our model was confirmed after projection in the validation area where the sensitivity score was also close to 1 given a prevalence value around 0.3 (Table SI2). Only 6 locations remained unfitted (highlighted by red crosses in Figure SI4) out of 561 occurrences in total. One unfitted occurrence was located in Algeria in a 10' grid cell where there was no evidence of irrigation ( Figure   SI4a). Four other unfitted points were located in Spain and Italy at the fringe of the suitable area which could be explained by local dispersion of the insect. The last unfitted occurrence represented the capture of one C. imicola in a single trap in Switzerland [75], probably the result of a rare event of long distance dispersal. The figure SI4b also highlights 10' cells in which only absences of C. imicola were reported [47,48,75,63,51,52,87,76,34,65,66,69,88,89,70,9,21] and illustrates the sampling effort in Europe. The European and North African absence data must be interpreted in the context of an active biological invasion. This means that some of the apparent absences may be suitable for C.
imicola but have not yet trapped this species due to an invasion lag or the existence of an unsuitable non-climatic habitat factor [90]. Without additional information, it is impossible to distinguish such locations from those that are climatically unsuitable. The difference between the modelled potential distribution and observed distribution is discussed in the main manuscript.  Figure SI6: Forecast variation of ecoclimatic suitability between historical climate (current EI) and 2030 or 2070 projected climate (future EI) from A1B or A2 climate scenarios using CSIRO-MK3.0 or Miroc-h climate models (respectively abbreviated as CS and MR). Quantitative variation illustrated climatic suitability increases or decreases. Small and large EI variations represented respectively differences in the range 3-30 and 31-100. The forecast of quantitative variation for 2070/A2/CS future climate is presented in Figure 3b.