Spatial variation in anuran richness, diversity, and abundance across montane wetland habitat in Volcanoes National Park, Rwanda

Abstract The spatial distribution of species has long sparked interest among ecologists and biogeographers, increasingly so in studies of species responses to climate change. However, field studies on spatial patterns of distribution, useful to inform conservation actions at local scales, are still lacking for many regions, especially the tropics. We studied elevational trends and species‐area relationships among anurans in wetland habitats within Volcanoes National Park (VNP) in Rwanda, part of the biodiverse Albertine Rift region. In VNP, wetlands are key sites for anuran reproduction, and anurans are likely threatened by wetland desiccation which has occurred for the last few decades. Between 2012 and 2017, we sampled anuran communities in ten VNP wetlands located along an elevational gradient of c. 600 m (from 2,546 to 3,188 m a.s.l.) and found at least eight species, including at least two Albertine Rift Endemics. We show that species richness, diversity, and abundance likely decline with a decrease in wetland size and with an increase in elevation, though additional sampling (e.g., at night) might be needed to derive definite conclusions. Larger wetlands at lower elevations contained most species and individuals, which indicates the potential threat of wetland size reduction (through desiccation) for anuran conservation. However, we also found that wetlands differed in species composition and that some species (e.g., Sclerophrys kisoloensis) were likely restricted in distribution to only a few of the smaller wetlands—suggesting that the conservation of each individual wetland should be prioritized, regardless of size. We propose that all wetlands in VNP require additional conservation measures, which should be based on knowledge gathered through long‐term monitoring of anuran communities and research on drivers of wetland decline. Only such extended research will allow us to understand the response of anurans in VNP to threats such as climate change and wetland desiccation.


| INTRODUC TI ON
An understanding of the spatial distribution of biodiversity allows us to address how evolutionary and ecological processes shape ecosystems and guides conservation efforts, such as the spatial planning of protected areas. Furthermore, only with an understanding of current spatial patterns will we be able to assess species responses to environmental change (e.g., climate or land cover change) or predict future distributions of species under different scenarios of change.
Patterns of elevational variation in species richness and species-area relationships are two especially well-studied topics in biogeography, enabling us to increasingly understand how historic, evolutionary, climatic, and geometric factors are linked to the distribution of biota (Dengler, 2009;Laiolo, Pato, & Obeso, 2018;Lomolino, 2001).
In certain regions, and for some taxa, species richness peaks at low elevations or plateaus at mid elevations due to geometric constraint known as the mid-domain effect (Colwell, Rahbek, & Gotelli, 2004). Alternatively, richness may decline monotonically with an increase in elevation, a pattern that closely resembles latitudinal trends and is considered to be predominantly driven by environmental correlates of elevation, such as a decline in temperature (Rahbek, 1995). Which pattern dominates in a given ecosystem or locality depends on many factors, such as the scale of study (Rahbek, 2005) and local and regional climatic conditions (Fu et al., 2006). In fact, elevational patterns may even differ between southern and northern slopes of the same mountain (Shuai, Ren, Yan, Song, & Zeng, 2017).
Thus, despite the wealth of existent literature, we cannot extrapolate from one study to the next. It, therefore, remains imperative to extend studies of elevational species distribution to fragile or threatened ecosystems in understudied regions, such as Central and East Africa, where these studies may inform conservation at local scales.
Similar to the exploration of elevational patterns, we have advanced considerably in our understanding of the relationships between species diversity and area; relationships that were first given context by the theory of island biogeography (MacArthur & Wilson, 1967). Although a general tendency of increasing species richness with increases in area of investigation remains a common result of studies in biogeography and landscape ecology, there are now many alternative explanations for underlying mechanisms (for an overview of some of the most important, long-standing, hypotheses see McGuinness (1984) and Connor and McCoy (1979)). Given this complexity of underlying mechanisms and associated describing models, it is unsurprising that it is hard to generalize species-area relationships, even within one taxon, such as amphibians.
Anurans and other amphibians play many intrinsic roles in ecosystems (Hocking & Babbitt, 2014), and local extirpation of amphibians may lead to severe alterations of ecosystem structure and functioning (Whiles et al., 2006). Unfortunately, many amphibian populations are experiencing declines because of a host of interacting threats, including climate change, chemical contamination, disease, the spread of invasive species, and habitat loss and degradation (Blaustein et al., 2011;Kiesecker, Blaustein, & Belden, 2001).
Here, we explore the existence of elevational gradients and species-area relationships in species richness, diversity, and abundance among anurans (i.e., frogs and toads) in the Virunga Massif, which is part of the Albertine Rift biodiversity hotspot in Central Africa.
We know relatively little of anuran ecology in the Virunga Massif, despite our recognition of the wealth of biodiversity, high levels of endemism, and considerable conservation challenges in the wider region (Plumptre et al., 2007). However, we do know that amphibians of the Virunga Massif are likely threatened by climate change and habitat loss, degradation, and fragmentation (Plumptre et al., 2007;Salerno et al., 2017). We can also predict that anuran distribution in the Virunga Massif is at least partially determined by elevational gradients, as is the case in nearby regions in East Africa (Malonza & Veith, 2011;Poynton, Loader, Sherratt, & Clarke, 2006;Zancolli et al., 2014). And finally, we predict to find positive effects of habitat size on anuran species richness, diversity, and abundance, as has been shown for other anuran communities in the tropics (Almeida-Gomes, Vieira, Rocha, Metzger, & Coster, 2016).
Volcanoes National Park (VNP) in Rwanda, which protects a part of the Virunga Massif, experiences enormous anthropogenic pressures. Most impacts stem from a particularly high human population density near the park (Bush, Ikirezi, Daconto, Gray, & Fawcett, 2010). The volcanic soil found in and near the park allows for highly productive agriculture, which drives local people to cultivate up to the edge of the park (Weber, 1987). In turn, direct water extraction from the park, alterations of ground water levels through agriculture activities, widespread use of nonnative plant species that may affect hydrological regimes (e.g., Eucalyptus spp.; Dye, 2013), and altered frequencies of rainfall (Campos et al., 2017), possibly contributed to observed reductions of water availability inside wetlands (The Dian Fossey Gorilla Fund International, unpublished data). Over the last decades, previously permanently wet ponds in certain wetlands now remain dry for parts of the year, while other wetlands have dried up completely and are experiencing associated vegetational succession toward shrubbery and forest. This process of wetland desiccation is likely to have a negative impact on the anurans particular to this habitat (McMenamin, Hadly, & Wright, 2008), especially in combination with changes in climatic factors, such as altered frequency and periodicity of rainfall (Hartter et al., 2012).
Previous studies in VNP showed the presence of at least nine different anuran species (Roelke & Smith, 2016). Many, if not all, of these species have some association with wetland habitat, especially with regards to reproduction. Of these species, the Karisimbi tree frog Leptopelis karissimbensis is currently considered Vulnerable (IUCN SSC Amphibian Specialist Group, 2016).
We investigated the distribution of these species in wetlands of varying sizes, found along an elevational gradient of c. 600 m that spans all vegetation zones in VNP. In contrast to previous studies conducted in the Albertine Rift region, which address broadscale distribution patterns across vegetation types (Behangana, Kasoma, & Luiselli, 2009;Malonza & Veith, 2011;Zancolli et al., 2014), we focus on fine-scale patterns within broadly similar vegetation types (wetlands). In particular, we provide an overview of the anuran community composition of each wetland and test for the separate effects of elevation and wetland size on anuran distribution. This information provides a baseline to test what drives anuran distribution and anuran responses to increasing threats in this biodiverse, but fragile, region of Central Africa. 9°C at lower elevations to 0°C at higher elevations, and monthly variation in rainfall can be substantial (Campos et al., 2017).

| Data collection and study area
The vegetation in VNP includes mixed forest, bamboo, Hagenia woodland, herbaceous, brush ridge, subalpine and alpine zones (Plumptre, 1991), whereas the vegetation inside the focal wetlands is characterized by species such as Juncus oxycarpus and Mariscus karisimbiensis, with encroachment by the dryland species Hypericum revolutum and Yushania alpina at many of the study sites. Focal wetlands ranged in size from approximately 99 to 45,355 m 2 and were found between 2,546 and 3,188 m.a.s.l. (Table 1).
Every study year, we visited each wetland two to five times between November and March and searched for anurans using the Visual Encounter Survey method (Crump & Scott, 1994) while walking 5 m wide strip transects that were spaced 15 m apart. During each visit, three or four observers slowly walked each transect parallel to each other, gently separating vegetation (e.g., clusters of grass) and looking at all sides of leaves. We did not monitor waterbodies for tadpoles, juveniles or adults, and though we used vocal cues to guide us to the location of individuals we did not abandon our gradual search method to record calling individuals (instead recording them as we reached them). The length and number of transects were consistent for each wetland across years but differed between wetlands to account for wetland size. To adjust for these varying sampling efforts across wetlands, we opted to use "site visits" as the sampling units in further statistical analyses, while sampled area (as a proxy of sampling effort) was used as an offset in our models (see below). Finally, we note that we were obliged to conduct all visits during daylight hours, despite the potential underestimates of abundance and richness this may have generated, because of restrictions on evening use of the park.
F I G U R E 1 Location of focal wetlands (example top left photo) in which we sampled anurans (example Leptopelis karissimbensis bottom left) within Volcanoes National Park, Rwanda. The map shows the boundaries of the park as well as the adjacent Virunga National Park (Democratic Republic of Congo) and Mgahinga National Park (Uganda), which together protect a large part of the Virunga Massif Unfortunately, we were unable to collect tissues and adopt molecular techniques for specimen identification. That said we did aim to collect photographic evidence to maintain a base of material for further clarification of identities (e.g., by consulting external experts). For the first three years of our study, we aimed to photograph the dorsal side of each specimen, a method we had to abandon in 2015 due to logistic restrictions. We thereafter aimed to photograph at least one specimen of each species found during every visit to a wetland. Despite these efforts, we acknowledge that both confusing taxonomy and morphological similarity between various congeneric species may have introduced some inaccuracy regarding the identities of certain specimen, see Discussion. For taxonomy, we followed Frost (2018), whereas we used information from for example, Zimkus, Rodel, and Hillers (2010), Roelke, Greenbaum, Kusamba, Aristote, andSmith (2011), andChanning, Dehling, Lotters, andErnst (2016) to confirm that our species records were plausible in terms of geographic distribution.

| Statistical analyses
To estimate species richness for each wetland and to assess sampling efficiency, we generated sample-based observed species richness rarefaction curves with 95% confidence intervals, using the "vegan" (Oksanen et al., 2007) package in R (R Core Team, 2014).
For this, we used the accumulative data for all years and site visits.
In addition, we calculated estimated species richness as the mean of four abundance-based species richness estimators (CHAO, JACK1, JACK2 and Bootstrap), following 999 randomizations of observed species richness. We subsequently calculated sampling completeness as the ratio between observed and estimated species richness.
We estimated species diversity using the Shannon-Wiener index, which takes species richness and the relative proportions of species found in each wetland into account.
After testing for potential collinearity among variables using Pearson's correlations, and finding none, we tested for the effects of wetland size using generalized linear mixed effects models (glmm) constructed with the "lme4" package (Bates, Maechler, Bolker, & Walker, 2015) for models with abundance and species richness as response variables. We used "glmmTMB" (Brooks et al., 2017) for zero-inflated glmms with Shannon-Wiener diversity as response variable. In these models, we treated wetland size as a fixed effect, sampling area (as a proxy of sampling effort) and elevation as offsets, and survey year and wetland identity as random effects. Next, we looked at the separate effect of elevation, for which we constructed similar models with abundance, species richness and Shannon-Wiener diversity as response variables, but with elevation as a fixed effect, survey year and site identity as random effect terms, and sampling area and wetland size as offsets. Appropriate error structures were applied for all models, and we used Chi-square to compare models to null models.
To assess the similarity of wetlands in terms of their species composition, we used a nonmetric multidimensional scaling (NMDS) ordination analysis using Bray-Curtis pairwise distances based on standardized, square-root-transformed abundance data (to reduce the influence of the most dominant species). We visualized results of this analysis in an ordination plot and used a permutational multivariate analysis of variance (adonis in the "vegan" package) to test for differences in Bray-Curtis similarity between wetlands.

| RE SULTS
Between 2012 and 2017, we recorded 2,638 specimens of at least eight anuran species (Tables 2 and 3), of which at least two are

Albertine Rift Endemics (ARE). The most abundant species was
Hyperolius castaneus (ARE), which was found at all sites, whereas we recorded only two specimens of Sclerophrys kisoloensis. Given the morphological similarity between Leptopelis karissimbensis (ARE) and L. kivuensis (ARE), we acknowledge that certain individuals that we identified as L. karissimbensis in the field could actually be L. kivuensis, and as such we opted use L. karissimbensis/kivuensis for further analyses. In addition, although we are certain that we have recorded both Phrynobatrachus graueri and P. parvulus, it was difficult to identify certain specimen of Phrynobatrachus to species-level in the field.
We thus also merged these two species as P. graueri/parvulus for further analyses. For reference, a list of all species records can be found  Figure 2).

| D ISCUSS I ON
Although we urge that our results should be interpreted with caution, due to biases introduced by skewed sampling efforts, such as a lack of nighttime sampling (see below for further discussion), we emphasize that there is considerable spatial variation among the anuran communities of the study wetlands in VNP. We found we confirm that anuran richness, diversity, and abundance are likely to decline with elevation in the wetlands of VNP, which shows that these trends exist even within one relatively homogeneous habitat type (wetlands) across a relatively small elevational range. Anuran distributions are thus not only determined by a turn-over between habitats, but potentially also by the effects of elevation, a proxy for climatic differences, per se.
We did not test for geographic and mid-domain effects on the elevational distribution of anurans as we focused on fine-scale trends among populations sampled across a relatively short gradient (~600 m) at mid-to-high elevations (i.e., we had "truncated sampling," TA B L E 2 Summary of anuran sampling in ten focal wetlands in Volcanoes National Park, Rwanda, between 2012 and 2017. Given are the total abundance relative abundance (per hectare sampled), observed (Sobs) and estimated (Sest) species richness (± SD), and proportion of estimated species detected (Sobs/ Sest) in each wetland.
richness and abundance with an increase in elevation (Behangana et al., 2009).
Few studies have shown patterns of decline in anuran abundance within species across elevational gradients (but see e.g., Behangana et al., 2009), but they logically extend from patterns observed for species richness and diversity (for a review of potential causes for these patterns see McCain & Grytnes, 2010). As species approach the limits of their elevational range, reproductive capacity tends to decline, despite remarkable physiological adaptations and forms of acclimations of many high-elevation anurans (Brattstrom, 1968;Christian, Nunez, Clos, & Diaz, 1988;Navas, Carvajalino-Fernández, Saboyá-Acosta, Rueda-Solano, & Carvajalino-Fernández, 2013). At lower temperatures and arguably under the influence of changes in rainfall and water availability, anurans tend to be developmentally, physiologically and behaviorally limited (Navas et al., 2013;Putnam & Bennett, 1981). As a result, they have shorter breeding seasons, take a relatively long time to mature (Amat & Meiri, 2017), and produce fewer (but larger) eggs (Morrison & Hero, 2003) with low abundances. In turn, this leads to potentially high vulnerability (i.e., higher local extirpation probability) of populations living at high elevations (Morrison & Hero, 2003).
We show that larger wetlands contain higher abundances of anurans and host more species, which strengthens the notion that desiccation and wetland size reduction will be detrimental to anuran survival in VNP. Not all of these species are restricted to wetland habitats, but most will likely depend on these wetlands as spawn sites and for early life stages (Roelke & Smith, 2016). More research will be needed on the nature of species-area relationships, the natural history, and breeding micro-habitat requirements (i.e., where the various species spawn) of anurans in VNP, if we are to understand future responses of these species to declines in wetland size through desiccation. For example, we recommend exploring whether habitat heterogeneity is higher in larger wetlands, which may drive higher species richness, or whether larger wetlands contain a larger random sample of species from the regional species pool (Lomolino, 2000).
Although we deem it probable that both elevation and wetland size play a role in observed differences among anuran communities across VNP wetlands, we acknowledge that various additional factors could also have contributed to these patterns. For one, there is a host of both local (e.g., species interactions) and landscape-level factors (e.g., characteristics of the surrounding vegetation) that can influence differences in anuran distribution across our focal wetlands (see e.g., Van Buskirk, 2005). Isolation of populations (Scheffer et al., 2006) and stochastic processes (reviewed in Marsh & Trenham, 2001) could similarly affect dispersal, local population persistence, and local extinction, and thus shape local anuran communities.
We also acknowledge that some of our estimates of richness and diversity may be confounded due to the lack of clarity on the taxonomy and phylogeny of East African anurans (see e.g., Channing et al., 2016;Roelke et al., 2011;Zimkus & Schick, 2010). This problem is exacerbated by the uncertainty regarding the identity of at least some of the individuals we found. Thus, we emphasize that our estimates of abundance, richness, and diversity TA B L E 3 Anuran species and number of specimens found in each focal wetland in Volcanoes National Park, Rwanda, between 2012 and 2017 During field work, we identified that both P. graueri and P. parvulus are likely to be present in eight wetlands, thus there is likely at least one more species in these wetlands than presented here.

Species
F I G U R E 2 Species accumulation curves constructed using sample-based rarefaction curves for ten focal wetland sites F I G U R E 3 Effects of elevation and wetland size on anuran abundance, species richness and Shannon-Wiener diversity, showing a linear regression (black line) ± SE (gray shade) should be interpreted in the light of the spatial patterns that we explored and should not be seen as exact estimates.
Although Roelke et al. (2011) provided an overview of the anurans present in VNP, and others have addressed aspects of anuran community composition of either the larger Albertine Rift region (Plumptre et al., 2007) or a specific wetland (Lac Ngezi) in VNP (Sinsch, Greenbaum, Kusamba, & Lehr, 2011), our study is the first to provide data on anuran community composition across different wetlands in the park. Our results suggest that species composition is highly variable across wetlands (Table 3), illustrated by the fact that we were unable to detect all nine species in one single wetland.
Although these differences in local community composition and relative abundances of varying species could stem from sampling biases, they also suggest that every single wetland in VNP is worth conserving. Indeed, we found that wetlands differed significantly in anuran community composition, and at least some species are likely to be found only in a few, small wetlands and are uncommon any-  (Sinsch et al., 2012), but is now limited due to logistic restrictions. One possible solution to logistic limitations to nocturnal surveys would be to conduct automated acoustic sampling, which would not require our presence in the field at night (Acevedo & Villanueva-Rivera, 2006).
The main challenge for conservation of montane anurans will be to understand the response of individual species to climate change (Beniston, 2003;Duarte et al., 2012), drying of wetlands (McMenamin et al., 2008), and a host of other threats. Certain species may be able to adapt to novel environmental conditions (Duarte F I G U R E 4 Nonmetric multidimensional scaling (NMDS) ordination of anuran community assemblages at different focal wetland sites, based on square-root transformed, standardized abundance data (stress 0.14). Each point represents accumulative data for the 2012-2017 study period for a single transect, and the number of transects (i.e., points) ranges between one and four for different focal wetland sites Hoffmann, Chown, & Clusella-Trullas, 2013), whereas others may require more dramatic conservation actions, such as habitat restoration or translocations (Carvalho, Brito, Crespo, & Possingham, 2010;Shoo et al., 2011). A first step in setting the stage for adequate protection would be the implementation of long-term monitoring protocols. Here, we report how such long-term monitoring reveals spatial patterns in anuran species distribution in VNP wetlands, findings which serve as a baseline for future monitoring in a region where even basic descriptions of the natural history of most species is still lacking.

ACK N OWLED G M ENTS
We thank the Rwanda Development Board (RDB) for granting us access to conduct fieldwork in Volcanoes National Park. We are very grateful to Biodiversity Program field team of The Dian Fossey Gorilla Fund's Karisoke Research Center for their hard work in the long-term anuran monitoring program. We thank C. Roelke and two anonymous reviewers for their many useful recommendations on initial drafts of this manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
DT conceived of the study, designed the study, collected the data, and helped draft the manuscript; WE contributed to the design of the study, data interpretation and manuscript writing; MAD performed data analysis, interpreted analyses, and helped draft the manuscript; NG aided with data formatting, analysis, and interpretation as well as drafting the manuscript; YvdH interpreted data analyses, and drafted the final manuscript. All authors provided comments and approved the final manuscript.

DATA ACCE SS I B I LIT Y
Species abundances and richness data, used in this manuscript, are available in