The effect of urbanisation and local environmental heterogeneity on phenotypic variability of a tropical treefrog

Urbanisation reduces species richness and change community composition. However, little is known on how the phenotype of organisms with low dispersal ability respond to environmental changes associated with urbanisation in fast urbanizing centres, such as those in the Global South.

Here, we tested how urbanisation rate, local environmental heterogeneity, land surface temperature, and spatial gradients affect phenotypic traits associated with dispersal, resource acquisition, and performance, namely: body size, head shape, and leg length of the Dwarf Treefrog (Dendropsophus nanus) using a space-for-time substitution approach.

We took linear measurements from 768 individuals in 21 ponds along an urban gradient in central Brazil. We also measured local environmental variables and summarized them using Hill-Smith Principal Component Analysis. The spatial arrangement of ponds at multiple scales was described using Moran Eigenvector Maps. Those variables were entered into a Structural Equation Model to test their direct and indirect effects on the mean and coefficient of variation (CV) of phenotypic traits. Additionally, we calculated the Scaled Mass Index as a proxy for fitness and estimated the adaptive landscape for body size, size-free leg length, and head shape. We also tested for spatial autocorrelation in traits.

Body size decreased from the periphery to the urban centre, whereas CV of body size and head shape had the opposite pattern. Body size increased, whereas CV of body size and head shape decreased in man-made ponds. The CV of leg length decreased with increasing land surface temperature. The remaining traits were not affected by any predictor variable. None of the traits were spatially autocorrelated. Both body size and head shape were under weak directional selection, but in opposite directions.

Our results suggest that the lack of a clear spatial variation in phenotypic traits can be due to a weak selection, due to a recent, although intense, urbanisation process. In conclusion, eco-evolutionary dynamics in tropical cities seem to have a different pace compared to temperate ones. Our results can contribute to building urban ecological theory that explicitly includes city age, their development, growth rate, and history.

Methods
Study site and sampling design

We conducted fieldwork in Campo Grande, Mato Grosso do Sul, central Brazil (Projection: UTM–21S; Coordinates: 751114, 7731659; DATUM=SIRGAS/2000; 529 m a.s.l.; Figure 2). Campo Grande has been founded 127 years ago (Arruda 2006; PLANURB 2019) and has about 897,938 citizens (IBGE, 2023). Although recent, urbanisation has been fast, with the city gaining ~ 13,000 ha of urban infrastructure from 1985 to 2021 (MRS, pers. obs. based on data from MapBiomas, 2023). The urban area has about 35,903 ha with most human population living in the urban centre (PLANURB, 2019). The climate is Equatorial Savanna or Köppen’s Aw, with dry winters and rainy summers (Kottek et al., 2006). Vegetation is composed mainly by Cerradão, Cerrado, Campo Sujo, and Vereda vegetation types (Ribeiro & Walter, 2008).

We sampled adults of the hylid frog Dendropsophus nanus in 21 ponds (Figure 2; Figure S1) along an urban gradient through surveys at breeding sites (Scott Jr. et al., 1994) and visual encounter surveys (Crump & Scott Jr., 1994). Field work was conducted between 17:30 and midnight from November 2021 to April 2022 and from November 2022 to January 2023 visiting two ponds per day. We did not visit all ponds in all months; instead, each pond was visited until we obtained 30 individuals to guarantee the same sampling effort.

Phenotypic traits

We took the following linear measurements in the field using hand digital calliper (MTX 150 mm) to the nearest 0.01 mm: body size (Snout-Vent Length – SVL), head width, head length, and leg length (foot, thigh, and tibiofibula) following Watters et al. (2016). Relative leg length is related to vulnerability to predation (Citadini et al., 2018; Emerson, 1985a) and dispersal ability (Phillips et al., 2006), while head length and width is related to size and variety of feeding resources (Emerson, 1985b; Parmelee, 1999). Body size is related to thermoregulation, dispersal, and desiccation resistance (Wells, 2007). After measurements, every frog was tagged with visible implant elastomer (VIE) and released back.

Prior to analysis, we transformed head width, head length, and leg length to remove the effect of body size by applying a linear regression on each trait as a function of body size and then retained the residuals. To obtain head shape, we took the log-shape ratio (Mosimann, 1970) by computing the geometric mean of head length and head width to obtain their sizes. Then, we divided each value by their geometric means and log-transformed the result (Claude, 2013). Finally, we performed a Principal Component Analysis and took the first axis to represent head shape. PC1 retained 92.5% of the variance and was negatively correlated with head length (-0.74) and width (-0.67). Negative scores along PC1 represent individuals with wider and longer heads, while positive ones represent those with shorter and narrower heads. Analysis was performed in R v. 4.2.3 (R Core Team, 2023).

Fitness proxy

We weighted all individuals using a Pesola (Lightline 10 g #10010) to calculate the Scaled Mass Index (SMI), following Peig & Green (2009). This index relates body mass and size using a Major Axis Regression to obtain a measure of body condition, which was used here as a fitness proxy. Analysis was conducted in lmodel2 R package (Legendre, 2018). The SMI assumes a non-linear growth relationship that removes the effect of body size. This index has been used to access body condition of amphibians, especially in urbanized environments (Franco-Belussi et al., 2024), since it provides information about animal health status and is positively related with lipid and protein content (MacCracken & Stebbings, 2012). Lipid content in turn is correlated with survival in amphibians (Scott et al. 2007) but body condition also seems directly correlated with survival in a wild toad population (Reading 2007). Field measurements of fitness components are challenging for many taxa, including frogs. We acknowledge that using SMI as a fitness proxy can be problematic (Wilder et al., 2016). However, evidence exists linking lipidic storage to energetic reserves and ultimately to physiological status in frogs. In this context, SMI can be considered a performance index (but see Franklin & Morrissey, 2017), since animals with good body condition that have energetic reserves provided by a past successful feeding are more likely to survive (Jakob et al., 1996; Wells, 2007).

Urbanisation and local environmental heterogeneity

We measured urbanisation as the percentage of building area and roads (Szulkin, Garroway, et al., 2020) in a 500-m buffer around each pond (Figure S2). We also extracted this same variable from raster layers from 1985 to 2021 provided by MapBiomas (2023). To obtain urbanisation rate, we subtracted the current (2021) urban area from that of 1985 and divide it by 36 (number of years in the period). We repeated all analyses using a 1-km buffer and results did not differ. Spatial data handling was conducted in the R packages terra (Hijmans, 2023) and landscapemetrics (Hesselbarth et al., 2019).

To quantify local environmental heterogeneity, we measured the following variables at pond level: pond area (m2; using Google Earth®), pond depth (mean of five points per pond; in m), pond hydroperiod (permanent or temporary – binary variable), presence or absence of cattle (binary variable), aquatic predators (vertebrate or both invertebrate and vertebrate – categorical variable), margin profile (excavated, flat or sloping – categorical variable), water temperature (using a thermometer submerged at 1 m for 1 min; in ºC), and percentage of floating vegetation (visual estimation %). Ponds with excavated margins were artificial, man-made ponds. Likewise, ponds with vertebrate predators and cattle were more associated with rural environments (Figure S1). The type of pond margins limits the amount and kinds of perches available for frogs, influencing how they use the environment (Vasconcelos et al. 2009), likely constraining the range of phenotypic traits associated with them. Those were later summarized using a Hill-Smith Principal Component Analysis (Hill & Smith, 1976), which allows combining discrete and quantitative variables. The PC1 retained ~ 63% of the variance in the data. The environmental variables that most contributed to PC1 (correlation > 0.6) were: excavated margin (positive), temporary hydroperiod (negative), invertebrate predators (negative), and cattle presence (negative; Table S1, Figure S3 and S4). Analysis was conducted in the ade4 R package (Dray & Dufour, 2007).

Lastly, we calculated land surface temperature (LST) and urban heat island (UHI; Figure S5). We used the following satellite images: LANDSAT 8 thermal band 10 for 2022 (07 November 2022) and LANDSAT 5 thermal band 6 for 1985 (17 November 1985). All images were obtained from Earth Explorer (United States Geological Survey – USGS: https://earthexplorer.usgs.gov/). The pixels of these images contain solar light variables, such as reflectance and radiance. To extract land surface temperature, we used the equation proposed by Coelho (2013) and Coelho & Correa (2013), which converts digital numbers to solar radiance, radiance to brightness temperature, and then to temperature in Kelvin or Celsius. Then, we used land surface temperature to calculate urban heat island (UHI) (see Ahmed et al., 2013; Faisal et al., 2021) for the whole urban area, as follows:

UHI=(Ts-Tm)/SD

where, Ts is land surface temperature; tm is the mean of land surface temperature, and SD is its standard deviation. Finally, we built an urban heat island map and calculated maximum, minimum, variance, and mean temperature in 500-m buffers around each pond. Mean land surface temperature was later used as predictor variable in the path analysis. Data handling and extraction were done in QGIS v. 3.22.1 (QGIS.org, 2020).

Building spatial predictors and spatial autocorrelation

To describe the spatial arrangement of ponds at multiple scales, we used Moran’s Eigenvector Maps (MEMs). Firstly, we used the geographical coordinates of ponds to build a spatial neighbourhood network. We tested different types of networks that represent different hypotheses of connections between ponds (i.e., Delaunay triangulation, Gabriel graphs, Minimum spanning trees, and a custom-built connectivity network). Then, we built a spatial weighting matrix in which weighting was a linear function of the inverse of distance. Afterwards, we generated a set of MEMs for each spatial neighbourhood network. Finally, we performed a selection of spatial matrices in the adespatial R package (Dray et al., 2023). We retained the best subset of MEMs that represents the autocorrelation of body size, built from the best spatial neighbourhood, which was selected based on the Akaike's Information Criterion corrected for small sample sizes (AICc). Then, we computed the corresponding Moran’s I statistic of each MEM and calculated its significance using a Monte Carlo permutation test to retain only eigenvectors with positive and significant autocorrelation.

The best spatial model was the custom neighbour network with six positive and significant MEMs (Figure S6). MEM3 had the lowest AICc and was used to represent pond spatial arrangement in the path analysis. Positive scores along MEM3 were arranged in a Northeast-Southwest direction and comprised most of the urbanized ponds (Figure S2 and S6). Finally, to test for spatial autocorrelation in traits and environmental variables, we used Moran's I correlograms (Legendre & Legendre, 2012) using the same custom-built spatial neighbourhood network.

Space-for-time substitution approach

We use this approach to make inference about microevolutionary processes. By using space as a substitute for time, we can use present environmental conditions to understand patterns produced by a past process (Wogan & Wang, 2018). Here, we sampled ponds along a rural-urban gradient (space) and calculated their rate of change in urbanisation from 1985 to 2021 (36 years) (time). Urbanisation rate varied from 0% to 2.11% per year. Ponds in the core urban had high urbanisation rates and had already undergone a moderate building process before 1985, providing evidence that ponds in urban core experienced a long lasting change. Thus, we assume frogs living in highly urbanised ponds have experiencing the effect of urbanisation for more time than those in peri-urban or rural sites. Provided tropical treefrogs live approximately 12 years (Stark & Meiri, 2018), at last three generations have passed since 1985.

Data analysis

We checked for multicollinearity of all environmental variables using Variance Inflation Factor (VIF) in the R package usdm (Naimi et al., 2014), and standardized the data to zero mean and unit variance in the R package vegan (Oksanen et al., 2018). All variables had VIF < 6 and were retained for further analysis. Finally, the PC1 of Hill-Smith PCA (local environmental variables), MEM3 (spatial arrangement), urbanisation rate, and land surface temperature were entered into a Structural Equation Model to test for their direct and indirect effects on the mean and coefficient of variation of each phenotypic trait.

We calculate the mean as a sum of each phenotypic trait divided by the number of individuals in each pond. Also, we calculated the Kvålseth coefficient of variation (KCV; see Kvålseth, 2017) for each phenotypic trait as the standard deviation divided by the square root of the mean of squared values, with an R function provided by Lobry et al., (2023). Path analysis was performed in the R package lavaan (Rosseel, 2012). Figures were made in the R package semPlot (Epskamp, 2022).

To estimate the adaptive landscape, we separately regressed the Scaled Mass Index as a function of body size, size-free leg length, and head shape. We used a non-parametric regression technique called projection-pursuit regression implemented in the gsg R package (Morrissey, 2014). This approach assumes that the partial regression coefficient represents the selection gradient as a way to estimate the adaptive landscape (Lande & Arnold, 1983). The data and associated R code are available in Dryad (Severgnini & Provete 2024).