Evaluation of Risk Factors Associated With Herds With an Increased Duration of Bovine Tuberculosis Breakdowns in Castilla y Leon, Spain (2010–2017)

The persistence of bovine tuberculosis (bTB) in certain cattle herds is a major concern in countries pursuing disease eradication worldwide. The chronic nature of the disease, the lack of performance of diagnostic tools, and the presence of wildlife reservoirs may lead infected herds to require longer periods to achieve the officially tuberculosis-free (OTF) status. Here, we evaluated the impact of farm and breakdown characteristics on the probability of disease persistence in infected farms in Castilla y Leon, a bTB-endemic region of Spain, using survival and logistic regression models. Data from bTB breakdowns occurring in 3,550 bTB-positive herds detected in 2010–2017 were analyzed. A multivariable Cox proportional hazards model was fitted using time to recover OTF status as the response variable, and a multivariable logistic regression model using the chronic status (yes/no) for herds experiencing particularly long breakdowns as the outcome variable was also used. Both analyses revealed that county-level bTB herd prevalence, herd size, number of incoming animals in the previous 3 years, number of skin test reactors in the disclosing test, and number of days between the disclosing and follow-up tests were associated with increased breakdown duration. Production type was not consistently associated with chronic infection, suggesting that once infected, it is not a significant predictor of outbreak duration beyond the initial stages of the breakdown. Province-level location and number of animals that are bacteriology-positive also affected significantly the expected herd breakdown duration, but their effect became less significant over time. Risk factors identified in this study may help to identify herds more prone to suffer chronic bTB infection that may require additional control measures early on in a breakdown.

The persistence of bovine tuberculosis (bTB) in certain cattle herds is a major concern in countries pursuing disease eradication worldwide. The chronic nature of the disease, the lack of performance of diagnostic tools, and the presence of wildlife reservoirs may lead infected herds to require longer periods to achieve the officially tuberculosis-free (OTF) status. Here, we evaluated the impact of farm and breakdown characteristics on the probability of disease persistence in infected farms in Castilla y Leon, a bTB-endemic region of Spain, using survival and logistic regression models. Data from bTB breakdowns occurring in 3,550 bTB-positive herds detected in 2010-2017 were analyzed. A multivariable Cox proportional hazards model was fitted using time to recover OTF status as the response variable, and a multivariable logistic regression model using the chronic status (yes/no) for herds experiencing particularly long breakdowns as the outcome variable was also used. Both analyses revealed that county-level bTB herd prevalence, herd size, number of incoming animals in the previous 3 years, number of skin test reactors in the disclosing test, and number of days between the disclosing and follow-up tests were associated with increased breakdown duration. Production type was not consistently associated with chronic infection, suggesting that once infected, it is not a significant predictor of outbreak duration beyond the initial stages of the breakdown. Province-level location and number of animals that are bacteriology-positive also affected significantly the expected herd breakdown duration, but their effect became less significant over time. Risk factors identified in this study may help to identify herds more prone to suffer chronic bTB infection that may require additional control measures early on in a breakdown.
Keywords: bovine tuberculosis, Mycobacterium bovis, chronic breakdowns, risk factors, cattle, case-control study, survival analysis INTRODUCTION Bovine tuberculosis (bTB) is a zoonotic disease affecting cattle caused by members of the Mycobacterium tuberculosis complex, mainly Mycobacterium bovis and Mycobacterium caprae, which has a major impact in the economy of affected countries due to its effect on trade. Despite being eradicated in several countries (1,2), bTB is still a challenge in many others (3). In Spain, bTB is still present in several regions, and although the implementation of a national eradication program has led to a decrease in the herd prevalence in the last decades, current levels are still similar to those recorded 17 years ago (2.2% in 2002 and 2.3% in 2018) (4).
Accurate diagnosis of bTB in live animals is often difficult, and several factors that influence test performance and can hence lead to possible diagnostic failures have been identified. These include the intrinsic limited sensitivity of currently available tests, the choice of the diagnostic cutoff, the test procedure, the disease stage of infected animals, the possible desensitization to the test in the case of the skin test, the existence of host or pathogen genetic variations, the occurrence of cross-reactions, and the effect of concurrent infections (5)(6)(7)(8). Other factors that impede disease eradication are survival and persistence of M. bovis in the environment and the aggregation at communal water and food sources that can promote closer contacts between cattle and may increase the likelihood of contact with infected wildlife reservoirs (9).
An extensively assessed feature of bTB is its persistence in certain herds, in terms of either herd recurrence or prolonged periods of restriction (10,11). These may imply reinfection, which could be attributed to a local source (such as a contaminated environment, infected wildlife, or farm-to-farm contacts with infected neighboring herds) or to the entry of undetected infected animals and/or ongoing transmission due to residual (persistent but undetected) infection (12). Infection persistence in a herd may indicate a test failure, thereby allowing false-negative infected animals not only to remain in the population but to potentially act as an ongoing source of infection to other herds and wildlife (13,14).
A large proportion of bTB-positive herds in Spain is located in areas where conditions that favor disease persistence are relatively common, such as the presence of potential wildlife reservoirs (mostly red deer and wild boar) and the predominance of extensively managed bullfighting and beef herds (15,16). The Castilla y Leon autonomous community, which holds 20% of Spain's total cattle population, is classified as a high-prevalence region (>1%) in the country (4) and has areas that meet the requirements described above to hamper the progress of the eradication program. Out of 3,550 positive herds detected during the period 2010-2017, bTB infection was confirmed through bacteriological culture in 41% of them. Of these, >70% tested bTB positive for 2 or more years, showing that a major part of the bTB burden detected in the region is concentrated on a subset of herds. A study performed on the cattle movement network demonstrated that infected herds in this region were clustered in space but not in the movement network, suggesting that factors other than movements may be related to disease introduction and maintenance in at least a proportion of positive farms (17).
Although several studies have been conducted on herdlevel risk factors for bTB in Spain (18,19), to date, factors contributing to these persistently infected herds experiencing particularly long outbreaks have not been clarified. To improve our understanding on why certain herds have a higher risk of experiencing prolonged bTB breakdowns, we examined the impact of farm characteristics (such as production type, herd size, and animal trade flow) and bTB breakdown-specific variables (such as results in the initial bTB tests and bTB prevalence in the region) on the duration of bTB breakdowns during the 2010-2017 period in Castilla y Leon. This study may provide useful information on the characteristics that influence the persistence of bTB in highly prevalent areas and help in the design of targeted strategies to manage bTB-positive herds in Spain.

Bovine Tuberculosis Eradication Program in Castilla y Leon
In accordance with EU Directive 64/432/EEC (20), the Spanish bTB eradication program is based on test and slaughter surveillance and live animal testing using in Castilla y Leon the single intradermal test (SIT) in all herds (except certain fattening herds in which compulsory testing may not be required if these only move cattle to the slaughterhouse) plus the interferongamma (IFN-γ) assay as a complementary test in infected herds to maximize the diagnostic sensitivity. Briefly, all animals above 6 weeks old in tested herds are subjected to routine SIT herd tests by intradermal inoculation of 0.1 ml of bovine purified protein derivative (PPD) (Cz Veterinaria, Porriño, Spain) in the anterior neck area with a frequency dependent on the bTB prevalence in their local area (ranging from 1-2 every year); after 72 h, animals with a >2-mm increase in skinfold thickness and/or presence of necrosis, edema, exudation, or inflammation of lymph nodes peripheral to the inoculation site are considered reactors (severe interpretation), culled within the following 15 days and subjected to postmortem analysis. Subsequently, OTF status is either suspended or withdrawn, and movement restrictions are applied. Farms confirmed as infected (≥1 positive animal in bacteriological culture and/or epidemiological evidence such as forming a single epidemiological unit sharing facilities with other herds where bTB has been confirmed) are then subjected to follow-up tests using the SIT and IFN-γ tests that must be conducted within the following 2-6 months until they recover the OTF status (two consecutive negative herd tests separated by at least 60 days; Supplementary Figure 1). During the study period, two versions of the IFN-γ assay (using bovine and avian PPDs, CZ Veterinaria) were authorized for application in the region according to the Spanish National Eradication Program: during 2010-2016, the Bovigam IFN-γ test (cutoff value of 0.05, Bovigam R , Thermo Fisher Scientific, Waltham, MA, USA), and in 2016-2017, the IDvet IFN-γ test (manufacturer-recommended cutoff value of 35, ID Screen R Ruminant IFN-γ, IDvet, Grabels, France) (21,22). Additionally, pre-movement tests within 30 days prior to the animal movement are routinely performed in all herds (4). When bTB is suspected in a herd (i.e., when reactors are found), the official veterinary services conduct an epidemiological investigation in which they collect information about the herds where the reactor(s) resided and its movements in the 2 years prior to bTB confirmation, the neighboring herds, and other possible sources of infection.

Data Sources
The primary population was comprised of all cattle herds in Castilla y Leon subjected to bTB testing during 2010-2017. Information on herd characteristics, cattle movements, and bTB status of farms was collected through the SITRAN Information System (23). Demographic information available for each herd included its unique identification number, herd production type (beef, fattening, bullfighting, breeding heifers, dairy or mixedbeef, and dairy), and number of animals present at the beginning of each year. Considering the similar management in dairymixed herds and the limited number of fattening units compared to beef herds (see Results), these were grouped into two categories for simplicity (dairy/mixed and beef/fattening, respectively). Additionally, data about the number of movements (contacts) and animals received by each farm during 2007-2016 were also available. Finally, information on the date and type of bTB test (routine testing, pre-movement test, or follow-up tests in bTB-positive farms), number of animals tested, and number of reactors found in the SIT and/or IFN-γ assay for herd tests performed in the frame of the bTB eradication program during 2007-2018 was collected. For herds subjected to whole-herd depopulation, information on the date of depopulation was also obtained. The total number of cattle older than 6 weeks tested per herd in routine and follow-up tests was used to calculate herd size.

Study Definitions and Duration of Bovine Tuberculosis Episodes
A bTB-positive herd test was defined as a herd test with at least one positive animal in the skin test, IFN-γ assay, or bacteriology. When infection was detected through passive surveillance at the slaughterhouse, this was also considered a bTB-positive event. Two bTB-positive herd tests were considered related if they were separated by <18 months regardless of the number of negative herd tests that could take place in the meantime (Figure 1).
Eighteen months was selected as a conservative threshold to increase the power to detect epidemiologically related positive tests since it was close to the median bTB breakdown duration in herds with ≥3 bTB-positive herd tests over the study period (see Results). A bTB breakdown was formed by all related bTBpositive herd tests. The bTB breakdown duration was defined as the period elapsed from the first bTB-positive herd test until the first negative herd test following the last related bTB-positive herd test (Figure 1). A bTB breakdown was defined as resolved when the herd did not experience subsequent bTB-positive tests within the next 18 months. For herds experiencing more than one bTB breakdown during the study period, only the longest bTB breakdown was kept in the database. Depopulated herds during a bTB breakdown were excluded to calculate the median bTB breakdown duration.
Available explanatory variables for each bTB breakdown were province of the herd, bTB prevalence in the county the year before the start of the bTB breakdown (extracted from official records provided by the regional veterinary services and available for the period 2011-2017), herd type, median herd size in the year the bTB breakdown started, relative change in herd size prior to the disclosing bTB-positive herd test (percentage difference between the herd size in the first positive herd test of the bTB breakdown and the median herd size the previous year), number of herds sending animals to the bTB breakdown herd (in-degree) and number of incoming animals in the previous 3 years (for breakdowns occurring in 2010-2017), number of days between the disclosing and follow-up tests, and number and proportion of positive animals to the SIT and bacteriology in the disclosing test. For herds in which bTB was disclosed through pre-movement tests or abattoir inspection, the number of animals tested in the previous routine test was used to calculate the relative change in herd size at the disclosing test. Models considering certain risk factors (number of incoming contacts and animals in the 3 years prior to the start of the bTB breakdown and countylevel herd prevalence) did not include breakdowns starting before 2010/2011 as no information for those factors was available. Herd size was categorized into three categories based on terciles (small = ≤41.1 animals, medium = >41.1-<134 animals, and large = ≥134 animals), relative change in herd size at the disclosing bTB-positive herd test, number and proportion of SIT reactors in the disclosing test, number and proportion of bacteriologyconfirmed animals in the disclosing test, in-degree in the 3 years prior to the start of the bTB breakdown, and number of days between the disclosing and follow-up tests were categorized into quartiles, while county-level herd prevalence was categorized into quintiles.

Survival Analysis
The association between each available predictor variable and the outcome variable "duration of the bTB breakdown" was explored in a survival analysis. The survival distributions for the categories of each of the considered variables were plotted and compared using the log-rank test. These analyses accounted for left truncation (i.e., incomplete information for the start date, for herds with bTB breakdowns already started before January 2010) and right censoring (i.e., incomplete information for the end date, for herds with ongoing bTB breakdowns at the end of the study period, December 2017) as seen elsewhere (24). Here, the included truncation time was days between the bTB disclosure date (in 2007-2009) and the first test in 2010. Depopulated herds were also considered right-censored observations. Univariable Cox proportional hazards models were fitted to compare the hazard ratios of resolving an outbreak for herds depending on the covariates Z 1 . . . Z k , so that: were t represents the survival time, λ 0 (t) is the baseline hazard function, and β i is the regression coefficient of covariate Z i (25). Variables considered in the analyses were production type, number of SIT reactors, and number of confirmed animals through bacteriology in the disclosing test, herd size, relative change in herd size, location at a province level, in-degree, number of incoming animals, and county-level herd prevalence. Variables with p ≤ 0.2 were then considered in a multivariable Cox proportional hazards model. Categorical variables were preferred over continuous based on Akaike's information criteria (AIC) scores. The choice of variable to include in the model for correlated variables (in-degree and number of incoming animals) was based on the AIC (26). The assumption of proportional hazards was evaluated by computing Schoenfeld residuals for each of the study variables. For variables that did not satisfy the proportional hazards assumption (i.e., the true hazard ratio does not change over time), we introduced time-dependent effects β i (t) using a parametric continuous function to construct timedependent covariates (27,28).

Case-Control Study
A case control study to identify risk factors associated with herds with an increased duration of bTB breakdowns ("chronic herds") was also carried out. For this, we defined a case as any herd with a bTB breakdown duration equal to or greater than 784 days. This 784-day threshold was selected as our aim was to identify bTB breakdowns with an unusually long duration for Castilla y Leon, and this included only the top 25% longer resolved bTB breakdown of herds with ≥3 bTB-positive herd tests (see Results). A control was any herd with a bTB breakdown duration less than or equal to the median bTB breakdown duration in the whole study population, i.e., irrespective of the number of positive herd tests (133 days, see Results) as previously performed (10), that was resolved without recurring to depopulation. In this case, calculation of the bTB breakdown duration for breakdowns starting in 2007-2009 or finishing in 2018 included the time in those years. For case herds that underwent depopulation during the bTB breakdown, the date of the clearance/end of the bTB breakdown was set to the date of the depopulation. When herds tested positive following depopulation, this was considered a new bTB breakdown.
Each of the potential risk factors listed above was then tested in a univariable logistic regression model using the chronic status (case/control) as the outcome variable. Risk factors that were significant in the univariable model at a liberal p < 0.20 were considered for inclusion in a multivariable model. Multicollinearity between potential covariables was assessed using the variance inflation factor (VIF) to ensure a mean VIF of <5 among the variables (29) before being offered to the multivariable model. For correlated variables, AIC was used again to perform variable selection (26). The final model considered the selected risk factors along with significant two-way biologically plausible interactions and was built using a backward selection procedure based on a likelihood ratio test (p ≥ 0.05). Results in the model were expressed as odds ratios (ORs) and 95% confidence intervals (CIs). The Hosmer-Lemeshow statistic was used to test the goodness of fit of the model (30).
Statistical analyses were performed using R version 3.5.0 (31). Survival analyses and Cox proportional hazards models were built using the survival (32) and survminer (33) packages, and model fitting assessment was performed through the ResourceSelection (34) package in R.

RESULTS
According to the definition used in this study, there were 3,550 (20%) bTB-positive herds out of 17,793 tested herds in Castilla y Leon at some point during the period 2010-2017. Beef was the predominant production type among the bTB-positive herds (3,026/3,550, 85.2%), followed by dairy (274/3,550, 7.7%), mixed (beef and dairy) (148/3,550, 4.2%), bullfighting (63/3,550, 1.8%), fattening (38/3,550, 1.1%), and raising heifer herds (1/3,550, <0.1%, excluded from further analyses so that the total number of positive herds considered is 3,549). Throughout the study period, the majority of bTB-positive herds was found in the province of Salamanca (n = 1,452/3,550, 40.9%), followed by Avila (n = 603/3,550, 17%), Leon (n = 348, 9.8%), and Segovia (n = 299, 8.4%), whereas the rest of the provinces accounted for the remaining 23.9% (n = 848) of the total number of bTBpositive herds in Castilla y Leon (Supplementary Figure 2). Out of these bTB-positive herds, 49.5% (n = 1,758) experienced only one bTB-positive herd test, 18.1% (n = 642) two, and 32.4% (n = 1,149) of the herds were positive in three or more herd tests. Four hundred nine herds tested negative at least once between bTB-positive tests within the same breakdown. The median bTB breakdown duration in herds with ≥3 bTB-positive herd tests FIGURE 2 | Kaplan-Meier survival estimates of bovine tuberculosis (bTB) breakdown duration by production type. Crosses indicate censored observations. Number of censoring table shows the number of herds that did not clear the infection during the study period or were depopulated during the study period at regular time intervals.
after the exclusion of 49 herds that were subjected to wholeherd depopulation during the bTB breakdown was 511 days [interquartile range (IQR) = 284-784].

Survival Analysis
Ninety-five out of the 3,549 bTB-positive herds were excluded from the survival analyses as they were tested only once during the study period. Up to 321 (9.3%) out of the remaining 3,454 bTB-positive herds detected in the period 2010-2017 did not clear the infection during the study period or were depopulated and were therefore considered censored observations in the analysis. In addition, 247 (7.2% out of 3,454 herds) of the bTB breakdowns had started before 2010 and were therefore left-truncated.
Kaplan-Meier curves for the 10 predictor variables are shown in Figure 2 and Supplementary Figures 3-5. Results of the univariable analyses are shown in Table 1. No data on the number of incoming contacts (and animals) in the 3 years prior to the start of the bTB breakdown and county-level herd prevalence was available for breakdowns starting before 2010 (258 herds) and 2011 (437 herds), respectively ( Table 1). Additionally, no information about relative change in herd size was found for 83/3,454 herds, as no routine tests were performed prior to the disclosing test. For the remaining 3,371 herds, relative change in herd size (median = 2.3%, IQR = −7.2-14.3) was not associated with the length of bTB breakdown duration ( Table 1 and Supplementary Figure 4). The number of SIT and bacteriology-positive animals in the disclosing test were retained in the models over the proportion of the herd positive to SIT Frontiers in Veterinary Science | www.frontiersin.org  and bacteriology based on better AIC. All the remaining variables (increasing herd size, in-degree, number of incoming animals in the 3 years prior to the start of the breakdown, county-level herd prevalence, number of days between the disclosing and followup tests, production type, province, number of SIT reactors and positive to bacteriology animals in the disclosing test) were significantly (p < 0.001) associated with an increasing time to recover OTF status ( Table 1 and Supplementary Figures 3-5). However, the last five variables had non-constant Schoenfeld residuals across time with a large departure of the proportional hazards assumption (Supplementary Figure 6).
The multivariable Cox regression model, run on 3,001 herds with complete information, contained all predictor variables except relative change in herd size (%) and in-degree. The final Cox model was fitted with both time-fixed (production type and number of SIT reactors in the disclosing test) and time-varying coefficients (province, number of days between the disclosing and follow-up tests, and number of positive to bacteriology in the disclosing test), and obtained coefficients did not differ substantially (<19% change) compared with those estimated in the univariable models except for production type, province, and number of SIT reactors in the disclosing test (Tables 1, 2). Time to recover OTF status was not significantly shorter in dairy than in beef/fattening (and bullfighting) herds (p = 0.555, HR = 1.04, 95% CI 0.9-1.2; Table 2). The probability of resolving the bTB breakdown over time was lower in herds with no reactors in the disclosing test (thus detected by other means such as passive surveillance) and herds with ≥5 reactors compared with those with just one reactor ( Table 2). Province also affected significantly the expected herd breakdown duration ( Table 2). Time required to recover OTF status was longer with increasing herd size, number of incoming animals (although not linearly), number of days between the disclosing and follow-up tests, and county-level herd prevalence. There was also a trend of increasing hazard of longer bTB breakdowns associated with the number of bacteriology-positive animals in the disclosing test, although differences with the baseline category (no reactors) decreased over time ( Table 2).

Case-Control Study
A total of 347 herds with a bTB breakdown duration equal to or greater than 784 days (cases) and 1,443 herds suffering a bTB breakdown with a duration ≤133 days (controls) were considered in this analysis. No herd had more than one chronic bTB breakdown. Spatial distribution of these chronic and nonchronic herds in Castilla y Leon from 2010 to 2017 is represented in Figure 3. Valladolid was the province with the highest proportion of chronically infected herds among all its herds included in this analysis (11/26, 42.3%), followed by Soria (19/71, 26.8%) (Figure 3 and Table 3). Fourteen case herds underwent depopulation during the bTB breakdown. Forty-seven case herds with ongoing bTB breakdowns at the start of the study period (January 2010) and 30 case herds with ongoing bTB breakdowns after its end (December 2017) were included in the study since the breakdown they experienced already classified them as a case herd.
The majority (>80%) of both case and control herds were beef/fattening ( Table 3). Overall median herd size in the year

P-value
Median herd size in the year the bTB breakdown started (n = 3,447) SIT, single intradermal test. The proportion of herds with missing information for any given variable was always <8%; out of the 1,790 herds considered in the case-control study, no information for number of incoming contacts (and animals) in the 3 years prior to the start of the bTB breakdown and county-level herd prevalence was available for breakdowns starting before 2010 (51 herds) and 2011 (130 herds), respectively, and were subsequently removed from models considering these variables (Table 3). Moreover, no information about relative change in herd size was found for 49 herds, as no routine tests were performed prior to the disclosing test. In the univariable analysis, nine out of the 10 evaluated variables (all except relative change in herd size at the disclosing bTB-positive test; Table 3) were potentially associated with experiencing a chronic bTB infection. In-degree and number of incoming animals were highly correlated (ρ = 0.74, p < 0.001), and the latter was selected for the multivariable model based on better AIC. According to the final multivariable model including 1,657 herds with complete information on all covariates considered, herd size, number of incoming animals in the 3 years before the outbreak, county-level prevalence before the outbreak, the province where the herd was located, number of SIT-positive and bacteriology-confirmed animals in the disclosing test, and an increasing number of days between the disclosing and follow-up tests were all associated with being a chronically infected bTB herd ( Table 3), while production type was not. Infected herds with a higher probability of suffering a chronic breakdown were larger (≥134 animals, OR = 9.45 CI = 4.8-19.8), located in Soria and Valladolid (OR = 4.73, CI = 1.9-11.7 and OR = 5.32, CI = 1.5-17.4, respectively), had either Frontiers in Veterinary Science | www.frontiersin.org no (OR = 1.7, CI = 1.1-3) or above 5 (OR = 2.75, CI = 1.2-5.9) SIT reactors in the disclosing test and at least three (OR = 89.96, CI = 28.1-408) animals confirmed through bacteriology, and were located in counties with a herd prevalence above 0.2% (OR = 2.78, CI = 1-7.4) in the year preceding the start of the bTB breakdown (Table 3). There were no significant interaction terms. The model fitted well the data (Hosmer-Lemeshow test, χ 2 = 10.29, p = 0.245).

DISCUSSION
The persistence of bTB in certain cattle herds in terms of either herd recurrence or prolonged periods of restriction is a major problem in countries pursuing disease eradication worldwide. Among the factors that may substantially extend the time to recover OTF status and thus hamper eradication programs are the chronic nature of the disease, the presence of wildlife reservoirs, and the lack of performance of diagnostic tools. Several studies have conducted risk factor analyses for chronic bTB herd breakdowns in Europe (10,(35)(36)(37)(38), but factors associated with disease persistence in infected herds had not been characterized at this level yet in Spain. Survival and case-control analyses performed in this study revealed the impact of both farm and breakdown characteristics on the probability of disease persistence in infected farms.
In this study, we defined two positive bTB herd tests as potentially related if they were separated by no more than 18 months, the median duration of bTB outbreaks in herds with ≥3 bTB-positive herd tests, which represented 32.4% of the total bTB-positive herds during 2010-2017. Herds with ≤2 bTBpositive herd tests were not considered to select the threshold in order to focus on problematic herds. The median bTB breakdown for all bTB-positive herds in Castilla y Leon duration was 133 days, which is in agreement with values reported in Northern Ireland (10,36).
The univariable survival analyses revealed that production type was associated with longer bTB breakdown durations: unsurprisingly, outbreaks in bullfighting herds were significantly longer compared to beef/fattening and dairy/mixed herds in agreement with previous studies conducted in a different region in Spain (39). Some tests may have differing performances depending on breed type, what could impact the time to recover OTF status (40): in this sense, both the skin test and the IFNγ assay can have a lower performance in bullfighting cattle due to a combination of stress, difficulties for performing correctly the test due to the temper of the animals, and the lower response to tuberculins/decreased IFN-γ production, what could lead to a higher proportion of false-negative animals remaining in the herd and thus further contribute to the increased outbreak duration (15,41). In addition, bullfighting herds are extensively managed and therefore can have an increased risk of contact with infected wildlife (42), which could lead to frequent reinfections, thus increasing the length of the outbreak. However, when other variables were considered in the multivariable model, differences were not significant. The absence of a significant effect of the bullfighting production type in outbreak duration may be due to the low sample size of bullfighting herds (<2% of the total herds included in the survival analysis). Interestingly, herd type was identified as a risk factor only in the univariable models (survival analysis and case-control study) but not in the multivariable, thus suggesting that other factors (potentially associated with production type such as province, county-level herd prevalence, and herd size) could be explaining at least part of the apparent effect of this variable in the probability of a herd being classified as a case.
The number of SIT reactors at the disclosing test was significantly associated with both outbreak duration and the odds of a chronic infection: as expected, a higher number of reactors (≥5), suggestive of active circulation of the disease in the herd, was associated with increased duration/risk of being a case herd compared with the reference category (one reactor). However, herds in which no SIT reactors had been found in the disclosing test were also at higher risk (Tables 1, 3) compared to herds with one SIT reactor. These herds (n = 777) included herds in which IFN-γ was being used (n = 321, 41.3%), what indicates that even in the absence of SIT reactors, there was conclusive evidence of the presence of disease, since the IFN-γ test is only used on units confirmed by bacteriology or based on epidemiological grounds (i.e., multiple herds managed as a single epidemiological unit with one or more having a confirmed bTB infection or that had animals that were in shared pastures with other positive herds) (Supplementary Figure 1). During 2016, both the Bovigam R and the IDvet test versions of the IFNγ were applied. However, the potential effect of the different versions of the IFN-γ test was not considered in this study, as no information on which version of the IFN-γ assay was used on specific herds was available. In the remaining 456 herds, IFN-γ was not being applied, and infection was therefore found through slaughterhouse surveillance (detection of lesions) or, for OTF herds with a previous history of bTB, bacteriology performed routinely on older animals sent to the abattoir (performed in <1% of the herds of this study). In this case, therefore, the higher duration of breakdowns and increased risk of becoming a chronic herd could be associated with the presence of anergic animals with visible lesions, especially in extensively managed herds, that could be infecting other animals while remaining undetected for some time during the breakdown. When herds in which IFN-γ was applied were removed, the same result was obtained (change in coefficient <9%, data not shown), thus suggesting that herds found through passive surveillance are in fact more prone to become problematic from an eradication standpoint. A median of 5.9%  breakdowns per year were detected via postmortem analysis (data not shown), although this value is well below the estimates of 27%−46% of all new bTB breakdowns per year reported elsewhere (43,44), results obtained here highlight the usefulness of abattoir surveillance to complement antemortem tests (and partially overcome limitations in their sensitivity). This limitation in the sensitivity of skin tests could also be related to an increased duration of outbreaks due to the difficulties in the removal of all infected animals in one single herd test and would be partially compensated by the repeated application of the tests (45). Additionally, and to maximize the sensitivity of the diagnostic tests, mandatory training courses for veterinarians conducting SITs organized by the Spanish Ministry of Agriculture and Fisheries, Food and Environment have been in place since 2012 (4).
We also found evidence of time-varying effects for the number of bacteriology-positive animals in the disclosing test, province, and number of days between the disclosing and follow-up tests in the survival analyses, as the mean estimates did not distribute evenly throughout the study time. The interval of time between the disclosing test and the follow-up test was significantly associated with both outbreak duration and the odds of a chronic infection, with herds tested at increased intervals (>124 days; Tables 2, 3) being at higher risk compared with herds subjected to follow-up tests within 2 months after the disclosing test. These results are suggestive of the potential benefits of short retesting intervals, which would allow the early removal of positive animals from the herd to reduce the risk of bTB spread.
In the case-control study, cases were selected among herds experiencing bTB outbreaks in the top 25% of the distribution of durations (>784 days, ∼26 months) in order to focus on outlier herds in the Castilla y Leon context. This threshold is well above the 1-year threshold used to differentiate chronic infections arbitrarily selected in previous studies in Ireland (46) and Northern Ireland (10,37). Other studies conducted in the UK used periods of >240 days (twice the minimum restriction period for a confirmed bTB breakdown) (35) and >6 months (38).
The case-control study further revealed the effect of farm-level management and outbreak characteristics with bTB persistence in infected herds: an increased number of bacteriology-positive animals in the disclosing test were significantly associated with increased odds of experiencing long-duration breakdowns. These findings are consistent with a higher risk of chronic infection in herds with multiple SIT reactors, as these herds are more likely to have a higher proportion of animals confirmed through bacteriology. This result may also be related to a higher bacterial excretion in animals with an advanced disease stage and thus leading to a higher risk of bTB persistence due to the increase of the time of contact with susceptible animals (8).
Even though chronic case herds were widespread throughout Castilla y Leon (Figure 3), infected herds located in specific provinces were subjected to longer breakdown durations and a higher probability of becoming chronically infected (Table 3). Interestingly, the two provinces with the highest indication of an increased risk, Soria and Valladolid, hold a relatively low proportion of the population [2.1% (369/17,793) and 1.9% (331/17,793) of the herds considered in the study period, respectively] and had relatively high herd level prevalences (1.5 and 3.7% in 2016, respectively) (47). To date, there is no substantial evidence of disease spillover from wildlife to cattle in these provinces, and further investigations are needed to clarify the reasons for this increased risk.
The association between an increased herd size and the risk of experiencing a chronic outbreak could be explained with the increased odds of finding reactors (both true and false positive) in larger herds (45,48), since herd size is a known risk factor for bTB detection (19,37). According to our study definition of a case herd (experiencing a bTB outbreak with duration ≥784 days), the likelihood of case herds not being truly bTB infected can be considered minimal, and therefore this finding may suggest an increased herd sensitivity of the tests in the eradication program in larger herds, as previously reported in the same region (45). Still, bTB infection was not confirmed through bacteriology in 102/347 herds, some of which could have even recovered the OTF status while experiencing what was defined here as a bTB outbreak. Lack of bacterial isolation in those herds could be related to a very small number of culled animals that would be subjected to postmortem analyses due to a very low level of disease, although the presence of repeated cross-reactions due to non-tuberculous bacteria cannot be ruled out.
The number of incoming animals in the 3 years preceding the start of the bTB breakdown was retained in the case-control logistic model based on AIC and in agreement with findings in the survival analysis, even though no significant differences between the variable categories were observed. Still, several studies have consistently evidenced that cattle movements may be important in bTB transmission (49,50), and even if positive farms were not clustered in the movement network in the region under study here, there was an association between increased connectivity and positivity at the farm level in a previous study (17).
Both survival and case-control studies showed that a countylevel herd prevalence above 0.2% in the year prior to the start of the bTB breakdown was significantly associated with an increased risk of experiencing chronic bTB breakdowns. This result is in agreement with previous studies that identified proximity to infected neighbors as a risk factor for persistent bTB infection in Spain (18) and Ireland (14) or prolonged outbreak duration in Northern Ireland (10,36). Local sources of infection that could contribute to this effect include local movements (51), contact with infected cattle and wildlife reservoirs (52), or environmental sources of M. bovis (53), thus contributing to disease recurrence or local persistence (6,12,54).
The inclusion of these variables (incoming animals in the 3 previous years and county-level herd prevalence in the previous year) in the final models forced the elimination of outbreaks starting before 2011 (n = 437 outbreaks), since no information on them was available for these outbreaks. These outbreaks had been considered to estimate the overall median duration of outbreaks in the region and the threshold for selection of chronic herds. Nevertheless, if these are excluded for these calculations, no major differences were observed (median duration of outbreaks remains at 113 days, and the threshold used to define chronic herds increases from 784 to 844 days), suggesting that they are not significantly affecting the study definitions.
Other factors not explicitly included in the study that could be related to increased duration of outbreaks include wildlife variables or the presence of concurrent infections compromising the diagnostic sensitivity. The effect of paratuberculosis in bTB diagnostic tests has been described in the past in Spain (5,55,56) and elsewhere (57)(58)(59). Similarly, cross-reactivity with environmental saprophytic mycobacteria could also artificially increase breakdown duration due to the presence of unspecific SIT reactions (60,61). However, by selecting herds with ≥3 bTBpositive herd tests in which the disease is typically confirmed through bacteriology (∼74% of all bTB breakdowns in herds with ≥3 bTB-positive herd tests were confirmed through bacteriology, and 72% of the 347 case herds) should have minimized the presence of false-positive herds in our study population. Ideally, the epidemiological relatedness between positive herd tests should be demonstrated by characterizing the M. bovis strains circulating in the farm using molecular typing techniques [spoligotyping, Variable number of Tandem Repeat (VNTR), and whole-genome sequencing (WGS)], but this information was not available for a large proportion of the herds evaluated here.
Time-varying effect emerges when the proportional hazards assumption is not fulfilled (62), and for this reason, we applied a multivariate Cox model with a mixture of time-varying and timeindependent parameters as seen elsewhere (63). This allowed accommodating the time-varying effect observed in two out of the nine variables included in the analyses. In all of them, a decrease in the hazards ratio estimates over time was observed, suggesting that the longer the breakdown lasted, the least important their contribution was. In the survival analyses, left truncation was accounted for 7.2% (out of 3,454) of these herds with already started bTB breakdowns before 2010. However, the fraction truncated was not high enough so that estimates became unstable, as the amount of truncation did not approach or exceed the 50% threshold reported elsewhere (24).
Our findings should be interpreted with caution, as associations found here may not reflect causation. Still, our results demonstrate that certain farm and breakdown characteristics (available at the beginning of the breakdown) can help to predict the probability of experiencing a chronic bTB infection among those already infected, what may help to implement specific measures aimed at preventing reinfection (e.g., increased biosecurity) or decrease residual infection (e.g., early detection and removal to slaughter). Further research is needed to better understand the mechanisms behind persistence of bTB in the region, including the risk factors associated with bTB recurrence and the role of wildlife reservoirs in bTB reinfection and/or maintenance in the herd.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/restrictions: datasets included in this study are owned by the Official Veterinary Services in Castilla y Leon and include information that cannot be made publicly available. Reasonable data requests can be submitted to the authors of the study and will considered individually. Summary data is provided in the article. Requests to access these datasets should be directed to mingonol@jcyl.es.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because data from animals used in this study was generated through the normal functioning of the tuberculosis eradication campaigns in Spain and therefore no additional procedures were performed as a result of this study.

AUTHOR CONTRIBUTIONS
PP and JA designed the study with the help of AG, JN, and OM, who also provided the data. PP analyzed the data with the help of JA, BR, JB, and JLS in data interpretation. PP and JA wrote the initial manuscript and received substantial feedback from all other authors.