Quantifying the potential for bluetongue virus transmission in Danish cattle farms

We used a mechanistic transmission model to estimate the number of infectious bites (IBs) generated per bluetongue virus (BTV) infected host (cattle) using estimated hourly microclimatic temperatures at 22,004 Danish cattle farms for the period 2000–2016, and Culicoides midge abundance based on 1,453 light-trap collections during 2007–2016. We used a range of published estimates of the duration of the hosts’ infectious period and equations for the relationship between temperature and four key transmission parameters: extrinsic incubation period, daily vector survival rate, daily vector biting rate and host-to-vector transmission rate resulting in 147,456 combinations of daily IBs. More than 82% combinations of the parameter values predicted > 1 IBs per host. The mean IBs (10–90th percentiles) for BTV per infectious host were 59 (0–73) during the transmission period. We estimated a maximum of 14,954 IBs per infectious host at some farms, while a best-case scenario suggested transmission was never possible at some farms. The use of different equations for the vector survival rate and host-to-vector transmission rates resulted in large uncertainty in the predictions. If BTV is introduced in Denmark, local transmission is very likely to occur. Vectors infected as late as mid-September (early autumn) can successfully transmit BTV to a new host until mid-November (late autumn).

The objective of this study was to use a mechanistic transmission model to: i) estimate the potential number of IBs per infectious host resulting from an introduction of BTV at farms in Denmark at any given day by including known parameter estimates and equations together with national vector surveillance and meteorological data; ii) quantify the uncertainty associated with each parameter.

Methods
Estimating microclimatic temperatures in Denmark, 2000-2016. We obtained hourly meteorological parameters (temperature, solar radiation, humidity, wind speed) for the period 2000-2016 from the Danish Meteorological Institute (DMI) at 320 grid points across Denmark. We considered the nearest grid point to each of the 22,004 cattle farms to be the temperature of that farm. The type of land cover within a 500 m radius of each cattle farm was quantified using CORINE Land Cover 37 . The land cover was reclassified as dry meadow (83%), hedges (6%), wet meadow (3%), and forest (3%) 26 . Using a previously published microclimatic model for Denmark 25,26 , these meteorological parameters were converted to estimates of hourly microclimatic temperatures for the period of April 1 st to December 31 st for each of the four different microclimatic habitats. This resulted in four different hourly temperatures for the 22,004 cattle farms over the 17 years. From each of these four microclimatic series we calculated the hourly mean of the first quantiles (>minimum and <25 th percentile), mean of the second quantiles (≥25 th percentile <50 th percentile), mean of the third quantiles (≥50 th percentile <75 th percentile) and mean of the fourth quantiles (≥75 th percentile) over the 17-year period. To quantify the impact of temperature on IBs, we ran our model with all four series of hourly temperatures for each of the four habitats in all combinations with the other model parameter settings (different equations of EIPs, daily survival rate, biting rate, host-to-vector and vector-to-host transmission rate). In addition, to quantify the worst-and best-case scenario, we also identified the maximum and minimum temperature of all farms each hour and for each microclimatic habitat, respectively from this 17-year hourly temperature dataset. An annual average of the six temperature series from each habitat data is presented in Fig. 1. Culicoides density data. Danish

Wet Meadow
Hours (1st April−Dec 31st)  38 , while surveillance collections are only identified to the two ensemble levels. In this study, we assumed the used Ondersteport trap acts like a host and would attract the same number of Culicoides as one host will attract as also assumed by Guis et al. 39 . Thus, we prepare a dataset with counts of biting midges for each week of different years. From the obtained Culicoides abundance data, we calculated the mean of the first quantiles (>minimum and <25 th percentile), mean of the second quantiles (≥25 th percentile <50 th percentile), mean of the third quantiles (≥50 th percentile <75 th percentile) and mean of the fourth quantiles (≥75 th percentile) of the number of Obsoletus and Pulicaris ensembles separately from the daily trap data. We then generated a daily abundance for each of the four estimates for each of the two species ensembles by assuming the abundance would be the same on each day of the week. Finally, we smoothed each of these eight daily Culicoides abundance series by a 15-day running average and used these daily vector abundances on cattle farms in Denmark as estimates of abundance for an average vector season. The resulting estimates of IBs are therefore averages for Denmark ignoring spatial variation as well as year-to-year variation in both microclimatic temperatures and vector abundance. In addition, to quantify the worst and best case scenario, we identified the maximum and minimum number of Obsoletus and Pulicaris ensembles each week from the 10-years Culicoides dataset and smoothed them with 15-day running averages. We present summaries of vector abundance data in Fig. 2. The equations used. The identified equations used in the transmission model and their references are listed in Table 1   www.nature.com/scientificreports www.nature.com/scientificreports/ replication, which is 10.4, 11.0, 14.1 and 13.3 °C respectively. We refer to the three different daily insect survival rate equations as Survival rate equation I 16 , Survival rate equation II 18 and Survival rate equation III 18 . Survival rate equation I shows a very high survival rate (e.g. 90% at 30 °C) whereas survival rate equation III shows a low survival rate (e.g. 6% at 30 °C). All the survival rates are based on the assumption that midge daily survival rate is independent of age. We used an equation for frequency of blood feeding by the biting midges on the host 7 .
Most BTV transmission models use a fixed rate of transmission from host to vector 1,14,23 . In this manuscript, we used two temperature dependent equations described by Turner et al. 14 from the study of Paweska et al. 13 and one fixed rate of transmission (median value from a range of distributions) from host to vector used by Gubbins et al. 1 and Szmaragd et al. 23 (Table 1). We used these three equations for both Obsoletus and Pulicaris ensembles considering that the rate of BTV transmission from host to vector is the same in both species. We refer to the rate of BTV transmission from host to vector as Host-vector transmission I (developed for C. imicola) 13,14 , Host-vector transmission II (developed for C. bolitinos) 13,14 and Host-vector transmission III (fixed rate) ( Table 1). We define the transmission season as the period when a vector ingesting an infected blood meal will survive to become infectious and infect a new host.

Mechanistic model for estimating the number of daily IBs by Culicoides insects.
We used a mechanistic model to estimate the potential number of IBs originating from one infectious host via Culicoides. This is a biological process-driven model based on the above-mentioned parameters. The model has previously been described 40 , and a similar model for another VBD, Setaria tundra (Filarioidea: Onchocercidae), has been described by Haider et al. 41 .
The model was designed to follow daily cohorts of biting midges throughout the season at hourly temperatures estimated for the four habitats: dry meadow, wet meadow, hedges and forest. In the model, biting midges took a blood meal infected with BTV, rested until the gonotrophic cycle was completed and then successfully took a new blood meal. Completion of the EIP was solely dependent on the hourly temperature experienced by the cohort of midges on each consequtive day. After the EIP was complete, we assumed that all subsequent bites by the biting midges were infectious until all vectors in the cohort were dead. We used a rate summation model in which the EIP or blood meal digestion rate was calculated hourly and summed up daily until the virus development/blood meal digestionwas complete.
The steps in the model are described below 41 .
1. The daily survival rates for Culicoides midges were calcuated using the daily mean temperaure recorded/ predicted by the Danish Meteorological Institute and the equations listed in Table 1. We assumed a maximum survival time of 60 d for Culicodes biting midges in Denmark with a daily maximum survival rate of 90% and minimum survival rate of 1% (Table 1) www.nature.com/scientificreports www.nature.com/scientificreports/ based on the successive hourly temperatures for each daily cohort. It was assumed that the biting midges would take a new blood meal immediately after the gonotrophic cycle was complete. 4. The model identified the date of the IBs in each cohort as the date when vectors bite after the EIP was complete. This date was then merged with information on the survival rate of biting midges to calculate how many vectors of the original cohort have survived until that day. 5. The model estimated the proportion of vectors that became infected while taking a blood meal. This was done seperately for the Obsoletus and Pulicaris ensembles by using two temperature-dependent formulae as well as one fixed rate (Table 1). This proportion was multiplied by the number of surviving vectors to estimate the total number of IBs produced by the cohorts. 6. The model assumed 80% of IBs would successfully infect the host by multiplying the total number of IBs by a factor of 0.80 (the rate of BTV transmission from vector to host, Table 1 25 on average. We considered that biting midges would rest randomly in the different habitats surrounding the farm (the larger proportion a specific habitat constitutes around a farm, the larger the proportion of vectors will be resting in that habitat). To adjust for different areas covered by each habitat we multiplied our model output by the IBs estimates (28 times, as there were approx. 28 times more dry meadow than forest or wet meadow) by those estimated from dry meadow temperature. Hedges were two times more abundant than forest or wet meadow and we multiplied the IBs 2 times with those estimated by hedges, and used the complete dataset including the IBs estimated with wet meadow and forest temperature in summary analysis. Finally, our daily estimated IBs resulted in 147,456 IBs.
To identify the worst-and best-case scenarios, we also ran the model with the maximum hourly temperatures and the daily maximum abundance of biting midges and the minimum hourly temperatures and the daily minimum abundance of biting midges, respectively. This gave another 144 daily estimates (4 EIP × 3 Survival rate × 3 rate of transmission from host to vector × 4 habitats) for each of the worst and best-case scenarios for both the Obsoletus ensample and the Pulicaris ensemble. We then summed the IBs estimated for the Obsoletus and Pulicaris ensembles to get the total IB from Culicoides vectors per day. We identified the worst-case scenario as the maximum value of IBs for each day from 144 different estimates. Likewise, we identified the best-case scenario as the daily minimum value of IBs for each day from 144 different estimates. Data analysis. We calculated summary statistics to report the daily mean temperature and 10-90 th percentiles at each of the four microclimatic habitats as well as the DMI weather stations. We also did this for the number of biting midges per day. We calculated the mean IBs and the 10-90 th percentiles from the 147,456 daily IBs. To quantify the impact of the variation in each of the three different input parameters (temperature, vector abundance, habitats etc.), we calculated the mean IBs for each category of that particular parameter while allowing all other parameters to vary in all remaining combinations. To quantify the uncertainty in each of the parameter equations (EIPs, survival rate, host to vector transmission rate), we ran the entire model for each equation one at the time and present the mean IBs for each equation. The mechanistic model for estimating the IBs per host was developed in SAS version 9.4 42 and all summary analyses and plots were performed in R version 3.4.0 43 . Figure 1 summarizes the microclimatic temperature data and data from DMI.  (Fig. 2). We found four or five generations of Obsoletus using the annual average, with the first peak seen in May, the second in June, third in July, fourth in August and the fifth and final peak was surprisingly seen between late October and the beginning of November. Although not as clear as for the Obsoletus ensemble, five or six generations of Pulicaris ensembles were observed between April and November (Fig. S2). Figure 3 summarizes the estimated number of IBs per host with BTV. The mean (10-90 th percentiles) number of IBs per host estimated for BTV in Denmark was 59 (0-73) per infectious host for the period April 1 st to October 31 st , taking into consideration all possible combinations of parameter values. In the worst-case scenario, the daily maximum number of IBs per host was 14,954 IB, per infectious host, while no transmission was detected in the best-case scenario. More than 82% of the estimated values showed more than one IBs per host.

Seasonality. While considering the median IBs of all combinations, the period when vectors could become
infected and subsequently successfully transmit BTV lasted almost 4 months, starting in the first week of May and ending in the third week of August (Fig. 3). According to the median value of all combinations (50 th percentiles), the earliest possible day a cohort of Culicoides could become infected and successfully transmit BTV was in the first week of May with the resulting IBs being delivered up to 80.6 d later and the last day was in the third week of September (or second week of October based on the worst-case scenario; Fig. 3).
We observed three peaks of IBs per host: at the end of May, middle of June, and the final and largest peak was observed between the second and fourth week of July. These peaks of transmission from host to vector correlated with the number of Obsoletus ensemble during the same period. The mean number of IBs with BTV per infectious host was 27 in April, 67 in May, 125 in June, 160 in July, 34 in August, and 3 in September (Fig. S2). We did not find a large discrepancy in the start or end of the season as estimated by the different equations (EIP, survival rate, host-to-vector transmission). However, different temperatures and vector abundances resulted in different durations of the transmission season. While the worst-case scenario showed almost 6 months of transmission, first quantile mean temperature showed around 3 months of transmission (mid-May to mid-August) and first quantile mean Obsoletus abundance showed less than 3 months of transmission (mid-May to the end of August).

Uncertainty of parameter estimates in BTV transmission. The mean IBs per infectious host during
the transmission season (April to October) was estimated as 26 using EIP equation I, 42 using EIP equation II, 41 using EIP equation III, and 77 using EIP equation IV (Fig. 4). Survival rate equation I estimated much higher IBs (mean: 103) compared to survival rate equation II (mean: 23) and survival rate equation III (mean: 13) (Fig. 4).
The mean number of IBs by host-to-vector transmission estimated with equation II was more than ten-fold higher than with equation I (53.5 vs. 4.7 for Obsoletus ensembles and 33.9 vs. 3.1 for Pulicaris ensembles). The  Difference in BTV transmission potential due to natural variation in observed temperature, vector abundance and resting habitats. The number of IBs estimated from the fourth quantile mean temperature was three times higher than when estimated from the third quantile mean temperature, 52 times higher on average compared to estimates from the second quantile temperature, and 128 times higher compared to the first quantile temperature (Fig. 5).
In presence of favorable temperatures (when temperature is high enough to accomplish the EIP in 60 d), variation in vector abundance was the most influential parameter driving the IBs. The number of IBs estimated from the fourth quantile mean Obsoletus ensemble abundance was on average 13 times higher than the IBs estimated from the third quantile mean, 110 times higher than estimates from the second quantile, and 173 times higher than estimates from the first quantile (Fig. 5). For the abundance of the Pulicaris ensemble, the number of IBs estimated from the fourth quantile mean was on average 13 times higher than the IBs estimated from the third quantile mean, 346 times higher than estimates from the second quantile, and 653 times higher than estimates from the first quantile mean (Fig. 5).
The mean number of IBs per host was 49.1 when insects rested at dry meadow temperature, 40.7 at hedge temperature, 20.2 at forest temperature and 13.7 at wet meadow temperature (Fig. 5).

Discussion
Our estimates of BTV IBs per host are based on the assumption that the identified equations used in BTV modelling are equally likely to be correct at estimating the number of IBs. Our model generated 147,456 different IBs by combining different parameter estimates of BTV in Denmark for each day of the potential transmission season (April 1 st to October 31 st ). Instead of calculating point estimates of IBs, we estimated a distribution of IBs per host for each day by combining all the natural variations (temperature, vector abundance and resting habitats types) and uncertainty associated with the different parameter equations (EIP, the daily survival rate, the host-to-vector transmission rate, biting rates and rate of transmission from vector-to-host). This allowed us to predict a wide range of IBs per host each day during the transmission season (April to October). Our estimates showed that only 18% of the combinations had less than one IBs per host during the entire season. The worst-case scenario predicted a very high average of 2,861 IBs for the entire season (April-October), with a maximum daily value of 14,954. The worst-case scenario was identified as the highest number of IBs estimated from the highest temperature recorded every hour over 17 years from any part of the country along with the highest number of biting midges recorded that week. The best-case scenario showed that no transmission would be possible, and  www.nature.com/scientificreports www.nature.com/scientificreports/ was identified as the lowest number of IBs estimated with the lowest hourly temperatures over 17 years from any part of the country, leading to a prediction where virus development was never possible within the lifespan of the vector. However, both best-case and worst-case scenarios are applicable only to the farms with lowest and highest temperature each hour and vector abundance each week and should be considered only as extreme prediction scenarios for the country. The mean number of IBs per host for all combinations was 59 for the entire season, indicating that BTV has a high potential to spread if introduced 44 .
The number of IBs per host refers to the ability of the vector population to transmit a pathogen to the host population from an infectious host 7 and is different from the basic reproduction rate or R 0 in that the number of susceptible animals available is not accounted for. IBs are bites where infectious vira are transmitted to a host. However, a host may not be susceptible to BTV, or may already be infected or may have become immune. How IB translates into R 0 on the given farm depends both on the availability of alternative hosts and on the proportion of vectors that disperse from the farm. The period of viremia in naturally infected animals varies for different strains of BTV with a range of 14-63 d and a mean of 20.6 d 44,45 . Our estimates therefore show that if BTV is introduced into the country, local spread is likely to occur. This finding is supported by empirical observations of BTV outbreaks in 2008 in Denmark, where a large number of animals were infected with BTV and local transmission was observed 6 .
We observed a large difference in estimates of BTV IBs per host due to natural variation in temperature, vector abundance and types of resting habitats. An increase of one quantile mean number of Obsoletus ensemble resulted in IBs that were at least 13 times higher than another quantiles. The 4 th quantiles mean Pulicaris ensemble estimated IBs that were 653 times higher than IBs estimated using the first quantile mean Pulicaris number. This is due to the large variation in Culicoides abundance observed between traps each week in different years or at different locations in Denmark. We found several generations of Obsoletus and Pulicaris ensembles during an average season. The peak of each vector generation coincided with the peak of IBs. In earlier studies, peaks of estimated R 0 based on observed and predicted Obsoletus complex also coincided with a peak in the population of vectors 3 . The Pulicaris ensemble population peaks earlier in the spring, and there was a small peak in the daily IBs in the beginning of May due to vectors in the Pulicaris ensemble. Pulicaris ensemble has previously not been considered to be a competent vector in BTV transmission models 1,34 , although an experimental study found the Pulicaris ensemble to have similar vector competence as the Obsoletus ensemble 30 . The field collected specimens have showed higher prevalence of BTV in Obsoletus ensemble than in Pulicaris ensemble 29 . Pulicaris are abundant in the "non-imicola" region where BTV was detected 46 , and recent studies indicate that they play a role in BTV transmission 29,30 . In studies from southern and central Europe 3,31 , the Pulicaris ensemble populations only constitutes a small fraction of the biting midges populations, yet we found that in Denmark the Pulicaris ensemble  www.nature.com/scientificreports www.nature.com/scientificreports/ populations were around half that of the total population of Obsoletus over the entire period. In general, Pulicaris are considered to be more abundant in Scandinavia than in central and southern Europe 47 and if species in the Pulicaris ensemble are competent vectors, their high abundance makes it necessary to include them in transmission models for Denmark. Our model shows that Pulicaris ensemble could influence the infectious bites per host or R 0 if the vectors are competent for BTV transmission. It is important to note that the host biting rates used to calculate IBs are based on the assumption that the number of female Culicoides collected in Onderstepoort light traps, reflects the biting rate, but light traps may both under-and overestimate actual biting rates 48 .
Temperature plays a critical role in driving many parameters including the EIP, survival rate, biting rate of vectors and host-to-vector transmission rate. In the best-case scenario, we used the minimum temperature among all the farms each hour to estimate the IBs of BTV which predicted transmission was not possible on any day. There was no transmission because virus development was not completed within the maximum lifespan of the vectors. In a situation like this, the number of insects becomes irrelevant; BTV transmission will not be possible despite high abundance of vectors. The model that used the first quantile mean temperature estimated a very low number of IBs per host because virus development was very slow, resulting in only few insects surviving to transmit the infections to new hosts when EIP was completed. One of the survival rates proposed by Gerry and Mullens (2000) 18,49 (Survival rate III) estimated a very low number of IBs, because only a small proportion of vectors survived until the EIP was complete. We estimated higher numbers of IBs per host in dry meadow compared to other habitats, simply because the temperature on average was higher than in any other resting site.
In an earlier study from Denmark, Graesbøll et al. (2012) showed that the temperature and seasonality of vectors determined the period during which an incursion of BTV could lead to epidemic spread 45 . Furthermore, the authors concluded that within the transmission season, the number of affected animals will depend on the temperature and vector abundance, vector behavior and vector ability to locate hosts 45 . If the temperature remains favorable, the size of outbreaks will depend on the vector population. A generation of Culicoides will generally take 4-5 weeks to reach peak abundance. Our findings show that peaks of vector abundance coincide with peak IBs estimates, which indicate that seasonality of vectors drives the IBs.
During the 2008 BTV outbreak in Denmark, most of the cases were identified in late autumn (generally between September and October, while the last case was detected on November 17 th ) 6 . The date of some infections could possibly be earlier due to a delay in outbreak identification, but it has been unclear how BTV is transmitted at the low autumn temperatures found in Denmark. Our model used estimated microclimatic temperatures of potential vector habitats and showed that the last date when a cohort of vectors could be infected and be able to infect new hosts after completion of the EIP was in the second week of September (early autumn). Considering the maximum lifespan of Culicoides is 60 d, our model therefore showed that a successful transmission from vector to host is possible even in mid-November (late autumn) in Denmark.
We found a large variation in the estimated number of IB when using temperature-dependent equations for different parameters. We observed the largest variation in IBs per host between the equations used for the rate of BTV transmission from host to vector. The mean number of IBs estimated from the Host-vector transmission II equation was 10 fold higher IB than for equation I. The mean number of IBs estimated from survival rate equation I was 3.6 times higher than for equation II and 6.3 times higher than for survival rate equation III. The mean number of IBs estimated from EIP equation IV was 3.8 times higher than for EIP equation I, 1.9 times higher than for equation II and 1.8 times higher than for equation III. Such large uncertainty for a parameter will lead to substantial uncertainty in the outputs of a BTV R 0 model, making it difficult to use model predictions to plan prevention and control strategies.
Although host to vector transmission is an important parameter, most models use a fixed number or a range of values for BTV transmission models. We used both temperatures dependent equations and a fixed rate for estimation of BTV transmission. The time it takes from when the virus enters the Culicoides to when it reaches a cell where it can replicate is unknown. However, virus attachment to and entry of a midgut cell within the Culicoides must occur before the peritrophic matrix (an acellular layer) form around the virus if a successful infection is to take place 50 . Therefore, we used the mean temperature of the day when the insect took the blood meal to calculate the temperature dependent probability of successful transmission of virus from host to vector. In total, we identified and used four equations for the EIP, three equations for the daily survival rate, two equations for the host to vector transmission rate and one equation for biting rates. However, there is a need to identify more precise parameters estimates for virus transmission in European Culicoides. Denmark experienced BTV outbreaks in two consecutive years in 2007 and 2008 6 . The IBs per host estimated from our model showed that only 18% of all estimates resulted in IBs per host less than 1. We suggest the quantitative predictions from modelling the transmission of BTV in northern Europe could be improved if more efforts were put into identifying and quantifying the correct relationships between temperature and Culicoides transmission parameters for BTV in European climates.

conclusion
We estimated 147,456 different IBs by combining different parameter estimates of BTV in Denmark for each day of the potential transmission season of which more than 82% of the estimated values showed more than one IBs per host. The mean (and 10-90 th percentiles) number of IBs was 59 (0-73) per infectious host over the transmission period. The best-case scenario in our model showed transmission was not possible. In the worst-case scenario, the transmission season lasted around 6 months (mid-April to mid-October), with a maximum number of IBs per host of 14,954. Therefore, it is likely that local spread will occur if BTV is introduced in Denmark. We identified a large uncertainty associated with the number of IBs estimated by different equations, including those for daily survival rate of Culicoides and host-to-vector transmission rates. We found temperature and the number of vectors to be the most influential factors in determining the BTV transmission (the daily number of IBs). Our model showed that the effective BTV transmission period is long in Denmark, and vectors infected as late as mid-September (early autumn) can successfully transmit BTV to a new host in mid-November (late autumn).