Lurking in the depth: Pond depth predicts microhabitat selection by Rhinella icterica (Anura: Bufonidae) tadpoles at two different sampling scales

Habitat selection has long been a central theme in ecology and has historically considered both physiological responses and ecological factors affecting species establishment. Investigating habitat selection patterns at different scales can provide important information on the relative roles of the environmental factors influencing the organisms’ abilities to use their surrounding habitat. This work aimed at investigating which environmental factors determine habitat selection by Rhinella icterica tadpoles, and also took the opportunity to investigate how the scale in which tadpoles and environmental data are sampled might influence the habitat selection results. A total of 2.240 tadpoles were counted in the whole sampling area, and while substrate cover and depth were the variables that better explained the abundance of tadpoles at the larger scale (plot level), depth and water turbidity better explained tadpoles’ abundance at the smaller scale (subplot level). The results suggest that avoiding predation by matching the background color is a likely process explaining tadpoles’ occupancy at both scales. Depth is known to influence tadpole habitat use in the tropics, and although its combination with turbidity and substrate cover varied between scales, our study suggests that sampling at different scales might not affect the inferred ecological processes driving habitat selection. This information might also be useful to predict tadpoles’ responses to micro-environmental perturbations and help in guiding the choice of parameters that should be taken into account when analyzing the effects of habitat degradation in Atlantic Forest amphibian populations.


Introduction
There is currently little doubt that habitat loss and degradation are among the major threats to biodiversity worldwide (Venter et al. 2006;Betts et al. 2017). Habitat selection has long been a central theme in ecology and has historically considered physiological requirements and other limiting factors, such as competition and environmental conditions, as important factors affecting species establishment (Huey 1991). Therefore, investigating habitat selection patterns can help biologists to more efficiently determine how these abiotic and biotic characteristics will influence the individual's abilities to use their surrounding habitat (MacArthur and Levins 1964;Whitham 1978), and evaluate the effects of habitat degradation on these patterns (Serrano and Astrain 2005;Campbell et al. 2018).
In ectotherms, the biotic and environmental factors affecting habitat selection, such as predation pressure and temperature, will certainly influence population dynamics over time (Huey 1991;Bradford et al. 1992). For amphibians, and specifically their aquatic tadpoles, habitat occupancy depends on many environmental characteristics of the water bodies they inhabit, such as temperature, oxygen concentration, depth, substrate coverage and surrounding aquatic vegetation (Laurila 1998;Nie et al. 1999;Torres-Orozco et al. 2002;Provete et al. 2014;Clevenot et al. 2018). These factors will directly affect the capacity of tadpoles to thermoregulate, which is linked to their abilities to select the microhabitat that provides better foraging opportunities (Noland and Ultsch 1981), and also influence their abilities to colonize other habitats (Clevenot et al. 2018). Furthermore, tadpoles' ecomorphological traits variation within a population can be largely explained by the availability of microhabitats (Jordani et al. 2019) and resource partitioning (Eterovick and Fernandes 2001). Hence, intra-specific phenotypic plasticity in response to micro-environmental variations will also influence tadpole habitat use (Fatorelli and Rocha 2008;Jordani et al. 2019). Taken together, these patterns indicate that understanding the mechanisms behind tadpole habitat selection may help explain the factors that determine amphibian population structure and fluctuations, especially in the megadiverse Neotropical region.
The Atlantic Forest is one of the most diverse and understudied ecosystems in the world (Morellato and Haddad 2000;Ribeiro et al. 2009). This diversity is mirrored in amphibian reproductive patterns, and the Atlantic Forest harbours the largest number of amphibian reproductive modes known to date (Haddad and Prado 2005). Our understanding of amphibian reproductive modes and its implication for conservation, population ecology, and community structure has advanced in the last decades (Crump 2015), and this development was accompanied by a continuous growth of research on tadpole ecology (Altig 2018). Only approximately 30% of the anuran species in the world even had their tadpoles described  and, although ~60% of Brazilian tadpoles have already been described (Provete et al. 2012), most of what we know on habitat selection involves choices made by reproducing adults (Buxton and Sperry 2017).
Researchers working on tadpole ecology and habitat selection must always take practical decisions on how to effectively sample tadpoles in the field. In most cases, individuals are sampled near the stream or pond edges using dip-nets (Laurila 1998;Both et al. 2009;Provete et al. 2014), or using localized traps (Torres-Orozco et al. 2002), but efforts to sample the entire extension of the water bodies are rarely employed (Jordani et al. 2019). Many factors might influence these practical decisions, such as the particular characteristics of field sites and the local water bodies, as well as the time and resources available for sampling. However, sampling at the edges of the water body might bias habitat selection inferences, since space occupancy by tadpoles might be reflecting the particular characteristics of these places, and not taking into account how tadpoles effectively use the habitat considering the whole range of environmental conditions and biotic interactions. For example, in ponds or slow flowing streams, edges might be more shaded and the substrate might be more covered with leaves because of the surrounding vegetation. Since the results found in habitat selection studies could be influenced by the scale on which the data was collected (Morris 2003;McGarigal et al. 2016), sampling tadpoles and environmental data at different scales might uncover different ecological patterns. At larger scales, we can expect that environmental variables that are more related to tadpole space use, such as depth and water flow, will be more important to explain the distribution of individuals in the habitat. On the other hand, at smaller scales, especially close to the edge of the ponds, predation might be a more important factor affecting the distribution of individuals (Heyer et al. 1975), and habitat structure such as the presence of surrounding vegetation and leaves in the substrate will determine the presence of such individuals (Kopp et al. 2006). Hence, habitat selection patterns at larger scales will reflect the processes to which these individuals are responding in order to forage (Pearman 1993;Both et al. 2009;Domingos et al. 2015) whereas, at smaller scales, will reflect individuals behavior in order to avoid predation (Bridges 2002).
The toad Rhinella icterica (Spix, 1824) is distributed throughout southern South America, from eastern Paraguay to northern Argentina, in southern and southeastern Brazil, and northwards to the state of Bahia (Frost, 2019). This species differs from other Rhinella marina group species mainly by the presence of a subtriangular parotoid and tibial glands, as well as less developed foot webs, and tadpoles with isometric growth (Kwet et al. 2006;Lima and Pederassi 2012). Besides being a large-bodied anuran, surprisingly little is known about this species in terms of reproduction (Pereyra et al. 2015) and habitat selection (Guix et al. 1998). In the Atlantic Forest, this species can be found from 30 to 1870 meters above sea level, and adults were found occupying a wide range of natural and anthropized habitats from streams and ponds to cleared and native forests (Guix et al. 1998). Its reproductive period varies according to the local climate and, although calling males and amplectant couples can be found throughout the year, reproduction usually peaks from July to February (Bertoluci 1998;Rodrigues and Bertoluci 2002;Lima et al. 2010). Oviposition and tadpole development occur in ponds or slow-flowing streams (Dixo and Verdade 2006;Hartmann et al. 2010), and there is anecdotal evidence that, in each particular year, R. icterica tries to avoid interspecific competition by reproducing earlier than other Rhinella species (Lima et al. 2010).
In this paper, we investigated which environmental factors determine microhabitat selection of R. icterica tadpoles. In addition, we examined how spatial scale might influence habitat selection results. We hypothesized that the environmental variables explaining tadpoles' abundance will differ according to the sampling scale: variables related to the physical characteristics of the pond will be more important at the larger scale, whereas variables that provide better camouflage opportunities will be more important at the smaller scale. Thus, at the larger scale, we predict tadpole abundance will be greater in shallow and larger habitats, with higher substrate cover and slower water flow. These expectations are based on the assumption that in these conditions there are more opportunities for shelter and less risk of dragging. At the smaller scale, near to the pond edge, higher levels of vegetation cover, water turbidity, as well as substrate leaf content, will possibly provide better camouflage opportunities and a stronger role in the abundance of tadpoles.
The sampling site was a 30 m long temporary pond with width varying from 60 cm to approximately 6 m and maximum depth of approximately 30 cm. This pond was surrounded by vegetation, predominantly secondary forest, and emergent macrophytes in the pond margins. Since we sampled during the rainy season (August), the pond had a constant influx of water from a nearby creek, and the excess of water would escape from the opposite part of the pond through a small water channel excavated by the rain. Ultimately, this means that the pond had constant running water throughout its whole extension.

Data collection
With the help of labelled stakes we delimited 20 plots along the pond extension where it was possible to visually count the tadpoles from the margin and without entering the pond. This condition was necessary to avoid biases due to the possibility of disturbing the tadpoles if an observer entered the pond. Each plot was 1 m long, consisting of the whole rectangular area across the pond, and this delimitation comprised our larger scale data collection. For the smaller scale data collection, each plot was subdivided in two equal parts considering a central imaginary line parallel to the pond edge (Fig. 2). In other words, the larger scale measurements considered the whole area of each of these one-meter long plots within the pond, while the smaller scale considered only one side of each plot, from the center to the pond edge. Hereafter, we refer to the larger scale sampling units as 'plots', and to the smaller scale as 'subplots' (Fig. 2).
Within each plot and subplot we counted the total number of tadpole individuals (abundance), measured pond width, depth, water flow velocity, and visually estimated water turbidity, surrounding vegetation cover percentage (including emergent macrophytes present in pond margins), and percentage of the sand-soil covered by leaves (hereafter substrate cover). Tadpoles were counted by a single observer, and always considering one subplot and then the other (i.e., the plot count was the sum of both subplots). Tadpole counting was performed from the outer edge of the pound before all other measurements to avoid disturbing the tadpoles. Hence, after counting the tadpoles, water flow velocity was measured in meters per second (m/s) by the formula D/T, where D is the distance (meters) and T is time (seconds). Distance was measured using a polystyrene foam ball and a 30 cm scale, where the ball was gently placed on the water surface, parallel to the streamflow direction, and next to the scale for visualization. Pond width and depth were measured (in cm) with a measuring tape. Since the pond was very shallow, it was not possible to use the standard Secchi disk protocol to estimate water turbidity. Thus, we used a Secchi disk as a proxy to estimate turbidity by visually estimating how well it could be seen (similar to a 'percentage' of how well the disc could be visualized). Turbidity, surrounding vegetation cover, and substrate cover were estimated by the same observer.
The width was measured only once (for the plot), and the subplot width was half of the plot width. To replicate the usual data collection in habitat selection studies, all other measurements were taken twice (depth, velocity, turbidity, vegetation cover, and substrate cover), one in the plot center and another one in the subplot center. Since Rhinella tadpoles, if undisturbed, tend to remain still for long periods Perotti 2009, 2010), we believe our tadpole abundance counts were accurate, and we did not observe any tadpole disturbance before entering the pond to measure physical variables. Data were collected in a single day on 22 August 2019. After all data collection procedures, a few tadpoles were collected from each plot to confirm identification. We based our identification in the comparative data provided in Tolledo and Toledo (2010).

Data Analysis
We separately analyzed the data of the plots (n = 20) and the subplots (n = 20) to compare the habitat selection results inferred in these two different ecological scales represented by the different sampling strategies. Although each plot consists of two subplots, we only analyzed the data of one randomly chosen subplot per plot, to try to re-flect the most commonly adopted tadpole sampling strategy used in ecological studies.
We checked for multicollinearity of our data by calculating Variance Inflation Factors (VIF) and using a threshold of VIF > 5 to exclude variables (Sheather 2009). Since no variable exceeded this threshold at the plot or subplot scale we ran all analyses using our full dataset. Considering the geographic structure of our sampling strategy, before proceeding with the analyses we tested for spatial autocorrelation in the data using Moran's I statistic by implementing a permutation test with 199 simulations in the package spdep (Bivand et al. 2008). Because abundance data at both scales showed spatial autocorrelation (pseudo-p values < 0.05), we used statistical procedures to remove such effects from the residuals of our multiple regressions analyses. Thus, for each scale, we performed Wavelet-Revised Models (WRM) with tadpole abundance as the dependent variable and environmental data as independent variables (width, depth, water velocity, turbidity, surrounding vegetation, and substrate cover). The WRM is an extension of Generalized Linear Models (GLM) developed to remove the effect of spatial autocorrelation in multiple regressions (Carl and Kühn 2010). This method is particularly suited for the analyses of our data for a multitude of reasons: first, GLM methods such as WRM do not require data to be normally distributed and are also robust to deviations from other parametric assumptions; second, the calculations can be fitted to data in different statistical distributions and; third, the WRM was developed to accommodate grid-based datasets, which is in line with the structure of our data. All WRM analyses were performed using package spind (Carl et al. 2018). In addition, to build the best possible models, we investigated the best fit of our data to a statistical distribution so we could use it in WRM model inference. Us- ing package fitdistrplus (Delignette-Muller and Dutang 2015) we visually inspected our data and compared the fit to the three different distributions available to run the WRM in package spind (normal, binomial or poisson).
After building the full WRM model, we performed an automated backward stepwise model selection using the Akaike Information Criterion (AIC) to select the environmental variables that better explain tadpole abundance. In short, the analysis starts with a previously created full WRM model where abundance is explained by all environmental variables, and variables are removed in a stepwise fashion while the AIC value of the models are compared until the best model fit with the lowest AIC value is reached. Thus, for each scale, the best model explaining the abundance of tadpoles was selected by AIC, and then the significance of the variables of the best model was calculated. All statistical procedures were performed in R v3.6.0 (R Core Team, 2019).

Results
In total, we counted 2.240 tadpoles in the whole sampling area. Tadpole abundance at the plot scale was almost twice that of the subplot scale (Table 1). In general, environmental variables had higher values at the plot scale, except for substrate cover that was slightly higher at the subplot scale. Turbidity did not vary between scales (Table 1). Visual inspections suggested that the best distribution that should fit our plot abundance data was a normal distribution. We compared the normal, binomial and poisson distributions using the Kolmogorov-Smirnov statistic (respectively 0.213, 0.232 and 0.450) and also Cramér-Von Mises statistic (respectively 0.120, 0.179, and 1.33). For our subplot data, visual inspections also suggested that the best distribution was a normal distribution. Again, we compared the normal, binomial and poisson distributions using the Kolmogorov-Smirnov statistic (respectively 0.249, 0.254 and 0.494) and the Cramér-Von Mises statistic (respectively 0.191, 0.267 and 1.426).
Hence, we fitted both our WRM models using a normal distribution. All R codes and statistical results, and the visual inspection graphs, can be found at Suppl. Material 1: R Markdown file.
At the plot scale, the stepwise model selection procedure kept substrate cover and depth as the variables that better explained abundance, while depth and turbidity were the ones selected at the subplot level (Table 2). Contrary to our expectations, the relationship of tadpole abundance and depth was positive, since there was an increase in the number of tadpoles in deeper plots (Fig. 3A). Nonetheless, the increase in tadpole abundance in plots with higher substrate cover percentage is in agreement with our hypothesis for the plot scale (Fig. 3A). In the subplot scale, tadpoles' abundance increased with turbidity and depth (Fig. 3B). While the relationship with water turbidity is in agreement with our initial hypothesis, we did not predict that depth would influence abundance at this scale.

Discussion
Habitat selection theory predicts that individual choices will influence foraging efficiency and survivorship, which will be directly related to individual fitness (Rosenzweig 1981;Chesson 2000;Morris 2003). Understanding tadpole habitat use might help to evaluate how habitat change will influence the underlying ecological processes related to individuals' choices (Holbrook and Schmitt 1988;Carey et al. 1992) and, thus, provide important information on how environmental disturbance might directly affect survivorship. Even though this work focuses on only one species, our results reveal useful patterns that help to fill the knowledge gaps we still have on Atlantic Forest tadpole ecology and habitat selection in particular.
Many environmental variables have been found to influence tadpole habitat selection such as water flow and substrate composition (Hoff et al. 1999;Haramura 2006). Nonetheless, water depth is usually a very strong component influencing habitat occupancy by tadpoles both at population (Pearman 1993) and community scales (Both et al. 2009). For our population of R. icterica tadpoles, depth was positively related to abundance both at the plot and subplot scales. The other two variables that better explained tadpole abundance were substrate cover (plot scale, Fig. 3A) and water turbidity (subplot scale, Fig. 3B). Taken together, all these variables seem to be essentially related to predation avoidance by the tadpoles which, in our particular system, can be at least partially explained by the tadpoles' dark body color. The sand-soil at the bottom of our sampled pond was clear, in stark contrast to the Rhinella tadpole color. However, by staying on top of dark decomposing leaves at the bottom (the main substrate cover), the tadpoles would be less conspicuous to visually oriented predators (Nomura et al. 2013;Gontijo et al. 2018). In fact, while measuring our environmental variables we noticed that, when disturbed, tadpoles would quickly swim through the clear sand-soil until other leaves in the nearby substrate were reached, or they would also swim to deeper parts of the pond. In line with these strategies, more turbid waters also provide better camouflaging opportunities from visually oriented predators (Gardner, 1981;Minello et al., 1987). Since predation pressure has an important role in tadpole habitat use (Eterovick and Barata, 2006), being in the deep, in more turbid water, and on top of leaves, could be directly influencing tadpole predation avoidance.
There is no information on the predation of R. icterica tadpoles in the wild, but it can be inferred from results of other species of the genus Rhinella that they should be mainly consumed by invertebrate predators such as dragonfly larvae and waterbugs (Jara and Perotti 2009), fishes (Nomura et al. 2011), and birds (Beckmann et al. 2011). In laboratory experimental conditions, Rhinella tadpoles only slightly change their activity patterns in response to the presence of predators (Jara and Perotti 2010), which suggests that they may rely on other ecological strategies to escape predation in the wild, such as camouflage.
Another experiment showed that although matching the background coloration did not affect tadpole mortality rate, tadpoles also decreased their activity in response to the presence of predators, and were less active in white background (Nomura et al. 2013). These laboratory habitat selection results depict a generalization of the individual's behavioral changes, and seem to partially corroborate our conclusions that adjusting their coloration to the background is the process driving habitat selection in R. icterica tadpoles.
Although the selection of both substrate cover and turbidity by our analyses can be attributed to protection against predators (Semlitsch 1993), it is worth noting that these characteristics might also be correlated with food availability (Skelly et al. 2002;Schiesari 2006;Provete et al. 2014). Rhinella tadpoles primarily eat algae scraped from the bottom of the water body (Hinckley 1963;Rossa-Feres et al. 2004) and the presence of leaves in the bottom is also linked to higher primary productivity and decomposition rates (Skelly 1996;Skelly et al. 2002), which provide food resources to tadpoles. Pond depth may also be determinant in tadpoles' survival and growth, since deeper ponds tend to retain water for longer periods, allowing tadpoles to grow bigger until metamorphosis (Peltzer and Lajmanovich 2004). However, the different combination of these variables in the plot and subplot level, and the fact that depth was selected in both scales, suggests that the more proximal explanation, i.e., decreasing predation risk, might be a better explanation for the process affecting R. icterica tadpoles' habitat selection. Nonetheless, given the lack of information on tadpole habitat use in the Atlantic Forest, further studies aiming at evaluating these and other factors influencing tadpoles' habitat use at different spatial scales are highly warranted.
According to our results, water flow was not important for tadpoles' abundance. However, this result needs to be interpreted with caution because of the low variation of this data observed among the replicates (Table 1), and the fact that flow was not particularly high throughout the whole pond. Not surprisingly, this observation supports the idea that R. icterica adults prefer slow-flowing water for reproduction, reducing the risk of egg dragging Figure 3. Relationship between tadpole abundance and pond depth at both sampling scales. In the plot level (A), the substrate cover percentage was also selected as an important variable determining tadpole abundance and is represented by the gradient color scale from light to dark brown. In the subplot level (B), water turbidity was also selected as an important variable determining tadpole abundance and is represented by the gradient color scale from light to dark blue.  (Lima et al. 2010). During fieldwork, we also sought for tadpoles in all the nearby creeks and other water bodies in the reserve, but with no success. We noted that these other water bodies were clearer, had fast-flowing water, and were mostly free of substrate leaves, reinforcing the mating site preferences of R. icterica adults. In tropical ponds, characteristics such as water flow and pond size have distinct effects on tadpole distribution (Borges Júnior and Rocha 2013), but tadpoles that inhabit ponds with slow water flow generally use sites with vegetation cover (Eterovick and Barata 2006). Predation rates are supposedly lower when the habitat has more hiding opportunities (Kopp et al. 2006), but vegetation cover did not influence the abundance of R. icterica tadpoles in our study site, also supporting the idea that they might predominantly rely on camouflage to avoid predation. Even though we were not able to test it, another factor that could influence tadpole microhabitat use is their toxicity. Previous studies detected toxins in the eggs and tadpoles of Rhinella marina, although in lower concentrations compared to adults (Flier et al. 1980;Akizawa et at. 1994;Hayes et al 2009). For instance, R. marina tadpoles are toxic to some native predators and other tadpole species in Australia, where it is an invasive species (Covacevich and Archer 1975;Hamley and Georges 1985;Crossland and Alford 1998;Crossland and Azevedo-Ramos 1999;Crossland 2000). The toxins in the tadpoles' bodies make them unpalatable, which, in turn, may favor a gregarious behavior by inclusive fitness (Waldman 1982;Griffiths and Denton 1992). Besides, the black coloration in unpalatable tadpoles might also be related to aposematism (D'Heursel and Haddad 1999). However, we did not find other tadpole species in our study site, which prevented any comparison with co-occurring species. Unfortunately, we cannot infer the impact of R. icterica toxicity in tadpoles' spatial distribution and microhabitat use from our data and observations, and the currently available evidence indicates that camouflage is possibly the main factor driving habitat selection in this case. Nonetheless, R. icterica tadpoles do have gregarious behavior and a body coloration that is consistent with aposematism. From an ecological perspective, we suggest that further studies should focus on tadpole individuals' spatial distribution, group sizes, and the comparison of such characteristics with sympatric species.
Although our two different sampling scales might reflect common sampling practices in tadpole ecology studies, they were probably not divergent enough to capture different ecological properties of tadpole habitat selection. Contrary to our hypotheses, all three variables selected at the plot and subplot scales are putatively related to predation avoidance, even though at least one environmental variable differed. In this regard, we suggest that sampling at a smaller scale should be enough to capture the ecological processes behind habitat selection, and future studies in tadpole habitat selection might benefit from this information during the research planning stage. Nonetheless, given that different variables were selected in each scale, our results also indicate that tadpoles are sensitive to micro-environmental changes, which suggests that small scale habitat perturbations might influence the behavior of these animals. Tadpoles might be highly sensitive to environmental disturbances (Babini et al. 2018), and our results highlight the need to consider these putative influences when planning conservation actions for Atlantic Forest amphibians.
In conclusion, although variables that predict R. icterica tadpoles' habitat selection in the Atlantic Forest might be explained by the different ecological processes discussed above, avoiding predation seems to be a reasonable explanation for our results, considering that water depth was selected at both scales and correlated with substrate cover and turbidity. Our study also indicates that sampling at different scales might not affect the inference of the ecological processes behind habitat selection. In face of the concerning conservation status of the Atlantic Forest, and considering the lack of data on the ecology of many amphibian species, the conservation of anurans reproductive sites might be one of the most important conservation measures to maintain viable populations in the Atlantic Forest (Gonçalves et al. 2015). In this context, our study can also be helpful in guiding the choice of parameters that should be taken into account when analyzing the effects of habitat degradation in amphibian populations.