Bovine anaplasmosis, caused by the rickettsia Anaplasma marginale, is an economically important tick- and blood-borne disease of cattle found worldwide. Its clinical effects of severe anemia, weight loss, spontaneous abortions, and death negatively impact cattle welfare and create a significant economic burden for cattle producers (1). A vaccine that is consistently safe and effective has proved elusive (2). Treatment and control usually depends on administration of oral and/or injectable tetracyclines; however, this approach is not always effective (3) and may be associated with the development of antimicrobial resistance (4). Cattle that survive acute infection are lifelong subclinical carriers of A. marginale and represent a key reservoir of infection for other cattle through tick transmission, inadvertent mechanical transmission via blood-contaminated needles (5) and veterinary instruments, and mechanical transmission via certain biting flies (6, 7).

Evidence suggests that bovine anaplasmosis is becoming more prevalent and more geographically widespread (8–13). Spatiotemporal modeling demonstrated increases of anaplasmosis in Kansas since 2005 (14), where over 50% of beef herds are infected (15). Human tick-borne illnesses follow a similar trend, having increased significantly over the last decade (16), particularly in the midwestern U.S. Concurrently, expansion of the geographic distributions and active seasons of medically important ticks in North America has been well documented (17–20), while human encroachment on tick habitats increases the likelihood of encounters (21).

In the U.S., bovine anaplasmosis is vectored by Dermacentor spp. including the American dog tick (D. variabilis), the Rocky Mountain wood tick (D. andersoni), the winter tick (D. albipictus), and the western dog tick, D. occidentalis (22). Both D. variabilis and D. albipictus are established in Missouri, although D. albipictus preferentially feeds on cervids (23). Distribution of tick-borne diseases, such as bovine anaplasmosis are often concordant with the spatial distribution of their transmitting vectors (17, 24, 25), and are therefore indirectly affected by the natural fluctuations in climate and habitat availability for ticks and their hosts. In the case of bovine anaplasmosis, additional factors that may affect spatiotemporal distribution include management practices (5, 26) and cattle movements leading to comingling of naïve and carrier animals.

Spatiotemporal disease mapping models are a useful tool to describe such disease patterns, more specifically to identify the presence of any clusters of incidences over space, time, or both (27). The presence of such clusters often points to underlying management, environmental, or other factors that could be potentially driving the disease incidence, which can then be targeted for interventions, and resource allocation (28, 29). Making inferences based on crude values of disease incidence over geographic space and time is problematic as the underlying spatiotemporal structure of the data is unaccounted for (27, 30). Bayesian spatiotemporal models provide a robust and flexible platform for space-time analysis of disease incidences and have been utilized to explore the potential effects abiotic covariates such as climate and land cover features in the disease ecology. Such models have been applied to the study of bovine anaplasmosis and other tick-borne diseases (14, 24, 25).

The objective of this study was to evaluate the spatial, temporal, spatiotemporal risk patterns of bovine anaplasmosis in the Midwestern state of Missouri using a Bayesian approach with counts of bovine serum samples that tested positive via competitive ELISA (cELISA) at the University of Missouri Veterinary Medical Diagnostic Laboratory (MU VMDL) between the years 2010–2021.

Bovine anaplasmosis is endemic in Missouri, with an estimated statewide seroprevalence of 46.5%. (37). This study analyzed the current spatial and spatiotemporal patterns in the disease's distribution across the state using retrospective diagnostic data routinely recorded in a laboratory information management system (LIMS), often a source of information for making epidemiological and disease management decisions (38, 39). The Bayesian models considered in this study were data-driven, and the better performing model (Model-3) utilizes the CAR autocorrelation, which is effective in borrowing information from neighboring space-time units (27, 40–42). A strength of the models in this study is their ability to accommodate missing counts as well as missing sampling sizes, which alleviates the missing data problem often encountered in passive surveillance datasets. Such model structure is widely applied in site-occupancy models for accommodating local population uncertainties (33), and could be useful in situations similar to ours wherein LIMS data are assessed. The risk calculated in this study, i.e., the inverse logit of Bayesian posterior means for the purely spatial, purely temporal, and spatiotemporal random effects fitted in Model-3 reveal the spatial, temporal, and spatiotemporal trends for bovine anaplasmosis seropositivity over the 12-year study period.

The overall structured spatial heterogeneity (Figure 4), when classified in three categories (low, medium, and high risk), shows clusters of high-risk counties for bovine anaplasmosis seropositivity along the eastern half of the state, in the north, central and southeast, while the medium and low risk counties are found throughout, and are also found in contiguous groups of counties. This suggests that the underlying spatial processes that drive such clustering are likely to be similar or potentially arise from the same source, i.e., environmental, disease management practices by producers and veterinarians, or other similar factors. The overall unstructured spatial heterogeneity on the other hand (Figure 5) does not reveal any clustering, suggesting that there is no consistency in the effect of any local (or county-level) factors that play a major role in the spatial dynamics of bovine anaplasmosis seroprevalence in Missouri. The temporally structured and unstructured heterogeneity that accounted for the temporal effect in the retrospective dataset are presented in Figure 6, revealing that the risk of observing bovine anaplasmosis seroprevalence in Missouri has changed over the study period, with some years recording higher risk than others. While the trend is non-linear, it appears the disease burden is more or less stable over the years in Missouri, with some noticeable interannual variability.

Choropleth plots showing posterior risk of structured spatiotemporal heterogeneity? (Figure 7) indicates that the variations in presence of seropositive cases across Missouri counties are not independent but are spatially correlated, and shows a varying trend, with the disease risk appearing to expand and contract over the years. Disease risk increased significantly between 2010 and 2012, peaked around 2016, and again in 2019. The northeastern and southeastern regions consistently show higher risk throughout the years. After 2019, a decline in risk is noticeable; however, this trend could be likely due to less information available to the model that could be derived from the dataset for the latter years. Additionally, values outside these clusters change abruptly, with high and low-value counties scattered throughout the study area, sometimes in proximity. This pattern indicates that while the model successfully identifies regions potentially with similar characteristics affecting presence of seropositive cases, other unknown or unexamined factors may also contribute to the spatial distribution of the disease risk.

The risk maps of structured spatiotemporal heterogeneity further reveal that certain counties remain high-risk over multiple years, indicating persistent high-risk areas. These regions likely have stable environmental factors conducive to disease transmission, such as tick-favorable habitats or frequent livestock movement. Conversely, some counties consistently appear low-risk, representing stable low-risk areas where environmental or management conditions may naturally limit anaplasmosis transmission. Additionally, several counties exhibit fluctuating risk levels over time. Some areas increase in risk in later years, while others decline in risk. These shifts may result from interannual variations such as El Niño Southern Oscillation effects, seasonal variations, land use changes, or changes in disease management practices. Furthermore, emerging high-risk zones are evident, as some counties that were previously lower-risk become high-risk over time. This expansion of high-risk areas may be driven by factors such as climate change, shifts in livestock distribution, or increasing tick populations.

In contrast, the unstructured or spatially unstructured heterogeneity (Figure 8) captures variations in anaplasmosis prevalence that are independent across counties, without spatial correlation to neighboring areas. Higher values suggest that anaplasmosis risk in those counties is driven by unique factors rather than broader spatial processes that operate at larger geographic scales spanning multiple contiguous counties. While values vary widely across the study area, no distinct clusters emerge. This indicates that although global trends exist, local deviations also occur, where specific local level factors may influence individual counties differently. Finally, the spatiotemporal plot of Bayesian p-values for exceedance probability (Figure 9) indicates that while Model-3 performs well for a vast majority of the counties in Missouri, there is a consistent set of counties in Missouri where the model underperforms. One of the reasons for this underperformance could be related to fewer submissions from these counties to the VMDL, which are not adequately addressed by the missing data imputation method used. Additionally, this is also potentially due to the existence of unique factors that drive seropositivity risk in these counties, which are not adequately captured by the random effect terms included in Model-3. This can be alleviated by including biological (environmental, physical), and disease management-related variables that may capture the underlying spatiotemporal process more efficiently. Although our models did not directly evaluate the effects of environmental or management covariates, the spatial and temporal patterns observed in Missouri suggest that factors such as the distribution of wooded vs. open pastureland, regional variation in climate and vector habitat suitability, and differences in herd management practices (e.g., grazing intensity, herd movement, and vector control) may help explain variation in anaplasmosis risk. Incorporating these data into future analyses would allow for more precise assessments of how such factors shape disease dynamics and could guide the development of region-specific management strategies to reduce the burden of bovine anaplasmosis.

Clusters of high-risk counties scattered throughout the state (Figure 4) suggest the presence of similar underlying processes that operate across local and/or regional scales and which are not captured in the present model. Possible explanatory variables include local variations in microhabitat and host abundance, which can influence populations of vector ticks, such as soil characteristics (43), encroachment by woody plants such as invasive red cedar (44), and variations in population density of white-tailed deer (45). Local differences in cattle population and movement are also likely to be involved. The lack of any statewide trend in plots of spatially unstructured heterogeneity (Figures 5, 8) further illustrates the importance of local risk factors to anaplasmosis epidemiology in Missouri.

While the overall risk of bovine anaplasmosis seropositivity in Missouri appears relatively steady during the 12-year study period (Figure 6), some fluctuation is nonetheless apparent. While some of this fluctuation may be attributable to random “noise” in the available data, other potential explanatory variables not captured in the present model include temporal fluctuations in weather patterns, which in turn can influence cattle population dynamics through their effects on market pressures. A specific example is the marked reduction in beef cow inventories prompted by drought in 2011–2013, which reached its nadir in 2014 (46, 47). This was eventually followed by restocking of herds, often with replacement heifers shipped in from out-of-state, which increases the opportunities for comingling of naïve and carrier cattle (48). Temporal fluctuations in weather patterns are also likely to affect tick populations and host-seeking behaviors (49, 50).

Plots of county- and year-specific spatially structured heterogeneity (Figure 7) reveal four distinct situations: (1) counties that are consistently higher-risk, (2) counties that are consistently lower-risk, (3) counties that transition from lower to higher risk, and (4) counties that transition from higher to lower risk. Illustrated differently, if the faceted plot in Figure 7 were animated as a time series, the color gradient of certain counties would become progressively darker (increasing risk) or lighter (decreasing risk), sometimes multiple times during the study period. This is uncharacteristic of biological diffusion, i.e., “natural” spread of disease cases over space and time (51) and instead suggests introduction of cases via cattle movements, similar to the manner in which human diseases can appear in new locations via travel (52). Cattle movement, generally through markets, is strongly suspected to be an important source of variation in anaplasmosis risk around Missouri.

Finally, spatiotemporal plots of Bayesian p-values for exceedance probability (Figure 9) indicate 13 counties where the model performs poorly. Some are likely attributable to insufficient submissions from these counties to the VMDL; the independent city of St. Louis, which, as expected for beef cattle, has no submissions during the study period, is perhaps the best example. Nonetheless, several other counties (such as Dunklin and New Madrid, in the “bootheel” of southeastern Missouri) are also without submissions during the study period and are not flagged as underperforming via exceedance probabilities. Further work is required to understand why Model-3 performs poorly in some counties. One potential way to improve the model performance is by including fixed effect terms that may have further explanatory potential, such as those in the categories of environmental, herd management, and tick distribution. The availability of such data at the scale of present analysis and their suitable spatial transformation may however be limited.