Detecting spatial ontogenetic niche shifts in complex dendritic ecological networks

Ontogenetic niche shifts (ONS) are important drivers of population and community dynamics, but they can be difficult to identify for species with prolonged larval or juvenile stages, or for species that inhabit continuous habitats. Most studies of ONS focus on single transitions among discrete habitat patches at local scales. However, for species with long larval or juvenile periods, affinity for particular locations within connected habitat networks may differ among cohorts. The resulting spatial patterns of distribution can result from a combination of landscape-scale habitat structure, position of a habitat patch within a network, and local habitat characteristics—all of which may interact and change as individuals grow. We estimated such spatial ONS for spring salamanders (Gyrinophilus porphyriticus), which have a larval period that can last 4 years or more. Using mixture models to identify larval cohorts from size frequency data, we fit occupancy models for each age class using two measures of the branching structure of stream networks and three measures of stream network position. Larval salamander cohorts showed different preferences for the position of a site within the stream network, and the strength of these responses depended on the basin-wide spatial structure of the stream network. The isolation of a site had a stronger effect on occupancy in watersheds with more isolated headwater streams, while the catchment area, which is associated with gradients in stream habitat, had a stronger effect on occupancy in watersheds with more paired headwater streams. Our results show that considering the spatial structure of habitat networks can provide new insights on ONS in long-lived species.


INTRODUCTION
Life history responses to environmental variation include changes in niche position as individuals mature, and these ontogenetic niche shifts (ONS) contribute to regulation of population growth and community structure (Werner and Gilliam 1984, Claessen and Dieckmann 2002, Nakazawa 2014, van Leeuwen et al. 2014, especially for species with complex life cycles or delayed maturity. For these species, consequences of habitat choice in early life stages may influence population growth rates by allowing access to differentially available resources (e.g., Mittelbach et al. 1988) or via differential sensitivity of life stages (e.g., Crouse et al. 1987, Rudolf and Lafferty 2011, Regnier et al. 2012, Burton and Metcalfe 2014. Resource levels within different habitats may depend on the position of sites within the landscape (Turner et al. 1997, Laurance 2008, complicating the assessment of ONS at small spatial scales. Most studies of ONS focus on changes in diet or habitat use during a single transition within the lifespan of an organism, such as changes associated with maturation or achieving a size threshold (Graham et al. 2007, Ayll on et al. 2010). Some studies have included variation in local resources among sites, but they generally do not describe how different cohorts vary in their response to broader patterns of habitat structure and isolation within a landscape. Specifically, although studies often compare resource levels associated with ONS, how landscape structure mediates variation in resources is poorly understood. A new focus on spatial ONS, which accounts for the position of a habitat patch within landscapes as well as the geometric configuration of landscapes, can provide new insights into the mechanisms affecting population distribution and persistence.
The geometry of stream ecosystems drives resource distributions (Benda et al. 2004), results in predictable gradients of habitat isolation, and constrains movement for aquatic-obligate species, making these good systems for studying spatial ONS. Unlike terrestrial environments, streams are formed by geomorphic processes that result in predictable gradients in habitat and resource distributions (Vannote et al. 1980, Rodriguez-Iturbe and Rinaldo 2001, Rosi-Marshall and Wallace 2002, Raymond et al. 2016. For example, resource distributions may change at confluences, and species that inhabit small streams may respond to this spatial variation in resources by moving to areas with higher resource levels (Uno and Power 2015). The ability to take advantage of these spatially distributed resources depends on the degree to which species are restricted to moving within or along stream channels (Vannote et al. 1980, Heino et al. 2005. Further, the overall structure of the landscape changes the relative importance of isolation (Andren 1994, Carrara et al. 2012. As individuals age and require different resources for survival, growth, or reproduction, their ability to access complementary habitats can be strongly influenced by the stream network geometry. Therefore, considering both the structure of stream networks and the isolation of habitat within networks is important to describe spatial ONS. We test for spatial ONS by characterizing shifts in distributions of different cohorts of aquaticobligate larval spring salamanders (Gyrinophilus porphyriticus) in headwater stream networks. Spring salamanders have a prolonged larval period that can last 4 years or more (Bruce andCastanet 2006, Resetarits 1995), have continuous growth, and multiple age classes may be present in a given stream segment. We identify age cohorts based on size frequency data of captured individuals, and use the occurrence patterns of each age class to test two predictions regarding spatial ONS in spring salamanders. First, we predict that occupancy patterns will differ among cohorts and be related to measures of network geometry. Second, we predict that any response to stream network position will be mediated by the overall structure of the network. The role of space in population dynamics has long been appreciated; the study of ONS can also benefit from new approaches that account for the isolation of resources in landscapes.

Field surveys
To detect spatial ONS, we sampled 59 stream segments across 10 randomly selected watersheds within the Shenandoah National Park, Virginia, United States. The park ranges in elevation from 400 to 1200 m above sea level along a north-south ridge that runs over 150 km, and is primarily covered by hardwood forest. We selected five watersheds with outlets at the park boundary that were drained by at least a second-order stream. We also sampled an adjacent watershed on the opposite side of the main ridge that divides the park to control for any effect of aspect on Gyrinophilus porphyriticus distribution.
Stream segments were sampled for spring salamander larvae during June and July 2012. We sampled each 20-m stream segment with three consecutive passes, and individual larvae were placed into a plastic bag at the time of their capture. Stream segments were located 10 m and 150 m upstream of confluences on the main channel and tributary, and they contained a range of stream habitats (e.g., riffles, runs, pools). Larval salamanders were measured to the nearest millimeter (snout-vent and total lengths, Heyer 1994) before being released near their point of capture.
We quantified in-stream habitat at the time of surveys by measuring stream width and stream depth at three random locations along the sampling transect. We also estimated percent of each transect covered by silt, sand, pebble, gravel, cobble, and boulder in the streambed (Hauer and Lamberti 2007). We log-transformed stream width and stream depth and standardized variables prior to conducting a principal components analysis to describe major gradients in stream habitat.

Characterizing stream network topology
To characterize spatial ONS of larval salamanders, we calculated three metrics that describe the geometry of stream networks (i.e., network topology). First, the catchment area of a stream segment (the total area drained by a stream segment) is correlated with stream size and associated gradients in stream habitat (Vannote et al. 1980). We also derived two centrality metrics that characterize the position of stream segments within a network relative to potential dispersal behavior of animals (Estrada and Bodin 2008). First, betweenness centrality reflects the degree to which a location acts as a central connection for pathways through a network. Betweenness centrality quantifies the proportion of pathways between all the sites in a network that pass through a given site. We hypothesized that betweenness could be related to spatial ONS because downstream movements would put individuals in areas with higher betweenness centrality, which may be advantageous for exploration of the network (Estrada andBodin 2008, Altermatt 2013). Second, farness centrality reflects the average distance from one stream segment to all other stream segments within the network (Altermatt 2013). This metric may characterize spatial ONS because downstream movements out of first-order streams would see individuals moving into stream segments with lower farness centrality. Farness centrality can be calculated with weighted distances to allow for upstream or downstream bias in movement, and we calculated it by weighting downstream movements twice as costly as upstream movements based on a prior study of dispersal behavior of this species (Lowe 2002). Both centrality metrics were calculated after deriving stream networks from a 10-m resolution digital elevation model (www.nationalmap.gov; accessed 2 March 2015), splitting networks into 20-m stream segments, and normalizing values by dividing centrality by the total number of stream segments within a network. We used the r.watershed, r.stream.stats, and v.net.centrality functions in GRASS GIS for these analyses (version 7.0, GRASS Development Team 2012).
To allow for an effect of the overall watershed structure, we also calculated two metrics of stream branching structure at the level of the entire network (Dodds and Rothman 1999): Horton's length ratio (or the ratio of the length of tributaries to main channels, R T ), and the average number of major tributaries within a stream network (T 1 ), which is defined as being one Strahler order less than the main channel (e.g., where major tributaries are defined as being first-order streams that flow into a second-order stream; Fagan et al. 2010). These metrics are calculated based on all streams within a watershed, but they are appropriate for characterizing aspects of first-order streams because of the fractal nature of rivers (Rodriguez-Iturbe and Rinaldo 2001). For species that inhabit first-order streams, lower R T and higher T 1 reflect more paired first-order streams, which are associated with higher probabilities of occupancy and abundance for headwater-associated species Bolger 2002, Grant et al. 2009).

Niche modeling
Because the sampled population contains individuals of multiple age classes, we fit mixtures of normal distributions to size frequency data from individuals captured in occupancy surveys to identify age cohorts (MacDonald and Pitcher 1979). We verified the size distribution data using annual growth rate estimates from markrecapture data (Appendix S2). The presence of each age class was related to covariates using occupancy models (MacKenzie et al. 2006), fit with JAGS (http://mcmc-jags.sourceforge.net/) using the jag-sUI library in R (Kellner 2015, R Core Team 2015. These models account explicitly for the probability of detecting an individual of a given age class, given it is present during a sampling occasion, and thus provide unbiased estimation of occupancy and the relationship between occupancy and our covariates (Royle and Dorazio 2008). We fit models for each age class separately, and related occupancy to R T , T 1 , betweenness centrality, farness centrality, catchment area, and the first three principal components of stream habitat (Appendix S1), which accounted for 76% of the variation in stream habitat. We also included the first three principal components of stream habitat as covariates for detection probability because we believed larval salamanders were more difficult to capture in certain microhabitats (Crocker et al. 2007).
We visually inspected trace plots for evidence of good mixing among chains, and we examined theR statistic to evaluate model convergence (Kery and Schaub 2012). We ran three chains for 100,000 iterations with a burn-in period of 20,000 iterations, and we did not thin the results (Eaton and Link 2011). We evaluated covariate effects on occupancy by examining credible intervals of posterior distributions. More precise estimates have narrower credible intervals. We considered effects to be strong when 95% of the mass of a posterior distribution does not include 0.

RESULTS
Mixture models identified four age classes of larvae present across our 59 stream networks. Growth rates estimated from mark-recapture data for this species (Appendix S2) confirmed that these mixture models provided an accurate representation of the age cohorts within the population. The probability of detecting each age class during one sampling pass varied: P class1 = 0.05 (95% credible interval [0.02-0.15]); P class2 = 0.13 (95% credible interval [0.06-0.27]); P class3 = 0.37 (95% credible interval [0.30-0.45]); P class4 = 0.34 (95% credible interval [0.25-0.50]). In general, principal components of stream habitat (Appendix S1) were weakly related to detection probability. The only exception was the second principal component, which was associated with shallow, narrow stream reaches and which had a negative effect on the detection probability of age 2 larvae (À0.95, 95% credible interval [À1.65, À0.32]).
Age classes varied in their response to stream habitat (Fig. 1). Age classes 1 and 3 did not show a strong response to the first principal component of stream habitat (58.2% and 18.1% of posterior mass <0, respectively), but age classes 2 and 4 were less likely to occur in areas with higher amounts of cobble (99.4% and 95.9% of posterior mass <0, respectively). Age classes 2, 3, and 4 were more likely to occur in shallow, narrow reaches (99.9%, 73.4%, and 95.5% of posterior mass >0), but age class 1 was less likely to occur in this habitat (17.4% of posterior mass >0). Age classes 1 and 3 did not show a strong response to the third principal component, which was associated with reaches dominated by gravel and pebble streambeds (25.7% and 59.3% of posterior mass >0, respectively), while age class 4 was more likely to occur in this habitat (97.8% of posterior mass >0) and age class 2 was less likely to occur in this habitat (0.3% of posterior mass >0). Results for age classes 1 and 2 should be interpreted with caution because of the low numbers of individuals that were captured from these age classes. Age classes varied in their response to stream network structure and within-network position. Catchment area had a weak negative relationship with occupancy of age classes 3 and 4 (80.0% of posterior mass <0 for class 3, 92.5% of posterior mass <0 for class 4). Betweenness centrality was unrelated to ❖ www.esajournals.org the presence of age 3 individuals, but positively associated with the presence of age class 4 (58.2% of posterior mass >0 for age class 3, 90.0% of posterior mass >0 for age class 4), while farness centrality showed a positive relationship with the presence of age classes 3 and 4 (93.4% of posterior mass >0 for age class 3, 98.6% of posterior mass >0 for age class 4). R T was negatively associated with the occupancy of both age class 3 (5.3% of posterior mass >0) and age class 4 (92.6% of posterior mass <0), but T 1 did not have a strong relationship with either age class 3 (71.6% of posterior mass <0) or age class 4 (67.2% of posterior mass <0).
We created two sets of maps to illustrate the variation in patterns of occupancy for different age classes. First, we mapped the predicted occupancy for each age class based on models with average stream habitat conditions in the watersheds that had the highest and lowest R T (Fig. 2). Age classes 3 and 4 had similar occupancy patterns in a watershed with the lowest R T , but age class 4 was less likely than age class 3 to occur in tributaries near the outlet of the watershed that had the highest R T .
Second, we compared maps of the predicted occupancy for age class 4 in these two watersheds with maps of predicted occupancy that included one measure of network position and both measures of network structure, while fixing the other two measures of network position and measures of stream habitat at their average values (Fig. 3). This was done because the relative isolation of a site within a habitat network might be less important in a landscape where paired headwaters are more common. When a variable has a strong effect on occupancy patterns, a map showing the predicted effect of that variable with other predictors held at average conditions will resemble the map produced with all terms within a model, which is similar to calculating partial correlation coefficients for parameters of a regression model. The predicted pattern of occupancy from the full model was most closely matched by the effect of catchment area in the watershed with the lowest R T and the effect of farness centrality in the watershed with the highest R T .

DISCUSSION
Habitat choice during early life stages has implications for overall population dynamics and trends, though detecting these preferred niches using static distribution data may be difficult for species with overlapping generations, long larval life stages, and for species living in continuous habitat networks. Theoretical models for ONS often focus on the role of body size and transitions between different habitat types (Nakazawa 2014), but network structure may place constraints on these optimal ONS by restricting movement among habitat types (Carrara et al. 2014, Bertuzzo et al. 2015. Although the spatial scale at which populations respond to habitat has been examined for some species (Ayll on et al. 2010), the importance of the topological relationships among patches within habitat networks remains poorly understood. We used size frequency data to identify age classes in populations of the stream-associated spring salamander, which allowed us to identify spatial ONS related to local-and landscape-scale factors. Our findings illustrate that ONS can be driven by habitat structure at large spatial scales.
The location of sites within stream networks was an important predictor of ONS in larval spring salamanders. While betweenness centrality and farness centrality were more strongly associated with the presence of the fourth age class than with that of the third age class, catchment area had a consistent negative effect on the presence of these older ages. When predicted occupancy is mapped across networks, it is apparent that the branching structure of stream networks mediates the effect of the position of stream segments within networks (Fig. 3). In a watershed with more paired headwaters (Big Run, Fig. 3), occupancy patterns are well represented by the response to catchment area, but in a watershed with less complex topology (Jeremy's Run, Fig. 3), farness centrality reflects the variation in predicted occupancy. This indicates that the isolation of stream segments matters more in a watershed that contains more isolated firstorder streams. While a previous study of stream fish found that extinction risk declined as Horton's length ratio, R T , (the ratio of the length of tributaries to main channels) increased (Fagan et al. 2010), we found lower occupancy rates for stream salamander larvae in watersheds with higher R T . This difference may result from different preferences for habitat position (e.g., main channel vs. tributary) among species, suggesting that habitat isolation is more important for headwater species. In watersheds with higher R T , first-order streams ❖ www.esajournals.org are more isolated from one another, while largerorder streams are always adjacent to one another. This difference is compounded by the upstream bias in movement of many species that inhabit first-order streams (Grant et al. 2007). All age classes of spring salamanders have an upstream bias in movement, especially the adults (Lowe 2003), so selection for network positions with higher centrality allows for better access to headwater habitats, which are the main sites for reproduction. Occupancy site selection of younger age classes may indicate that those positions in the ❖ www.esajournals.org 6 February 2017 ❖ Volume 8(2) ❖ Article e01662 stream network are better for survival and growth (Cecala et al. 2009).
Our modeling approach shows that larval salamanders choose sites based on local habitat after accounting for network topology. Both the third and fourth age classes were associated with shallower, narrower reaches. Although the third age class was not associated with reaches having Fig. 3. Predicted occupancy for age class 4 of larval spring salamanders (Gyrinophilus porphyriticus) in Big Run (panels a-d, lowest Horton's length ratio) and Jeremy's Run (panels e-h, highest Horton's length ratio) from the full model (top) and predicted occupancy (w) when considering only network structure and catchment area (second row), network structure and farness centrality (third row), and network structure and betweenness centrality (fourth row).
❖ www.esajournals.org 7 February 2017 ❖ Volume 8(2) ❖ Article e01662 more cobble, gravel, or pebble substrates, the fourth age class was positively associated with these conditions. Indeed, we find that habitat characteristics of stream segments were unrelated to their position within stream networks (Appendix S1). The vegetation structure of terrestrial habitats may be affected by patch area and patch isolation via factors like edge effects or seed dispersal limitation (Rodewald 2003, Laurance 2008. In contrast to this, stream habitats are often the result of stream flow and other fluvial processes, so the range of habitats associated with smaller streams can be similar in watersheds that vary in their size or the relative isolation of headwater streams (Benda et al. 2004). Understanding spatial ONS improves our ability to identify intraspecific variation in distribution, can elucidate cryptic spatial dynamics in populations with multiple life history stages, and presents new opportunities for empirical tests of ecological theory. The complexity of a habitat network is related to metapopulation capacity (Bertuzzo et al. 2015) and persistence (Grant 2011) in theoretical landscapes. In species with prolonged larval or juvenile life stages, spatial ONS may further enhance metapopulation capacity and persistence by maximizing survival in complementary locations in the habitat network. Our results show that larval salamanders appear to have distributions that shift upstream as they age, which is likely to put them in areas with fewer predators and lower resource levels.