Study area and study design

A cross-sectional study was conducted on the Potohar plateau including Islamabad Capital Territory (ICT), Rawalpindi and Attock districts of Pakistan (Fig. 1). The Potohar plateau is a hilly area having a great diversity of fauna and flora. The area is located in north eastern Pakistan with an elevation of 575 m between the northern part of the Punjab and the western part of Azad Kashmir. Rain water is the main source of irrigation of agricultural land. Other parameters related to the sampling sites have previously been described [28].

This area has all major breeds of buffaloes and cattle of Pakistan which are reared under extensive and semi-extensive grazing systems. According to the 2014–2015 provincial livestock population survey, the number of cattle in this area was estimated to be 19.4 million (49% of total cattle in Pakistan), with 22.5 million buffaloes (65% of total buffaloes in Pakistan), providing more than 67% of the total milk produced in the country [11].

Buffaloes and cattle for blood/milk sampling were selected randomly from eight major sampling locations [Ahmadal (Latitude 33°17′ N; Longitude 72°29′ E), Attock (Latitude 33°46′ N; Longitude 72°21′ E), Chak Shahzad (Latitude 33°39′ N; Longitude 73°8′ E), Chauntra (Latitude 33°30′ N; Longitude 72°22′ E), Kahuta (Latitude 33°34′ N; Longitude 73°22′ E), Kallar (Latitude 33°24′ N; Longitude 73°22′ E), Kherimurat (Latitude 33°30′ N; Longitude 72°52′ E) and Rawat (Latitude 33°29′ N; Longitude 73°11′ E)] located in ICT, Rawalpindi and Attock districts from 2009 to 2011 (Fig. 1). A sample size of 202 herds was calculated expecting a herd seroprevalence of 15.6%, a confidence level of 95% and a desired absolute precision (d) of 0.05. Contingencies were taken into account by adding another 25% of animals and herds leading to a total of 253 herds. The 253 herds were randomly selected from the 8 sampling sites due to the lack of a detailed herd and cattle/buffalo identification system. The number of herds was estimated using the formula n = (1.96)² p(1 − p)/d ² [29, 30]. The herds were divided into two categories on the basis of the median value of their sizes; below the median value, the herd was considered as “a small holding cluster” (≤8) and above the median value (≥9) as “a large holding cluster”. Herds were of three types, those having only cattle, only buffaloes and those with both cattle and buffaloes (mixed type). Blood/milk samples were collected from 50% of the animals of a herd, for most small holdings, all animals were sampled. To avoid false positives due to the presence of maternal antibodies, only cattle older than 1 year were sampled.

The questionnaire was distributed paper-based through face-to-face interviews (Additional file 1). Data related to age, sex, urbanity, districts/territory, sampling sites, animal species (cattle or buffalo), abortions in third trimester, metritis, herd size, insemination method, source of replacement of animals and body condition of animals were collected at the sampling day. All data were kept for further assessment or if requested.

Sample collection

A total of 2709 serum samples were randomly collected [1462 buffaloes (53.97%) and 1247 cattle (46.03%)]. Moreover, 2330 milk samples were collected from 1168 cattle and 1162 buffaloes. Approximately 10 ml of blood was collected aseptically from the jugular vein of each animal according to standard procedure [31]. These samples were immediately stored at 4 °C. Samples were then transported to the laboratory. Sera were separated and stored at −20 °C while milk samples were stored at 4 °C.

Serology

Serum samples were initially screened with RBPT antigens (Institute Pourquier, France). Samples positive to RBPT were confirmed with the serum agglutination test (SAT) (Veterinary Research Institute, Pakistan). All serological tests were performed and results were interpreted according to standard procedures [7, 31, 32]. Briefly, 25 µl of serum were mixed with an equal volume of antigen preparation on a glass plate; the plate was agitated gently for 4 min. A serum sample was considered positive if agglutination occurred. A serum sample positive in RBPT as well as in SAT was considered as positive at the animal level.

SAT was carried out with ethylene diamine tetra acetic acid (EDTA) as described previously [32]. The Brucella antigen used in this study was purchased from Immunostics, Inc., USA. One hundred and sixty-eight microliters of Serum Agglutination de Wright (SAW) buffer were added to the first well and 100 μl to the second and third well of a 96-well microtiter plate. 32 µl of test serum was added to the 1st well to reach dilution of 1/6.25. After adequate mixing, 100 μl from the 1st well were transferred to the 2nd well to reach dilution of 1/12.5. Similar to the previous method 100 μl were transferred from the 2nd to the 3rd well to reach dilution of 1/25 and 100 μl discarded from the 3rd well. Then in each well 100 μl of standardized SAW antigen was added giving the serial serum dilutions of 1/12.5, 1/25 and 1/50. The plate contents were thoroughly mixed and incubated for 20–24 h at 37 °C. The value reading was done according to the degree of agglutination [33].

Milk ring test (MRT)

Milk samples were initially screened by MRT. As per manufacturer recommendations, the MRT antigen was kept at room temperature before use. One milliliter milk sample was added to the test tube. Then 30–40 µL of antigen were added, mixed and incubated at 37 °C for 1 h. A sample having a change in color at the top of the milk was considered positive [7, 31].

Isolation and identification of Brucella

Milk samples considered as positive by MRT were used for isolation of Brucella. Isolation was conducted on modified Farrells serum dextrose agar according to standard procedures [31, 34]. Identification and biotyping of these isolates was done according to standard procedures [7, 31, 35].

DNA extraction and qRT-PCR

Serum samples that tested positive in serology were further subjected to DNA extraction using the High Pure PCR Template preparation kit (Roche Diagnostic, Germany). Purity and concentration of DNA was checked by Nano-Drop ND-1000 UV–Vis spectrophotometer (Nano-Drop technologies, USA). DNA samples were stored at −20 °C until further analysis. A Brucella genus-specific (31-kDa salt-extractable immunogenic protein gene, bcsp31) qRT-PCR assay was used for further screening of seropositive samples [36]. Primers and probes were purchased from TIB MOLBIOL (Berlin, Germany). The reactions were conducted in duplicate in microtiter plates (Applied Biosystem, Germany) using M×3000P thermocycler platform (Stratagene, La Jolla, Canada). The thermal profile for assays was 1 cycle of 50 °C for decontamination for 2 min, 1 cycle of 95 °C for initial denaturing for 10 min, 50 cycle of 95 °C for denaturing for 25 s and 1 min for annealing at 57 °C. Cut-off value of cycle threshold (Ct) for a positive sample was ≤40 for Brucella genus specific qRT-PCR being automatically generated by the instrument. Herds with at least one animal positive in qRT-PCR were considered as positive.

Statistical analysis

Se = Se _RBPT  * Se _SAT Whereas the combined specificity is given by Sp = Sp _RBPT  + (1 − Sp _RBPT ) * Sp _SAT .

The data extracted from each selected study included the Se, Sp, total number of subjects considered from which the number of true positives, true negatives, false positives and false negatives were calculated. These data were then analyzed using “metandi” in Stata 12.1 [38]. The outcome of this analysis is a synthesized estimate (and 95% confidence interval) of the sensitivity and specificity of each test adjusted for the total number of subjects in each of the studies included. The true individual animal level prevalence was estimated across the different potential risk/indicator factors.

Screening of the different potential risk factors related to brucellosis seropositivity was done using univariate random effects logistic regression analysis. Sampling location and herd were both used as random effects to interpret potential clustering of animals within herds and for the differences in herd sizes for the animal level analysis whereas only sampling location was used as a random effect for the herd level analysis to account for clustering of herds within sampling sites. This also accounted for the differences in number of animals within herds and number of herds within sampling regions, respectively. Variables with a p value <0.25 in the univariable analysis were further analysed in a multivariable random effects logistic regression model. Manual forward stepwise selection was applied to select the final model using the Akaike’s information criterion (AIC) as the calibrating parameter.

When the removal of a non-significant variable led to a change of more than 25% of the estimated odds ratio, that variable was considered a confounder and was kept in the final model. Multicollinearity was assessed among the independent variables using the Cramer’s phi prime statistic with values >0.7 indicative of co-linearity [39].

All two-way interaction terms of the variables remaining in the final model were assessed for significance based on the likelihood ratio test comparing the model with the desired interaction term and the corresponding model with no interaction terms.

The models were built using the xtmelogit function in STATA, version 12.1, software (SataCorp LP, College station, Texas). Model selection was done using Laplacian approximation and the robustness of the final model was assessed by increasing the number of Quadrature (integration) points and monitoring changes in parameter estimates [41, 42].

True animal and herd level seroprevalence

The results of the meta-analysis yielded an estimated Se and Sp of RBT of 0.93 and 0.98, respectively whereas for SAT the estimated Se and Sp were 0.67 and 0.99, respectively. The overall estimated Se and Sp based on a serial interpretation of both tests were 0.63 and 1.0, respectively. These values were used to compute the adjusted (true) prevalence.

Out of 170 serological positive samples, 89 (52.4%) were positive using the Brucella genus specific qRT-PCR.

Moreover, a total of five isolates were recovered from 156 (6.7%) milk samples of positive animals and these isolates were identified as B. abortus biovar 1 according to standard biotyping procedures.

Risk factors associated with animal and herd level prevalence

The results of the univariable analysis which was based on a random effects model correcting for animal and herd-level clustering indicated that at the individual animal level, sex (cows versus bulls), animal species (cattle or buffaloes) and stock replacement (self-reared versus purchased) were significantly associated with seropositivity at the animal level (p < 0.05) (Table 2). At the herd level, herd size and insemination method were significantly associated with seropositivity at the herd level (p < 0.05) (Table 2). In addition, the animal level factors i.e. district and age and the herd level factors i.e. presence of animals with metritis were not significant at the 5% level but since their p values were <0.25 (Tables 2, 3), they were considered as potential risk factors and thus subjected to the multivariable random effects logistic regression analysis.

Based on the final animal level model, the odds of brucellosis seropositivity were found to be 2.4 (95% CI 1.1–5.5%) times higher among older animals compared to those of the younger animals. In addition, the odds of brucellosis seropositivity were 5.6 (95% CI 2.6–12.0%) times higher among bulls compared to cows. Lastly, animals that were self-reared were 1.6 (95% CI 1.0–2.4%) times more likely to be seropositive for brucellosis compared to those that were purchased. At the herd level, the odds of brucellosis seropositivity were found to be higher for large herds, OR = 5.0 (95% CI 1.7–14.6%) compared to small herds and for herds with occurrence of abortion in third trimester, OR = 17.4 (95% CI 1.4–214.1%) compared to herds with no occurrence of abortion in third trimester, respectively (Table 5). In addition, for herds in which both, AI and natural insemination methods applied, the odds of brucellosis seropositivity were 4.7 times higher compared to those of herds in which natural insemination was practiced. Moreover, the use of AI was found to be a protective factor against brucellosis seropositivity since herds which applied AI were less likely to be brucellosis seropositive, OR = 0.2 (95% CI 0.1–0.8%) compared to those in which only natural insemination was applied.

The value of σ _INT:Herd ² , which is the proportion of the total variance in the model explained by the variance between herds, was 0.63 corresponding to an ICC of 0.16 and thus indicating that 16% of the variability in bovine brucellosis seroprevalence occurrs between herds whereas the rest is due to differences between animals within herds. The ICC for sampling location was close to zero indicating that there were no differences in seropositivity of herds based on their sampling sites.