Habitat amount and ambient temperature dictate patterns of anuran diversity along a subtropical elevational gradient

Patterns of diversity along elevational gradients are driven by species characteristics but remain poorly understood. Filling this gap is imperative given the deteriorating conservation status of anurans worldwide. Here, we examine frog diversity and species composition along a sharp subtropical elevational gradient and assess the degree to which these are determined by environmental and spatial predictors.

Proposed explanations generally invoke multiple drivers that can potentially act and interact at different scales (Szewczyk & McCain, 2019;Willing & Presley, 2016). Variation in contemporary climatic conditions, such as temperature and precipitation, may directly filter out organisms given their physiological tolerance, leading to changes in species richness and composition (i.e. a positive climate-diversity relationship; Fu et al., 2006;McCain, 2009).
Warmer climates mostly lead to higher primary productivity, resulting in more available resources, which in turn ensure coexistence of more individuals (McCain, 2007). Therefore, productivity can indirectly promote diversity by maintaining more viable populations and reducing local extinction rates, also known as the more-individuals-hypothesis (MIH: Binkenstein et al., 2018;Srivastava & Lawton, 1998). Diversity also scales to the available habitat area, and montane regions covering larger areas or providing more available habitat should also harbour more species and individuals (Fahrig, 2013;Rosenzweig, 1995). This area-diversity relationship is based on the assumption that larger habitat areas reduce extinction rates and are more likely to encompass more habitat types (MacArthur & Wilson, 1967;Rosenzweig, 1995). Biotic processes such as competition, habitat complexity and habitat heterogeneity across ecotones have also been proposed to explain elevational patterns of diversity (Beck et al., 2017;McCain & Beck, 2016;McCain et al., 2018). However, the contribution of these factors to montane diversity remains unclear given the difficulty in measuring these processes along natural elevational gradients.
However, the AF is a highly threatened global biodiversity hotspot with only 26% of its original forest cover remaining, most of which scattered across thousands of fragments <50 ha (Rezende et al., 2018;Ribeiro et al., 2009). Virtually all extensive tracts of AF are in montane regions, rendering them a conservation priority (Tabarelli et al., 2010). However, there have been few attempts to elucidate the elevational patterns of diversity throughout the AF (e.g. Almeida-Neto et al., 2006;Martinelli, 2007;da Silva et al., 2018).
Atlantic Forest anurans are especially diverse, with 625 recorded species, 85% of which are endemics .
This high level of diversity has been attributed to the exceptional heterogeneity of macro-and microhabitats, topography and climatic variation across this biome (Brown & Brown, 1992;Haddad & Prado, 2005). However, this frog diversity is threatened by the relentless impacts of habitat fragmentation, habitat loss, climate change and infectious diseases, all of which elevate anurans to the most threatened vertebrate taxon worldwide (Becker et al., 2007;IUCN, 2019;Stuart et al., 2004). Despite the vast remaining extent of montane AF habitat, investigations on the variation in frog diversity along elevational gradients remain scant, and available data are restricted to a few sites and narrow elevational ranges (350-750 masl; Giaretta et al., 1997;Giaretta et al., 1999;Goyannes-Araújo et al., 2015). Given the rapid pace of habitat conversion and degradation, basic information on AF frog diversity is therefore crucial to help develop evidence-based conservation strategies (Pereira et al., 2013).
Here, we investigate the patterns and potential drivers of frog diversity along an extensive subtropical elevational gradient in the southern Atlantic Forest. First, we describe how frog species richness and abundance changes with elevation. We expected mid-elevation peaks in frog diversity, which is consistent with the trend for all major vertebrate taxa (McCain & Grytnes, 2010). Second, we assess the role of environmental and spatial variables on frog species richness and abundance by simultaneously testing predictions from six major hypotheses of diversity: area, climate, habitat complexity, habitat amount, MIH and productivity (Table 1). Finally, we describe the compositional variation in frog assemblages along the elevational gradient and assess the contribution of environmental predictors and species relationships on the structure of these assemblages.

| Study area
The largest remaining AF montane forest habitat span the Serra do Mar, Serra da Mantiqueira and Serra Geral ridges of southeastern and southern Brazil, where elevations range from sea level to ~2,900 m (Gontijo-Pascutti et al., 2012). This study was carried out within the Serra Geral, a complex montane landscape formed K E Y W O R D S altitudinal gradient, Atlantic Forest, elevational patterns, environmental filter, frogs, habitat integrity mainly by Early Cretaceous lava flows, which are associated with the break-up of Gondwana, separating the modern continents of South America and Africa (Hartmann, 2014). The ~920,000 km 2 Serra Geral region extends from Argentina through Uruguay, Paraguay and southern to south-eastern Brazil (Frank et al., 2009). Specifically, we focused on the ~49,800 ha São Joaquim National Park (hereafter, SJNP), a strictly protected area located in the state of Santa Catarina, Brazil ( Figure 1). The study area spans approximately 1,200 km 2 of dissected terrain and altitudes ranging from 200 to 1,822 masl (Vianna et al., 2015). The vegetation comprises of a mosaic of dense evergreen forest largely occurring up to 800 masl, replaced by mixed Araucaria forest and high-altitude grasslands at higher altitudes (IBGE, 2012). The climate is humid subtropical without a dry season with some variation depending on elevation: hot summers and mean annual temperature of 19°C below 800 masl, and cold winters with frequent frosts and a mean annual temperature of 13.5°C above 800 masl (Alvares et al., 2014).

| Sampling design
We sampled frogs along 38 discrete lentic aquatic environments (ponds, humid areas, marshes and small dams; hereafter, ponds) located between 312 and 1796 masl (Figure 1). We acknowledge that some specialist anuran species (i.e. ground-dwelling, stream-dwelling and bromeliad breeders) may have been undersampled, but ponds are considered suitable environments given that most AF frogs use them and their vicinities to reproduce and forage (Crump, 2015;Haddad et al., 2013). We employed the following steps to select survey ponds. We first subdivided the elevational gradient into eight 200-m elevational bands, from 200 to 1,800 masl, which ensured that pond selection was not biased for or against any single region of the elevational gradient (McCain & Grytnes, 2010 (Valeriano & Rossetti, 2012). We randomly selected five ponds for each elevational band that (a) were at least 1 km apart; (b) were surrounded by negligible anthropogenic habitat disturbance; and (c) were no farther than 5 km from the nearest access road.
Prior to sampling, we visited all selected ponds to confirm their presence and accessibility. Due to restricted availability, only four ponds were selected between 1,000 and 1,200 masl and between 1,600 and 1,800 masl. A detailed description of each pond can be found in Appendix S1.1 in Supporting Information.

| Frog survey
We quantified frog species richness and abundance at each pond based on 10 monthly samples from September 2017 to August 2018, except for December 2017 and April 2018. We combined visual and acoustic encounter surveys, within a 5-m buffer along the pond perimeter and recorded all frogs on the water surface, above and underneath any vegetation, rocks and logs (Scott-Jr & Woodward, 1994). A 5-m survey buffer was selected do maximize visual and acoustic anuran detectability, avoiding double counting.
Sampling initiated at least one hour after dusk (17:30 hr-19:18 hr) and ended no later than 01:00 hr. Sampling effort was standardized by one person-hour per pond per month. Each pond therefore accumulated a sampling effort of 10 person-hours, amounting to an overall effort of 380 person-hours of active search time. Up TA B L E 1 Major hypotheses and predictors of frog diversity along elevational gradients considered in this study

Hypothesis Theory Predictor Scale
Area Diversity is positively related to area (MacArthur & Wilson, 1967;Rosenzweig, 1995) Pond surface area (PSA) Local

Elevational band area (Area) Landscape
Climate Diversity is positively related to temperature (Allen et al., 2002) Mean annual temperature (MAT) Landscape Habitat complexity Diversity is positively related to habitat complexity (Rosenzweig, 1995) Water surface cover (WSC) Local

Pond vegetation structure (PVS) Local
Habitat amount Diversity is positively related to total habitat amount (Fahrig, 2013) Pond habitat amount (PHA) Local Band habitat amount (BHA) Landscape More-individuals-hypothesis Abundance is positively related to productivity, and therefore, species richness is positively related to abundance (Binkenstein et al., 2018;Srivastava & Lawton, 1998)

| Environmental and spatial predictors
We measured environmental and spatial variables at both local and landscape scale. Ponds and their immediate surroundings (within a 500-m radial buffer) were defined as the local scale. At the landscape scale, we opted for a different approach from most previous altitudinal gradient studies because our sampling protocol was not restricted to a single survey site per evenly spaced elevational band. We therefore followed Peters et al. (2016) in  Habitat complexity was measured at the local scale, using two widely used proxies based on previous studies on neotropical frog diversity (e.g. Melchior et al., 2017;da Silva et al., 2012). The first was the proportion of pond water surface covered by aquatic and/or emergent vegetation (WSC). The second was the proportion of pond margins covered by shrubs and/or trees, which we used as a proxy of pond vegetation structure (PVS), with higher values indicating more complex habitats. WSC and PVS values were determined using aerial photographs (summer 2018) taken by a rotary-wing drone (Phantom 3 Standard, DJI, China) and subsequently analysed using the ImageJ software (Rasband, 1997). Habitat amount was considered as the percentage of natural vegetation cover measured at both the pond (PHA) and the elevational band (BHA) scales. We included these predictors because even though SJNP is a protected area, private agricultural landholdings are still active inside the park boundaries.
Additionally, some ponds were also located outside the SJNP area and had been exposed to some anthropogenic activities. PHA and BHA values were obtained using previously classified, georefer-

| Data analysis
Species were assumed to occur at all elevations between the lowest and highest elevations at which they were detected, a common range interpolation approach in elevational gradient studies (Rowe & Lidgard, 2009). This can be considered the most realistic approach to understand diversity patterns as species tend to be distributed continuously (Rowe & Lidgard, 2009). We also used total species richness and abundance based on all surveys at each pond, given that overall sampling effort was equal at each pond. We acknowledge that local diversity may change through time and data pooling may lead to overestimated species richness and abundance. However, few individuals are usually found more than once in a reproductive event, which is associated with the generally low recapture rates within year compared to between years (e.g. Nomura et al., 2012;Tucker, 1995). Therefore, we considered this a valid approach in representing frog diversity at the pond scale, and any temporal variation in frog diversity will be investigated elsewhere.
We assessed our data reliability by evaluating sample completeness by means of sample coverage (Chao et al., 2014). We also estimated asymptotic species richness using the species abundance data at each pond using the iNEXT R package (Hsieh et al., 2016). We then used Pearson correlations to examine observed and estimated species richness values.
Prior to analysis, we examined a predictor correlation matrix to reduce issues of multicollinearity. Mean annual temperature, observed temperature and precipitation were highly correlated (Pearson r = +0.93, -0.98 and -0.90, respectively; Appendix S2.1) so we opted to retain the former as a proxy of climatic conditions. Precipitation values are relatively high along the SJNP elevational gradient (range = 1,463-1,817 mm/yr; see Appendix S1.1 for predictor values along the elevational gradient), and given the marked precipitation aseasonality (absence of a dry season), we assumed that rainwater is not a limiting factor in the study region. Additionally, MAT was also highly negatively correlated with PHA and BHA (r = -0.75 and r = -0.81, respectively; Appendix S2.1), but we opted to retain MAT and both habitat amount variables in order to assess the purposed hypotheses.
We examined the distribution of frog species richness and abundance along the elevational gradient using generalized additive models (GAMs), which incorporate more flexible smoothing functions in modelling nonlinear relationships compared to generalized linear models (Hastie & Tibshirani, 1990). We used a thin plate regression spline on elevation, set the data family to Poisson for richness and negative binomial for abundance (to account for overdispersion) and allowed scaling parameters to be estimated in model fitting.
Smoothing parameters were estimated using restricted maximum likelihood (REML) and basis dimension (k) set to five to prevent overparameterization (Peters et al., 2016). GAMs were computed using the mgcv R package (Wood, 2011). The presence of any residual spatial autocorrelation in both species richness and abundance was tested using Moran's eigenvector maps based on 42 different spatial weighting matrices (Bauman et al., 2018), but none could be found (see details in Appendix S3).
We used structural equation modelling (SEM; Grace, 2006) to simultaneously test the direct and indirect effects of multiple hypothesized predictors of frog species richness and abundance. SEM constructs a network of causal relationships, where pathways and their directions indicate previously formulated hypotheses (Grace et al., 2010). Using this framework, we developed an initial full model encompassing the hypothetical a priori causal relationships among the measured predictors and response variables (Appendix S4.2).
First, we hypothesized that productivity (NPP) would be directly and positively influenced by temperature (MAT). We would also expect that frog abundance was directly influenced positively by habitat complexity (WSC and PVS), area (PSA and Area), habitat amount (PHA and BHA) and productivity, but also via a direct or indirect effect of temperature (mediated by productivity). Lastly, we expected that these same predictors would directly influence positively frog species richness, but also indirectly, via indirect abundance-mediated effects.
SEM was applied in a Bayesian framework with local estimation of parameters to account for specific error distribution on each response variable using the brms R package (Carpenter et al., 2017). Therefore, productivity, abundance and species richness were modelled with gamma, negative binomial and Poisson error distribution (all log link), respectively. For each model, we ran four chains with 50,000 iterations each, with the first 25,000 discarded as burn-in and the remaining thinned by 100, so that 1,000 samples from the posterior distributions were retained for analysis. We used brms default weekly informative priors for parameters estimate (Lemoine, 2019). Chains convergence was assessed by ensuring that potential scale reduction factor (Rhat) values were less than or equal to 1.01 for all parameters (Bürkner, 2017). We considered relationships important when predictors' highest 95% posterior density intervals (HPDIs) exclude zero. In addition to raw estimates, we also computed paths coefficients standardized on relative ranges (Grace et al., 2018) and calculated the Bayesian R 2 for each response variable (Gelman et al., 2019).
Classic approaches to analyse multivariate species data do not account for proper mean-variance relationships of abundance, nor provide sufficient quantitative information on the contributions of different community assembly processes, which may lead to misleading interpretations Warton et al., 2012).
Therefore, we employed a joint species distribution model (JSDM; Pollock et al., 2014) with the inclusion of latent variables (LVs) to simultaneously describe the pattern of compositional variation and estimate the contribution of predictors on assembly processes.
JSDMs extend multivariate generalized linear models by using LVs, which ensure separating species co-occurrences given (dis)similar responses to known environmental predictors (e.g. environmental filtering), to unexplained species associations ("residual correlation"; Ovaskainen et al., 2017;Warton et al., 2015). These residual correlations can either represent responses to unmeasured predictors or the outcome of biotic interactions (D'Amen et al., 2018;Warton et al., 2015).
We fitted JSDMs using the boral R package (Hui, 2016), which uses Bayesian Markov chain Monte Carlo (MCMC) sampling to estimate model parameters. First, we fitted a pure latent variable model (LVM) including only two LVs, which can be understood as a model-based unconstrained ordination Warton et al., 2015). From this model, we extracted the posterior median LVs values and used as unconstrained ordination axis coordinates to visualize the variation in species composition across ponds. We then fitted a second model by adding the same covariates used to model species richness and abundance (MAT, NPP, BHA, PHA, PVS, Area, WSC and PSA), but also included vegetation type (either dense evergreen forest or mixed Araucaria forest) where ponds were located, plus two LVs. We also extracted the posterior median LVs values to visualize pond-scale species composition after controlling for environmental predictors (residual ordination). We then quantified correlations in frog species' abundance due to shared environmental responses and their respective residual correlations (responses to unknown covariates or biotic interactions) (Ovaskainen et al., 2016;Warton et al., 2015). Lastly, we performed variation partitioning among the included predictors to quantify their relative contribution in structuring frog assemblages (Ovaskainen et al., 2017).

Predictors were centred and standardized prior to inclusion into
JSDMs. For this analysis, we only included species that occurred in at least three ponds to avoid misinterpretations (n = 28; Warton et al., 2015). We ran models using the default uninformative normal priors for model parameters. The number of iterations on MCMC chains was set to 200,000 with a 100,000 burn-in period and a thinning factor of 100. A negative binomial error distribution was used to account for overdispersion in modelling species abundance simultaneously. We also included a site random effect to account for any possible spatial autocorrelation (Hui, 2016). Model convergence was assessed using a combination of the Geweke diagnostic (Geweke, 1992) and visualization of trace plots, while Dunn-Smyth residuals plots were used to inspect violation of model assumptions (Hui, 2016). Correlations were considered important whenever their corresponding 95% HPDIs excluded zero.

| RE SULTS
We recorded a total of 12,636 frogs across the 38 surveyed ponds, In the mechanistic SEM proposed with the direct effects of temperature, productivity, area, habitat amount and habitat complexity on frog species richness, plus their indirect effects though frog abundance, only the proportion of natural habitat remnants within elevational bands (BHA) was considered an important predictor of richness (i.e. 95% HPDIs excluded zero; Figure 4).
Standardized path coefficient indicates that local frog species richness increases in 45% of its observed range values when moving from areas with minimum to maximum BHA (observed range: 57% -99%). Similarly, mean annual temperature (MAT) was the only important predictor of frog abundance (Figure 4). Moving from colder to hotter areas (range: 11.7°C -18.5°C), frog abundance increases by 43% of its observed range. Lastly, even though we found that mean annual temperature positively influences net primary productivity, our data provide no support for the more-individuals-hypothesis. Neither the direct nor indirect effect of productivity through abundance was important in predicting local  Pond sites showed a clear elevational gradient in species composition ( Figure 5a). The unconstrained ordination showed a separation along the second latent variable axis, with a tight cluster of ponds at the bottom left and another sparser cluster spanning the top along the first latent variable. The first cluster was formed mainly by ponds at lower elevations (300-800 masl), whereas the second consisted of mid-to high-elevation ponds (>800-1,800 masl), thereby reflecting the observed altitudinal gradient (Figure 5a). The residual ordination indicates that most of the variability among ponds were no longer detected once the environmental predictors were included ( Figure 5b). Therefore, it was not possible to identify any distinct clusters based on elevational differences, but a few ponds remained more segregated from the core of the ordination (Figure 5b).  0 0 1 8 0  We detected both positive and negative important correlations (i.e. 95% HPDIs exclude zero; n = 106) due to similar responses to measured predictors for several species (Figure 5c), indicating that environmental filtering was important in producing correlations in species abundance. Of these, 57% were positive (n = 60) and 43% (n = 46) were negative, indicating that observed correlations in species abundance were mostly attributed to convergences rather than divergences in species-specific responses to predictors. In contrast, we observed a greatly reduced small set of important residual correlations after accounting for environmental predictors (n = 16), all of which were positive (Figure 5d). The joint species distribution model indicated that all measured predictors combined accounted for 77% of the variation in species abundance (averaged over species; Figure 5e). Of these, mean annual temperature was most influential in frog species composition, as on average 22% of the variance in species abundance was attributed to this covariate (Figure 5e).
Other predictors contributed individually to less than 10% of the variance.

| D ISCUSS I ON
To the best of our knowledge, this is the first comprehensive study documenting transitions in anuran diversity along a steep elevation   Elevation gradient (range of 1,500 m) of the Brazilian Atlantic Forest that simultaneously tests multiple hypotheses explaining these patterns.
In this study, frog species richness showed a hump-shaped distribution while frog abundance was high in the lowlands but decreased at mid-to high elevations. The variation in frog species richness was best explained by the available area of natural habitats within elevational bands, whereas changes in mean annual temperature best explained frog abundance along the elevational gradient. In addition, we observed marked compositional divergence between frog assemblages from lowland to mid-to-high-elevation sites, which was also best explained by changes in temperature.
The mid-elevational peak in species richness we found here is consistent with previous studies on vertebrate elevational ranges, for which this pattern is considered to be the most pervasive (e.g. Yet, there is still no consensus on the operational mechanisms promoting this pattern, which appears to be driven by a complex interplay of different predictors (Beck et al., 2017;Lomolino, 2001).
We evaluated the support for most frequently proposed hypotheses explaining elevational patterns of species richness (e.g. Laiolo et al., 2018;McCain et al., 2018;Wu et al., 2013). Surprisingly, frog species richness was positively related to only landscape-scale habitat amount. This indicates that ponds located in areas containing more native vegetation should also harbour more frog species. This provides direct support for the habitat amount hypothesis (HAH; Fahrig, 2013), rather than the more general species-area relationship. The HAH was also found to be important in predicting taxonomic diversity in small mammals and functional diversity in frogs in highly fragmented Atlantic Forest landscapes (Almeida-Gomes et al., 2019;Vieira et al., 2018).
Areas corresponding to more extensive natural vegetation are more likely to harbour a higher diversity of suitable (micro)habitats and thereby confer a higher chance of accommodating habitat specialist species (Henneron et al., 2019;Leibold et al., 2004). Many frog species are considered habitat specialists and their local occupancy depends exclusively on some environmental conditions. For instance, Ischnocnema henselii exhibits direct development and depends on moist leaf-litter to lay their eggs (Haddad & Prado, 2005).
Similarly, Fritziana mitus is found in epiphytic bromeliads, which they use for shelter, foraging, ovipositing and raising their tadpoles (Walker et al., 2018). At SJNP, these two species were only found where their specific microhabitat requirements were present.
This highlights the strict dependence of some species on specific microhabitat types, which can partly explain the mid-elevational species richness peak documented here. At the elevational band scale, overall natural vegetation cover ranged from 57% to 99% (mean = 70±12%, SD; Appendix S1.1) and was positively correlated with elevation (r = 0.81, p < .001). This positive effect therefore indicates sufficiently strong environmental filtering in excluding species from unsuitable habitat areas along the elevation gradient.
Like most topographically dissected landmasses worldwide, lowland Atlantic Forest areas succumbed to the heaviest impacts of agricultural land use, particularly compared to adjacent high-elevation forest lands (Tabarelli et al., 2010), and this reflects the strong positive relationship uncovered here between elevation and primary habitat amount. Moreover, we cannot rule out the possibility that lowland areas in the past retained higher anuran species diversity.
For instance, water body density (e.g. ponds) throughout the entire 785,400-hectare planimetric area of our study region decreases at increasingly higher elevations (see Appendix S9.3). Water body density can also be seen as proxy of suitable habitat conditions for anurans, given that many species use them to reproduce (Crump, 2015;Haddad et al., 2013). We would therefore expect that higher lowland density of water bodies could also lead to higher species packing. Wachlevski et al. (2014) recorded twice as many frog species (n = 32) within similar lowland ranges (240 m -490 masl) at Serra do Tabuleiro State Park, a protected area ~75 km east of our study landscape. Serra do Tabuleiro is one of the largest AF remnants in southern Brazil and most of its lowlands still remain relatively intact (MMA, 2000). In this context, historical land use trajectories-that are stacked against the integrity of contemporary lowland habitatslikely contribute to the observed pattern of frog species richness.
In general, larger areas tend to harbour more individuals, thereby reducing local extinction rates, and promote higher levels of habitat diversity and heterogeneity (Rosenzweig, 1995). However, habitat area per se (i.e. pond size and area of elevational bands) was a poor predictor of frog species richness, which is consistent with montane regions in China (Fu et al., 2006;Hu et al., 2011;see also Bueno et al. (2019) for a recent discussion on island species-area relationships). We found no support for the habitat complexity hypothesis even by taking into account many distinct montane habitats spanning a broad spectrum of habitat complexity (Figure 4; Table S1.1).
Although local habitat complexity was not an important predictor of frog diversity, previous studies have found that higher local habitat complexity promotes species diversity, which is attributed to greater numbers of suitable sites for anuran egg-laying, vocal activity and shelter against predators (e.g. Vasconcelos et al., 2014). Species inhabiting low-complexity habitats may exhibit temporal segregation in pond use (Vasconcelos & Rossa-Feres, 2005), and this resource partitioning could partly explain the co-occurrence of similar frog diversity compared to more complex habitats. In summary, local habitat complexity seems to be species and/or system dependent and is not an essential element for pond-dwelling frog diversity along the elevational gradient.
Our data also indicate that neither temperature nor a proxy of habitat productivity was important predictors of frog species richness. This implies that, although anurans exhibit a critical thermal minimum lower than other ectothermic vertebrates (Buckley et al., 2012), their physiological limitations do not necessarily restrict their distribution along the elevational gradient. Indeed, some subtropical Brazilian frog species are known to be reproductively active all year-round, even showing physiological and behavioural adaptations to exploit cold environments (Kiss et al., 2008;Rossa-Feres & Jim, 1994). Should these limitations hold true, we would expect a monotonic upslope decline in frog species richness, which does not reflect our data. Productivity was high at low elevations but was not followed by a monotonic upslope decline, even though ~40% of the variance in productivity was explained by ambient temperature. This suggests that available potential energy may not be limiting along the elevation gradient, which could lead to similar levels of food availability for frogs. In addition, freshwater pond environments are likely to buffer upslope declines in local productivity and none of the ponds we surveyed ever freeze even during severe winters, even though widespread frosts are common above 800 masl. We acknowledge that we did not quantify food resource availability directly. However, both arthropod species richness and abundance, which comprise the bulk of anuran diets (Carvalho-Rocha et al., 2018;Solé & Rödder, 2010), are positively related to productivity (Schuldt et al., 2019), which may support this hypothesis. Frog species richness was apparently decoupled from frog abundance and productivity, lending little if any support to the "more-individuals-hypothesis" in our study system, although further studies on food availability along elevational gradients are needed.
Frog abundance was also decoupled from productivity but was positively affected by temperature. Lower temperatures may impose prohibitive limitations on anuran reproductive success, even in coldadapted species (Kiss et al., 2008;Navas et al., 2013), which may result in small populations at high-elevation sites. For instance, lower temperatures can impose changes in species phenology and restrict breeding seasons to shorter periods, but can also increase the time to complete larval metamorphosis and affect tadpole survival (Amat & Meiri, 2017;Bernal & Lynch, 2013). While these mechanisms can explain the lower frog abundance at lower temperatures, we highlight that more physiological and behavioural studies are needed to help us better understand the mechanisms behind these patterns.
The latent variable model showed that lowland frog assemblages were very distinct compared to those in the highlands.
Differentiation among frog assemblages along the elevational gradient was also observed in palaeotropical montane areas in Tanzania (Zancolli et al., 2014), China (Hu et al., 2011) and Borneo (Menegon & Salvidio, 2005). This observed segregation is consistent with variation in vegetation macromosaics across our study region. This suggests that each vegetation formation retains its own pool of species adapted to local environmental conditions and exploiting specific microhabitats (Crump, 2015;Haddad & Prado, 2005 Our results confirm that local abiotic conditions, mainly temperature, play an important role in predicting frog assembly composition. For example, Scinax tymbamirim, Leptodactylus gracilis and Rhinella abei were common sympatric species at warmer sites (lower altitudes; Figure 2) and were strongly positively correlated with one another and this climatic predictor. These species were also negatively correlated with species that are also typically sympatric at cold, high-elevation sites, such as Boana leptolineata, S.
granulatus and Physalaemus gracilis (Figure 2). We detected just a few positive relationships between frog species that could be associated with positive interactions. Competition among frog species, specially noncongeners, is often acknowledged to be rare (Cloyed & Eason, 2017;Duellman, 1999), and general assembly process is more related to responses to environmental conditions (e.g. Ernst & Rödel, 2006;Indermaur et al., 2010), corresponding to a Gleasonian community worldview (Gleason, 1926). This strongly suggests that frog assemblages along the SJNP elevational gradient are formed by species with similar thresholds of thermal tolerance.
Although frog species richness and assemblages were associated with their own key predictors, both were conditional on effects of environmental filtering, providing further evidence that this process is the main driver of different levels of frog diversity (Werner et al., 2007). Several studies have shown that climatic gradients are important drivers of the spatial distribution of neotropical anurans at large biogeographic scales (Luiz et al., 2016;Silva et al., 2014;Vasconcelos et al., 2014). However, we have shown that steep climatic variation imposed by elevational settings is also important in shaping frog assemblages, even at a relatively narrow spatial scale.

| CON CLUS IONS
Patterns of heterotherm diversity along tropical elevational gradients remain poorly explored. Our study shows that different components of anuran species diversity are influenced by mechanisms related to environmental filtering along the elevational gradient.
Specifically, frog species richness was strongly associated with the proportion of natural habitats remaining within each elevational zone, while frog abundance was highly associated with mean annual temperature. This reinforces the detrimental effects of the widespread historical degradation of lowland areas in the Atlantic Forest, which is directly associated with the leading cause of declines in frog diversity worldwide. Climatic conditions were strong determinants of frog assemblage structure along the elevational gradient. Lowland assemblages clearly diverge from those at midto high elevations, leading to two distinct groups of frog species.
Nonetheless, future studies considering species functional traits and phylogenetic relationships may help elucidate if climatic restrictions are due to niche conservatism or species-specific environmental filtering operating on life history characteristics.
Lastly, the entire elevational gradient we studied was clearly important in retaining the overall regional scale (gamma) diversity.
This indicates that anthropogenic habitat disturbance has exerted an enormous overarching influence on montane frog diversity, reinforcing the need for effectively protected areas that are well distributed across all elevational zones.

ACK N OWLED G EM ENTS
We thank the staff at São Joaquim National Park for logistical support. We are grateful to Silvia Onofre and all members of

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13187.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw data are provided in Appendix S1.1 and S5.3.

B I OS K E TCH
Vítor Carvalho-Rocha is interest in diversity patterns and biodi-