Fine-scale spatial variation in fitness is comparable to disturbance-induced fluctuations in a fire-adapted species.

The spatial scale at which demographic performance (e.g. net reproductive output) varies can profoundly influence landscape-level population growth and persistence, and many demographically pertinent processes such as species interactions and resource acquisition vary at fine-scales. We compared the magnitude of demographic variation associated with fine-scale heterogeneity (< 10 m), with variation due to larger-scale (> 1 ha) fluctuations associated with fire disturbance. We used a spatially explicit model within an IPM modeling frame-work to evaluate the demographic importance of fine-scale variation. We used a measure of expected lifetime fruit production, EF , that is assumed to be proportional to lifetime fitness. Demographic differences and their effects on EF , were assessed in a population of the herbaceous perennial Hypericum cumulicola (~ 2600 individuals), within a patch of Florida rosemary scrub (400 m x 80 m). We compared demographic variation over fine spatial scales to demographic variation between years across six years after a fire. Values of EF changed by orders of magnitude over < 10 m. This variation in fitness over fine spatial scales (< 10 m) is commensurate to post-fire changes in fitness for this fire-adapted perennial. A Life Table Response Experiment indicated that fine-scale spatial variation in vital rates, especially survival, explains as much change in EF as demographic changes caused by time-since-fire, a key driver in this system. Our findings show that environmental changes over a few tens of meters can have ecologically meaningful implications for population growth and extinction.


INTRODUCTION
The spatial scales at which plant population performance is estimated, and the extent to which vital rates (survival, growth, reproduction) vary through space, have important ramifications for the robustness of demographic estimates (Gurevitch et al. 2016). At regional scales, the aggregate behavior of collections of populations can vary widely from any individual population within that collection (Gurevitch et al. 2016, Hui et al. 2017, and the overall extinction risk may be much lower than that indicated by the average population performance (Abbott et al. 2017, Hui et al. 2017, Dibner et al. 2019. Moreover, when scaling from individuals to small patches and subpopulations, up to the entire species range, there is often no natural scale marking the end of one population and the beginning of another. Identifying the scale at which different subpopulations behave independently of each other is important because averaging over small-scale variation may mask important relationships with environmental drivers, and change our understanding of the population dynamics (Clark 2003, Diez et al. 2014, Vindenes and Langangen 2015. For example, a few fast-growing individuals can disproportionately contribute to population-level growth rates (Zuidema et al. 2009), and even if the broad-scale climate is unsuitable for a species, populations may still persist in favorable microhabitats (Csergő et al. 2017).
Many drivers of plant demographic performance vary at fine scales. The intensity of interactions with other organisms is a function of the immediate composition of neighboring vegetation that determines the strength of competition or facilitation (Casper and Jackson 1997). Neighboring associations also affect the abundance and composition of seed predators, herbivores, and microbiota (García and Chacoff 2007). Nutrient and groundwater availability can also differ between locations at short distances (Jackson and Caldwell 1993). Demographic variation at these small scales can have profound consequences on overall species performance at regional or landscape levels (Levin 1992).
There is, however, a trade-off in determining the appropriate scale to measure and model natural populations. Working at fine scales, ideally at the individual level, can substantially increase the cost and complexity of both data collection and analysis (Diez et al. 2014, Gurevitch et al. 2016. Further, using increasingly finer scales is likely to yield diminishing returns in terms of the insights and predictive performance gained because of spatial autocorrelation (Legendre 1993). In most cases, the trade-off between efficient sampling and averaging over important variation (Diez et al. 2014) is set a priori through a combination of study-system knowledge, goals, and practicality. Surprisingly, despite the potential consequences for demographic estimation and population forecasting, the effects of these decisions on our understanding of demographic processes are rarely tested.
Not all vital rates contribute equally to demographic performance (de Kroon et al. 1986), and they may not vary equally over a given spatial scale (Jongejans et al. 2010). Thus, even knowing how vital rates vary over space is not sufficient, because integrative measures of performance such as population growth rate (λ) or lifetime reproductive output (R 0 ) may show a different spatial pattern to that of the underlying vital rates. Attributing variation in population performance to variation in vital rates over space can identify relevant population drivers, highlight which vital rates drive differences in fitness at different scales, identify why spatial mismatches between population performance and its components occur, and show the relative importance of spatial vs. temporal variation (Jongejans et al. 2010).
The spatial covariance between vital rates can also be important for population-level demographic performance. To maintain demographic performance over environmental gradients, individuals can trade off lower performance in one vital rate for higher performance in another, a strategy known as demographic compensation (sensu Villellas et al. 2015). Demographic compensation results in negative covariance between vital rates across space, and can stabilize reproductive performance across a landscape Morris 2010, Villellas et al. 2015). Although demographic compensation is often thought of as occurring over larger spatial scales Morris 2010, Villellas et al. 2015), mechanistically it happens at the individual level, and so might, in principle, act over fine spatial scales.
To capture spatial variation in demographic performance adequately, we must estimate both the integrative metrics of population performance and the spatial scales over which they vary, rather than the scale being assumed a priori (Clark 2003, Ibáñez et al. 2007, Diez et al. 2014. Here, we estimate the spatial scale at which individual-level demographic performance and expected fitness varied using an individual-level data set for a single population of H. cumulicola, listed as federally endangered in the United States. We used a spatial randomization of this data to test if this population of H. cumulicola exhibits demographic compensation over small spatial scales. We use E F , expected lifetime fruit production, as our measure of integrative demographic performance. We expected to identify significant demographic variation in H. cumulicola at relatively fine scales of 10-20 m because of high heterogeneity in plant associations in the Florida scrub. In this ecosystem, it is common to transition from wetlands to dry vegetation in <50 m. To assess its importance, we compare population variation at these fine scales to the annual demographic variation associated with fire, an extremely important driver in this ecosystem (Quintana-Ascencio et al. 1998, 2003, Dolan et al. 2008.

Study site and species
Hypericum cumulicola is a herbaceous perennial plant species endemic to Polk and Highland counties, Florida, USA. This short-lived plant is mostly limited to open sandy areas (hereafter gaps) between shrubs within Florida rosemary scrub vegetation, composed mainly of Florida rosemary (Ceratiola ericoides), small-stature oaks (Quercus inopina, Quercus geminata, and Quercu myrtifolia) and palmettos (Sabal etonia and Serenoa repens). We defined gaps as continuous open areas >1 m 2 without shrubs >0.5 m in height (as in Menges et al. 2008Menges et al. , 2017a. Gap boundaries were defined by shrubs <0.5 m apart. Hypericum cumulicola flowers and fruits mainly from June to September. Fruits mature as red capsules with many seeds (Quintana-Ascencio and Morales-Hernández 1997). Primary seed dispersal is by gravity, leading to low levels of migration between gaps and high levels of inbreeding (Dolan et al. 1999).
The Florida rosemary vegetation (hereafter Florida rosemary patches) and its immersed gaps constitute a recognizable unit within the Florida scrub ecosystems. The Florida rosemary patches occur within a matrix of less xeric and more continuous and dense vegetation. The numerous gaps within the shrub matrix, which are expanded by fire and contract with time-since-fire, offer habitat for many endemic herbaceous and subshrub species and ground lichens (Menges et al. , 2017a. Hypericum cumulicola is one of several species specializing in Florida rosemary scrub that has greater demographic performance and less chance of extinction in regularly burned landscapes (Quintana-Ascencio et al. 2003, Menges andQuintana-Ascencio 2004).
The study site is located at Archbold Biological Station (270°10'50'' N, 81°21'0'' W, 42 m above sea level [a.s.l.]), Highlands County, Florida. All sampling occurred in a single Florida rosemary scrub patch on a xeric ridge, extending ca. 400 m from north to south and 80 m at its widest (see Appendix S1: Fig. S1). Prior to a fire in 2001, the site was previously burned in July 1986.
Fire is the dominant ecological disturbance in this community (Menges et al. 2017b, Quintana-Ascencio et al. 2018 and has a strong effect on H. cumulicola demography. Fire kills individual H. cumulicola, although the species has a persistent seedbank that allows it to survive fires. Fire also releases surviving and germinating individuals from competition with the surrounding woody matrix (Quintana-Ascencio and Morales-Hernández 1997). Consequently, recruitment and population growth are highest during the first years postfire, and then decrease with time-since-fire (Quintana-Ascencio et al. 1998, 2003, Dolan et al. 2008

Demographic and spatial data
Following the 2001 fire, we located, tagged, and mapped every individual H. cumulicola within the study site between January and August of 2002. The location of every individual was mapped with a laser (Impulse, Laser Technology Inc., Englewood, Colorado, USA, 1-cm accuracy). Every new recruit was also mapped and recorded in August of 2002August of , 2003August of , and 2004. We recorded the annual survival, reproductive condition (vegetative or reproductive) and plant height (in millimeters) of all H. cumulicola individuals between 2002 and 2004. Between 2005 and 2007 we randomly sampled gaps, and mapped and recorded the annual survival, reproductive condition (vegetative or reproductive) and plant height of all individuals within the chosen gaps (Appendix S1: Tables S1-S4, Appendix S1: Figs S2-S4).

Statistical analyses
Our measure of integrated fitness is asymptotic expected lifetime fruit production of a new plant, assuming the environment is constant (E F ). E F depends directly on survival and fruit production. In addition, both survival and fruit production are affected by plant size in H. cumulicola (Quintana-Ascencio et al. 2018), and so growth also indirectly affects reproductive output. To model survival, growth, and fruit production, with the ultimate goal of estimating E F , we used a generalized linear, spatial errors modeling framework. There were slight differences between the modeling approached for each vital rate, because survival (s) and probability of reproduction (r) are binary variables and growth (g) is a continuous variable. We first describe the linear predictor, which is common to all three vital rates and includes the spatial error term. The linear predictor has a year-specific intercept for vital rate v ∈ s,r, g ½ , β v 0,j , where each observation i was made in year j ∈ 2002,2003,2004,2005,2006,2007 f g , that includes the effects of both year and time-since-fire (Quintana-Ascencio and Morales-Hernández 1997, Quintana-Ascencio et al. 1998, 2003, Dolan et al. 2008). In our data set, year and time-since-fire are confounded. The linear effect of an individual's height on vital rate v is β v , and h j i is the height of individual i in the previous year for survival and growth, and the current year for reproduction. For survival and reproduction, we assume a common slope on h j i across years. In the growth models, the response is the height of an individual, predicted by its height in the previous year; therefore, the slope terms on the previous years' heights β g j can be interpreted as the average growth between surveys. To allow individuals to grow at different rates in different years, we let the slope terms on height in the previous year, β g j , vary between years. A fire in 2001 killed most standing individuals, opening up the study site, favoring large recruitment and reducing competition (Quintana-Ascencio et al. 1998, 2003, Dolan et al. 2008. As a result, growth in the year immediately following the fire was much higher (see Fig. 1, "Growth").
The spatial error term for vital rate v, W v q i ð Þ, can be thought of as a random effect that varies continuously over space. The spatial error term is computationally challenging to fit to the full set of observation locations.
To make the set of locations (q i ∈ q) smaller and more computationally tractable, we combine nearby observations so that densely packed individuals share a location, and thus the estimate for the spatial error term. To start, the x-y coordinate position of each observation is designated as a location, the set of all locations is q. To shrink this set: 1. Each q i ∈ q acted as the target location in turn. 2. The average position of all locations in q within 50 cm of the target location were added to q as a new aggregate location.

We then removed the original locations aggregated in
Step 2 from q.
Because this process was iterative as it progressed, the locations being aggregated could be the x-y coordinates of individual plants, or previously aggregated locations. In constructing the set of locations, q, our aim was to have a smaller number of locations, where the distance between every observation and its nearest neighboring location was small, and with small distances between locations. We used 552 locations in q, with good coverage of the data. All observations had a location in q within 56 cm, and 75% had a location within 25 cm. Seventy-five percent of locations had another location within 82 cm, and 50% had at least one other location within 61 cm; thus after combining locations we still had submeter resolution; see Appendix S2 for more details.
The spatial error term is drawn from a multivariatenormal distribution that lets nearby locations influence the estimate of location q i , while simultaneously estimating the scale of that influence: 0 is a vector of 0's the same length as q and the covariance matrix between locations in q is Ω q;α s , σ 2 s À Á . The covariance is modeled as an exponentially decreasing function of distance so that the covariance at the mth row and kth column is where d m,k ð Þ is the Euclidean distance between locations m and k, σ 2 v is the variance, and α v controls the rate that covariance between locations decays with distance.
For the binary variables, survival (s) and probability of flowering (r), we assumed the data are drawn from the Bernoulli distribution B Á ð Þ and the likelihood is where η v i,j is given in Eq. 1. For growth (g), we assume height in the current year is drawn from a normal distribution and the link function is identity where σ g is the standard deviation on the error term.
For the year effect of binary variables, we use a weak prior, β r,s 0,j ∼ N 0, 20 ð Þ. For all other β parameters in the linear predictor (Eq. 1) and σ g we use a flat prior. The variance component of the spatial error term is drawn from a weakly informative Cauchy prior, σ 2 v ∼ Cauchy 0, 5 ð Þ, truncated at 0.00001. The correlation component of the spatial error term α v is sampled on the inverse scale to improve sampling efficiency; α v ¼ 1=α * , where α * is sampled from the weak prior α * ∼ Cauchy 0, 5 ð Þ, truncated at 0.1 and 1000. We truncate the prior just above 0 to improve numerical stability of the sampling. The lower limit of truncation results in α v ¼ 10, a distance decay rate where points just 30 cm apart are independent. Because this distance is smaller than the resolution of our combined location points (Appendix S2) it is pointless to sample at smaller spatial scales. The upper truncation limit (α v ¼ 0:001) implies that all knot locations co-vary, that is, a site level effect, and so there is little point in sampling higher.
Fruit production.-The number of fruiting bodies (assumed to be proportional to seed production; Quintana-Ascencio et al. 2018) produced by each individual was not observed at the study site during the study period. Instead, we use data from 15 sites (including the study site), located in the same general area at the Archbold Biological Station (Quintana-Ascencio et al. 2018) and with similar ecological conditions, across 9 yr (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) to estimate the number of fruiting bodies based on plant height, with a separate intercept chosen for each study site (common slope on height estimated across sites).
where f i is the predicted number of fruiting bodies produced by individual i, β f is the effect of height (h i ) on fruit number, and β f 0 is the intercept for fruit study site f .

Height effect
Year effects Height effect Year effects

Height effects
Year effects 1. Effect sizes of year (β v 0j ) and height (β v and β g j ) on survival (s), probability of flowering (r), and growth (g). Note that in order to show the effect of height on the same scale as the intercepts for survival and flowering probability the slopes are the increase in probability for an individual with a height of the upper quantile (41 cm) compared with the mean height (32 cm), and for growth this is the effect of height for an average-sized individual. Note that the effect sizes for survival and reproduction are on the logit scale and on the arithmetic scale for growth. Article e03287; page 4 S. R. COUTTS ET AL. Ecology,Vol. 102,No. 4 We use the intercept fit to our study site in all analysis. Fruit numbers are counts, and we use a Poisson error distribution; μ f i is the expected number of fruits for individual i. All parameters are drawn from flat priors.
Initial size distribution.-To calculate integrated lifetime fruit production of a new plant, an initial distribution of new plant sizes is needed (Ellner et al. 2016). In 2001, a fire killed all aboveground individuals, and so we assume all individuals in 2002 were less than 1 yr old. We use a log-normal distribution fit to the observed 2002 size distribution.
All models were fit using Hamiltonian MCMC in stan, using the 'rstan' interface (Carpenter et al. 2017). Trace plots and effective sample size, along with residual plots and fitted versus predicted plots were checked for all models (not shown).

Individual-level performance
One of the most widely used measures of individual demographic performance is reproductive output (Caswell 2001, Ellner et al. 2016). We do not have data on seed demography over the study period in this population, so we use cohort-specific expected fruit production of a new plant in year j at location q i , E ij F , as an integrative measure of reproductive performance. We assume this is proportional to lifetime seed production. We constructed an integral projection model (IPM; Easterling et al. 2000) to estimate per-individual fruit productivity (E F ). Briefly, an IPM is a demographic model where a population distribution is projected through discrete time intervals as a function of one (or more) continuous state variables, in this case height (Easterling et al. 2000). Because we calculate expected E F from this IPM at an individual level (as opposed to the population level), distributions of E F over size arise from the distribution of future sizes and demographic fates predicted by the model, rather than variation within a population.
We use asymptotic E ij F , which assumes that the environment, and its effects on vital rates, at location i stays as it was in year j over the long term, an assumption we know is not supported in this system (Quintana-Ascencio et al. 2018). Thus, when interpreting E ij F one must note that it is based on a theoretical construction of demography, and not a literal prediction of the number of fruits produced per capita.
Adapting the calculation of asymptotic R 0 in Ellner et al. (2016:58-61) to fruit production, I is an identity matrix, S z, j, q i ð Þis the survival kernel, and G z, z 0 ,j,q i ð Þis the growth kernel from size z to size z 0 in the next year, both over size z in year j and location q i . The term I À S z, j, q i ð ÞG z, z 0 , j, q i ð Þ ð Þ À1 is the fundamental operator, and it calculates the expected amount of time an individual will spend at each size over their life, conditional on starting life at size z (see Ellner et al. 2016:61 for a more in-depth explanation). The expected number of fruits produced by an individual of size z in year j, location q i is F z, j, q i ð Þ . The initial size distribution of an individual was We used the same initial size distribution at all 552 locations, so we can attribute any spatial structure in E ji F to spatial structure in the vital rates.
In Appendix S3, we describe the fecundity, survival, and growth kernels in more detail.

Life table response experiment
To examine which vital rates drive differences in E ij F , and whether changes in vital rates between years or over space were most responsible for change in E ij F , we apply a fixed-effects two-way life table response (LTRE) analysis (Caswell 2001) to estimate how changes in the parameters of vital rate models contribute to differences in ln (E ij F ). We use a fixed-effect LTRE to visualize and interpret the contributions of locations and specific years to deviations in ln(E ij F ). We use a linear approximation to decompose differences in ln(E ij F ) over year and location relative to the mean year and location (Caswell 1989). Decomposing the differences in E ij F on the log scale greatly improved the accuracy of the linear approximation to our nonlinear IPM.
Following the notation of Caswell (1989) ln E ij rate v, for year j, averaged over location, and K •j ð Þ are the kernels built under those location-averaged parameters.
The total contribution of change in the parameters of the vital rate models at location q i to change in E F isθ i , and the total contribution of the change in parameters of the vital rate models in year j to change in E F isγ j . d θ i ,γ j À Á is the interactive effect of changes to parameters of the vital rate models over both year and location, over and above the additive effects. Finally δ ln EF is the sensitivity of ln E F ð Þ to vital rate parameter β vÃ , under kernels built with vital rate parameters at the midpoint between the values averaged over year and location and the value at location q i and year j. We used a finite difference approximation to calculate partial derivatives (Ellner et al. 2016:219).
We used a randomization procedure to test if the distributions of contributions to changes in ln E F ð Þ of location, time, and the interaction between the two, had different dispersion (i.e., were some of the distributions "wider" than others?). We used median absolute deviation (MAD) as the measure of dispersion. MAD is robust to extreme values in the tails and differences in the number of samples used to estimate it. This is an important feature in our case, because there are 546 locations for which contributions are estimated, but only 5 yr. The randomization tests the null hypothesis that both sets of contributions are drawn from the same distribution. The randomization ignores parameter uncertainty and uses the mean of the posterior for estimated contributions for each location, year, and location:year combination (for the interaction). The randomization procedure combined two sets contributions (pairs of location, time, or the interaction) into a joint vector. This joint vector was resampled with replacement 10,000 times into two vectors (C 1 and C 2 ) the same size as the original two vectors. The test statistic D rand ¼ abs MAD C 1 ð ÞÀMAD C 2 ð Þ ð Þ is calculated for each randomly generated pair of contribution vectors. The resulting distribution of differences in dispersion under the null model is compared to D obs , the test statistic calculated for the two observed sets of contributions. To reduce the number of pairwise comparisons, we only test the "full" contribution of location, time, and their interaction (θ,γ and d θ i , γ j À Á , respectively). To simplify the notation in the results and discussion, we drop the year and location superscripts from E i,j F and refer to reproductive performance as E F .

RESULTS
Year (and by proxy time-since-fire) clearly affected both survival and growth, with the years immediately following the 2001 fire being the best for growth and survival in Hypericum cumulicola (Fig. 1).
The scale of covariance in the spatial error term (W v q i ð Þ in Eq. 1) was similar for all three vital rates (survival, s, growth g, and probability of flowering, r), with correlation between locations falling to 0.5 (moderately correlated) at roughly 5 m and is less than 0.2 (nearly independent) at 10-m distance (Fig. 2). This suggests that similar, small-scale, environmental processes are influencing survival, growth, and reproduction. Because 69% of the occupied gaps were larger than 10 m 2 and 88% of the total plants occurred in gaps larger than 10 m 2 , the observed scale of spatial correlation indicates that vital rates can vary considerably within the sandy gaps in which H. cumulicola occurs (Fig. 3d-f).
Although each vital rates varies over a similar spatial scale, they do so in different ways (Fig. 3). In the "example gap" at the top left of Fig. 3d, survival is better than average after year and height effects are controlled for. The spatial structure of growth and probability of flowering in that same example gap is more complex (Fig. 3, e&thinsp;f respectively). This gap is <10 m in length and at one end there is a higher-than-average probability of flowering (yellows) and a lower-than-average growth rate (blues), whereas at the other end of the gap probability of flowering is lower than average and growth is higher than average.
Reproductive output (E F ) varied more than an order of magnitude over less than 10 m (Fig. 4). This level of spatial variation suggests that most of the reproductive output in our study population of H. cumulicola comes from a few small but highly productive areas within gaps. Recall that E F is the asymptotic seed production assuming vital rates remain constant at a given location, in the case of Fig. 4 vital rates in the average year for which data were available (1-6 yr postfire). In this system, biological conditions for H. cumulicola tend to worsen with increasing time-since-fire (Quintana-Ascencio et al. 1998, 2003, Dolan et al. 2008. Therefore, in this study E F is a model-based, relative measure we expect to be proportional to individual fitness, not a literal prediction of lifetime fitness. The contribution of changes in vital rates over space to changes in E F was of similar magnitude to the contribution of changes in vital rates between years ( Fig. 5a "full" compared to Fig. 5b "full"). The randomization test showed that the difference in the dispersion (as measured by MAD) of these two distributions (full contribution of location and time) was not different to that expected under the null hypothesis that both sets of contributions were drawn from the same distribution (P value = 0.21). The contribution to changes in E F of the interaction between location and year, over and above the additive effect of each, was modest. The dispersion of the distribution of contributions of the interaction term, location:time was significantly smaller than that of both location and time (P values from the randomization were <0.0001 in both cases). There were larger differences between locations in the contribution of the interaction effect in the immediate postfire years 2003 and 2004, with the effect fading in subsequent years Article e03287; page 6 S. R. COUTTS ET AL. Ecology,Vol. 102,No. 4 ( Fig. 5c "full" distributions become more concentrated around 0 with increasing time-since-fire). Changes in survival (particularly its spatial component) were the most important source of deviation in E F from the year and location-averaged reference, followed by changes in growth. Variation in probability of flowering over either space or between years had no substantial effect on E F (Fig. 5). Height had a large effect on probability of flowering (Fig. 1), and most H. cumulicola quickly achieved a height that made flowering likely (~15 cm), even in locations that had a worse than average probability of flowering.
We used a randomization test to show there was no evidence of demographic compensation (Villellas et al. 2015) at our study site (Appendix S4: Fig. S1). Although there was negative covariance between the spatial error terms for probability of flowering and both survival (covariance −0.27) and growth (covariance −0.35), this did not affect the variation of E F over locations because variation in probability of flowering had a negligible effect on E F (Fig. 5). Further, the two vital rates that did drive the observed spatial variation in E F , growth and survival, did not co-vary spatially, and so there was no potential for demographic compensation.

DISCUSSION
We show that fine-scale spatial variation in vital rates leads to large variation in an integrative measure of reproductive success over small spatial scales. A few locations inside a few semi-independent gaps (168 gaps), within a single continuous Florida rosemary scrub patch had reproductive outputs (measured by E F ) an order of magnitude larger than locations less than 10 m away. We use an LTRE to show that these spatial effects are at least as big as the effect of a good or bad year. This is an especially important finding in our data, which covers a period immediately following a fire. Postfire demographic dynamics are known to have a large influence on the population performance of H. cumulicola and its persistence (Quintana-Ascencio and Morales-Hernández 1997, Quintana-Ascencio et al. 2003, Dolan et al. 2008, Quintana-Ascencio et al. 2018). Our results suggest that the magnitude of spatial variation in vital rates and fitness over short spatial scales (<10 m), within gaps, is similar to that of postfire demographic changes in this population. Although mostly devoid of vegetation aboveground, these environments may have strong interaction and resource gradients associated with the belowground distribution of the roots of the dominant shrub species in the gap boundaries.
We currently have a poor understanding of how finescale demographic variation translates to landscape-level dynamics like extinction and population growth (Fieberg and Ellner 2001). Because such small areas make a disproportionately large contribution to the reproductive output of the whole landscape, we might expect landscape-level reproductive output to be highly variable and unpredictable. Our study species exists in discrete sandy gaps that are opened up by fire, and if not killed by fire the surrounding shrubs encroach and outcompete H. cumulicola Morales-Hernández 1997, Dolan et al. 2008). Our work shows that as well as the fire-return interval driving temporal population dynamics, the locations of the gaps that the fire opens up is a key driver of performance at the landscape scale. This spatial-temporal interaction may influence how population performance changes with time-since-fire. For example, if one of the few highly productive locations is encroached by shrubs, the overall reproductive output will be greatly reduced.
In the case of H. cumulicola, almost all seeds are dispersed within the parent gap (Dolan et al. 1999, Quintana-Ascencio et al. 2019. The limited dispersal of H. cumulicola means highly productive regions may result in small areas with very large seed banks, which could make them key sites for postfire regeneration (Quintana-Ascencio et al. 2018). These dynamics would require spatially referenced data that spans several fire cycles to unpick fully. An important caveat to these results is that we did not have data for seed production per fruit, seed bank mortality, and germination rate, and so could not close the life cycle. Thus, although we estimate that a H. cumulicola individual in the most productive areas can produce >20,000 fruit over its lifetime (assuming favorable postfire conditions would persist), we do not have the data to determine how many new individuals each H. cumulicola can be expected to produce, that is, R 0 (Ellner et al. 2016). It may be that there is important covariance between the vital rates that contribute to fruit production and seed mortality and germination. If that  Ecology,Vol. 102,No. 4 covariance is negative (areas of high fruit production are worse than average for seed survival and germination), then the spatial variation in R 0 will be less extreme than that seen in E F . On the other hand, if regions that are better than average for growth and survival are also better than average for germination and seed survival, then fine spatial variation in R 0 will be even more extreme than that shown for E F . Important population processes, such as widespread population declines, extinctions, and population spread, are most meaningful at landscape and regional scales. At these scales, the overall landscape growth rate tends to increase, and variance in that growth tends to decrease, as the number of populations in that landscape increases, due to the diversified portfolio effect (Dolan et al. 2008, Hui et al. 2017. Under the diversified portfolio effect the more independent populations there are, the greater the chance that at least some of them will display high positive growth rates in any given year (Dibner et al. 2019), and poor growth in a few populations can be offset by a larger number of average populations. The diversified portfolio effect depends on both spatial variation in habitat quality, which drives differences in population performance across the landscape, and how the spatial pattern of that habitat quality changes over time. We identify a landscape where there is potential for portfolio effects to operate because the environment varies over very small spatial scales, driving large differences in individual-level reproductive performance. However, we also Median E F (x1,000) find that the contribution of spatial variation in vital rates to changes in reproductive performance attenuates with time-since-fire. This suggests that in this system the potential for the portfolio effect changes over time and is greatest in immediate postfire years. Demographic compensation has been shown to stabilize and maintain reproductive performance across environmental gradients over larger spatial scales Morris 2010, Villellas et al. 2015). In principal, demographic compensation might also act over demographically meaningful environmental gradients at smaller spatial scales, such as those found in this study. However, our spatial randomization test shows no evidence of demographic compensation in this population. If vital rate covariance (the basis of demographic compensation) is a heritable trait (Villellas et al. 2015), then demographic compensation may not occur at such small scales, because gene flow is likely to be too high to allow differentiated demographic strategies.
Demographic models based on the spatial errors framework, like those we develop, will need further developments in demographic analysis to exploit. For example, extending the fixed-effects LTRE analysis we performed to decompose the variance in reproductive performance or fitness (a spatial random effects LTRE) directly would allow the estimation of how demographic contributions of vital rates shift in importance at different distances. , and the interaction between location and time (c). "Full" panels show the total variation in ln(E F ) attributable to variation in all vital rates, with location, time, and interaction given byθ i ,γ j , and d θ i ,γ j À Á , respectively, in Eqs. 10-12. (a) The contribution of each location to change in ln(E F ) relative to the year and location-averaged reference. For location variation in the vital rates survival, growth, and probability of flowering arises from the spatial error term W v q i ð Þ, v ∈ s,g,r ½ , respectively (Eq. 2). Violin plots show the total variation, including both the spatial variation and parameter uncertainty from the posteriors of the vital rate parameters. Golden points show the mean for each location. (b) The contributions of each year to change in ln(E F ) relative to the year and location-averaged reference. For time differences in the vital rates survival, growth, and probability of flowering arise from the year-specific intercepts β v 0,j , v∈ s,g,r ½ , and differences in growth slope come from the yearspecific slopes in the growth model β g j (Eq. 1). Uncertainty in the contribution of each year is from uncertainty in vital rate parameters. (c) The contribution of the interaction between location and year to change in ln(E F ), over and above the additive contributions of location and year, relative to the year and location-averaged reference. Violin plots in (c) show total variation including both spatial variation and parameter uncertainty from the posteriors of the vital rate parameters. Golden points show the mean contribution of the interaction between location and year for each location within each year (i.e., without the effect of parameter uncertainty). Article e03287; page 10 S. R. COUTTS ET AL. Ecology,Vol. 102,No. 4 The scales over which demographic rates and performance have important implications for what inferences we can draw about the landscape-level performance we are often concerned with, such as regional scale invasions and extinction risk (Clark 2003, Dolan et al. 2008, Diez et al. 2014, Abbott et al. 2017, Hui et al. 2017. Determining the right scale to model demography at, and collecting the necessary spatially explicit data over large enough extents, at a fine-enough resolution, is a nontrivial hurdle, and is a task made more difficult by the need for long temporal data sets to describe fully how spatial variation changes over time. Even with five annual transitions (6 yr of data) our results are a snapshot of the spatial structure that has emerged in this population. However, more such snapshots across systems and species would allow us to estimate what scales different types of species typically vary over. Spatial information is routinely collected in demographic studies, often in the form of nested study site, quadrats, and subquadrats locations. This spatial information is also routinely discarded in demographic analysis. Tools such as the spatial errors framework presented here, or those used by Clark (2003) (which is less demanding of the data), can be used to exploit this spatial information. These approaches will allow for more robust estimates of population dynamics (Clark 2003) and vulnerability (Hui et al. 2017), and insights only possible with a spatially explicit analysis, such as the identification of highly productive locations (as shown here).