Recent changes in populations of Critically Endangered Gyps vultures in India

Summary

Populations of the White-rumped Vulture Gyps bengalensis, Indian Vulture G. indicus and Slender-billed Vulture G. tenuirostris declined rapidly during the mid-1990s all over their ranges in the Indian subcontinent because of poisoning due to veterinary use of the non-steroidal anti-inflammatory drug diclofenac. This paper reports results from the latest in a series of road transect surveys conducted across northern, central, western and north-eastern India since the early 1990s. Results from the seven comparable surveys now available were analysed to estimate recent population trends. Populations of all three species of vulture remained at a low level. The previously rapid decline of White-rumped Vulture has slowed and may have reversed since the ban on veterinary use of diclofenac in India in 2006. A few thousand of this species, possibly up to the low tens of thousands, remained in India in 2015. The population of Indian Vulture continued to decline, though probably at a much slower rate than in the 1990s. This remains the most numerous of the three species in India with about 12,000 individuals in 2015 and a confidence interval ranging from a few thousands to a few tens of thousands. The trend in the rarest species, Slender-billed Vulture, which probably numbers not much more than 1,000 individuals in India, cannot be determined reliably.

Introduction

A ban on the veterinary use of the non-steroidal anti-inflammatory drug (NSAID) diclofenac in India was announced in 2006 and the official completion of the banning process was an extraordinary gazette notification in 2008 (Gazette of India Notification No. GSR 499(E)). The ban was an attempt to halt the precipitous decline in three species of ‘Critically Endangered’ vultures endemic to South Asia: White-rumped Vulture Gyps bengalensis, Indian Vulture G. indicus and Slender-billed Vulture G. tenuirostris. Veterinary use of diclofenac was the main and probably the only cause of these population declines. Evidence concerning the importance of diclofenac relative to that of other postulated causes of the decline has been presented in detail elsewhere (Oaks et al. 2004, Green et al. 2004, Shultz et al. 2004). Vultures die from diclofenac-induced kidney failure if they consume sufficient tissue from the carcass of an ungulate that has died within a few days of treatment with the drug. In the early 2000s, before the ban, the proportion of carcasses of domesticated ungulates in India contaminated with diclofenac and the concentration of the drug in their tissues were sufficient to have caused vulture declines at the observed rates without the involvement of any other factor (Green et al. 2007). Prakash et al. (2012) reported results from counts of the three Gyps vulture species on road transects in northern India in six comparable surveys between 1992 and 2011. They found that the rapid population declines of all three species up to 2002 had slowed and, in the case of White-rumped Vulture, possibly even reversed by 2011. In this paper, we report results from the latest of this series of counts: the seventh survey conducted in 2015.

Methods

Survey method and data limitations

Vultures were counted in 2015 on road transects distributed across northern (Haryana, Punjab, Uttar Pradesh, Uttarakhand, Bihar), central (Madhya Pradesh), western, (Rajasthan, Gujarat, Maharashtra) and north-eastern (West Bengal, Assam, Meghalaya, Arunachal Pradesh) India. Transect locations and methods followed those of similar surveys conducted in 1991–1993, 2000, 2002, 2003, 2007 and 2011 (Prakash et al. 2012). Transects covered in 2015 were located in and near protected areas (99 transects and total length 5,221 km) and along roads between protected areas (55 transects and total length 10,296 km). All transects had been covered in one or more of the previous surveys. The initial surveys in this series were conducted in one year of a three-year period (1991–1993). For the purpose of analysis we treated them as having been conducted in 1992, the mid-point of the period. In all years, surveys analysed here were conducted between March and July. This period was chosen because it includes the end of, and a period after, the vulture breeding season, so adult birds were not unavailable for survey because they were incubating eggs or brooding small nestlings. In addition, this period is in the early part of the monsoon season when road travel and observation were unlikely to be hampered by heavy rain. In 2015, four teams, each consisting of one observer and one driver, surveyed the four regions described above. Transects were driven starting between 07h00 and 11h00 and finishing between 15h00 and 19h00 local time at 10–20 km/h in and near protected areas and ∼50 km/h between protected areas. Vultures observed on the ground, in trees, on cliffs, flying and soaring within 500 m on either side of the transect were identified and recorded. All were fully-grown birds and our counts did not include nestlings at breeding sites. Because they are large birds (c.5 kg body weight), vultures were easy to detect within this distance without optical equipment, but identification of species was done using binoculars. Estimated distances to individuals were not recorded, so detection probability in relation to perpendicular distance to the transect line cannot be estimated to adjust the counts for detection probability. Observations were made from a vehicle and there was therefore little or no disturbance to the vultures. Vultures were identified as White-rumped Vulture, Indian Vulture and Slender-billed Vulture from 2002 onwards. Before 2002, Indian Vulture and Slender-billed Vulture were not identified or recorded separately because they were only recognised as two separate species in 2001 (Rasmussen and Parry 2001). For that reason, we analysed combined counts of Indian Vulture and Slender-billed Vulture together for two periods (1992–2015 and 2000–2015) and analysed counts for the two species separately in 2002–2015. Vultures were so numerous in 1991–1993 that only groups of five or more were recorded. Hence, differences in counts between 1992 and all subsequent years are underestimated, to an unknown extent. However, we believe that most vultures seen in 1992 were in groups of five or more, so this negative bias is probably small. No specific permits were required for the surveys, but permission was obtained to enter all protected areas. Further details of previous surveys are given in Prakash et al. (2007, 2012). A map showing the locations of transects is in Green et al. (2007).

Calculation of annual population indices

Not all transects were covered in every survey. Some were added to the set after 1991–1993, whilst others were temporarily or permanently omitted from the survey. To allow for the turnover and missing values, we fitted regression models that allowed for the effects of the changing composition of the sample of transects. We only included data from transects that were surveyed more than once in the study period and on which vultures of the focal species or species group had been recorded at least once. We called these informative transects. In these regression models, ‘count’ was the dependent variable and ‘transect’ and ‘survey year’ were fixed effect factors. Models were fitted in R, with a Poisson error term and a logarithmic link function. The form of the model was

$${C_{ij}} = \exp \left( {{k_i} + {p_j}} \right),$$

where C _ij is the count for the j ^th transect in the i ^th year. Site effects are represented by the regression coefficients p _j. The coefficients k _i represent the year effects and are the logarithms of the abundance of vultures in i ^th year, allowing for site effects, expressed as a proportion of the abundance of vultures in the first year of the series in the study period. Hence, exp(k _i) provides an index of population density in the i ^th year, relative to that in the first year. We obtained 95% confidence intervals for the population index values using a bootstrap method. In a period in which there were m informative transects eligible for analysis for a species or species group, we took random bootstrap samples of m transects, with replacement, from the m transects available. We then fitted the log-linear Poisson regression model for this bootstrap sample and recorded the value of exp(k _i) for each of the survey years. This procedure was repeated 1,000 times, the bootstrap estimates ranked, and the bounds of the central set of 950 estimates taken to define the 95% confidence interval of each of the population indices.

Calculation of mean annual population multiplication rate and changes in population trend over time

We estimated the mean annual rate of population change by fitting a Poisson regression model with a logarithmic link function and transect as a fixed factor, as before, but with the effect of year modelled as a continuous explanatory variable t, the number of years elapsed since the first survey of the series being used. Hence, t = 0 for 2000 and t = 15 for 2015. The form of the model was

$${C_{ij}} = \exp \left( {{b_0} + {b_1}\,t} \right),{\rm{\ }}$$

where C _ij is the count for the j ^th transect in the i ^th year, which is t years after the initial year of the series. We did this only for the period 2000–2015 because we considered it unwise to estimate the average annual rate of population decline over the earlier period 1992–2015, given that the rapid vulture population decline began at an uncertain time, but probably in the 1990s. The regression coefficient from this model b ₁ provides the mean annual population multiplication rate λ = exp(b ₁). To examine changes in population trend over time, we tested whether λ had altered significantly over time using a bootstrap method. We fitted Poisson regression models with a logarithmic link function and transect as a factor, similar to the previous model, but with the effect of the quadratic and cubic terms t ² and t ³ added. The forms of these models were

$${C_{ij}} = \exp \left( {{b_0} + {b_1}\,t + {b_2}\,{t^2}} \right),$$

and

$${C_i}_j = \exp \left( {{b_0} + {b_1}\,t + {b_2}\,{t^2} + {b_3}\,{t^3}} \right).$$

If the inclusion of the higher order polynomial terms significantly improved the fit of the regression, this would indicate that the population multiplication rate changed significantly over time. We tested this possibility using a backwards elimination bootstrap procedure. We drew 1,000 bootstrap samples of data, as described above, and fitted the cubic Poisson regression model to each bootstrap sample. We took the central 950 values of b ₃ as defining its 95% confidence limits and counted the number of bootstrap samples in which the coefficient was of opposite sign to that calculated from the full dataset. If the 95% confidence limits for b ₃ overlapped zero, we eliminated the cubic term and conducted the equivalent procedure for the quadratic term b ₂. If the 95% confidence limits for b ₂ overlapped zero, we eliminated the quadratic term and conducted the equivalent procedure for b ₁. We stopped this backwards elimination procedure if the confidence interval of the highest order regression coefficient remaining in the model did not overlap zero and accepted that model as the minimal adequate model.

Estimates of total population size

We used a regression model fitted to road transect survey data to estimate vulture density in relation to survey year and covariates and then used this model to estimate the total vulture population size for India. We analysed counts from 159 road transects in the years 2003, 2007, 2011 and 2015 for which we had information on the length of transect driven. We modelled the density of vultures recorded on the road transect surveys in relation to survey year, the geographical position of the centroid of the transect in India and the distance of the centroid of the transect from centroid of the nearest protected area. Protected areas varied considerably in extent (National Parks were up to 3,350 km² in extent and Wildlife Sanctuaries up to 8,500 km²), but we used the distance of the transect centroid to the centroid of the nearest protected area in our analysis for simplicity. The largest National Parks have an average diameter of about 60 km. We fitted Poisson regression models with a logarithmic link function. The dependent variable in the regression was the count of vultures of a particular species on each transect in one of the four survey years. We included the natural logarithm of the length of each transect in kilometres in the regression as an offset. This makes the model equivalent to one in which the dependent variable is the density of vultures per square kilometre, because the strip of land covered by each transect was one kilometre wide. Hence, each kilometre of transect driven represents a survey of one square kilometre. The effect of survey year was modelled by including it in the regression as a factor with four levels. It was necessary to take into account the geographical position of the transects because the geographic distributions of two vulture species (Indian Vulture and Slender-billed Vulture) do not extend to all parts of mainland India and the abundance of all three species is thought to vary geographically. We modelled the effect of transect position by including the latitude and longitude of the transect centroid in decimal degrees as continuous variables together with the squares of each of latitude and longitude. Hence, both latitude and longitude were modelled as having a quadratic effect on abundance. This allows the density of vultures potentially to have a hump-shaped relationship to latitude and longitude. In order to avoid the fitting of large numbers of regression parameters, we assumed that the coefficients of the functions relating density to latitude and longitude varied among vulture species, but were the same for a given species in all survey years. Confidence limits of regression coefficients were obtained by bootstrapping, with transects being used as the bootstrap units. To obtain each bootstrap sample, we drew 159 sets of count data by selecting results for transects at random from the original data, with replacement. We fitted the regression model to each of 1,000 bootstrap samples obtained in this way and took the central 950 of the bootstrap regression coefficient estimates as 95% confidence limits.

Results from analysis of a previous survey in 2011 indicated that most vultures were located in or near National Parks (Prakash et al. 2012), so we modelled vulture densities in relation to the proximity of the transect to protected areas. We used the 2014 United Nations List of Protected Areas of India (Deguignet et al. 2014), supplemented by internet searches, to obtain the centroids in decimal degrees of latitude and longitude of all 79 National Parks and all 338 wildlife sanctuaries larger in extent than 10 km² in mainland India. We calculated the geodesic distance in kilometres between the centroid of each transect and the centroids of all protected areas; and then found the distance from each transect centroid to that of the nearest National Park (NPD) and the distance to the centroid of the nearest Wildlife Sanctuary (WSD). National Parks and Wildlife Sanctuaries are both types of protected areas. Although actual levels of protection of National Parks and Wildlife Sanctuaries vary considerably, National Parks tend to have greater emphasis on restrictions on human activities and maintenance of natural ecosystem function than Wildlife Sanctuaries. National Parks are accorded a higher status than Wildlife Sanctuaries (Category II vs Category IV) in IUCN’s global classification of types of protected areas (Deguignet et al. 2014). We included NPD and WSD as continuous variables in the regression models.

We considered possible methods to allow for the effects of spatial autocorrelation in the transect data used to fit the regression models. Statistical methods are widely used for this purpose for models with a normally distributed continuous response variable or counts and data from evenly distributed grids of sampling points (Dormann et al. 2007), but appropriate methods for the data such as ours with irregularly distributed sampling sites, a Poisson dependent variable with many zeros and offsets are less easily implemented and less thoroughly tested. We therefore fitted the Poisson models as described, without allowing for spatial autocorrelation, and then performed a global Moran’s I test on the residuals from the selected final model (see Results). The residuals were the differences between the observed mean density across the survey years in the period 2003 – 2015 and the expected mean density from the regression model for that period. We used the reciprocal of the geodesic distance between transect centroids as weights in the calculation of Moran’s I.

We used the regression model fitted to combined data for all three species, with the effect of NPD included, to estimate the total numbers of vultures of each species in mainland India. To do this, we obtained the latitude and longitude of the centroids of all 3,278,983 1-km squares in mainland India and the geodesic distance, in kilometres, of each 1-km square centroid to the centroid of the nearest National Park. Using the parameter estimates from the regression model, we calculated the expected number of vultures in each square from its latitude, longitude and distance from its centroid to that of the nearest National Park and summed the expected numbers across all 1-km squares to give a total for India for each species in each of the four survey years. To obtain confidence limits for these estimates, we used the 1,000 sets of bootstrap estimates of the parameters of the regression model and used the method described above to calculate estimated vulture populations from each set. We took the central 950 of the bootstrap population estimates for each species and survey year as the 95% confidence limits of the population estimates.

Results

Annual population indices

The total numbers of White-rumped Vultures, Indian Vultures and Slender-billed Vultures counted in 2015 were 102, 139 and 12 respectively, compared with 99, 299 and 15 in 2011. The annual indices of population density differed little between 2011 and 2015 for White-rumped Vulture (Table 1, Figure 1), but the 2015 index for Indian Vulture and Slender-billed Vulture combined was about half of that in 2011, after being approximately stable since 2003 (Table 1, Figure 2). Populations of both these species groups in 2015 remained low relative to the 1992 level: about one five-hundredth of the 1992 level for White-rumped Vulture and about one-hundredth of the 1992 level for Indian Vulture and Slender-billed Vulture combined. Too few Slender-billed Vultures have been counted per survey to quantify a reliable trend for this rare species separately, but the index values obtained since they were first counted separately in 2002 suggest an initial decline between 2002 and 2003 and no consistent trend since then (Table 1).

Table 1. Indices of population size for White-rumped Gyps bengalensis, Indian G. indicus and Slender-billed G. tenuirostris Vultures in India across three periods. Indices are population densities, relative to those of the first year of the respective period indicated, estimated by log-linear Poisson regression performed on data from road transect surveys in northern India. Also shown are 95% bootstrap confidence intervals for each index (in brackets) and the number of informative transects used in each analysis.

Figure 1. Population indices and trend of White-rumped Vulture in India. Points show indices of population density, relative to that in 1992, estimated by log-linear Poisson regression performed on data from seven road transect surveys in northern India. Vertical lines show 95% bootstrap confidence intervals. The curve is the quadratic log-linear population trend fitted to data for the period 2000–2015.

Figure 2. Population indices and trend of Long-billed Vulture and Slender-billed Vulture combined in India. Points show indices of population density, relative to that in 1992, estimated by log-linear Poisson regression performed on data from road transect surveys in northern India. Vertical lines show 95% bootstrap confidence intervals. The line is the log-linear population trend fitted to data for the period 2000–2015.

Changes over time in annual population multiplication rate

Bootstrap tests on cubic and quadratic regression models of population density on year were used to determine whether the annual population multiplication rate has changed significantly since 2000 (see Methods). For White-rumped Vulture, in the model with both quadratic and cubic terms, the 95% confidence interval of the regression coefficient for the cube of years elapsed since 2000 overlapped zero by a wide margin (coefficient = +0.000952, 95% CL -0.004420 to +0.006852) and the sign of the coefficient was opposite to that fitted to the full dataset for a large proportion (0.336) of bootstrap samples. We concluded that the data do not justify the inclusion of the cubic term in this model and it was deleted. However, in the model with the quadratic term, the 95% confidence interval of the regression coefficient for the square of years elapsed did not overlap zero (coefficient = +0.02904, 95% CL +0.01271 to +0.04480) and the sign of the coefficient was opposite to that fitted to the full dataset for a very small proportion (0.001) of bootstrap samples. We therefore concluded that the inclusion of the quadratic term was justified. The fitted regression model for the relationship between population index relative to that in 2000 and years since 2000 was index = exp(-0.6524 years + 0.02904 years²). The significantly positive quadratic regression coefficient indicates a departure from continuous exponential population decline at a constant proportion per year for White-rumped Vulture. Instead, the rate of decline has slowed significantly since 2000 and the population has stabilised and may be increasing (Figure 1).

For Indian Vulture and Slender-billed Vulture combined, in the model with both quadratic and cubic terms, the 95% confidence interval of the regression coefficient for the cube of years elapsed since 2000, overlapped zero substantially (coefficient = -0.00162, 95% CL -0.00583 to +0.00324) and the sign of the coefficient was opposite to that fitted to the full dataset for a large proportion (0.242) of bootstrap samples. We concluded that the available data do not justify the inclusion of the cubic term in this model and it was deleted. The equivalent analysis for the quadratic term in the quadratic model also indicated that its inclusion in the regression model was not justified by the data (quadratic coefficient = +0.00573, 95% C.L. -0.0081 to +0.01973) and it was deleted. The sign of the coefficient was opposite to that fitted to the full dataset for a large proportion (0.202) of bootstrap samples. However, the bootstrap test on the regression coefficient for the first-degree term b ₁, in the model containing just this term, indicated strong evidence for a negative trend in population index since 2000 (first-order coefficient = -0.1182, 95% CL -0.1821 to -0.0652). The fitted regression model for the relationship between population index relative to that in 2000 and years since 2000 was index = exp(-0.1182 years). Hence, the analysis of survey results for Indian Vulture and Slender-billed Vulture combined indicates a continuous exponential population decline since 2000 at a constant proportion per year with no indication that the decline has slowed (Figure 2).

Estimates of total population size

Regression analysis indicated that there was a consistent negative effect of increasing distance to the centroid of the nearest National Park on the density of vultures (Table 2). The fitted regression coefficients for this variable (NPD) were negative and of similar magnitude for all three species, whether NPD was fitted on its own or in models that also included distance to centroid of the nearest Wildlife Sanctuary (WSD). The bootstrap 95% confidence limits for the regression coefficient for NPD did not overlap zero for any of the models. The coefficient was also negative for the model fitted to data for all three species. In contrast, the regression coefficients for WSD were not consistent in sign. In the models with effects of both NPD and WSD, the bootstrap 95% confidence limits for the regression coefficient for WSD overlapped zero, except for the model for White-rumped Vulture, where the coefficient was positive and almost overlapped zero. Because any possible effect of WSD was weak and inconsistent, we used only the regression model with the effect of NPD alone for further analyses. The coefficient for the effect of NPD on vulture density appeared to be similar for all species, so we used the model with a coefficient common to all species, fitted to the combined data for all three species for population estimation. There was little evidence of spatial autocorrelation of the residuals of mean vulture density from this model (Moran’s I = -0.00497, standard deviate -0.193). The relationship between vulture population density and NPD is illustrated in Figure 3.

Table 2. Poisson regression models of the effects of distance from the transect centroid to the centroid of the nearest National Park (NPD) and distance to the centroid of the nearest Wildlife Sanctuary (WSD), in kilometres, on the density of vultures per square kilometre observed on road transects in the years 2003, 2007, 2011 and 2015. All analyses include data for all four survey years. All models include effects on density of survey year (as a factor), latitude and longitude (both as quadratic models: coefficients not shown). The first three models in each section of the table were fitted separately for each species. The fourth model is for all species combined with the main effect of species and two-way interactions of species with survey year, latitude and longitude included. The upper part of the table shows results for models with both NPD and WSD and the lower part shows results with only NPD included. 95% confidence limits were obtained by bootstrapping, but could not be calculated for the models for Slender-billed Vulture.

Figure 3. Numbers of vultures recorded per square kilometre of road transect surveyed in relation to the distance of the centroid of the transect to the centroid of the nearest National Park. Transects were grouped into three categories: less than 50 km from a National Park, 50–100 km and over 100 km. The total number on vultures seen in all four survey years was divided by the number of square kilometres surveyed and the resulting densities are plotted against the mean distance from a National Park of the transects in each distance category. Results are shown separately for White-rumped Vulture (diamonds), Indian Vulture (squares) and Slender-billed Vulture (triangles).

Total vulture populations estimated for 2003, 2007, 2011 and 2015 from the regression model are shown in Table 3. Estimated numbers changed between years in the way expected from the population indices shown in Table 1. The confidence limits of estimated population sizes were wide, with the upper 95% limit being about three times the estimate and the lower limit being about one-third of the estimate, even for the most abundant of the three species (Indian Vulture). Hence, vulture population sizes are estimated only crudely by this method to about one order of magnitude. We were unable to calculate confidence limits for the rarest species Slender-billed Vulture because the small number of transects upon which it was recorded prevented the reliable implementation of the bootstrap procedure.

Table 3. Estimates of population size in each survey year for three species of Gyps vultures in mainland India calculated from a regression model of density in relation to survey year, latitude, longitude and the distance to the centroid of the nearest National Park. 95% confidence limits (CL) were obtained by bootstrapping, but could not be calculated for Slender-billed Vulture.

Discussion

For White-rumped Vulture, our latest update in 2015 of a previous series of road transect surveys (Prakash et al. 2012) indicates that the rapid decline in numbers of this species, which began in the mid-1990s, stopped in about 2010. The population has stabilised since then or may be increasing slowly. However, the total population of this species in India is precariously small. Our estimate based upon a regression model is that there were about 6,000 individuals in 2015, with a confidence interval ranging from less than one thousand to a few tens of thousands.

Indian Vulture and Slender-billed Vulture were not considered to be different species until 2001 and we do not have separate information on population trends of these species until after 2003. Previous indications that the population index values for Indian Vulture and Slender-billed Vultures combined had stabilised between 2003 and 2011 (Prakash et al. 2012) are not confirmed by our latest results. Addition of the new survey results for 2015 suggests instead that populations of Indian Vulture and Slender-billed Vulture have been continuing to decline, albeit at a much slower rate than was the case for White-rumped Vulture up to about 2010. Counts of White-rumped Vulture nests at Keoladeo National Park suggest that the decline of that species began in 1994, which was also the median year for first veterinary use of diclofenac reported by Indian veterinary professionals (Cuthbert et al. 2014). Assuming that the rapid declines of Indian Vulture and Slender-billed Vulture also began in 1994, the mean annual rate of decline rate of these species between 1994 and 2000 was about 35% per year (100 (1-0.0751^1/6)), compared with a mean rate of decline for these two species combined from 2000 to 2015 of 11% per year. We estimated the total population of Indian Vulture in India at about 12,000 individuals in 2015, with a confidence interval ranging from a few thousands to a few tens of thousands. Our survey data were too sparse to estimate a confidence interval for the population of Slender-billed Vulture, but our best estimate is that there were a little over one thousand individuals in India in 2015.

A ban on veterinary use of diclofenac in India was first announced in 2006. Repeated surveys of the prevalence and concentration of diclofenac in tissues from carcasses of domesticated ungulates available to vultures in India showed that they both declined after the ban. The expected risk of death from diclofenac poisoning per meal for White-rumped Vulture, calculated from these ungulate survey data, had fallen to one-third of its 2006 level by 2009 (Cuthbert et al. 2014), but post mortems and tissue analyses showed that wild Gyps vultures in India continued to die from diclofenac poisoning, though probably at a lower rate than before the ban (Cuthbert et al. 2016). Simulation models of the Indian population of the White-rumped Vulture Population models indicate that the observed cessation of the decline for this species is in accord with the change in vulture population trend expected from data on diclofenac contamination of ungulate carcasses (Prakash et al. 2012). However, these findings do not throw any light on why the decline in the combined populations of Indian Vulture and Slender-billed Vulture did not cease or slow significantly after the diclofenac ban. Much of the continuing exposure of vultures to diclofenac is attributable to the illegal sale for veterinary use of diclofenac formulated for use in human medicine in large multiple-dose vials (Cuthbert et al. 2011). The large vials are convenient for injecting large-bodied domesticated ungulates. In 2015, the Government of India banned the manufacture of human formulations of diclofenac in multiple-dose vials (Gazette of India Notification GSR 503(E)), and it is hoped that this will further reduce exposure of vultures to diclofenac.

In addition to the continuing threat from diclofenac, other veterinary NSAIDs that are toxic to Gyps vultures are approved for legal use in India and are likely to be causing mortality. These include ketoprofen, for which there is experimental evidence of toxicity to vultures below the maximum level of exposure for White-rumped Vulture (Naidoo et al. 2010) and aceclofenac, which is largely metabolised to diclofenac within cattle (Galligan et al. 2016). In addition, nimesulide residues have been found associated with visceral gout in vultures found dead in the wild in India (Cuthbert et al. 2016). Although experimental tests of nimesulide on captive vultures have not yet been done, the co-occurrence of nimesulide residues and visceral gout in dead vultures makes it probable that nimesulide is nephrotoxic to vultures. At present, meloxicam is the only NSAID known not to be toxic to vultures and other scavengers at levels up to the maximum likely level of exposure (Swan et al. 2006, Swarup et al. 2007).

Other actual and potential threats to vulture populations in India and changes in their prevalence are poorly quantified. Poisoning is a frequent cause of death of vultures throughout the Old World, including Europe, South East Asia and Africa, where poison baits that are usually set to kill other species kill vultures incidentally (Hernández and Margalida 2008, Clements et al. 2013, Ogada et al. 2016). The baits, which often use widely available agricultural pesticides, also kill vultures that scavenge the carcass. Poison baits are set in India at carcasses of domesticated ungulates killed by mammalian carnivores, such as feral dogs and jackals, to kill them. It seems likely that a vicious circle has occurred in which populations of feral dogs have increased because of the increased cattle carrion food supply no longer consumed by vultures (Markandya et al. 2008). This may have led to more killing of livestock by dogs and other scavenging mammals and more reprisal poisoning. However, this is conjecture. The numbers of vultures of these three species reported dead from this cause annually in India is small, but the degree to which instances of it are detected, reported and correctly attributed is uncertain and difficult to estimate. Similar lack of robust quantification applies to other causes of death. It is hoped that future recovery for post-mortem studies of carcasses of wild vultures fitted with GPS tags will allow the estimation and comparison of per capita annual death rates from NSAID poisoning, poison baits and other causes. However, such studies have yet to be conducted. Nonetheless, estimates of per capita additional mortality rates of vultures due to diclofenac poisoning have been made based upon two types of data: (1) proportions of dead vultures with diclofenac residues and visceral gout (Green et al. 2004; Cuthbert et al. 2016), and (2) surveys of diclofenac prevalence and concentration in carcasses of cattle available to scavengers (Green et al. 2007, Prakash et al. 2012). Both of these sets of results indicate a high level of additional mortality of vultures due to diclofenac which was sufficient to account for the observed rate of population decline without the involvement of other causes. In addition, recent changes in diclofenac prevalence after the ban on its veterinary use were sufficient to account for changes in the observed rate of population decline of vultures (Prakash et al. 2012).

Our estimates of total populations of vultures in India in 2007 were smaller than those made by Prakash et al. (2007) for the same year. This difference occurred despite the fact that Prakash et al. (2007) only calculated total population size for part of India (about 80% of the land area, excluding Goa, Andhra Pradesh, Teleganga, Karnataka, Kerala and Tamil Nadu), whereas we did so for the whole of the Indian mainland. The explanation for this difference in estimates for 2007 is that the method of Prakash et al. (2007) assumed that transects were randomly placed and did not take into account distance from National Parks. Road transects were not located randomly with respect to the distance from National Parks. More transects were positioned with their centroids close to the centroid of the nearest National Park than would be expected by chance (Figure 4), because the survey was designed to repeat, in part, surveys of all raptors conducted in the early 1990s, before the vulture population decline began. In these initial surveys, many transects were deliberately placed in and near protected areas so as to increase survey coverage of scarce raptor species, some of which are reliant on natural ecosystems protected in National Parks.

Figure 4. Proportions of 1-km squares in mainland India (light grey bars) and road transect surveys (dark grey bars) in 50-km categories of distance between the centroid of the square or transect and the centroid of the nearest National Park.

Our estimates of total population size are subject to several caveats because of limitations in the data available. The first caveat is that, we estimated populations for the whole of mainland India, but did not conduct surveys in every state. We have survey data from 13 states of mainland India, which comprise 58% of its land area. Sampled states were widely distributed in the northern two thirds of India by latitude, which comprise about 80% of the India’s land area. We suggest that extrapolation of our regression model of population density to the unsampled northern states may be quite accurate, given that we allowed for geographical effects by including quadratic effects of latitude and longitude in our regression models. However, no surveys were done in any of the southern states of Goa, Karnataka, Andhra Pradesh, Teleganga, Kerala and Tamil Nadu and parts of that region are about 1,000 km from the nearest transect. Hence, extrapolation to that region is less secure. However, we believe that errors introduced by this extrapolation to the south are unlikely to be large because there are probably relatively few vultures there and this is reflected in our models. Our opportunistic observations suggest that average densities of two of the three vulture species are much lower in the south than in the north. For the third species, the Slender-billed Vulture, the breeding range does not include the south of India (del Hoyo and Collar 2014). This conclusion is reflected in results from our regression models of the effects of latitude and longitude within the sampled region. These models predict densities of all three species at a typical latitude of the unsampled southern region (13°N) less than 1% of the density at a typical latitude of the sampled northern region (25°N) because of marked north-south negative trends in density within the sampled area. In addition to these low expected densities in the south, the southern region comprises only 19% of the area of mainland India. Hence, we believe that total numbers of vultures in the unsampled southern region are likely to be small and that errors in the predicted numbers due to extrapolation are unlikely to cause substantial bias in the total population estimates.

The second caveat about our population estimates is that they are based on data from road transects. Roads are not placed in representative parts of the landscape and therefore average vulture densities along roads might not be representative of those in India as a whole. In the absence of comparable density estimates away from roads, which are not practical to collect, we cannot evaluate this possible effect. However, we note that people on foot or in vehicles do not usually attempt to kill or disturb vultures in India and the birds are quite tame and appear not to be afraid of humans and built infrastructure. Hence, we suggest that it is unlikely that there was underestimation of population size due to vultures avoiding roads because of fear.

A third caveat is that our regression analysis did not allow for possible effects of spatial autocorrelation for technical reasons (see Methods). However, we consider that this is unlikely to have had a large effect on the regression results or the population estimates based upon them because the regressions included quadratic effects of latitude and longitude and the degree of spatial autocorrelation in the density residuals from the fitted model was slight.

Before diclofenac came into widespread veterinary use in India, the millions of tonnes of carrion from cattle carcasses discarded annually provided a safe and widely-distributed food supply for vultures, in addition to the less plentiful carcasses of wild ungulates. Spatial variation in the occurrence of wild ungulates in India is positively correlated with forest cover and, additionally, with the presence of protected areas. This indicates that both the area of natural habitats and protection from hunting have important effects on wild ungulates (Karanth et al. 2009). It seems likely that vultures have declined less in and near National Parks than far from them at least partly because a greater proportion of the food of birds foraging to some extent in National Parks consists of carcasses of wild ungulates that are more abundant there than outside and are never contaminated with NSAIDs. In addition, the health hazard and nuisance arising from cattle carcasses not being rapidly eaten by vultures has led to a proportion of them being disposed of by methods such as burial. This may have resulted in carrion from wild ungulates now being a larger proportion of the total available food supply than it was before the vultures declined.

However, populations of vultures living in and near National Parks have also declined, though not by as much as those elsewhere (Prakash et al. 2012). Vultures range over long distances from their breeding and roosting sites whilst foraging. Gilbert et al. (2007) found that five adult male White-rumped Vulture, satellite tagged in Pakistan, ranged up to 316 km from their breeding or roosting sites (mean convex polygon range area 24,155 km²), even though plentiful supplementary food was provided near these sites during part of the period. This mean foraging range is about fifty times the mean area of National Parks in India (490 km²) and seven times larger than the largest park. Hence, there is likely to be a risk of exposure to diclofenac for vultures breeding in National Parks from contaminated carcasses of domesticated ungulates well beyond their boundaries, even though feeding from carcasses of uncontaminated wild ungulates in the parks may reduce it. In addition, Gyps vultures may be more numerous in National Parks because of the greater availability of nesting and roosting sites, such as trees or cliffs, in the relatively undisturbed forests, woodlands and mountains within the parks.

Although we found that vulture densities in 2003–2015 were higher near to National Parks than distant from them, we did not find a similar effect of proximity to Wildlife Sanctuaries. If the explanation of the effect of proximity to National Parks is the safe food supply provided by carcasses of wild ungulates (see above), it might be that the abundance of wild ungulates is lower in Wildlife Sanctuaries than in National Parks leading to a smaller and undetectable effect on the level of exposure of the vultures to diclofenac. Densities of wild ungulates per unit area of natural habitat in a sample of 11 protected areas in India were found to vary by more than a factor of 10, with differences in the level of protection of ungulates from hunting being one of the principal variables affecting density (Karanth et al. 2004). If Wildlife Sanctuaries tend to have lower, and perhaps more variable, levels of protection of wild ungulates than National Parks, this might account for our failure to find robust evidence for an effect on vulture density of proximity to Wildlife Sanctuaries.

The future persistence of wild Gyps vulture populations in India will depend upon effective implementation of the existing regulation of the veterinary use of diclofenac and measures to prevent the use of other veterinary drugs with similar effects. However, our findings also imply that measures to maintain or improve the effectiveness of the protection of wild ungulate populations and habitats within National Parks and Wildlife Sanctuaries also have a role to play in slowing or reversing vulture declines.