Long‐term shifts in the seasonal abundance of adult Culicoides biting midges and their impact on potential arbovirus outbreaks

Abstract Surveillance of adult Culicoides biting midge flight activity is used as an applied ecological method to guide the management of arbovirus incursions on livestock production in Europe and Australia. To date the impact of changes in the phenology of adult vector activity on arbovirus transmission has not been defined. We investigated this at two sites in the UK, identifying 150,000 Culicoides biting midges taken from 2867 collections over a nearly 40 year timescale. Whilst we recorded no change in seasonal activity at one site, shifts in first adult appearance and last adult appearance increased the seasonal activity period of Culicoides species at the other site by 40 days over the time period. Lengthening of the adult activity season was driven by an increase in abundance of Culicoides and correlated with local increases in temperature and precipitation. This diversity in responses poses significant challenges for predicting future transmission and overwintering risk. Policy implications. Our analysis not only shows a dramatic and consistent increase in the adult active period of Culicoides biting midges, but also that this varies significantly between sites. This suggests broad‐scale analyses alone are insufficient to understand the potential impacts of changes in climate on arbovirus vector populations. Understanding the impact of climate change on adult Culicoides seasonality and transmission of arboviruses requires the context of changes in a range of other local ecological drivers.


| INTRODUC TI ON
Arthropod-borne viruses (arboviruses) include some of the most important emerging and re-emerging pathogens of humans, livestock, and wildlife worldwide (Gould, Pettersson, Higgs, Charrel, & Lamballerie, 2017;Liang, Gao, & Gould, 2015;Weaver & Reisen, 2010). The seasonal incidence and abundance of arthropod vectors capable of transmitting arboviruses, alongside environmental temperatures enabling replication of the arbovirus in these vectors, are key determinants of the timing, intensity, and spread of outbreaks (Lafferty, 2009;Purse, Carpenter, Venter, Bellis, & Mullens, 2015;Rogers & Randolph, 2006). The potential influence of climate on the distribution and dynamics of both arthropod vectors and the arboviruses they transmit, both historically and in prediction under future climate change scenarios, has been a key area of debate (Gould & Higgs, 2009;Kovats, Campbell-Lendrum, McMichael, Woodward, & Cox, 2001;Lafferty, 2009;Tabachnick, 2016).
The unprecedented incursion and establishment of multiple strains and species of arboviruses transmitted by Culicoides biting midges (Diptera: Ceratopogonidae) across Western Europe is a spectacular example of a shift in global pathogen distribution and has been suggested to be driven by changing climate (Elbers, Koenraadt, & Meiswinkel, 2015;Gubbins, Carpenter, Baylis, Wood, & Mellor, 2008;Purse et al., 2005). This hypothesis is underpinned by the fact that Culicoides are among the vector groups most likely to be affected by changes in temperature and precipitation, being smallbodied (<1.5 mm body length) and entirely reliant on the presence of suitable semi-aquatic habitats for larval development (Purse et al., 2005). Competing hypotheses for these incursions, including changes in livestock husbandry, distribution, and abundance, have previously been discounted (Carpenter, Wilson, & Mellor, 2009;MacLachlan & Guthrie, 2010;Purse et al., 2005).
Insect phenology is a strong biological indicator for the impacts of climate change (Cleland, Chuine, Menzel, Mooney, & Schwartz, 2007;Forrest & Miller-Rushing, 2010;Root et al., 2003) and plays a key role in our understanding potential changes in ecosystem processes (Diez et al., 2012;Rafferty, CaraDonna, Burkle, Iler, & Bronstein, 2013). The magnitude and direction of response to shifts in climate varies between species, often depending on ecological traits, but there is a general trend for spring-time events such as adult emergence, to take place earlier in response to shortening winters in temperate regions (Menzel et al., 2006;Thackeray et al., 2016Thackeray et al., , 2010. The seasonal timing of interactions between vectors, pathogens, and hosts may also shift in ways that may promote or hinder pathogen transmission (Altizer et al., 2006).
Surveillance of adult Culicoides flight periods is currently used in Europe and Australia as a means of demarcating a period during which the probability of arbovirus transmission is negligible (Searle et al., 2014). In northern Europe this enables the lifting of movement restrictions on susceptible livestock during winter months during outbreaks, reducing the economic impact of bluetongue virus (BTV) which is a notifiable arbovirus of ruminants (Carpenter et al., 2009).
Anecdotally, the seasonal period in the UK within which adult flight activity occurs has been lengthening for the past two decades (J. Boorman & P.S. Mellor, Personal Communication). This could not only reduce the utility of approaches that rely on reduced Culicoides adult activity, but also increase the probability of the arboviruses they transmit overwintering in the UK and becoming established, with significant consequences for livestock production.
Established in the 1960s, the Rothamsted Insect Survey network is the longest continuous survey of insects in the world (Harrington et al., 2007). The 12.2 m high suction-traps collect a daily sample of aerial invertebrates that have been used to investigate changes in the population dynamics of pest aphid species (Harrington et al., 2007). Records of the first annual flight of pest aphid species show a significant change in the timing of first flight, with an average of 0.6 days per year earlier whilst last flight and abundance remained relatively constant (Bell et al., 2015).
The Rothamsted suction-traps have been previously used to monitor Culicoides activity (Fassotte et al., 2008) and describe spatial and seasonal patterns of Culicoides flight in relation to meteorological conditions over a single season (Sanders et al., 2011). These collections provide a unique opportunity to examine the seasonal dynamics and abundance of adult vectors over the period of time relevant to that at which global climate has been changing. Here we use this resource to describe changes in the seasonality and abundance of Culicoides vectors of BTV in the Palearctic, at two sites over a period spanning nearly 40 years and investigate the role of local climate conditions and assess the likely impact of these changes on the transmission of arboviruses.

| Trap collections, climatic variables, and livestock density data
The Rothamsted suction-trap network sites are situated in arable areas where Culicoides abundance is low and abundance in collections highly variable (Fassotte et al., 2008;Sanders et al., 2011). Two of the 12 UK sites, however, are situated in pastoral-dominated habitat and had sufficient Culicoides abundance and temporal coverage to enable population changes over a period from 1974 to 2012 to be examined. The Preston, Lancashire (53°51'16"N, 2°45'48"W) and Starcross, Devon (50°37'44"N, 3°27'13"W) trap sites were selected in the north-west and south-west of England respectively, where different changes in climatic conditions have been experienced ( Figure S1). The Preston trap is located within a diverse mixture of broad leaved woodland (27%), heathland (21%), and improved grassland (18%) with smaller proportions of sub-urban (14%) and arable (10%) cover, whilst the area around Starcross is largely a mix of arable (66%) and suburban (20%) land cover with small amounts of improved grassland (10%).
Suction-traps at both sites collected large numbers of Culicoides in 2008 (Sanders et al., 2011) and had near-complete sample records from 1974 to 2012 from which Culicoides were identified. Samples were examined for Culicoides from every fourth day of every even year for Preston and every fourth day of every fourth year, with the addition of 1980, 2008, and 2012 for Starcross. The daily samples for both sites from 2008 (Sanders et al., 2011) were also included in the analysis.
Culicoides were counted and identified according to morphological keys (Campbell & Pelham-Clinton, 1960;The Pirbright Institute, 2007) to species or species group level. Female Culicoides of the Avaritia subgenus were identified to the level of the polyphyletic Obsoletus group, described here as Culicoides obsoletus Meigen, Culicoides scoticus Downes and Kettle, Culicoides dewulfi Goetghebuer, and Culicoides chiopterus Meigen. The long-term storage of these samples precluded molecular analysis to species level and separation of C. dewulfi and C. chiopterus by wing-pattern morphology in older samples was considered to be unreliable. Males of the Obsoletus group were identified to species-level and have been used previously as a proxy for the activity of the females of each species (Sanders et al., 2011;Searle et al., 2014). In addition, age grading from pigmentation of female abdomens was also precluded by the age and storage method of specimens.
The number of Culicoides present in large samples of more than 500 individual Culicoides was estimated using a randomised grid sampling method as described previously (Sanders et al., 2011).
Temperature and precipitation data for 1973-2012 were obtained from the UK Climate Projections (UKCP09) gridded observation data-sets. These cover the UK at 5 km × 5 km resolution with the data for each trap site extracted for the grid square in which it is located ( Figure S2a (1976, 1979, 1981, 1988, 1993-1997, 2000, 2003-2004, and 2010) and a linear interpolation or extrapolation was used for non-census years ( Figure S3). Fine scale spatial information (25 m 2 land parcels) on land cover during the latter half of the study period was derived from CEH Land Cover Maps (Appendix S1; Figure S4).  (Sanders et al., 2011)). In the model, the trap catch on day t was assumed to follow a Poisson distribution with mean μ(t) given by, where b 0 scales the abundance and the a n s and b n s are the seasonality parameters (i.e., they control phenology). The model parameters were allowed to vary between years (to allow for variation in seasonality and abundance between years) and the model also included temporal autocorrelation (to allow for dependence between observations) (see Appendix S2 for details). The fitted models were then used to generate replicated simulated phenology and abundance measures with the trends in these replicated measures analysed by linear regression (see Appendix S2 for details). In addition to assessing the robustness of the trends, this approach also allows a more detailed exploration of the underlying trends in the Rothamsted suction-trap data and, in particular, to disentangle the roles of abundance and phenology from one another.

| Statistical methods
The relationships between the five phenology and abundance measures and climate were assessed for 19 summary climate variables and for density of two livestock species (Table S1). Following exploratory data analysis to assess the form of potential relationships, linear regression was used to fit a straight line relationship between each measure and the variable of interest and testing the significance of the slope. As with the annual trends, this was done using the measures computed from the data and using replicated simulated measures to assess the robustness of any relationships.
Subsequently, multiple regression was carried out using a model including all those variables for which there was a significant univariable relationship with a measure. Model selection started from a model including all variables and proceeded by stepwise deletion of non-significant (p > 0.05) terms. Again, this was done using the replicated simulated measures to assess the robustness of any relationships.

| Implications for arbovirus transmission
To explore the possible consequences of changing Culicoides phenology on the transmission of Culicoides-borne viruses we used a simple model for viral replication inside the adult midge vectors based on cumulative thermal time Wilson, Carpenter, Gloster, & Mellor, 2007). In the model a Culicoides midge is able to transmit BTV once the accumulated degree-days since it was infected equal or exceed the level required for completion of the extrinsic incubation period, with threshold temperature for viral replication and virus replication rate above threshold estimated from outbreaks of BTV-8 in Great Britain (Sumner, Orton, Green, Kao, & Gubbins, 2017).
Assuming no infected midges survive from one season to the next (Wilson, Darpel, & Mellor, 2008), we used the model to estimate the earliest date at which newly infected midges would first become infectious (because they have accumulated sufficient thermal time) and so the earliest time at which transmission could occur. than female Culicoides (31.6%) (Figure 1).

| Annual trends in phenology and abundance
The full results of the statistical modelling are presented in the supporting information, including assessment of the fit of the model to data (see Appendix S1; Figures S5-S26). Model checking indicated that the models adequately captured the data in terms of overall fit, total catch and maximum daily catch.
Between 1974 and 2012 at the Preston site, dates of first appearance were earlier, dates of last appearance were later, seasons were longer and abundance (maximum and mean daily catch) greater over time, though these trends were not significant for all 11 Culicoides species/groups (Table 1; Figure 2). By contrast, no significant trends were identified in any of the measures of phenology or abundance for any of the 11 species/groups at Starcross over the same time period (Table S2; Figure S27). These conclusions are robust to uncertainty in the phenology and abundance measures ( Figure 3).
At Preston, the date of first appearance was 0.5-1.5 days per year earlier for total Culicoides, C. obsoletus group females, C.
Furthermore, the date of last appearance was delayed by a similar amount for these species/groups (Table 1). Accordingly, the length of the active season for these species/groups increased by around 1-3 days per year (Table 1). In addition, the date of last appearance was delayed by around two days per year and season length increased by about three days per year for C. punctatus females (Table 1). There were increases in abundance of between 5% and 13% per year for seven (out of 11) species/groups as measured by maximum daily catch and of between 2% and 10% per year for ten (out of 11) species/groups as measured by mean daily catch (Table 1).
In the underlying statistical models used to infer the annual trends in phenology and abundance measures, there were no significant trends in any of the seasonality parameters (i.e., the a n s and b n s) for any of the species at either Preston or Starcross, nor were there significant trends in the abundance parameter (i.e., b 0 ) for any of the species at Starcross ( Figure S28). By contrast, the abundance parameters increased significantly over time for nine (out of 11) species at Preston ( Figure S28) indicating that the earlier dates of first appearance and later dates of last appearance (Table 1) are a result of increased abundance (and, hence, an increased chance of being trapped) rather than to changes in phenology.

| Climate, host, and land use variables
There was a significant increase in annual mean temperature at Preston (b = 0.03; p < 0.001) between 1974 and 2012 ( Figure   S2a), but no significant trends in either annual total precipitation

| Variables driving the trends in phenology and abundance measures
At Preston, higher annual mean temperatures, annual total precipitation, and annual mean soil moisture were associated with later dates of last appearance and longer seasons, as well as with higher maximum daily catches and mean daily catches ( Figure 4). Furthermore, these associations were significant and consistent across a number of Culicoides species/groups (Figure 4).
At Starcross, there were far fewer significant relationships than at Preston. Annual total precipitation was associated with an increase in mean catch and annual mean soil moisture was associated with an increase in season length, maximum catch, or mean catch for a small number (≤4) of species/groups ( Figure S29). For the remaining climate variables (Table S1) Figure S30). The only exception was for C. dewulfi males, where the maximum and mean daily catches were influenced by both annual mean temperature and annual mean soil moisture in combination ( Figure S30).

| Implications for transmission of Culicoidesborne viruses
Using the model for viral replication there was no significant annual trend in the date at which newly infected Culicoides would F I G U R E 2 Annual trends in Culicoides phenology and abundance at Preston, 1974Preston, -2012. Results are presented for five measures: day (from 1 January) of first appearance (d. first app.); day (from 1 January) of last appearance (d. last app.); season length (seas. length; in days); log 10 maximum daily catch (log 10 max. catch); and log 10 mean daily catch (log 10 mean catch). In each figure the observed measures calculated from the trap catches are shown as red circles. The simulated measures generated from the generalised linear models fitted to the RST data are shown as box and whisker plots (median: horizontal black line; interquartile range: grey box; and 2.5th and 97.5th percentiles: whiskers) Accordingly, a later date of last appearance results in a longer transmission season ( Figure 5).   . Results are presented for five measures: day (from 1 January) of first appearance (d. first app.); day (from 1 January) of last appearance (d. last app.); season length (seas. length); log 10 maximum daily catch (log 10 max. catch); and log 10 mean daily catch (log 10 mean catch). Each figure shows the regression coefficients for the trend line for the measure estimated in two ways. First, the black circles and error bars show the estimate and 95% confidence interval, respectively, based on the measures computed directly from the Rothamsted Insect Survey suction trap data. Second, violin plots show the posterior density (shape), median (circle), and interquartile range (line) for the regression coefficient inferred from the posterior predictive distribution for the data. Plots are coloured red where evidence for the trend is robust (median posterior p-value <0.05) and blue where it is not (median posterior p-value >0.05). Where a plot is missing for Starcross, insufficient individuals were caught to be able to define the dates of first and last appearance (through replacement with new designs over time as they became available), or have been conducted in parallel with attempts at control (e.g., long term monitoring of mosquito populations in Singapore (Ooi, Goh, & Gubler, 2006)).

| D ISCUSS I ON
The increase in the length of the adult flight activity season was related to a significant increase over time in Culicoides abundance (driven by temperature, precipitation, and soil moisture), that could not be explained by other changes in the environment and would be likely to increase the probability of arbovirus transmission and overwintering at the site. The study also emphasises heterogeneity in response that will be challenging to capture as part of applied attempts to predict future risks of arbovirus transmission and overwintering. At −0.5-−1.5 days per year, the observed rate of advancement in the first capture of adult Culicoides at Preston are greater than the reported average for phenological events of terrestrial invertebrates of −0.4 days per year (Thackeray et al., 2016(Thackeray et al., , 2010 and similar to that observed in some Lepidoptera and Homoptera (Bell et al., 2015;Roy & Sparks, 2000). ing site and host availability for Culicoides (Purse et al., 2015(Purse et al., , 2005. Although livestock density, in particular that of sheep and cattle, is a key driver of vector Culicoides abundance in providing hosts and associated breeding habitats (Searle et al., 2014), changes in host livestock density were not associated with the changes observed.
Soil type and soil moisture determine Culicoides larval habitat suitability and availability (Mellor, Boorman, & Baylis, 2000;Purse et al., 2015). The sandy soil at Starcross retains lower soil moisture than the clay loam at Preston. Ephemeral larval habitats may dry more quickly and therefore habitat suitability and availability dependent on precipitation rather than temperature may be limiting the Culicoides popu- The increase in abundance of Culicoides at Preston was driven by the impact of warmer temperatures and increased precipitation experienced at the site over time. The RSTs collect Culicoides undergoing active dispersal flight at 12.2 m, with a male bias and lower trapping efficiency than the standard UV light-suction traps (Sanders et al., 2011). In the absence of comparative studies, there is no evidence that these collections are not proportional to the population engaged in host-seeking behaviour closer to the ground.  (Bessell et al., 2013;Napp et al., 2011;Searle et al., 2014).
This study is the first to examine long-term changes in the phenology and abundance of arbovirus vectors, using a unique data set of standardised trapping over 40 years with observed climatic changes.
The study demonstrates trends in abundance and phenology of Culicoides vectors that vary between species and sites are sensitive to local variation in climate change. Rather than use a simple correlation with climate changes we examined other changes in the potential driv-

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section at the end of the article.