Farm-scale risk factors for bovine tuberculosis incidence in cattle herds during the Randomized Badger Culling Trial

Summary

We analysed the incidence of cattle herd breakdowns due to bovine tuberculosis (Mycobacterium bovis) in relation to experimental badger culling, badger populations and farm characteristics during the Randomized Badger Culling Trial (RBCT). Mixed modelling and event history analysis were used to examine the individual risk factors. The interdependencies of covariates were examined using structural equation modelling. There were consistent findings among the different analyses demonstrating that during a badger culling programme farms experiencing: reactive culling, larger herd sizes, larger holdings and holdings with multiple parcels of land were all at greater risk of a herd breakdown. Proactive culling reduced risks within the culling area, but we did not assess any potential effects in the periphery of the treatment area. Badger-related variables measured prior to the start of culling (number of social groups and length of badger territorial boundaries) did not consistently point to an increase in risk, when set against a background of ongoing badger culling. This could be because (1) the collected variables were not important to risk in cattle, or (2) there were insufficient data to demonstrate their importance. Our findings highlight the difficulty in identifying simple predictors of spatial variation in transmission risks from badger populations and the consequent challenge of tailoring management actions to any such field data.

INTRODUCTION

The management of disease in livestock is particularly challenging when a reservoir of infection exists in wildlife [1]. Management interventions in wildlife populations can be logistically demanding, while the behavioural and demographic responses of wild animals can give rise to unexpected outcomes [2]. Furthermore, as the behaviour of wild hosts will be modified by local ecological conditions, in heterogeneous environments their relative importance as sources of infection for domestic animals may vary in time and space [1]. Similarly, local variations in farm characteristics may influence both the risks of disease transmission among livestock and between wildlife and livestock. An understanding of such variation will be important in the development of sustainable approaches to the control of disease at the interface between wildlife and livestock. In order to be able to effectively target management efforts effectively it is necessary to identify predictors of disease risk that operate at a local scale, and to differentiate wildlife-related risk factors from farm-related factors.

Bovine tuberculosis (bTB) caused by Mycobacterium bovis is a serious disease of cattle, and constitutes a significant economic burden for the UK cattle industry and taxpayer. Eurasian badgers (Meles meles) are susceptible to infection and are recognized as a significant reservoir of infection for cattle. Since the 1970s badgers have been culled under a variety of strategies. However, the effectiveness of badger culling as a tool for controlling the risks of transmission to cattle remained unclear. In 1998, following the recommendations of an independent review [3] the UK Government initiated the Randomized Badger Culling Trial (RBCT). The trial had a simple design with 10 ‘triplets’ each of which comprised three 100 km² treatment areas. In ‘Proactive’ treatment areas badgers were culled approximately once per year on all accessible land, in ‘Reactive’ treatment areas badger culls only took place on farm holdings in response to the detection of infection in cattle, and in ‘Survey-only’ control areas no culling took place [4].

Analyses of the results of the RBCT at the ‘triplet’ scale, i.e. each 100 km² area constituted a single replicate, confirmed that disease incidence in cattle was significantly affected by badger culling. However, the outcomes of culling were complex. Within proactive culling areas, disease incidence in cattle was reduced; while in the immediate perimeter and in reactive culling areas, the incidence of infection in cattle increased [5]. These patterns changed over time, such that the positive effects of culling lasted up to 42 months and the negative effects diminished [6]. These outcomes were associated with perturbation of the social structure of badger populations [7] which was in turn associated with increased risks of transmission among badgers [8, 9] and, it is hypothesized, increased risks of onward infection of cattle [10]. Although local predictors of the magnitude of such effects have yet to be described in any detail, they are likely to relate to prevailing ecological conditions and the characteristics of the culling operation [11].

Previous analyses of data arising from the RBCT have focused on identifying differences in the incidence of infection in cattle in relation to badger culling at relatively large spatial and temporal scales, i.e. treatment triplets (i.e. 100 km²) and over the duration of the trial (8 years) and thereafter. While this is appropriate given the original experimental design and the aim of evaluating badger culling as a disease control option for potential application at the landscape scale and in the long term, it may be of less value in understanding the epidemiology of the disease in the more local context of the herd or farm holding and in the shorter term. Indeed, other analyses of RBCT data have revealed complex patterns of spatial and temporal effects on trial outcomes [5, 12, 13]. A 100 km² sampling unit is too coarse to capture ecological processes that may operate at the scale of the individual badger social group, or farm-holding-related characteristics that may influence local transmission risks.

The RBCT took place over a period of 8 years during which testing of herds on farms in the triplet areas was undertaken on a regular basis. This provides a longitudinal record of outbreaks of disease at the level of the farm holding, coupled with information on the distribution of badgers. Analysis of disease incidence in cattle herds in relation to risk factors measured at the same scale might provide further insights into the role of farm-holding characteristics, local landscape and badger distribution. Furthermore, the longitudinal nature of the study allows partition of fixed effects from random unmeasured effects. Hence we carried out analyses on a subset of RBCT data at the scale of the County Parish Holding (CPH or ‘holding’), which is the unit at which infection in cattle herds is recorded. We specifically investigated events in the trial areas, i.e. we did not analyse events in the periphery of culling areas, and restricted our analyses to the period during culling, i.e. not before or afterwards. We considered the roles of farm-holding and cattle herd characteristics, habitat composition and the distribution of badger populations at the start of the trial as risk factors. Although some work has been carried out on herd-based risk factors for bTB incidence in cattle, the present study represents the first attempt to simultaneously investigate the role of factors relating to cattle and to wildlife at the scale of individual farm holdings. Further understanding of risk factors that operate at the scale of the farm holding is important for the development of approaches to manage M. bovis transmission at the interface between cattle and badger populations.

MATERIALS AND METHODS

Data

The data used in the present analyses were derived from the RBCT and were provided by the Veterinary Laboratories Agency (VLA). Owing to various limitations on the availability of data at the farm-holding scale (see Supplementary material, available online), we confined our analyses to data from four of the 10 triplets. The selected triplets were located on the Devon/Cornwall border (Triplet B), in East (C) and West Cornwall (F), and on the Somerset/Devon border (H). For the same reasons analyses were restricted to cattle herd breakdowns (definition below) that started between 1 January 2002 and 31 December 2005.

Infection in cattle

The incidence of bTB infection in cattle is determined routinely in the UK by the application of the tuberculin skin test [14]. Regular testing of cattle herds occurred on all CPHs on an annual basis during the RBCT. The presence of at least one confirmed positive test result (a reactor) is indicative of an incident event, known as a cattle herd breakdown (CHB). If a breakdown was recorded, then that herd was subsequently tested more frequently, as were those herds on any contiguous holdings. This gave rise to a sequence of test results associated with each CHB event. In our analyses we only included the first positive test of any given confirmed CHB event (i.e. the index case). Negative herd tests were only considered for a holding if they occurred outside the start and end dates of the sequence of tests associated with a given index case.

Covariates

Most of the covariates in the analyses originated directly from the available RBCT data for each CPH, although some were derived from manipulations of spatial data (see Supplementary material for details). Spatial data, stored in ESRI shape file format, were interrogated to provide measures of the number of parcels of land associated with each farm holding, the number of neighbouring holdings and the distance to the nearest herd placed under restriction due to an ongoing breakdown. In addition, we used spatial data on the approximated distribution of badger social group territories measured at the beginning of the trial (and prior to the initiation of culling treatments) in relation to each farm holding. Each badger social group generally occupies a territory, with boundaries between groups being characterized by the presence of latrines [15] containing potentially infectious excretory products. For the data examined here, the extent and configuration of badger social group territories was determined by a combination of bait-marking [16] and expert opinion based on field signs, and these home ranges were stored in a GIS. In our study, we estimated the total number of badger social group territories and the total length of territorial boundary present on each CPH by overlaying the badger territories and CPH shape files. UK Landcover 2000 data [17] was used to derive measures of habitat composition for each CPH (see Supplementary material). Since woodland and scrub is not generally used by cattle but is an important habitat for badger foraging and for sett location in particular [15], it was hypothesized that this may relate to the number of badger social groups present on a holding. Hence we estimated the percentage of woodland and scrub cover present in each CPH. Moreover, as managed grasslands provide excellent foraging habitat for badgers (see e.g. [18]) and hence opportunities for direct and indirect contact with grazing cattle, we also estimated the percentage of grassland present in each CPH. There was a clear association between these variables, so only the proportion of woodland was included as a covariate in the models.

Data analysis

General

A progressive modelling strategy was used to initially identify risk factors for disease at the level of the farm holding, then to investigate spatio-temporal variations in herd breakdown events, and finally to describe the interactions among potential risk factors. Generalized linear mixed models (GLMM) were used to analyse the influence of variation in treatments and covariates and their interactions on the subsequent incidence of infection in cattle. Then to investigate the time to detection of infection in cattle herds on holdings under each treatment we used event history analyses. While event and mixed-effects modelling approaches are useful for analysing the relative contribution of putative risk factors for disease incidents individually and through time, they ignore the fact that risk factors themselves may be inter-dependent. Hence, structural equation modelling (SEM) was used to investigate the relationships among different processes. SEM and path analysis have been previously used to investigate complex ecological systems [19, 20] and disease transmission in social environments [21]. In this study SEM was used to develop a conceptual model of how each of the risk factors considered in the mixed modelling analyses may interact and precipitate disease incidents. This model was then populated with the RBCT data to investigate the incidence of bTB in cattle herds during the trial.

CHB incidence constituted the response variable for our analyses. For mixed-effects modelling it was assumed that each CHB represented the result of a binomial test of whether the herd had been detected positive at the time of testing or not. This was also assumed for the event history analyses. The response variable for SEM was the proportion of positive tests for a CPH tested throughout the study (2002–2005).

Mixed-effects modelling

GLMMs were used to investigate variations in CHB risk on individual holdings in relation to farm-level covariates and triplet-level treatments (culling strategies). The farm holding was modelled as a random effect, initially nested within triplet. The underlying hypothesis was that measures of farm, landscape and badger distribution characteristics, proximity to other CHB events and season were drivers of CHB risk on individual holdings (Table 1). We also investigated the interactions between triplet and treatment. To investigate the extent of seasonality in the data, harmonic covariates were included within the model. Since CHB was effectively a binomial response models were fitted using penalized-quasi-likelihood (glmmPQL) in the MASS library [22] within the statistical package R [23].

Table 1. Descriptions of spatial covariates measured at the level of the County Parish Holding (CPH) used in each analysis

Event history analysis

The impacts of treatment, habitat, badger-related covariates, and farm-holding characteristics on the incidence of CHBs were investigated using event history analysis. We defined a CHB event as being test-positive to bTB. Here, it was assumed that the hazard for CHBs was proportional and dependent on measured covariates. We developed models that assumed a baseline hazard common to all farms in all triplets and used relevant covariates in Cox proportional hazards models. Herd breakdowns were not recorded at the same time that an individual animal became infected, instead herds were tested at set intervals. This means that the events could not be ascribed to the point at which they occurred, rather the data were effectively interval-censored. Interval censoring leads to bias which impacts on both the direction and magnitude of the effects of covariates in the Cox model, the magnitude of which cannot be predicted a priori. We used the iterative convex minorant algorithm (ICM) as implemented in the intcox library [24] in R to estimate the regression coefficients of the Cox model for the full dataset. Since this procedure cannot estimate standard errors for the regression coefficients, we used a bootstrapping procedure to create upper and lower bounds on 999 samples of the original data. For this we bootstrapped by farm holding rather than individual breakdown event, following the rationale in Therneau & Grambsch [25]. In order to assess the impact of each covariate we counted the number of samples for which the estimated coefficient was less than zero if the mean for all coefficients was greater than zero and conversely the number of samples greater than zero if the mean was less than zero. This effectively estimates the probability that a regression coefficient would encompass zero, where it could be assumed that the covariate was not having a significant effect on the hazard, in analogy with the Wald test.

SEM

Path analysis was used within a SEM framework to investigate postulated hypothetical models of the causal relationships between landscape, badger and holding-related characteristics and the occurrence of CHBs on individual farms in each triplet. The proportion of cattle tests on each CPH that were positive over the sampling period was used as the response variable. The full conceptual model (Fig. 1) had three exogenous variables, which were the treatment (as imposed in the RBCT), the area of a CPH, and the proportion of woodland present. Triplet was not included in this analysis because it would have led to a further four categorical variables and it was hypothesized that any differences between triplets should have been encapsulated in the landscape characteristics which were more proximal drivers of the CHB process. In effect, the exogenous variables represent the underlying differences in the landscape and the treatment in each triplet, and uncontrollable aspects of holding size. It was hypothesized that these were drivers of other variables, specifically farm characteristics relating to the size of herd present and the presence of contiguous holdings. It was proposed that proportion of woodland, holding area and the number of land parcels would determine the number of badger social group territories recorded on the holding and that the more groups there were on a holding the more boundaries (between groups) would be present. We hypothesized that the risk of a CHB on an individual holding would be dependent on the number of badger groups present and the length of territorial boundaries present, since this would lead to increased contact between livestock and an environment potentially contaminated by infected badgers. In addition, it was assumed that the greater the number of holdings with which any CPH was contiguous, the greater the risk of infection from cattle-to-cattle transmission at their boundaries. Treatment would impact directly on CHB risk through an effect on the total number of badgers in overlapping social groups.

Fig. 1. Conceptual model of the relationships between potential risk factors and cattle herd breakdowns. CHB, Cattle herd breakdown.

Models were assessed by comparing comparative fit index (CFI) criteria, χ² measures of association and the root mean square error of approximation (RMSEA) for each model, we identified the most parsimonious model for the available data [26, 27]. Models were fitted using maximum likelihood in Mplus [28] following the approach used by Rushton et al. [21] where the formulation assumes conditional normality rather than the more restrictive assumptions of multivariate normality. The conditional normality assumption allows non-normality for the response variables [28] as might be expected when modelling proportions or categorical responses.

RESULTS

General

There were a total of 1309 holdings in the four triplets sampled over the 4-year period. Cattle herds were tested on 10994 occasions, and of these 561 resulted in a confirmed CHB. Summary statistics for the covariates in the treatment areas of each triplet that were used in the modelling are shown in Table 2. There was considerable variation in the number of holdings, their size, and number of land parcels across the four triplets, with holdings being smallest in Triplet F.

Table 2. Characteristics for County Parish Holding (CPH) level covariates across the four triplets in this study

CHB, Cattle herd breakdown.

Mixed-effects modelling

The parsimonious model (Table 3) indicated that culling treatment, CPH area, contiguity with other holdings and the presence of a breakdown within 1 km were significant predictors of CHB risk on individual CPHs. Interactions between triplet and treatment were not significant. Holdings subjected to reactive culling and survey-only treatments experienced higher CHB rates (23% and 18%, respectively) than those subjected to proactive culling (Table 3). There was also evidence that CHB risk diminished over the period of study by about 2% per year. The number of badger social groups, length of badger boundaries and the proportion of woodland occurring on a CPH were not found to be significant risk factors. The number of parcels of land comprising a CPH was also not significant.

Table 3. Binomial mixed-effect model fitted by maximum likelihood investigating cattle herd breakdown risk on individual holdings in relation to badger culling strategy and holding-level covariates

Event history analysis

A sample of 438 farm holdings that experienced a CHB event during the period 2002–2005 were included in these analyses. However, key covariates were not available for all CPHs so sample sizes were reduced accordingly.

Triplet was not a significant predictor of risk of breakdown, with the exception of Triplet F, which had a lower hazard rate (Table 4). The parsimonious model for first CHB events suggested that holdings in reactive-culling and survey-only treatment areas experienced breakdown risks that were slightly lower (4% for both) than those in proactively culled areas (Table 5); however, these differences were non-significant. Of the farm-scale covariates herd size slightly reduced risk of breakdown (Table 5), the proportion of woodland on a holding was shown to reduce the risk of CHB by 40%

Table 4. Interval-censored event analysis, mean regression coefficient and s.d. from 999 bootstrapped samples

Column n<>0 is the number of samples for which the coefficient exceeded or was exceeded by zero. A value of 50 would approximate to 95% upper or lower intervals on the mean.

Table 5. Mean regression coefficient and s.d. from 999 bootstrapped samples for a parsimonious interval-censored event model from which covariates have been removed

Column n<>0 is the number of samples for which the coefficient exceeded or was exceeded by zero. A value of 50 would approximate to 95% upper or lower intervals on the mean.

SEM

An initial χ² test was used to compare deviations of the coefficients of the model from those expected from the covariance matrix. There was no evidence that the full model deviated from the expectations of the data as assessed by the χ² statistic but there is considerable doubt as to the validity of the model when considered in the context of other statistics. Not all of the hypothesized drivers of CHBs were shown to be statistically significant when assessed using the approximate standard error for each parameter, suggesting no causal link. The most parsimonious model (Fig. 2) indicated that CHB risk was related more strongly to farm-holding attributes than to landscape or badger-related parameters (Table 6). CHB risk was positively related to farm-holding size, the number of contiguous neighbouring holdings and cattle herd size. The number of badger group territories on holdings was, as hypothesized, greater on farms with a higher proportion of woodland. However, neither the treatment (proactive, reactive, survey) nor the numbers of badger social groups present on each holding were drivers of CHB risk in this model.

Fig. 2. Parsimonious structured equation model and related risks. Arrows indicate the direction of significant paths and associated coefficients and error. Model fit: root mean square error of approximation=0·05; comparative fit index=0·994. CHB, Cattle herd breakdown.

Table 6. Parsimonious structural equation model diagnostics

RMSEA, Root mean square error of approximation.

DISCUSSION

Bovine TB infection in UK cattle herds has a complex aetiology. Potential risk factors associated with cattle, badgers and their interactions may vary in space and time, and in relation to management interventions. The results of the RBCT and associated research highlighted the complexity of this multi-factorial problem and the challenges of developing sustainable management solutions [29]. In the present study we used data originating from the RBCT to conduct multivariate analyses of farm and wildlife-related risk factors that were associated with the incidence of bTB in cattle during experimental culling and which were acting at the scale of the farm holding.

As in all complex statistical studies, many of the approaches adopted make assumptions of the data that are not fully met and we have to consider how these assumptions impact on our analyses. In such cases the significance of the results can only be considered indicative. Many statistical tests assume that the data are normal, and for biological systems this assumption is rarely met. Nonetheless, the results from such testing are considered adequate, as the most likely consequence is that the ‘real’ P value is smaller than that calculated. In the mixed-effects modelling the data had non-normal error structures and the penalized quasi-likelihood approach we used has been considered to be unreliable and potentially biased for large variance [30]. There are also likely to have been both spatial and temporal components to the error structure that should be considered in the modelling, given the landscape context. Trivially, badger social groups were found on more than one CPH, so some of the predictors were necessarily non-independent. More critically, we might expect spatial dependence in CHB events between CPH (e.g. through the mechanism of cross-farm boundary transmission and/or farm–badger–farm transmission) and temporal correlation between events (disease persists on the farm although not in the cattle). These issues cannot be addressed in the mixed-effects modelling framework developed here.

In the event history analyses, events were recorded at rather wide intervals relative to the time domain of the study. Holdings were tested for breakdowns at regular intervals but the exact date of an event occurring was not known hence the data had to be analysed taking account of interval censoring. The events as recorded were unlikely to be independent in the first place, since the occurrence of a positive test result for M. bovis would trigger a series of subsequent tests on both the farm and locally. While we did not include repeat testing of the same CHB event within the analyses, occurrence of an event on one CPH would have led to increased sampling in local farms, potentially emphasizing the relative importance of local CHB events as predictors of events on any of the local CPHs. One aspect that the analyses highlighted was the complexity of the disease process. The SEM approach allowed us to investigate the multiple dependencies between covariates. This dependency could lead to confounding in a mixed-effects or event modelling context, and should be incorporated in any future models if we are to understand the pattern of CHBs in the RBCT.

These concerns notwithstanding, there was a high degree of concordance between the results of the different analytical approaches used. Badger culling treatment was shown to have an effect in both the mixed modelling and the event analysis. In the former, risk was shown to be greater in the reactive and survey-only areas, which confirms that a proactive culling regimen reduced CHB events for individual farms within the culling area in the years 2002–2005. This finding is entirely consistent with earlier landscape-scale analyses both in Ireland [31, 32] and in the RBCT [6, 10, 12], although clearly does not take into account the counteracting detrimental effects detected in the periphery of proactive culling areas [6]. The background risk was lower in Triplet F in the mixed-effects models and the event analyses, but no obvious reason for these differences was identified.

The most consistent feature of the study was that farm characteristics, specifically herd and farm size, number of land parcels and contiguity were significant risk factors in determining CHB. This confirms the results of previous studies (e.g. herd size [33], multiple premises [34], holding size [35]). These relationships may also be indicative of indirect effects arising from other risk factors related to cattle purchase, housing type and grazing systems, which have also previously been found to have a significant effect on herd breakdown rate [34] and will be linked to farm size.

No consistent association with CHB risk was found for habitat and badger-related variables at the level of the individual premises, although deciduous woodland has previously been found to be significant [5]. There are several possible explanations for this. First, the data used may not have been ideal for testing the underlying hypotheses relating CHB to badgers, since badger social group data were collected at the entry of each RBCT triplet and CHB events were recorded up to 7 years later. It is well established that culling caused substantial perturbation to the social structure of badger populations [35–37], and that the suspension of cattle controls during a foot-and-mouth disease epizootic increased bTB prevalence in badgers [10]. Therefore it is likely that the similarity between the badger covariates measured at the start and those in place at the time CHBs occurred diminished as time passed and that the available badger data did not represent a strong ‘surrogate’ measure of the true state of the badger population and hence of the risk of bTB in badgers. Second, this lack of any clear relationship may be due to limited sample size, suggesting that the relationship between CHBs and badger variables is weaker than that with farm characteristics. Third, it is possible that the badger variables were poor representations of CHB risk, since there is an implicit assumption in using the number of badger social groups as a predictor that the disease threat posed by badgers would be higher in higher density badger populations. However, bTB has been shown to persist in spatial clusters in badger populations [38] and disease prevalence may be higher in smaller social groups [39]. It is important to note that results from the mixed modelling suggest that badger culling had an effect on CHB rates for some time after implementation, thus supporting similar findings in previous analyses of the RBCT [40].

In conclusion, our farm-scale analyses were in agreement with earlier landscape-scale work indicating the detrimental effects on cattle disease risk associated with localized, reactive badger culling. We have also confirmed the positive effects of proactive badger culling at the farm scale within the culling area, although it is important to recall that we did not analyse the experience of farms in the periphery of the culling areas. We have confirmed earlier work showing that farm-holding characteristics were important factors in the occurrence of CHBs. The results of the event analyses, mixed-effects modelling and SEM were consistent in suggesting that large holdings with multiple parcels or large herds were more likely to suffer CHB than smaller holdings [33, 41]. No such consistency was found in the badger-related covariates as predictors of CHBs. This highlights the difficulty in identifying simple predictors of spatial variation in transmission risks from badger populations and the consequent challenge of tailoring management actions to any such field data.