Composition and seasonality of Culicoides in three host environments in Rabat region (Morocco)

The Culicoides Latreille biting midges (Diptera: Ceratopogonidae) are small hematophagous insects, which are biological vectors of viruses responsible for major livestock diseases, such as bluetongue (BT), African horse sickness (AHS), epizootic hemorrhagic disease (EHD) or Schmallenberg disease (SB) (Purse et al., 2015). Several incursions of AHS and EHD have occurred in recent decades in the Mediterranean region, though outbreaks were limited geographically and lasted only a few years. In contrast, after sporadic incursions, the Mediterranean BT epidemiology changed radically in the last two decades. What happened in Morocco illustrates this epidemiological drift in the BT situation compared to other Culicoides-borne viruses. ■ S A N TÉ A N IM A LE E T ÉP ID ÉM IO LO G IE


Summary
Morocco has suffered several outbreaks of Culicoides-borne viruses in recent decades and most studies have focused on Culicoides imicola, considered for a long time as the only important vector. The change in bluetongue (BT) epidemiology in the Mediterranean Basin and Europe over the past two decades has highlighted the role of other Culicoides species in BT virus transmission. The objective of this study was to provide new insights on the Culicoides species composition and seasonality in three different host environments (a horse-riding center, a goat farm and a cattle farm) around Rabat, the capital of Morocco, where BT has been endemic since 2004. Light / suction trap collections were carried out on two consecutive nights at fortnight intervals from May 2016 to May 2017. Culicoides were identified morphologically at the species level when possible. Multivariate analyses were used to compare the impact of the site / vertebrate species, and the collection month on the species communities. In addition, statistical modeling was used to identify environmental drivers of the Culicoides seasonality. A total of 12,460 Culicoides individuals belonging to at least 15 different species were collected during the survey. Culicoides imicola was by far the most abundant species (71.4% of total catches). The site location, and thus the vertebrate species, did not influence the species composition, which was mainly impacted by the month of collection. Surprisingly, the atmospheric pressure was the environmental parameter the most frequently selected in seasonal models. The potential impact of this meteorological parameter along with the other selected variables is discussed. Identifying the environmental parameters driving Culicoides seasonal abundance is the first step to implementing robust Culicoides dynamic models that could later be used in transmission risk modeling.
Morocco experienced AHS outbreaks in the 1960s (serotype 9) and in the late 1980s (serotype 4) (Mellor and Hamblin, 2004), and EHD outbreaks (serotype 6) in 2006 (Savini et al., 2011). During the 1980s AHS outbreak, observed circulation of AHS virus (AHSV) was reported for five years before returning to an absence of transmission (Mellor and Hamblin, 2004). No EHD clinical cases were reported after 2006. Bluetongue (serotype 10, BTV-10) was described for the first time in Morocco in 1956. The virus spread was controlled thanks to the implementation of a vaccination campaign and a temperature drop in the fall (Lhor, 2016). Meanwhile, other BTV serotypes were extensively transmitted in the Mediterranean Basin from 1998, but the disease had apparently remained absent in Morocco until BTV-4 emergence in 2004 (Lhor, 2016). In 2006, BTV-1 was reported in Algeria, then it spread to Morocco in 2006-2007(Lhor, 2016. In 2008, sanitary policies implemented by veterinary services stopped these outbreaks, but the disease reappeared in 2009, with BTV-1 and BTV-4 co-circulating during Culicoides activity periods in both northern and southern parts of the country (Lhor, 2016). The last census reported 305 BT outbreaks across the country with 1456 cases leading to 449 livestock deaths in 2017. The BT endemic situation is weighing on the sheep industry, which increased from 1.6 billion to 2.4 billion euros between 2008 and 2017. The national herd is currently estimated at 20.6 million head, including 2.5 million of the Sardi breed, which are reared for the Eid al-Adha sacrifice fest. The economic impacts of BT include i) direct losses due to mortality, weight loss, abortion, reduction in the fertility rate, and decreased meat production efficiency, and ii) indirect costs due to vaccination, monitoring, surveillance measures, and trade restrictions limiting the access to higher-value markets (Garros and Balenghien, 2017).
A total of 54 Culicoides species have been recorded in Morocco. This list has been mainly established in the 1970s after the 1965-1967 AHS outbreaks and has been comprehensively reviewed and updated recently (Bourquia et al., 2019).
Culicoides distributions were investigated in Morocco in 1994-1995 after the 1989-1991 AHS outbreaks (Baylis et al., 1997), and in 2000-2003 just before the 2004 BT emergence (Lhor, 2016). These studies mainly focused on Culicoides imicola Kieffer, as a proven AHSV and BTV vector, which was, at that time, considered the only significant vector in the Mediterranean region. The rest of the individuals have been identified at group level, mainly Obsoletus and Pulicaris groups (Bourquia et al., 2019). C. imicola abundance (i.e. the annually-averaged mean daily trap catch) has been found to be positively correlated with the Normalized Difference Vegetation Index (NDVI), and negatively with the wind speed . NDVI was suggested to be related to soil moisture and thus to the suitability for C. imicola breeding sites . Wind speed is suggested to affect the abundance of C. imicola by increasing adult mortality rates through desiccation or, more likely according to authors, through the loss of adults due to wind dispersal . The seasonal pattern of C. imicola populations was described in Morocco during the 1994-1995-2010surveys (Baylis et al., 1997Lhor, 2016) but was not analyzed to identify climatic drivers. C. imicola dynamic is described as mainly unimodal, and the catch of C. imicola peaks in general in late summer and fall, with a smaller peak in spring. In Senegal and South Africa, this species is associated with rainfall, NDVI, temperature, wind speed and percentage cover of water bodies (Diarra et al., 2015). Other studies including environmental parameters (related to land use, landscape and climate) conducted in the Mediterranean Basin suggest that these parameters improve the prediction of habitat suitability occupied by C. imicola as well as the accuracy of abundance estimates at local scale (Ippoliti et al., 2013).
The species composition and seasonality of Culicoides in Rabat region were thus investigated to provide new insights on the role of Culicoides in BTV transmission by highlighting the environmental drivers of the seasonality.

Study sites and entomological collections
This study was conducted in the region of the Moroccan capital, Rabat. The city is located on the Atlantic coast in northwestern Morocco, with a Mediterranean climate (characterized by hot, dry summers, and mild, wet winters). The average annual precipitation reaches 492 mm and the average annual temperature is 17.1°C (1970-2000 seasonal averages, WorldClim, http://worldclim.org/version2). The presence of the Atlantic Ocean is conducive to a high relative humidity, which fluctuates depending on the distance from the sea and time of the year.
Culicoides midges were collected around Rabat at two-week intervals from May to October 2016 in two sites, a cattle farm ( 33° 57' 44.1" N; 6° 48' 11.7" W) and a goat farm (33° 51' 36.0" N; 6° 51' 00.0" W), and from May 2016 to April 2017 in a horse-riding center (33° 53' 46.0" N; 6° 49' 53.3" W), using a UV light / suction trap (OVI type) manufactured by the Onderstepoort Veterinary Institute (Venter et al., 2009). The traps were placed near animals and were operated from dusk to dawn for two consecutive nights fortnightly. The Culicoides collected were transferred to 70% ethanol, then, for most individuals, identified morphologically to species-level under a stereomicroscope using the key of Mathieu et al. (2012). Additionally, molecular identification was performed i) to differentiate morphologically close species (Culicoides obsoletus Meigen, Culicoides scoticus Downes and Kettle, or Culicoides montanus Shakirzjanova), ii) to confirm morphological identifications if specimens were damaged, and iii) to confirm new species records (Bourquia et al., 2019). Females of all species were age-graded as nulliparous, parous, gravid or freshly blood-fed after abdominal examination of pigmentation.

Environmental data
Daily rainfall estimates were obtained from 1983 to 2017 from the database of Tropical Applications of Meteorology using SATellite data and ground-based observations (TAMSAT, www.tamsat.org.uk/data/ archive) (Maidment et al., 2014). Day and night land surface temperatures (LSTday and LSTnight), vegetation indices (Enhanced Vegetation Index [EVI] and NDVI) from 2000 to 2017 were extracted from version 6 MODIS MOD11A1 product (Wan et al., 2015) and from version 6 MODIS MOD13A1 (Didan, 2015) for each site location and sample date. Hourly climatic data, including 2-m temperature, 2-m dewpoint temperature, surface pressure, 10-m u-component of wind, and 10-m v-component of wind were extracted for 2016 and 2017 from the fifth generation of the European Centre for Medium-range Weather Forecasts (ECMWF) atmospheric reanalysis (ERA5) of the global climate (Dee et al., 2011). The saturated vapor pressure and the environmental vapor pressure, then the relative and the absolute humidity (see Shaman and Kohn [2009] for calculation details and definitions) were calculated from the 2-m temperature and the 2-m dewpoint temperature. Then, meteorological parameters from ECMWF-ERA5 were rescaled from hourly data to daily data (corresponding to 24 hours which included night collection) by computing the minimum, maximum and mean of each parameter. To determine broadly the diversity per site, species communities were described using a principal component analysis (PCA) on the log-transformed abundances. Then, a between-class analysis was carried out on this PCA to assess the importance of both the site and the collection month. Between-class analysis is a particular case of a PCA, where the variability between groups is optimized.

Statistical analyses
A correlation matrix was produced for all environmental variables to select non-correlated variables to be included in the modeling of Culicoides seasonal abundance. Then, cross correlation maps (CCM) were used to assess the correlation between the Culicoides abundance averaged between sites and the selected environmental variables at different time lags. CCM were used to assess average or accumulated meteorological quantities over a period beginning at a first time lag and ending at a second time lag (Brugger and Rubel, 2013). Analyzing CCM determined which time lags of environmental variables may have had an impact on Culicoides abundance.
Finally, abundance was modeled with these environmental parameters and time lags (after centring and scaling) using a generalized linear model (GLM) or a negative-binomial log linear model to account for any over-dispersion of insect collection data. The model selection was based on corrected Akaike's (cAIC) and Bayesian information criteria (BIC). The validity of the selected GLM was assessed by i) plotting the observed versus the fitted values and testing Pearson's product moment correlation coefficient, ii) graphically checking the normality of the residuals, iii) graphically testing the linearity hypothesis (random distribution of residuals around 0 after they had been plotted by fitted values), and iv) graphically confirming the homogeneity of residuals.
All statistical analyses and graphs were implemented with version 3.4.2 of R software, using ade4 package for PCA and between class analysis, fields for CCM, MuMIn for calculating the corrected Akaike and Bayesian information criteria, and aod for analysis of over-dispersed data.

Culicoides collections
The goat and cattle farms were sampled fortnightly over a period of 22 weeks (12 and 11 collections, respectively, the last sampling was not considered due to electrical failure) during most of the Culicoides activity period (from mid-May 2016 to mid-October 2016). The horse-riding center was sampled fortnightly for 48 weeks (22 collections) to measure the population seasonality over an entire year (from mid-May 2016 to the end of April 2017).
During sampling periods, the temperature conditions were similar to the seasonal normal, with an average annual temperature of 20.5°C in 2016 and 20.6°C in 2017, compared to 20.1°C for the 2000-2015 period ( Figure 1). The 2016 collections began with a relatively dry period with a total of 125 mm rainfall from January to May 2016, compared to 193 mm on average for the same months in the 1986-2015 period (Figure 1). In contrast, the winter of 2016-2017 was wetter than normal with 440 mm rainfall from September 2016 to May 2017 compared to 389 mm on average for the same months in the 1986-2015 period ( Figure 1).
The structure of the PCA on log-transformed abundances was driven by the seven most abundant species (or pair of species), leading to a cumulative projected inertia of 70.4% on the first three axes. The site (or the dominant domestic vertebrate species) explained 17.1% of PCA inertia (p = 0.001, with a permutation test), whereas the collection month explained 26.7% of PCA inertia (p = 0.001). The 'site/host effect' was mainly due to a higher abundance of C. circumscriptus in the cattle farm than in the other two sites (Figure 2A), and, but less significantly, to a higher abundance of C. imicola in the horse-riding center and of C. puncticollis and C. paolae in the goat farm (see   * Data from the weeks when all the sites were sampled so as to compute a global abundance rank for each species; Indiv.: individuals positive and obtained with the wind speed 36 days before collection, which seemed incidental and without any biological sense (the loss of adults is generally due to wind dispersal the day of collection). Thus, only the wind speed on the day of collection was included as a variable for modeling procedure. Finally, the EVI of the week of collection was negatively correlated with C. imicola abundances (Figure 3). Table II shows  The selected GLM included atmospheric pressures (27-to-43-day average before collections) and mean temperatures (25-to-27-day average before collections). It correctly predicted the seasonal pattern of C. imicola populations (R² = 0.778, Table III), i.e. a slow increase in population from March to September when it reached a maximum, then a rapid decrease from October to February (Figure 4).
The selected models (Table III) correctly predicted the seasonal pattern of C. circumscriptus (R² = 0.752), C. cataneii / C. gejgelensis (R² = 0.754) and C. puncticollis (R² = 0.670), with the lowest accuracy for the latter (see Suppl. Mat. X for details on model selection, validation and prediction for all species). The selected models did not successfully predict the bimodal seasonal patterns of C. newsteadi (R² = 0.463), and C. obsoletus / C. scoticus (R² = 0.700). GLM was able to predict the abundance peaks of C. newsteadi populations in fall 2016 and spring 2017, but not in spring 2016 (Suppl. Mat. X). Similarly, GLM was able to predict abundance peaks of C. obsoletus / C. scoticus populations in spring 2016 and 2017, but not fully in fall 2016 (Suppl. Mat. X). Finally, the selected model predicted the general seasonal pattern of C. kingi populations, but failed to reproduce the November peak of abundance (Suppl. Mat. X). The atmospheric pressure was selected as a significant predictor in 5 of 7 models, always with a negative correlation, rainfall in 4 models but only once with a significant positive effect, humidity in 2 models, wind speed in 2 models with a positive correlation, and EVI and temperature only in 1 each (Table III).
Supplementary Material I for abundance per species and site). The between class analysis on the collection month highlighted the higher abundances of C. imicola and C. kingi in August/September, of C. obsoletus / C. scoticus and C. newsteadi in May/June and October, and the lesser abundance of C. puncticollis in June/July ( Figure 2B).
C. imicola populations increased progressively from spring to peak in September (Suppl. Mat. II), whereas C. kingi was abundant only in August/September. Population densities of C. circumscriptus, C. cataneii / C. gejgelensis and C. puncticollis were unimodal with a maximal abundance in August for the first species and in July for the latter two (Suppl. Mat. II). Finally, both C. obsoletus / C. scoticus and C. newsteadi populations showed two peaks of abundance, the first in June, the second in October (Suppl. Mat. II).

Influence of environmental parameters on Culicoides abundance
Supplementary Material III details the correlation matrices produced with all environmental parameters, including EVI and NDVI MODIS parameters, TAMSAT daily rainfall and ECMWF-ERA5 meteorological parameters. Comparison of the day and night land surface temperatures from MODIS products and the daily 2-m temperature from ECMWF-ERA5 showed high correlations (Suppl. Mat. III). Therefore, only the temperature ECMWF-ERA5 data was used for the rest of the analysis. The temperature, which was highly correlated with the environmental vapor pressure, the relative and the absolute humidity, the atmospheric pressure, the wind speed, the rainfall and EVI, which was highly correlated with NDVI, were also retained.
Bivariate correlations between C. imicola abundance and environmental variables highlighted a positive impact of the temperature and absolute humidity for a large range of time lags (Figure 3). The time lags that had the highest correlation coefficients with abundance were selected (Pearson's product moment correlation coefficient ρ = 0.80 for 25-27 days as time lags for temperature and ρ = 0.75 for 23-26 days for absolute humidity). In contrast, the atmospheric pressure was negatively correlated with C. imicola abundance (Figure 3, ρ = 0.83 for 27-43 days). CCM did not indicate any specific relation between the wind speed and the number of C. imicola collected. The best correlation was   are known to breed in specific site types. C. imicola larvae develop in moist clay mud rich in nutrients exposed to sunlight or in moist or water saturated soils rich in organic matter (Braverman et al., 1974). Species such as C. newsteadi and C. puncticollis are found near substrates rich in water-saturated organic matter (Braverman et al., 1974). C. kingi larvae are more likely to be found in sunny and very salty mud (Cornet and Brunhes, 1994), whereas C. cataneii and C. gejgelensis larvae may be found in rivers or pond edges in wet meadows (Garros and Balenghien, 2017). C. obsoletus and C. scoticus may be considered ubiquitous as they develop in forest litter, tree holes, corn silage residues or composting manure (Garros and Balenghien, 2017). Finally, C. paolae larvae are considered specific to decaying prickly pear trees, which are common in the Mediterranean Basin.
The seasonal pattern of C. imicola populations observed around Rabat is consistent with previous studies carried out in Morocco, with a peak in late summer / early fall (Baylis et al., 1997), which is associated with BTV transmission in Northwestern Morocco (Lhor., 2016). This seasonal pattern is usually observed in the Mediterranean Basin (Garros and Balenghien, 2017). The peak of abundance of C. imicola populations is observed at the same period (September/October) in Senegal with a tropical climate, where populations are highest during the rainy season (Diarra et al., 2014). Seasonal patterns of C. obsoletus and C. scoticus vary widely depending on the climate (Garros and

■ DISCUSSION
Culicoides collections carried out in three different host environments around Rabat confirmed the presence of at least 15 Culicoides species, among which were several proven or probable vectors of arboviruses of veterinary interest. C. imicola is a proven BTV and AHSV vector and a suspected EHDV vector (Purse et al., 2015), C. obsoletus and C. scoticus are probable BTV vectors (Purse et al., 2015), whereas C. newsteadi and C. paolae are suspected BTV vectors (Foxi et al., 2016;Foxi et al., 2019).
C. imicola was by far the dominant species in the three sites combined. It has been reported as the most abundant and frequent species in Morocco (Baylis et al., 1997;Lhor, 2016). The species composition was similar between sites, the month of collection having more impact than the site / host species. Even if the site and host effects were intertwined, it is likely that species composition is driven first by the site location because of the location of breeding sites. For instance, Larska et al. (2017) highlight no difference in species composition between cattle and horse farms. The main difference in species composition was due to higher abundance of C. circumscriptus, known as a coastal species (Garros and Balenghien, 2017), in the cattle farm, which was the closest to the coast. Other abundant species  impact on abundance by reducing flight activity (Sanders et al., 2011;Scolamacchia et al., 2014).  add that wind speed negatively affects the abundance of C. imicola in distribution modeling, through the loss of adults caused by wind dispersal. The positive correlation shown in our study should thus be considered as incidental.
It is worth noting that including wind speed on the day of collection in the model did not change the outcome (data not shown). Finally, the temperature, absolute humidity and EVI were rarely selected in the models. This result contrasts with those from other studies where temperature is considered a main positive driver of Culicoides seasonal abundance (Sanders et al., 2011;Rigot et al., 2012;Brugger and Rubel, 2013;Scolamacchia et al., 2014;Diarra et al., 2015). The positive impact of humidity has been less often highlighted in other studies (Diarra et al., 2015) though it is known to impact adult survival (Purse et al., 2015). NDVI, which is correlated to EVI, is associated with high abundance of C. imicola in distribution models (Baylis and Rawlings, 1998;Tatem et al., 2003;Acevedo et al., 2010), but rarely to abundance in seasonality modeling (Diarra et al., 2015). In this latter study, NDVI was higher during the rainy season, and the association with Culicoides abundances may be coincidental.

■ CONCLUSION
The primary objective of this study was to provide new insights on the Culicoides species composition and seasonality around Rabat, to understand better the role of Culicoides in BTV transmission in Morocco, where BT is currently endemic. Neither the site location (except for the cattle farm where the coastal species C. circumscriptus was more abundant than the other ones), nor the main vertebrate species influenced the global species composition. The seasonal pattern of Culicoides described was typical of the Mediterranean climate. Finally, the impact of environmental parameters which may drive Culicoides abundance was investigated, questioning the potential role of the atmospheric pressure.
Balenghien, 2017). In the Mediterranean Basin, as observed in Spain or Sardinia, they exhibit a first main peak of abundance in May/June, and later a secondary peak in October (Lucientes and Alarcón-Elbal, 2016;Foxi et al., 2016). Our observations agreed with these previous reports. The seasonal pattern of C. newsteadi populations highlighted the same bimodal distribution. This species is abundant during spring in Morocco and Sardinia (Lhor, 2016;Foxi et al., 2016). Finally, C. kingi populations were mainly abundant around Rabat in November. On the contrary, in tropical regions of Sudan and Senegal, this species is most abundant in July/August during the rainy season, and also in April when temperature is highest (El Sinnary et al., 1985;Diarra et al., 2014).
Determining the environmental factors driving these seasonal patterns is useful not only to understand better the variations in abundance, but also to be able to predict abundance as a first step for transmission risk modeling. Other factors are of course involved. Adult population seasonality is also the consequence of the long-term impact of environmental factors on the different steps of the life cycle, e.g. the duration of larval development, the longevity of both adult and immature stages, the size and frequency of egg laying (Purse et al., 2015). The same meteorological parameters, such as temperature, rainfall, or wind, may also have a short-term impact on the Culicoides flight activity leading to important daily variation of the proportion of the Culicoides population which is active and can be collected. Moreover, a single meteorological variable, such as the temperature, may have a mainly positive impact on population dynamics for a given range of values, but mainly negative impacts for another range, leading to non-linear effects on population abundance. These complexities may explain why many models have been developed to predict the presence and distribution of Culicoides species, in particular C. imicola in Europe, using climatic factors (Wittmann et al., 2001), satellite imagery (Tatem et al., 2003) or a combination of both (Baylis and Rawlings, 1998), but only a few have described the influence of meteorological and environmental parameters on Culicoides populations using statistical (Sanders et al., 2011;Rigot et al., 2012;Searle et al., 2013;Brugger and Rubel, 2013;Scolamacchia et al., 2014;Diarra et al., 2015) or mechanistic (White et al., 2017) modeling.
The influence of environmental parameters on the abundance was explored for species collected in 2016-2017 around Rabat. The atmospheric pressure was negatively correlated with the abundance of C. imicola, C. newsteadi, C. circumscriptus, C. kingi and C. obsoletus / C. scoticus at different time scales. Although never assessed in Culicoides seasonality modeling before, the atmospheric pressure has long been known to impact insect populations (Wellington, 1946), especially mating or phototaxis behavior (Pellegrino et al., 2013;Zagvazdina et al., 2015). At this stage, it is not possible to conclude if the atmospheric pressure has a real impact on Culicoides population abundance or if this correlation is incidental (a seasonal pattern similar to Culicoides seasonal pattern) or indirect (through another meteorological parameter). This is the main limitation of statistical modeling. The rainfall was positively correlated to C. newsteadi and C. cataneii / C. gejgelensis abundances with time lags corresponding to two weeks before collections, but negatively to C. circumscriptus and C. kingi with longer time lags (about 40 days before collection). Rainfall may have a direct negative impact on Culicoides activity (Murray, 1991), and long-term effects on Culicoides abundance by increasing the availability of breeding sites or perhaps drowning nymphs (Nevill, 1967). Long-term positive effects have been highlighted in areas with temperate climates (Brugger and Rubel, 2013) and short-term negative effects in both temperate (Sanders et al., 2011) and tropical climates (Diarra et al., 2015). The wind speed was positively correlated with the abundance of C. cataneii / C. gejgelensis and C. obsoletus / C. scoticus with a large time lag (up to 45 days before collection). Wind speed on the day of collection has often been reported as having a negative