Indicators of wood frog (
 Lithobates sylvaticus
 ) condition in a suburbanizing landscape

Wood frogs (Lithobates sylvaticus) are threatened by habitat degradation associated with suburbanization. Suburban development near pools may affect larvae with ramifications for population persistence and vitality. We studied larval development and documented key vernal pool and terrestrial characteristics within 1000 m of 43 pools across the suburban development gradient near Bangor, Maine, USA. Specifically, we examined how survival and morphological characteristics (e.g., developmental phenology, condition, body length, and tail length and shape) varied with characteristics at pool and landscape scales. Secondarily, we explored associations between morphology and survival. Differences in tadpole morphology were associated with suburban land development, hydrology, within-pool vegetation indicative of light availability at the water’s surface, and density of pool-breeding amphibian egg masses. Across all pools, tadpoles had delayed development and were larger in pools with longer hydroperiods, while tadpoles in suburban pools developed earlier and were larger than those in rural pools with comparable hydroperiods. Tadpoles developed later and had longer tails in pools with greater canopy cover. Morphology profiles also differed between rural and suburban sites and among years. Survival in suburban pools was predicted to be 15 times greater than in rural pools, but across all pools (including those at intermediate intensities of suburbanization), survival was not predicted to vary with either morphology or site characteristics. No strong relationship existed between developmental phenology and any condition or size metric. Because rural and suburban tadpoles responded similarly to within-pool conditions, our results support the need to maintain natural hydrology and vegetation conditions in pools even in developing areas. Although we detected benefits to tadpoles with increasing suburbanization, suburbanization is well known to extirpate breeding populations; thus, it is likely that wood frog population declines associated with suburbanization are responding to stressors beyond the pool at terrestrial life stages.


INTRODUCTION
Pool-breeding amphibians are often threatened by habitat loss and degradation as land development associated with suburbanization encroaches near and even eliminates intermittently inundated pools (Zacharias et al. 2007, Windmiller and Calhoun 2008, Baldwin and deMaynadier 2009). More specifically, suburban development has been associated with both extirpation and reduced occupancy of amphibians (Homan et al. 2004, Rubbo and Kiesecker 2005, Clark et al. 2008. Draining or impounding wetlands for development can eliminate breeding pools (Preisser et al. 2000, Beja and Alcazar 2003, Windmiller and Calhoun 2008, and adult road mortality and loss of critical postbreeding habitat via deforestation also likely contribute to population declines (Fahrig et al. 1995, Homan et al. 2004, Eigenbrod et al. 2008. Even when breeding populations persist in urbanizing areas (Le Viol et al. 2012, Scheffers andPaszkowski 2013), factors common to suburban landscapes (e.g., road runoff) can have sublethal influences (Bommarito et al. 2010, Brand et al. 2010) that may not be detectable in demographics until decades later (Griffis-Kyle 2007, Blaustein et al. 2011. Accounting for sublethal effects of suburbanization may be especially important in assessing population condition of species, such as wood frog (Lithobates sylvaticus), that are commonly surveyed by egg mass counts that can have high interannual variability (Berven 1990, Crouch and Paton 2000, Capps et al. 2015. Sublethal responses to suburbanization may also serve as early indicators of demographic declines that would not otherwise be detected for decades (L€ ofvenhaft et al. 2004, Gagn e andFahrig 2010). Measuring multiple responses (e.g., condition, survival) during postembryonic life stages may provide better indicators of the impact of suburbanization. For wood frogs, assessing larval responses may be appropriate as environmental conditions experienced during larval development may influence terrestrial stage morphology and performance (Berven 1990, Boes andBenard 2013). Larval wood frogs are sensitive to pool conditions that can be influenced by nearby vegetation and land use (Skelly et al. 2002, Watkins and Vraspir 2006, Karraker et al. 2008, and their survival may decline with road salt (Sanzo and Hecnar 2006, Karraker et al. 2008, Hall et al. 2017) and herbicide contamination (Relyea 2012) and light pollution (Perry et al. 2008). Few studies have examined the effects of natural and suburban-associated land cover types on amphibian larvae. Some studies have compared larval responses between forest and roadside pools (Karraker et al. 2008, Brady 2013, Hall et al. 2017 and examined the relationship between skeletal abnormalities and road proximity (Reeves et al. 2008). Although Shepack et al. (2017) experimentally introduced larvae into pools where breeding populations had been extirpated by urbanization and examined their survival and morphology, to our knowledge, no study has examined the relative influences of suburban-associated land cover types at multiple spatial scales on the larval survival and development of a vernal pool-breeding amphibian.
Shifts in within-pool vegetation with suburbanization (Ehrenfeld 2000) may also influence wood frog larvae as they respond to vegetation structure (Skelly et al. 2002) and leaf litter composition (Rubbo and Kiesecker 2004). For example, fluctuations in canopy cover and leaf litter composition may influence larval development, growth, and survival (Werner and Glennemeier 1999, Skelly et al. 2002, Halverson et al. 2003, Rubbo and Kiesecker 2004. Changes in vegetation may also alter community interactions involving tadpoles. Submerged vegetation may reduce predation risk by providing cover for larvae (Formanowicz andBobka 1989, Kopp et al. 2006), and reductions in canopy cover may increase the amphibian species richness (Skelly 2014), potentially changing competition and predation pressures.
Understanding the influence of suburban development on wood frog persistence can be particularly challenging because a response to one condition can be mediated by another condition through synergistic or antagonistic effects (Relyea 2004, Marino et al. 2014. For example, tail morphology in wood frogs can respond to aquatic insect predators (Relyea 2012) as well as pool temperature during incubation (Watkins and Vraspir 2006). However, examining multiple metrics can provide a more comprehensive population status assessment than any single measure (Welsh et al. 2008) and can reduce potential error associated with extrapolating single measures, especially from early life stages (e.g., size at metamorphosis) to fitness (Earl and Whiteman 2015).
In this paper, we examined larval wood frog responses to land cover type at the landscape scale (within 1000 m), water quality, hydrology, within-pool vegetation, and indicators of amphibian community competition and predation across a gradient of low-intensity suburbanization (0-38% impervious cover within 1000 m). First, we examined relationships among developmental phenology, condition, body length, tail length, and tail shape within individuals and cohorts (i.e. site-year level) and determined if morphology varied among years and sites or between the least-and most-suburbanized (i.e., rural and suburban) pools. We then identified (1) which morphological and/or site characteristics were important predictors of survival; (2) which site characteristics were important predictors of larval morphology; and (3) if the influence of these predictors varied with suburbanization.

Study area
The greater Bangor, Maine, area is located in the glaciated northeastern United States, a region historically dominated by mixed coniferoushardwood forest (Eakin et al. 2018). Bangor encompasses 90 km 2 with a population of 33,000 (U.S. Census Bureau 2011). The 200-km 2 study area also included Orono, Hampden, and Old Town (populations of~7000-10,000; U.S. Census Bureau 2011). Developed land uses are primarily residential and commercial with development intensity extremes of nearly 100% impervious surface in downtown Bangor (44°48 0 8″ N, 68°46 0 13″ W) to <1% impervious surface (e.g., in the City Forest; Fry et al. 2011).

Site description
Each site consisted of a wood frog breeding pool and the area within 1000 m of the pool's high-water mark. We studied 36, 30, and 13 sites (43 total) in 2014, 2015, and 2016, respectively. We selected sites to represent the range of development intensity at which wood frog breeding occurred using percent impervious cover within 1000 m as an index of development. Selected sites had 0-34% impervious cover within 100 m and 0-38% within 1000 m, 16-82% tree cover within 1000 m, 0-100% above-pool tree canopy density, and 63-9978 m 2 pool area at spring high water.

Measures of morphology
We conducted dip net surveys for wood frog larvae throughout larval development at 25, 30, and 12 pools (39 total) in 2014, 2015, and 2016, respectively. Surveys were conducted once every two weeks 10 June-19 August in 2014, and weekly 15 May-15 September in 2015 and 18 May-13 September in 2016. Up to 20 wood frog larvae were measured per survey. Within each pool, we sampled larvae from multiple locations, representing available vegetation and light conditions, to account for potential larval response to spatial variation within a pool. We recorded Gosner (1960) stage, mass, snout-vent length (SVL), tail length, and tail fin depth (following Relyea 2000).
To account for shifts in morphology with developmental stage and to satisfy assumptions of linearity in the relationship between morphological variables (Green 2001), we calculated individual morphological responses as the residuals of ANCOVA models where Gosner developmental stage was a factorial covariate (Packard andBoardman 1988, Garc ıa-Berthou 2001). We calculated indices of relative developmental phenology (hereafter considered a morphological response, Dev), condition (Cond), body length (Len), and tail length and depth (TailL, TailD; Table 1). Negative developmental phenology values represent delayed phenology for a given developmental stage, and conversely, positive values represent advanced phenology. Condition is used to indicate amount of fats and other energetic reserves relative to body size (reviewed by Green 2001). Because tail fin depth relative to body size (residuals of fin depth regressed against SVL and developmental stage) and fin depth relative to tail length (TailD) were highly correlated (r = 0.83, Pearson's correlation coefficient), we only used TailD as this provided a metric of tail shape and TailL already represented tail size.
For pool-level analyses, we extracted aggregate morphology metrics: the 10%, 50%, and 90% quantiles and standard deviation (SD) from each site-year. Quantiles and SD better represent the variability within a population compared to the mean response (Cade and Noon 2003) and can be used to address the possibility that different segments within one population may be limited by different factors. We calculated aggregate measures for site-years with ≥20 tadpoles total collected during ≥2 site visits in a year. All animal handling was conducted under the University of Maine Institutional Animal Care and Use Committee Protocols A2014-05-04 and A2015-03-03.

Egg mass surveys
We conducted egg mass surveys at 43 pools to provide an index of the initial abundance of individuals in a pool as a baseline for estimating survival through metamorphosis. We also used the density of egg masses (number 9 m À2 pool area) to indicate conspecific competition (wood frog eggs) and predatory pressure from larval salamanders (spotted salamander and blue-spotted salamander eggs, which included the unisexual complex, Ambystoma laterale-jeffersonianum). We counted these pool-breeding amphibian egg masses following the apparent peak of breeding using a dependent double-observer method to increase detection (Grant et al. 2005). If eggs had been recently laid (within approximately two days), we revisited pools and counted new eggs.

Estimation of survival
We used repeated counts of late stage (Gosner stage ≥41) tadpoles to calculate an index of survival (hereto referred to as survival). We conducted tadpole counts using dip net surveys at 42 sites (35, 26, and 12 pools in 2014, 2015, and 2016). The number of net sweeps was based on pool size at time of sampling. And surveys began at a site when we detected a tadpole developed to Gosner stage ≥41 (approximate stage indicating pro-metamorphosis and thus a useful proxy of survival through metamorphosis) during morphology observations. We conducted 3 surveys within 3-4 d at each pool, and tadpoles were sampled without replacement. Additional details of survey methods including calculation of the number of sweeps are in Appendix S1.
We used package unmarked (Fiske et al. 2017) to fit negative-binomial N-mixture models (MacKenzie et al. 2006) that estimated abundance (k) of tadpoles 9 m À2 within sampled areas by site while accounting for imperfect detection (p; Royle 2004). We extrapolated the estimated tadpoles 9 m À2 to the entire sampled pool area to estimate the abundance of tadpoles in a pool. We then used the residuals of total abundance regressed against the number of egg masses to estimate survival by pool. We recognize that this method does not incorporate error in abundance estimates and egg mass counts. However, markrecapture methods, which would theoretically provide more robust estimates of survival, are not feasible because embryos cannot be marked and marking techniques are unreliable for tadpoles (Grant 2008, Carlson andLangkilde 2013) or require a tissue sample at each capture for genetic analysis (Ringler et al. 2015), potentially affecting survival. Further details of abundance estimation are in Appendix S1.

Site characteristics
We used ArcView GIS10.2 and the Maine Land Cover Dataset (2004 all land use; 2011 impervious surface) to quantify the percent impervious and tree cover within 100, 300, 600, and 1000 m from pool high-water marks. We used impervious cover within 300 m (based on scales of land cover type identified as important during random forest analyses, RFA) as an indicator of suburbanization intensity. We assigned the 10 pools with the most impervious cover at 300 m (16.5-36.8%) to a suburban group and the 10 pools with the least impervious cover at 300 m (0-1.9%) to a rural group.
We measured pool area at spring high water using a sub-decimeter accuracy GPS (Trimble Notes: Missing data for some variables for some individuals and site-years account for differences in the number of siteyear, site, and individual observations. Dev, residuals of the inverse of Julian day of measurement regressed on developmental stage; Cond, residuals of square-root-transformed mass (g) regressed on SVL and developmental stage; Len, residual of body length regressed on developmental stage; TailD, residuals of tail depth regressed against tail length and developmental stage; TailL, residual of tail length regressed on SVL and developmental stage.
† Responses only calculated for site-year. ‡ 10, 50, 90, and standard deviation indicate 10, 50, and 90% quantiles and standard deviation of responses within a site-year.
❖ www.esajournals.org Geo 7x, Westminster, Colorado, USA) and maximum spring-high-water depth (17 April-22 May) using a pole marked in centimeter increments. Hydroperiod was determined by the Julian day that standing water was no longer present as observed during field visits or by trail cameras placed at pools for a separate study. Because tadpoles can survive in small volumes of water (<1 L) for a short time, even small puddles were used to indicate standing water. To treat hydroperiod as a continuous variable in subsequent analyses, pools that did not dry in a given year were assigned Julian day 280 if they dried to within 10 cm or Julian day 300 if they had ≥10 cm of water at the deepest point (values only determine rank order of continuous explanatory variables used in RFA; De'ath and Fabricius 2000). These values were selected because no pool dried during the study after Julian day equivalent 266 (23 September).
We measured vegetation in spring and summer within-pool basins at 25, 36, and 11 sites in 2014, 2015, and 2016, respectively. Spring surveys were conducted 2-9 June 2014, 2-20 May 2015, and 19 May-7 June 2016 when vernal poolbreeding amphibian eggs were present. Summer surveys were conducted 24 July-21 August 2014, 22 July-17 August 2015, and 19 July-26 July 2016 at late summer dry down. We visually estimated percent cover of shrub, emergent vegetation, and submerged vegetation to the nearest 10%. We separately estimated duckweed (subfamily Lemnoideae) because we suspected it may be related to anthropogenic disturbance. Duckweed was estimated to the nearest 10%, using 1% to indicate any presence. Woody vegetation canopy density over pools was measured~1 m above the ground near the pool center using a spherical convex densitometer.
We used water probes (Hach, Loveland, Colorado, USA) at 43 pools to sample specific conductance, dissolved oxygen (DO), and water temperature, more specifically 34 pools (31 pools once and 3 pools twice) 3 May-10 June 2014; 37 pools (2 pools once and 35 pools twice) 2 May-12 June 2015; and 12 pools (4 pools once, 6 pools twice, and 2 pools thrice) 6 May-16 June 2016. On each date a pool was sampled, we collected and tested 1 L of surface water~1 m from the water edge at each of three equidistant points around the perimeter. Only one sample was taken at pools that were almost dry and <2 m 2 . All testing was conducted at the pool edge within minutes of sample collection. Each metric was averaged by day and then year to calculate values used in analyses.
Some missing values occurred in site characteristic data: ≤5% in datasets used for analysis of tadpole responses. We used k-nearest neighbor (kNN) imputation (package DMwR; Torgo 2010) to replace missing values with a weighted mean of the 10 nearest neighbors. Single imputation using kNN provides robust missing value estimates for datasets with ≤20% missing values (Troyanskaya et al. 2001).

Morphology profile analysis
We conducted principal components analyses (PCAs) using the vegan package (Oksanen et al. 2017) in R version 3.3.1 (R Development Core Team 2016) to identify dominant gradients of variation in tadpole morphology and development within site-years and individuals (Borcard et al. 2011). To examine site-year profiles, we created a correlation matrix of 10%, 50%, and 90% quantiles and SD of developmental phenology, condition, body length, and tail length and depth for site-years with complete observations across these variables (64 site-years; 38 sites). To examine individual morphology profiles, we created a correlation matrix of developmental phenology, condition, body length, tail length, and tail shape variables for individuals with complete observations from site-years used in the site-year PCA. Variables were standardized to have a mean of zero and standard deviation of one to avoid unequal weighting among variables within the individual and site-year datasets, respectively. We used a Monte Carlo randomization test (1000 permutations) to select significant axes (P < 0.01) for further examination.
We compared the individual and site-year ordinations to examine the possibility that the tadpoles in a pool may display multiple, distinct profile types (i.e., multiple morphology strategies within one pool). We assumed that the individual ordination was a fair representation of the relationships among morphology characteristics, even at the site-year level, as the relative eigenvector contributions of site-year variables to the correlation matrix were similar to those of the same variable type in the individual ❖ www.esajournals.org ordination (e.g., the relationships of Cond10, Cond50, and Cond90 relative to other site-year variables were similar to the relationship of Cond relative to other individual variables). Because individual and site-year profiles were similar, we used individual profiles to examine differences in morphological response (site and year effects and between rural and suburban groups). Because survival may not depend on all tadpoles in a pool (e.g., the smallest or slowdeveloping) and/or may be influenced by within-pool variance in responses, we used siteyear profiles to examine the relative importance of aggregate responses (including within siteyear variance) to predict survival.
To test for Site, Year, and Site 9 Year effects and suburbanization level (Suburban), Year, and Year 9 Suburban effects in individual profiles, we conducted permutations of multiple analysis of variance (PERMANOVA, package vegan) for a site-year model and a suburban-year model. The suburban-year model was conducted for the subset of individuals included in rural and suburban pools. We ran 5000 permutations of PERMA-NOVA using a Bray-Curtis similarity matrix (Anderson and Walsh 2013) and nested individuals by site-year for permutations (Anderson 2001). We used the F-value as a signal-to-noise ratio to indicate effect size (McCune et al. 1997, Short andMorris 2016). If there were significant differences in multivariate profiles identified in site-year or suburban-year models, we conducted separate ANOVAs to examine morphology response among sites and years (Site, Year, and Site 9 Year effects) and among suburbanization levels and years (Suburban, Year, and Suburban 9 Year as fixed effects; nested by Site; package lme4; Bates et al. 2017). We created a set of five site-year models and a set of five suburban-year models, each including all possible combinations of fixed effects including a null (intercept only) model and used Akaike's information criterion (AIC) to select the top model from each set. We considered models DAIC < 2 that ranked above the null model to be plausible (Burnham and Anderson 2002).
For the site-year dataset, we used classification and regression trees analyses, specifically RFA, to identify the relative importance of aggregate tadpole morphology variables in explaining survival. We conducted RFA using package randomForest (Liaw and Wiener 2002) in R version 3.3.1 (R Core Team 2016). We bootstrapped with replacement to build 10,000 regression trees (random forest error stabilized at~1000-2000 trees for each response variable), using 2/3 of the data at each iteration. We calculated explanatory variable importance using the mean percent decrease in accuracy resulting from removal of each variable to rank the importance of explanatory variables. We used package randomForestSRC (Ishwaran and Kogalur 2014) to create partial dependence plots (PDPs) that examine the marginal effects of predictor variables while holding all other predictors at average values (Friedman 2001, Cutler et al. 2007). Because PDPs display general trends, all reported values are approximate.

Modeling tadpole response on site characteristics
We conducted ANCOVAs to examine if suburbanization level (rural, suburban) or year changed the relationship between number of egg masses and estimated abundance of tadpoles that survived to metamorphosis (i.e., survival) between rural and suburban pools. We used etasquared (g 2 ) to indicate effect size (i.e., % change in the response accounted for by the predictor; Levine and Hullett 2002). Survival was examined for rural and suburban pools with <70 egg masses; 4 suburban site-years with 100-402 egg masses were dropped from analyses to maintain comparable x-axis values between rural and suburban sites.
We used RFA following the same methods described above to identify the importance of site characteristics ( Table 2) that may affect survival and aggregate measures of developmental phenology, condition, and tail length and shape (Table 1). We included the respective sample size (DEV.N, COND.N, LEN.N, TAILL.N, TAILS.N) as a predictor in each model because scatterplot trends suggested that variance decreased (condition, body length, and tail length) or increased (developmental phenology) with sample size. A decrease may be an artifact of larger sample size leading to increasingly precise estimates of group means, or a decrease or increase may be a result of site characteristics that affect the availability of tadpoles for sampling throughout the season (therefore the sample size) also affecting variance. Because we were interested in the influence of site conditions at the pool scale, we modeled annual (site-year) responses.
For each survival or aggregate response with variance explained by RFA, we conducted a separate ANCOVA for each of the best predictors identified during RFA (those for which we present partial dependence plots). Global model terms included suburbanization level and an RFA predictor. We dropped non-significant interaction terms (F test; P > 0.05) and kept all main effects contained in significant interactions. We disregarded models with statistically significant RFA predictors that did not have a similar range in values (x-axis distribution) across suburbanization levels.

Morphological profile
Individual morphology (Dev, Cond, Len, TailL, TailD) was influenced by the combined effect of site and year in site-year models and by the combined effect of suburbanization level and year in suburban-year models. In both morphology profile models as well as in models of single morphology responses, all predictors (main effects and interactions) were significant (PER-MANOVA, P < 0.001; full models were the top model for each morphology response; Appendix S2). For all models, the combined effects of the main effects (interactions) prevented the interpretation of main effects across all factor levels (e.g., when considering 95% confidence intervals, condition was higher in suburban pools in 2014 but not necessarily in 2015 and 2016). Developmental rate was not strongly associated with condition, body length, or tail morphology, but tail length and condition were positively associated, as were body length and tail depth (Fig. 1). No distinct profile types (i.e., clusters) emerged in morphology profiles. Tree cover within 100 m of pool (%) TREE1000 Tree cover within 1000 m of pool (%) TREE300 Tree cover within 300 m of pool (%) TREE600 Tree cover within 600 m of pool (%) Notes: All variables except for impervious and tree cover were measured by site-year; impervious and tree cover were measured by site.
† Julian day of drying adjusted for pools that did not dry (details in Methods).

Survival
Survival in suburban pools was predicted to be~15 times greater (1530% greater, b = 2.73, standard error = 1.12) than in rural pools (F 1,20 = 5.93, P = 0.02, Adj. R 2 = 0.36; Fig. 2). Random forest analyses survival models using morphology or site characteristic predictors did not explain any variation in the response (i.e., not different from random); thus, predictor importance was not examined.

Developmental phenology
Tadpole phenology was predicted to be earlier in warmer pools with a shorter duration and less tree cover. Later phenology was associated with longer hydroperiods, especially for pools drying after Julian day 200, 19 July, and in pools with >60% canopy cover (Fig. 3a, b and Appendix S4: Fig. S1a-d; see Appendix S3 for ranking of predictors). Conversely, earlier developmental phenology was predicted by warmer water (Appendix S4: Fig. S1e) with the greatest increases in the lower third of the range of measured temperatures. Finally, greater tree cover was associated with later phenology, with the greatest decrease at >60% tree cover within 100 m of pools (Fig. 3c).

Condition
Tadpoles in pools with characteristics associated with a longer hydroperiod and deeper water, including less understory vegetation, were predicted to be in better condition. Emergent vegetation and shrub cover were important predictors of condition, as was greater depth and duration of water (positively associated with depths >0.4 m and longer hydroperiods; Figure 3d-g; Appendix S4: Fig. S1g-i). Impervious cover and sample size were also high-ranking predictors positively associated with condition (Appendix S4: Fig. S1f, h, and j).

Body length
Site characteristics indicative of more water (depth, lower conspecific density, longer duration of water) were important predictors of greater body length. Pool depth and hydroperiod Fig. 1. Principal components analysis (PCA) of individual wood frog tadpole morphology profiles. The first and second axes (PC1 and PC2) were significant (P < 0.001) and account for 34.7 and 27.6% of the total variation, respectively. The first PC was positively described by relatively equal contributions from tadpole condition and tail length (0.64 and 0.62 respective eigenvalues). The second PC was negatively described by relatively equal contributions from body size and tail shape (À0.69 and À0.63 eigenvalues). Developmental phenology is poorly represented (communality value = 0.14 compared to 0.70-0.79 for other metrics). were important predictors of body length, with deeper pools and pools drying after Julian day 190, 9 July, associated with longer bodies (Figure 3i and k; Appendix S4: Fig. S1l). Greater density of wood frog eggs (up to 0.2 egg masses 9 m À2 ) and sample size (up to 100 individuals) were associated with shorter bodies (Fig. 4h, j; Appendix S4: Figure S1n). Although canopy was a top-ranking predictor, lowess smoothing of predicted trends suggests body length is relatively invariant to canopy cover (Appendix S4: Figure S1k, m).

Tail length
Site characteristics related to low-light conditions or greater suburbanization were important predictors of greater tail length. The greatest increases in tail length were where summer canopy cover >90% ( Fig. 3l; Appendix S4: Figure S1o,q), at sites with 0-10% impervious cover within 300 m (Figure 3o; Appendix 4Ss), and at sites with < 40% spring emergent vegetation and/or <30% summer shrub cover ( Fig. 3m; Appendix S4: Figure S1p). Greater tree cover at 300 and 1000 m was associated with shorter tails (Figure 3n; Appendix S4: Figure S1r).

Tail shape
Pool area was by far the most important predictor for tail shape (pools <50 m 2 predicted to have tadpoles with broader fins; Appendix S4: Figure S1t, u). However, as tail shape was predicted to be relatively invariant to pool area, we hesitate to interpret the direction of association between area and tail shape.

Morphological variation within pools
Site characteristics indicative of more water (duration or depth) and higher light conditions (less canopy, greater submerged vegetation) were important predictors of variance in developmental phenology (positively associated at 0-20% spring submerged vegetation and >0.5 m maximum depth; negatively associated with canopy cover >60%; Fig. 3p-r), condition and body length (both positively associated with pool drying between Julian day 200 and day 230; 19 July-18 August; Fig. 3s; Appendix S4: Figure S1z), and tail length (negatively associated with canopy cover >80%; Appendix S4: Figure S1a1, a2). Amphibian egg mass densities were also important predictors of variance; variance in developmental phenology was positively associated with 0-3 blue-spotted salamander eggs 9 m À2 , and variance in tail length was negatively associated with wood frog egg density (Appendix S4: Fig. S1w and a3). Sample size-which was likely responsive to duration of larval presence, abundance, and within-pool detection probability -was also an important predictor of variation in developmental phenology (roughly linearly, positively correlated), condition (negatively associated with increased sample size up to 80 individuals), and body length (negatively associated with increased sample size up to 75 individuals; Appendix S4: Fig. S1v, x, and y). Impervious cover was an important predictor of variance in tail length (positively associated with increases in impervious cover within 300 m between 0% and 10%; Appendix S4: Figure S1a3).

Differences in response to site characteristics between suburbanization levels
Developmental phenology, condition, and body length were greater in suburban than rural pools when accounting for the influence of important site characteristics ( Fig. 4; Appendix S5). The effect of these site characteristics was the same across suburbanization levels (P > 0.05 for interaction terms).

DISCUSSION
This study provides evidence that human land disturbance surrounding pools (within 1000 m) can influence tadpole survival and morphology, an issue that can influence subsequent terrestrial stages and potentially population vitality (Berven 1990, Boes andBenard 2013). Moreover, we compared the relative importance of land cover and site conditions within vernal pools and demonstrated that conditions both within and surrounding pools are likely to influence larval survival and morphology. The importance of site characteristics on multiple spatial scales emphasizes the complexity of wood frog larval response and the need to consider multiple variables when assessing population response to suburbanization. Tadpole responses primarily varied with variables related to suburbanization, pool hydrology, and within-pool vegetation, each covered in the next three sections. Fig. 3. Partial dependence plots from random forest analysis predictions of relative developmental phenology for larval wood frog. Responses (y-axis) are 10%, 50%, and 90% quantiles of developmental phenology (Dev; ac), condition (Cond; d-g), body length (Len; h-k), tail length (TailL; l-o), and tail depth (TailD; p-s) and standard deviation (SD; v-a3) for each morphology characteristic. Predictors are hydroperiod indicated by Julian day of pool drying (HYDRO); impervious cover within 300 m of pool (%; IMP300); L. sylvaticus egg masses 9 m À2 (LSDEN); maximum pool depth at spring high water (m; MAXDEPTH); within-pool spring tree canopy, emergent vegetation, shrub, and submerged vegetation cover (%; SPCAN, SPEMERG, SPSHRUB, SPSUBMG);

Urbanization
Suburbanization was positively associated with earlier developmental phenology (both individuals and cohorts) and greater cohort survival. Additionally, changes in land cover type consistent with suburbanization (less tree cover and greater impervious cover) were associated with advanced developmental phenology and increased condition, tail length, and greater variation in tail length. These survival and condition responses seem counterintuitive given (1) documented negative impacts of similar intensities of suburbanization on adult wood frogs (Homan et al. 2004, Rubbo and Kiesecker 2005, Skidds et al. 2007, Clark et al. 2008) and (2) evidence that wood frog eggs and larvae are harmed by elevated road salt concentrations (Brady 2013). However, larval wood frogs have been observed to be larger at metamorphosis in stormwater wetlands (Scheffers and Paszkowski 2016) and have greater survival through metamorphosis after the forest surrounding pools was cut ). Furthermore, some other amphibian species, including common green frog (Pelophylax perezi), southern leopard frog (Rana utricularia), pig frog (R. grylio), American bullfrog (R. catesbeiana), and eastern narrowmouth toad (Gastrophryne carolinensis), have also shown positive morphological and survival responses in suburban areas Paszkowski 2011, Iglesias-Carrasco et al. 2017) and in areas with less terrestrial habitat (Salice et al. 2011). Earlier developmental phenology paired with greater survival in areas of greater suburbanization suggests that advanced phenology is advantageous at these pools. Thus, it is more likely that advanced development could be prompted by the well-demonstrated effects of low conspecific density and high food availability for anurans, including wood frog (Licht 1967, Berven and Chadra 1988, Murray 1990; but see Wilbur 1977) rather than by accelerated pool drying, which results in tadpoles completing metamorphosis earlier and at a smaller size (Crump 1989, Dimauro and Hunter 2002, Loman and Claesson 2003. Completing metamorphosis early may also confer further benefits to froglets by providing more time for froglets to develop metabolic reserves before hibernation. Longer tails and greater variation in tail length in suburban areas could indicate relatively large and variable predator and competition pressures across a gradient of increasing suburbanization. Longer tails may provide a competitive advantage, potentially by increasing the force with which larvae scrape food (i.e., periphyton), thus increasing foraging efficiency (Relyea 2000). Longer and broader tails are associated with greater survival in the presence of invertebrate predators (Van Buskirk and Relyea 1998) and can distract invertebrate predator attacks away from the body and head of a tadpole, reducing the lethality of strikes (Van Buskrik et al. 2003). Tadpoles with this predator-induced tail shape may be less vulnerable to predation but suffer lower growth rates and smaller bodies (Van Buskirk and Relyea 1998, Relyea 2004. Our finding that larval wood frogs respond to differences in conditions at the landscape scale differs from that of Shepack et al. (2017). However, this difference may be because our study was observational and thus did not control for sources of variation (genetic, predatory pressures, or resource quality and availability within pools) that could be influenced by land cover type and likely affects tadpoles.
We also demonstrate that within a landscape with relatively low-intensity suburbanization (0-38% impervious surface within 1000 m), larvae are likely influenced by land cover types indicative of suburbanization at multiple spatial scales. For example, tadpole development rate was sensitive to reductions in tree cover nearer to pools, with a steep change in rates predicted with initial reductions in tree cover (up to 30%) within within-pool summer tree canopy and shrub cover (%; SUCAN, SUSHRUB); tree cover within 100 m and 300 m of pool (%; TREE100, TREE300); and sample size for body length measurements (LEN.N). Dashed lines correspond to lowess-smoothed lines representing the partial dependence between an explanatory variable and response. The solid lines indicate a smoothed error bar of AE two standard errors. The dots indicate the partial values used to fit the lowess function. (Fig. 3. Continued) ❖ www.esajournals.org 100 m of pools, whereas greater tree cover within 1000 m was correlated with shorter tails (Fig. 3; Appendix S4). Additionally, even slight increases in impervious cover (0-10%) within 300 m were associated with shifts in body condition. Other studies have found adult wood frogs to respond to land cover types at similar scales Kiesecker 2005, Skidds et al. 2007) and with similar sensitivity (e.g., 88% habitat cover threshold within 30 m of pools, 44% cover threshold at 1000 m; Homan et al. 2004); thus, our findings suggest that larvae respond to land cover at scales similar to those affecting adults.
Because survival differed between suburbanization levels but the RFA using site characteristics as predictors did not explain any variance in survival, some unmeasured difference between rural and suburban pools may be responsible for differences in survival. Moreover, it is possible that positive tadpole responses to suburbanization may be driven by local adaptation to aquatic and/or terrestrial environments, as wood frog can rapidly adapt to local conditions (pool and landscape) on the order of decades (Relyea 2002, Skelly 2004. Unmeasured factors may contribute to greater tadpole survival as well as size in suburban areas. For example, we did not quantify invasive vegetation, but glossy buckthorn (Rhamnus frangula), an invasive plant in our study area, is associated with larger and faster developing wood frog tadpoles with higher survival (Stephens et al. 2013). Pathogens could also explain greater survival and larger morphology in suburban pools given that Batrachochytrium dendrobatidis (Bd) and ranavirus (Rv) have higher rates of infection and occurrence in rural pools than suburban pools (Crespi et al. 2015, Saenz et al. 2015. Additionally, advanced developmental phenology increases survival to emergence in the presence of Rv (Reeve et al. 2013), and this may explain greater survival even in infected suburban populations. Suburbanization may also alter aquatic insect communities and decrease predation pressure on wood frog larvae. For example, insecticides can reduce insect predator abundance, thereby increasing tadpole survival (Relyea 2005) and potentially increasing feeding (Petranka and Hayes 2011, Cothran et al. 2013).
Greater larval survival and individual size may translate to advantages during terrestrial life stages in suburban areas (Scheffers and Fig. 4. Associations between wood frog tadpole cohort responses and important site characteristics for which there were significant differences between rural and suburban pools (Dev50: F urban:1,31 = 15.57, P < 0. 001; Cond90: F urban:1,29 = 4.79, P = 0.006; Len90: F urban:1,28 = 4.81, P = 0.004). Responses are the 50% quantile of relative developmental phenology (Dev50); the 90% quantile of relative condition (Cond90); and the 90% quantile of relative body length (Len90). Site characteristics are hydroperiod indicated by Julian day of pool drying (HYDRO) and sample size for measurements of relative body length (LEN.N). Points represent pool-year responses. Shaded areas represent parameter 95% confidence intervals.
Paszkowski 2016), as has been noted in rural areas (Berven 1990). However, because juvenile wood frog growth can respond to terrestrial conditions (Boone 2005, Dananay et al. 2015, larval benefits from suburbanization may be overshadowed by degradation of terrestrial habitat. Removal of tree cover, which protects from desiccation, can increase adult mortality (Rittenhouse et al. 2009). Other novel threats, such as roads, sewer grates, predatory pets, and lawn mowers, may reduce abundance and survival of terrestrial wood frogs (Eigenbrod et al. 2008. Alternatively, larval benefits from suburbanization may persist beyond metamorphosis and lengthen the time breeding populations persist in suburbanized areas and thus contribute to the multi-decadal time lag noted for amphibians in urbanizing landscapes (L€ ofvenhaft et al. 2004, Metzger et al. 2009, Gagn e and Fahrig 2010. Even if suburbanization benefits larvae and adequate terrestrial habitat remains near pools, reduced habitat connectivity associated with human land use conversion can increase the risk of wood frog extirpation (Harper et al. 2008).

Hydrology
Tadpole responses to hydroperiod and pool depth across a suburbanization gradient were similar to those found by others in more natural landscapes. More specifically, our results of delayed developmental phenology and greater body length and condition with longer hydroperiod align with other studies noting positive correlations between hydroperiod and wood frog larvae size (Karraker and Gibbs 2009) and larval period (Dimauro and Hunter 2002, Skelly 2004, Gervasi and Foufopoulos 2008. Although developmental phenology and condition in suburban and rural pools had similar rates of change with increases in hydroperiod (Figure 4; Appendix S5), suburban pools were predicted to develop earlier (all quantiles), have greater condition (90% quantiles), and show more variance in condition than rural pools. This suggests that although hydroperiod is an important predictor of developmental phenology and condition, other factors related to impervious cover within 300 m (our metric for delineating rural and suburban categories) influence these responses. For example, predator pressure, which leads to slower development (Relyea 2002), may decrease as suburbanization levels increase.
Additionally, we found better condition and greater length associated with deeper pools perhaps because deeper pools have more volume and lower conspecific density, and thus competition, relative to pool area (Brooks and Hayashi 2002). Pool depth also may correlate with increased condition and body length because deeper pools are likely cooler at the bottom and cool water is linked to greater body size in wood frog tadpoles (Berven 1982). The influence of pool depth on condition and length may affect adult stages because larger tadpoles produce larger adults (Werner 1986, Berven 1990) and larger females produce more eggs (Berven 1988). This might explain the increase in reproductive effort with pool depth that has been observed by others (Dimauro andHunter 2002, Calhoun et al. 2003).

Vegetation
Across the suburbanization gradient, withinpool vegetation conditions indicative of low-light conditions at a pool's surface were associated with greater condition, and tail length and slower development (i.e., positive associations with canopy cover and negative associations with emergent vegetation or shrub cover). Our results align with Boes and Benard (2013), who observed larger wood frog metamorphs in closed canopy pools, and Skelly (2004), who observed later development in closed canopy pools. Lowlight conditions could also facilitate increased condition and size by maintaining lower water temperatures which can increase size of wood frog larvae (Berven 1982).

Morphology
We did not detect strong relationships among larval developmental phenology, length, and condition ( Fig. 1) or between morphology and survival. Although we did not necessarily expect consistent (i.e., strong and directional) relationships among morphology metrics (Berven 1987, Berven andChadra 1988), we did expect survival to increase with measures indicating size (Berven 1990). The lack of relationship between survival and morphology metrics suggests that unmeasured factors may determine survival or that a wide range of morphological responses may result in similar survival depending on the local context. Although larval survival was not well predicted by morphology, other size metrics may also be consequential for populations (i.e., breeding population size is limited by terrestrial density-dependent factors ;Berven 1990). Given that larval body size is usually correlated with postmetamorphic body size (Boes andBenard 2013, Berven 1990; but see Earl and Semlitsch 2013), tadpole size may translate into fitness advantages associated with greater body size (Berven 1988).

CONCLUSION
Conserving amphibians in developing landscapes is challenging given the complexities of responses to suburbanization, which can introduce a multitude of novel factors that may influence amphibians in both aquatic and terrestrial stages and habitats (Windmiller and Calhoun 2008). Our results indicate that larvae respond to conditions within 1000 m of pools as well as within-pool site conditions. Vegetation conditions consistent with greater canopy cover benefit larval development, and canopy cover is greater in rural than suburban pools. Notably, other vegetation metrics that correspond to canopy density (e.g., lesser emergent and shrub cover) as well as pool depth and duration, which also benefit larval development, do not differ between rural and suburban pools. Additionally, the site characteristics most closely associated with suburbanization (e.g., impervious cover, road salt contamination) were not necessarily the most important predictors of tadpole responses. This suggests that differences in response may be influenced by the confluence of changes across many site characteristics resulting from suburbanization. Although we found responses in suburban pools were more beneficial than in rural pools (e.g., greater survival), it is well documented that ultimately suburbanization extirpates breeding populations. Given the suburban-associated benefits to larval stages, it is unlikely that wood frog population declines stem from larval responses to suburbanization. We suspect that any suburban-associated benefit to larval stages is outweighed by negative impacts to terrestrial life stages, and thus, we suggest future study focus on this issue. We also recommend examining larval and terrestrial stage responses to greater intensities of suburbanization than were represented in our study area because response strength may increase with suburbanization.