Tracking Community Timing: Pattern and Determinants of Seasonality in Culicoides (Diptera: Ceratopogonidae) in Northern Florida

Community dynamics are embedded in hierarchical spatial–temporal scales that connect environmental drivers with species assembly processes. Culicoides species are hematophagous arthropod vectors of orbiviruses that impact wild and domestic ruminants. A better sense of Culicoides dynamics over time is important because sympatric species can lengthen the seasonality of virus transmission. We tested a putative departure from the four seasons calendar in the phenology of Culicoides and the vector subassemblage in the Florida panhandle. Two years of weekly abundance data, temporal scales, persistence and environmental thresholds were analyzed using a tripartite Culicoides β-diversity based modeling approach. Culicoides phenology followed a two-season regime and was explained by stream flow and temperature, but not rainfall. Species richness fit a nested pattern where the species recruitment was maximized during spring months. Midges were active year-round, and two suspected vectors species, Culicoides venustus and Culicoides stellifer, were able to sustain and connect the seasonal modules. Persistence suggests that Orbivirus maintenance does not rely on overwintering and that viruses are maintained year-round, with the seasonal dynamics resembling subtropical Culicoides communities with temporal-overlapping between multivoltine species. Viewing Culicoides-borne orbiviruses as a time-sensitive community-based issue, our results help to recommend when management operations should be delivered.


Introduction
Multihost disease systems are influenced by the community dynamics of hosts, vectors and the environment in which they live to facilitate parasite transmission that can vary over space and time [1][2][3][4][5][6][7]. In particular, spatio-temporal dynamics drive the phenological patterns of multihost parasite transmission, as defined by the space-time duality theory [8][9][10]. Influential temporal scales may include short-term processes such as flowering time, while species distribution patterns may be driven by long-term biogeographic processes such as the timing of the last glacial maxima. When modeling host-vector interactions that can be influenced by these varying spatio-temporal scales, we create a risk of mismatching these phenomena when scales are arbitrarily set [11]. Episystem the Culicoides community and the spatial distribution of the species in the preserve. The minimum distance between the traps was 328.7 m, and the minimum distance that connected all the traps was 605.6 m (Supplementary File 2; Figure S2.5). Between July and November, 10 more traps were placed in the property but not included in this study to ensure a balanced sampling effort between months. Timers in traps were set to begin collections one hour prior to sunset and turn off one hour after sunrise. The Culicoides collected were identified using morphological keys [30] and sorted by physiological status (nulliparous, parous, gravid or blood fed) following published methods [78,79].

Environmental Data
Daily temperature (average, maximum and minimum) and precipitation were obtained from the National Oceanic and Atmospheric Administration climate online datasets [80] (station USW00093805), and total stream discharge and average gage height were obtained from the United States Geographical Survey National Water Information System [81] (station USGS2330000). Lags were determined for 15, 30, 45 and 60 days prior to each trap date and included temperature, precipitation and stream data (environmental submodels).

Culicoides Data
Because day-to-day variability in Culicoides count is not relevant to the recruitment process, we used the day with the maximum trap catch each week in our analyses [47]. Data were consolidated into trap-week by a species abundance matrix (1010 observations: 101 weeks × 10 traps) where only the count of total identifiable females (i.e., irrespective of physiological stage) was used.

Culicoides β-Diversity Temporal Profile
We used a multivariate time series analysis of community composition to assess the assemblage of putative HD vectors across an annual cycle. A β-diversity analysis was implemented to test the departure of the Culicoides community from the four season temperate regimes (winter, spring, summer and fall), how environmental fluctuation constrained the assembling process and what other wider temporal scales (≥12 months) were observed. In addition, the distribution of the species among the space and habitat can be a confounder of the temporal change in community composition [9]. We did not observe independent temporal patterns in species composition among traps (Supplementary File 2). Summaries of species distribution among habitats were not reported here since the scope of this study was not to analyze the pattern of diversity partitioning among habitats.
We employed four steps in the time series modeling framework ( Figure 1). First, we averaged the trap-week abundance matrix (n >10) to a weekly abundance species matrix since spatial structures were not significant. This was the response matrix for the model (Figure 1, Step 1). Second, we derived self-informed temporal submodels from the species composition data and, as a regressor, the sampling temporal frame by asymmetric eigenvectors maps (AEM) (Figure 1, Step 2). Third, we assessed the fit of the temporal submodels operating on the observed species abundance matrix by redundancy analysis (RDA) to better understand how temporal scales influenced the models (Figure 1, Step 3). Fourth, we used hierarchical variation partitioning analysis (HVP) [82,83] to assess how the temporal scales and environmental variables explained species composition (Figure 1, Step 4).

Figure 1.
Workflow to analyze Culicoides temporal dynamic and its drivers: the sequence of steps from collecting Culicoides field abundance data, drawing temporal scales and identifying β-diversity scales and their relatedness with environmental variability. White filled boxes outline the aim for each analysis, and dark gray boxes mention the statistical tools used.
Step 1 shows the temporal dynamic of the observed abundance in a nine-species virtual community; therefore, it was the response matrix informing weekly community composition.
Step 2: the sampling temporal profile was thrown to spectral decomposition by asymmetric eigenvectors maps (AEMs) after accounting for the dependency between weeks (autocorrelation threshold). The resulting AEMs were half splitting in positive and negative correlation structures, and a subset was selected by their suitability to model temporal patterns in the observed abundance matrix. Selected subsets were arranged in groups representing broad to narrow temporal waves, hereafter referred to as temporal submodels.
Step 3: the community composition change was modeled with each temporal submodel to assess the fitting of the temporal scale on the observed species abundance matrix. The assemblage of species contributing on each significant temporal scale was untangled by the correlation between each species abundance and Lc for each RDAaxis.
Step 4: the contribution of the temporal submodels and environmental variables to explain the amount of variation was visualized by Venn diagrams after hierarchical variation partitioning analysis. Since temporal scales were nested, even when there were orthogonal structures, both diagrams show intersection areas between temporal submodels as a hierarchical overlapping pattern. Letters inside diagrams are R 2 adj for each unique and intersection fraction in the variation partitioning table (as were introduced in section: 2.4. Culicoides β-diversity temporal profile; equations for the fractions are detailed in [84]). AEM: asymmetric eigenvectors maps; AEMs: AEM eigenfunctions; FWD: forward selection; Hd: Hellinger's distance; Lc: linear combinator scores from an RDAaxis; MMC: multivariate Mantel correlogram; RDA: redundancy analysis; Rho: Rank Spearman correlation; R 2 adj: adjusted determination coefficient. This scheme was adapted considering [9,83] (Figure 1).
A detailed explanation of the statistical toolbox and its expedience has been published elsewhere [9][10][11][85][86][87] and expanded on Supplementary File 2. Asymmetric eigenvectors maps were chosen Workflow to analyze Culicoides temporal dynamic and its drivers: the sequence of steps from collecting Culicoides field abundance data, drawing temporal scales and identifying β-diversity scales and their relatedness with environmental variability. White filled boxes outline the aim for each analysis, and dark gray boxes mention the statistical tools used.
Step 1 shows the temporal dynamic of the observed abundance in a nine-species virtual community; therefore, it was the response matrix informing weekly community composition.
Step 2: the sampling temporal profile was thrown to spectral decomposition by asymmetric eigenvectors maps (AEMs) after accounting for the dependency between weeks (autocorrelation threshold). The resulting AEMs were half splitting in positive and negative correlation structures, and a subset was selected by their suitability to model temporal patterns in the observed abundance matrix. Selected subsets were arranged in groups representing broad to narrow temporal waves, hereafter referred to as temporal submodels.
Step 3: the community composition change was modeled with each temporal submodel to assess the fitting of the temporal scale on the observed species abundance matrix. The assemblage of species contributing on each significant temporal scale was untangled by the correlation between each species abundance and Lc for each RDA axis .
Step 4: the contribution of the temporal submodels and environmental variables to explain the amount of variation was visualized by Venn diagrams after hierarchical variation partitioning analysis. Since temporal scales were nested, even when there were orthogonal structures, both diagrams show intersection areas between temporal submodels as a hierarchical overlapping pattern. Letters inside diagrams are R 2 adj for each unique and intersection fraction in the variation partitioning table (as were introduced in Section 2.4: Culicoides β-diversity temporal profile; equations for the fractions are detailed in [84]). AEM: asymmetric eigenvectors maps; AEMs: AEM eigenfunctions; FWD: forward selection; Hd: Hellinger's distance; Lc: linear combinator scores from an RDA axis ; MMC: multivariate Mantel correlogram; RDA: redundancy analysis; Rho: Rank Spearman correlation; R 2 adj : adjusted determination coefficient. This scheme was adapted considering [9,83] (Figure 1).
A detailed explanation of the statistical toolbox and its expedience has been published elsewhere [9][10][11][85][86][87] and expanded on Supplementary File 2. Asymmetric eigenvectors maps were Viruses 2020, 12, 931 6 of 25 chosen because the temporal dependency is a forward directionally process that can be modeled by monotonic oscillation, and inferences about scales are independent [10,11,85]. The shortest independent time interval for community composition change was computed by multivariate Mantel correlogram (mMC) [87]. Positive (abiotic interactions) and negative (biotic interactions) eigenfunctions (AEMs) were selected separately by forward selection (FW) [88] and grouped representing wide (≥12 months), medium (≥6 months) and narrow temporal scales (1 week to 6 months). The fit of a scale's submodel as a predictor of the Hellinger transformed [89] weekly abundance species matrix was tested by RDA. Variation within scales was judged by the retained profiles (RDA axis ) after FW. With HVP, the relative importance of the environmental dimensions of the niche and neutral theory (birth, death and dispersal rates) drivers in Culicoides β-diversity was assessed as the unique fraction of variation from each environmental submodel (precipitation, temperature and stream level) and each temporal submodel, respectively. The additive fractions were used to judge how environmental constraints or demographic processes in the temporal scales could lead to the observed dynamics in Culicoides composition. The joint fractions explain the intercorrelation between scales and environment. The environmental submodels were selected by the FW using the HVP procedure [84] which allowed collinearity. The significance of the additive and unique fractions was tested by 1000 permutations, but there is no statistical procedure to test the joint fraction [87]. Negative values of R 2 adj lack ecological meaning and were reported as R 2 adj ≈ 0 [87]. Analyses were carried out in R 3.5.1 with vegan [90], adespatial [91] and adapted published codes [10,84,87].

Seasonality in Species Recruitment, Weekly Occurrence and Co-Occurrence Structure
We tested if spring months were a transition period for the coexistence of cold and warm tolerant Culicoides species by analyzing the relationship between diversity structure and species composition. Species recruitment across calendar seasons was analyzed using time to species saturation in Culicoides richness during spring with species-accumulation curves [92] in iNEXT [93]. As such, accumulation curves were fitted with rarefaction and extrapolation sampling curves for asymptotic Hill numbers (q 0 , Chao 2 ) with the species by week incidence matrix. Additionally, we explored if the vector assemblage (C. stellifer, C. venustus, C. debilipalpis, C. biguttatus and C. pallidicornis) had sufficient sample coverage by season. With sample completeness curves [92], the number of weeks necessary to complete 97.5% of the vector species diversity was compared with the total sampled weeks during each season in iNEXT [93]. The linkage between richness and species composition was assessed to explore if the pattern of diversity change was related to loss or replacement of species using nestedness analysis [94] in vegan [90]. A nested pattern meant that species composition was structured within the season, where smaller communities form species subsets within larger communities [95]. Then, observations were sorted by maximizing the richness gradient through time and nestedness calculated using NODF (nestedness metric based on overlap and decreasing fill) [96]. Departure from randomness was tested as the deviance from a quasiswap generated null model (n = 1000).
The hypothesis that vectors are able to support year-round virus transmission was tested with a temporally explicit co-occurrence network. This analysis was selected in virtue of modeling simultaneously the probability of sharing time between Culicoides species at community level but also the case of EHDV-BTV vectors. Two species were considered to co-occur when they were detected the same week. First, those relevant occurrence interactions (OR > 3) in species pairwise association analysis constrained by temporal autocorrelation [97] were retained in an interaction matrix from the weekly species presence data (weeks: n > 0; species: n > 5, incidence > 5%) with sppairs [98]. Modularity (M) summarized the degree of compartmentalization in the species co-occurrence network [99] (i.e., seasons) and was defined by Culicoides species that co-occurred significantly more often with each other than with other species in the community and was estimated according to [100] in igraph [101]. In order to identify the role of each Culicoides species in bridging across seasonal modules, we identified four species role parameters that considered the persistence of a species in a seasonal module (within-module degree, z) and its ability to connect among them (among-module connectivity, c) [102] in bipartite [103]. The cut-off to define species role was fixed as [104]. Culicoides were defined as connectors (z < 1.3; c > 0.6) if they sustained network coherence by interacting among modules. Culicoides species that sustained the coherence of a module were considered a module hub (z > 1.3; c < 0.6). If a midge supported a module and increased network coherence, it was defined as a network hub (z > 1.3; c > 0.6). Peripheral Culicoides (z < 1.39; c < 0.6) co-occurred irregularly with species within a defined module and uncommonly interacted with other modules.

Environmental Thresholds in the Temporal Dynamics of the Vector Assemblage
Environmental thresholds in the temporal partitioning of putative HD vectors C. stellifer, C. venustus, C. debilipalpis, C. biguttatus and C. pallidicornis [61,62] were predicted by the monothetic clustering method multivariate regression tree (MVRT) [105] in mvpart [106]. This analysis allowed us to double-check if seasonality in the vector subassemblage and at community level (Section 2.4) overlapped. Culicoides insignis was not considered due to its scarcity on this property. Abundance data were Chord transformed [87,89] and regressed with the standardized temperature, precipitation and stream variables previously described. Clusters informed epidemiological times, and tree branch configuration was shaped by environmental thresholds. The output of the MVRT showed the sampled weeks organized in clusters that informed about epidemiological times as vectors abundance were regressed. Furthermore, the configuration of the tree relied on how the detected environmental thresholds explained the partitioning of the sampled weeks in clusters. To compute the environmental thresholds, the weeks frame was iteratively split into n clusters, and the value of the environmental variables that minimized the Euclidean distance within a cluster was selected. For the comparison of data between both years, weeks were standardized as epidemiological weeks. The first epidemiological week ended on the first Saturday of January, as long as it fell at least four days into the month. Cross validation (CV) of the tree was performed after 1000 iterations, and tree size was judged following [87]. Indicator species analysis (id) selected species which were statistically representative of each branch node [107].
Tracking variation in the community composition across time revealed that both environmental and temporal constraints emerged to mold the pattern of the Culicoides community, explaining 64% of the observed variation ( Figure 3A). Half of the overall variation accounted for the intercorrelation between neutral and niche drivers (environmental and temporal submodels joint effect; [j]R 2 adj = 0.26; [k]R 2 adj = 0.07). Precipitation did not contribute in the environmental submodel ( Figure 3B) Table S1.1) encompassed the 66% of the community composition (|Rho| < 0.25), but environmental data did not fit the model ( Figure S1. 3), and almost all the species were involved in this profile ( Figure 2D). The compositional changes that occurred in the temporal windows of nine or fewer weeks were autocorrelated ( Figure S1.2), and biotic interactions were not observed (negative AEMs; Table S1.1). In general, Culicoides β-diversity was composed of a series of phenological modules which described sequences of nonoverlapping temporal scales and species composition.   Table S1.1 summarizes temporal scales and their week range, AEMs and statistical significance. The x-axis is shared between panels (A,C). Panel A also shows the rainy (light blue) and dry (light yellow) seasons in Florida. (B,D) Panels indicate compromised species for the changing composition at each temporal profile (RDA) after Spearman rank correlation. Only significant |Rho| showed where size and color followed the magnitude of the species-temporal structure association (Table S1.2). Larval habitat traits are indicated as wet soils and tree holes with shapes. The question mark indicates those species where the larval habitat was not described. Panel (E): Venn diagram represents the amount of variation in the community change composition at each temporal scale by hierarchical variance partitioning. Fractions were described in Culicoides β-diversity temporal profile section and their R 2 adj printed.

Seasonality in Species Recruitment, Weekly Occurrence and Co-Occurrence Structure
Species recruitment (defined as the first appearance of a species) occurred during each seasonal compartment, and the highest recruitment occurred during the transition from the dry into the wet season (Figures 4 and 5). Over the two-year study period, the species inventory was fully represented (Chao2 = 29.49 (29.03; 37.37); sampling coverage (sc = 0.99)). Spring months were key to shaping Culicoides community structure, as interpolation-extrapolation curves show (Figure 4)

Seasonality in Species Recruitment, Weekly Occurrence and Co-Occurrence Structure
Species recruitment (defined as the first appearance of a species) occurred during each seasonal compartment, and the highest recruitment occurred during the transition from the dry into the wet season (Figures 4 and 5). Over the two-year study period, the species inventory was fully represented  (11.25-29.05); sc = 0.94). The rate of increase in the accumulated species across sampling weeks was the steepest during spring months (interpolation curve). During spring, richness was saturated before rarefaction ended. However, in summer months, the sampling effort had to increase 17%, 31% in fall and 32% in winter to reach 98.5% sampling coverage.
The sample completeness curve showed that in each season, 97.5% of the sampling coverage was reached as soon the field sampling effort reached 27 weeks ( Figure S1.7). Consequently, the deployed sampling effort was enough to analyze the seasonal dynamics in the vector assemblage.
Culicoides richness and composition were linked through a pattern of nestedness, with the assemblage organized in a series of nested time-derived subassemblages ( Figure 5A). Assemblages in weeks with low richness were formed with a limited set of species also occurring in high-richness weeks. Nestedness (NODF = 65.45; p < 0.01; matrix fill = 21%; Figure 5A) was more likely to be driven by richness and composition arrangement (NODF row = 69.25; p < 0.01; Figure 5C) than the pattern of species occurrence (NODF column = 53.80; p = 0.65; Figure 5B). Moreover, week order in the nested matrix depicted seasonality because species richness and composition followed a winter-fall-summer-spring gradient (subassemblage to full assemblage).
The Culicoides community was active year-round and connected by the pattern of vector co-occurrence ( Figure 6A). The weekly occurrence network supported the two season modules (M = 0.24; p < 0.01): winter-spring and summer-fall ( Figure 6A). Species displayed most of their interactions with species within the same module ( Figure 6A; lighter shade links), e.g., winter prevalent species and spring prevalent species overlapped their activity during several weeks. In addition, species did not co-occur equally. For example, species occurrence and degree (C d ) were heterogeneously distributed as C. venustus (C d = 15), C. stellifer (C d = 13) and C. haematopotus (C d = 7) and occupied 80% of surveyed weeks in both years, while the remainder were recorded in less than 40% of the weeks ( Figure S1.4, Table S1.4). A total of 90% of species also co-occurred with at least one species out of its module ( Figure 7A; dark shade links). This result emphasizes how connectivity among modules expanded co-occurrence beyond traditional seasons ( Figure 6B). Culicoides venustus played a connector role among the two modules by overlapping time, environment and HD transmission role between species. In addition, C. stellifer, as a module hub, connected the summer-fall module. All other Culicoides were peripheral, as they co-occurred more commonly with species displaying similar seasonality. Midge activity was not detected in five weeks in 2016 and just one week in 2017 (January and February, Table S1.5), when average minimum temperatures (2.16 ± 0.38 • C) were below the winter months' first quartile (Table S1.3, Figure S1.5). Finally, occurrence length and co-occurrence pattern among Culicoides linked vector persistence across the entire year, producing a suitable scenario to locally support the maintenance of arbovirus transmission.
Viruses 2020, 12, x FOR PEER REVIEW 11 of 26 was the steepest during spring months (interpolation curve). During spring, richness was saturated before rarefaction ended. However, in summer months, the sampling effort had to increase 17%, 31% in fall and 32% in winter to reach 98.5% sampling coverage. The 95% confidence interval is the shaded area, and significance between profiles was the inference of their nonoverlapping. Observed q0 and the reference or complete sample size (t) are indicated with The 95% confidence interval is the shaded area, and significance between profiles was the inference of their nonoverlapping. Observed q0 and the reference or complete sample size (t) are indicated with round dots and values detailed with bold legends. Diamonds, dots and grey legends show q0 and t necessary to attain 98.5% of the expected richness after coverage-based rarefaction and extrapolation analysis. displaying similar seasonality. Midge activity was not detected in five weeks in 2016 and just one week in 2017 (January and February, Table S1.5), when average minimum temperatures (2.16 ± 0.38 C°) were below the winter months' first quartile (Table S1.3, Figure S1.5). Finally, occurrence length and co-occurrence pattern among Culicoides linked vector persistence across the entire year, producing a suitable scenario to locally support the maintenance of arbovirus transmission. winter-spring: blue). Co-occurrence interactions between modules are drawn by a stronger opaque links pattern (violet) than within module interactions (transparent violet). Label size is relative to species degree (Table S1.4); vectors are red printed and plotted with circleplot [108]. (B) Species roles in the co-occurrence network sustaining and connecting modular configuration are shown in dimensions defined by z (within-module degree) and c (among-module connectivity) scores (

Discussion
The community assembly in Culicoides species followed several independent temporal profiles that were linked by environmental variables, neutral drivers and historical process. Importantly, the seasonal dynamics of Culicoides relied on a biannual phenology, highlighting that overimposed calendar schemes do not necessary fit the community responses [44]. For instance, northern Florida falls within the temperate domain of southeastern USA making it tempting to analyze communities using a four-season regime. The ordination pattern of Culicoides that we observed, however, is more consistent with a biannual seasonality structure commonly seen in subtropical and tropical biomes [39,[109][110][111]. This pattern of subtropical seasonality was consistently detected with our different modeling approaches (diversity composition, diversity structure, co-occurrence network and vector subassemblage partitioning among environmental thresholds).

Environmental Thresholds in the Temporal Dynamic of the Vector Assemblage
Clustering analysis defined weather thresholds as an important driver of vector assemblage ( Figure 7A). The resulting tree predicted 61% of the variation in the abundance with four informative nodes (Root node error = 0.46 (47.03/101); CV error = 0.67; SE = 0.09) and three environmental dimensions that were defined by stream condition and temperature lags (Figure 7). Precipitation was not selected by MVRT as a relevant source to predict variation in the vector assemblage. Nodes or clusters represented groups of independent epidemiological weeks which were allocated between nodes in the two sampling years. First, a minimum temperature of 13.6 • C (T min30 ) was able to explain 44% of the variation in the vector assemblage ( Figure 7A) and was explained by a warm node that spanned between the last half of the spring to the first half of fall. This node absorbed 97% of the abundance for C. debilipalpis (id = 0.6; p < 0.01) and 58% for C. stellifer (id = 0.57; p < 0.01; Figure 7A). Second, a 10% improvement in the clustering was explained by stream flow level (gauge 60 = 69.1 ft), at the end of fall, and represented 57% of C. venustus abundance (id = 0.57; p < 0.01). Finally, maximum temperature (T max45 = 23.8 • C) explained 7% of the vector temporal dynamic by a node related with some early spring weeks (97% in C. biguttatus abundance; id = 0.97; p < 0.01) and a consistent winter node (43% in C. pallidicornis abundance; id = 0.54; p < 0.01). Additionally, epidemiological weeks mismatched with calendar season after node allocation ( Figure 7B). Phenological dynamics of the vector community were dominated by arrays of species associated with defined ecological thresholds that were not congruent with calendar seasons.

Discussion
The community assembly in Culicoides species followed several independent temporal profiles that were linked by environmental variables, neutral drivers and historical process. Importantly, the seasonal dynamics of Culicoides relied on a biannual phenology, highlighting that overimposed calendar schemes do not necessary fit the community responses [44]. For instance, northern Florida falls within the temperate domain of southeastern USA making it tempting to analyze communities using a four-season regime. The ordination pattern of Culicoides that we observed, however, is more consistent with a biannual seasonality structure commonly seen in subtropical and tropical biomes [39,[109][110][111]. This pattern of subtropical seasonality was consistently detected with our different modeling approaches (diversity composition, diversity structure, co-occurrence network and vector subassemblage partitioning among environmental thresholds).

Seasonality of Culicoides β-Diversity
The winter-spring and summer-fall modules of species co-occurrence and niche partitioning of the vector assemblage coincide with the dry-wet scheme of the Oceanic and Tropical-Subtropical climate regimes of Florida [39,111]. Our data indicate that Culicoides species have a sequential and overlapping emergence that ensures temporal connectivity within the vector community. The early season, cold-tolerant species, transitional species and warm weather species may provide seasonal persistence of pathogen vectors, as supported by the harmonic overlap of species demographies in our models which identified recruitment waves of vector assemblages rather than any single dominant vector species. Multivoltinism in subtropical-tropical Culicoides [48] should increase overlap. Even though large pulses during the summer-fall months were observed (Figure S1.6, Table S1.5), we did not find biotic interactions, supporting the assertion that species coexistence relies on competitive processes [112].
Both temperature and stream flow shaped the likelihood of Culicoides species co-occurrence ( Figure 7). Gerry and Mullens [113] identified temperature rather than rainfall as a driver of seasonal dynamics of Culicoides in temperate regions. However, temperature usually interacts with rainfall in Culicoides dynamics [39,110,114,115]. Therefore, the distribution of the thermal niche between species subassemblages appears to contribute to the stability of BTV-EHDV transmission. The distribution of the thermal niche among species in the subassemblages is relevant (Figures 3 and 7, Figure S1.3) as both longevity and larval developmental times in Culicoides species are modulated by the species thermal sensitivity [116][117][118]. We propose that the fluctuations in temperature associated with the temporal change in Culicoides composition could also be affecting the overall load of EHDV-BTV transmission. For instance, the outcome of the virus-Culicoides interaction in the episystem (i.e., vector capacity, [23]) is sensitive to temperature because abundance, survival (demography) and extrinsic incubation period (virogenesis) covary with the temperature a Culicoides species experiences [116,[119][120][121][122]. Consequently, temperature has a strong net effect on Culicoides-borne Orbivirus.
The finding that stream dynamics and temperature drive seasonality in Culicoides β-diversity, but precipitation does not, indicates that nonprecipitation sources of water shape the moisture conditions driving larval development and population pulses of adults [37,[123][124][125]. In the panhandle of Florida, water from a superficial aquifer emerges in the slopes as seepages and runs off into streams, and clay soils cause surface pooling [51]. Consequently, the aquifer and stream interact with the Culicoides' niche. For instance, C. haematopotus, C. stellifer and C. venustus develop at stream borders, puddles and seepage habitats, respectively [126], and moisture and nutrient availability in these habitats are subject to the dynamics of stream flooding and aquifer surfacing [51].
The lack of significance of rainfall is surprising, given the intrinsic link between water and semiaquatic breeding sites of Culicoides [127] and the previously documented association between precipitation and Culicoides seasonality [39,110,114,115,128]. Nevertheless, precipitation was not related to Culicoides abundance in some studies [36,114,[129][130][131]. For instance, rainfall might be patchily distributed, but this spatial variability is expected to underestimate daily precipitation rather than the lag explored here. Importantly, our finding does not exclude the effect of rainfall on Culicoides but helps to explain how weather may constrain each population, community and functional group differently. As such, phenology proved less sensitive to rainfall, specifically, than to other measures of hydrology (stream flow).
The co-effect of temperature and stream flow configure the timing for community seasonal dynamics given the entire Culicoides species deploys a range of responses to moisture and temperature. For example, the least diverse assemblage (RDA 1 ) arose early in the epidemiological calendar ( Figure 2B), was mainly composed of winter-early spring species and was highly correlated with all of the selected stream and temperature variables ( Figure S1.3); meanwhile, the most diverse assemblage (RDA 2 ) arose later in the calendar, and was composed of a mixture of spring and summer-fall species exclusively correlated with temperature. We hypothesize that flooding early in the wet season resupplies the soil moisture lost during dry fall-winter months and nutrients drained by precipitation and recurrent flooding events during the wet season. In addition, it is interesting how the dynamics of the soil moisture and HD are linked. In the panhandle, the pattern of HD occurrence is supported by the enzootic stability hypothesis [69] as the relationship between high BTV/EHDV transmission rate and herd immunity [27] is associated with low mortality [69]. Christensen et al. [132] shows that the increase in drought severity has a protective effect on HD mortality occurrence within the counties in the top 20% of wetland cover and bellow 30 • latitude, which define our study location. Therefore, drought-flooding dynamics may lead to the enzootic stability of HD in association with Culicoides community change.
Our analysis showed that niche and neutral drivers partitioned equally Culicoides β-diversity seasonality, not from variability in flying activity and null noise in field data, as previously suggested [36,47,50,131,133]. Neutral drivers could be overestimated, as other variables (i.e., wind patterns, humidity, atmospheric pressure and daylight duration) were not measured but are expected to correlate with the environmental conditions explored here. It seems less likely that confounding neutral drivers, unexplored environmental layers and noisy data are major drivers of Culicoides β-diversity seasonality. Importantly, this so-called black box of "unexplained" variation in populations and community of Culicoides are susceptible to being biologically explained after stochastic processes are considered.

Seasonal Culicoides Diversity Structure
Organization in Culicoides diversity aligns with both nested and modular temporal topologies, relying upon the persistence of a few species (Figures 5 and 6, Figure S1.4). For example, we observed that 90% of Culicoides species occurred in weeks when species of another seasonal module were present (Figures 5 and 6A). As such, Culicoides richness and composition are temporally nested ( Figure 5) and may follow similar ecological rules. It also means that changes in diversity depend on how Culicoides species are lost over time but not by the replacement of some species by others. Community climax occurred in the spring weeks ( Figure 4) even though only 14% of the species displayed unique springtime phenology. In temperate regions, univoltine and rare species drive the high Culicoides richness in spring months [48]. Some Culicoides species might be ecologically redundant, as species exhibit temporal co-occurrence that significantly overlaps niche conditions, and they seem to share abundant resources [40,134]. Overall, this indicates that community structure is not organized by temperate calendar seasons, and the unified neutral theory of biodiversity is suitable to explain Culicoides diversity.

Persistence of Vector-Borne Pathogen in the Ecosystems Depends on Environmental Stability
Northern Florida provides a permissive environment to sustain the year-round transmission of Culicoides-borne Orbivirus as shown by vector persistence in the region. The longest observed period without host-seeking activity (5 weeks) does not signify arrested larval development, since C. pallidicornis adult emergence was continuous during winter weeks (Figure S1.6). The absence of Culicoides in traps during this period could be due to biological shifts in circadian rhythms, with normally nocturnal species shifting their host-seeking activity to warmer daylight hours [135][136][137].
Culicoides venustus and C. pallidicornis can be exposed to EHDV infectious blood meals as seroconversion and viremic resident white-tailed deer have been recorded in northern Florida during the winter. Alternatively, vertical and horizontal transmission, in addition to the persistence host infection, may facilitate the maintenance of EHDV-BTV during winter [26,[138][139][140]. Vector recruitment is nearly continuous for HD transmission, and our data suggest that C. venustus provides a link between fall and spring months. Nonetheless, no Culicoides species behaved as a network hub connecting all species and modules. Additional data on the infection prevalence of BTV and EHDV in Culicoides will help to understand the link between the continuous vector recruitment and the year-round persistence of HD. Operations such as animal movement and trade during winter months should proceed with caution as the vector-free period does not occur in northern Florida.

Is Culicoides β-Diversity Operating under Wider Temporal Scales?
Three wider temporal scales (≥12 months, Figure 2C) revealed that Culicoides communities experience annual and cyclic-periodic oscillation that could extend beyond two years. We propose that different processes operate between the wider scales. Here, we describe the annual variation at the Culicoides community level for the first time in Florida ( Figure 2; RDA 2,3 ) and that it is strongly related to the fluctuation in the environmental conditions. Importantly, the annual incidence of HD also correlates with environmental variables, especially drought severity, weather conditions and land cover change in the USA [25,132,133]. We speculate that the annual temporal scale is related to how species recruitment and the local environment are shaped by global-scale weather patterns. Both phases of El Niño Southern Oscillation (ENSO) took place during the study [141] and coincide with the observed effect of the temperature and stream level on annual Culicoides dynamics ( Figure 3A, Figure S1.3; RDA 2 ). Interestingly, ENSO has been previously implicated in Culicoides-borne episystems [142] where in southeastern USA, HD morbidity has occurred every 2-3 years [26,68,143].
When periodicity in community composition shows innate signals (RDA 1 ; Figure 2C,D, Figure S1.3) long-term biogeographical processes under wider spatial scales are expected to impact β-diversity. This emphasizes the trade-off between detrending data for satisfying mathematical assumptions and the risk of obscuring ecological patterns [11,87]. Altogether, Culicoides β-diversity agrees well with the space-time duality predictions [9]. In general, the wider scales are more dynamic and better suited to environmental filtering than seasonality. Moreover, temporal patterns in HD occurrence are strongly dependent on the spatial hierarchy since between years and periods, HD is nested at county and biogeographical scales [132,133]. Thus, we suggest that regional patterns in HD can also be predicted by surveying local scenarios of Orbivirus transmission in Culicoides community data.

Conclusions
The dynamics in the Culicoides community followed a subtropical seasonal regime driven by the atmospheric temperature and ground water early in the season and by temperature later in the season. The merit of this study relies on how we allowed the system to inform the phenological pattern instead of imposing a classical temperate seasonal scheme. Neutral drivers were relevant as niche drivers. Conversely, the relevance of the ecological drift on assembling was lacking, even though unexplained variation was related to noisy data niche-based studies.
The management of Culicoides-borne Orbivirus risk depends on disrupting the pathogen niche through seasonally-appropriate and targeted activities. In northern Florida, the disrupting of larval breeding during spring could be a valuable tool to interrupt HD transmission; but the use of stations for insecticide and repellent application may not be economically viable given the persistence and dispersal ability of Culicoides. Field evaluations are needed to assess these statements.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4915/12/9/931/s1. Supplementary File 1. Table S1.1: Temporal scales and their AEM eigenfunction associated in Culicoides composition change; Table S1.2: Species of Culicoides related with community change composition at each temporal scale; Table  S1.3: Weeks without Culicoides activity recorded and associated daily meteorology; Table S1.4: Temporal network of co-occurring Culicoides; Table S1.5: Monthly, year and total weekly Culicoides abundance, density and presence; Figure S1.1: Map of the study area and light traps distribution in Gadsden County, Florida; Figure S1.2: Culicoides composition time autocorrelation by multivariate Mantel's correlogram; Figure S1.3: Environmental constraints on Culicoides β-diversity temporal profiles; Figure S1.4: Culicoides season occurrence in the two modules derived from the weekly co-occurrence network; Figure S1