Assessing effects from four years of industry-led badger culling in England on the incidence of bovine tuberculosis in cattle, 2013–2017

The objective was to measure the association between badger culling and bovine tuberculosis (TB) incidents in cattle herds in three areas of England between 2013–2017 (Gloucestershire and Somerset) and 2015–2017 (Dorset). Farming industry-selected licensed culling areas were matched to comparison areas. A TB incident was detection of new Mycobacterium bovis infection (post-mortem confirmed) in at least one animal in a herd. Intervention and comparison area incidence rates were compared in central zones where culling was conducted and surrounding buffer zones, through multivariable Poisson regression analyses. Central zone incidence rates in Gloucestershire (Incidence rate ratio (IRR) 0.34 (95% CI 0.29 to 0.39, p < 0.001) and Somerset (IRR 0.63 (95% CI 0.58 to 0.69, p < 0.001) were lower and no different in Dorset (IRR 1.10, 95% CI 0.96 to 1.27, p = 0.168) than comparison central zone rates. The buffer zone incidence rate was lower for Gloucestershire (IRR 0.64, 95% CI 0.58 to 0.70, p < 0.001), no different for Somerset (IRR 0.97, 95% CI 0.80 to 1.16, p = 0.767) and lower for Dorset (IRR 0.45, 95% CI 0.37 to 0.54, p < 0.001) than comparison buffer zone rates. Industry-led culling was associated with reductions in cattle TB incidence rates after four years but there were variations in effects between areas.

Culling licence criteria were informed by results from the Randomised Badger Culling Trial (RBCT) conducted between 1998 and 2006 9 . The RBCT showed a statistically significant decrease in confirmed TB incidence (where M. bovis has been detected through post-mortem tests in at least one animal from each herd) of 23% (95% Confidence Interval (CI) 12 to 33%, p < 0.001) over four years in 100 km 2 culled areas compared to non-culled areas 6,10 . An increase of 25% (95% CI 1 to 56%, p = 0.057) in confirmed TB was also observed on two km wide land surrounding the culled areas relative to land surrounding non-culled areas, during the period culling was conducted 10 .
Licences for culling have been issued for two areas (in Gloucestershire and Somerset) from 2013, one area (in Dorset) from 2015, seven areas from 2016 and 11 areas in both 2017 and 2018. All areas are located within the HRA except one area in 2018, which is located in the Low Risk Area for TB in England. The culling areas are larger than the RBCT areas and also differ in a number of other important ways. For example, culling during the RBCT was carried out by government and involved cage-trapping badgers; with the whole area simultaneously trapped over a period of approximately 10 days. In contrast, current badger culls are carried out by the farming industry, and involve a combination of cage trapping and controlled shooting (where free roaming badgers are shot at night), with effort spread over a period of six weeks or more. Consequently, it is unclear whether current badger culls would produce similar changes in cattle TB incidence to those observed during the RBCT.
The Animal and Plant Health Agency (APHA) has been commissioned to monitor and evaluate the effects of culling on TB incidence in cattle herds in culling areas. The analysis of effects in cattle herds over two years since badger culling began in Gloucestershire and Somerset has already been reported 11 . This analysis showed a statistically significant decrease in Officially Tuberculosis Free Withdrawn (OTF-W) incidence in both culled areas compared to areas with no culling. An OTF-W incident is a new outbreak of TB in a herd disclosed by field testing or slaughterhouse surveillance and where M. bovis is detected through post-mortem tests in at least one animal in the herd. The OTF-W incidence rate ratio (IRR) was 0.79 for Somerset (95% CI 0.72 to 0.87) and 0.42 for Gloucestershire (95% CI 0.34 to 0.51) after adjustment for confounding factors. An increase in incidence was observed in the two km buffer zone around the Somerset culled area (IRR 1.38, 95% CI 1.09 to 1.75) but not in Gloucestershire (IRR 0.91 95% CI 0.77 to 1.07).
The statistically significant decrease in cattle TB incidence associated with culling in Gloucestershire and Somerset over two years was unexpected given the statistical power of the study 12 . The aims of this new analysis is to determine firstly whether a decrease in TB incidence rates in cattle associated with culling has been sustained in Gloucestershire and Somerset; secondly whether a similar effect can be detected after two years in Dorset; and thirdly whether an effect on TB incidence can be detected in two km wide buffer zones surrounding culled areas. The null hypothesis being tested is that TB incidence is the same in the intervention areas (where culling has been conducted) and their comparison areas in the years since badger culling began 11 .

Methods
Intervention area coverage, cohort of herds exposed to culling, follow-up period. The intervention areas in Gloucestershire, Somerset and Dorset selected by industry are 311, 256 and 223 km 2 in size respectively and located in the HRA of England (Fig. 1). Culling took place during a six-week period each autumn in each area. In this analysis, the culling areas are referred to as the central zones. Each central zone has been allocated a two km wide buffer zone from its boundary using ArcGIS (ESRI Releases 10.0-10.3 Redlands, CA, USA). Cattle herds in existence in the APHA surveillance database (Sam) in the central and buffer zones when culling started (the baseline date, Table S1) formed the cohort of herds exposed to badger culling and monitored for changes in TB incidence 11 . The follow-up period in this analysis was four years from autumn 2013 to autumn 2017 in Gloucestershire and Somerset and two years from autumn 2015 to autumn 2017 in Dorset (Table S1).
Selection of comparison areas. The methodology and criteria for the selection of comparison areas for the Gloucestershire and Somerset intervention areas have been reported in detail in Brunton et al. 11 . The same comparison areas for Gloucestershire and Somerset are used in the current analysis. Comparison areas for Dorset, which was licenced two years later, have been selected for the current analysis using the methodology reported previously 11 . Geographical Information System programming was used to generate a population of potential comparison areas within the HRA of the same size and shape as the Dorset area. Each area was ranked by degree of similarity to the Dorset area on the basis of the number of OTF-W incidents in the previous three years, the number of herds in the area, the median herd size, the proportion of the area (if any) that was in a RBCT proactive cull area and distance to the Dorset intervention area. The 10 areas that most closely matched the Dorset intervention area and did not contain land already within Gloucestershire or Somerset intervention areas or land within two km of any intervention area boundary (at the time of the baseline date) were selected. Cattle herds in existence in comparison area central and buffer zones on the baseline date formed the cohort of herds compared to herds in the intervention area central and buffer zones.
Main outcome variable and time at risk. An area level analysis was conducted. The main outcome variable was the OTF-W incidence rate in cattle herds for each area. Effects were also estimated for all TB incidents (OTF-W plus OTFS (Officially Tuberculosis Free Suspended)), where reactors to the Single Intradermal Comparative Cervical Tuberculin (SICCT) test had been detected in a herd but infection had not been confirmed by post-mortem tests. Herd years at risk (HYR) were calculated using results from whole herd tests as the sum of the time cohort herds were unrestricted and therefore at risk of new infection (a new TB incident) during the period of interest 13 . HYR data used previously 11 were updated to remove implausible values for time at risk. This meant that the sum of time at risk for all herds in an area as opposed to the median estimate of time at risk could be used thereby bringing the methodology in line with national statistics reporting 14 . Herds could move out of an area during the follow-up period because either the farmer moved the herd off the land in the area to which the farm was originally assigned to land outside of that area, or the original geo-reference for the location of the herd was found to be incorrect and the subsequent correction placed the herd outside of the area. Herds that moved out of an area but were still in existence during the follow-up period contributed incidents and HYR for the calculation of area-level incidence rates for the area in which they were originally located on the baseline date. Similar to an "intention to treat" analysis, the purpose was to reduce any bias due to differences in the number of herds leaving areas between intervention and comparison areas.
Confounding factors. Comparison areas, although selected based on their similarity to intervention areas will vary to a greater and lesser extent from the intervention areas in factors associated with TB incidence in cattle 14,15 . Factors (attributes of the herds and the area) that could be associated with TB incidence and with selection of an area for badger culling were extracted from the APHA cattle surveillance database and other sources and summarised for each area. These factors included all those in Brunton et al. 11 plus a new variable describing badger density that included new survey data 16,17 (see Supplement). As a result of data checks, corrections were made to two variables used previously; to the proportion of land in flood zone 3 and the proportion of land that had been in a proactive area of the RBCT.
Overlap between areas. Early in the project design phase, it was apparent that comparison areas could be overlapped by intervention areas licensed in later years and that comparison areas to different intervention areas might also overlap 18 . Rules were developed, prior to any analyses, to address allocation of herds between overlapping areas and for summarising confounder data for overlapping areas (Fig. S1, Tables S2, S3). In the case where more than 25% of a previously identified comparison area was overlapped by a more recently licensed intervention area (central and/or buffer zone), that comparison area was removed from the analysis. In the case where a comparison area was overlapped by an intervention area but at least 75% of it remained as one continuous area, the comparison area boundaries (central and/or buffer zone) were redrawn so that its central zone was at least two km from the central zone of the intervention area. In the case of overlapping comparison areas central zones, herds were randomly allocated between the comparison areas and confounder data associated with the herds in each area recalculated. In the case of a comparison area buffer overlapping a comparison area central zone, the herds remained in the central zone and the boundary to the buffer was redrawn. In the case where buffer zones to different comparison areas overlapped, herds were randomly allocated between the buffers and confounder data associated with herds recalculated. Statistical analysis. The analytical approach was the same as that followed in Brunton et al. 11 . Effects in central zones and buffer zones were estimated separately because previous research has suggested that different effects in cattle TB in the areas surrounding the culling areas compared to areas where culling is conducted 10,11 . Crude OTF-W and all TB incidence rate ratios (IRRs) comparing incidence rates between intervention and comparison areas were calculated for three 12-month periods prior to the baseline date for all areas; and 12-month periods after the baseline date until the end of the follow-up (four years for Gloucestershire and Somerset and two years for Dorset).
Multivariable Poisson regression was conducted to measure changes in TB incidence over four or two years in central and buffer zones whilst controlling for confounding factors. IRRs were estimated for all variables as the exponent of the estimated Poisson regression coefficient. Matching of comparison areas to either Gloucestershire or Somerset was achieved through the inclusion of an area variable. The mathematical description of the model is as follows: where μ i = the expected number of OTF-W or TB incidents for the areas during the time period under analysis, with the observed number being Poisson distributed with mean μ i , β 0 = the intercept. Coefficients β 1 to β p are multiplied by explanatory variables x 1 to x p including an intervention effect variable matching each intervention area (Gloucestershire, Somerset and Dorset) to its comparison areas.
The natural logarithm was taken of explanatory variables that were counts such as herd years of risk, herd size and number of badgers removed as part of historical badger control operations. Where a value was zero, 0.5 was added to the count before taking the logarithm. Effects by culling area (Gloucestershire, Somerset or Dorset) were measured by including a single intervention status variable and variables representing the geographical area (coded as Gloucestershire, Somerset or Dorset) which included both the intervention areas and their matched comparison areas. Statistical interaction by area was tested in the final models.
Using the dataset compiled for the current follow-up, OTF-W incidence rates were initially estimated in models including the same explanatory factors as in Brunton et al. 11 . Models were then rebuilt to investigate whether a better fit might be achieved with different covariates. Initial models included all explanatory factors shown to differ between intervention areas in the descriptive analyses. Predictors not statistically significantly associated with OTF-W incidents in the initial model were then eliminated. The subsequent model included terms for area, herd years at risk, OTF-W incidence rate in the three years up the baseline date when culling started and predictors that had been statistically significantly associated with OTF-W incidents in the initial model. Other predictors were then re-inserted one at time in a forward selection procedure whilst model fit was assessed using deviance 19 and Pearson statistics 20 and Akaike's Information Criterion 21 (AIC). Herd size and OTF-W incidence rate over the three years prior to baseline when culling started were retained in all models. Residuals were examined in the best fitting models and model stability examined by removing comparison areas that had the largest leverage statistics. Separate models were built to estimate effects for the individual year's one, two, three and four since culling began. All TB incidence rates (OTF-W plus OTF-S incidents) were estimated in models including the same explanatory variables used in the models for OTF-W incidence rates.
Sensitivity analyses included (1) rerunning models reported in Brunton et al. 11 for two years follow-up in Gloucestershire and Somerset with the new HYR, badger density, flood zone 3 and RBCT proactive area data (2) examining effects in final models with the removal of one comparison area at a time and (3) comparing effects in herds in the central zones subdivided into outer zone herds located within two km of the boundary and inner zone herds located two or more km inside the boundary.
The statistical analyses were conducted using Stata release 14.1 (Stata Corp., Texas, USA). Standard errors on effect estimates were calculated controlling for any non-Poisson variation (greater or smaller variance in the counts of OTF-W incidents between areas than would be expected by chance) using a robust estimator of variance in Stata 22,23 . Probability values (p values) of less than five percent (p < 0.05) were interpreted as statistically significant.

Overlap between comparison areas and intervention areas licensed in 2016. Four comparison
areas (including central and buffer zones) were excluded from the analyses because more than 25% of their land was overlapped by land in one of the eight intervention areas licenced between 2015 and 2016 (Tables S2, S3). Gloucestershire and Somerset intervention areas lost one comparison area each and Dorset lost two. Additionally, the land area and number of herds in the central zones of two comparison areas for Gloucestershire and one for Dorset were reduced due to intervention area overlap (Table S2). Land area and the number of herds in 21 buffer zones were also reduced due to intervention area or comparison area central zone overlap (Table S3) www.nature.com/scientificreports www.nature.com/scientificreports/ Baseline characteristics in intervention and comparison areas. Dorset intervention area herds were larger and more likely to be dairy than Gloucestershire and Somerset herds (Tables 1, 2). The estimated badger density prior to industry-led culling was almost twice as high in the Dorset intervention area compared to the Somerset and Gloucestershire areas and fewer badgers had been removed in historical operations. More of the Gloucestershire intervention area land was classified as flood zone 3 than in Somerset or Dorset. Somerset comparison areas were on average located further from the intervention area than Gloucestershire and Dorset comparison areas.
Comparison areas were fairly similar to their matched intervention areas (Tables 1, 2). However, a higher proportion of herds were dairy in Dorset and there was almost twice as much land in flood zone 3 in Gloucestershire compared to their respective comparison area central zone averages. A lower proportion of Gloucestershire intervention area central zone farms had land parcels outside the area compared to its comparison central zone average. There was also twice as much motorway in the Gloucestershire intervention area buffer zone compared to its comparison area buffer zone average.

Crude OTF-W incidence rates in intervention and comparison areas.
In the four years since the baseline date when culling started there were 64 and 80 OTF-W incidents in the Gloucestershire and Somerset intervention central zones respectively. In the two years since the baseline date there were 44 OTF-W incidents in the Dorset intervention central zone. There were 34 and 40 OTF-W incidents respectively during four years in the Gloucestershire and Somerset intervention buffer zones and 13 incidents during two years in the Dorset intervention buffer zone.
Incidence rates over one and three years before the baseline date in each intervention area were correlated, particularly OTF-W incidence rates in the Dorset central zone (Spearman correlation coefficient 0.867, Spearman p value = 0.003). Annual crude OTF-W incidence rates were lower in the Gloucestershire intervention central zone compared to comparison areas from two years prior to the baseline date to four years post baseline when culling started (Table 3); with the strongest relative decline in the intervention area central zone in the fourth year after culling started. Rates in the Somerset intervention central zone were higher than the mean for comparison zones until the second year after culling started. Other differences between OTF-W between intervention and comparison area central and buffer zones were not as strong. Rates in comparison buffers were lower than rates in the intervention area buffers throughout the follow-up period in Gloucestershire and Dorset. Differences in rates for all TB incidents showed a similar pattern to the OTF-W rates (Table S4).
After four years follow-up, 19.1% of the Gloucestershire and 15.6% of the Somerset herds in the central zone were no longer located in their respective areas. The percentage of Dorset herds that had moved after two years follow-up was 3.8% (Table S5). www.nature.com/scientificreports www.nature.com/scientificreports/ Adjusted OTF-W incidence rates in central zones of intervention and comparison areas. The fit of models for OTF-W incidence rates with the confounding factors reported in Brunton et al. 11 was good for effects in the central zones over four years since baseline (deviance goodness-of-fit (gof) p = 0.506, Pearson gof p = 0.511) but poor for effects over two years (deviance gof p = 0.001, Pearson gof p = 0.001).
The best fitting model for effects in the central zone over four years since baseline showed statistically significantly lower OTF-W incidence rates in cattle herds in both Gloucestershire and Somerset intervention central zones relative to comparison central zones (Table 4, model A). This model had an AIC of 149.492 compared to an AIC of 152.847 for a model fitted to the same dataset using the same confounding factors as in Brunton et al. 11 . The effect was strongest in Gloucestershire where the central estimate was 66% lower than in comparison areas compared to 37% lower in Somerset (p = 0.038 for interaction by area). Comparison area labelled WS03 had the greatest leverage but its removal had little effect on estimates (Table S6).
The best fitting model for effects in the central zones over two years since baseline included a variable with quintiles of the distribution of the total number of badgers culled historically (Table 5, model C). This model had an AIC of 197.112 compared to an AIC of 221.735 for a model fitted to the same dataset using the same confounding factors as in Brunton et al. 11 . OTF-W incidence rates for Gloucestershire and Somerset intervention central zones were statistically significantly lower than in comparison central zones. There was no difference between OTF-W incidence rates in the Dorset intervention central zone compared to comparison central zones (p = 0.001 for interaction by area). Comparison area labelled WS01 had the greatest leverage but its removal had little effect on estimates (Table S6).
The central estimates for annual IRRs declined in Somerset and Gloucestershire with years since culling started (Fig. 2). The OTF-W incidence rates for the Gloucestershire intervention were lower than in comparison central zones each of the four years since baseline (p < 0.001). Rates for the Somerset intervention were no different to comparison rates in the first year since baseline but were statistically significantly lower in subsequent years (p < 0.001). Differences between annual rates in the Dorset intervention central zone and comparison central zones were not statistically significant.

Adjusted OTF-W incidence rates in buffer zones of intervention and comparison areas. The fit
of models for OTF-W incidence rates with the confounding factors reported in Brunton et al. 11 was not good for effects in the buffer zones over four years since culling started (deviance gof p = 0.035, Pearson gof p = 0.037) but good for effects over two years (deviance gof p = 0.443, Pearson gof p = 0.394).
The best fitting model for effects in the buffer zones over four years since baseline showed that rates were around 36% lower in the Gloucestershire intervention buffer zone relative to comparison buffer zones (p < 0.001, Table 4, model B). This model had an AIC of 139.210 compared to an AIC of 147.634 for a model fitted to the same dataset using the same confounding factors as in Brunton et al. 11 .There was no difference between incidence rates in the Somerset intervention and comparison buffer zones but also no evidence for interaction by intervention area (p = 0.396). The best fitting model for effects in the buffer zones over two years showed that Dorset rates were around 55% lower than in comparison zones (p < 0.001, Table 5, model D). This model had an AIC of 170.884 compared to an AIC of 176.801 for a model fitted to the same dataset using the same confounding factors as in Brunton et al. 11 . The incidence rate was lower for the Gloucestershire buffer compared to comparison buffers but the difference was not statistically significant. There was no statistically significant difference in the rates in the Somerset buffer zone relative to comparison buffer zones. The interaction p value was highly statistically significant showing effects differed by intervention area (p < 0.001). Comparison area labelled WS03B had the greatest leverage in both the two year and the four year model for buffer zones but its removal did not change the interpretation of effects (Table S7).
Incidence rates in the Gloucestershire and Dorset intervention buffer zones were lower than in comparison buffer zones each year since culling started (Fig. 3). Incidence rates in the Somerset intervention buffer zone were higher than in comparison areas (p = 0.001) in the first year after culling started but were similar to comparison zones in years two, three and four.

Adjusted all TB (OTF-W plus OTF-S) incidence rates in central and buffer zones. The adjusted
IRRs for all TB incidents in the Somerset and Gloucestershire intervention central zones over four years were statistically significantly lower than in comparison areas zones (Table S8) but not over two years (Table S9). The Dorset intervention central zone had a statistically significantly higher incidence rate over two years than comparison areas zones (Table S9). Rates in Gloucestershire intervention buffer zone were statistically significantly lower than in comparison areas over four years, but not over two years. There were no differences between rates in the Somerset and Dorset intervention buffer zones to comparison area buffer zones.
Sensitivity analyses. IRR estimates for the effects over two years in Somerset and Gloucestershire using the updated estimates for HYR, percentage of land exposed to RBCT proactive culling and badger density as opposed to badger sett density were very similar to those previously reported for the central zone in Brunton et al. 11 (Table S10). However, the increase in incidence in the Somerset buffer zone was smaller and no longer statistically significant.
Removal of comparison areas from the final four and two year central and buffer zone models (  www.nature.com/scientificreports www.nature.com/scientificreports/ The Somerset intervention central zone had a higher proportion of herds (68.2%) located within two km of the boundary to the zone than Gloucestershire (53.5%) or Dorset (51.6%). There was no obvious pattern to differences in OTF-W incidence rates between inner and outer zones to the central zones (Fig. 4).
The central estimates for IRRs were the same but confidence intervals were generally wider in models without control for non-Poisson variation compared to models with control. The OTF-W incidence rates over four years in central zones were still statistically significantly lower in Gloucestershire and Somerset relative to comparison area zones as were OTF-W incidence rates in the Gloucestershire buffer zone (compare Table 4 to Table S11). Central estimates for OTF-W rates in Gloucestershire and Somerset central zones over two years were lower than comparison central zones but not statistically significantly so; there was no change to the interpretation of effects in Dorset (compare Table 5 to Table S12).

Discussion
The results from this study showed that there were statistically significant decreases in cattle TB incidence in the Gloucestershire and Somerset intervention areas after four years of culling, consistent with an earlier analysis based on two years of badger culling 11 . The decreases in confirmed TB (OTF-W) incidence rates relative to comparison areas observed were 66% (95% CI 61 to 71%) and 37% (95% CI 31 to 42%) in Gloucestershire and Somerset respectively. However, there was no change in OTF-W incidence after two years of culling in Dorset central zone. Decreases in OTF-W incidence rates in the Gloucestershire and Dorset buffer zones relative to comparison area buffer zones were also observed, which were unexpected. Other research suggests that an increase in cattle TB may occur in areas surrounding culling areas due to increased M. bovis transmission caused by perturbation (increased mixing) of badger populations 24 .
Establishing causality between an intervention and a disease that involves transmission between two animal species is challenging. Each component of the causal pathway is affected by environmental factors that differ between areas such as the spatial distribution of badger habitat, the location of cattle and the effectiveness of the culling operations. Furthermore, evaluating control policies that change over time in response to new information and political climate adds additional challenges.
Since two original pilot culls of 2013, a further 29 licences for culling in other areas of the HRA have been issued. We included 10 comparison areas per intervention area in the study design because we had no control over the selection of future culling areas and anticipated comparison area land would be lost to newly licenced areas over time 18 . Four comparison areas have been lost and the land area of 21 buffer zones reduced because of overlap by culling areas licensed after the culling areas in Gloucestershire, Somerset and Dorset. Land that became part of a new culling area had to be excluded entirely from our analysis because herds on that land became exposed to culling. Incidence rates throughout the follow-up period were only calculated from the herds in existence on comparison central and buffer zone land when culling started (the baseline date) that also did not become licensed for culling at any point in the follow-up period. This was to ensure that the TB incidence rates  www.nature.com/scientificreports www.nature.com/scientificreports/ for comparison zones were from herds that were not directly exposed to culling at any point. The total number of TB incidents recorded on comparison area land will be lower than if none of the land had been lost to culling. However the incidence rates (incidents per herd years at risk) should be an accurate reflection of disease rates in herds on comparison area land.
The RBCT was a rigorously conducted scientific trial and the best evidence of what might be achievable from badger removal on a large scale in England 9 . As a randomised controlled trial, it is less vulnerable to confounding by differences in the distribution of TB risk factors between areas, than the current study. We attempted to control for effects from confounding by adjusting for factors known to be associated with TB risk. However, these factors were derived from routinely collected surveillance data and may be subject to misclassification biases. There may have been factors associated with the granting of culling licenses e.g. greater uptake of biosecurity that could affect TB risk. We could not adjust for unknown confounders nor factors where we did not have information from both intervention and comparison areas. A reduction of 23% (95% CI 12 to 33%) in confirmed TB incidents in areas subject to at least four years of widespread systematic culling compared to non-culled areas was observed in the RBCT 6,10 . Reductions over four years in Gloucestershire and Somerset were larger. However, this increase may be because industry-led culling was conducted over a longer period each year and the areas culled were larger than in the RBCT.
Some heterogeneity in effects is to be expected due to chance and other factors e.g. differences in cattle distribution, the burden of infection in badger populations and culling coverage. There were some inconsistencies between our results and those of the RBCT as well as between areas. We did not find evidence for a stronger beneficial effect on cattle TB from culling with increasing distance from the boundary of the central zones. There was initially evidence for a larger beneficial effect with increasing distance inside culling areas in the RBCT, although the trend was non-statistically significant in later analyses 10,25 . Buffer zone incidence rates were statistically significantly lower in the Gloucestershire and Dorset buffer zones compared to comparison buffer zones throughout the follow-up period but statistically significantly higher in the Somerset buffer during the first year of the follow-up. In the RBCT, a 25% increase (95% CI -1% to 56%) in incidents was initially observed in the buffer zones around culled areas, which reduced over time 10 .  One third of buffer zones to Gloucestershire, Somerset and Dorset culling areas were reduced in size to take account of culling areas licensed later. It is likely that both cull coverage of areas and boundary permeability will have affected badger ranging behaviour. Up to 20% of intervention area cohort herds were no longer in existence in the area after four years follow-up and other herds had moved into the areas. However, trends in crude incidence rates in the cohort herds are similar to those for herds in existence in the areas on the baseline date and anniversaries of the baseline date 26 .
There has been debate about the delay that might be expected between culling badgers and an observable effect on the incidence of TB in cattle TB 27 . Behavioural data show that local reductions in badger density can cause badgers to alter their ranging behaviour within a few weeks 24 . Cattle most likely acquire infection from badgers from indirect contact, e.g. from pasture or contaminated feed 28 . Experimental work suggests pulmonary exposure of cattle to very low numbers of colony forming units of M. bovis will result in a positive SICCT test response within 12 weeks 29 . However detection is also dependent on the frequency of TB surveillance tests in the cattle herds (annual in the HRA during the follow-up period). There was no decline in TB incidence in cattle, between the first and second year of follow-up in Gloucestershire and Dorset central zones although there was in Somerset. This may reflect area Figure 2. Adjusted incidence rate ratios (central estimates and 95% confidence intervals) of Officially Tuberculosis Free Withdrawn (OTF-W) incidence in cattle herds in the central zones of intervention areas compared to herds in the central zones of comparison areas. Annual incidence rate ratios (IRRs) for Gloucestershire and Somerset compared to comparison areas were estimated in Poisson regression models adjusting for the OTF-W incidence rate in the three years prior to baseline, median herd size, percentage of herds that were dairy, distance between intervention areas and comparison areas, the log transformed total number of badgers removed between 1972 and 2006 and the percentage of farms with at least one land parcel in the central zone. The annual IRRs for Dorset compared to comparison areas were estimated in models adjusting for effects in Somerset and Gloucestershire, the OTF-W incidence rate in the three years prior to baseline, median herd size, percentage of land in flood zone 3, distance between intervention areas and comparison areas and quintiles for total numbers of badgers removed between 1972 and 2006. www.nature.com/scientificreports www.nature.com/scientificreports/ differences in reduction in infection transmission and detection of infection. Furthermore, infection prevalence in the wildlife reservoir will vary between areas, affecting the relative impact of the badger removal policy.
As more areas are licensed for culling, they are likely to be less similar to RBCT areas and consequently the RBCT results may be less predictive. In contrast to Gloucestershire and Somerset, none of the Dorset area was within the RBCT. Dorset also contained larger herds, more dairy herds and had a higher estimated baseline density of badgers. The most recent badger population estimates imply that the proportion of the Dorset badger population culled in the first year at least, may have been lower than in Somerset and Gloucestershire (Table S13). For these reasons a reduction in the incidence of TB in cattle may take longer to emerge in Dorset. A longer follow-up increases power to detect effects 12 and supports delaying analyses to a time point where there is strong evidence for sufficient power for robust evaluation of effects. The decline in cattle TB for the Gloucestershire and Somerset central zones over four years was slightly stronger and more robust than after two years of follow-up.
Martin et al. 30 concluded that infected badgers explained 9-19% of cattle TB incidents in the East Offaly trial in Ireland. A more recent analysis of RBCT data estimated that 5.7% (95% CI 0.9-25%) of transmission to cattle herds is from badgers 31 . This alongside the level of reduction in OTF-W incidence rates observed in the current analysis suggests that there are other mechanisms at play that amplify effects associated with badger controls. Implementing culling may lead to greater focus on cattle controls, TB testing quality and implementation of biosecurity. . Adjusted incidence rate ratios (central estimates and 95% confidence intervals) for Officially Tuberculosis Free-Withdrawn (OTF-W) incidence in cattle herds in the buffer zones of intervention areas compared to herds in the buffer zones of comparison areas. Annual incidence rate ratios (IRRs) for Gloucestershire and Somerset compared to comparison areas were estimated in Poisson regression models adjusting for the OTF-W incidence rate in the three years prior to baseline, median herd size, percentage of urban land, the log transformed total number of badgers removed between 1972 and 2006 and the percentage of farms with all land in the buffer zone. Annual IRRs for Dorset compared to comparison areas were estimated in models adjusting for effects in Somerset and Gloucestershire, the OTF-W incidence rate in the three years prior to baseline, median herd size, percentage of urban land, length of motorway, distance between intervention areas and comparison areas and the percentage of farms with all land in the central or buffer zone.
The results presented here, which were reasonably consistent with the RBCT, show that a culling policy implemented by the farming industry can result in statistically significant reductions in the incidence of cattle TB. However, given the observational nature of the study we cannot exclude entirely biases in our results due to for example, unknown or unmeasured confounding. We recommend that evaluation of the effects from culling continues. We need to know whether the beneficial effects that have been observed on cattle TB continue and can also be observed in other culled areas. We also need to understand why an increase in TB incidence rates in cattle has not been detected in buffer zones surrounding culling areas. The analysis highlights the difficulties in predicting effects from large scale interventions aimed to reduce infection transmission between animal species. . Adjusted incidence rate ratios (central estimates and 95% confidence intervals) for Officially Tuberculosis Free-Withdrawn (OTF-W) incidence in cattle herds in inner and outer zones of intervention central zones compared to in the inner and outer zones of the comparison area central zones. Inner zone herds were from farms with a geolocation two or more km inside the boundary to the central zone. Outer zone herds were from farms with a geolocation within two km of the boundary to the central zones. Annual incidence rate ratios (IRRs) for inner zone herds over four years were estimated in models adjusting for intervention area, OTF-W incidence rate in the three years prior to baseline, median herd size, percentage of flood zone 3 land and quintiles for total numbers of badgers removed between 1972 and 2006. Annual incidence rate ratios (IRRs) for inner zone herds over two years were estimated in models adjusting for intervention area, OTF-W incidence rate in the three years prior to baseline, median herd size, percentage of flood zone 3 land, percentage of dairy herds, badger density and the percentage of farms with all land in the central zone. Annual incidence rate ratios (IRRs) for outer zone herds over two years were estimated in models adjusting for intervention area, OTF-W incidence rate in the three years prior to baseline, median herd size, percentage of urban land, percentage of farms with at least one land parcel in area and percentage of land that was in a Randomised Badger Culling Trial (RBCT) Proactive area. Annual incidence rate ratios (IRRs) for outer zone herds over four years were estimated in models adjusting for intervention area, OTF-W incidence rate in the three years prior to baseline, median herd size, percentage of dairy herds, percentage of farms with at least one land parcel in area and percentage of land that was in a RBCT proactive area.