Epidemiology of Mycobacterium bovis infection in free-ranging rhinoceros in Kruger National Park, South Africa.

Mycobacterium bovis infection, which is a prominent cause of bovine tuberculosis, has been confirmed by mycobacterial culture in African rhinoceros species in Kruger National Park (KNP), South Africa. In this population-based study of the epidemiology of M. bovis in 437 African rhinoceros (Diceros bicornis, Ceratotherium simum), we report an estimated prevalence of 15.4% (95% CI: 10.4 to 21.0%), based on results from mycobacterial culture and an antigen-specific interferon gamma release assay from animals sampled between 2016 and 2020. A significant spatial cluster of cases was detected near the southwestern park border, although infection was widely distributed. Multivariable logistic regression models, including demographic and spatiotemporal variables, showed a significant, increasing probability of M. bovis infection in white rhinoceros based on increased numbers of African buffalo (Syncerus caffer) herds in the vicinity of the rhinoceros sampling location. Since African buffaloes are important maintenance hosts for M. bovis in KNP, spillover of infection from these hosts to white rhinoceros sharing the environment is suspected. There was also a significantly higher proportion of M. bovis infection in black rhinoceros in the early years of the study (2016–2018) than in 2019 and 2020, which coincided with periods of intense drought, although other temporal factors could be implicated. Species of rhinoceros, age, and sex were not identified as risk factors for M. bovis infection. These study findings provide a foundation for further epidemiological investigation of M. bovis, a multihost pathogen, in a complex ecosystem that includes susceptible species that are threatened and endangered.

African rhinoceros (Diceros bicornis, Ceratotherium simum) are currently under threat due to poaching activity and habitat destruction, as well as the underrecognized threat of infectious diseases (1, 2). Mycobacterium bovis infection has been confirmed in the African rhinoceros population in Kruger National Park (KNP), South Africa (3–6). The discovery of M. bovis infection in this population has led to a quarantine of rhinoceros intended for translocation from the park to other protected areas, which has significant conservation consequences (1). The paucity of knowledge regarding the epidemiology and risk of transmission from infected rhinoceros has become evident when assessing impact on the KNP population and potential for spread to other populations (1).

Because M. bovis infection is chronic and may not cause clinical signs of disease for months to years, its presence in an ecosystem with multiple susceptible host species may not be recognized for decades, as has been documented in several bovine tuberculosis (bTB) afflicted wildlife populations worldwide (7–9). KNP is considered endemic for bTB, with African buffaloes (Syncerus caffer) being the key maintenance hosts (8–10). Historically, M. bovis is believed to have originated from infected cattle adjacent to the park boundaries in the 1960s and 1980s, but was not detected until the 1990s in infected buffalo herds (11). Since then, 15 additional wildlife species in KNP have been documented with infection (12), including rhinoceros, in which infection was confirmed using mycobacterial culture and M. bovis species confirmation (3–6).

The epidemiology of bTB in a complex system that contains multiple hosts with varying susceptibility results in an array of opportunities for infection spread. Black and white rhinoceros in KNP share environmental resources (including browse/grazing, and water sources) with potentially M. bovis–infected African buffaloes (8–10), greater kudu (Tragelaphus strepsiceros) (13, 14), warthogs (Phacochoerus africanus) (15), and other species (12). A recent review of potential scenarios for interspecies transmission has suggested that rhinoceros may become infected with M. bovis in ecosystems containing other infected hosts (1). Interspecies spread has been demonstrated in other systems, including badgers and cattle in the United Kingdom (16), wild boars, deer, and cattle in Spain (17), and deer and cattle in the United States (18). Mechanisms of transmission between herbivores are poorly understood but have been attributed to indirect interaction through shared resources such as pastures, feed, or water holes that are contaminated by M. bovis–shedding hosts (17, 19, 20). The potential for intraspecies transmission between rhinoceros is also plausible based on antemortem detection of mycobacteria in respiratory secretions (1, 21).

For species (like white and black rhinoceros) that are considered threatened or endangered, the presence of a controlled infectious disease can significantly hamper conservation efforts and potentially impact population health and survival (22–24). Regulations imposed by the Department of Agriculture, Land Reform and Rural Development in KNP, due to the diagnosis of M. bovis in rhinoceros and other species, are an additional barrier to the movement of rhinoceros from the park to other national or private reserves. This can have a significant impact on conservation of the species, as KNP has historically been an important source population of rhinoceros for other conservation strongholds in South Africa and other African countries. In order for a captured rhinoceros to be moved out of KNP, it must first be placed in a quarantine facility for 3 mo, and must test negative for M. bovis infection during repeated testing events.

It is therefore critical to be able to assess infection status in these populations. Recent advances in the development of diagnostic tests (12, 25) for M. bovis infection in wildlife beyond conventional mycobacterial culture and M. bovis species confirmation has allowed for antemortem testing. The QuantiFERON TB Gold In-Tube Plus (QFT) (Qiagen)–interferon gamma release assay (QFT–IGRA) was recently validated for use in white rhinoceros (4), and has been used for testing KNP rhinoceros for M. bovis infection. These results provide the opportunity to generate an understanding of epidemiological determinants and risk factors for infection and disease transmission within the rhinoceros population.

Here, we report a population-based study on the epidemiology of M. bovis infection in free-ranging African rhinoceros. We investigated the distribution of M. bovis infection in rhinoceros over the KNP landscape and identified the demographic, spatial, and temporal factors that may drive infection in this population (1). Our findings begin to uncover the complex epidemiology of bTB for rhinoceros in a multihost system where bTB is endemic. Results from this study emphasize the importance of disease surveillance in managed wildlife systems and support current quarantine and testing requirements for rhinoceros in KNP. These findings are important for preventing the spread of M. bovis infection to other rhinoceros populations (2, 3). In a broader sense, this study reflects the complexity of investigating a multihost pathogen that has been introduced into a previously naïve system. It provides an opportunity to increase awareness of the global impact that TB and other zoonotic pathogens can have on domestic and wild animal populations, food security, and conservation of species and ecosystems.

Results

Prevalence of M. bovis Infection and Spatial Clustering.

The study population consisted of 475 free-ranging African rhinoceros that were opportunistically sampled in KNP from 2016 to 2020, as described in . The M. bovis infection status could be determined for 437 (92%) of 475 rhinoceros, largely based on antemortem test results from IGRA (n = 428) (4, 5), with a few individuals’ case statuses confirmed using (postmortem) mycobacterial culture of tissues with M. bovis species confirmation using a rapid diagnostic multiplex PCR (RD-PCR) (n = 9) (3, 26) (Table 1). Fifty-eight of the 437 study population individuals with a known infection status were tested after recent translocation out of the park to quarantine areas. Of these, only five were M. bovis positive according to the IGRA result. All five of these individuals were sampled for testing within hours of their translocation out of the park; therefore, their positive status reflects infection acquired in KNP.

Table 1.

Mycobacterial QFT-IGRA test and culture results for 475 African rhinoceros in KNP, South Africa (2016–2020)

M. bovis status	Test result for IGRA and culture*	White rhinoceros	Black rhinoceros	Total
M. bovis positive (n = 82)	IGRA+/no culture performed	55	19	74
	IGRA+/culture+	3	0	3
	IGRA−/culture+	5	0	5
M. bovis negative (n = 355)	IGRA−/no culture performed	253	101	354
	IGRA−/culture-	1	0	1
Unknown (n = 38)	IGRA status undetermined/no culture performed	30	8	38
	Total	347	128	475

IGRAs were completed on plasma from QuantiFERON-stimulated rhinoceros whole blood (4, 5) to detect immune sensitization to mycobacteria. Rhinoceros were classified as IGRA positive, negative, or unknown for M. bovis infection based on criteria outlined in .

*M. bovis was isolated in tissues obtained at necropsy from eight of nine white rhinoceros through conventional mycobacterial culture using BACTECMGIT platform, with M. bovis species confirmation using a rapid diagnostic multiplex PCR (RD-PCR) (26). Culture positive animals included three adult males (one necropsied in 2016, two in 2018), one subadult male (necropsied in 2016), two adult females (necropsied in 2016 and 2018), and two subadult females (necropsied in 2016 and 2018).

The apparent M. bovis prevalence was estimated and adjusted according to the sensitivity and specificity of the IGRA assay, as described in (27). The overall adjusted prevalence (based on the IGRA) during the study period was 15.4% (82/437: 95% CI: 10.4 to 21.0%). The majority of M. bovis–positive cases were considered clinically normal (83%, 68/82), based on veterinary clinical assessment at the time of immobilization. Species-specific prevalence was 17.0% (63/317; 95% CI: 11.0 to 23.9%) for white rhinoceros and 11.2% (19/120; 95% CI: 3.1 to 22.2%) for black rhinoceros, and were not statistically different (Fisher’s exact P = 0.41). The 38 rhinoceros that were classified as having an “unknown” infection status were excluded from prevalence calculations and further analyses. Further description of demographic characteristics of the study population, according to M. bovis status, species, sex, and age, is shown in .

Prevalence of M. bovis infection (adjusted based on IGRA sensitivity and specificity) in rhinoceros according to ranger area and ecozone is shown in Fig. 1. In total, 420 individuals were included in prevalence calculations across the different areas. For ranger area, the highest M. bovis prevalence was in Pretoriuskop (28.1%), and the lowest was in Tshokwane (9.7%; Fig. 1); however, no significant differences in prevalence across ranger areas were identified (Fisher’s exact P > 0.05 for all comparisons). Similarly, the highest prevalence by ecozone was in the Pretoriuskop Sourveld (36.4%), and the lowest prevalence was in the Lebombo Mountain Bushveld (13.3%; Fig. 1). The Pretoriuskop Sourveld ecozone had a significantly higher prevalence of M. bovis infection than the Sabie/Crocodile Thorn Thickets (P < 0.001), the Mixed Bushwillow Woodlands (P = 0.02), and the Lebombo Mountain Bushveld (P = 0.01) ecozones. No other significant differences were detected (Fisher’s exact P > 0.05 for all other area comparisons). Ecozone and ranger area were not evaluated in the univariate and multivariable analyses, due to the numerous categories and the potential for missing covariate patterns.

Fig. 1.

Prevalence (percent) of M. bovis infection in rhinoceros in KNP, South Africa, 2016–2020 (n = 420). Prevalence estimates are reported in the key in descending order. Areas with insufficient data (n ≤10 sampled animals) are shown in white. (A) The prevalence of M. bovis in the study population according to ranger management area. No significant differences in M. bovis prevalence according to ranger area were identified (Fisher’s exact P > 0.05 for all comparisons). A single, statistically significant cluster of radius 6.5 km is depicted by a hatched yellow circle based on Kulldorff’s spatial scan statistic (28); 12 cases of M. bovis were identified in this cluster, whereas the model predicted only 4 (relative risk = 3.5, P = 0.036). (B) The prevalence of M. bovis according to ecozone (76). Significant differences in M. bovis infection prevalence were detected between the ecozones with the same superscript letter. Prevalence in the Pretoriuskop Sourveld ecozone was significantly different from in Sabie/Crocodile Thorn Thickets (P < 0.001), Mixed Bushwillow Woodlands (P = 0.02), and Lebombo Mountain Bushveld (P = 0.01). Fisher’s exact P > 0.05 for all other area comparisons.

Further exploration of differences in geographical distribution of M. bovis infection using Kulldorff’s spatial scan statistic (28) showed significant spatial clustering of infection in white rhinoceros. A single, statistically significant cluster of radius 6.5 km was detected toward the northern border of the Pretoriuskop ranger area (Fig. 1). Twelve cases of M. bovis were identified in this cluster, whereas the model predicted only four (relative risk = 3.5, P = 0.036). No other significant spatial clustering was detected.

Univariate Analyses.

Since environmental risk factors were hypothesized to be similar for black and white rhinoceros due to sharing of resources, the initial univariate analyses were performed with data from the two rhinoceros species (black rhinoceros, or D. bicornis, and white rhinoceros, or C. simum) combined. In total, 13 variables were evaluated in the univariate analysis using logistic regression (as outlined in Tables 2 and 3). This included three spatial variables: number of buffalo herds, number of kudu herds, and buffalo density. Because the spatial scale for these variables was unknown, we used our data from repeated captures to estimate a plausible spatial scale over which rhinoceros in the study population could potentially move. We then evaluated the three spatial variables across different potential circular home ranges, with radii of 5.75, 11.5, 17.25, and 23 km. For each variable, a single home range size was ultimately chosen for further statistical evaluation based on the best-fitting single-variable logistic regression model. Chosen home ranges included 17.25 km for number of buffalo herds and 5.75 km for both number of kudu herds and buffalo density (Table 3). show all results from evaluations using home ranges with various distances; methodology is described in detail in .

Table 2.

Frequency distributions and univariate analyses of potential risk factors (measured as categorical variables) for M. bovis infection in African rhinoceroses in KNP, South Africa, 2016–2020 (n = 437)

Risk factor	Number of M. bovis–positive rhinoceros n = 82 (percent of total in category)	Number of M. bovis–negative rhinoceros n = 355 (percent of total in category)	OR (95% CI)	P
Species				0.33*
White rhinoceros (C. simum)	63 (77)	254 (72)	1.3 (0.8–2.3)
Black rhinoceros (D. bicornis)	19 (23)	101 (28)	Reference
Sex				0.82*
Female	48 (59)	203 (57)	1.1 (0.6–1.7)
Male	34 (41)	152 (43)	Reference
Age				0.27*
Adult	52 (63)	192 (54)	1.7 (0.8–3.6)
Subadult	21 (26)	107 (30)	1.2 (0.5–2.8)
Calf	9 (11)	56 (16)	Reference
Orphan status (calves only, n = 65)				0.61
Orphaned	2 (22)	17 (30)	0.7 (0.1–3.5)
With mother	7 (78)	39 (70)	Reference
Health status				0.38
Injured/abnormal health	14 (17)	47 (13)	1.3 (0.7–2.6)
Appear healthy	68 (83)	307 (87)
Sampling year				0.02*,†
2016	10 (12)	19 (5)	4.5 (1.7–12.1)
2017	15 (18)	42 (12)	3.1 (1.3–7.2)
2018	13 (16)	59 (16)	1.9 (0.8–4.5)
2019	33 (40)	141 (40)	2.0 (1.0–4.2)
2020	11 (13)	94 (26)	Reference
Season				0.61
Dry	43 (52)	197 (55)	0.9 (0.5–1.4)
Wet	39 (48)	158 (45)	Reference
Nearest permanent water source type‡				0.93*
Waterhole	36 (44)	156 (44)	1.0 (0.6–1.6)
River	46 (56)	195 (56)	Reference
Number of kudu herds nearby‡,§,¶				0.30
8–14	24 (29)	95 (27)	1.1 (0.5–2.3)
5–7	23 (28)	132 (37)	0.7 (0.4–1.6)
3–4	21 (26)	64 (18)	1.4 (0.7–3.0)
0–2	14 (17)	60 (17)	Reference

*Selected for inclusion in multivariable model.

†Met screening criteria.

‡Here n = 433, and excludes individuals with unknown capture locations.

§Within home range of 5.75 km radius.

¶Categories created according to quartile distribution of measure over the study population.

Table 3.

Frequency distributions and univariate analyses of potential risk factors (measured as continuous variables) for M. bovis infection in African rhinoceroses in KNP, South Africa, 2016–2020 (n = 437)

Continuous variables	Median for M. bovis–positive rhinoceros n = 82 (IQR)	Median for M. bovis–negative rhinoceros n = 355 (IQR)	OR (95% CI)	P
Distance to nearest water source (km)*	2.34 (1.33–3.62)	2.41 (1.32–4.08)	0.9 (0.8–1.0)	0.15†,‡
Distance to nearest M. bovis infected rhinoceros (km)*,§,¶	2.76 (1.00–4.81)	2.54 (1.35–4.89)	0.8 (0.6–1.1)	0.11†,‡
Number of buffalo herds nearby*,§,#	86 (46–111)	64 (41–102)	1.8 (1.2–2.8)	0.006†,‡
Surrounding buffalo density (estimated buffalo per square kilometer)*,§,¶,ǁ	1.21 (0.55–2.30)	1.13 (0.67–2.60)	0.9 (0.6–1.3)	0.41

IQR, interquartile range.

*Here n = 433, and excludes individuals with unknown capture locations.

†Selected for inclusion in multivariable model.

‡Met screening criteria.

§Odds ratio and CI calculated with a log transformation of the associated variable.

¶Added one before log transformation of the measured variable for the study population as variable is equal to zero for at least one of the included individuals.

#Within home range of 17.25 km radius.

ǁWithin home range of 5.75 km radius.

Four of the 13 variables met screening criteria (P < 0.2; Tables 2 and 3). Two of these factors were statistically significant (P < 0.05) in the univariate analyses, including sampling year and number of buffalo herds within the rhinoceros home range (circular buffer; radius = 17.25 km). Rhinoceros species, age group, sex, and nearest permanent water source type were also considered important demographic covariates or plausible effect modifiers and were evaluated in the final multivariable model(s).

Even though 82 of the 437 study individuals tested positive for M. bovis infection, no significant association between M. bovis infection and (apparent) clinical health status was detected. The majority of the test-positive individuals appeared clinically normal (83% of infected animals were clinically normal, 68/82), and all but 2/14 M. bovis–positive individuals with recorded clinical abnormalities had afflictions that were associated with poaching or fighting injuries, rather than evidence of infection.

Multivariable Analyses.

The final, combined-species model with all rhinoceros included 433 individuals, and consisted of five variables, namely, species, sex, age, sampling year, and number of buffalo herds within a 17.25-km radius of the capture location ().

This model indicated that the year of sampling and the number of buffalo herds within a 17.25-km radius of the capture location were significantly associated with M. bovis infection in rhinoceros. Specifically, rhinoceros sampled in years 2016 (odds ratio [OR] = 4.4; 95% CI: 1.6 to 12.3), 2017 (OR = 3.4; 95% CI: 1.4 to 8.1), and 2019 (OR = 2.2; 95% CI: 1.0 to 4.6) had higher odds of infection compared to the reference year 2020 (P = 0.01). Additionally, for each additional log-transformed buffalo herd in the rhinoceros home range, the odds of M. bovis infection increased by 75% (OR = 1.75; 95% CI: 1.1 to 2.8). However, there was significant effect modification identified across the sampling year by species (). Therefore, we also constructed species-specific models with the same variables. The final species-specific models are reported in Fig. 2, and in greater detail (including model fit parameters) in . Additional effect modification by species or other factors were not identified, and other variable combinations did not improve the fit of the model or indicate additional sources of confounding.

Fig. 2.

Forest plots depicting species-specific multivariable models of factors associated with M. bovis infection African rhinoceros in KNP, South Africa (2016–2020). Parameters for both models, including coefficients, SEs, and fit statistics, are reported in detail in .

Species-Specific Models.

Variables found to be significantly associated with M. bovis infection differed between white and black rhinoceros. For the white rhinoceros, each additional log-transformed buffalo herd in the home range resulted in an increase in odds of infection by 77% (OR = 1.77; 95% CI: 1.07 to 2.92, P = 0.02; Fig. 2 and ). This corresponds to a probability of M. bovis infection of 0.057 when the number of buffalo herds in the white rhinoceros home range is at the minimum (number of buffalo herds = 6) and all other factors in the model are held at their mean. This probability increases to 0.192 when the number of buffalo herds is at the median (number of buffalo herds = 66) and all other factors are held at their mean.

Importantly, the numbers of buffalo herds should be considered a relative, rather than an absolute, measure of exposure to buffaloes, since we do not have a precise measure of exposure to buffalo herds for each rhinoceros; however, our data support the hypothesis that white rhinoceros in areas with more buffalo herds are at an increased risk of M. bovis infection compared to those in areas with fewer buffalo herds, while controlling for other factors. Sampling year was not significantly associated (P = 0.15) with M. bovis infection in white rhinoceros.

Conversely, for the black rhinoceros, the year of sampling was significantly associated with M. bovis infection (P = 0.01; Fig. 2 and ), while controlling for other factors in the model; individuals sampled in years 2016 (OR = 18.11; 95% CI: 2.09 to 157.15), 2017 (OR = 4.36; 95% CI: 0.98 to 19.41), and 2018 (OR = 4.60; 95% CI: 1.11 to 19.05), compared to years 2019 and 2020 (note that the years 2019 and 2020 were combined for black rhinoceros, due to small numbers of animals in those categories). However, the number of buffalo herds nearby (within 17.25 km) was not significantly associated (P = 0.68) with M. bovis infection in the black rhinoceros (Fig. 2 and ).

ORs from final models reporting adjusted associations for M. bovis infection among white rhinoceros and black rhinoceros, separately, are reported graphically in Fig. 2. Species-specific model estimates, including coefficients and SEs, are included in .

Discussion

A considerable and widespread M. bovis infection burden was reported for the KNP rhinoceros population (15.4%), with similar rates of infection found in males and females of all age groups and in both black and white rhinoceros. Although demographic factors were not associated with risk, an increasing number of buffalo herds in the white rhinoceros home range, and year of sampling in black rhinoceros, increased the risk of M. bovis infection in this population. The KNP rhinoceros are central to the “Integrated Strategic Management of Rhinoceros” plan introduced by the South African Department of Environmental Affairs (29, 30). This strategy relies, in part, on translocation of individuals from the KNP population to newly developed safeguarding strongholds around the country. Therefore, the findings in this study support the decision to impose quarantine (31) on all rhinoceros (regardless of demographics) prior to translocation, in order to mitigate the risk for inadvertent M. bovis spread to other ecosystems outside KNP.

The distribution of M. bovis infection in KNP rhinoceros is similar to that reported for other species in the park. For example, a 1991–1992 survey of bTB in 1,122 African buffaloes in KNP showed widespread bTB in the central and southern regions of the park (including Houtboschrand and regions to the south of it), with individual herd bTB prevalence up to 67% (10). A later study showed spread of M. bovis infection to African lions (n = 70) sampled in 2012/2013 in the same areas of KNP, with an overall infection prevalence of 44% (32). Such extensive infection is increasingly observed in additional species in KNP, including warthogs (33), African wild dogs (Lycaon pictus) (23), and African elephants (34), with cases identified in more than 15 species in the park, to date (12). Taken together, these findings suggest that spillover of bTB is not a new occurrence and support the need for ongoing bTB surveillance across species to continuously assess disease risk and conservation impact, and to better understand transmission within and from the multihost system in KNP.

The detection of a single statistically significant M. bovis infection cluster, with a 6.5-km radius, in white rhinoceros toward the northern border of the Pretoriuskop ranger area (Fig. 1) is in concordance with the higher infection prevalence in the Pretoriuskop Sourveld ecozone, compared to other areas in the east and the north (Fig. 1). The identified cluster is in close proximity to the KNP border with the surrounding Mpumalanga province. Importantly, the region outside of the southern KNP borders is primarily farmland and home to livestock herds, specifically, cattle. Livestock in areas around the southern border of the park have historically been implicated in spillover of M. bovis to wildlife in KNP, including African buffaloes (10, 35). Now recognized as maintenance hosts for M. bovis in KNP (10, 35–37), African buffaloes share similar vegetation preferences to that of white rhinoceros. Interestingly, favored vegetation is abundant in the Pretoriuskop Sourveld ecozone where the highest M. bovis prevalence in rhinoceros, as well as the only significant case cluster, occurred.

Results from this study suggest a role for buffalo in M. bovis infection, specifically, in white rhinoceros. Adjusted associations showed an increasing risk of M. bovis infection in white rhinoceros with increasing numbers of nearby buffalo herds. This suggests that the African buffalo in KNP may serve as a potential source for spillover of M. bovis infection into white rhinoceros, as observed in other species (1, 10, 14, 38, 39). The importance of buffalo as a predictor of infection in white rhinoceros as compared to black rhinoceros was expected, as the landscape and preferred vegetation of white rhinoceros closely mirror that of buffaloes (40, 41); both are grazing species, and, therefore, white rhinoceros are more likely to share habitat with buffaloes than are black rhinoceros.

Interestingly, the association was only significant when examining buffalo herds—the number of individual buffaloes nearby and buffalo density were not significantly associated with M. bovis infection. Importantly, the crude nature of the available data makes it difficult to define the exact relationship of the buffalo variable with M. bovis infection. It is possible that the identification of a buffalo “herd” is a more precise measure compared to estimating total number of buffalo in a herd from aerial surveys. Alternatively, this finding could simply be due to variation in individual versus herd-level prevalence in KNP buffaloes (37). Our findings could also suggest that effective contact rates between infected buffalo and susceptible rhinoceros that lead to transmission do not depend on buffalo herd size.

Presence of infected buffalo herds was likely to increase environmental contamination with M. bovis due to shedding, consequently increasing exposure of rhinoceros to the pathogen. This phenomenon is observed in other bTB multihost systems, which involve wild boar, red deer, and cattle populations, across continental Europe (17, 19, 20). Mechanisms of transmission between herbivores are unclear but have been attributed to indirect interaction through shared resources such as pastures, feed, or water holes that are contaminated by M. bovis–shedding hosts (19, 20). In the current study, we did not find any associations with distance to nearby water source or water type. Prospective evaluations with refined measures of frequency of individuals at particular water sources may further elucidate exposure to contaminated environments.

Potential transmission of TB between rhinoceros was also plausible (1, 21, 42), but we did not find an association in adjusted models when evaluating distance to a nearby infected rhinoceros. For our initial exploratory analyses, our models included only crude evaluations of landscape-level effects and did not include the potential for infection in more than one nearby rhinoceros. We also did not have information on infection status of all (unsampled) individuals within a rhinoceros group, which may further lead to misclassification of exposure that could mask associations (43). More refined, individual-level measures that include longitudinal social network effects may improve understanding of the potential for intraspecies transmission.

Our findings also show a strong temporal association with risk of M. bovis infection, especially among black rhinoceros, albeit data were sparse for this group. Infection prevalence and odds of M. bovis infection were significantly higher in the earlier years of the study compared to 2019 and 2020. The magnitude of the adjusted OR was substantial in 2016 (OR = 18.1) and high in 2017 and 2018 (OR = 4.4 and OR = 4.6, respectively), compared to the most recent years of 2019 and 2020 (combined). It is unclear whether the observed decreasing prevalence over the study period was a result of changes in infection incidence or infection clearance rates, changes in contact patterns, an artifact of sampling bias, or another unobserved cause. All of these scenarios could be considered plausible.

Importantly, severe drought occurred in the year preceding and the first year of this study (2015–2016) (44). Associations between drought and M. bovis infection or disease have been detected in other multihost systems. In a large study of Mediterranean wild boars (n = 3,923), which are known bTB reservoirs, the occurrence of drought and increasing drought severity were significantly associated with increased occurrence of TB-like lesions (45). Another study investigating risk factors for bTB in cattle in Great Britain showed a positive association between atmospheric dryness and areas of high risk for M. bovis infection (46). Drought may impact rhinoceros body condition, with adverse consequences for immune responses (47, 48), and susceptibility to TB (49–51). Alternatively, reduced availability of water sources during periods of drought could lead to greater congregation of animals at the limited available drinking or wallowing sites. This could increase the risk of transmission through indirect interaction via M. bovis–contaminated water sources or surrounding resources, as is posited to occur in other multihost systems (19, 20). Multiyear, longitudinal investigations may elucidate potential interactions between disease, climate change, and other factors that vary over time, on rhinoceros’ susceptibility and exposure to M. bovis infection.

Of the 18.8% of rhinoceros that tested positive for M. bovis infection, most appeared clinically normal (83% of infected individuals were clinically normal, 68/82, and only 2/14 M. bovis–positive individuals with clinical abnormalities showed possible evidence of infection); this is consistent with reports that M. bovis infection does not always lead to clinical disease in rhinoceros. Findings from a study of experimentally infected white rhinoceros suggest that individuals that are healthy may be able to contain and clear the infection before developing active or overt disease (21, 52). However, it was difficult to confirm that there was no association between M. bovis infection and clinical disease or death in KNP rhinoceros; if present, this association could be obscured by the currently high mortality rate due to poaching (53).

Whether or not M. bovis infection is likely to progress to disease in KNP rhinoceros, the substantial infection burden in this population should raise concern that changes in (unknown) factors impacting disease progression could lead to increased TB-related morbidity and mortality in this population. This could have further negative consequences for survival of this population, which is already experiencing pressures associated with habitat loss, climate change, and poaching (44, 53, 54). Overall, more research that utilizes sensitive indicators of rhinoceros clinical health is required to improve understanding of M. bovis infection and pathogenesis.

In recent years, there have been numerous important diagnostic advances for the detection and characterization of M. bovis infection and disease in wildlife species, including free-ranging populations (12). These tools have been used in bTB surveillance efforts in wildlife species and, importantly, provide platforms to investigate bTB at the population level. Mycobacterial culture of tissue or secretions, followed by M. bovis confirmation using PCR, is a highly specific, gold standard diagnostic method, which can confirm true infection. However, this technique is most accurate when applied postmortem, and tissues from mortalities in wildlife populations are often unavailable. For this reason, most of our study population was defined according to IGRA (5) results, which were used to classify rhinoceros infection status antemortem. Importantly, there is a potential for misclassification of M. bovis infection status, due to the sole reliance on this single available diagnostic platform. The IGRA is a standard method for diagnosis of active and latent Mycobacterium tuberculosis infection in humans (55, 56), and for active M. tuberculosis complex (MTBC) infection in animals (25, 57), including rhinoceros (4, 5, 21, 52). A limitation of this test is that it may not always detect cases of recently acquired infection (52). A study using M. bovis experimentally infected rhinoceros showed that an immune response (detected by IGRA) was measurable within 1 mo to 2 mo after exposure by airway inoculation, decreasing gradually between 5 and 12 mo postinfection and reverting to negative results between 12 and 16 mo postinfection, indicating clearance of active infection (52). In general, the chronic nature of M. bovis infection makes it difficult to reliably estimate infection incidence, although the presence of a positive IGRA result signifies current infection and, therefore, can be used to estimate prevalence. To minimize misclassification of infection status, 38 individuals that could not be defined as positive or negative were removed (refer to ). Efforts focused on improving M. bovis diagnosis in rhinoceros are currently ongoing and may lead to more accurate and reliable case classification in future studies.

The spatial analyses in this study relied on capture location data for the study population as proxy for the location or area where each rhinoceros may have been exposed during the study period, since information on the individual ranges was not available. Our rationale for using our own location data for individuals captured multiple times over the study was to provide contemporaneous and contextual estimates of the spatial scale in KNP over which study individuals could move (and therefore be exposed to spatial factors). The chosen sizes for potential home ranges around each individual’s capture location were approximated from our data based on the distances that rhinoceros traveled, as indicated by repeated captures, with 95% of point-to-point distances traveled falling within 23 km. These were consistent with limited reports of rhinoceros movement within home ranges, typically reported to have sizes between 5 and 65 km2, in the literature (41, 58–68). However, the single-time-point capture locations do not reflect the movement patterns of rhinoceros and provide only a crude representation of the area that the rhinoceros regularly inhabited. Similarly, movements of other rhinoceros, and potential maintenance hosts (buffalo and kudu), were based on single-time-point location data for a cross-sectional sample of individuals, and the true rhinoceros home ranges are likely to vary over time according to the landscape and the resource availability, and by individual (40, 41, 59–68). Additionally, information on the infection status of these potential maintenance hosts was also not available; therefore, the applied model assumed a uniform risk over space and time, which is unlikely to be the case. Lack of precision in these measures would be expected to misclassify exposure variables, often biasing associations toward the null (43). The fact that we detected associations between numbers of nearby buffalo herds and risk of M. bovis infection in rhinoceros may suggest that the true magnitude of the association is higher than our estimate. Future studies could circumvent these challenges by tracking each rhinoceros over the course of the study, using satellite trackers to record their movement over time and space.

Conclusion

This study examines the epidemiology of M. bovis in a free-ranging population of rhinoceros and includes a large sample population from what is historically the world’s largest population of free-ranging rhinoceros in KNP. We detected evidence of widespread M. bovis infection in African rhinoceros in KNP, with a substantial infection burden (M. bovis prevalence was 17.0% [63/317; 95% CI: 11.0 to 23.9%] for white rhinoceros and 11.2% [19/120; 95% CI: 3.1 to 22.2%] for black rhinoceros), the extent of which was previously unknown. This emphasizes the importance of cross-species surveillance in bTB-afflicted multihost systems. Since bTB can affect wildlife, domestic animals, and humans, its spread to different areas could have serious consequences for human and animal health and, consequently, the agriculture and tourism industries in southern Africa. For rhinoceros specifically, translocation to other populations is an integral part of conservation strategies but may be accompanied by the risk of introducing novel pathogens, including M. bovis, into other ecosystems. Due to the presence and widespread impact of M. bovis in KNP rhinoceros, imposed quarantine and testing requirements prior to translocation are warranted across both rhinoceros species and all age groups.

Results from this study also highlight the potential role of different factors in infection risk for each species. In black rhinoceros, we found temporal associations, with a higher risk of infection in individuals sampled in the early years of the study compared to the final years (2019-2020). This may suggest the involvement of drought and changing climatic conditions in infection risk in black rhinoceros. In white rhinoceros, we found an association with distribution of buffalo in the individual’s surrounding vicinity. This highlights the potential role of buffalo, a recognized bTB maintenance host, in infection of white rhinoceros. Further study of M. bovis risk in these populations is warranted.

Future related work should focus on the development of diagnostic tools that may improve surveillance in these species. These techniques would also enhance the ability to classify cases and improve resolution to understand the epidemiology of bTB in complex systems. In particular, a targeted cohort study that tracks individual rhinoceros longitudinally and measures their resource usage in relation to their infection status could aid in further developing the understanding of infection and transmission risk. Characterizing threats to the survival of these species in the KNP ecosystem is vitally important for conservation and protecting other vulnerable populations.

Materials and Methods

Source Population and Data Collection.

Black and white rhinoceros populations in the KNP were sampled opportunistically during postmortem examinations (n = 9) or immobilizations performed as part of management and veterinary activities conducted in 2016–2020. In total, 528 rhinoceros (130 black rhinoceros and 398 white rhinoceros) were sampled and considered for inclusion in this study. Data collected for individual rhinoceros included date of sample collection, GPS coordinates for capture locations, demographic characteristics (sex, species, and age class), and general health status of the animals prior to immobilization.

During immobilization, whole blood was collected from the auricular or radial vein of the rhinoceros in lithium heparinized vacutainer tubes (BD Biosciences), as previously described (3). Postmortem tissue samples were collected from nine white rhinoceros during necropsy. These samples included submandibular, retropharyngeal, cervical, prescapular, axillary, inguinal, mediastinal, tracheobronchial, and mesenteric lymph nodes and lung, which were frozen at −20 °C for transport to Stellenbosch University for further laboratory testing under biosafety level 3 conditions, as previously described (3).

All living animals (n = 519) were immobilized by wildlife veterinarians for management, or by other approved procedures according to KNP’s Wildlife Veterinary Services’ standard operating procedures for the capture, transportation, and maintenance in holding facilities of wildlife (South African National Parks). Ethical approval for this project was granted by the Stellenbosch University Animal Care and Use Committee (ACU‐2020‐19019), and a section 20 research permit was issued by the Department of Agriculture, Land Reform and Rural Development (DALRRD; 12/11/1/7/2).

Study Design and Study Population.

A cross-sectional retrospective study design was used to identify factors associated with M. bovis infection in rhinoceros from KNP. An individual rhinoceros from the sampled source population (described above) was included in this study if 1) a blood sample was obtained from the individual during immobilization, or tissues were sampled from the individual at necropsy, and 2) it was sampled while free ranging in KNP or within 2 mo [this is the expected length of time after infection with M. bovis within which an immune response is detectable using QFT-IGRA (52)] after translocation out of the park to a quarantine location. A total of 475 out of the 528 rhinoceros (90%) from the source population met these inclusion criteria.

Determination of M. bovis infection status.

The M. bovis infection status of individual rhinoceros was determined using one of the previously described methods: 1) QFT-IGRA (4, 5) (sensitivity = 78%; 95% CI: 52.3 to 93.5%: specificity = 92%; 95% CI: 63.9 to 99.8%) or 2) mycobacterial culture for M. bovis isolation from a (postmortem) tissue (3, 69) with RD-PCR for M. bovis species confirmation (26). QFT-IGRA is a standard method for diagnosis of active or latent M. tuberculosis infection in humans (56, 70), and for MTBC infection in animals, including rhinoceros (4, 5, 21, 52). A rhinoceros was classified as M. bovis infected if it had a positive IGRA result (antigen-specific TB response ≥ 21 pg/mL) or positive BACTEC mycobacterial growth indicator tube (MGIT) culture result with subsequent PCR identification of M. bovis. An individual was assigned a negative infection status if it had a negative IGRA result (antigen-specific TB response ≤ 21 pg/mL, mitogen response ≥ 21 pg/mL, nil response ≤ 21 pg/mL) and, if conducted, any MGIT culture result with subsequent PCR that did not identify the presence of M. bovis. Individuals for whom the M. bovis status could not be defined according to these criteria were classified as “unknown” infection status.

A small number of rhinoceros (n = 38) had multiple immobilizations and QFT-IGRA results during the study period. For these individuals, a negative infection status was assigned if all test results (QFT-IGRA and mycobacterial culture, if applicable) were negative. A positive infection status was assigned if any of the tests (QFT-IGRA, and/or mycobacterial culture with RD-PCR confirmation of M. bovis) were positive. Rhinoceros were considered positive for M. bovis infection on the date of the first positive test result. Importantly, data assigned to each of the individuals in this subset were associated with a single GPS point location at which the animal was sampled—either the point of capture at which an M. bovis–infected individual first tested positive for M. bovis, or, for individuals that consistently tested negative for M. bovis, a randomly selected point from their multiple capture locations.

Evaluated Risk Factors.

Risk factors hypothesized to be associated with M. bovis infection in black and white rhinoceros populations from KNP (1) were evaluated in this study based on availability of data at each individual’s sampling event. For M. bovis–positive rhinoceros that were sampled multiple times, the immobilization date (and corresponding data) associated with the first positive sample was included. For M. bovis–negative rhinoceros that were sampled multiple times, a single date (and corresponding data) was randomly selected from among all of that individual’s capture dates. Evaluated factors are further described below and include the following 13 variables: species, sex, age class, orphan status, health status at time of sampling, sampling year, sampling season, nearest permanent water source type, distance to nearest permanent water source, distance to the nearest M. bovis–infected rhinoceros, number of nearby kudu herds, number of nearby buffalo herds, and buffalo density.

Demographic, health, and temporal risk factors.

Demographic factors selected and evaluated in this study included species (white or black rhinoceros), sex (male, female), and age. Age was estimated by field veterinary staff and categorized as follows: adult (>7 y), subadult (>2 y to 7 y), and calf (0 y to 2 y). Calves were further classified as orphaned or with their mother at the time of sample collection. Health status at the time of sampling was assessed by veterinary staff as normal or abnormal. Examples of conditions associated with abnormal health status included poor body condition, visible injuries, or any treatment undergone for illness or injury at the time of sampling. The health status variable was used to determine whether there was an association between M. bovis infection and health status. Temporal factors like year of sampling (2016–2020) together with rainfall season, that is, dry (March–August) or wet (September–February), were evaluated.

Spatial risk factors.

Spatial risk factors were based on a single GPS point for each rhinoceros’ immobilization location plotted onto a map using geographic information system software (GIS; ArcGIS Pro, version 2.8; Environmental Systems Research Institute). Spatial data were further processed and evaluated in the GIS software. Rhinoceros that did not have a GPS point recorded (n = 4) were omitted from the spatial analyses.

Spatial data describing the distance between the rhinoceros capture point and the nearest water source, or other M. bovis–infected rhinoceros were summarized. Risk of M. bovis infection was then evaluated as a function of distance (continuous predictor).

A circular polygon buffer was placed around each rhinoceros capture location to approximate a crude home range, in which exposure to African buffalo and greater kudu, which are known bTB maintenance hosts (14), could occur. The radius of the home range was derived from a subset of rhinoceros (n = 38) with GPS coordinates from multiple immobilization events occurring within a maximum of 3 y of each other. The distribution of distances between pairs of capture points for the same individual (n = 70 total pairs of capture points for the 38 rhinoceros individuals; 4 individuals had four immobilization events, 7 individuals had three immobilization events, and 27 individuals had two immobilization events) was examined, and 95% of pairwise distance observations occurred within 23 km of each other. This distance served as an approximation of the upper limit of the distance that a rhinoceros would travel; however, it is assumed that most of the movement probably occurs within a smaller core area (3, 41, 68) of unknown size. Therefore, the size of the circular home ranges was varied to include radii at 75%, 50%, and 25% of the maximum, corresponding to 17.25, 11.5, and 5.75 km, respectively. Variables characterizing relative exposure to African buffalo and greater kudu (further described below) were summarized for each of the four circular home range sizes.

African buffaloes in home range.

Two different spatial data layers were available to estimate the presence and density of African buffaloes in each rhinoceros’ assigned home range. The first dataset was a zero-inflated Poisson model generated by Hughes et al. (71) for prediction of buffalo density per square kilometer across KNP. The predicted buffalo density map was overlaid with rhinoceros home ranges to derive an estimated buffalo density (per square kilometer) value for each rhinoceros as a proxy for M. bovis exposure due to the presence of these maintenance hosts (9, 10, 14, 72).

The second dataset was obtained with permission from SANParks GIS Scientific Services, and contained buffalo census data (describing distribution of herds and individuals) that were collected across KNP in 2015 and 2017 using aerial line transect sampling and distance analysis methods, as previously described (73). The census data for the 2 y were combined into a single mapped data layer, which was applied as a crude estimate of buffalo distribution in KNP. The combined census data were then overlaid with the rhinoceros' home ranges to estimate the number of buffalo individuals and herds within the rhinoceros home range.

Greater kudu in home range.

Overlap of rhinoceros distribution with kudu was evaluated because greater kudu are considered M. bovis maintenance hosts (14). Census data on greater kudu were collected in three separate years (2014, 2016, and 2017), using aerial line transect sampling and distance analysis methods as previously described (73), and made available by SANParks GIS Scientific Services. The three datasets were combined into a single mapped data layer, providing a crude measure of kudu distribution in KNP. The combined kudu census data were then overlaid with each rhinoceros’ home range to estimate the number of kudu individuals and herds that each rhinoceros may have been exposed to.

Proximity to water source(s).

Mapped datasets of the rivers and water holes in KNP were provided by SANParks GIS Scientific Services. The river dataset described main and secondary rivers and was compiled in 2018 using older data sources in combination with National Geo-spatial Information aerial imagery (74). The water hole dataset described the location of available water holes in KNP, updated through June 2016.

The distance (kilometers) between each rhinoceros’ capture location and the nearest water source (river or water hole) was determined in the GIS; this value was used to represent the proximity of each rhinoceros to water. This variable was tested based on the hypothesis that aggregation of infected hosts at available water sources may result in increased infection exposure of susceptible hosts living in close proximity to these water sources, either through direct interactions with infected hosts or due to mycobacterial loads shed into the environment. The nearest water source type (river or water hole) was also evaluated as an independent risk factor.

M. bovis status of nearby rhinoceros.

A continuous variable was created to evaluate whether the risk of M. bovis infection was a function of the infection status of other nearby rhinoceros. All rhinoceros capture location points were plotted in ArcGIS, and the infection status of each animal was determined. The distance (kilometers) from each study subject to the nearest M. bovis–positive rhinoceros was then determined with ArcGIS and ascribed to the study subject.

Data Analyses.

Apparent M. bovis infection prevalence was estimated for the full KNP study population (number of test-positive rhinoceros/total number of study rhinoceros), as well as within different ranger sections (75) and different ecozones (number of test-positive rhinoceros in specific area/total number sampled in specific area) (76). Within ranger sections and ecozones, prevalence calculations were limited to the areas where >10 animals were sampled. These apparent prevalence values were then adjusted to account for the sensitivity and specificity of the IGRA (sensitivity = 78%; 95% CI: 52.3 to 93.5%: specificity = 92%; 95% CI: 63.9 to 99.8%), thereby estimating true prevalence, using the following equation: Estimated true prevalence = (apparent prevalence + IGRA specificity – 1)/(IGRA sensitivity + IGRA specificity – 1) (based on formulas implemented in ref. 27). Prevalence values were compared across ranger areas and ecozones using Fisher’s exact tests.

Differences in geographical distribution of M. bovis infection were further explored using Kulldorff’s spatial scan statistic (28). The statistic was applied using a Bernoulli based model and SatScan software (version 10.0). SatScan is a trademark of Martin Kulldorff and developed under the joint auspices of Martin Kulldorff, the National Cancer Institute, and Farzad Mostashari of the New York City Department of Health and Mental Hygiene (28). The methods identify significant case clustering by moving a circular window over the geographic area; the maximum spatial cluster size was set to half of the population. For this statistic, the null hypothesis assumed that the relative risk of M. bovis infection is the same inside the geographic area compared to outside. Significance was determined by comparing likelihood ratio tests from 999 iterations of a Monte Carlo simulation. We performed this evaluation among all rhinoceros, and then among black and white rhinoceros separately.

Univariate logistic regression was used to screen for associations between each factor and M. bovis infection (Tables 2 and 3). Crude ORs, 95% CIs, and type III Wald’s P values were estimated. Evaluation of the association between M. bovis infection status and orphan status was completed only within the subset of calves. Functional forms of continuous variables were determined by fitting a logit-transformed Loess curve for single-variable models. Natural log transformations were used for covariates that were not normally distributed. If there was evidence of nonlinearity in the logit, then the variable was categorized into quartiles. Important demographic covariates, potential effect modifiers, and associations with P ≤ 0.2 were further examined in multivariable analyses.

The three factors derived from rhinoceros areas of exposure (i.e., number of kudu herds, number of buffalo herds, buffalo density) were further evaluated to determine the optimal spatial scale for each variable separately. For each variable, a single radius distance for area of rhinoceros exposure was selected for further evaluation in the multivariable model. The chosen distance was based on the strongest association and the best-fitting single-variable logistic regression using the Akaike information criterion (AIC) statistic (77).

Multivariable logistic regression analyses evaluated associations between multiple factors and M. bovis infection. A backward stepwise approach was used to fit models that included species and demographic factors (age and sex) as well as those that met inclusion criteria (P ≤ 0.2). Effect modification was evaluated by including an interaction term between plausible effect modifiers (species, age, and sex) and the other factors in the model. We also examined whether water type was an effect modifier of the association between distance to the nearest water source and M. bovis infection. Competing models with similar predictors were chosen based on the AIC. The final model included covariates (species, age, and sex) as well as other significant risk factors. Since significant interaction between species and year of sampling with M. bovis infection was identified in the final multivariable model with all rhinoceros, the same model was further explored for black and white rhinoceros separately, to describe potential differences in risk factors between species.

Spatial data processing and analyses were performed in ArcGIS (version 2.8), except for tests of spatial clustering performed with SatScan, as described above. Statistical analyses were performed in R (version 4.0; R Core Team); true prevalence was calculated with the package EpiR with the function epi.prev (27), and univariate and multivariable models were fit with the glm function (78). Associations with P < 0.05 were considered statistically significant.