Introduction
Populations of three species of vultures endemic to South and South-East Asia, the White-rumped Vulture Gyps bengalensis, Indian Vulture Gyps indicus and Slender-billed Vulture Gyps tenuirostris began to decline rapidly in the 1990s, leading to them to be classified as ‘Critically Endangered’ in the IUCN Red List in 2000; a status which has been maintained in all subsequent assessments (BirdLife International 2017). Until the beginning of these declines, the majority of the large populations of each of these three species were found in the Indian subcontinent. The most abundant of them, the White-rumped Vulture, was described in 1985 as possibly the most abundant large bird of prey in the world, with a global population of tens of millions (Houston Reference Houston, Newton and Chancellor1985). However, the population of this species in India declined by a factor of about 1,000 between 1992 and 2007, at rates of population decline of about 50% per year (Prakash et al. Reference Prakash, Green, Pain, Ranade, Saravanan, Prakash, Venkitachalam, Cuthbert, Rahmani and Cunningham2007). Indian and Slender-billed Vulture populations in India also suffered large declines at rates of about 35% per year (Prakash et al. Reference Prakash, Galligan, Chakraborty, Dave, Kulkarni, Prakash, Shringarpure, Ranade and Green2017). Declines of similar magnitude, rate and timing also occurred in the formerly large White-rumped and Indian Vulture populations in Pakistan (Gilbert et al. Reference Gilbert, Virani, Watson, Oaks, Benson, Khan, Ahmed, Chaudhry, Arshad, Mahmood and Shah2002) and there were also declines in White-rumped and Slender-billed Vulture populations during the same period in Nepal (Chaudhary et al. Reference Chaudhary, Subedi, Giri, Baral, Bidari, Subedi, Chaudhary, Chaudhary, Paudel and Cuthbert2012) and Bangladesh (Khan Reference Khan2013). Since these declines began in South Asia, evidence has also emerged of a wide range of other serious threats to Old World vultures in most parts of their global range (Safford et al. Reference Safford, Andevski, Botha and Bowden2019).
There is a large body of evidence that veterinary use on cattle of the non-steroidal anti-inflammatory drug (NSAID) diclofenac was the main cause of the rapid vulture population declines in the Indian subcontinent (Oaks et al. Reference Oaks, Gilbert, Virani, Watson, Meteyer, Rideout, Shivaprasad, Ahmed, Chaudry, Arshad, Mahmood, Ali and Khan2004, Green et al. Reference Green, Newton, Shultz, Cunningham, Gilbert, Pain and Prakash2004, Pain et al. Reference Pain, Bowden, Cunningham, Cuthbert, Das, Gilbert, Jakati, Jhala, Khan, Naidoo, Oaks, Parry-Jones, Prakash, Rahmani, Ranade, Baral, Senacha, Saravanan, Shah, Swan, Swarup, Taggart, Watson, Virani, Wolter and Green2008). Diclofenac causes kidney failure and death in vultures within a few days of them consuming tissue from the carcass of an ungulate that has died within a few hours or days after treatment with the drug (Oaks et al. Reference Oaks, Gilbert, Virani, Watson, Meteyer, Rideout, Shivaprasad, Ahmed, Chaudry, Arshad, Mahmood, Ali and Khan2004, Green et al. Reference Green, Taggart, Das, Pain, Sashi Kumar, Cunningham and Cuthbert2006, Swan et al. Reference Swan, Cuthbert, Quevedo, Green, Pain, Bartels, Cunningham, Duncan, Meharg, Oaks, Parry-Jones, Shultz, Taggart, Verdoorn and Wolter2006a). In the early 2000s, 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. Reference Green, Taggart, Senacha, Raghavan, Pain, Jhala and Cuthbert2007).
In response to evidence about the impact of veterinary diclofenac on vultures, the licence to manufacture veterinary formulations of diclofenac was withdrawn by the Governments of India, Nepal and Pakistan in 2006 and in Bangladesh in 2010. Repeated surveys of the prevalence and concentration of diclofenac in tissues from carcasses of domestic ungulates available to vultures in India after the ban showed that the use of the drug had declined, but that it was still present in a substantial proportion of carcasses in 2009, three years after the ban (Cuthbert et al. Reference Cuthbert, Taggart, Prakash, Chakraborty, Deori, Galligan, Kulkarni, Ranade, Saini, Sharma, Shringarpure and Green2014). Modelling based upon these domestic ungulate survey data indicated that the expected risk of death from diclofenac poisoning per meal for White-rumped Vultures had fallen markedly, but only partially, by 2009 in India; then being about one-third of its 2006 level (Cuthbert et al. Reference Cuthbert, Taggart, Prakash, Chakraborty, Deori, Galligan, Kulkarni, Ranade, Saini, Sharma, Shringarpure and Green2014). Post mortems and tissue analyses of Gyps vultures found dead in the wild in India after the ban showed that they continued to die from diclofenac poisoning, though probably at a lower rate than before it (Cuthbert et al. Reference Cuthbert, Taggart, Mohini, Sharma, Das, Kulkarni, Deori, Ranade, Shringarpure, Galligan and Green2016). Diclofenac formulated for human use and packaged in large vials was being sold illegally for veterinary use in India (Cuthbert et al. Reference Cuthbert, Dave, Chakraborty, Kumar, Prakash, Ranade and Prakash2011). In response to this misuse of human diclofenac, the Government of India banned the manufacture of all diclofenac products in vial sizes containing more than the volume of a single dose for humans (3 ml) in 2015. This was expected to make the drug too expensive and inconvenient to use on large-bodied animals and therefore reduce its misuse in livestock. These changes in India were also expected to affect neighbouring countries, including Nepal, because much of the diclofenac used there is imported from India. Since these actions were taken, Gyps vulture declines in India have slowed or ceased (Prakash et al. Reference Prakash, Bishwakarma, Chaudhary, Cuthbert, Dave, Kulkarni, Kumar, Paudel, Ranade, Shringarpure and Green2012), but there is no indication yet of measurable recovery there (Prakash et al. Reference Prakash, Galligan, Chakraborty, Dave, Kulkarni, Prakash, Shringarpure, Ranade and Green2017).
Although the introduction of drug-control regulations is an important component of vulture conservation programmes in South Asia, it was recognised that the bans were unlikely to solve the problem quickly on their own. In response, a network of Vulture Safe Zones (VSZ) has been established in South Asia (Mukherjee et al. Reference Mukherjee, Galligan, Prakash, Paudel, Khan, Prakash, Ranade, Shastri, Dave, Donald and Bowden2014). Each VSZ encompasses one or more nesting colonies of ‘Critically Endangered’ Gyps vultures and > 30,000 km2 of the foraging area surrounding them. VSZs involve communities in advocacy to conserve vultures and awareness-raising about the hazard posed by diclofenac and the availability of the vulture-safe alternative drug meloxicam (Swan et al. Reference Swan, Naidoo, Cuthbert, Green, Pain, Swarup, Prakash, Taggart, Bekker, Das, Diekmann, Diekmann, Killian, Meharg, Patra, Saini and Wolter2006b, Swarup et al. Reference Swarup, Patra, Prakash, Cuthbert, Das, Avari, Pain, Green, Sharma, Saini, Das and Taggart2007). VSZs also include monitoring of the availability and use of NSAIDs among farmers, veterinarians and pharmacists. By 2010–2011, Nepal’s VSZs had become effective in reducing the apparent availability of diclofenac to the owners of livestock (Prakash et al. Reference Prakash, Bishwakarma, Chaudhary, Cuthbert, Dave, Kulkarni, Kumar, Paudel, Ranade, Shringarpure and Green2012). However, it was unclear from the count data available at that time whether vulture populations had ceased to decline (Chaudhary et al. Reference Chaudhary, Subedi, Giri, Baral, Bidari, Subedi, Chaudhary, Chaudhary, Paudel and Cuthbert2012, Prakash et al. Reference Prakash, Bishwakarma, Chaudhary, Cuthbert, Dave, Kulkarni, Kumar, Paudel, Ranade, Shringarpure and Green2012). In this paper, we report a longer series of road transect results for White-rumped and Slender-billed Vultures, covering a large proportion of the range of these species in Nepal. We also present results of further open and undercover surveys of pharmaceuticals on sale in pharmacies for veterinary use. We compare the results of vulture surveys with those from the pharmacy surveys to assess the degree to which changes in population trend and the availability of the toxic NSAID diclofenac coincided. Our results indicate that the population of White-rumped Vultures in Nepal ceased to decline in the early 2010s and has partially recovered at as rapid a rate as could be expected. Similar changes may also be occurring in the Slender-billed Vulture population. These changes in population trend coincide approximately with the time when diclofenac availability fell to a low level.
Methods
Study area
We surveyed vultures in the Terai lowlands and Pahad hills of Nepal (Figure 1). This area lies between the Himalayan Mountains to the north and the Gangetic Plain in India to the south. Both the Terai and Pahad zones consist of fragmented patches of grassland, forest and urban areas within a matrix of agricultural land. The Terai and Pahad hold the largest densities of wild and domesticated ungulates in Nepal and therefore provide suitable foraging areas for vultures (Baral et al. Reference Baral, Giri, Virani, Chancellor and Meyburg2004). Two of the three ‘Critically Endangered’ Gyps vulture species endemic to South and South-East Asia, the White-rumped Vulture and Slender-billed Vulture, nest in these regions of Nepal. These species both nest in trees, of which considerable numbers are widely distributed in these regions. Surveys were conducted in four regions whose boundaries are shown in Figure 1. We refer to these as the Eastern Terai, Western Terai, Western Pahad and Far Western Pahad.
Figure 1. Map of Nepal showing the road transects surveyed (solid black lines) within the four survey regions (delineated by dashed white lines; Far Western Pahad (FWP); Western Terai (WT); Western Pahad (WP); and Eastern Terai (ET)). Also shown are the Vulture Safe Zone (dark grey area), Diclofenac-free District Zone (light grey area) and the area outside these vulture conservation zones (white area).
Deployment of VSZs, diclofenac-free districts; and vulture safe feeding sites in Nepal
The first VSZ in South Asia was initiated formally in the Western Terai lowlands of Nepal in 2009, but built upon advocacy and awareness-raising activities begun in 2006. The VSZ network in Nepal gradually expanded to include the entire Western Terai and Western Pahad mid-hills regions by 2010 (Figure 1). This VSZ network was developed and implemented alongside two other initiatives intended to conserve vultures: Diclofenac-free Districts; and Vulture Safe Feeding Sites. In Diclofenac-free Districts, district-level decision-makers pledge to make the food supply of vultures in their district free from diclofenac. This is often a preliminary step which prepares areas outside the existing VSZ network to be declared as additional VSZs. More than 80% of districts in Nepal are now VSZs or Diclofenac-free Districts (Figure 1). Vulture Safe Feeding Sites provide vultures with access at undisturbed locations to cattle carcasses free from diclofenac and other toxic NSAIDs. This is usually achieved by gathering old and ailing cattle from the surrounding region and caring for them without the use of NSAIDs until they die a natural death. Six Vulture Safe Feeding Sites are located inside VSZs, with one being in a Diclofenac-free District (Figure 2). Within the historical range of White-rumped and Slender-billed Vulture in Nepal, only parts of the Eastern Terai have had little or no vulture conservation activity so far.
Figure 2. Map of Nepal showing all sightings of White-rumped Vultures (white circles) and Slender-billed Vultures (black circles) recorded during annual road transects surveys between 2008 and 2018. Each circle represents one or more individuals. The black and white stars indicate where Nepal’s seven Vulture Safe Feeding Sites are located. Also shown are the Vulture Safe Zone (dark grey area), Diclofenac-free District Zone (light grey area) and the area outside these vulture conservation zones (white area).
Vulture survey method
We used Nepal’s East-West Highway to survey vultures. The Highway runs continuously through the Terai from the eastern to the western borders of Nepal with India (Figure 1). Because White-rumped Vultures are most abundant in the Western Terai (Baral et al. 2014), we surveyed most often on the western section of the Highway, which runs for 638 km through the districts of Nawalparasi, Rupandehi, Kapilbastu, Dang, Banke, Bardia, Kailali and Kapilbastu. We also surveyed less frequently from the eastern section of the Highway, which runs for 470 km through the districts of Chitwan, Makwanpur, Bara, Rautahat, Sarlahi, Mahottari, Dhanusa, Siraha, Saptari, Sunsari, Morang and Jhapa. Arterial roads lead off the Highway, heading north into the Pahad and then running parallel to the Highway. In the most recent surveys, we also surveyed an additional 1,108 km of these arterial roads through the districts of Salyan, Rukum, Rolpa, Pyuthan, Arghakhanchi, Gulmi, Palpa, Syanja, Kaski and Tanahu (Figure 1). The same roads were traversed in each survey in each of the regions mapped in Figure 1.
We undertook surveys in 14 of the 17 years between 2002 and 2018, but only covered the whole of the road network described above in 2012. We surveyed the Western Terai in 2002, 2003, 2006 and 2008–2018; the Eastern Terai in 2002, 2003, 2009, 2012, 2015 and 2018; the Western Pahad in 2012–2014; and the Far Western Pahad 2012, 2013 and 2014. In each survey, we covered regions in the same order: Eastern Terai, Western Terai, Far Western Pahad and Western Pahad. Surveys began in early May, except in 2002 and 2003, when the Eastern Terai was surveyed in early April and the Western Terai early May. April and May are the first and second months after the end of the breeding season of White-rumped Vultures (October–March), when numbers of this species are largest because of recruitment of juveniles. In most years, the survey took 12–19 consecutive days: (3-4 days in Eastern Terai; 4-5 days in Western Terai, 5-7 days in Far Western and Western Pahad and up to 3 non-survey days travelling between regions). However, civil unrest sometimes caused disruption of the survey schedule and lengthened the number of days required to complete the survey. There were 2-day interruptions in 2008 and 2011; and a 14-day interruption in 2012. In addition, there was an interval between survey stages of 21 and 30 days in 2002 and 2003, respectively (see above).
Surveys were conducted between 07h00 and 17h00. Two surveyors stood, facing opposite directions, in the cargo bay of a utility truck, which was driven at 20 km h-1. Surveyors scanned continuously for vultures, which they then viewed with 8 x 42 or 10 x 42 binoculars. In most years, all transects in the year were undertaken by the same surveyors. During the 15 years of this study, a total of nine surveyors was involved. In every year, there was at least one surveyor who had done the survey before. All surveyors were experienced in spotting vultures and identifying them to species. The driver did not count vultures. On spotting one or more vultures, the truck was stopped and the surveyors both counted the vultures and then agreed on their number and species. All vultures seen perching, flying and soaring were counted. Geographic co-ordinates of vulture sightings were recorded using a hand-held GPS device, except in 2002, 2003 and 2006.
Accurate estimation of the distance of vultures from the observers or transect line was considered to be impractical, so no distance measurements were attempted. Surveyors were asked to count vultures seen within 1 km on either side of the road. The degree to which this survey strip width was adhered to in practice is uncertain because of the lack of an accurate method for measuring distances, but we believe that mean detection distances are unlikely to have changed systematically during the survey period.
Surveys of NSAIDs offered for sale in veterinary pharmacies
Surveys of veterinary pharmacies were conducted in the Western Terai region, within the VSZ (Figure 1). Initial surveys by Bird Conservation Nepal (BCN) were undertaken in December 2006–January 2007 following the ban on veterinary diclofenac. In these surveys, the total number of pharmacies visited and the total stock (30 ml vials and equivalent strips of bolus formulations) of diclofenac or meloxicam per district were recorded. These surveys were repeated in July 2007, November–December 2007, March 2008, February–March 2010, April–May 2010, June–July 2010, May–June 2011 and December 2011. Not all districts surveyed in the first survey were surveyed again in each of these periods, but they were surveyed at least twice. The mean stock level, as items per pharmacy, in the districts surveyed was used a measure of the availability of these two drugs. Additional surveys by BCN were begun in 2007 in some districts, but in 2008 or 2009 in others, according to when Diclofenac-free District activities began in these districts. In these additional surveys, detailed data was collected from pharmacies, including the shop name, shop address, types of NSAIDs stocked, but not the quantity of stock. All districts were surveyed again in 2010 and 2011. Therefore, for analysis, we grouped data from the initial surveys (a single period 2007–2009) and compared these to data from the two later surveys. BCN staff and a district livestock official visited veterinary pharmacies and asked which NSAIDs were available for the treatment of cattle. Surveyors made it clear to the pharmacy staff that the visit was connected with vulture conservation and then discouraged veterinary use of diclofenac and encouraged veterinary use of meloxicam. Pharmacies were selected for survey opportunistically but were widely distributed within the region. Some pharmacies were visited in more than one survey period, but records are not available to indicate which ones were surveyed repeatedly.
In addition to these open surveys, in which pharmacy staff knew that the surveyors’ visit was connected with vulture conservation, there were also undercover surveys, designed with the intention that pharmacy staff would be unaware of the purpose of the visit. However, undercover surveys did not begin until 2012, so no data from them are available for the first six years of the implementation of the VSZ and Diclofenac-free District programme. Between 2012 and 2016 four sets of surveys of veterinary pharmacies were conducted in the Western Terai region, with the sets being spaced at approximately annual intervals. In each set, the survey team visited settlements selected at random from those mapped for the region, dressed in ordinary clothes and travelling in an unmarked vehicle. In each settlement, a local man was asked to purchase medicines from one veterinary pharmacy in exchange for a small payment. He was told that the medicines were for a sick cow and given adequate money to purchase the drugs. Vulture conservation was not discussed with him. The team took the man to the nearest pharmacy selling animal medicines and gave him a basic prescription written by a veterinarian for one 30 ml vial of diclofenac and five tablets of a deworming medication (Matemar). He was told that the prescription was written by an unqualified veterinarian, of which there are many in Nepal, and that he should purchase the first drugs he was offered, regardless of whether they were diclofenac and Matemar or not. On the day after this undercover survey, the team themselves conducted an open survey of several randomly selected veterinary pharmacies in each settlement, which sometimes included the same pharmacy as that surveyed undercover. During the open surveys, the surveyors wore clothes with vulture conservation logos and slogans, distributed vulture conservation posters and pamphlets and talked to pharmacy staff about the plight of vultures and the ban on diclofenac. They asked to see all veterinary NSAID products stocked and recorded the products they were shown. The reason for conducting the open survey in addition to the undercover survey was to allow direct comparison of the results from the two methods. This was useful because only open survey data were available until 2012.
Statistical analysis
We were unable to estimate population densities of vultures or allow for possible differences in detection of birds between regions, years or observers using a package such as DISTANCE (Buckland et al. Reference Buckland, Anderson, Burnham, Laake, Borchers and Thomas2001) because we lacked estimates of the distances of vulture records from the transect line (see above). We therefore fitted generalised linear models to the raw count data in Table 1. In all models, the count of vultures in a given year and region was the dependent variable. We fitted a Poisson model first and used a quasi-Poisson model instead if the overdispersion factor (residual deviance divided by residual degrees of freedom) for the Poisson model was > 2 (Crawley Reference Crawley2007, Thomas Reference Thomas2017). We also considered the use of zero-inflated Poisson regression, but did not pursue this option because there was only one count of zero in the data for White-rumped Vulture, which was the main focus of our analysis because of its more substantial dataset (Table 1). Although the data for Slender-billed Vulture included a higher proportion of zero counts, some of those were not included in the model fitting because it was necessary to include region as an explanatory variable in all models (see below). Hence, data for the two regions (Eastern Terai and Far Western Pahad) in which Slender-billed Vultures were never recorded on the transects, were not informative about changes in population of this species and were excluded. Eight of the 21 remaining counts in the other two regions were zeroes, which might justify the use of zero-inflated Poisson models. However, we considered that the dataset was too sparse to support their use, because the models require the estimation of additional parameters. Furthermore, Poisson models performed reasonably well for the Slender-billed Vulture data (see Results).
Table 1. Counts of White-rumped and Slender-billed Vultures in 2002–2018 on road transect surveys in four regions of Nepal.
We used two model structures to describe changes over time in an index of population. The first model structure was log-linear with fixed effects of region and year, both of which were modelled as categorical variables (factors). This model produces annual indices of population change which are assumed to be the same across all regions. Because not all regions were surveyed in every year, it was necessary to include region in the model as a factor to allow the estimation of a single annual population index for each year. Region was modelled as a fixed, rather than random, effect because of the small number of regions (four for White-rumped Vulture and two for Slender-billed Vulture) and because region effects would not necessarily be expected to be log-normally distributed. The model has a logarithmic link function, so the vulture count 
${C_{ij}}$ in region i and year j is modelled as,
$${\log _e}\left( {{C_{ij}}} \right) = {r_i} + {y_j},$$where the r i and y j are constants, representing the effects of region and year respectively, which are estimated so as to maximise the likelihood of the data under the model. The model is equivalent to
$${C_{ij}} = {e^{ri}}{e^{yi}},$$with the eyj being annual estimates of the population index relative to the population in the first year of the series (2002 in this case).
The second model structure was also log-linear, but with the effect of time being modelled as continuous rather than categorical. In this case, the model is
$$\log _e \left( {C_{ij} } \right) = r_i + b_1 t,$$ where t is the time elapsed, in years, since 2002. The fitted parameter b 1 is an estimate of the mean population growth rate, and 
${e^{b1}}$ estimates the mean population multiplication rate λ. Inspection of the data in Table 1 suggested an abrupt change in population trend during the study period, so we also fitted a piecewise version of the log-linear model as
$$\log _e \left( {C_{ij} } \right) = r_i + b_1 t + b_2 \left( {t - t^* } \right),$$where t* is the estimated breakpoint in years since 2002, at which time the mean population growth rate changed, b 1 is the mean population growth rate before the change and b 2 is the difference in mean population growth rate after the change, compared with that before it. The value of the variable (t-t*) was set to zero when (t-t*) < 0. We identified the maximum-likelihood value of the breakpoint t* using an iterative bisection search algorithm and determined the 95% confidence limits of this estimated t* value from its profile likelihood (Kalbfleisch Reference Kalbfleisch1979). We performed an F test of whether the piecewise log-linear model with a breakpoint performed significantly better than the model without a breakpoint. All regression models were implemented using the function glm in the package R, version 3.5.1 (R Development Core Team 2018).
We wished to compare the population multiplication rates λ, estimated from the regression analyses described above, with the maximum population multiplication rate λ max expected to be possible for vultures, given the most favourable likely values of their demographic rates. To do this, we used the Demographic Invariants Method (DIM) proposed by Niel and Lebreton (Reference Niel and Lebreton2005). The DIM requires values for the age at first reproduction α and adult survival probability s, both of which should be assessed under optimal conditions. Age of first reproduction has not been measured for wild Asian vultures, so we used values for the Eurasian Griffon Vulture Gyps fulvus of α = 4 years from Ferrière et al. (Reference Ferrière, Sarrazin, Legendre and Baron1996) and α = 5 years from Cramp and Simmons (Reference Cramp and Simmons1980). Green et al. (Reference Green, Newton, Shultz, Cunningham, Gilbert, Pain and Prakash2004) estimated adult survival s for Asian Gyps vultures when they are not impacted by the effects of diclofenac, by using combined data on White-rumped and Indian Vultures in India and Pakistan, as presented in their Figure 2(a). Their best estimate was s = 0.907. We used the plot of residual deviance in relation to adult survival shown in Figure 2(c) of Green et al. (Reference Green, Newton, Shultz, Cunningham, Gilbert, Pain and Prakash2004) and the profile likelihood method of Kalbfleisch (Reference Kalbfleisch1979) to estimate the 95% confidence interval of this estimate of s (0.867–0.937). We used the same values of α and s for both White-rumped and Slender-billed Vulture and estimated λ max from them using Niel and Lebreton’s equation 17. We calculated the 95% confidence interval for λ max from the confidence limits for s. We had no estimate of uncertainty for α and therefore did not allow for its effect on the confidence interval for λ max, except by calculating λ max for both of the alternative values of α. Noting that the 95% confidence intervals of loge(λ max) obtained by this method were symmetrical about the maximum-likelihood value of loge(λ max), we used a parametric bootstrap method, with 10,000 replicates and the estimated means and standard errors of loge(λ), to perform a one-tailed significance test of the difference between loge(λ) from the piecewise regression analysis of vulture counts and loge(λ max) from the DIM calculations.
Results
Distribution and counts of vultures
The summed total of all counts of the White-rumped Vulture was 1,632 whilst that of Slender-billed Vulture was much lower at 30 (Table 1). Inspection of the annual counts of White-rumped Vultures indicates a similar tendency across regions for numbers to decrease during the 2000s, followed by an increase in the next decade in the Eastern Terai, Western Terai and Western Pahad. In the Far Western Pahad, counts were only conducted in 2012–2014 and only a few White-rumped Vultures were counted there, so no long-term pattern can be detected in that region on its own. A map of White-rumped Vulture locations on transects conducted in 2008 – 2018 indicates that they were concentrated in the Western Terai and Western Pahad regions, where five of the seven Vulture-Safe Feeding Sites are located (Figure 2). All but one of the records of this species outside this region were close to the other two Vulture Safe Feeding Sites (Figure 2). Counts of Slender-billed Vultures also suggested a tendency to decrease over time, followed by an increase, with birds being detected in the period 2002–2009 and again in 2013–2018, but with none counted in the 2010, 2011 or 2012 surveys (Table 1). However, the data for this species are too sparse for any firm conclusions to be drawn. A map of Slender-billed Vulture locations on transects conducted in 2008–2018 indicates that they were concentrated in the Western Terai and Western Pahad regions with most records being close to Vulture Safe Feeding Sites (Figure 2).
Changes in vulture population index
The Poisson regression with year modelled as a categorical variable, used to estimate annual population index values from the White-rumped Vulture counts in Table 1, indicated substantial overdispersion (overdispersion factor = 3.534), so we used a quasi-Poisson model instead. This model gave the same index estimates as the Poisson, but the precision of the index values was lower (the standard errors of the y i values were larger). The dispersion parameter θ of the quasi-model was 3.259. When plotted against time, the index values showed a tendency to decline from 2002 until about 2013 and then to increase until 2018 (Figure 3).
Figure 3. Annual index values for populations of White-rumped (WRV: filled circles) and Slender-billed (SBV: open circles) Vultures in Nepal for 2002–2018, relative to 2002. Vertical lines for the WRV points are 95% confidence intervals from the quasi-Poisson model. Curves show results from the fitted piecewise log-linear regression models for WRV (solid line) and SBV (dashed line). Crosses in the upper part of the diagram show the estimated breakpoints and their 95% confidence intervals.
The Poisson regression model used to estimate annual population index values from the Slender-billed Vulture counts did not indicate overdispersion (overdispersion factor = 0.815). The data were so sparse that most of the standard errors of the index values for individual years were large and the individual values therefore did not represent robust estimates of population relative to that in 2002. Hence, we do not show confidence intervals for index values for this species in Figure 3. However, when plotted against time, the index values show a broadly similar pattern to that indicated more robustly for the White-rumped Vulture, with a tendency to decline and then to increase (Figure 3).
Changes in vulture population trend
The piecewise Poisson regression model used to estimate population growth rates and a breakpoint from the White-rumped Vulture counts indicated substantial overdispersion (overdispersion factor = 3.315), so we used a quasi-Poisson model instead. The dispersion parameter θ of the quasi-model was 3.109. The effect of including a breakpoint in the model was highly significant (F 1,23 = 76.25; P < 0.0001). The estimated breakpoint was in March 2013, with 95% confidence limits from January 2012 to January 2014 (Figure 3). In the period before the breakpoint, the estimated mean population multiplication rate λ was 0.863 (annual decline of 14% per year) with 95% confidence limits for λ of 0.840 to 0.886. In the period after the breakpoint, the estimated mean population multiplication rate λ was 1.217 (annual increase of 22% per year) with 95% confidence limits for λ of 1.148 to 1.289. Visual comparison of the curves fitted using the piecewise regression models with the annual population index values from models with survey year as a categorical variable indicate that the piecewise models describe the patterns in the data reasonably well for both species (Figure 3).
The piecewise Poisson regression model used to estimate population growth rates and a breakpoint from the Slender-billed Vulture counts did not indicate overdispersion (overdispersion factor = 0.969), so we used it as the final piecewise model and did not use a quasi-Poisson model. The effect of including a breakpoint in the model was significant (F 1,16 = 13.15; P = 0.002. The estimated breakpoint was in June 2012, with wide 95% confidence limits from February 2010 to April 2017 (Figure 3). In the period before the breakpoint, the estimated mean population multiplication rate λ was 0.752 (annual decline of 25% per year) with 95% confidence limits for λ of 0.659 to 0.857. In the period after the breakpoint, the estimated mean population multiplication rate λ was 1.413 (annual increase of 41% per year) with 95% confidence limits for λ of 1.118 to 1.786. 1 Hence, although the count data for Slender-billed Vultures are too sparse for any reliable conclusions to be drawn, our analyses suggest a broadly similar pattern of population change, in terms of timings and rates of change, to that for the White-rumped Vulture.
Comparison of recent vulture population multiplication rates with the maximum expected rate
The mean population multiplication rates λ of White-rumped and Slender-billed Vultures, estimated for recent periods of apparent increase, were somewhat higher than the expected maximum population multiplication rate λ max under optimal conditions obtained by the DIM method (Table 2). The comparison is most meaningful for the White-rumped Vulture because λ during the 2010s was estimated with greater precision for that species. Our best estimate of the recent mean rate of increase for the White-rumped Vulture was 22% per year, which is considerably more rapid than the expected maximum possible rates of 11–13% and significantly more rapid than the DIM estimates both for α = 4 years (P = 0.011) and α = 5 years (P = 0.002). The estimated recent mean rate of increase for the Slender-billed Vulture was 41% per year, which was also significantly more rapid than the DIM estimates for α = 4 years (P = 0.031) and α = 5 years (P = 0.023).
Table 2. Comparison between mean annual population multiplication rates λ of White-rumped and Slender-billed Vultures in Nepal, estimated for recent periods of apparent increase, and the expected maximum population multiplication rate λ max calculated by the Demographic Invariants Method (see text). λ max is shown for two alternative values of the age at first reproduction α.
Availability of diclofenac and meloxicam in veterinary pharmacy shops
The first open survey of veterinary pharmacies covered 1184 pharmacies across 9 districts in the Western Terai. The mean stock of diclofenac decreased by 99.9% between December 2006–January 2007 and July 2007; and, in survey periods between December 2007 and December 2011, the mean stock of diclofenac was 0–2 items per pharmacy per district (Figure 4). The mean stock of meloxicam found in open surveys fluctuated between 23 and 122 items per pharmacy per district, but increased by 99.8% between June–July 2010 and May–June 2011 and remained high in December 2011 (Figure 4). The proportion of veterinary pharmacies in the Western Terai region which offered diclofenac for veterinary use in open surveys fell from 24% in 2007–2009 to 1.6% in 2010. No pharmacies offered diclofenac in 2011, as was also the case for the open surveys in the period 2012–2016 (Table 3). A high proportion (73%) of shops offered the vulture-safe alternative NSAID meloxicam in open surveys during 2007–2009, but the proportions were even higher (93% and 87%) in 2010 and 2011. The proportion of veterinary pharmacies offering diclofenac was also low (2.3%) in the undercover surveys in 2012–2016, although slightly higher than in open surveys on the following day of shops in the same settlements (0%; Table 3). The proportion of shops offering meloxicam was high during 2012–2016 and was similar for the undercover and open surveys (93.5% and 93.2% respectively; Table 3).
Figure 4. Mean (± 1 SE) stock per pharmacy shop of the NSAIDs (A) diclofenac and (B) meloxicam in items per shop for areas of the Western Terai region surveyed at various times (months) since the ban in diclofenac was introduced in December 2006. Results are for open surveys only. The ordinate is on a logarithmic scale.
Table 3. Proportion of veterinary pharmacy survey shops in the Western Terai region at which diclofenac and meloxicam were offered for sale during open and undercover surveys.
Discussion
Our results indicate that the rapid declines of White-rumped and Slender-billed Vulture populations observed in Nepal in the 2000s (Chaudhary et al. Reference Chaudhary, Subedi, Giri, Baral, Bidari, Subedi, Chaudhary, Chaudhary, Paudel and Cuthbert2012) ceased in the early 2010s. A previous report that the decline in the White-rumped Vulture population in western Nepal had slowed or ceased during 2009–2011 was based on a suggestive, but non-significant, regression analysis (Prakash et al. Reference Prakash, Bishwakarma, Chaudhary, Cuthbert, Dave, Kulkarni, Kumar, Paudel, Ranade, Shringarpure and Green2012). Our analysis of counts from more areas of Nepal and over a longer period, provide more robust evidence that a change has occurred. For the White-rumped Vulture there is now strong evidence from our piecewise regression analysis that a decline up to about 2013 has given way to a rapid increase from about 2013 to 2018. The population count data for Slender-billed Vulture are too sparse to draw firm conclusions, but they nonetheless suggest that a decline gave way to increase or stability from about 2012 onwards. For both species, the mean annual rate of population increase estimated since the early 2010s was significantly greater than the maximum rate of increase calculated according to the DIM method of Niel and Lebreton (Reference Niel and Lebreton2005). This suggests that vulture populations may now be benefiting substantially from reduced competition for food or other resources arising from their low population size, which might also be attracting some birds from adjacent areas of India.
Veterinary use of diclofenac was banned in Nepal in 2006. If this had resulted in an immediate removal of diclofenac from all carrion available to vultures, the decline in vulture populations would have been expected to cease in the years immediately after this action. In contrast, it took a further 6–7 years after the ban was introduced before there was evidence of a change in population trend. However, when the change did occur, in 2012–2013, it was very marked. No direct measurements have been made of the prevalence and concentration of diclofenac in ungulate carcasses available to vultures in Nepal before and after veterinary use of the drug was banned. We therefore cannot assess whether, or how rapidly, reduced diclofenac contamination of the vultures’ food supply may have contributed to the apparent change in population trend. However, information on the timing of implementation of the community-based Vulture Safe Zone programme suggests that may have taken several years to become fully effective.
Widespread vulture conservation advocacy and awareness activities for vultures, as well as a programme of work to ensure the removal of veterinary stocks of diclofenac from pharmacies was initiated in several lowland districts in Nepal in 2006 and 2007 during the development of the Diclofenac-free District and VSZ programme. These activities spread to a programme across all western districts during 2008–2010 and have spread further since. An initial open survey of veterinary pharmacies found a sudden large reduction in diclofenac stock in the areas covered by advocacy work immediately after the ban in December 2006 and very little or no diclofenac stock present between December 2007 and December 2011. Another open survey of veterinary pharmacies undertaken by vulture conservationists accompanied by a district livestock official found nearly one-quarter of shops surveyed had diclofenac during the period 2007–2009 and very little and no diclofenac stock in 2010 and 2011. We do not suggest that these open pharmacy survey data necessarily reflect the true timing of changes in diclofenac availability to livestock owners in the entire range of these vulture species in Nepal because pharmacy staff had almost certainly become aware that veterinary use of diclofenac was being discouraged by state authorities, non-government conservation bodies and community leaders. Pharmacy staff may therefore have concealed the availability and stock of diclofenac from the surveyors, particularly in the initial survey when conservationists undertook surveys without district livestock officials, and done so to an increasing extent as the surveys proceeded. In addition, the area covered by of the Diclofenac-free District and VSZ programme expanded over time and monitoring of shops was concentrated in the areas where advocacy was conducted. Hence, the apparent decline in the availability of diclofenac in these areas may not have been representative of that in the vulture range in Nepal as a whole. However, undercover surveys of veterinary pharmacies in the Western Terai region in 2012 - 2016, during which pharmacy staff were not aware that the potential purchaser was associated with vulture conservation work, found that diclofenac was offered for sale only in a slightly higher proportion of shops (2.3%) than for open surveys conducted in the same settlements during the same period. Comparing the results from the open advocacy visits in 2007–2009 when 24% of shops offered diclofenac, but some may have concealed it, with undercover surveys in 2012–2016, when 2.3% of shops offered diclofenac and concealment is unlikely, suggests that a large change in the availability of diclofenac to livestock owners through pharmacy retail outlets occurred between 2009 and 2012. However, the lack of data from undercover surveys during 2009–2011 precludes a precise assessment of when during this period the change in the availability of diclofenac in veterinary pharmacies occurred. However, it is probable that the vulture population recovery began at some time during this period of uncertainty. The pharmacy survey results do not preclude the possibility that some veterinary use of diclofenac has continued in Nepal. Livestock owners may have found other ways to obtain the drug than through local pharmacy shops, such as sourcing it covertly from neighbouring areas of India. Some other veterinary drugs known or suspected to be toxic to Gyps vultures are used in Nepal (BCN/RSPB unpubl. data).
Most sightings of both vulture species were made close to Vulture Safe Feeding Sites: within 50 km of the nearest site for Slender-billed Vulture and within 100 km for White-rumped Vulture. This apparent association is partly because Vulture Safe Feeding Sites were mostly placed in parts of Nepal where vulture breeding sites were known to remain before the feeding sites were established. Rigorous testing of the effects of Vulture Safe Feeding Sites on population growth rates would require a comparison of population changes between regions with and without feeding sites, but with similar other conservation measures in place. Such a comparison is not possible using the available data. Gilbert et al. (Reference Gilbert, Watson, Ahmed, Asim and Johnson2007) found that satellite-tagged adult White-rumped Vultures in Pakistan often visited sites where diclofenac-free supplementary food was provided regularly and that the mean daily distance travelled during periods was reduced when food was provided. The mean number of vultures found dead per day in the Pakistan study area was also reduced during feeding periods, but it was still high enough to cause continuing rapid population decline. The Vulture Safe Feeding Sites in Nepal may have contributed to the observed increase in the numbers of vultures counted on transects by reducing the proportion of meals eaten by vultures which contained a lethal concentration of diclofenac. The cattle whose carcasses were provided at Vulture Safe Feeding Sites also did not receive other NSAIDs which are known (ketoprofen; Naidoo et al. Reference Naidoo, Wolter, Cromarty, Diekman, Duncan, Meharg, Taggart, Venter and Cuthbert2010) or suspected (nimesulide; Cuthbert et al. Reference Cuthbert, Taggart, Mohini, Sharma, Das, Kulkarni, Deori, Ranade, Shringarpure, Galligan and Green2016) to be nephrotoxic to vultures, but whose veterinary use has not yet been banned in Nepal. The programme took old and ailing cattle from surrounding towns and villages, many of which might otherwise have been treated with one of these NSAIDs before death. The Vulture Safe Feeding Sites diverted cattle which would die elsewhere in the region and therefore probably did not increase the total quantity of food available. However, provisioning at the sites was regular and predictable, so it is possible that some vultures from India may have been attracted to these sites in Nepal. Hence, immigration of vultures from India is a possible contributory factor to the rapid rate of population growth since 2012.
Our findings suggest that well-resourced, community-based advocacy and awareness-raising actions, combined with the provision of safe food, may be important in augmenting the effects of changes in the regulation of veterinary drugs toxic to vultures. Despite these encouraging results, a continuing lack of safety testing and regulation of other veterinary NSAIDs incurs the risk that vulture populations will again suffer rapid declines, but this time from a much lower starting point than was the case for the diclofenac-induced decline. The veterinary use of ketoprofen, which is known to be toxic to Gyps vultures: (Naidoo et al. Reference Naidoo, Wolter, Cromarty, Diekman, Duncan, Meharg, Taggart, Venter and Cuthbert2010), and of nimesulide, which is suspected to be toxic (Cuthbert et al. Reference Cuthbert, Taggart, Mohini, Sharma, Das, Kulkarni, Deori, Ranade, Shringarpure, Galligan and Green2016), have not yet been banned in Nepal.
Acknowledgements
This research was supported by Department of National Parks and Wildlife Conservation-Nepal and funded by the Government of the United Kingdom’s Darwin Initiative (18-008) the RSPB, and the Conservation Leadership Programme (03154313). We are grateful to Bird Conservation Nepal for logistical support and assistance that enabled these surveys to be undertaken. We thank Tulsi Ram Subedi, Jeet Bahadur Giri, Basu Bidari, Hem Subedi and Badri Chaudhary for their assistance during the fieldwork. We thank Ian Newton and Clare Stringer for constructive criticism of early drafts.
