Culicoides Midge Abundance across Years: Modeling Inter-Annual Variation for an Avian Feeder and a Candidate Vector of Hemorrhagic Diseases in Farmed Wildlife

(1) Background: Epizootic hemorrhagic disease virus (EHDV) and bluetongue virus (BTV) are orbiviruses that cause hemorrhagic disease (HD) with significant economic and population health impacts on domestic livestock and wildlife. In the United States, white-tailed deer (Odocoileus virginianus) are particularly susceptible to these viruses and are a frequent blood meal host for various species of Culicoides biting midges (Diptera: Ceratopogonidae) that transmit orbiviruses. The species of Culicoides that transmit EHDV and BTV vary between regions, and larval habitats can differ widely between vector species. Understanding how midges are distributed across landscapes can inform HD virus transmission risk on a local scale, allowing for improved animal management plans to avoid suspected high-risk areas or target these areas for insecticide control. (2) Methods: We used occupancy modeling to estimate the abundance of gravid (egg-laden) and parous (most likely to transmit the virus) females of two putative vector species, C. stellifer and C. venustus, and one species, C. haematopotus, that was not considered a putative vector. We developed a universal model to determine habitat preferences, then mapped a predicted weekly midge abundance during the HD transmission seasons in 2015 (July–October) and 2016 (May–October) in Florida. (3) Results: We found differences in habitat preferences and spatial distribution between the parous and gravid states for C. haematopotus and C. stellifer. Gravid midges preferred areas close to water on the border of well and poorly drained soil. They also preferred mixed bottomland hardwood habitats, whereas parous midges appeared less selective of habitat. (4) Conclusions: If C. stellifer is confirmed as an EHDV vector in this region, the distinct spatial and abundance patterns between species and physiological states suggest that the HD risk is non-random across the study area.


Introduction
Epizootic hemorrhagic disease virus (EHDV) and bluetongue virus (BTV) are globally distributed orbiviruses presenting important threats to livestock and wildlife health.Both viruses are known to infect cattle and sheep with outbreaks occurring in both North America and Europe [1,2].Typically, these viruses circulate in stable enzootic/epizootic patterns [3] with subclinical infections or mild-to-moderate febrile disease in domestic herds [1,4], but in North America, EHDV and BTV have been found to consistently cause disease in wild ruminants where outbreaks in white-tailed deer (Odocoileus virginianus) [5][6][7], pronghorn (Antilocapra americana) [8,9], and bighorn sheep (Ovis canadensis) [6,10] have resulted in high morbidity and mortality.
White-tailed deer (WTD), the most abundant wild ruminant in North America, are particularly susceptible to multiple serotypes of EHDV and BTV [3,[11][12][13][14].These viruses present similarly in WTD and together are referred to as hemorrhagic disease (HD) [3].Clinical signs of HD are characterized by high fever, respiratory distress, severe head and neck edema, hemorrhaging in body tissues, and cracked hooves.Acute infections can result in death within 8-36 h without clinical signs [4].It is currently estimated that WTD in risk areas have a 29% infection with a 20% mortality rate [4,15].WTD are also the most common species of farmed native and exotic wildlife [16,17].WTD farming is an important industry in rural areas in the United States like Texas where antibodies against HD were detected in greater than 80% of surveyed farmed deer [18].Given this high level of seroprevalence, WTD susceptibility to HD presents an ever-increasing economic and population health risk across North America and warrants research supporting disease management and prevention methods [1,3,4,7].
EHDV and BTV are vector-borne viruses transmitted by Culicoides (Diptera: Ceratopogonidae) midges with a few confirmed and several suspected competent Culicoides spp.vectors in different regions; virus distributions are linked to the geographic range of vectors (though disease distributions reflect multiple virus serotypes and vectors) [3].In Florida, where HD causes significant mortality in WTD, several Culicoides spp.are suspected of transmitting EHDV and BTV.The larval ecology of Culicoides species varies widely, from wet treeholes to seepages and stream margins [19,20], making it challenging to clearly delineate the fine scale and local distribution of these viruses based on the distribution of putative vectors.
The range of these viruses is expanding, potentially through a mix of additional competent Culicoides species transported from other regions, expanded Culicoides vector range, and/or altered vector ecology due to climate change [21,22].In North America, recent HD outbreaks have been reported in new regions of Canada that did not historically have HD [23,24].At the same time, exotic serotypes of EHDV and BTV have been introduced to native vector species.For example, a previously unknown serotype of EHDV was recently discovered in North America [12,25], and a sub-Saharan serotype of BTV was introduced into an indigenous vector species in northern Europe [26][27][28][29].How exotic serotypes are introduced to new regions remains unclear, although wind currents are one of the suspected mechanisms involved in midge dispersal [30].Combining these unknowns with the expanding distributions of these viruses puts both naïve domestic livestock and wildlife species in previously unaffected areas at risk of more severe and more frequent outbreaks.Current research efforts by ours and various other teams around the world [31][32][33] are working to define the ecology of Culicoides vectors, namely their seasonal abundance, host, and habitat preferences, as identifying these details can help understand the expansion of HD viruses [34].With this knowledge, we may be able to proactively manipulate herd habitat use to avoid specific areas at times of high vector abundance [35,36].
Here, we modeled spatially explicit Culicoides abundance on a local landscape over two years in the Florida panhandle.We evaluated how the site-specific abundance of several Culicoides species at the parous and gravid physiological states are related to different habitat types over repeated sampling events during the transmission seasons in 2015 and 2016 on a well-studied wildlife ranch.Our objectives were (1) to develop a universal model to predict the abundance of Culicoides spp.across years and states, (2) assess whether the physiological states of different species of Culicoides have contrasting habitat preferences, and (3) determine if abundance and habitat preferences are similar across multiple years.Lastly, we used the resulting model to predict and map the abundance of three Culicoides spp.during the HD transmission seasons to test our hypothesis that different species in different states prefer different habitats.

Study Area
This study was conducted on a 180 ha privately owned, high-fenced ranch in Gadsden County, Florida.The property was divided into a 172 ha hunting preserve and an 8 ha deer breeding facility.In the hunting preserve, there were between 130-150 free-ranging WTD and approximately 150 exotic cervids and bovids of 13 different species, resulting in an animal density of approximately 1.48 animals/ha that was managed with 12 supplementary protein feeders that were regularly filled by ranch personnel [37,38].The animals had access to seven food plots across the ranch, although the timing of planting and food availability varied between years, and palatable vegetation was often low during the summer months during the study period.Water features on the ranch included one large pond (2.3 ha), two small ponds (0.35 ha and 0.1 ha), and a permanent stream fed by spring seepages.There were also 10 double-fenced WTD breeding pens within a larger exclusion zone totaling ~9.3 ha [34]; hunting preserve animals were excluded from the entire breeding pen area.The primary habitat type of the ranch was hardwood hammock with upland short-leaf pine species, such as loblolly pine (Pinus taeda), growing throughout the ranch [34].

Entomological Sampling
Complete entomological trapping methods are detailed in McGregor et al. [38] and Dinh et al. [34] Briefly, 20 spatially random trap locations were selected to represent all habitat types present on the ranch.Each site was equipped with a miniature CDC light trap (Model 2836BQ, BioQuip, Rancho Dominguez, CA, USA) and a blacklight-emitting diode array (Model 2790 V390, BioQuip, Rancho Dominguez, CA, USA), and controlled by a timer to operate starting 1 h prior to sunset to 1 h after sunrise.Samples were collected twice weekly between July 2015-October 2015 and May 2016-October 2016.Trapped midges were identified to species according to Blanton and Wirth [39] and categorized as nulliparous (never bloodfed), bloodfed (engorged with host blood), gravid (laden with eggs), or parous (having previously laid at least one batch of eggs).Resulting counts for Culicoides species at each physiological state were grouped by sampling week, and then count data were filtered to only those samples collected during the suspected HD season (May-October).This limited the data for analysis in this study to samples collected from July to October 2015 and May to October 2016.Histograms were created to visualize filtered count data and identify the most abundant species collected during the 2015 and 2016 seasons for comparative analysis across years and species with large enough sample sizes.Data were then filtered to counts of parous (females seeking blood meals after completing a gonotrophic cycle) and gravid (females carrying eggs) midges, as these states are most likely to transmit viruses through subsequent blood meals [40].

Environmental Data
Seven environmental covariates were derived, rasterized, and resampled to a 10 m cell size and clipped to the extent of the ranch (Table 1, Figure 1).We used the week in the HD season to account for time when developing our models and accounted for space by including the standardized latitude-longitude coordinates of the center of each 10 m raster cell containing a trap site.Additional covariates were included to investigate if different Culicoides species preferred different environmental factors at different physiological states.Similar to previous work detailed in Dinh et al. [34], the Euclidian distances from the nearest feeder and the nearest water body were included as proxies for the availability of potential blood meal hosts, such as WTD, and potential oviposition sites, respectively.Habitat type was included as unordered factor levels based on land cover data reported on the Cooperative Land Cover map (version 3.2) from the Florida Fish and Wildlife Conservation Commission (FWC) and the Florida Natural Areas Inventory (FNAI) [41] and was reclassified as defined in Dinh et al. [34,42] as upland pine covering 35.29% of the ranch, mixed hardwood pine covering 6.91% of the ranch, mixed bottomland hardwood covering 43.50% of the ranch, or rural/developed/pasture covering 14.29% of the ranch.A soil survey of the ranch was downloaded from the US Department of Agriculture's Natural Resources Conservation Service (NRCS) Web Soil Survey application [43] and grouped by the Natural Drainage Class designation because Culicoides spp.larvae typically occur in moist and muddy substrates [39,44,45].Soil types with Map Unit Symbols (MUSYMs) 6, 9, 31, and 36 were classified as well-drained, and soil types with MUSYMs 66, 86, and 88 were classified as poorly drained.Lastly, we included weekly utilization distributions (UD) [46,47] from 15 WTD, 1 fallow deer (Dama dama), and 1 Père David's deer (Elaphurus davidianus) that were collared in a previous study on the same ranch during the midge sampling effort [42,48] to represent the probability of animal presence in the study environment regardless of the environmental characteristics or proximity to feeders [34].These species were confirmed as preferred bloodmeals for midges on this ranch [38].For this analysis, collared animals were resampled from 15 or 30 min intervals to 6 hours using the T-LoCoH R package.The kernelUD function from the adehabitatHR R package was then used to estimate a weekly UD for all movements combined.Kernel density estimates require the computation of a bandwidth.Here, we used least-squares cross-validation (LSCV) and calculated an average LSCV for all weeks as the bandwidth to create UDs per week.Kernel density estimates were output to the same 10 m raster surface as the environmental covariates.We assessed the potential correlation between variables to avoid multicollinearity and model overfit by computing Pearson correlation coefficients (r) between numerical variables and ANOVAs between pairs of numerical and categorical variables in R. We did not find any correlation between the numerical variables as none of the Pearson's r values were greater than |0.7| (Table S1), and there was no correlation between the numerical and categorical variables except for UD and soil (Table S2).

Model Construction
Using the workflow developed by Dinh et al. [34], environmental covariate values were extracted for each of the 20 trap sites.Continuous variables were standardized, and repeated count models were developed using the 'unmarked' R package version 1.1.1[49,50].Repeated count models were fitted with covariates as site-level covariates and a negative binomial prior mixing distribution because insects were not randomly distributed in space, and this distribution allows the density of animals to vary spatially [51].Models were constructed based on the N-mixture model presented by Royle [52], which reasons that individual midges are always available to collect and assumes that a lack of collection suggests non-detectability or an apparent absence.This was applicable to our study because midge abundance, and therefore the likelihood of detection, is likely unaffected by repeated sampling events.Within our models, we used static site-level covariates to define abundance because other time-specific survey-level covariates (such as local weather conditions) were unavailable at the study ranch for the fine spatial resolution used here (10 m 2 raster resolution).
First, we determined if counts were temporally and/or spatially dependent through null models that tested all variations of sampling week and trap coordinates as covariates.Null models also served to identify which distribution was most appropriate for the data: Poisson, negative binomial, or zero-inflated Poisson random variable.The most representative null model using the negative binomial distribution then served as the base for the development of 31 alternative models with different combinations of variables (week, latitude, longitude, feeder, water, habitat, soil, UD (Table S1)) for either the parous or gravid state.All 31 models were run for each species by physiological state for each year and were ranked separately by AIC.(Table 2).From the resulting list of best models for each species/state/year combination, the most common best model (GlobalE) was identified as the 'universal model' and used to predict the weekly abundance of biting midges of all three species for each life stage across the study ranch throughout the 2015 and 2016 transmission seasons, with additional predictions using the individual species/state/year best models provided in the supplement.Maps of predicted abundance were created in R and exported as JPEGs to create weekly GIF animations.Lastly, we visualized each

Model Construction
Using the workflow developed by Dinh et al. [34], environmental covariate values were extracted for each of the 20 trap sites.Continuous variables were standardized, and repeated count models were developed using the 'unmarked' R package version 1.1.1[49,50].Repeated count models were fitted with covariates as site-level covariates and a negative binomial prior mixing distribution because insects were not randomly distributed in space, and this distribution allows the density of animals to vary spatially [51].Models were constructed based on the N-mixture model presented by Royle [52], which reasons that individual midges are always available to collect and assumes that a lack of collection suggests non-detectability or an apparent absence.This was applicable to our study because midge abundance, and therefore the likelihood of detection, is likely unaffected by repeated sampling events.Within our models, we used static site-level covariates to define abundance because other time-specific survey-level covariates (such as local weather conditions) were unavailable at the study ranch for the fine spatial resolution used here (10 m 2 raster resolution).
First, we determined if counts were temporally and/or spatially dependent through null models that tested all variations of sampling week and trap coordinates as covariates.Null models also served to identify which distribution was most appropriate for the data: Poisson, negative binomial, or zero-inflated Poisson random variable.The most representative null model using the negative binomial distribution then served as the base for the development of 31 alternative models with different combinations of variables (week, latitude, longitude, feeder, water, habitat, soil, UD (Table S1)) for either the parous or gravid state.All 31 models were run for each species by physiological state for each year and were ranked separately by AIC.(Table 2).From the resulting list of best models for each species/state/year combination, the most common best model (GlobalE) was identified as the 'universal model' and used to predict the weekly abundance of biting midges of all three species for each life stage across the study ranch throughout the 2015 and 2016 transmission seasons, with additional predictions using the individual species/state/year best models provided in the supplement.Maps of predicted abundance were created in R and exported as JPEGs to create weekly GIF animations.Lastly, we visualized each model's goodness of fit by plotting the actual counts of each species at each state against the predicted counts at each trap location over the duration of the sampling periods.

Entomological Sampling
Entomological data collected between July-October 2015 and May-October 2016 are summarized in Figure S1 and Table S4.The most abundant species that were collected in both 2015 and 2016 are C. stellifer, C. haematopotus, and C. venustus (Figure S1).In 2015, a total of 28,887, and in 2016, a total of 29,296 midges were trapped and identified as either C. haematopotus, C. stellifer or C. venustus with counts subdivided by species and physiological status as parous, gravid, bloodfed, or nulliparous (Figure 2A).Here, only the parous and gravid categories are considered, with a relative abundance of each life stage for each species illustrated over time in Figure 2B.

Universal Model for Culicoides Species Abundance
Null models were run with different combinations of week and standardized latitude and longitude variables.For 2016, the best null model includes week and location, whereas the best null model for 2015 only includes location.Although time was not significant for the 2015 data, the null model that included week, latitude, and longitude was used for 31 alternative models to ensure we captured the variation of the HD season.
We then applied the 31 alternative models to evaluate three Culicoides species at the parous and gravid states for two years and encountered differences in the resulting best models for each of the 12 situations.With our priority being to draw comparisons among species, state, and year, we identified our universal model (GlobalE) as the most common,

Universal Model for Culicoides Species Abundance
Null models were run with different combinations of week and standardized latitude and longitude variables.For 2016, the best null model includes week and location, whereas the best null model for 2015 only includes location.Although time was not significant for the 2015 data, the null model that included week, latitude, and longitude was used for 31 alternative models to ensure we captured the variation of the HD season.
We then applied the 31 alternative models to evaluate three Culicoides species at the parous and gravid states for two years and encountered differences in the resulting best models for each of the 12 situations.With our priority being to draw comparisons among species, state, and year, we identified our universal model (GlobalE) as the most common, best model among the 12 best models identified, regardless of the ∆AIC > 2 for some situations (Table 2).
The universal model identified to predict C. haematopotus, C. stellifer, and C. venustus over two years at the parous and gravid states returned coefficient estimates for all covariates included in this study, as summarized in Table 3.Like predictions from running null models, all situations (C.haematopotus parous, C. haematopotus gravid, etc.) were significantly affected by week in 2016, with abundance decreasing throughout the season, while time had little effect on abundance in 2015, with only C. haematopotus gravid samples significantly decreasing over the season.When considering habitat, C. haematopotus gravid in 2016, C. stellifer gravid in 2015 and 2016, and C. venustus gravid in 2015 and 2016 have significant positive coefficients for mixed bottomland hardwood habitats.Across C. haematopotus and C. stellifer for both years, parous midge abundance appears to be marginally affected by mixed bottomland hardwood habitats, as only the C. stellifer parous samples in 2015 were significantly reduced in this habitat.In both years, across species or states, strong patterns were observed for distance to supplementary feeders and permanent water sources.Negative coefficients for both variables imply higher midge abundance with closer proximity to these features, particularly for gravid midges.Proximity to water may suggest higher midge abundance; however, there is also a positive correlation for well-drained soil for both parous and gravid C. haematopotus in 2015 and 2016, both parous and gravid C. stellifer in only 2016, and only gravid C. venustus in both 2015 and 2016.Lastly, only C. haematopotus parous in 2016 and, in 2015, C. stellifer and C. venustus gravid had significant correlations relating abundance to deer utilization of the area (UD).
Species-specific differences in covariates are more ambiguous, likely because of sample size differences among species.The most prominent species difference is with C. haematopotus for which both parous and gravid abundance for 2015 (β = 1.2017 and 0.4277, respectively) and 2016 (β = 1.0429 and 0.2375, respectively) is significantly increased with higher longitude, while C. stellifer parous abundance decreases with longitude for both years (β = −0.3199and −0.4059, respectively).Culicoides stellifer and C. venustus gravid abundances increase in both 2015 (β = 1.3422 and 1.4783) and 2016 (β = 2.5083 and 2.3591) in mixed bottomland hardwood habitats.Culicoides haematopotus, as well as C. stellifer gravid abundance, are also significantly reduced in rural/developed/pasture habitats in both years.In general, all three species in both states have increased abundance closer to supplementary feeders, closer to permanent bodies of water and in habitats with well-drained soil, while abundances did not seem to be strongly affected by UD (Table 3).
The universal model was considered the best model (∆AIC = 0) for 2015 C. haematopotus gravid, 2016 C. haematopotus parous and gravid, and 2016 C. stellifer gravid; however, it factored in all covariates (week, latitude, longitude, habitat, feeder, water, soil, and UD) when predicting Culicoides midge abundance.This resulted in higher ∆AICs for the universal model in situations with limited sample sizes, such as parous C. venustus samples in 2015, which had an ∆AIC of 8.3 (Table 2).For this reason, we developed secondary best models for the remaining combinations of year, species, and states for which the universal model was not the most parsimonious (Table S5).Furthermore, we opted to exclude all predictions for 2015 C. venustus parous from our analysis because the sample size was too small to generate representative models.

Spatial Predictions Using the Universal Model
For simplicity, we have displayed the spatial predictions for the 14th week of the 2015 and 2016 hemorrhagic disease seasons using the universal model to predict C. haematopotus, C. stellifer, and C. venustus abundance at the parous (Figure 3) and gravid (Figure 4) states (GIFs for the entire season are animated in Figure S2 for 2015 and Figure S3 for 2016).The 14th week was approximately representative of the average midge counts for each species for the entire season.The clearest spatial difference observed in predictions with the universal model is that parous C. haematopotus would be most abundant on the northwestern section of the ranch close to traps 5 and 15 (Figure 3A), particularly close to the largest permanent body of water, while parous C. stellifer abundance is consistently highest on the southeastern side of the ranch near trap 17 distributed along the creek boundary (Figure 3B).Similar spatial patterns are observed in the predicted abundance for C. haematopotus and C. stellifer in their gravid states, but both species are predicted to have a wider geographic range in these general areas than in their parous states (Figure 4).These differences in spatial patterns are similar in both years, with the range for C. haematopotus concentrated in the northwest corner of the ranch, while C. stellifer occupied the length of both creeks running approximately north to south.Parous midges of both C. haematopotus and C. stellifer were expected to populate areas with primarily upland pine habitats, with the greatest abundance predicted along the border of upland pine and mixed bottomland hardwood habitats.These bordering areas approximately overlap with the border between poorly drained soil and well-drained soil.In contrast, gravid midges of these two species demonstrate a strong preference for mixed bottomland hardwood habitats with poorly drained soil, expanding out into the mixed hardwood pine habitats, which were avoided by the parous states.Other parts of the ranch with rural/developed/pasture habitats were predicted to be most populated by C. stellifer parous midges, with these areas being avoided by C. stellifer gravid females, C. haematopotus parous females, and C. haematopotus gravid females.
The predicted abundance for parous C. haematopotus is comparatively low throughout the season and gradually decreases until week 26 in both years (Figure 5).In contrast, the abundance of parous C. stellifer is higher than C. haematopotus, and was predicted to slightly increase over the transmission season in 2015 but decrease over the transmission season in 2016 (Figure 5).Abundance of gravid individuals for both species is predicted to be higher and more variable than the abundance of parous individuals, but also gradually decreased throughout the transmission season (Figure 6).
The predicted abundance for parous C. haematopotus is comparatively low throughout the season and gradually decreases until week 26 in both years (Figure 5).In contrast, the abundance of parous C. stellifer is higher than C. haematopotus, and was predicted to slightly increase over the transmission season in 2015 but decrease over the transmission season in 2016 (Figure 5).Abundance of gravid individuals for both species is predicted to be higher and more variable than the abundance of parous individuals, but also gradually decreased throughout the transmission season (Figure 6).
Only 24 samples of parous C. venustus were collected in 2015, so spatial predictions for that year and states were not reliable and were excluded.We include predictions for parous C. venustus in 2016, when there were 96 individuals.While still a limited sample size compared to C. haematopotus and C. stellifer in 2016, which had 2403 and 7592 counts, respectively, the predicted abundance of parous C. venustus in 2016 appeared to be concentrated along the creek throughout the ranch and quickly decreased by week 14 of the transmission season to have a maximum abundance less than 200 midges throughout the ranch (Figure 5).An abundance of both stages of C. venustus is concentrated in the mixed bottomland hardwood habitats with poorly drained soil with a higher predicted abundance for gravid than parous, which also has an overall downward trend over the transmission season (Figure 6).

Discussion
Hemorrhagic disease, caused by EHDV and BTV, is one of the biggest challenges facing deer ranches across North America.Developing a better understanding of the ecology of the Culicoides vectors responsible for transmitting these viruses can help mitigate the threats they pose to deer ranch population management and economic success.The present analysis utilized occupancy modeling to predict the parous and gravid abundance of C. haematopotus, C. stellifer, and C. venustus on a wildlife ranch in the Florida panhandle for the duration of the HD transmission seasons in 2015 and 2016.By evaluating specific ecological criteria of Culicoides vectors, our modeling approach developed a universal model to predict spatially explicit counts of suspected vector abundance by species and physiological status for the most abundant species.The universal model allows for comparing species' abundance and distribution using the same covariates and model structure.These predicted distributions and counts for putative vector species are proxies for HD transmission risk for this landscape.Variation in universal model success led to developing species' specific models which are provided in the supplement.
We observed similar patterns in actual abundance in our collected samples from 2015 and 2016, with peaks in C. haematopotus and C. stellifer parous abundance around weeks 17-21, corresponding to late August-mid-September. Previous estimates suggest this is the high point of the EHDV transmission season [3,53], and this abundance peak agrees with previous work on the same ranch, which also reported a peak in EHDV deer mortalities during September [19].However, the University of Florida Cervidae Health Research Only 24 samples of parous C. venustus were collected in 2015, so spatial predictions for that year and states were not reliable and were excluded.We include predictions for parous C. venustus in 2016, when there were 96 individuals.While still a limited sample size compared to C. haematopotus and C. stellifer in 2016, which had 2403 and 7592 counts, respectively, the predicted abundance of parous C. venustus in 2016 appeared to be concentrated along the creek throughout the ranch and quickly decreased by week 14 of the transmission season to have a maximum abundance less than 200 midges throughout the ranch (Figure 5).An abundance of both stages of C. venustus is concentrated in the mixed bottomland hardwood habitats with poorly drained soil with a higher predicted abundance for gravid than parous, which also has an overall downward trend over the transmission season (Figure 6).

Discussion
Hemorrhagic disease, caused by EHDV and BTV, is one of the biggest challenges facing deer ranches across North America.Developing a better understanding of the ecology of the Culicoides vectors responsible for transmitting these viruses can help mitigate the threats they pose to deer ranch population management and economic success.The present analysis utilized occupancy modeling to predict the parous and gravid abundance of C. haematopotus, C. stellifer, and C. venustus on a wildlife ranch in the Florida panhandle for the duration of the HD transmission seasons in 2015 and 2016.By evaluating specific ecological criteria of Culicoides vectors, our modeling approach developed a universal model to predict spatially explicit counts of suspected vector abundance by species and physiological status for the most abundant species.The universal model allows for comparing species' abundance and distribution using the same covariates and model structure.These predicted distributions and counts for putative vector species are proxies for HD transmission risk for this landscape.Variation in universal model success led to developing species' specific models which are provided in the supplement.
We observed similar patterns in actual abundance in our collected samples from 2015 and 2016, with peaks in C. haematopotus and C. stellifer parous abundance around weeks 17-21, corresponding to late August-mid-September. Previous estimates suggest this is the high point of the EHDV transmission season [3,53], and this abundance peak agrees with previous work on the same ranch, which also reported a peak in EHDV deer mortalities during September [19].However, the University of Florida Cervidae Health Research Initiative has reported that farmed deer mortalities in the state are generally in late September to mid-October [54], which extends the season beyond the sampling of this study.While sampling may not have included abundance peaks related to these mortalities, C. stellifer has previously been linked with EHDV because of its high abundance during outbreaks [55,56] and was corroborated by a study implicating C. stellifer and C. venustus as vectors of EHDV in Florida [19].Our data represent a nearly 2.5-fold higher count of C. stellifer in 2015 than in 2016.Culicoides species are likely to be affected by changes in temperature [21], and on average temperatures in Gadsden County, Florida, were warmer in August 2015 (82.3 • F/27.9 • C) than in August 2016 (81.8 • F/27.7 • C) [57], providing a plausible explanation for the count difference and suggesting that Culicoides abundance and the timing of peaks may be more accurately predicted when detailed climate data are incorporated into a predictive model [58].It may also be that trapping success declines as vector activities shift from night to day as the year progresses [59].
Predicting the timing of peaks in Culicoides abundance would be advantageous for developing a vector management plan (such as targeted insecticide use or habitat modification to reduce midge larvae sites) or, more realistically, an animal management plan to avoid areas and times of high vector abundance (such as moving feeders or changing the time of food availability throughout the 24 h period).As a working wildlife ranch, feeder placement and timing of feeding may both be ways to reduce risk.Although we observed similar peaks in our collected samples around weeks 17-21 in both years, the resulting coefficient estimates from our universal model show time as a significant factor in 2016 but not 2015.Similarly, Culicoides studies in England found that emergence is continual and highly variable throughout the year [60], suggesting that time in the HD season may not be the most reliable factor for predicting abundance [61].Here, we show a more accurate predictor of biting midge abundance may be a combination of covariates such as habitat type, proximity to water, proximity to bloodmeals and soil drainage.For example, gravid midge abundance was significantly increased in mixed bottomland hardwood habitats for C. haematopotus in 2016 and both C. stellifer and C. venustus in 2015 and 2016.Over 75% of the mixed bottomland hardwood habitat on the study ranch was classified as having poorly drained soil and gravid midge abundance was also predicted to be significantly increased closer to water for all gravid situations, except for 2015 C. venustus.Because gravid midges are most likely to occur near potential oviposition sites, like muddy or sandy substrates on the shores of ponds and streams, that provide the required moisture for larval development [45], it can be reasoned that habitat type, soil type, and boundary areas between well-and poorly drained soil types can be a strong predictor of gravid midge abundance [62].
In contrast, our universal model did not demonstrate a specific habitat preference for parous midges, those that have previously taken a blood meal and, therefore, are potentially infectious.This indicates that parous midges of these three species may be equally adapted to live in upland pine, mixed hardwood pine, or mixed bottomland hardwood habitats.Previous work found that regardless of physiological status, adult C. stellifer females preferred habitats in close proximity to supplemental protein feeders [34], yet this preference was not evident when physiological status was factored in, as was also found in the current study.We observed stronger patterns for parous midge habitat use when applying the individual best models for each parous midge cohort investigated.In these models, all parous midge cohorts except for C. haematopotus in 2016 were found in areas that were closer to permanent bodies of water, perhaps because they had recently laid eggs in moist environments along water boundaries.Areas immediately surrounding permanent water bodies would also provide parous midges with new blood meals as WTD and other animals will frequent these habitats throughout the day for hydration [63,64].
Our analysis was focused on parous and gravid states because they are most likely to transmit viruses.Covariate estimates reported here demonstrate differences in used habitats between these life stages, indicating the risk of hemorrhagic disease across the study ranch is not random with clear spatial patterns.Our spatial predictions highlight the striking difference in geography between C. haematopotus and C. stellifer, particularly for the parous life stage.While the covariate estimates do not indicate any species-specific differences, the C. haematopotus and C. stellifer parous maps from 2015 and 2016 reveal that C. haematopotus was most abundant in the northwest corner of the ranch, while C. stellifer was most abundant in the southeast corner.
The highest abundance of C. haematopotus was congregated closest to trap 5, which was on the bank of the largest permanent body of water, and trap 15, which was close to a frequently used supplemental feeder in mixed bottomland hardwood habitat on the border of well and poorly drained soils.Conversely, C. stellifer abundance was greatest surrounding the southern end of the eastern creek on the ranch closest to trap 17, which was on the border of upland pine and rural/developed/pasture habitats and near two supplemental feeders, the nearest of which was in mixed bottomland hardwood habitat.Culicoides haematopotus is known to breed in pond margins and the margins of spring-fed streams with freshwater soil habitats [39,45,65].It is also hypothesized to prefer feeding at tree canopy levels [39,53], and prior bloodmeal analyses have mostly identified birds as host species [38,66], while bloodmeal analyses from C. stellifer frequently reveal cervids, particularly white-tailed deer (Odocoileus virginianus), as hosts [38,67,68].Previous work has also shown that, like C. haematopotus, C. stellifer will occupy the forest canopies for a portion of the day or night, descend for bloodmeals, and then return to the canopy [53].However, niche partitioning among preferred bloodmeal hosts may have induced analogous partitioning between these species.The environmental features associated with areas of abundance for these two species align well with the likely habitats of the preferred hosts, thus verifying the importance of proximity to bloodmeals for both species.
Collected samples for C. venustus were limited for parous midges in both 2015 and 2016, with only 24 and 96 samples, respectively.Count and spatial predictions for 2015 using any model were excluded from this study, and it is likely that the included predictions for 2016 are not reliable.In fact, prior research on the same study ranch found C. venustus at traps in the geographic range of C. haematopotus [38] rather than in the same area as C. stellifer, as shown on the current prediction maps.Culicoides venustus gravid counts were also much lower than counts for C. haematopotus and C. stellifer but were high enough to reasonably model with 331 samples in 2015 and 682 samples in 2016, respectively.Spatial predictions for C. venustus were similar for both years but predicted counts were more variable in 2015 than in 2016, suggesting that there may be a specific sample threshold to meet to generate reliable predictive models.
Small sample sizes, such as those for C. venustus, are a limitation of this study and affected the breadth of analysis possible.These sample sizes may have been a result of the singular trap height used throughout the entire study, as many species of Culicoides frequently move between the ground level and forest canopies throughout the day [53].However, because WTD is bitten near ground level (not in the canopy), the modeled distributions of midges presented here are relevant.A second limitation is the lack of high-resolution local climate and weather data.Specifically, future studies with this design should deploy wind, temperature, and precipitation measuring devices in proximity to traps to better capture local variation, such as wind or thermal refugia, that may increase catches at some sites.Identifying patterns in those factors would better explain relationships between vector abundance, the timing of vector habitat use, as well as host habitat use, which may better predict midge abundance.Lastly, the continuation of this study through to 2017 informed us that the duration of sampling at the start of this study was not conducted long enough [38].Because the peak of reported HD mortalities is not reflected in our 2015 or 2016 sampling, it is possible that abundance peaks do not match our modeled estimates [54].This identifies an area for further study and suggests that year-round trapping may be required to fully determine midge abundance.

Conclusions
Here, we estimated the abundances of three Culicoides species for two physiological states over two years to develop a universal model and better understand the risk of disease transmission during the HD season.Distinct spatial and abundance patterns were identified for both parous and gravid females of two species, C. haematopotus and C. stellifer.Given the putative vector status of C. stellifer for EHDV and lack of evidence for C. haematopotus as a putative vector, the differences in habitat use among these two species suggest that the HD risk to deer is not random across the landscape, rather it is correlated to species-specific biting midge abundance and their habitat use.Continued research to further define these patterns can lead to improved animal management strategies that intentionally concentrate animal use in different areas depending on the predicted HD risk.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v16050766/s1,Table S1.Pearson correlation r values for assessing correlation between continuous covariates used to model Culicoides spp.weekly abundance on a wildlife ranch in the Florida panhandle.Table S2.ANOVA p-values for testing correlation between continuous numerical and categorical variables used to model Culicoides spp.weekly abundance on a wildlife ranch in the Florida panhandle.Table S3.List of models that were run to compare how selected variables relate to Culicoides midge preferences by species, physiological state, and year.Table S4.Summary of midge samples identified to species and physiological state that were collected from a wildlife ranch in the Florida panhandle between July-October 2015 and May-October 2016.Table S5.Individual best models of abundance for three species of Culicoides in the parous and gravid physiological states based on sampling during the hemorrhagic disease (HD) season in 2015 and 2016.Covariate estimates with blank cells were not included in the selected individual best model for the specific year, species, and physiological state, and estimates with shaded cells are identical to those reported in Table S5.Counts for parous C. venustus were too low to use for predictions in 2015.Figure S1.The total number of Culicoides midges collected and identified by species during the 2015 (A) and 2016 (B) sampling seasons on a wildlife ranch in the Florida panhandle.S4 and S1, with covariate estimates displayed in Table S5.S4 and S1, with covariate estimates displayed in Table S5.

Figure 1 .
Figure 1.Map of study ranch in Gadsden County, Florida, illustrating the numbered trap locations, feeder locations, and all environmental variables considered in the current study.Poorly drained soil is identified by the lined area with the rest of the ranch categorized as well-drained soil.Numerals indicate trap site number.

Figure 1 .
Figure 1.Map of study ranch in Gadsden County, Florida, illustrating the numbered trap locations, feeder locations, and all environmental variables considered in the current study.Poorly drained soil is identified by the lined area with the rest of the ranch categorized as well-drained soil.Numerals indicate trap site number.

Figure 2 .
Figure 2. The total number of females of the three most abundant species (Culicoides haematopotus, C. stellifer, C. venustus) collected on a wildlife ranch in the Florida panhandle during 2015 (left) and 2016 (right) HD seasons at the parous, gravid, bloodfed, and nulliparous life stages (A) and the relative abundance of parous (B,C) and gravid (D,E) life stages throughout each season.Here, each week is reported from the first week of sampling in May forward to October.In 2015, sampling started in July (week 10), and in 2016, sampling started in May.

Figure 2 .
Figure 2. The total number of females of the three most abundant species (Culicoides haematopotus, C. stellifer, C. venustus) collected on a wildlife ranch in the Florida panhandle during 2015 (left) and 2016 (right) HD seasons at the parous, gravid, bloodfed, and nulliparous life stages (A) and the relative abundance of parous (B,C) and gravid (D,E) life stages throughout each season.Here, each week is reported from the first week of sampling in May forward to October.In 2015, sampling started in July (week 10), and in 2016, sampling started in May.

Figure 3 .
Figure 3. Predicted parous Culicoides abundance on a wildlife ranch in the Florida panhandle during the 14th week of the HD transmission season using the universal model to make predictions for C. haematopotus (A) and C. stellifer (B) in 2015 and predictions for C. haematopotus (C), C. stellifer (D), and C. venustus (E) in 2016.Also noted are the locations of CDC light traps 5 and 15 (A,C) and 17 (B,D).Counts for C. venustus were too low to use for spatial predictions in 2015.

Figure 4 .
Figure 4. Predicted gravid Culicoides abundance on the study ranch during the 14th week of the HD transmission season using the universal model to make predictions for C. haematopotus (A), C. stellifer (B), and C. venustus (C) in 2015, and predictions for C. haematopotus (D), C. stellifer (E), and C. venustus (F) in 2016.

Figure 3 .
Figure 3. Predicted parous Culicoides abundance on a wildlife ranch in the Florida panhandle during the 14th week of the HD transmission season using the universal model to make predictions for C. haematopotus (A) and C. stellifer (B) in 2015 and predictions for C. haematopotus (C), C. stellifer (D), and C. venustus (E) in 2016.Also noted are the locations of CDC light traps 5 and 15 (A,C) and 17 (B,D).Counts for C. venustus were too low to use for spatial predictions in 2015.

Figure 3 .
Figure 3. Predicted parous Culicoides abundance on a wildlife ranch in the Florida panhandle during the 14th week of the HD transmission season using the universal model to make predictions for C. haematopotus (A) and C. stellifer (B) in 2015 and predictions for C. haematopotus (C), C. stellifer (D), and C. venustus (E) in 2016.Also noted are the locations of CDC light traps 5 and 15 (A,C) and 17 (B,D).Counts for C. venustus were too low to use for spatial predictions in 2015.

Figure 4 .
Figure 4. Predicted gravid Culicoides abundance on the study ranch during the 14th week of the HD transmission season using the universal model to make predictions for C. haematopotus (A), C. stellifer (B), and C. venustus (C) in 2015, and predictions for C. haematopotus (D), C. stellifer (E), and C. venustus (F) in 2016.

Figure 4 .
Figure 4. Predicted gravid Culicoides abundance on the study ranch during the 14th week of the HD transmission season using the universal model to make predictions for C. haematopotus (A), C. stellifer (B), and C. venustus (C) in 2015, and predictions for C. haematopotus (D), C. stellifer (E), and C. venustus (F) in 2016.

Figure 5 .
Figure 5. Predicted parous midge counts according to the universal model compared to the actual midge counts observed during the hemorrhagic disease season on a wildlife ranch in the Florida panhandle presented by week and trap number with C. haematopotus (A) and C. stellifer (B) data from 2015, and C. haematopotus (C), C. stellifer (D), and C. venustus (E) data from 2016.Counts for C. venustus were too low to use for spatial predictions in 2015.

Figure 5 .
Figure 5. Predicted parous midge counts according to the universal model compared to the actual midge counts observed during the hemorrhagic disease season on a wildlife ranch in the Florida panhandle presented by week and trap number with C. haematopotus (A) and C. stellifer (B) data from 2015, and C. haematopotus (C), C. stellifer (D), and C. venustus (E) data from 2016.Counts for C. venustus were too low to use for spatial predictions in 2015.

Figure 6 .
Figure 6.Predicted gravid midge counts according to the universal model compared to the actual midge counts observed during the hemorrhagic disease season on a wildlife ranch in the Florida panhandle presented by week and trap number with C. haematopotus (A), C. stellifer (B), and C. venustus (C) data from 2015, and C. haematopotus (D), C. stellifer (E), and C. venustus (F) data from 2016.

Figure 6 .
Figure 6.Predicted gravid midge counts according to the universal model compared to the actual midge counts observed during the hemorrhagic disease season on a wildlife ranch in the Florida panhandle presented by week and trap number with C. haematopotus (A), C. stellifer (B), and C. venustus (C) data from 2015, and C. haematopotus (D), C. stellifer (E), and C. venustus (F) data from 2016.
Figure S2.Still image of GIFs included in the supplemental PowerPoint file illustrating the predicted parous and gravid Culicoides abundance on a wildlife ranch in the Florida panhandle during the 2015 hemorrhagic disease transmission season using the universal model to make predictions for parous C. haematopotus (A) and C. stellifer (B), and gravid C. haematopotus (C), C. stellifer (D), and C. venustus (E).

Figure S3 .
Still image of GIFs included in the supplemental PowerPoint file illustrating the predicted parous and gravid Culicoides abundance on a wildlife ranch in the Florida panhandle during the 2016 hemorrhagic disease transmission season using the universal model to make predictions for parous C. haematopotus (A), C. stellifer (B), and C. venustus (C), and gravid C. haematopotus (D), C. stellifer (E), and C. venustus (F).
Figure S4.Predicted parous Culicoides abundance on a wildlife ranch in the Florida panhandle during the 14th week (8 March 2015-8 July 2015, 8 January 2016-8 May 2016) of the sampling season using the individual best model for each year, species, and physiological state to make predictions for C. haematopotus (A) and C. stellifer (B) in 2015, and predictions for C. haematopotus (C), C. stellifer (D), and C. venustus (E) in 2016.Counts for C. venustus were too low to use for spatial predictions in 2015.Figure S5.Predicted gravid Culicoides abundance on a deer ranch in the Florida panhandle during the 14th week (8 March 2015-8 July 2015, 8 January 2016-8 May 2016) of the HD transmission season using the individual best model for each year, species and physiological state to make predictions for C. haematopotus (A), C. stellifer (B), and C. venustus (C) in 2015, and predictions for C. haematopotus (D), C. stellifer (E), and C. venustus (F) in 2016.Figure S6.Still image of GIFs included in the supplemental PowerPoint file illustrating the predicted parous and gravid Culicoides abundance on a wildlife ranch in the Florida panhandle during the 2015 hemorrhagic disease transmission season.Maps were created using the individual best model for each year, species, and physiological state to make predictions for parous C. haematopotus (A) and C. stellifer (B), and gravid C. haematopotus (C), C. stellifer (D), and C. venustus (E).Best model IDs are noted in each panel and correspond to Tables Figure S7.Still image of GIFs included in the supplemental PowerPoint file illustrating the predicted parous and gravid Culicoides abundance on a wildlife ranch in the Florida panhandle during the 2016 hemorrhagic disease transmission season.Maps were created using the individual best model for each year, species, and physiological state to make predictions for parous C. haematopotus (A), C. stellifer (B), and C. venustus (C), and gravid C. haematopotus (D), C. stellifer (E), and C. venustus (F).Best model IDs are noted in each panel and correspond to Tables Figure S8.Predicted parous midge counts according to the individual best model compared to the actual midge counts observed during the hemorrhagic disease season on a wildlife ranch in the Florida panhandle presented by week and trap number with C. haematopotus (A) and C. stellifer (B) data from 2015, and C. haematopotus (C), C. stellifer (D), and C. venustus (E) data from 2016.Counts for C. venustus were too low to use for spatial predictions in 2015.Figure S9.Predicted gravid midge counts according to the individual best model compared to the actual midge counts observed during the hemorrhagic disease season on a wildlife ranch in the Florida panhandle presented by week and trap number with C. haematopotus (A), C. stellifer (B), and C. venustus (C) data from 2015, and C. haematopotus (D), C. stellifer (E), and C. venustus (F) data from 2016.

Table 1 .
Covariates used to model weekly Culicoides spp.abundance on a wildlife ranch in the Florida panhandle.

Table 2 .
Corresponding universal model ∆AICs for each year, Culicoides species and state investigated on a wildlife ranch in the Florida panhandle, as well as the individual best model identified for each iteration and the covariates of importance.Gray boxes indicated variables removed from best models.

Table 3 .
Universal models of abundance for three species of Culicoides in the parous and gravid states based on sampling during the hemorrhagic disease season in 2015 and 2016.Counts for parous C. venustus were too low to use for predictions in 2015.