Bovine tuberculosis breakdown duration in cattle herds: an investigation of herd, host, pathogen and wildlife risk factors

Study area

Northern Ireland is approximately 14,000 km². The official bTB control programme is administered over ten Divisional Veterinary Office DVO areas, comprised of 123 “patches”; mean size 110 km² (SD ± 53); Fig. 1.

Study design

Two retrospective analyses were undertaken, firstly (i) quantifying the risk factors associated with bTB breakdown duration, using negative binomial (count) regression with the outcome measured in days, and (ii) quantifying the risk associated with bTB breakdown duration using ordinal regression, with the outcome modelled as a categorical ordered variable. This approach was considered necessary to account for bTB breakdown administration in NI. Herd breakdown duration measures arise as a result of a disease management process, and are not a wholly natural phenomenon. Generally, once bTB has been confirmed, the herd Officially Tuberculosis free status is Withdrawn (OTW). Usually, two clear herd tests are required to restore Officially Tuberculosis Free (OTF) status. Each herd-level test is scheduled to occur a minimum of 60 days apart (DAERA, 2018). For the ordinal regression therefore, breakdown duration was classified into four distinct categories based on multiples of 60 days (it should be noted, however, that breakdown length may not always correlate exactly with the number of tests done, as herds may delay testing). The first category contained breakdowns ≤180 days (approx. 6 months; 3 tests until OTF status restored), the second category included breakdowns which ended up to 120 days later; ≤300 days (approx. 9 months; 5 tests until OTF status restored), and the third category included breakdowns which ended up to 120 days after this; ≤420 days (approx. 13 months; 7 tests until OTF status restored). The final category included breakdowns which lasted longer than 421 days (8 or more tests until OTF status restored). Breakdown start dates are denoted by the date at which the first SICCT reactor or lesioned animal identified at slaughter was disclosed, and the breakdown end date was the test at which the last clear herd test was achieved.

Dataset creation

BTB breakdown data spanning January 2003 to December 2015 inclusive (n breakdowns = 27,718) were made available from the NI Department of Agriculture, Environment and Rural Affairs (DAERA) database, the Animal and Public Health Information System (APHIS) (Houston, 2001). This dataset was restricted to only include OTW breakdowns (n = 19,084; 8,634 breakdowns removed), which were defined by policy guidance at the time of study as the presence of more than five SICCT reactors, or two positive results to the four possible bTB tests; confirmation via histopathology, culture or spoligotyping, or the identification of a lesion at routine slaughter. Breakdowns with incomplete or erroneous information were also excluded (e.g., missing GIS information, MLVA information, or breakdowns lacking end dates (n = 17,114; 1,970 breakdowns removed). The dataset was further restricted to include breakdowns which started and ended between 01∕01∕2009 and 31∕12∕2014 (n = 7,478; 9,636 breakdowns removed). These dates were chosen because surveillance using M. bovis MLVA genotyping data occurred at the herd level between 2003 and 2008, but from 2009 onwards, all culture confirmed animal-level M. bovis isolates were genotyped. Finally, breakdowns which were recorded as lasting less than 60 days were excluded from the final dataset (n = 5 breakdowns removed), as 60 days is the minimum restriction period which may be permitted under some circumstances e.g., less than five positive SICCT animals with no post-mortem or laboratory confirmation (DAERA, 2019). The final dataset contained information on 7,473 breakdowns. All data were assembled and analysed using Microsoft Access 2007 (12.0.6735.5000) SP3 MSO and R Version 3.2.5 (R Core Team, 2013).

The fixed-effect variables considered in the analysis are shown in Table 1. They were derived and defined as follows; herd size (number of animals in the herd at the time of breakdown); outbreak reactors (the number of SICCT reactors present in the disclosing test); total reactors (the total number of SICCT reactors during a breakdown); yearly patch prevalence (herd level bTB prevalence for the year); mean patch prevalence (mean herd level bTB prevalence), outward moves year before (the number of outward cattle moves in the year prior to breakdown), and inward moves year before (the number of inward cattle moves in the year prior to breakdown). A categorical herd type variable was included (beef, dairy, other, or unknown). Binary variables were the presence or absence of a milk license, whether lesions consistent with tuberculosis were identified during routine slaughter (LRS), the presence or absence of associated herds (herds are “associated” via e.g., shared management, shared grazing, or shared family responsibilities), and whether the herd had any previous breakdowns during the study period. The herd DVO, the year of breakdown, and the herd unique identifier were included as random effect variables. The distribution of explanatory variables across each DVO in NI is illustrated in Fig. 1.

M. bovis MLVA genotype data

M. bovis MLVA genotype data were derived from isolates obtained from skin-test reactors, and from lesioned animals identified at routine slaughter. These animal-level data were then associated with bTB breakdown-level data. From this, breakdown-level metrics of MLVA genotype richness (number of different MLVA types) were calculated. The process of genotyping M. bovis isolates has been described more fully elsewhere (Kamerbeek et al., 1997; Skuce et al., 2010; Skuce et al., 2005). Briefly, all culture-confirmed bTB cases were sub-cultured to single colonies and heat-killed to create PCR-ready bacterial cell lysates. These were then used as PCR templates for molecular characterisation of pathogen variation. Eight VNTR loci across the M. bovis genome were genotyped; MV2163B/QUB11B, MV4052/QUB26A, MV2461/ETRB, MV1955/Mtub21, MV1895/QUB1895, MV2165/ETRA, MV2163/QUB11A and MV3232/QUB3232 (Durr, Hewinson & Clifton-Hadley, 2000).

Badger density

Badger main sett density was incorporated into models by using a data from the Northern Ireland Badger Survey 2007-08 (Reid et al., 2012). This enumerated and mapped badger main setts within 212 regularly spaced 1 km² squares throughout Northern Ireland, and subsequently spatially interpolated using the Kriging function of the ArcMap 10.5 (ESRI, California, USA) Spatial Analyst toolbox, providing a heat-map proxy of badger density throughout the region (Reid et al., 2012).

Data modelling

During the univariable stage of model fitting for both the count and ordinal models, predictor variables were explored using summary statistics and cross–tabulations with the outcome variable. The relationship between each predictor and the outcome was also visually scrutinised using ggplot2 (Wickham, 2009). Predictor variables were then considered individually for association with the outcome. Correlation coefficients between variables were determined. Variables with moderate or strong correlation ≥ 0.5 or ≤−0.5 were identified, and from these, only those variables with the strongest association with the outcome were retained, based on log-likelihood values. Following univariable assessment, generalised linear mixed models (GLMMs) were fitted. The count model was constructed using the package lme4 (Bates et al., 2015), and the ordinal model was constructed using the package ordinal (Christensen, 2019). Initial modelling of the count data using Poisson regression indicated the presence of over-dispersion (the variance was greater than the mean); a negative binomial model was instead found to be more suitable for these data (Zuur, Hilbe & Ieno, 2015; Zuur et al., 2009).

In both count and ordinal models, the DVO, breakdown year, and herd identifier were included as nested random effects (Zuur et al., 2009). Continuous variables were log-transformed in the final models for computational efficiency, and to improve the model fit (i.e., ensure all explanatory variables were on the same scale, to approximate a more linear relationship, reducing skew and to limit the influence of outliers). All predictors were initially included in the model, including biologically plausible two-way interactions. Final models were assembled using backwards stepwise selection routines; better fitting models were selected on the basis of likelihood ratio tests (Christensen, 2019; Zuur et al., 2009). At each stage, however, model coefficients were manually assessed for confounding (Dohoo, Martin & Stryhn, 2009). Once final models were constructed, excluded predictor variables were again offered to the model and the impact assessed using likelihood ratio tests. Final models were screened for correlations between fixed effects and random effects and were assessed by visual examination of residuals. Plots of residual versus fitted values were firstly explored; residuals were then plotted against all covariates included in the model, and also against the covariates which had been excluded during model fitting. Residuals were used to identify influential data-points, and models were re-run with these data removed for comparative purposes.

Ordinal regression assumes that the effects of explanatory variables are consistent across all outcome categories (i.e., the assumption of proportional odds). We firstly attempted to test this using the nominal_test function of the Ordinal package (Christensen, 2019). However, at the time of analysis, this function was not available for models with multiple random effects. Furthermore, it is presently not feasible to construct an ordered regression model with multiple random effects for which the assumption of proportional odds is also relaxed (Christensen, 2019). To overcome this, we constructed an initial ordered regression model including only fixed effects (via the clm function) and tested the assumption of proportional odds on this model (the nominal_test function). Explanatory variables which violated the assumption of proportional odds were identified and the model was re-ran, wherein the proportional odds assumption was relaxed for these variables. However, the final clmm model was further validated by comparing the model coefficients against those derived from three binary logistic GLMMs (Armstrong & Sloan, 1989; Ananth & Kleinbaum, 1997), with the binary outcome variable dichotomised at the same levels as in the ordinal regression. In these three models, the outcome (breakdown length) was dichotomised as follows: Model 1; ≤180 days (breakdowns 180 days or less classified as 0, all others classified as 1); Model 2 ≤300 days; (breakdowns 300 days or less classified as 0, all others classified as 1); and Model 3 ≤420 days (breakdowns 420 days or less classified as 0, all others classified as 1).

Summary data

The final dataset contained 7,473 breakdowns associated with 5,378 herds. The mean breakdown length was 226 days (SD ± 140 days; approx. seven months) and median breakdown length was 188 days (Inter Quartile Range (IQR): 140–260 days; approx. six months). The longest breakdown was recorded at 2,288 days (6 years). When classified into categories, almost half of all breakdowns (47.18%, n = 3,526) lasted less than 180 days. 34.86% (n = 2,605) of breakdowns were between 181 and 300 days in duration, 11.33% (n = 847) lasted between 301 and 420 days, whilst 6.62% of all bTB breakdowns (n = 495) lasted 421 days or longer (13 months; i.e., 8 or more tests were required to restore OTF status). The distribution of the breakdown length outcome variable is shown in Figs. 2A–2B. Mean breakdown duration varied across NI, from a minimum of 192 days in Derry/Londonderry DVO to a maximum of 266 days in Newry DVO (Fig. 1).

Count model results

The results of the count model of breakdown duration is shown in Table 2 (Table S1). The final model contained seven explanatory variables. The exponentiated results are reported here as Incidence Rate Ratios (IRR) with associated 95% upper and lower confidence intervals (CI). The variables log herd size (IRR: 1.05, 95%CI [1.04–1.06]), log outbreak reactors (IRR: 1.05, 95%CI [1.04–1.06]), log mean patch prevalence (IRR: 1.04, 95%CI [1.01–1.07]) and log MLVA richness (IRR: 1.62, 95%CI [1.58–1.67]) were positively associated with breakdown duration. The binary variables for presence of an LRS (IRR: 1.12, 95%CI [1.09–1.14]), presence of associated herds (IRR: 1.10, 95%CI [1.07–1.13]) and a previous breakdown (IRR: 1.04, 95%CI [1.02–1.07]) were positively associated with breakdown duration. Re-running the model with influential data removed resulted in only minimal change in parameter estimates when compared to the original model (<15% change). The addition of a quadratic term for log MLVA richness was also found to significantly lower log-likelihood; this model is shown in Table S2.

Ordinal model results

Six variables were identified as important predictors in the ordinal model. The parameter estimates of the final model are shown in Table 2 (Table S3). All six variables in the final model were found to be positively associated with the increasing breakdown duration; log herd size (OR: 1.26, 95%CI [1.20–1.32]), log outbreak reactors (OR: 1.34, 95%CI [1.26–1.43]), log mean patch prevalence (OR: 1.20, 95%CI [1.04–1.37]), log MLVA richness (OR: 7.06, 95%CI [6.04–8.24]), the presence of an LRS (OR: 1.79, 95%CI [1.59–2.01]) and the presence of associated herds (OR: 1.49, 95%CI [1.32–1.69]). The coefficients derived from this model were similar to a fixed-effect ordinal regression model, however the variables log herd size, log MLVA richness and log outbreak reactors violated the proportional odds assumption (p < 0.05), suggesting that the effect size is not the same across all three breakdown duration categories. As the assumption of proportional odds was not met for all variables, the coefficients from ordinal model were also compared to those derived from three binomial logistic GLMMs (Fig. 3). There was only limited evidence of the parameter estimates differing between ordinal and binomial models. The binomial model of breakdowns lasting 420 days or less returned a higher odds ratio associated with herd size (OR: 1.49, 95%CI [1.32–1.67]) than the ordinal model (OR: 1.26, 95%CI [1.20–1.32]). The parameter estimate for the number of outbreak reactors was elevated in the binomial model of breakdowns lasting less than 180 days (OR: 1.56, 95%CI [1.45–1.49]) compared to the ordinal model (OR: 1.34, 95%CI [1.26–1.43]), and was also diminished in the model of breakdowns lasting 301 days or more (OR: 1.23, 95%CI [1.03–1.23]) and 421 days or more (OR: 1.08, 95%CI [0.94–1.24]).

MLVA Genotype richness

MLVA genotype richness was the most important variable in both count and ordinal models, in terms of both effect size and decrease in model deviance. This was particularly observable in the ordinal regression model (Fig. 4A). The MLVA genotype richness variable was moderately correlated with the number of inwards moves in the year prior to breakdown (r = 0.33), outwards moves in the year prior to breakdown (r = 0.34) and the number of total reactors over the breakdown (r = 0.39). Further investigation into this “total reactors” variable revealed significantly more reactors in herds with a milk license (mean = 11) than herds without a milk license (mean = 6; Univariable Negative Binomial Regression, IRR: 1.70; 95%CI [1.62–1.79]; Fig. 4B). However, the presence of a milk license was only ‘marginally significant’ in a univariable analysis of breakdown length in both count (IRR: 1.02, 95%CI [1.00–1.05]) and ordinal models (OR: 1.16, 95%CI [1.01–1.27]), and was not retained as a predictor of breakdown length in the finals GLMMs after model building. Further analysis showed that whilst mean breakdown length in herds with a milk license was indeed marginally longer (230 days  ± 141) than in herds without a milk license (224 days ± 140), some of the longest breakdowns were found in herds without milk licenses. For example, there were 27 breakdowns lasting over 1000 days; 10 were in herds with milk licenses, and 17 were in herds without; Fig. 4C. It would therefore appear that whilst production type per-se is not a useful predictor of breakdown length, the results show that some variables which vary between production types—the number of reactors over a breakdown for example (here, confounded with MLVA genotype richness, Fig. 4D), are indeed important predictors of breakdown length.

Badger density results

When modelled using a univariate negative binomial GLM, badger main sett density was a significant predictor of breakdown length (IRR: 1.13, 95%CI [1.13–1.14]). However, this variable was not retained in the final GLMM. Further investigation found that main sett density was correlated with other explanatory variables. Thus, main sett density per-DVO was moderately correlated with breakdown length per-DVO (r = 0.57) and with breakdown length per-patch (r = 0.32), suggesting that the spatial variables already included in the model, notably DVO, captured the general positive relationship observed between main sett density and breakdown length; Fig. 5A. Furthermore, when compared to a fixed effects univariate GLMs where DVO was the sole predictor of breakdown length, the addition of the main sett density variable did not result in a better fitting model (χ² = 0.02, df = 1, p = 0.90). An interaction between DVO and sett density was, however, significant when compared to the fixed-effects model with non-interacting DVO and main sett variables (χ² = 24.24, df = 9, p = 0.004; Tables S4 and S5), suggesting a differential relationship between sett density and breakdown length on a per-DVO basis which was not immediately observable when data were not stratified by DVO. Figures 5C–5D illustrates this observation. Whilst a positive association was found between main sett density and breakdown length in Ballymena, Coleraine, Dungannon, Larne and Derry/Londonderry DVOs, a negative relationship between sett density and breakdown length was observed in Armagh, Enniskillen, Newry, Newtownards and Omagh DVOs (Table S6). To explore this further, we therefore present a second GLMM, (Table 3) in which main sett density was permitted to differ on a per-DVO basis (i.e., a random slopes and random intercepts model). It should be noted, however, that the inclusion of the random slopes term for main sett resulted in only marginally improvements, compared to the original GLMM (Table 2) which did not include a random slope for main sett density per-DVO (χ² = 4.61, df = 2, p = 0.099)

Further analysis also indicated that that the main sett density variable exhibited moderate correlation with mean patch prevalence (r = 0.40; Fig. 3C). To better understand the effect of main sett density on breakdown duration in the absence of spatial confounders, two further alternative models were constructed, both omitting DVO from the random effects component and including log main sett in the fixed effects component. These models also incorporated the other fixed-effect variables reported in Table 2, however, one of these models included patch prevalence in the fixed component, and the other did not. In the model which omitted both DVO and log patch prevalence, log main sett was a significant predictor of breakdown length (OR: 1.08, 95%CI [1.05–1.11]; Table S7). Log main sett was also found to be an important predictor of breakdown length when log patch prevalence was included (OR: 1.08, 95%CI [1.05–1.11]; Table S8), however in this model, log patch prevalence was no longer an important predictor of breakdown duration (OR: 1.01, 95%CI [0.99–1.11]; χ² = 1.02, df = 1, p = 0.31). Confounding between main sett density and DVO was also observed in the ordinal regression. Thus, main sett density was positively associated with increasing breakdown duration categories in a univariable GLM (OR: 1.59, 95%CI [1.42–1.78]), but the main sett variable was not recovered as an important predictor of breakdown length in the mixed model context. We constructed a univariable ordinal GLM with DVO as the sole predictor of breakdown duration category. The coefficients from this model (i.e., the “risk” associated with each DVO) was positively associated with mean sett density per DVO (r = 0.59). Additionally, the inclusion of the main sett variable in this model did not improve model fit (χ² = 0.04, df = 1, p = 0.84).

Wildlife and breakdown duration

In the univariable context, we identified a general positive relationship between breakdown duration and main sett density. Unsurprisingly, in the multivariate context, main sett density was confounded with other spatial and local variables (i.e., DVO and patch), making inferences on the contribution of badger density to infection prolongation less clear. It is not yet possible to conclusively distinguish between local sources of infection (which may include wildlife, contiguous herds and environmental contamination), but given our data, we cannot reject the hypothesis that badgers may be involved in the maintenance of local patch bTB prevalence via spillback infection to cattle. Whether infected badger presence has a greater risk of sporadic introduction of infection into herds (singular badger-cattle spillover), than longer-term maintenance within herds (explosive introduction of infection elevated with cattle-cattle transmission), remains to be determined. Whilst this study was unable to conclusively clarify the relationship between badger density and breakdown duration, our data nevertheless reveal important features that warrant further investigation in future studies.

Thus, despite DVO capturing the risk associated with main sett density and the general positive association between breakdown length and main sett density, there was some evidence of within-DVO effects. Within five DVO areas (Ballymena, Coleraine, Dungannon, Larne and Derry/Londonderry), increasing sett density was generally associated with longer breakdowns. In the other five DVO areas (Armagh, Enniskillen, Newry, Newtownards and Omagh), increasing sett density was generally associated with shorter breakdowns. The five DVOs with a positive association between sett density and breakdown duration were areas of generally lower badger sett densities (Reid et al., 2012). The DVOs with a negative association between sett density and breakdown duration were generally associated with higher badger sett densities. Whilst the interpretation of this is not straightforward, differences in farming practice (e.g., farm fragmentation) or differences in badger ecology (e.g., population context dependent badger dispersal; (Byrne et al., 2019) across the region could partially explain this observation. However, this does not preclude that the relationship between herd bTB risk from badgers may not be simply be dependent on wildlife density; the sett density data provides no insight regarding disease prevalence within the badger population. A spatially explicit model of disease prevalence in badgers may resolve this in future. Indeed, future research could investigate variation in wildlife TB transmission risk (LaHue et al., 2016) as a function of infection prevalence as well as density, and investigate how that could help to partially explain patterns within cattle data.

Conversely, infection risk in cattle may not be linked to badger disease prevalence or population density, but may instead related to the relative frequency of interactions between infected badgers and susceptible cattle (Böhm, Hutchings & White, 2009). Alternatively, it may be that indirect transmission of bTB via, for example, cattle accessing badger latrines, is more critically associated with chronic bTB breakdowns as opposed to wildlife population density per se (Campbell et al., 2019; Drewe et al., 2013). Furthermore, despite sett density being a convenient metric, we must be careful when inferring the relationship between sett density and population density, as the magnitude of the association can change depending on the local dynamics. For example, population density can increase without necessarily increasing the number of setts, via an increase in the mean group size (Judge et al., 2014). Alternatively, where badger population densities are depressed (e.g., though hunting, culling, or illegal disturbance), sett density metrics can overestimate true local density. In Ireland, sett density was found to be good predictor of increased herd breakdown risk early in a six-year study, but progressively became a weaker predictor as a program of targeted badger culling reduced population density (Byrne et al., 2014b). Therefore, investigating intricate relationships between wildlife and domestic hosts may well require even more detailed information around population abundance at large scales, in the Northern Ireland context this could include mark-recapture and/or the use of remote camera trapping technologies (Campbell et al., 2019).