Long-term population dynamics in a multi-species assemblage of large herbivores in East Africa

Wildlife population declines in Africa are widespread. However, species-specific population trends and dynamics in mammal community composition have rarely been described over long time periods. To describe population trends of 13 large herbivore species in Lake Manyara National Park (Tanzania) from 1959 to 2016 and to discover whether the herbivore community structure changed, we used general additive models and additional statistical methods to detect structural changes in the time series. Population dynamics were non-linear and population growth rates were not correlated with precipitation anomalies. Relatively steep population declines of three megaherbivores occurred during the 1980s and early 1990s, resulting in severe reductions in African elephant and buffalo populations and the local extinction of black rhinoceros. These declines coincided with reported peaks of illegal hunting of these species and expansion of agriculture at the periphery of the park. Population densities of elephant and buffalo seem to have stabilized in recent times, yet have not recovered to previous densities. In contrast, eight species (giraffe, zebra, waterbuck, wildebeest, warthog, impala, bushbuck, and baboon) have apparently fared well (similar or higher densities in most recent compared to first decade), despite having undergone substantial fluctuations over the past 58 yr. Population fluctuations in these species were likely caused by disease outbreaks, heavy bush encroachment, and reduced competition with buffalo. Possibly, declines in megaherbivore densities (mainly elephants) facilitated bush encroachment. Albeit grazers are still dominating in the herbivore community, the proportion of browsers is currently increasing, likely encouraged by dense vegetation in the shrub layer in large parts of the park. Overall, herbivore biomass density has declined by ~40% compared to the baseline estimate in the first decade of the time series. Our analyses and ancillary information provide evidence that this overall decline in the herbivore assemblage was triggered by human-induced reductions in megaherbivore population densities during the 1980s, either through excessive poaching, insularization of the park, or both. Likely, this had cascading and interacting effects on the vegetation structure and the herbivore assemblage. Thus, legacy effects of ineffective megaherbivore conservation efforts 30 yr ago are likely still affecting the ecology of this national park.

trends and dynamics in mammal community composition have rarely been described over long time periods. To describe population trends of 13 large herbivore species in Lake Manyara National Park (Tanzania) from 1959 to 2016 and to discover whether the herbivore community structure changed, we used general additive models and additional statistical methods to detect structural changes in the time series. Population dynamics were non-linear and population growth rates were not correlated with precipitation anomalies. Relatively steep population declines of three megaherbivores occurred during the 1980s and early 1990s, resulting in severe reductions in African elephant and buffalo populations and the local extinction of black rhinoceros. These declines coincided with reported peaks of illegal hunting of these species and expansion of agriculture at the periphery of the park. Population densities of elephant and buffalo seem to have stabilized in recent times, yet have not recovered to previous densities. In contrast, eight species (giraffe, zebra, waterbuck, wildebeest, warthog, impala, bushbuck, and baboon) have apparently fared well (similar or higher densities in most recent compared to first decade), despite having undergone substantial fluctuations over the past 58 yr. Population fluctuations in these species were likely caused by disease outbreaks, heavy bush encroachment, and reduced competition with buffalo. Possibly, declines in megaherbivore densities (mainly elephants) facilitated bush encroachment. Albeit grazers are still dominating in the herbivore community, the proportion of browsers is currently increasing, likely encouraged by dense vegetation in the shrub layer in large parts of the park. Overall, herbivore biomass density has declined by~40% compared to the baseline estimate in the first decade of the time series. Our analyses and ancillary information provide evidence that this overall decline in the herbivore assemblage was triggered by human-induced reductions in megaherbivore population densities during the 1980s, either through excessive poaching, insularization of the park, or both. Likely, this had cascading and interacting effects on the vegetation structure and the herbivore assemblage. Thus, legacy effects of ineffective megaherbivore conservation efforts 30 yr ago are likely still affecting the ecology of this national park.

INTRODUCTION
In mature ecosystems, key parameters such as primary productivity or biomass density of consumers are believed to fluctuate around a dynamic equilibrium (Scheffer and Carpenter 2003). However, a large number of aquatic, marine, and terrestrial ecosystems have steadily declined to lower levels during the past decades (Vitousek et al. 1997, Dobson et al. 2011 or have undergone catastrophic changes (Pauly 1995, Scheffer and Carpenter 2003, Folke et al. 2004, Andersen et al. 2009). For example, in East Africa, an index of overall large mammal abundance has declined gradually from 1970 to 2005 to about half of the baseline (Craigie et al. 2010). Assessments in Kenya and Tanzania have further provided evidence that wildlife populations are declining in strictly protected and multiple-use areas although wildlife populations usually fare better in strictly protected areas (Ottichilo et al. 2000, Stoner et al. 2007, Western et al. 2009). These large-scale assessments of wildlife population trends are important for informing policy but need to be supplemented with more detailed descriptions and analyses of species-specific population trends within particular areas in order to advance hypotheses on the controlling factors of wildlife population sizes (Borner et al. 1987, Kiffner et al. 2009, Ogutu et al. 2012. In essence, this approach follows the declining population paradigm which attempts to identify causes and mechanisms that lead to declines in animal populations (Caughley 1994). Such insights can then be used to adjust management strategies and to advance our overall understanding on how wildlife populations are affected by changes in their environment (Caughley 1994, Caro et al. 2013. Environmental or anthropogenic changes may cause gradual or even abrupt shifts in an ecosystem, which may manifest themselves in altered animal community structure (Vitousek et al. 1997, Scheffer andCarpenter 2003). Such structural changes may have strong cascading effects (Pringle et al. 2007, Hollings et al. 2014, Ripple et al. 2015a. For example, declines in browsing species (e.g., due to disease outbreaks or poaching) may release woody vegetation regeneration from strong herbivory and may thus facilitate bush encroachment. In turn, higher shrub cover possibly supports increased population densities of browsing species (Prins and van der Jeugd 1993). Such environmental dynamics are usually difficult to identify due to the lack of long-term data for entire species assemblages and vegetation structure, as well as finding appropriate analytical methods to best explore the dynamics of observed ecological dynamics (Andersen et al. 2009).
Unfortunately, reliable long-term and direct assessments of wildlife populations are generally scarce (Caro 2011). In Africa, most wildlife population trends are based on aerial surveys which may yield biased results, particularly for smaller species and in areas with dense vegetation (Jachmann 2002). An example for ground-based, longterm wildlife population monitoring comes from Lake Manyara National Park, Tanzania. Intermittent ground counts of the major herbivore species date back to 1959, the year before the area became a national park in 1960, and indicated that the park supports one of the highest recorded large mammal biomass densities in the world (Watson and Turner 1965, Coe et al. 1976, Prins and Douglas-Hamilton 1990. While several mammal species have changed in population size over a 26-yr period, the overall large herbivore biomass had remained stable from 1959 to 1984 (Prins and Douglas-Hamilton 1990). However, beyond variation in precipitation-a main environmental driver of ungulate population dynamics in African savanna ecosystems (Coe et al. 1976)-many changes have occurred both inside and outside the park in the last decades. For example, poaching for ivory had severely reduced the African elephant (Loxodonta africana) population and the black rhinoceros (Diceros bicornis) population had been locally extirpated by poachers in the late 1980s (Borner 1981, Newmark 2008). In addition, epidemics of rinderpest and anthrax have decimated populations of African buffalo (Syncerus caffer) and impala (Aepyceros melampus), respectively. Furthermore, the park is being encroached, particularly along its northern boundary by a fast-growing human population, leading to further insularization of the park (Mwalyosi 1991, Newmark 2008, Msoffe et al. 2011. In light of these manifold anthropogenic and environmental changes, we describe population trends of 13 large herbivore species over a time span of 58 yr  and ask the following question: How did herbivore populations and associated growth rates respond to potential causal factors such as precipitation, disease outbreaks, severe poaching, and changes in other species abundance? In addition, we analyzed how the herbivore community structure changed functionally, how bush cover in the park changed over time and discuss how vegetation, herbivore community, and other variables were interacting. Finally, we estimated the net effect of observed changes in species-specific population sizes and overall herbivore biomass to assess the overall conservation performance of the park (Craigie et al. 2010) and discuss the findings in the context of wildlife conservation efforts in East Africa.

Study area
Lake Manyara National Park (S 3.50°, E 35.77°) is located in the Tarangire-Manyara ecosystem of northern Tanzania (Fig. 1). Since its establishment, the park has undergone several boundary adjustments. Most notable are the extensions to the south and southeast. In addition, the former Marang Forest Reserve (~220 km 2 ) became part of the national park in 2009. Currently, a dry land area of 428 km 2 and about 220 km 2 of the alkaline lake are within the boundaries of the park. Located on the valley rim of the Eastern Rift Valley, elevation ranges from 960 to 2000 m above sea level (Prins and Weyerhaeuser 1987). The climate is tropical semi-arid to semi-humid with two distinct rainy periods (November-December and March-May), and a dry season in the remaining months (Prins and Loth 1988). From 1959 to 2016, annual precipitation ranged from 108 to 1203 mm (rainfall data collected by park management at the park headquarters). In addition to precipitation, the park receives freshwater influx from several seasonal and permanent streams that originate in the highlands (Prins 1996). The majority of the park has hard pan, plain soil, and clay, but there are also marshy areas in the park (Greenway andVesey-Fitzgerald 1969, Loth andPrins 1986a). Fertility of the soil and a high ground water table yield high vegetative productivity (Prins 1988(Prins , 1996. The main vegetation types are alkaline grassland, groundwater forest, riverine vegetation, Acacia woodlands, and escarpment bushlands (Loth and Prins 1986b).

Herbivore counts
We analyzed population trends of the 13 herbivore species that were initially analyzed by Prins and Douglas-Hamilton (1990): African elephant, hippopotamus (Hippopotamus amphibius), black rhinoceros, giraffe (Giraffa camelopardalis), African buffalo, Burchell's zebra (Equus burchellii), waterbuck (Kobus ellipsiprymnus), blue wildebeest (Connochaetes taurinus), warthog (Phacochoerus africanus), impala, southern (common) reedbuck (Redunca arundinum), bushbuck (Tragelaphus scriptus), and olive baboon (Papio anubis). Population sizes from 1959 to 1984 were obtained from Table 1 in Prins and Douglas-Hamilton (1990); population sizes were assessed by attempted total counts which were conducted by systematically driving through the park. During these earlier times , the small park size (~100 km 2 ) and rather open vegetation (and thus good visibility) provided suitable conditions for reliable total counts. Additional population size estimates derived from total ground counts conducted by the park management were obtained for the years 1985, 1992, 1996, 2007, 2008, and 2009. These counts followed the same protocol as the previous counts.
From 2011 to 2016, we conducted line transect surveys and applied distance sampling methodology to estimate herbivore density (Thomas et al. 2010). The first line transect survey was conducted in November 2011, and from 2012 to 2016, we conducted surveys three times per year (March/April, July, and October/November) to capture the main seasons. Because of the dense vegetation in most areas of the park (e.g., Table 1) and restrictions to off-road driving, line transects followed the road network in the park (Fig. 1). Although this design is not strictly systematic, the road network in this park covers all habitat types approximately in proportion to their availability given that the main road crosses the park from north to south. Line transects covered the core areas that were initially delineated as the national park. In addition, we covered areas to the south that had been annexed to the park in the 1990s, but not the Marang forest; the area sampled in the lowland area is~168 km 2 in extent. Line transects were usually 2 km in length and were separated by 500 m to allow for independence of transects. While driving along transects at slow speed (10-20 km/h), three to seven observers looked out for mammals. Upon detection of a mammal or mammal group, the vehicle was stopped and the animals were identified and counted, and the perpendicular distance between the transect and the animal (group) was measured with a laser rangefinder (Bushnell Elite 1500). During most surveys, we used multiple vehicles to cover the entire park at the same time in order to reduce potential bias due to animal movement. During some surveys, transects were repeated on consecutive days. In these cases, we combined the effort and sightings for each transect (Thomas et al. 2010). With slight modifications due to road constructions or inaccessibility, we used the same transects (mean = 33) during each of the 16 surveys. Estimating animal density Because the park boundaries (and hence the covered counting area) changed multiple times since the establishment of the national park (from 100 to 168 km 2 of lowland areas), we divided population sizes that were derived from total counts ) by the lowland surface area of the national park in the matching year and used animal density (individuals/km 2 ) for all consecutive analyses. Animal densities for 2011-2016 were estimated using line distance methodology, implemented in the software Distance 6.0 (Thomas et al. 2010). Line distance sampling takes a variation in detectability of target species into account and is thus a robust analysis technique to estimate animal density. For each species, we pooled observations of the 16 surveys (since 60-80 detections are required to model robust detection functions) and fitted six different detection functions to the data (Thomas et al. 2010). We considered four conventional detection models (all with cosine expansion) that describe detectability as a function of perpendicular distance: uniform, half-normal, negative exponential, and hazard rate. In addition, we considered that detectability may vary by season and fitted two additional detection models (half-normal and hazard rate) that describe detectability as a function of perpendicular distance and the threelevel covariate season. Based on the Akaike's information criterion (AIC) score of each model and the visual fit of the detection function near the transect line and goodness-of-fit statistics, we selected the best-fitting detection model for each species and used this model to estimate seasonand species-specific densities (Thomas et al. 2010). For 2012-2016, we averaged the means of the three-seasonal density estimates to generate yearly mean densities for each species.

Validating density estimates
Since the counting methods changed after 2009, we employed independent techniques to calibrate the data derived from line transect surveys for five out of ten species.
For wildebeest and zebra, we conducted an independent total count simultaneously with the line transect survey in November 2016. Observers in one vehicle started at the southernmost road in the park and counted all zebra and wildebeest. The vehicle systematically followed the road northward to the gate, completing each of the park's road circuits on the way to cover as much area as possible and to avoid double counting of moving individuals. Because both target species prefer open habitats and since these areas are largely observable from the road network, this method likely produces an accurate estimate. Total counts were repeated on two consecutive days, and we used the highest daily estimate as a reference point for the distance sampling estimate. To estimate hippopotamus densities, the park management conducted a total count in 2014. We assumed no population size change from 2011 to 2016 in this species and used the 2014 figure (divided by the year's park size) as population density estimates for the years 2011-2016.
For estimating buffalo and elephant abundance, we used a mark-re-sight approach (e.g., Morrison et al. 2016b). On day one, observers in two vehicles searched for buffalo and elephant while systematically driving along all roads. We recorded individual and herd characteristics of target individuals and, if possible, obtained digital pictures of individuals. On day two, the same observers traveled in the same manner as on day one. Upon encountering buffalo or elephant and using recorded descriptions and digital photographs from the previous day, observers assessed whether each Notes: Values from 1985 and 1991 are from Prins and van der Jeugd (1993), and for 2017, we included 95% confidence intervals. The depiction of vegetation communities followed Prins and van der Jeugd (1993) albeit scientific nomenclature has changed for some species. Croton macrostachyus-Phoenix reclinata communities described in Prins and van der Jeugd (1993) were not found in 2017.
❖ www.esajournals.org individual had been captured on the previous day. Individual identification of elephants was unambiguous based on photographs taken on the previous day. For buffalo, it was not possible to obtain pictures of each individual, particularly in large breeding herds. Thus, for individuals in large buffalo groups, we used herd location, demographic structure of the group, and group size as indicators of whether individuals had been seen on the previous day in addition to individual characteristics. Therefore, the buffalo estimate may be associated with some level of uncertainty. We created a capture history for each individual and estimated abundance using a closed population mark-recapture model implemented in the R package RMark and selected the model with the lowest AIC score (R Core Team 2012, Laake 2016). For giraffes, we used density estimates from a recently conducted mark-re-sight survey in which giraffes were individually identified using pattern-matching software (Lee 2015).
We calculated a correction factor (total count or mark-re-sight estimate divided by distance sampling estimate) for wildebeest (correction factor: 0.7560), zebra (0.7466), giraffe (0.6301), buffalo (1.6560), and elephant (0.1605), and used this factor to calibrate the yearly density estimates. Despite this calibration effort, estimated elephant densities in 2014, 2015, and 2016 were considered unrealistically high. This is probably due to a preference for elephants to use the roads, low horizontal visibility, and associated low-effective strip width of transects, which may cause high elephant density estimates based on distance sampling procedures in this park. We hence used the markre-sight abundance estimate of 2016 (266 elephants), the population size assessed by a total count by the park management in 2014 (250 elephants), and the average of the two for 2015 (258 elephants).
For other species, we used non-transformed line transect estimates since alternative methods appeared rather impractical and unreliable in dense habitats and for species that are very abundant (baboon, bushbuck, impala, southern reedbuck, warthog, waterbuck). Density estimates for all survey years can be found in Appendix S1.

Estimating herbivore consumption and herbivore biomass
To analyze whether the herbivore assemblage changed functionally, we calculated the amount of grazing, browsing, and total herbivore biomass consumption per unit area and year, as well as the ratio of grazing to browsing. We calculated dry matter (DM) consumption using allometric scaling theory [DM (kg/day 2 ) = a 9 Bodymass À0.75 ] with a = 46 g/kg 2 (M€ uller et al. 2013). For elephants (a = 127.5 g/kg 2 ) and black rhinoceros (a = 33.5 g/kg 2 ), we used species-specific scaling parameters because these species appear to differ in their feeding physiology from other herbivores (Clauss et al. 2007). We used body masses (Prins and Douglas-Hamilton 1990) and proportions of grass and browse in each species' diet [bovids (Gagnon and Chew 2000); rhinoceros (Malan et al. 2012); other species: (Prins and Douglas-Hamilton 1990)] and multiplied these with the corresponding species' densities to estimate grass, browse, and overall (grass, browse and all other food items) consumption per hectare (ha) and year.
We estimated overall herbivore biomass density (kg/km 2 ) by summing up the product of body mass and the corresponding annual density estimate to summarize all species' densities to a single number.

Bush encroachment
To assess how the savanna vegetation changed over time, we estimated horizontal visibility-a measure that is closely and positively correlated with shrub cover (Prins and van der Jeugd 1993) -in four distinct vegetation communities in February 2017. Using a laser range finder, we measured the distance that observers could see horizontally without being obstructed by the vegetation. In line with Prins and van der Jeugd (1993), we measured horizontal visibility ten times along 48 500-m transects in each vegetation community in the area of the 1990 park boundaries (a total of 480 visibility scans). We compared the 2017 values with baseline horizontal visibility data that were assessed in a similar manner in 1985 and 1991 (Prins and van der Jeugd 1993). Since baseline data (1985 and 1991) were only available as mean values for each vegetation stratum, we calculated 95% confidence intervals for the 2017 measurements and assessed whether the stratum-specific mean values of 1985 and 1991 differed from the corresponding 2017 estimates. We did not make any vegetation assessments in the alkaline grasslands because we did not have baseline data and ❖ www.esajournals.org because we did not assume major changes in this particular vegetation stratum.

Time series of wildlife populations
Visual analysis of the time series suggested that linear analyses are inappropriate and we hence used statistical methods that allow for detecting and modeling non-linear changes (Andersen et al. 2009). To describe the time series of the sporadic census data, we used general additive models implemented in the package mgcv within the statistical software R (R Core Team 2012, Wood 2016). The variable year was modeled as a cubic spline. The degrees of freedom of the cubic spline smoother were chosen based on the AIC score of models with variable degrees of freedom of the smoothing term. To account for relatively small sample size (n = 21 yr of census data), we used a penalty of k = 4 for calculating the AIC score. To detect structural changes in the time series, we calculated breakpoints using the strucchange package in R (Zeileis et al. 2002). These breakpoints indicate the year and associated confidence intervals when the data were structurally different in a linear regression relationship and thus indicate when abrupt changes in the time series occurred (Zeileis et al. 2002).

Population growth rate and causes of population growth
We used the selected non-linear models to predict yearly population densities for each species, grass and browse consumption, grazing/browsing ratio, and herbivore biomass. Based on this complete, yet interpolated time series, we estimated the finite population growth rate k = D t+1 / D t , where D t and D t+1 are the densities of the preceding and following year, respectively. We tested whether species-specific population growth rates were correlated with deviations from the average rainfall (mean 1958-2015: 609 mm) of the same and the preceding year (since population growth may respond in a timelagged fashion) using non-parametric Kendall's correlation coefficient. Rainfall data were consistently collected at the park headquarters from 1953 to 2015. Since other potential causal factors of population growth (poaching, disease outbreaks, human development outside the park) were non-continuous and of unknown duration, we were unable to employ a strong inference approach (e.g., a multivariate model). Instead, we plotted these events-gathered from published articles (Prins and Weyerhaeuser 1987, Prins and van der Jeugd 1993, Prins 1996, Kiffner et al. 2015b) and from the park management-on the timeline of population growth rates and identified whether the timing of these events coincided with the identified breakpoints. This analysis was restricted to species with obvious breakpoints and well-documented historic events that may have caused population fluctuations.
Since we observed a drastic decline in buffalo densities and given previous research suggesting that competition with buffalo may limit populations of other grazing species such as zebra, wildebeest, waterbuck, and warthog (Sinclair 1985, Ogutu et al. 2012, we tested whether population densities of grazing herbivores were correlated with buffalo densities using Kendall's correlation test.

Net changes: comparing the first and last decade
To assess the net effect of changes in the herbivore community, we used the general additive models and predicted the mean population density of each species and overall herbivore biomass of the first (1959)(1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968) and the last decade (2007-2016) and compared mean predicted densities of the two decades with each other. While this analysis blurs the non-linear population dynamics of most species, it provides an overview of net changes in the animal community and allows for an assessment of the relative conservation performance of the park.

Bush encroachment
Compared to 1985 and 1991, mean horizontal visibility in 2017 was substantially reduced in all four considered savanna vegetation communities. Using 1985 data as baseline, horizontal visibility had been reduced to 43-95% in 1991 and to 9-38% in 2017, suggesting a considerable increase in shrub cover over time (Table 1).

Population trends and associated drivers
Population densities of the 13 analyzed species showed strong fluctuations from 1959 to 2016 (Fig. 2). Population growth rates (k) of all species except wildebeest were not significantly correlated ❖ www.esajournals.org with deviations from the mean precipitation (P ≥ 0.17). Population growth of wildebeest was weakly and positively correlated with rainfall anomaly (s = 0.189, P = 0.038). Rainfall anomalies in the previous year were not significantly correlated with population growth rates of any species (P ≥ 0.07; Appendix S2). Overall, precipitation did not change unidirectional over the assessed time span (F 1,55 = 0.9947; P = 0.323). Taken together, these findings suggest that population dynamics were not sufficiently explained by variation in precipitation (except partially for wildebeest) and that other causes were likely to be more influential in governing population dynamics of the herbivore assemblage (Appendix S3).
The elephant population density remained fairly constant in the first decades, declined during the 1980s, and seemed to have remained fairly stable since 2000 (Fig. 2). The main turning point was 1984, which coincides with reported high levels of poaching for this species (Appendix S3). Yet, elephant growth rates appear to have dropped below 1 earlier than 1984 (Fig. 3), coinciding with a Mycoplasma pneumonia outbreak in 1977 and substantial human development to the north of the park (Appendix S3). The hippopotamus density appears to have slightly declined over time (mean k = 0.992) with indicated turning points in 1975, 1984. Black rhinoceros were locally extinct by the late 1980s or early 1990s (Fig. 2), whereas the main breakpoints were 1965, 1975, and 1984. The population density declined particularly strongly in the early 1980s, coinciding with severe poaching (Fig. 3). Giraffe population density steadily increased over most of the time series, whereas a declining trend in more recent years becomes evident (Fig. 2) but no structural change in the time series was detected (Appendix S3). African buffalo density was high initially, and the population recovered apparently quickly from the rinderpest outbreak in 1959 (Appendix S3). The population density declined sharply during the early 1980s, coinciding with local development (expansion of settlements, conversion of swamps to irrigated agriculture just outside the national park) and high poaching pressure. In the last few years, the population density seems to have stabilized and shows signs of population increase (Fig. 3). Population trajectories of zebra density show an almost linear, positive trend (mean k = 1.010; Fig. 2), whereas the population growth seems to have slowed down in the mid-1980s (Fig. 3), coinciding with agricultural development outside the national park (Appendix S3). The wildebeest population density showed substantial fluctuations over time; after the Rinderpest outbreak in 1959, population growth was negative, around 1970 it was stable, and it was positive in the 1980s. From 1992 to 2010, population growth was negative, whereas the wildebeest population growth has been positive since 2010 (Fig. 3). The waterbuck population density followed a humpshaped curve with an increasing phase from 1984 to 2009 and a decreasing slope from 2009 to 2013 (Fig. 2); the (alkaline) lake dried up several times during both periods of waterbuck population density increase and decrease (Appendix S3; fresh water flow data from rivers and streams are not available). The warthog population density showed asymptotic growth. Interestingly, population densities of all grazers were significantly and negatively correlated with buffalo densities (s zebra = À0.664; s wildebeest = À0.312, s waterbuck = À0.443, s warthog = À0.505; all P-values < 0.001).
Bushbuck densities fluctuated over time and seem to have increased substantially in the last years of the time series. Similar to bushbuck, the impala population density fluctuated over time; from 1981 to 2009, it showed a slightly negative trend (possibly triggered by earlier anthrax outbreaks; Appendix S3), whereas it increased from 2009 onward (Fig 2); in the past few years, impala population densities were higher than from 1959 to 2000. The reedbuck population density showed a declining trend initiated already in the late 1960s; since 1976, reedbuck occur at very low population density (Fig. 2). Baboon population densities remained at relatively constant levels and then showed an exponential increase, starting around 2008 (Fig. 2).

Herbivore community structure
Grass consumption, and thus density of grazing animals, of the herbivore assemblage was at high and relatively constant levels from 1959 to 1984, then dropped, and has stabilized at a lower level since 2008 ( Fig. 4; Appendix S4). After a long period of relative constancy, the consumption of browse declined gradually beginning in the late 1970s and then more steeply in the mid-1980s (coinciding with severe poaching, local ❖ www.esajournals.org development, and disease outbreaks); in the last few years of the time series, the amount browsing appears to have stabilized and shows an increasing trend ( Fig. 4; Appendix S4).
Overall, herbivore consumption dropped substantially since 1984 (coinciding with severe poaching, local development, and disease outbreaks) and appears to have reached a new, yet much lower level ( Fig. 4; Appendix S4). The grazing/browsing ratio increased in the early 1980s and then decreased since the late 1990s (Appendix S4). In recent years (around 2009-2013), the herbivore community structure has shifted more toward browsing species, yet grazers are still dominating the overall herbivore community.

Net effects of population changes
Overall, four species (elephant, hippopotamus, African buffalo, and southern reedbuck) had substantially lower densities in the last decade compared to the first decade of data collection, and one species-black rhinoceros-is locally extinct (Table 2). Compared to the first decade of the time series, densities of eight species (giraffe, Burchell's zebra, blue wildebeest, waterbuck, impala, bushbuck, warthog, and olive baboon) increased substantially (Table 2). Despite population density increases of the majority of the species (8/13), the biomass density of the herbivore community has declined substantially and is recently only 58% compared to the first decade of the time series (Table 2). Biomass density declined sharply in the mid-1980s and appears to have increased slightly in more recent years (Fig. 5), whereas the main biomass decline (which occurred from the late 1970s to the late 1990s) can be mainly attributed to reduced elephant and buffalo densities (Fig. 6).

DISCUSSION
Lake Manyara National Park in Tanzania is near unique because of the very long record of direct and ground-based animal censuses. Population trends derived from intermittent surveys of large herbivores over a 58-yr period showed that the decline of elephant and buffalo population densities and the local extinction of black rhinoceros during the 1980s led to an overall 42% reduction in herbivore biomass density and a proportional change toward browsers in the herbivore assemblage. These results are in stark contrast with wildlife population trends from 1959 to 1984, which suggested overall stability in Fig. 4. Estimated consumption (kgÁha À1 Áyear À1 ) of vegetation by 13 large mammal species in Lake Manyara National Park over time, separated by grazing, browsing, and total consumption combined as well as the ratio of grazing to browsing over time. The temporal trends and associated confidence limits were estimated with a cubic spline in a general additive model environment.
this herbivore assemblage (Prins and Douglas-Hamilton 1990) and thus highlight the importance of very long-term monitoring and comprehensive data in order to understand ecological dynamics (Peters 2010). The observed patterns and timing of population dynamics in combination with ancillary information may allow for the formulation of hypotheses about the causality underlying the observed population trends (Caughley 1994, Venter et al. 2014, Bragina et al. 2015, Shoemaker and Akc ßakaya 2015. We are very much aware of the disparity of the data and the data collection methods that were used over the last 60-odd years. We thus realize that the data may not be sufficient to test causal ideas about population dynamics and herbivorevegetation interactions, but our long data series is at the same time perhaps one of the best available long-term datasets of terrestrial large mammal populations.

Population trends and associated conjectured factors
The massive elephant decline during the 1980s coincided with the estimated 60% poachingrelated mortality in the elephant population (Prins et al. 1994). Indeed, elephant poaching peaked during this time across East Africa (Caughley et al. 1990). An alternative proposition-conversion of swamps at the north of the park to agriculture as a main cause of elephant decline-is thus not very plausible but may perhaps have acted as a synergistic factor (Brook et al. 2008). Recent research suggests that the park's elephant population presently has a high proportion of young and sub-adult elephants, as Table 2. Predicted mean densities and associated 95% confidence limits of herbivore species (individuals/km 2 ) and overall herbivore biomass density (kg/km 2 ) in Lake Manyara National Park 1959-1968. Historic: 1959-1968Recent: 2007   well as a high female-to-male ratio, indicative of an increasing population. However, there is still evidence that occasional poaching of elephants occurs, mostly taking place outside park boundaries (Kioko et al. 2013). For example, between 1997 and 2002, 43 cases of elephant poaching were reported with only one taking place inside the park (Kaswamila 2009), and in 2011-2012, nearby Manyara Ranch experienced substantial elephant poaching (Kioko et al. 2013). In contrast to other protected areas in Africa where massive declines occurred since~2009, the elephant population in Lake Manyara National Park seems at present to be rather well protected from the more recent and continent-wide poaching crisis (Chase et al. 2016). In terms of absolute population size, the current (2014-2016) elephant population (~260 individuals) is less numerous than from 1959 to 1984 (~400 individuals).
Hippopotamus have seemingly declined linearly over time but have not undergone drastic population crashes. Since the hippopotamus population did not spread substantially from its original location in the north of the park, the observed decline in density (from about 1.5 ind/ km 2 to about 1.2 ind/km 2 ) is likely the result of the park doubling in extent and actually reflects an increase in overall numbers (from about 150 to about 200). However, park extension is not a valid explanation of density declines in other terrestrial megaherbivores since the areas annexed to the park contained suitable habitat, those areas receive higher amounts of rainfall and are thus likely to even support higher herbivore populations. In addition, if park boundary had substantial effects on herbivore population sizes, we would have expected to detect structural changes in the years of (or directly following) park extensions. This was not the case.
In contrast, the local extinction of black rhinoceros can unambiguously be ascribed to excessive poaching to feed the demand for rhinoceros horn during the 1980s, which eradicated not only the rhinoceros population in this park but indeed led to the extinction of this species in several protected areas of Africa (Borner 1981, Leader-Williams et al. 1990, Newmark 2008). Among the megaherbivores, giraffes have maintained relative constant population sizes since 1959. However, the more recent negative population growth has also been linked to mortality in adult females, likely caused by unsustainable and illegal hunting that is mainly taking place outside the park (Kiffner et al. 2015a, Lee 2015, Lee et al. 2016. Population declines of buffalo during the 1980s and 1990s were very apparent and contributed Fig. 6. Estimated biomass density (kg/km 2 ) for 13 large herbivore species in Lake Manyara National Park for years from 1959 to 2016. markedly to the decline of grass consumption and overall herbivore biomass (Fig. 6). The park once held exceptionally high buffalo densities, and buffalo quickly recovered after a severe rinderpest outbreak in 1959 Weyerhaeuser 1987, Prins 1996). According to Prins (1996: 11), the apparent reduction in water influx from the escarpment and the resulting temporal desiccation of streams during the 1980s can be excluded as a limiting factor for the resident buffalo population as surface water was (and still is) readily available. Prins (1996: 4) reported that hunting for bushmeat at the periphery of the park since 1991 decreased the buffalo population to 700 individuals in 1994 compared to c. 2300 individuals in the 1960s and 1980s. This timing of documented and excessive hunting falls within the confidence limit of structural changes in the time series (Appendix S3) and thus provides strong evidence that buffalo decline during the 1980s and 1990s was primarily driven by illegal hunting. In the same period, rinderpest had been eradicated and no additional outbreaks of this once decimating disease were reported. The recent annual population growth of buffalo can possibly be attributed to sustained anti-poaching efforts and a series of lake drying events in the recent decade that also coincided with a population breakpoint (Table 1). Lake drying allows for the expansion of more available grassland, which could increase food availability for buffalo that now share grasslands with much higher densities of zebra and wildebeest (Table 2). Thus, it seems that buffalo populations can cope relatively well with (moderate) droughts, as also exemplified by the buffalo population in Lake Nakuru National Park (Kenya), which increased significantly from the 1970s to 2012 despite frequent droughts (Ogutu et al. 2012). The anthrax outbreaks in 1961The anthrax outbreaks in , 1977The anthrax outbreaks in , and 1984 likely contributed to the 1985 breakpoint in the zebra time series and apparently slowed down the population growth (Prins and Weyerhaeuser 1987). The outbreak of rinderpest in 1959 also likely explains the initial small population size of wildebeest. However, release from this disease through extensive cattle vaccinations likely led to the density increases since the 1980s, and coincides with other cases of wildebeest population recovery after being released from this virus . The more recent increase in wildebeest and zebra populations is possibly linked to decreased competition over grass, coinciding with the substantial decline in buffalo (Sinclair 1974, Ogutu et al. 2012. Indeed, population densities of zebra and wildebeest were negatively correlated with buffalo densities, lending empirical support for the hypothesis that populations of these grazers are partially controlled by competition with buffalo as has been observed elsewhere in East Africa , Prins 2000, Ogutu et al. 2012). An additional, mutually non-exclusive explanation may be immigration of these two species to within the park boundaries to avoid hunting and human encroachment outside park boundaries (Msoffe et al. 2011, Kiffner et al. 2014. Since buffalo declines occurred in this same period of human development, it is unclear which of these factors (release from competition or immigration or both) led to the increased density of these two grazing species. The brief disappearance of the park's wildebeest population in the late 1960s and early 1970s is most likely attributed to a temporary migration out of the park due to flooding of the grasslands in that period (Prins 1987). However, the wildebeest population in the park is considered to be rather resident in recent times (Morrison and Bolger 2012). Compared to the migratory wildebeest population in nearby Tarangire National Park which has experienced a substantial population decline compared to historic times, the resident wildebeest population in Lake Manyara National Park showed higher recruitment in recent years (Morrison et al. 2016a).
In light of the quantitative evidence for competition with buffalo (see also Mwasi et al. 2013), it seems plausible that the relative high density of waterbuck during the 1990s and 2000s and the higher levels of warthog numbers since the 1980s were facilitated by reduced competition for grass. However, additional research is required to understand why the waterbuck population has declined in more recent years, with bush encroachment being a plausible hypothesis explaining this decrease. Indeed, bush encroachment particularly set in after the collapse of the impala populations which preceded the strong decline in elephant numbers (Prins and van der Jeugd 1993). The decrease in horizontal visibility suggests that bush cover in the park has substantially increased in the past 30 yr throughout all major vegetation types (Table 1).
Given the strong preferences of the bushbuck for dense habitats and browse (MacLeod et al. 1996, Brock et al. 2003, bush encroachment is thus a promising candidate to explain the overall increase in bushbuck densities. Anthrax outbreaks certainly caused fluctuations in the impala population, killing 1400 individuals in a year during the 1961 outbreak and even more during the outbreak of 1985 (Prins and Weyerhaeuser 1987). However, in line with the hypothesis for strong density-dependent population growth in smaller mammals (Sinclair 2003), impala populations recovered quickly after disruptions and recently occurred at densities higher than ever historically sustained by the park (Fig. 2). While we cannot unambiguously assign causality between herbivore crashes and bush encroachment, reductions in impala and elephants numbers likely facilitated bush encroachment due to reduced browsing pressure (Prins and van der Jeugd 1993). While many factors (e.g., soil, precipitation, herbivory, and fire) may facilitate bush encroachment (Ward 2005), there is substantial evidence that elephants play a key role in determining mortality of woody vegetation (Morrison et al. 2016a). Since elephants also exhibit a disproportionally high food intake (Clauss et al. 2007), it is very likely that a decline in the elephant density promoted bush encroachment. As mixed feeders, impala are exceptionally adaptable to changing environments and can utilize habitats that have become encroached by shrubs and bushes and this adaptability may have in turn facilitated the more recent high population density of impala.
The population of southern (common) reedbuck has historically been small (Appendix S1) and reductions in swamp vegetation caused by drying up of seasonal streams and an increase in alkalinity may have further reduced their population size (Loth and Prins 1986b). Indeed, from 2011 to 2016, only two reedbuck sightings have been recorded by the first author, and the fate of this small population is likely subject to stochastic events (Caughley 1994).
Bush encroachment may have facilitated the substantial increase in baboon numbers in recent years. Changes in methodology (total counts vs. line distance sampling) may partially be responsible for this increase in estimated baboon population densities, whereas density estimates before 2011 were likely underestimates. While we have no validation of absolute baboon density (total counts were attempted in 2016 but were considered not feasible in the dense vegetation), baboons are considered to be superabundant in the park (D. Zinner and S. Knauf, personal communication). Similar increases in baboon densities have been observed in other African protected areas and have been linked to reductions or extinctions of top predators (Prugh et al. 2009, Brashares et al. 2010). However, lion (Panthera leo) and leopard (Panthera pardus) populations still thrive in the park (C. Kiffner and J. Kioko, personal observation), and thus, the causation of the apparent steep baboon population increase remains unclear.

Dynamic of the animal community structure
Our analyses suggest that a rather steep decrease in the overall animal community occurred in the 1980s due to the major decline of dominating herbivore species (namely elephant, black rhinoceros, and African buffalo). Most of the evidence supports the notion that high levels of poaching for ivory and bushmeat initiated this drop in elephant and buffalo densities (Prins et al. 1994, Prins 1996. While local development outside the park also occurred during the same time and may have reduced animal movement between protected areas, it is unlikely that it leads to massive reductions in recruitment or increased mortality of elephant or buffalo populations. Instead, it is plausible that human population growth may have facilitated unsustainable illegal hunting. In combination with natural dieoffs (Young 1994) that led to short-term reductions in impala densities, this altered herbivore assemblage may have facilitated bush encroachment (Prins and van der Jeugd 1993), which in turn allowed a time-lagged recovery of browsing species. The proportional succession toward browsers mainly manifested itself around 2009 and thus several years after the first signs of bush encroachment in 1991 (Fig. 3). Indeed, such shifts in herbivore assemblages from grazers to browsers have been documented elsewhere in Africa following woody vegetation encroachment (Smit and Prins 2015). Ironically, in the 1970s, elephant impact on the vegetation was a major concern (Vesey-FitzGerald 1973), whereas some 30 yr later the opposite pattern seems to be the case. In addition, lowered competition with buffalo likely facilitated population growth of smaller grazing species.
More broadly, our findings are consistent with research showing that the likelihood of substantial ecological change increases with human disturbances (Folke et al. 2004, Craigie et al. 2010. Our time series supports notions that ecosystem dynamics are typically the result of non-linear dynamics and that ecosystems are naturally changing over long time periods (Andersen et al. 2009, Peters 2010. Moreover, we provide evidence for biotic homogenization, that is, that local anthropogenic interventions mainly affect large-bodied (elephant, black rhinoceros, buffalo) and specialized species (reedbuck) and favor smaller and usually generalist species such as impala and baboons (McKinney andLockwood 1999, Le Viol et al. 2012), which may be followed by a time-lagged loss of ecosystem services (Dobson et al. 2011, Ripple et al. 2015b).

Implications for management and conservation
Ideally, addressing the declining population paradigm (Caughley 1994) requires diagnosing causes of processes such as death rates and reproduction (Peery et al. 2004) rather than correlating population sizes with hypothesized drivers such as precipitation, competition, or anthropogenic threats (Shoemaker and Akc ßakaya 2015). Such demographic data are, however, rarely available for an entire animal assemblage. With due caution associated with correlation analyses, our time series provides support for the idea that humancaused top-down control of megaherbivores (which are usually not sensitive to natural predation due to their large body mass, except African buffalo in the present study) as well as natural die-offs can have cascading effects via changes in the vegetation structure and species assemblage in a terrestrial ecosystem and should thus be considered keystone species (Sinclair 2003, Ripple et al. 2015b). However, because of manifold interactions between these species (either competitive or facilitative, or by means of apparent competition), one cannot easily identify particular species with the keystone role. While smaller species recover faster from population crashes due to higher plasticity in growth rates, populations of larger-bodied mammals take long time periods to recover (Sinclair 2003). For practical conservation efforts, this means that conservation efforts need to be sustained for several decades until measurable success becomes apparent (Gaston et al. 2008, Lindsey et al. 2014, Coetzee 2017 and that direct and indirect interspecific and species-habitat interactions need to be considered more strongly when conserving or actively managing herbivore assemblages (Mwasi et al. 2013).
Essentially, our study shows that ineffective conservation some 30 yr ago-namely not effectively protecting black rhinoceros, African elephant, and African buffalo-reverberates for decades in an ecosystem. It is plausible that reductions in these megaherbivore populations set off a number of causal chains leading to severe shrub encroachment and changes in the herbivore assemblage. Even though illegal hunting has likely been much reduced compared to mid-1980 levels (poaching still occurs outside and inside the park), it is unlikely that the present ecosystem can easily revert back to the state that Lake Manyara National Park was in before the mid-1980s, and that present biomass density can double to revert to its original state. A perilous concoction has arisen because the declining population paradigm and the small population paradigm (Caughley 1994) now occur together since declining populations become small populations and thus may now be subject to stochasticitythe southern reedbuck as a case in point.
The declining population paradigm optimistically suggests that if the causes of the decline are removed, then the populations may revert to their original numbers (Caughley 1994). However, that may seem logical from a theoretical perspective without taking into account the reality of an ecosystem that now likely functions differently from the way it did when herbivore densities were nearly twice as high. Indeed, the conditions for the present-day wildlife assemblage are different than conditions were before the 1980s. The savanna habitats are much denser in woody undergrowth and have fewer grass cover, and animal movement in and out of the park (though still feasible) is substantially affected by human infrastructure and activities (Morrison and Bolger 2014). Fundamentally, the park faces an absence of the black rhinoceros, fewer habitat-engineering elephants, much more shrub and bush cover, and fewer buffalo to maintain the grasslands. In sum, these legacy effects likely militate against a full recovery of a system that once held the highest recorded mammalian biomass density on earth (Coe et al. 1976, Prins andDouglas-Hamilton 1990).