Research article | Open | Published:
Variation in Mycobacterium bovis genetic richness suggests that inwards cattle movements are a more important source of infection in beef herds than in dairy herds
BMC Microbiologyvolume 19, Article number: 154 (2019)
We used genetic Multi-Locus VNTR Analysis (MLVA) data gathered from surveillance efforts to better understand the ongoing bovine tuberculosis (bTB) epidemic in Northern Irish cattle herds. We modelled the factors associated with Mycobacterium bovis MLVA genotype richness at three analytical scales; breakdown level, herd level, and patch level, and compared the results between dairy and non-dairy production types.
In 83% of breakdowns and in 63% of herds, a single MLVA genotype was isolated. Five or more MLVA genotypes were found in less than 3 % of herds. Herd size and the total number of reactors were important explanatory variables, suggesting that increasing MLVA genotype richness was positively related to increases in the number of host animals. Despite their smaller relative size, however, the highest MLVA genotype richness values were observed in non-dairy herds. Increasing inwards cattle movements were important positive predictors of MLVA genotype richness, but mainly in non-dairy settings.
The principal finding is that low MLVA genotype richness indicates that small-scale epidemics, e.g. wildlife, contiguous farms, and within-herd recrudescence, are important routes of M. bovis infection in cattle herds. We hypothesise that these mechanisms will maintain, but may not explicitly increase, MLVA genotype richness. The presence of elevated MLVA richness is relatively rare and likely indicates beef fattening enterprises, which purchase cattle from over long distances. Cattle movements were furthermore an important predictor of MLVA genotype richness in non-dairy herds, but not in dairy herds; this may represent reduced cattle purchasing levels in dairy enterprises, compared to beef. These observations allude to the relative contribution of different routes of bTB infection between production types; we posit that infection associated with local factors may be more evident in dairy herds than beef herds, however in beef herds, inwards movements offer additional opportunities for introducing M. bovis into the herd.
Genetic approaches such as spoligotyping, MLVA genotyping, and more recently, whole genome sequencing (WGS) are routinely used to help answer epidemiological questions [1,2,3,4]. The molecular characterisation of infectious agents can be used to better understand disease outbreak, spread and maintenance, and to better inform interventions [5,6,7]. In Northern Ireland (NI), MLVA genotyping [1, 8] is a molecular epidemiological tool deployed as a surveillence measure to investigate the ongoing bovine tuberculosis (bTB) epidemic. Mycobacterium bovis is the causative agent of bTB, a serious bacterial disease of cattle and livestock. The animal and human health implications of bTB infection are considerable, and far-reaching economic costs are associated with disease control. As yet, however, eradication has remained elusive throughout much of the United Kingdom (UK) and Republic of Ireland (ROI) [9, 10].
MLVA genotyping is carried out on M. bovis isolates from culture-confirmed bTB cases, including reactor cattle which test positive to the Single Intradermal Comparative Cervical Tuberculin (SICCT) test, and lesioned animals identified during meat inspection . The information gained from these surveillance activities is used to undertake epidemiological investigations. For example, different M. bovis genotypes are associated with relatively different proportions of lesioned reactors in cattle  and increasing within-herd M. bovis genotype richness (i.e. numbers of different MLVA genotypes) is linked to prolonged and recurrent bTB breakdowns .
The M. bovis MLVA genotype data also reveals that the genetic structure of the M. bovis population in NI is spatially aggregated [14, 15] a phenomenon that has also been observed in M. bovis populations in cattle elsewhere [16,17,18,19]. Furthermore, M. bovis MLVA genotypes are co-localised between cattle hosts and a known wildlife reservoir, the European badger (Meles meles) . This spatial clustering suggests that bTB infection may be associated with small scale, local epidemics which maintain (but do not necessarily increase) M. bovis MLVA genotype richness. Such transmission routes may include spillback from infected wildlife in the area [4, 20], recrudescence of M. bovis MLVA genotypes within-herd , infection from a contaminated environment , or from contiguous herds [23, 24]. Notwithstanding these processes, increases in M. bovis MLVA genotype richness have also been observed within herds , alluding to other potential routes of infection. Long-range cattle movements [25, 26], or spreading of contaminated slurry sourced from different areas [27, 28] also present sources of infection, potentially associated with the introduction of additional MLVA genotypes. It is possible that on rare occasions, new M. bovis MLVA genotypes may also be introduced via infected wildlife ; badgers occasionally undertake long-range movements, which may facilitate wider dissemination of infection . However, the different routes of infection are not mutually exclusive, and their relative contribution to the infection burden is unknown in the context of NI. In GB, however, it is estimated that 16% of herd infections arise from cattle movements, with 75% of infection arising from “local effects”, including badgers and contiguous spread .
Different mechanisms of infection may furthermore be associated with different herd management practices. In GB, there is evidence that beef fattening farms may be at more risk of acquiring infection via inwards cattle purchases, whereas dairy farms are more likely to spread infection via the onwards sale of animals . In NI, the purchase of beef cattle is associated with increased risk of bTB breakdown , and thus in beef settings, inwards cattle movements may be a relatively more important source of infection than in dairy settings. Furthermore, Lahuerta-Marin et al.,  found evidence to suggest that in chronically infected herd in NI, the performance characteristics of the skin test and the interferon-gamma test may be poorer in dairy herds relative to beef herds . In GB and NI, dairy reactor animals were less likely to have confirmed infection at post-mortem than non-dairy animals [33, 34]. Dairy herds may therefore be at increased risk of persistent and recurrent breakdowns [13, 35, 36], a phenomenon associated with the problem of fully clearing infection using ante-mortem diagnostics . Together, this uggests a reduced ability to eradicate infection in dairy settings compared to beef settings and thus, within-herd recrudescence may be a more evident route of infection in dairy herds. As yet, however, the relative contribution of different mechanisms of bTB infection in different herd types in NI is not well understood. What follows, therefore, is the quantification of M. bovis MLVA genotype richness across multiple analytical scales in NI, in both dairy and non-dairy settings. We determine the factors associated with this observed richness, and use the results to better understand infection dynamics across both production types.
Final datasets were derived from 29,473 animal-level tests, associated with 7478 bTB breakdowns from 5378 herds. There was moderate correlation between breakdown level MLVA richness and herd level MLVA richness (Spearman’s Rank (Rs) = 0.52, p < 0.05) and every unit increase in breakdown level MLVA richness was associated with an increase in herd level MLVA richness of 38% (IRR:1.38, 95%CI: 1.36–1.40, Fig. 1a). Only weak correlation (Rs = 0.08, p < 0.05) was identified between individual herd level richness and patch level MLVA richness, with a 2.1% increase in patch level richness for a unit increases in herd level richness (Univariable Poisson Regression, IRR: 1.02, 95%CI: 1.01–1.03, Fig. 1b). Instead, the mean herd level MLVA richness value exhibited a stronger association with patch level metrics of richness (Rs = 0.40, p < 0.05) and a greater effect size (Univariable Poisson Regression, IRR: 1.35, 95%CI: 1.22–1.49). Measures of MLVA genotype richness at the breakdown, herd and patch level varied spatially across NI. Fig. 1c illustrates the broad variation in patch level MLVA genotypic richness in each DVO. Armagh DVO was associated with the highest mean breakdown level and herd level MLVA genotypic richness. However, the highest patch level MLVA genotypic richness was observed in Newtownards DVO. Londonderry DVO was associated with the lowest breakdown level MLVA genotypic richness, herd-level genotypic richness, and patch-level genotypic richness (Table 2).
Exploratory analysis revealed differences in the distribution of explanatory variables between herds with milk licences, and herds without. Whilst the mean herd size at time of breakdown was 141 animals (Standard Deviation (SD) ±146), herds with milk licences (242 animals, SD ± 171) were significantly larger than beef herds (95 animals, SD ± 104; Univariable Poisson Regression, IRR: 2.55, 95%CI: 2.54–2.56). In 40% of breakdowns (n = 2964), there were two or fewer reactors disclosed; the majority of these breakdowns (78%, n = 2300) occurred in herds without a milk licence. However, the total number of reactors variable did not exhibit the same distribution between dairy and non-dairy herds; increases in the number of reactors disclosed during a breakdown had a positive association with the presence of a milk licence (Univariable Poisson Regression, IRR: 1.7, 95%CI: 1.68–1.73). The presence of a milk licence was also significantly negatively associated with inwards movement intensity (Univariable Poisson Regression, IRR: 0.21, 95%CI: 0.20–0.22).
MLVA genotype richness at the breakdown level
In 82.5% of all breakdowns (n = 6168), only one MLVA type was isolated. The maximum MLVA genotype richness observed was 12 (n = 1), and 49 breakdowns (0.66%) yielded 5 or more different MLVA types. Whilst the presence of a milk licence was a positive predictor of MLVA richness, it was not statistically significant (Univariable Poisson Regression: IRR: 1.01, 95%CI: 0.96–1.05). Furthermore, non-dairy herds exhibited greater accumulation of MLVA genotypes. For example, the maximum number of different MLVA genotypes isolated from a breakdown in a non-dairy herd was 12, with 41 breakdowns yielding five or more different MLVA genotypes. In breakdowns in dairy herds however, the maximum number of MLVA genotypes isolated was seven, and only eight breakdowns yielded five or more different MLVA genotypes (Fig. 2a).
All models of breakdown level MLVA genotype richness performed better than “null models” containing the variable for the total number of reactors disclosed during a breakdown (Additional file 1: Table S3). Five main effect explanatory variables were included in the final model for MLVA genotype richness at the breakdown level (Table 3 and Additional file 1: Table S4). The total number of reactors over the breakdown, herd size at the time of breakdown, inwards movement intensity, and breakdown length in days were positively associated with MLVA genotype richness. Latitude was negatively associated with MLVA genotype richness. An interaction term between herd size and inwards movement intensity showed that MLVA genotype richness was elevated in breakdowns in larger herds with large numbers of cattle purchases.
When bTB breakdown data from herds with a milk licence were analysed independently, three explanatory variables were included in the final model (Table 3 and Additional file 1: Table S5); the total number of bTB reactors over the breakdown, herd size at the time of breakdown, and breakdown length in days. The model constructed using data from breakdowns in non-dairy herds (Table 3 and Additional file 1: Table S6) contained four fixed effect explanatory variables and one interaction term. The total number of reactors over the breakdown, herd size at the time of breakdown, inwards movement intensity, and the breakdown length in days were positively associated with MLVA genotype richness. An interaction between herd size and inwards movement intensity showed that larger breakdown herds, with more inwards cattle movements before the bTB breakdown, were furthermore associated with elevated MLVA genotype richness.
MLVA genotype richness at the herd level
After data were aggregated to herd level, no increases in MLVA genotype richness were observed in 70.6% of herds (n = 3795). However only one breakdown occurred in the majority of herds in the dataset (70.9%; n = 3803). In those herds with more than one breakdown (n = 1575), no increase in MLVA genotype richness was observed in 26.3% of cases (n = 414). In 57 herds (3.62%) which experienced more than one breakdown, the number of MLVA genotypes in the breakdown level analysis and the herd level analysis increased by five or more.
Single MLVA isolates were observed in 62.9% of herds (n = 3382) and the maximum number of MLVA genotypes isolated was 19 (n = 1), with 159 herds (2.96%) yielding five or more MLVA genotypes. In this analysis, the presence of a milk licence also exhibited a slight association with increasing MLVA genotype richness (Univariable Poisson Regression, IRR: 1.06, 95%CI: 1.02–1.1). Despite this, there was greater accumulation in the overall number of MLVA genotypes in non-dairy herds, as 133 non-dairy herds (3.5%) were associated with between five and 19 different MLVA genotypes, compared to the 26 dairy herds (1.7%) that were associated with between five and 10 MLVA genotypes (Fig. 2b).
The final model for MLVA genotype richness at the herd level for all herds in the dataset contained six main effect variables and two interactions (Table 3 and Additional file 1: Table S7). Thus, the total number of reactors over the breakdown, herd size at the time of breakdown, the count of bTB breakdowns, inwards movement intensity, and the mean breakdown length in days were all positively associated with MLVA genotype richness at the herd level. The number of animal-level M. bovis isolates of the most common MLVA genotype found in the herd was instead negatively associated with MLVA genotype richness. A positive interaction between mean herd size and inwards movement intensity indicated that MLVA genotype richness is amplified in larger herds which are also engaged with high levels of inwards cattle purchases. A negative interaction between breakdown length and herd size suggested that larger herds which experienced longer breakdowns were associated with fewer M. bovis genotypes than smaller herds which experienced longer breakdowns.
The final model of MLVA genotype richness for herds with a milk licence consisted of six main effect variables and one interaction (Table 3 & Additional file 1: Table S8). Here, the total number of reactors over the breakdown, herd size at the time of breakdown, the count of bTB breakdowns, inwards movement intensity, and the mean breakdown length in days were positively associated with MLVA genotype richness. As with the full model, the number of animal-level M. bovis isolates of the most common MLVA genotype was negatively associated with MLVA genotype richness in these herds. Larger herds with higher levels of inwards movement intensity were also associated with elevated herd-level MLVA genotype richness.
In non-dairy herds (Table 3 and Additional file 1: Table S9), MLVA genotype richness was positively associated with the total number of reactors over the breakdown, herd size at the time of breakdown, the count of bTB breakdowns, inwards movement intensity, and the mean breakdown length in days. Two interactions were identified; one interaction showed that larger herds with higher levels of inwards movement intensity were associated with higher herd-level MLVA genotype richness. The second interaction suggested that longer breakdowns in larger herds were associated with less MLVA genotype richness than in smaller herds.
MLVA genotype richness at the patch level
There were no patches identified with single isolates, however 23.6% of patches (n = 29) were associated with 10 or fewer different MLVA genotype isolates, and 82.1% of patches (n = 101) were associated with 20 or fewer different MLVA genotype isolates. Comparing patch level MLVA richness between dairy and non-dairy herds (Fig. 3a-c) revealed a lower median number of isolates when assessing richness in the dairy herds only (median = 8, IQR: 5–12) compared to the non-dairy herds (median = 14, IQR: 10–18). Patch level richness in dairy herds was also found to be significantly lower than patch level richness in non-dairy herds (Univariable Poisson Regression: IRR: 0.58, 95%CI: 0.52–0.64).
The explanatory model for patch-level MLVA genotype richness contained three fixed effect variables (Table 3 & Additional file 1: Table S10); the number of reactors in the patch, total inwards movement intensity in the patch, and mean herd level MLVA genotype richness. When patch-level MLVA genotype richness was derived using only data from those herds with milk licences, the final model consisted of three explanatory variables (Table 3 & Additional file 1: Table S11); the number of reactors in the patch, the total number of cattle on the patch, and the mean herd-level MLVA genotype richness value of the patch. Metrics of patch-level MLVA genotype richness from herds without milk-licences were also associated with three variables (Table 3 and Additional file 1: Table S12); the number of reactors in the patch, total patch-level inwards movement intensity, and the mean herd-level MLVA genotype richness value of the patch.
The results reveal that a single MLVA genotype was isolated in the majority of cases (i.e. 83% of breakdowns, and 63% of herds). The accumulation of MLVA genotype richness at the breakdown level and at the herd level was less common; less than 1 % of breakdowns and 3 % of herds were associated with five or more MLVA genotypes, respectively. Our principal finding is, therefore, that M. bovis infection is largely related to a single source of infection, or to multiple local sources of infection associated with the same MLVA genotype. These may include nearby infected wildlife, infected cattle on neighbouring contiguous farms , contamination from the environment, or recrudescence of infection associated with MLVA genotypes circulating within the herd [13, 21, 35, 36, 47]. In GB, the majority of M. bovis infection was also related to “local” factors such as wildlife and contiguous spread ; our results would appear to be in agreement with these findings. We also found that the number of breakdowns was positively associated with MLVA genotype richness, indicating that in some cases, M. bovis infection does appear to be associated with the introduction of new MLVA genotypes. However, the vast majority of herds in our study (71%) experienced only one breakdown during the study period. It would therefore appear that the occurrence of breakdowns associated with newly introduced MLVA genotypes is less common. Furthermore, in 26% of herds which did experience more than one breakdown, there were no increases in MLVA genotype richness over successive restriction periods. These results further highlight the relative importance of small-scale, local effects in maintaining the epidemic.
The M. bovis population in NI is spatially clustered, with MVLA genotypes localised to geographic areas . Consequently, the accumulation of MLVA genotype richness in herds is indicative of purchasing cattle from a wide geographic extent , or possibly the spreading of non-local contaminated slurry . However, elevated MLVA richness was observed in only a limited number of cases, suggesting that a process maintains richness in some herds without the wider dissemination of M. bovis genotypes. We posit that elevated MLVA genotype richness within-herd implies that these herds are beef fattening enterprises. These herds operate via the purchase of a large number of cattle, which are then raised for sale, straight to slaughter. This affords the opportunity for the accumulation of MLVA genotype richness with less risk of spreading infection to other herds via the onwards sale of cattle. In GB, beef fattening herds are indeed shown to have a higher turnover of animals than dairy herds  which may facilitate the introduction of novel MLVA genotypes. In beef fattening enterprises, the purchase of infected cattle therefore appears to serve as a particularly important source of M. bovis infection, in addition to wildlife, nearby farms, within herd recrudescence and environmental contamination.
The low MLVA genotype richness observed in most herds does not mean that cattle movements do not contribute more widely to the epidemic; cattle movements are a known risk factor for bTB in NI . Indeed, a main finding from this work is that higher levels of inwards movement intensity were associated with increasing MLVA genotype richness in non-dairy herds. The impact of inwards cattle movements on MLVA genotype richness in dairy herds was less evident. Previous work from GB revealed differences between purchasing patterns in beef and dairy settings, with dairy herds receiving fewer inwards movements than beef herds ; our findings would therefore appear to be in accord with these results. We argue that cattle movements may increase the risk of bTB breakdown, but may not always be important for introducing new MLVA genotypes in dairy herds. Thus, if most cattle purchases consist of short range movements (as opposed to buying and selling cattle from further afield), the same MLVA genotypes will circulate between herds, therefore maintaining current herd level MLVA genotype richness (but not increasing it). We further hypothesise that the spatially structured nature of the M. bovis population  means that disentangling the relative contributions of local infection routes is exceptionally challenging; a higher resolution phylodynamics approach such as whole genome sequencing may shed more light on this issue .
Whilst there was a moderate correlation between breakdown level and patch level metrics of MLVA genotype richness, individual herd level MLVA richness values exhibited only weak correlation with the overall patch level MLVA richness. Higher patch level MLVA richness will indeed partially reflect the larger host population size in patches compared to herds, but it also further shows that herds within a patch share some, but not all, of the same MLVA genotypes. This further hints at the importance of localised infection mechanisms in maintaining MLVA richness within herds. In addition, we found that the mean herd level MLVA richness in the patch did indeed exhibit a stronger association with patch level MLVA richness. This is likely due to the influence of herds with elevated MLVA type richness (i.e. beef fattening herds) impacting on the mean.
In all models, the total number of reactors disclosed during a breakdown was an important explanatory variable for MLVA genotype richness. Thus, at least in part, the accumulation of different M. bovis MLVA genotypes is likely to be related to increases in the number of hosts. Increases in herd size were also positively associated with MLVA genotype richness; this may be a further manifestation of the positive relationship between MLVA genotype richness and increasing cattle numbers. However, large herds may also be linked to other factors which afford more opportunities for the accumulation of new M. bovis MLVA genotypes e.g. inwards cattle movements, or a larger geographic footprint.
The results furthermore demonstrate that dairy herds were both larger, and contained more reactors, than herds without a milk licence. There was also evidence that at the herd level, the presence of a milk licence had a slight positive association with increased MVLA genotype richness. In spite of this, the highest MLVA genotype richness values were not found in dairy settings. Previous studies show that dairy herds in NI are associated with elevated risk of bTB breakdown compared to other herd types  and the diagnostic SICCT can perform relatively poorly in NI dairy settings . Together, this suggests a failure to clear infected animals within dairy herds, thus increasing the likelihood of within herd recrudesce of the same MLVA genotypes. Other factors associated with dairy herds may also limit the opportunity for MLVA genotype accumulation, for example, if dairy farms are less expansive, less fragmented, utilise less shared grazing, or share fewer boundaries with neighbouring farms. As yet however, how farm fragmentation factors and grazing practices influence bTB risk across different herd types is presently unknown in the NI context.
The herd-level analysis suggested that increasing numbers of the most common MLVA genotype found in dairy herds was negatively associated with MLVA genotype richness at the herd level. We consider that rapid propagation of a single MLVA genotype within-herd could arise if within-herd transmission dynamics were different between beef and dairy herds . In dairy settings, factors such as differences in stress profiles , higher stocking density, or increased access to shared housing or feed [24, 52], may provide opportunity for the rapid amplification of within-herd bTB infection. However, rates of within-herd transmission, and the factors influencing this metric in cattle herds in NI are not yet known.
This work alludes to the fact that different routes of infection may be of greater relative importance to different production types in NI. In the majority of herds, there was relatively low MLVA genotype richness, indicating that highly localised factors operating over a small geographical extent were primarily influencing the epidemic, e.g. infection from neighbouring herds, wildlife or within herd recrudescence. Only a limited number of herds were identified with elevated MLVA richness, and these were likely beef fattening enterprises, wherein long range cattle purchases may also represent an important source of new MLVA genotypes. Indeed, inwards movement intensity was an important predictor of increasing MLVA genotype richness in non-dairy settings. Cattle movements appear to be relatively less important for introducing MLVA genotypes into dairy herds, which may be a consequence of less movement intensity associated with dairy herds more generally, and/or the purchase of cattle from the local area which are more likely to share the same MLVA genotypes. Despite the larger number of potential M. bovis hosts, the maximum MLVA genotype accumulation was not observed in dairy herds. We therefore posit that in dairy herds, localised infection routes which do not explicitly introduce MLVA genotypes from outside the area, e.g. within herd recrudescence, are more evident. Non dairy herds, however, may also be exposed to additional infection via cattle purchasing. Future eradication should focus on understanding the contribution of infection pathways in different production types. An enhanced understanding of the NI cattle movement network will provide a better understanding of the risk associated with cattle purchasing in different settings. In addition, we advocate that future work should focus on the infection risk associated with farm fragmentation, contiguous herds and shared grazing between different herd types in the Northern Irish context.
This retrospective analysis used data collected between 2009 and 2016 inclusive, and was conducted at three analytical scales; breakdown level, herd level and patch level. Cattle herds in NI reside within one of 123 administrative patches, which are situated within 10 Divisional Veterinary Offices (DVOs; Additional file 1: Figure S1). The area of NI is approximately 14,000 km2 and mean patch size is approximately 110 km2 (SD ± 53). Data on confirmed breakdowns which occurred between January 2003 to April 2016 were made available from the Animal and Public Health Information System (APHIS) database . BTB breakdowns with missing or erroneous data (e.g. location co-ordinates, start and end dates) were removed. Breakdowns were only included if the start and end date (defined by the date of disclosing test, and the date of at which Officially Tuberculosis Free status was restored) was inclusive of the study period. These data were then associated with animal-level M. bovis genotype surveillance data (MLVA) from skin-test reactors and lesioned non-reactor animals identified at slaughter . Briefly, all culture-confirmed bTB cases were further sub-cultured with single colonies, then heat killed to create bacterial cell lysates, which were then used directly as Polymerase Chain Reaction (PCR) templates for all molecular characterisation of pathogen variation. Spoligotyping and genotyping of variable number of tandem repeat (VNTR) loci were undertaken using established methods [14, 39]. Standardised international nomenclature defined at https://www.mbovis.org/ was applied for all spoligotypes in the dataset. 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. Prior to 2009, MLVA and spoligotyping were usually only carried out on the first confirmed reactor or lesioned case, whereas post-2009, molecular typing was carried out at animal-level, on all culture-confirmed cases.
The outcome variable in each analysis was MLVA genotype richness, defined as the number of different MLVA genotypes present at each of the three analytical scales. In this study, reactor cattle were treated as individuals, and different M. bovis MLVA genotypes were treated as species. At all three scales, three different models were constructed to reflect different enterprise types. Firstly, data on all herds were analysed together. Next, data from herds with a milk-licence (i.e. dairy units) were analysed separately from herds without a milk-licence (non-dairy herds; e.g. beef herds). Nine models were therefore constructed in total. Seven variables were considered for inclusion in the breakdown level analysis (Table 1). These variables were; herd size at the time of breakdown, bTB breakdown length in days, the total number of SICCT reactor animals identified during the breakdown, and herd latitude and longitude in the form of the Irish six figure grid reference. Inwards movement intensity was defined as the number of animals moved into the herd in the six months prior to breakdown, as a proportion of herd size. The number of isolates from animals with the most common MLVA genotype from the breakdown was also derived (i.e. counts of the most common MLVA genotype). For the herd level analysis, data on bTB breakdowns were aggregated to the herd unit, resulting in ten herd level variables (Additional file 1: Table S1). In brief, the outcome variable was calculated as the total number of different MLVA genotypes present in the herd, throughout all breakdowns. The total number of reactors over the breakdown was summed to give the total number of reactors in the herd. The total number of breakdowns per-herd, and the number of isolates from reactors associated with the most common MLVA genotype were also derived. The mean values for herd size, inwards movement intensity, and breakdown length were also quantified. The maximum inward movement intensity, and breakdown length, were also calculated and included as potential explanatory variables. The herd latitude and longitude were as before. For the patch level analysis, breakdown level data were aggregated to the patch unit; 15 variables were included in this analysis (Additional file 1: Table S2). The patch-level MLVA richness variable was derived from the number of unique MLVA genotypes present in the patch. The total number of herds, the total number of herds with a milk licence, the total number of cattle in breakdown herds, the total number of reactors, the total number of breakdowns, the total inwards movement intensity (total number of inwards movements into the patch, as a proportion of cattle in breakdown herds), and the total number of breakdown days per-patch were calculated. The total number of isolates from reactors associated with the most common MLVA genotype in the patch was also quantified. The mean number of breakdowns per-herd, the mean herd-size, the mean inwards movement intensity, and the mean breakdown length per-herd were also included. The latitude and longitude of each patch centroid was calculated. All analyses were carried out in R  and Excel.
Models were constructed using the lme4 package . The outcome (counts of MLVA genotypes) was modelled using Poisson regression with a log link. The same modelling procedure was applied to each of the nine models. Firstly, univariable relationships between the outcome and each of the predictor variables were explored. Continuous predictor variables were assessed for correlation as a potential source of multi-colinearity. Variables exhibiting a correlation coefficient > 0.5 or < − 0.5 were considered for removal, whilst retaining predictors with the strongest association with the outcome, as determined by Akaike Information Criterion (AIC) values. The remaining covariates (including biologically plausible interactions) were considered for inclusion in final models. As these data were hierarchical, the most appropriate random effects structure in each model was firstly derived. For the breakdown level analysis, random effects included the herd, the patch and the DVO. For the herd level analysis, potential random effects included the patch and DVO. For the patch level analysis, DVO was the only random effect considered.
Following the procedures outlined in Zurr et al. (2009) , initial models were constructed with all potential covariates included in the fixed effects component, along with the appropriate random effects for the given model. Models were also constructed which included fixed effects only. The decision to use random effects was informed by comparing AIC values between models with and without random effects, using both AIC values and by assessing the variance associated with each random effect. Once the random effects structure was determined, AIC values were then used as a guide for model selection following a backwards stepwise process, with smaller AIC values indicating better fitting, more parsimonious models. During model construction, variables were also assessed for confounding effects. Finally, any variables omitted from the modelling process were reintroduced and the impact on AIC values was determined, in order to ensure that all important predictors were included. Models were checked to ensure that over-dispersion was not present. After model construction, models were validating by plotting covariates against residuals to assess for patterns. Residuals were also plotted against covariates which were not included in the model. Finally, models were re-run, excluding data associated with large residual values, and the impact on covariates was assessed. Spatial autocorrelation in the residuals was checked using the GeoR package  and Moran’s I , however no significant evidence of spatial autocorrelation was detected. All predictor variables were scaled and centred at the mean prior to analysis, to permit comparisons of the relative effects of covariates. The exponentiated parameters of the resulting models can be interpreted as Incidence Rate Ratios (IRR), which are reported along with the 95% upper and lower confidence intervals (CI). Final models were also compared to “null” models which contained only the number of reactors variable as the sole predictor (Additional file 1: Table S3).
Availability of data and materials
The datasets generated and analysed during the current study are not publicly available by request of the data controller (the project funder) who deem these data regarding disease status as sensitive. Any data related enquires should be made to the Department of Agriculture, Environment and Rural Affairs (DAERA; https://www.daera-ni.gov.uk).
Akaike Information Criterion
Animal and Public Health Information System
Divisional Veterinary Offices
Incidence Rate Ratio
Multi-Locus VNTR Analysis
Polymerase Chain Reaction
Republic of Ireland
- Rs :
Spearman’s Rank Correlation Coefficient
- SICCT test:
Single Intradermal Comparative Cervical Tuberculin test
Variable Number of Tandem Repeat
Whole Genome Sequencing
Roring S, Scott A, Brittain D, Walker I, Hewinson G, Neill S, Skuce R. Development of variable-number tandem repeat typing of Mycobacterium bovis: comparison of results with those obtained by using existing exact tandem repeats and Spoligotyping. J Clin Microbiol. 2002;40:2126–33. https://doi.org/10.1128/JCM.40.6.2126-2133.2002.
Skuce RA, McDowell SW, Mallon TR, Luke B, Breadon EL, Lagan PL, McCormick CM, McBride SH, Pollock JM. Discrimination of isolates of Mycobacterium bovis in Northern Ireland on the basis of variable numbers of tandem repeats (VNTRs). Vet Rec. 2005;157:501–4. https://doi.org/10.1136/vr.157.17.501.
Sola C, Filliol I, Legrand E, Lesjean S, Locht C, Supply P, Rastogi N. Genotyping of the Mycobacterium tuberculosis complex using MIRUs: association with VNTR and spoligotyping for molecular epidemiology and evolutionary genetics. Infect Genet Evol. 2003;3:125–33. https://doi.org/10.1016/S1567-1348(03)00011-X.
Biek R, O'Hare A, Wright D, Mallon T, McCormick C, Orton RJ, McDowell S, Trewby H, Skuce RA, Kao RR. Whole genome sequencing reveals local transmission patterns of Mycobacterium bovis in sympatric cattle and badger populations. PLoS Pathog. 2012;8:e1003008. https://doi.org/10.1371/journal.ppat.1003008.
van Belkum A, Struelens M, de Visser A, Verbrugh H, Tibayrenc M. Role of genomic typing in taxonomy, evolutionary genetics, and microbial epidemiology. Clin Microbiol Rev. 2001;14:547–60. https://doi.org/10.1128/cmr.14.3.547-560.2001.
Urwin R, Maiden MCJ. Multi-locus sequence typing: a tool for global epidemiology. Trends Microbiol. 2003;11:479–87. https://doi.org/10.1016/j.tim.2003.08.006.
Smith NH, Dale J, Inwald J, Palmer S, Gordon SV, Hewinson RG, Smith JM. The population structure of Mycobacterium bovis in Great Britain: clonal expansion. Proc Natl Acad Sci U S A. 2003;100:15271–5. https://doi.org/10.1073/pnas.2036554100.
Skuce RA, McCorry TP, McCarroll JF, Roring SMM, Scott AN, Brittain D, Hughes SL, Hewinson RG, Neill SD. Discrimination of Mycobacterium tuberculosis complex bacteria using novel VNTR-PCR targets. Microbiology. 2002;148:519–28. https://doi.org/10.1099/00221287-148-2-519.
Waters WR, Palmer MV, Buddle BM, Vordermeier HM. Bovine tuberculosis vaccine research: historical perspectives and recent advances. Vaccine. 2012;30:2611–22. https://doi.org/10.1016/j.vaccine.2012.02.018.
Allen AR, Skuce RA, Byrne AW. Bovine tuberculosis in Britain and Ireland – a perfect storm? The confluence of potential ecological and epidemiological impediments to controlling a chronic infectious disease. Front Vet Sci. 2018;5. https://doi.org/10.3389/fvets.2018.00109.
Abernethy DA, Denny GO, Menzies FD, McGuckian P, Honhold N, Roberts AR. The Northern Ireland programme for the control and eradication of Mycobacterium bovis. Vet Microbiol. 2006;112:231–7. https://doi.org/10.1016/j.vetmic.2005.11.023.
Wright DM, Allen AR, Mallon TR, McDowell SWJ, Bishop SC, Glass EJ, Bermingham ML, Woolliams JA, Skuce RA. Field-isolated genotypes of Mycobacterium bovis vary in virulence and influence case pathology but do not affect outbreak size. PLoS One. 2013;8:e74503. https://doi.org/10.1371/journal.pone.0074503.
Milne MG, Graham J, Allen AR, Lahuerta-Marin A, McCormick C, Presho E, Skuce RA, Byrne AW (2018) Herd characteristics, wildlife risk and bacterial strain genotypes in persistent breakdowns of bovine tuberculosis in northern Irish cattle herds. Conference: Society for Veterinary Epidemiology and Preventive Medicine.
Skuce RA, Mallon TR, McCormick CM, McBride SH, Clarke G, Thompson A, Couzens C, Gordon AW, McDowell SWJ. Mycobacterium bovis genotypes in Northern Ireland: herd-level surveillance (2003 to 2008). Vet Rec. 2010;167:684–9. https://doi.org/10.1136/vr.c5108.
Trewby H (2016) The genetic and spatial epidemiology of bovine tuberculosis in the UK : from molecular typing to bacterial whole genome sequencing University of Glasgow.
Woodroffe R, Donnelly CA, Johnston WT, Bourne FJ, Cheeseman CL, Clifton-Hadley RS, Cox DR, Gettinby G, Hewinson RG, Le Fevre AM, McInerney JP, Morrison WI. Spatial association of Mycobacterium bovis infection in cattle and badgers Meles meles. J Appl Ecol. 2005;42:852–62. https://doi.org/10.1111/j.1365-2664.2005.01081.x.
Furphy C, Costello E, Murphy D, Corner LAL, Gormley E. DNA typing of Mycobacterium bovis isolates from badgers (Meles meles) culled from areas in Ireland with different levels of tuberculosis prevalence. Vet Med Int. 2012.
Goodchild AV, Watkins GH, Sayers AR, Jones JR, Clifton-Hadley RS (2012) Geographical association between the genotype of bovine tuberculosis in found dead badgers and in cattle herds. Vet Rec 170: 259. doi: https://doi.org/10.1136/vr.100193.
Carvalho RCT, Vasconcellos SEG, Issa MA, Filho PMS, Mota PMPC, FRd A, Carvalho ACS, Gomes HM, Suffys PN, Figueiredo EES, Paschoalin VMF. Molecular typing of Mycobacterium bovis from cattle reared in Midwest Brazil. PLoS One. 2016;11:e0162459. https://doi.org/10.1371/journal.pone.0162459.
Byrne AW, White PW, McGrath G, James O, Martin SW. Risk of tuberculosis cattle herd breakdowns in Ireland: effects of badger culling effort, density and historic large-scale interventions. Vet Res. 2014;45:109.
Doyle LP, Gordon AW, Abernethy DA, Stevens K. Bovine tuberculosis in Northern Ireland: risk factors associated with time from post-outbreak test to subsequent herd breakdown. Prev Vet Med. 2014;116:47–55. https://doi.org/10.1016/j.prevetmed.2014.06.010.
Menzies FD, Neill SD. Cattle-to-cattle transmission of bovine tuberculosis. Vet J. 2000;160:92–106. https://doi.org/10.1053/tvjl.2000.0482.
White PW, Martin SW, De Jong MCM, O’Keeffe JJ, More SJ, Frankena K. The importance of ‘neighbourhood’ in the persistence of bovine tuberculosis in Irish cattle herds. Prev Vet Med. 2013;110:346–55. https://doi.org/10.1016/j.prevetmed.2013.02.012.
Johnston WT, Gettinby G, Cox DR, Donnelly CA, Bourne J, Clifton-Hadley R, Le Fevre AM, McInerney JP, Mitchell A, Morrison WI, Woodroffe R. Herd-level risk factors associated with tuberculosis breakdowns among cattle herds in England before the 2001 foot-and-mouth disease epidemic. Biol Lett. 2005;1:53–6. https://doi.org/10.1098/rsbl.2004.0249.
Gopal R, Goodchild A, Hewinson G, Domenech RR, Clifton-Hadley R. Introduction of bovine tuberculosis to north-East England by bought-in cattle. Vet Rec. 2006;159:265–71. https://doi.org/10.1136/vr.159.9.265.
Fielding HR, McKinley TJ, Silk MJ, Delahay RJ, McDonald RA. Contact chains of cattle farms in Great Britain. R Soc Open Sci. 2019;6:180719. https://doi.org/10.1098/rsos.180719.
Wolfe DM, Berke O, Kelton DF, White PW, More SJ, O’Keeffe J, Martin SW. From explanation to prediction: a model for recurrent bovine tuberculosis in Irish cattle herds. Prev Vet Med. 2010;94:170–7. https://doi.org/10.1016/j.prevetmed.2010.02.010.
O'Hagan MJH, Matthews DI, Laird C, McDowell SWJ. Herd-level risk factors for bovine tuberculosis and adoption of related biosecurity measures in Northern Ireland: a case-control study. Vet J. 2016;213:26–32. https://doi.org/10.1016/j.tvjl.2016.03.021.
Fitzgerald SD, Kaneene JB. Wildlife reservoirs of bovine tuberculosis worldwide: hosts, pathology, surveillance, and control. Vet Pathol. 2013;50:488–99. https://doi.org/10.1177/0300985812467472.
Byrne AW, Quinn JL, O'Keeffe JJ, Green S, Paddy Sleeman D, Wayne Martin S, Davenport J. Large-scale movements in European badgers: has the tail of the movement kernel been underestimated. J Anim Ecol. 2014;83:991–1001. https://doi.org/10.1111/1365-2656.12197.
Green DM, Kiss IZ, Mitchell AP, Kao RR. Estimates for local and movement-based transmission of bovine tuberculosis in British cattle. Proc R Soc Lond B Biol Sci. 2008;275:1001–5. https://doi.org/10.1098/rspb.2007.1601.
Lahuerta-Marin A, Milne MG, McNair J, Skuce RA, McBride SH, Menzies FD, McDowell SJW, Byrne AW, Handel IG, de C. Bronsvoort BM (2018) Bayesian latent class estimation of sensitivity and specificity parameters of diagnostic tests for bovine tuberculosis in chronically infected herds in Northern Ireland. Vet J 238: 15–21. doi: https://doi.org/10.1016/j.tvjl.2018.04.019
Downs SH, Broughan JM, Goodchild AV, Upton PA, Durr PA. Responses to diagnostic tests for bovine tuberculosis in dairy and non-dairy cattle naturally exposed to Mycobacterium bovis in Great Britain. Vet J. 2016;216:8–17. https://doi.org/10.1016/j.tvjl.2016.06.010.
Byrne AW, Graham J, Brown C, Donaghy A, Guelbenzu-Gonzalo M, McNair J, Skuce R, Allen A, McDowell S. Bovine tuberculosis visible lesions in cattle culled during herd breakdowns: the effects of individual characteristics, trade movement and co-infection. BMC Vet Res. 2017;13:400. https://doi.org/10.1186/s12917-017-1321-z.
Karolemeas K, Mkinley TJ, Clifton-Hadley RS, Goodchild AV, Mitchell A, Johnston WT, Conlan AJK, Donnelly CA, Wood JLN. Recurrence of bovine tuberculosis breakdowns in Great Britain: risk factors and prediction. Prev Vet Med. 2011;102:22–9. https://doi.org/10.1016/j.prevetmed.2011.06.004.
Doyle LP, Courcier EA, Gordon AW, O’Hagan MJH, Menzies FD. Bovine tuberculosis in Northern Ireland: risk factors associated with duration and recurrence of chronic herd breakdowns. Prev Vet Med. 2016;131:1–7. https://doi.org/10.1016/j.prevetmed.2016.06.016.
Conlan AJK, McKinley TJ, Karolemeas K, Pollock EB, Goodchild AV, Mitchell AP, Birch CPD, Clifton-Hadley RS, Wood JLN. Estimating the hidden burden of bovine tuberculosis in Great Britain. PLoS Comput Biol. 2012;8:e1002730. https://doi.org/10.1371/journal.pcbi.1002730.
Houston R. A computerised database system for bovine traceability. Rev Sci Tech. 2001;20:652.
Kamerbeek J, Schouls L, Kolk A, Mv A, Dv S, Kuijper S, Bunschoten A, Molhuizen H, Shaw R, Goyal M, Jv E. Simultaneous detection and strain differentiation of Mycobacterium tuberculosis for diagnosis and epidemiology. J Clin Microbiol. 1997;35:907.
Durr PA, Hewinson RG, Clifton-Hadley RS. Molecular epidemiology of bovine tuberculosis. I. Mycobacterium bovis genotyping. Rev - Off Int Epizoot. 2000;19:675–88.
Team RC. R: a language and environment for statistical computing. R Foundation for statistical computing. Vienna, Austria; 2013.
Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015;67:1–48. https://doi.org/10.18637/jss.v067.i01.
Zuur AI, E.N. ; Walker, N.; Saveliev, A.A.; Smith G.M. (2009) Mixed effects models and extensions in ecology with R. Springer978.
Ribeiro Jr. PJ, Diggle PJ (2001) GeoR: a package for geostatistical analysis. R-NEWS Vol 1: 15–18.
Bivand RW, David WS. Comparing implementations of global and local indicators of spatial association. TEST. 2018;27:716–48.
Broughan JM, Maye D, Carmody P, Brunton LA, Ashton A, Wint W, Alexander N, Naylor R, Ward K, Goodchild AV, Hinchliffe S, Eglin RD, Upton P, Nicholson R, Enticott G. Farm characteristics and farmer perceptions associated with bovine tuberculosis incidents in areas of emerging endemic spread. Prev Vet Med. 2016;129:88–98. https://doi.org/10.1016/j.prevetmed.2016.05.007.
Clegg TA, Good M, More SJ. Future risk of bovine tuberculosis recurrence among higher risk herds in Ireland. Prev Vet Med. 2015;118:71–9. https://doi.org/10.1016/j.prevetmed.2014.11.013.
Doyle LP, Courcier EA, Gordon AW, O'Hagan MJH, Stegeman JA, Menzies FD. Bovine tuberculosis in Northern Ireland: quantification of the population disease-level effect from cattle leaving herds detected as a source of infection. Epidemiol Infect. 2017;145:3505–15. https://doi.org/10.1017/S0950268817002424.
Lahuerta-Marin A, McNair J, Skuce R, McBride S, Allen M, Strain SAJ, Menzies FD, McDowell SJW, Byrne AW. Risk factors for failure to detect bovine tuberculosis in cattle from infected herds across Northern Ireland (2004-2010). Res Vet Sci. 2016;107:233–9. https://doi.org/10.1016/j.rvsc.2016.06.014.
Alvarez J, Perez AM, Bezos J, Casal C, Romero B, Rodriguez-Campos S, Saez-Llorente JL, Diaz R, Carpintero J, de Juan L, Domínguez L. Eradication of bovine tuberculosis at a herd-level in Madrid, Spain: study of within-herd transmission dynamics over a 12 year period. BMC Vet Res. 2012;8:100. https://doi.org/10.1186/1746-6148-8-100.
Carroll JA, Forsberg NE. Influence of stress and nutrition on cattle immunity. Vet Clin North Am Food Anim Pract. 2007;23:105–49. https://doi.org/10.1016/j.cvfa.2007.01.003.
Gallagher MJ, Higgins IM, Clegg TA, Williams DH, More SJ. Comparison of bovine tuberculosis recurrence in Irish herds between 1998 and 2008. Prev Vet Med. 2013;111:237–44. https://doi.org/10.1016/j.prevetmed.2013.05.004.
We wish to thank DAERA, the project funders. We also wish to extend our gratitude to all of the staff in the Veterinary Sciences Division laboratories who undertook the M. bovis genotyping work used herein.
This work was by the Department of Agriculture, Environment and Rural Affairs (DAERA; https://www.daera-ni.gov.uk/) under study 16/3/06 – Bovine TB Molecular Epidemiology – Analysis of Cattle Movements and Optimisation of Epidemiological Investigations. The project funders had no role in study design, data collection and analysis, interpretation of data, writing the manuscript or in decision to publish. CMcC was an AFBI employee involved in collecting the genetic data. CMcC is now an employee of DAERA (the project funder), however CMcC’s contribution occurred when they were an AFBI staff member, and they have had no further role in the project.
Ethics approval and consent to participate
Ethical approval was not required as the data used in this study is routinely collected as part of the NI bTB control programme.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Patches and Divisional Veterinary Offices (DVOs) within NI. Table S1. Summary statistics for herd level variables. Table S2. Summary statistics for patch level variables. Table S3. AIC values for each of the models, compared to null models containing only the “total reactors” variable. Table S4. Final full model for the breakdown level analysis constructed using all data. Table S5. Final full model for the breakdown level analysis constructed using only data from herds with milk licences. Table S6. Final full model for the breakdown level analysis constructed using only data from herds without milk licences. Table S7. Final full model for the herd level analysis constructed using all data. Table S8. Final full model for the herd level analysis constructed using data from herds with milk licences. Table S9. Final full model for the herd level analysis constructed using data from herds without milk licences. Table S10. Final full model for the patch level analysis constructed using data from all herds. Table S11. Final full model for the patch level analysis constructed using only data from herds with milk licences. Table S12. Final full model for the patch level analysis constructed using only data from herds without milk licences. (DOCX 187 kb)