# Effects of environmental and individual variation on time patterns of extinction in small experimental populations

Souleyman Bakker, Tom J M van Dooren, Romain Peronnet, Thomas Tully

► To cite this version:

Souleyman Bakker, Tom J M van Dooren, Romain Peronnet, Thomas Tully. Effects of environmental and individual variation on time patterns of extinction in small experimental populations. Oikos, In press, 10.1002/oik.11280. hal-05286951

HAL Id: hal-05286951

https://hal.science/hal-05286951v1

Submitted on 27 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research article

# Effects of environmental and individual variation on time patterns of extinction in small experimental populations

Souleyman Bakker<fnref n="1" />, Tom J. M. Van Dooren<fnref n="2" />, Romain Peronnet<fnref n="2" /> and Thomas Tully<fnref n="1" />

Correspondence: Tom J. M. Van Dooren (tvdooren@gmail.com)

Oikos

2025: e11280

doi: 10.1002/oik.11280

Subject Editor:

Paulo R. Guimaraes

Editor-in-Chief:

Paulo R. Guimaraes

Accepted 31 July 2025

<figure><img src="image_2.png" /></figure>

Predicting the fate of small populations is essential in ecology, epidemiology and conservation biology. Small populations can go extinct quickly or can develop into large established populations. These can still go extinct through demographic stochasticity, especially when declines in mean population size or large size fluctuations drive them into an extinction vortex. Individual differences and environmental conditions are two factors expected to determine extinction probabilities in recently founded populations. In populations surviving this initial phase, effects of initial conditions are expected to wash out, unless these differences in conditions persist systematically. We used the parthenogenetic Collembola *Folsomia candida* to investigate such effects in a crossed treatment on 290 replicate populations, each initiated with a single individual. Populations were initiated with founders of three different life stages and subjected to four levels of culling. Their extinction times were recorded in a period of up to 41 weeks after initiation. We fitted the parameters of an age-structured multi-type branching process model to extinction rates observed in the first weeks after initiation. As predicted by the model, extinction probabilities increased with the level of culling. Populations started by non-reproducing individuals had larger extinction probabilities which decreased over time. In established populations, extinction probabilities increased with culling levels, and the effects of the founder stage disappeared, except for a late increase in extinction probability in populations founded by newborns. This increase was not expected and forced us to reconsider when transient effects of the initial state would end. We hypothesize that size-structured competition in this system causes long-term transients. The results show that predictors of extinction probability can be shared between small and established populations. Transient dynamics of extinction risk can be long, with phases even mimicking a stationarity followed by the onset of an extinction vortex.

Keywords: branching processes, Collembola, declining populations, environmental effects, extinction vortex, *Folsomia candida*, long transients, small population paradigm

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Institut d'Écologie et des Sciences de l'Environnement de Paris (iEES-Paris), Sorbonne Université, CNRS, INRAe, IRD, Université Paris Créteil, Université Paris Cité, Paris, France</footnote>

<footnote n="2">²Archéozoologie, Archéobotanique : Sociétés, Pratiques et Environnements (AASPE), Muséum National D'Histoire Naturelle, CNRS, Paris, France</footnote>
# Introduction

A better understanding of the dynamics and fate of small pop-
ulations is central in conservation biology (Caughley 1994).
Small populations can pose environmental threats, as in the
case of invasive species which require control (Mack et al.
2000, IPBES 2023), or present conservation challenges, as
in the case of endangered or reintroduced species (Caughley
1994, Armstrong and Seddon 2008, Ralls et al. 2018).

Numerous environmental factors can affect the establish-
ment or long-term persistence of both small and large popu-
lations (Schwartz et al. 2006, Chevin et al. 2010, Primack
2012). They include habitat quality and size (Griffen and
Drake 2008), predation (Holyoak et al. 2000), or environ-
mental stochasticity (Stacey and Taper 1992, Halley and
Iwasa 1998, Lande et al. 2017). For example, the establish-
ment success of small introduced populations in Australia was
determined by predation pressure (Short 2009, Moseby et al.
2011). In established populations, these factors can reduce
population size or cause fluctuations in population density
and the increased demographic stochasticity can drive them
into an extinction vortex (Gilpin and Soulé 1986).

Individual variation in a small population can also affect extinction risk, for example through processes like inbreeding depression (Shaffer 1981, Stacey and Taper 1992). Other factors related to individual characteristics, such as age, size or life stage tend to have more transient effects, as they are not inherited across generations. For example, Sarrazin and Legendre (2000) used individual-based simulations to compare the efficiency of releasing juvenile versus adults in a griffon vulture reintroduction program in the French Alps. Their simulations indicated that releasing juveniles increased extinction probabilities relative to releasing adults. In the field, however, the survival cost of release disappeared after one generation (Sarrazin et al. 1994, 1996). Conditional on non-extinction, and if population regulation occurs as populations become large, it is expected that the numbers of individuals in different individual states will converge to a quasi-stationary distribution (QSD) and that the probability of extinction per unit of time will become constant (Gosselin and Lebreton 2000). The same holds for total population size (Belovsky et al. 1999, Gosselin and Lebreton 2000, Drake et al. 2011). The initial transients in the population dynamics are no longer relevant in larger, stationary non-extinct populations. Multi-type branching processes (Athreya and Ney 1972, Caswell 2001, Schreiber and Ross 2016) without density-dependence can serve as models for these transient phases where population regulation has not yet taken effect or remains negligible. Small and newly established populations often do not experience strong density dependence, but over time, density dependence will lead to changes in mean fecundity and survival rates (Coulson et al. 2000), eventually bringing the population growth rate closer to an average of zero (Sibly and Hone 2002), potentially increasing its variance in population size (Ferrer and Donazar 1996). Population regulation can already occur during the transient phase. For example, in springtails, the monopolisation of resources by large adults

can inhibit the maturation of juveniles (Le Bourlot et al. 2014). When juveniles fail to reach maturity, extinction risks can be extended over time.

Caughley (1994) already called for an integration of approaches for small populations that may become larger and large ones at risk of becoming small. Multi-type branching process modelling of regulated populations (Gosselin and Lebreton 2000, Lebreton et al. 2007) provides a theoretical backbone for this. Extinction probabilities are larger during transient phases starting from small populations and decrease to constant and lower extinction rates for the surviving populations reaching QSD.

In this study, we aim to test whether extinction proba-
bilities during a transient phase are accurately predicted by
multi-type branching processes without density dependence.
For presumably established populations, we further explore
whether the effects of individual state (such as age or size)
wash out over time, while the effects of systematic differences
in environmental conditions between populations persist.
Additionally, population declines driven by environmental
change (Caughley 1994), for example by the accumulation
of waste products or mould in experimental settings, could
undermine non-extinct populations and alter the expected
distribution of non-extinct population states. If this degrada-
tion occurs, there might be more than two distinct phases,
each with its own pattern of extinction probabilities. The
third phase might then correspond to an extinction vor-
tex (Gilpin and Soulé 1986) with extinction probabilities
increasing over time.

We extended a previous experiment on the extinction
of small populations (Van Dooren et al. 2024), using the
Collembola species *Folsomia candida* (Willem 1902), to a
much longer duration. We monitored replicated popula-
tions over nearly ten months, allowing to examine extinction
probability patterns over both short and longer time frames.
We developed a new multi-type branching process model to
obtain predictions tailored to this experiment. The duration
of this experiment allowed us to explore whether the tran-
sient phase, where this model is expected to apply, lead to
the establishment of the QSD regime. Each population was
founded by a single individual from one of three different age
classes, ensuring that a transient phase towards a stationary
state would occur. The populations were subjected to differ-
ent levels of culling, for which we anticipated different popu-
lation growth rates in small populations and distinct average
population sizes at stationarity.

Material and methods

Experimental design

We used the parthenogenetic springtail *Folsomia candida*
(Collembola, Isotomidae) as our experimental model organ-
ism (Fountain and Hopkin 2005) to study how founder age
and culling intensity affected extinction probability patterns
over time. *Folsomia candida* is a widespread soil-dwelling
arthropod which also inhabits caves and mines and which
can be found both at very low and high population densities. It shows indeterminate growth without clear metamorphosis (ametabolism). Life history statistics can be found in Fountain and Hopkin (2005) and Tully (2023). Average lifespan and total fecundity are temperature-dependent and reach over 200 days and 1000 eggs in favourable conditions. Egg production increases with body length, up to around 100 eggs per clutch on average. At 21°C, newborns measure approximately 0.1 mm in length, attain maturity in about ten days after hatching and reach up to 3 mm. This species comprises two distinct clades, A and B, which exhibit genetic differences (Tully et al. 2006), morphological variation (Tully and Potapov 2015), and divergent life-history traits (Tully and Ferrière 2008, Mallard et al. 2015, 2020, Tully 2023). The extinction experiment was carried out using a single clonal strain, 'GM' from clade B (Tully et al. 2006), known for its shorter lifespan and increased fecundity in comparison to clade A (Tully 2023). There could therefore not be any effects of genetic variability on the results.

We maintained the *F. candida* populations in standardised rearing containers made from polyethylene vials (52 mm in diameter, 65 mm in height) filled with a 30 mm layer of plaster of Paris mixed with 600 µl of Pebeo graphic Indian ink. The moistened plaster ensured nearly 100% relative humidity, while the ink enhanced the visibility of the white individuals and their clutches of eggs. All populations were kept at a constant temperature of 21°C in an incubator. Dried pellets of yeast (Tully and Ferrière 2008) were provided as food and regularly replaced to ensure ad libitum feeding conditions throughout the experiment.

To have founders of the required age classes, we isolated eggs from stock populations and raised cohorts by providing food immediately after they hatched. We used three types of founders: 1) newborns, which were small individuals that had hatched less than two days before using them, but had not been fed beforehand and therefore had not been able to start growing, 2) small adults, two weeks after hatching and capable of reproducing, and 3) large adults, fully grown but not yet showing signs of senescence, at least five weeks post-hatching (Fig. 1A). Small and large adults are both fertile and lay a clutch approximately every week (Fig. 1A). However, small adults produce smaller clutches than large adults, as fecundity increases with body size in *F. candida* (Fig. 1A, Supporting information in Tully 2023).

We founded populations at two time points, two weeks apart, allowing us to monitor 290 populations. At the first time point, we founded populations with single small or large adults, at the second with single newborns. The populations in a specific treatment group had to go through a transient phase before potentially reaching a quasi-stationary distribution across different population states. Our founder age classes were confounded with date. Newborn founders might have differed from adult founders not just in their age, but in other individual properties due to a different history before the start of the experiment. The populations founded by newborns shared the environment during the experiment with all other populations founded two weeks earlier.

Populations were either left untreated (control) or subjected to one of three culling treatments applied at each census. For culling, a plastic template was used to randomly isolate a quarter, half, or three quarters of the plaster surface area, and all individuals and visible eggs within these sections were removed. Every week, all populations were inspected to determine their extinction status, culled according to treatment, observed again to determine their extinction status after culling, and then provided with fresh food before being returned to the incubator.

This experimental protocol began one week after the populations were founded and continued until the end of the experiment, 41 weeks after the first populations were founded. We analysed the population states as binary extinct/present (0/1) response variables and calculated extinction times based on these observations. Among the explanatory variables, we also recorded the time intervals between censuses, these were not always exactly seven days due to logistical constraints. Occasionally, populations that initially appeared extinct following culling were later found to be non-extinct, due to the hatching of missed eggs. In these cases, we corrected the response profile by replacing zero values with ones for time points before the last individuals were observed alive.

## Age-structured branching process model with imperfect culling

We used an age-structured multi-type branching process (MBP) model, to predict extinction probabilities in our experiment (Athreya and Ney 1972, Van Dooren et al. 2024, Supporting information). Parameter estimates were derived from a previous experiment by Tully (2023, Fig. 1A, Supporting information). The MBP model was based on an age-structured matrix population dynamical model (Caswell 2001, Van Dooren et al. 2024). In our experiment, we mainly selected founder individuals according to life stage. We therefore also labelled the age classes with life stages and used the predictions of the age classes thus labelled ‘newborns’, ‘small adults’ and ‘large adults’. Note that these labels ignore within age-class heterogeneity in life stages, and that individuals can proceed to a different life stage before the next census. Individuals were categorised into six distinct age classes, each with its own fecundity ($f_j$) and survival ($s_j$) parameters. The life cycle diagram corresponding to the Leslie matrix (Leslie 1945) of this model for the treatment without culling is given in Fig. 1B. The model assumed that the time interval between observations is one week. Individuals observed for the first time make up the first age class. They are eggs of between zero and seven days old. At each census from then on, individuals in any age class contributed the eggs laid since the previous census, if any, to the first age class and advanced themselves to the next age class. Figure 1A shows the ages in days of the individuals in each age class at census. Following Van Dooren et al. (2024), fecundities in the first two age classes were set to zero in a first parameterisation (Supporting information), as reproductive maturity in *F. candida* at 21°C occurs from about 22 days since birth, based on observations by Tully (2023). In a re-estimated parameterisation
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 1. (A) Growth curves for individuals of the GM clone and other clonal lines belonging to the same clade. Age is given in days since egg laying and the age boundaries of each age class at census in the model of panel (B) are drawn. Thin grey lines show individual growth trajectory, the black line shows the average growth trajectory. The dots placed on a growth trajectories represent the clutches of that individual and their area is proportional to the clutch size. The colours are used to distinguish successive clutches. Data from Tully (2023). (B) Life cycle graph for *F. candida* corresponding to the age-structured matrix model in the absence of culling (Supporting information). The time interval between censuses is one week. This life cycle shows the different ages classes at census, of which the boundaries are given in (A) and a life stage description of each of them, which we use to denote the different types of founders. Individuals can change life stage in between censuses. The arrows express the contributions of each age class at census to the numbers of individuals in each age class at the next census (survival in blue, fecundity in grey). The thickness and intensities of the fecundity arrows represent the increase in average fecundity with the age and size of the individuals. The age classes used as founder individuals in our experiment are given a different coloured background. The fecundity in each age class is a weighted sum of the numbers of eggs produced by individuals in that age class between successive censuses. Therefore, fecundity *f*₃ is certainly non-zero because juvenile individuals in age class three mature before the next census. Individuals in age class two (labelled 'newborns' but in fact a mixture of eggs and newborns), are between 7 and 14 days old at census. At the next census, they are between 14 and 21 days old, which is one day away from the earliest observed age at maturity in (A). A fecundity *f*₂ is drawn for this age class, because we used the data of the experiment to estimate it anew and assess if maturation might have occurred earlier.</figcaption></figure>
below, we did allow for a non-zero fecundity in the second age class as the boundary of the juvenile age class and the earliest age at maturity are little different. When individuals reached the sixth age class, they remained in that class from then on (Fig. 1B). In this model, we ignored actuarial senescence (mortality increases with age) because its impact on mortality is negligible compared to the effects of culling treatments. Additionally, reproductive senescence is minimal in the clonal line used in the experiment (Tully 2023). Thus, we assumed constant mortality and fecundity rates for older age classes.

The model proposed by Van Dooren et al. (2024) assumed that all mature individuals produce a clutch during each interval between two culling events. However, after observing that some intervals between clutches exceeded one week (Fig. 1A), we introduced additional parameters $p_j$, representing the probability than an individual lay a clutch within a given week ($0 \le p_j \le 1$). If an individual produces a clutch, the clutch size is modelled as a Poisson-distributed variable with mean $F_j$. The mean fecundity across all individuals in age class $j$ is therefore $f_j = p_j F_j$. In the model parameter estimation detailed below, we adjusted $F_j$ in response to changes in $p_j$ such that the expected $f_j$ remained as in the Supporting information. This allowed for a zero-inflated fecundity without a change in the mean and thus no effect on the population growth rate.

Culling was quantified by the probability that an individual escapes culling $\gamma$, which took values of 1 (no culling), 0.25 (75% culling), 0.5 (50% culling), or 0.75 (25% culling). Given that some clutches were overlooked during culling, we introduced parameter $\delta$, corresponding to the proportion of clutches missed during culling ($0 \le \delta \le 1$). The probability that a clutch survived culling, $\gamma'$, was thus modified to $\gamma' = \gamma + \gamma(1-\gamma)$. The matrix population dynamical model $M$, resulting from the multi-type branching process, is the sum of two matrices (Supporting information). The logarithm of the dominant eigenvalue $\lambda$ of this matrix is the asymptotic population growth rate $r$, and predicts whether populations will eventually become extinct with probability one ($r \le 0$) or have a probability of escaping extinction indefinitely ($r > 0$, Supporting information). In MBP models, population growth rates decrease with higher levels of culling applied to all life stages (Van Dooren et al. 2024). Additionally, we expected populations to converge to smaller sizes as $\gamma$ decreases. The model assumed that survival and reproduction were independent within each time interval and across individuals, that individuals survived and were culled according to Bernoulli distributions, and that clutch sizes followed Poisson distributions (Supporting information). We used the branching process model to predict extinction probabilities at each census for populations founded by a single individual, accounting for founder age, culling levels and census number. The predictions were not adjusted to account for variable durations between censuses. As time proceeds, it can be expected that the assumed independence between individuals breaks down when population sizes increase after founding due to density-dependence.

A mismatch between model predictions and actual extinction patterns thus becomes indicative of density dependence and an approach to regulated population states.

## Statistical analysis

For all statistical analyses, we used R ver. 4.1.2 (www.r-project.org), with scripts made available at the Zenodo Digital Repository: https://zenodo.org/records/14179307. The extinction times observed were interval-censored, as data collection occurred at specific time points.

## Extinction time distributions and survivorship curves

We first explored an approach proposed by Drake et al. (2011) to detect transient phases in extinction probabilities. We plotted histograms of extinction times on a logarithmic scale (i.e. based on the census times when extinctions were observed). In some cases, these can show bimodality even while extinction probabilities are expected to decrease continually (Drake et al. 2011). As proposed, we fitted a Cox proportional hazard model to the observed and right-censored extinction times to test the assumption that extinction hazards remain constant over time and are proportional between treatments. We did this using the tests proposed by Grambsch and Therneau (1994) and used by Drake et al. (2011), using the `cox.zph()` function from the 'survival' package (Therneau and Grambsch 2000). This analysis ignored the interval censoring of the data. To test for the presence of multimodality, we also applied Hartigans' dip tests (Hartigan and Hartigan 1985) to the uncensored extinction times per treatment group and corrected tail probabilities for multiple comparisons using Holm's method (Holm 1979).

We subsequently compared survivorship curves between treatments using nonparametric generalised log-rank tests for interval-censored data (Sun 1996), using the `ictest()` function from the 'interval' package (Fay and Shaw 2010). We used the data from censuses after culling only, integrating survival in between and at culling, because we wanted to compare these curves to MBP predictions. We tested the overall null hypothesis that no differences exist between any of the survivorship curves, and complemented this with tests for each founder stage within a culling treatment and each culling treatment within a founder stage.

## Fitting branching process parameters using the data

To improve the accuracy of the predictions from our branching process model, we first conducted a sensitivity analysis of ultimate extinction probabilities for all estimated model parameters except the culling survival probability, at parameter values used in Van Dooren et al. (2024, Supporting information). For the sensitivities to parameters $p_j$, the probabilities to reproduce, we systematically adjusted $F_j$ in response to changes in $p_j$ such that mean fecundity $f_j$ was not affected. For the populations founded by individuals of age classes two, four and six, the average fecundity for age class two ($f_2$), the survival parameters for age classes two and four ($s_2, s_4$), the proportion of eggs surviving culling ($\delta$) and clutch laying probabilities $p_4$ and $p_6$ had the largest
absolute values of sensitivities (Supporting information). Consequently, we focused on these six parameters to maximise the likelihood of the data as a function of the parameter values in the set we allowed to change, while maintaining the other parameters at the values provided in the Supporting information.

For each observation, we calculated the probability of observing an extinction, assuming that the MBP prediction with a given set of parameter values and for that ordered culling event was true. We maximised the likelihood (ML) of the data as a function of the free parameters which could change, using the *optim* function in R with the limited-memory modification of the BFGS quasi-Newton method (Byrd et al. 1995).

Because we didn't expect assumptions of this model to be satisfied throughout the experiment, we restricted the likelihood optimisation to the first ten culling events. This maximisation was repeated for each subset of the six parameters and we compared model fit among all subsets using AICc (Hurvich and Tsai 1989). We report the ML parameter estimates and their SD from the best model with the lowest AICc. We compared the predictions of the model with the lowest AICc with actual survivorship curves along the entire experiment, to see whether the predictions of the branching process model failed at some point, indicating that the assumptions of the model would break down as time proceeds. When the best parameter subset predicted negligible extinction probabilities beyond a certain census for a treatment combination, we considered the populations established. Populations can only establish according a MBP model when their population growth rate is above zero. Conversely, for subcritical regimes where the population growth rate is below zero, an approach to a quasi-stationary distribution (QSD) with a constant extinction probability over time is expected.

## Treatment and time effects on extinction probabilities

The approaches using Cox proportional hazards models or non-parametric tests above are not modelling the time dependence of extinction probabilities explicitly. We did so, using binomial generalized additive models (GAM; Wood 2017) and generalized linear models (GLM; Collett 2023). We exploited the fact that population survival status was observed both right before and right after each census.

For a population *k* observed after culling at times *t* and *t+1*, the probability that it survives until *t+1* when it is non-extinct at time *t* equals $p_{c,k}(t+1)p_{\Delta,k}(t)$, where $p_{\Delta,k}(t)$ is the probability of surviving the interval of length $\Delta$ between these two observations and $p_{c,k}(t+1)$ is the probability of surviving the culling event at the end of the interval. In survival analysis, the survival probability across an interval until culling is given by $p_{\Delta,k}(t) = e^{-S_{t,t+\Delta,k}}$, where $S_{t,t+\Delta,k}$ is the cumulative extinction hazard for population *k* between times *t* and *t+Δ*. Approximating this cumulative hazard by $S_{t,t+\Delta,k} \approx \Delta(h_k(t+\Delta/2))$, with $h_k(t)$ the time-dependent instantaneous hazard at the interval midpoint, we obtain:

$$ \ln(-\ln(p_{\Delta,k}(t))) = \ln(\Delta) + \ln\left(h_k\left(t + \frac{\Delta}{2}\right)\right). \quad (1) $$

Using extinction probability $p_{c,\Delta,k}(t)$ per interval as the response

$$ p_{\Delta,k}(t) = 1 - p_{c,\Delta,k}(t), \quad (2) $$

this approximation produces the model form of a binomial generalized additive model (GAM) or generalized linear model (GLM) with complementary log-log link which estimates hazards $h_k$, which depend on time and on other explanatory variables characterising population *k*. The log interval length must be used as an offset in these models. At culling, survival probability $p_{c,k}(t+1)$ can be modelled with GAM or GLM with any link function. This is an instantaneous event for which no smooth hazard can be integrated over time. We analysed the extinction probabilities at culling $p_{c,c,k}(t+1)$ with binomial GAM and GLM and using a logit link.

To compare time patterns of extinction probabilities between the early and late phases of the experiment, we introduced an explanatory dummy variable, named ‘period’. This variable takes the value one starting from the twentieth culling event on, dividing the data into ‘Early’ and ‘Late’ periods. This division allowed us to examine changes in extinction probabilities over time. The other explanatory variables included were the culling treatment, founder stage and the interval midpoints or the times at culling.

We initially modelled time-dependent effects using splines to allow for non-linear effects. Observing that the spline functions were approximately linear within periods, we present GLM analysis with linear time effects (of time at culling, interval midpoints). A full model with all explanatory variables and interactions was constructed and simplified, first by eliminating effects that increased the AIC, then by removing non-significant effects according to a likelihood ratio test. This resulted in the selection of a minimum adequate model (MAM). We inspected the parameter estimates and generated model predictions for the MAM, including SDs of each prediction, to calculate predicted extinction probabilities between one observation after culling and the next $1 - p_{c,k}(t+1)p_{\Delta,k}(t)$.

## Results

During the experiment, 290 populations were monitored. The numbers per treatment combination are given in Fig. 2.

### Extinction time distributions and survivorship curves

When we plotted the observed extinction times on a logarithmic scale, bimodality was only suggested in some cases (non-zero culling treatments with adult founders), even though all populations had to go through a transient phase. In agreement with this, Hartigans’ dip tests found departures
<figure><img src="image_6.png" /><figcaption>Figure 2. Histograms of the times at which extinctions were observed. Extinction times are plotted on a logarithmic scale for each treatment combination. Non-extinct populations censored at the end of the experiment are shown in blue and observed extinctions are shown in red. Sample sizes are shown per panel. They vary between culling treatments and slightly so between founder stages.</figcaption></figure>

from unimodality in five out of eleven treatment groups with extinction events: the populations which were 50% culled with large adult and newborn founders (tail probabilities both < 0.001) and populations which were 75% culled (tail probabilities < 0.001 for adult founders and 0.008 for newborn founders). However, the treatment combinations with significant dip tests did not always show a clear departure from unimodality, indicating that the test might pick up a signal due to the discreteness of the distribution at small extinction times (Fig. 2). The limited detection of departures from unimodality can further be due to small sample sizes or most extinctions from QSD might only occur after the end of the experiment. Populations might not reach QSD, which can occur when newborns founded populations and with large culling probabilities. Also, when there is no culling, extinctions can be absent in the transient phase (Fig. 2). The tests on Schoenfeld residuals of a Cox model for extinction times with an interaction between culling levels and founder age classes rejected the proportional hazards assumption (overall p < 0.0001), but plots of the time dependence didn't show anything else than a constant pattern (Supporting information). Therefore, this approach provided inconclusive evidence of different phases with specific time-dependent patterns of extinction risks.

A visual inspection of the extinction times and survivorship curves (Fig. 2–3) revealed that extinctions occurred until

the end of the experiment in any treatment. This was different from our expectations based on Van Dooren et al. (2024, Supporting information). There, the predicted population growth rate from the branching process model was below zero in the severe culling treatment (75% culled on average per week), meaning that culling should drive nearly all populations to extinction by the end of the experiment, regardless of founder age. For the other culling treatments, extinctions were expected to occur in the first month only (Fig. 3, solid lines).

A global test comparing all survivorship curves revealed significant differences between them ($\chi^2_{11} = 131.4, p < 0.0001$), indicating that extinction dynamics varied across treatments. Tests examining the effects of culling levels within founder stage and founder stage within culling levels were also significant ($p \le 0.001$), except for founder stages within the control treatment, which is likely due to the smaller sample sizes in the control treatment. These findings demonstrate that both the environmental treatment (culling) and differences at founding significantly influenced population extinction risks.

## Branching process modelling

When we compared MBP models with different subsets of parameters which were ML optimized, we found that adjusting all six parameters resulted in the lowest AICc (AICc=718). The second-best subset had AICc=722 ($s_2$,
<figure><img src="image_7.png" /><figcaption>Figure 3. Survivorship curves (proportions of non-extinct populations, coloured lines) with 95% confidence intervals (coloured bands) are plotted for each combination of founder stage (rows) and culling treatment (columns). Black lines represent the predictions of the different branching process models. Solid lines correspond to the MBP model of Van Dooren et al. (2024) with $\delta = 0$, dotted lines to the best model where the likelihood of the data was optimised as a function of parameters $f_2, p_4, p_6, s_2, s_4$ and $\delta$.</figcaption></figure>

$s_4, f_2, \delta$). During the first ten weeks of the experiment, the predictions of survivorship based on the subset with lowest AICc stayed within the confidence bands for most treatments (Fig. 3). Without any optimised parameters, this was less so (Fig. 3). From ten weeks after founding on, the predictions over- or underestimated survival probabilities. In the absence of culling, Van Dooren et al. (2024) calculated an asymptotic population growth rate of 1.17 per week; here we found 0.76 for the best MBP model (Supporting information). Population growth rates decreased with culling intensity and became subcritical below $= 0.31$ in Van Dooren et al. (2024) and below $= 0.40$ for the model with lowest AICc here. MBP model predictions converged to a flat line above probability zero, indicating that no further extinctions were expected, while in reality, extinctions continued to occur. This suggests that model assumptions were violated as time since founding progressed. The prime candidate for that seems population regulation due to interactions between individuals.

The parameter values estimated for the subset with the lowest AICc are: $\delta = 0.61$ (confidence interval CI [0.43–0.77]), $s_2 = 0.30$ (CI [0.22–0.40]), $s_4 = 0.94$ (CI [0.01–1.00]),

$f_2 = 0.45$ (CI [0.26–0.78]), $p_4 = 0.93$ (CI [0.50–0.99]), $p_6 = 0.79$ (CI [0.62–0.90]). We conclude that a significant fraction of clutches has been missed during culling, that the survival of newborns was lower than in the previous experiment (as compared to the Supporting information) and that the interval between successive clutches was longer than a week, especially for large adults.

## Statistical modelling of extinction probabilities

For extinctions occurring between culling events, we found a significant three-way interaction between founder stage, period and time ($\chi^2 = 46.0, p < 0.0001$). Inspection of the parameter estimates (Supporting information) showed that, during the early period, populations founded by age class two had an increased extinction probability intercept (parameter estimate 12.28, SE 4.79, $p=0.01$) and a strongly decreasing extinction probability with time, with a slope difference estimate of –0.29 (SE 0.079, $p=0.0003$). This slope was not significantly different from zero for the other age classes (Supporting information). On top of that, for populations founded by age class two, there was a gradually increasing
<figure><img src="image_8.png" /><figcaption>Figure 4. Estimated extinction probabilities and number of populations at risk over time for each founder stage and culling treatment (points). Probabilities were estimated after culling for each weekly interval and treatment interaction and combine extinctions during the interval and by culling. Solid lines correspond to the combination of inter- and post-culling extinctions predicted using the GLMs described above and summarised in the Supplementary Information. The bands correspond to the upper and lower bounds of the extinction probabilities (95% confidence intervals). Point sizes are proportional to the numbers of populations at risk.</figcaption></figure>

extinction probability (Fig. 4) with time throughout the experiment (slope parameter difference estimate 0.041, SE 0.021, p = 0.05). There were no significant effects of culling level on extinction probabilities in between culling events (Supporting information).

For extinctions occurring immediately after culling events, there were significant differences between culling levels ($\chi^2 = 133.33$, p < 0.0001), and a three-way interaction was identified between founder stage, period and time effect ($\chi^2 = 25.94$, p < 0.0001). The effects of culling levels on extinction probability were all significantly different from each other (Tukey post hoc pairwise comparisons), and extinction probabilities increased with the intensity of culling (Supporting information). Thus, the extinction probability increased with culling intensity in both the early and late periods, in line with the model-based population growth rates and expected population sizes. For populations founded by age class two in the early period, the extinction probability intercept was increased (parameter estimate 7.79, SE 2.77, P(>|z|) = 0.005). In the early period, there was an additional strongly negative slope of time in populations founded by age class two (parameter difference estimate -0.068, SE 0.016, P(>|z|) < 0.0001). Extinction probability increased over time in populations founded by age class two throughout the experiment (slope parameter difference estimate 0.027, SE 0.013, P(>|z|) = 0.033). The significant

effects involving populations founded with newborns were concordant for extinctions in between and at culling. They support an early phase with an initially high and then strongly decreasing extinction probability. The same populations showed an unexpected gradually increasing extinction probability towards the end of the experiment. Figure 4 shows the overall pattern of extinction probabilities based on the combined predictions of the models in the Supporting information.

## Discussion

Our study demonstrates that environmental differences which persist over time (culling levels) and differences between populations at the moment of founding can significantly influence population extinctions. These effects were observed in both the first and second half of the experiment. For culling levels, this was as expected as it affects population growth rates and population sizes. For effects caused by population differences at founding, we expected that these would be due to individual variation and would wash out over time as populations became established. This was not fully the case here: for populations founded by newborns, extinction probabilities slightly increased with time late in the experiment, while for other founder stages, they did not.
A multi-type branching process model adjusted to the experimental protocol managed to predict extinction probabilities well for some time, after re-estimating parameters based on the data likelihood. Some time after the start of the experiment, the trajectories predicted by the branching process moved out of the confidence bands of survivorship curves. We attributed this to violations of the assumptions of the MBP used, in particular population regulation.

## Environmental effects

The parthenogenetic reproduction of *F. candida*, its limited age at maturity and substantial clutch sizes should reduce the effects of demographic stochasticity. The remaining expected lifetime number of offspring of a young adult can exceed 350 offspring (Mallard et al. 2015). In the rearing containers we used, long-term population sizes of unculled populations typically range from a few dozen to about 200 adults and several hundred juveniles, with numbers rarely reaching critically low densities (Mallard et al. 2020). With culling and assuming no population regulation, the predicted MBP population growth rates decrease with culling intensity and growth rates become subcritical when $\gamma$ decreases. We also expected that stationary population sizes near QSD would be on average smaller with more culling, as a fraction of the population was removed each week. The constant probabilities of extinction at QSD (Gosselin and Lebreton 2000) should therefore be larger. When a population occasionally drifts to low numbers, the lowered growth rates due to culling imply that they will remain there longer. Together, this implies that such environmental treatments will affect extinctions and population survival during both transient and stationary phases of population establishment and that these will covary. In the absence of culling, we expected a negligible risk of extinction but found that about 21% of populations nevertheless went extinct by the end of the experiment. Our analysis based on the MBP predictions estimated a reduced newborn survival in comparison to previous experiments and that not all adults lay a clutch every week. This does not explain that some populations founded by large adults went extinct late in the experiment. It might be that the mechanism which increased extinction probabilities in populations founded by newborns late in the experiment is at play in all populations, but not as easy to detect in all of them due to differences in the strengths of effects. However, the (negative) non-significant parameter estimates in the Supporting information are not suggesting any increase and therefore do not support this.

## Individual variation

The analysis of extinction probabilities found that these were larger in the early period for populations founded by age class two, and that the difference with adult age classes diminished over time, as expected for factors with transient effects. The differences between the two adult stages were always smaller than those between adults and newborns. Larger extinction probabilities for newborn founders short after founding is similar to what Sarrazin and Legendre (2000) found in their simulations: extinction probabilities were larger when

releasing juveniles and not adults. In a founding phase, release and handling effects might increase extinction risk temporarily. Next to that, juveniles might suffer more from a lack of experience in the new environment or be more vulnerable in general. The expected duration between founding and the first new clutch or offspring which survives independently of the parent must be important as well in the case of an age effect. However, we note that in our experiment founder age was confounded with the date where populations were founded, so that we cannot unambiguously assign effects to founder age. Newborn founders might have differed in other traits from adult founders than just age. By re-estimating parameters, we found that some individuals observed in age class two might reproduce before the following census, in contrast to the conclusion of Van Dooren et al. (2024). The time it took for age class two founders to lay a clutch was less than expected but the fecundity estimate was low, less than one egg per parent on average. This and our reduced estimate of newborn survival in comparison Van Dooren et al. (2024, Supporting information) might have been driven by specific properties of the individuals used as newborn founders and indicative of a date effect due to a different individual history from the other founders instead of an age effect. Such a date effect which also corresponds to individual differences, if it has occurred, it has enlarged the age effect expected from the multi-type branching process model. Unexpectedly, we found an increasing extinction probability late in the experiment for populations founded by newborns. Here, one would expect that also date effects due to individual history would gradually wash out like founder age effects. If non-extinct populations remain near QSD after a transient phase, a late increase is not expected. The populations founded with newborns might differ in other ways than just by individual history up to founding, for example when systematically different rearing containers were used, but that was not the case.

## Time patterns of extinction probabilities

Newly founded populations go through an initial phase characterised by small population sizes and a larger vulnerability to demographic stochasticity. Our results suggest that the initial phase lasts about the duration of the entire experiment of Van Dooren et al. (2024), which didn't allow a study of the approach to QSD. In a second, later and longer-lasting phase, extinctions are less frequent and affect populations that have grown, appear established and have reached a QSD (Gosselin and Lebreton 2000). The tests proposed by Drake et al. (2011) to detect the existence of a transient period where populations approach a QSD were significant, but the results did not show a clear time pattern. A direct analysis using binomial GLM allowed modelling time effects explicitly and found similar results in the extinction probabilities in between culling events and at culling events. All populations had to go through a transient phase because they were all initiated from a single individual. However, when adults founded a population, few extinctions were expected during this transient such that it might go undetected when focusing on extinction times alone. There are even situations
imaginable where extinction probabilities during a transient
are smaller than at QSD. The late increase in extinction
probability suggested a driven population decline (Caughley
1994) and no stationarity. It should correspond to average
population sizes gradually decreasing again or increasing fluc-
tuations. This slowly appearing change could be due to non-
linear population dynamics in structured populations, but it
could also have been caused by an environmental degrada-
tion of some sort. Here, it only appeared for a single category
of founder individuals. This affects which mechanisms are
plausible.

Population regulation

Competition may decrease early-stage survival and fecundity,
thereby reduce population growth rates and consequently
increase the probability of extinction for populations which
experience density-dependence. An experiment similar to
this one conducted by Pike et al. (2004), demonstrated that
a heavily culled population (25% every two days) remained
at substantial risk of extinction regardless of its age and state.
This seems to correspond to what is called a subcritical regime
in branching processes, which occurred for the strongest cull-
ing regime here. Then non-extinct populations do converge
to a QSD, but they are probably almost never large enough
for density dependence to make much of a difference to pre-
dictions. We observed that the multi-type branching pro-
cess predictions did move outside of the confidence band of
observed extinction probabilities relatively quickly in some
cases and could reduce this to some extent by re-estimating
model parameters. For weaker culling regimes where the
MBP predicts supercriticality, populations which don't go
extinct during the initial critical period are expected to grow
and to reach a size sufficiently large to become more resistant
to extinction. Branching processes then predict a survivor-
ship curve converging to a horizontal asymptote and there-
fore no further extinctions. However, extinctions continued
to occur. In regulated populations, this is expected when
populations are at QSD with population growth rates near
zero and our results on populations founded by adults sug-
gest that they often reached it because extinction probabilities
did not change significantly with time. The lagged time effect
in populations founded by newborns was not in agreement
with being at a QSD. Here, we explore the hypothesis that it
might be due to a long transient in the population dynamics,
which is a regime in the dynamics occurring without a clear
trend in imposed environmental conditions, persisting for
longer than a few generation times but not or only partially
representative of the stable long-term dynamics which would
eventually appear (Hastings et al. 2018). Hastings et al.
(2018) proposed several broad mechanisms which can cause
long transients. Trajectories of the population dynamics can
be in the vicinity of an unstable equilibrium for a long while
(a ‘crawl-by’) which is known from predator–prey interac-
tions or interactions between competitors (Hastings et al.
2018). This transient can be affected by initial conditions,
which can initiate or prevent the convergence to the unstable
equilibrium. Le Bourlot et al. (2014, 2015) and Mallard et al.

(2020) observed that phases in the population dynamics of
*F. candida* kept in similar conditions can occur which last
around 200 days (almost 30 weeks), where adults prevent
the maturation of juveniles due to asymmetric size compe-
tition. Such size-structured competition (Le Bourlot et al.
2014) can explain the effect of date/founder age when this
shifts the long phase in time or changes its nature. Late-
occurring extinctions in springtail populations may also be
linked to the progressive accumulation of organic matter
in rearing containers. Incomplete food consumption and
accumulating faeces in a humid environment can facilitate
the development of mould. Personal observations suggested
competition between springtails and mould, with a competi-
tive disadvantage against mould at low springtail densities. It
might have caused additional mortality, especially in treat-
ments with high culling levels where the springtail popula-
tion density is expected to remain low. In the longer run,
a competitive disadvantage could drive, small populations
into an extinction vortex (Courchamp et al. 1999). This is
comparable to extinctions induced by a predator-driven
Allee effect in a *Daphnia–Chaoborus* predator–prey system
by Kramer and Drake (2010). It can also cause long tran-
sients due to time scale differences between different actors
in the system (Hastings et al. 2018), such as in the case of
viruses and hosts. For this to play a role in the late increase in
extinctions, this interspecific competition should be affected
by the initial transient phase specific to newborn founders
or by a date effect. This seems less plausible, as populations
founded by newborns lag only two weeks behind in their
dynamics behind those founded by small adults, but can't
be excluded in particular in combination with a date effect.
Similarly, long transients due to environmental and demo-
graphic stochasticity (Hastings et al. 2018) cannot explain
the lasting founder effect well as they are expected to wash
out across replicate populations but cannot be excluded.
Another mechanism proposed by Hastings et al. (2018) for
long transients is high dimensionality of the dynamical sys-
tem. Age-structured models with six age classes are not high-
dimensional, but size-structured competitive interactions can
be high-dimensional, as the number of possible individual
states and sizes, is much larger as in Le Bourlot et al. (2014).
MBP for size-structured regulated populations might be able
to predict long transients and detailed time series of popula-
tion size distributions seem required to understand and pre-
dict such long-term effects. The detection of lasting transient
effects might call for a reconsideration of the time scale on
which QSD are expected to be reached or remain stable.

**Conclusion**

Patterns of extinction in newly founded populations are
expected to show a transient phase and a stationary phase
with a fixed probability of extinction. Environmental effects
and individual variation can affect extinction probabilities.
These effects are in part predictable by multi-type age-struc-
tured branching processes and can be studied using statistical
modelling. Our results demonstrate that the transient phase
can lack a substantial extinction probability and that increases
in extinction probability over time can occur when popula-
tions are expected to be at a quasi-stationary distribution.
This increase can be due to long transients, which can mimic
QSD followed by the onset of an extinction vortex.

*Funding* – This study was financed by recurrent funding from the French government allocated to CNRS and Sorbonne University.

Author contributions

**Souleyman Bakker:** Formal analysis (equal); Visualization (lead); Writing – original draft (lead); Writing – review and editing (equal). **Tom J. M. Van Dooren:** Formal analysis (equal); Methodology (lead); Supervision (equal); Writing – review and editing (equal). **Romain Peronnet:** Data curation (equal). **Thomas Tully:** Conceptualization (lead); Data curation (equal); Formal analysis (equal); Methodology (equal); Visualization (supporting); Writing – review and editing (equal).

Data availability statement

Data are available from the Zenodo Digital Repository: https://zenodo.org/records/14179307 (Bakker et al. 2025).

Supporting information

The Supporting information associated with this article is
available with the online version.

References

Armstrong, D. P. and Seddon, P. J. 2008. Directions in reintroduc-
tion biology. – Trends Ecol. Evol. 23: 20–25.

Athreya, K. B. and Ney, P. E. 1972. Branching processes. – Springer

Bakker, S., Van Dooren, T. J. M., Peronnet, R. and Tully, T. 2025. Data from: Effects of environmental and individual variation on time patterns of extinction in small experimental populations. – Zenodo Digital Repository https://zenodo.org/records/14179307.

Belovsky, G. E., Mellison, C., Larson, C. and Van Zandt, P. A. 1999. Experimental studies of extinction dynamics. – Science 286: 1175–1177.

Byrd, R. H., Lu, P., Nocedal, J. and Zhu, C. 1995. A limited memory algorithm for bound constrained optimization. – SIAM J. Sci. Comput. 16: 1190–1208.

Caswell, H. 2001. Matrix population models, 2nd edn. – John
Wiley & Sons.

Caughley, G. 1994. Directions in conservation biology. – J. Anim. Ecol. 63: 215–244.

Chevin, L.-M., Lande, R. and Mace, G. M. 2010. Adaptation, plasticity and extinction in a changing environment: towards a predictive theory. – PLoS Biol. 8: e1000357.

Collett, D. 2023. Modelling survival data in medical research, 4th edn. – Chapman & Hall/CRC Texts in Statistical Science,

Coulson, T., Milner-Gulland, E. J. and Clutton-Brock, T. 2000. The relative roles of density and climatic variation on population dynamics and fecundity rates in three contrasting ungulate species. – Proc. R. Soc. B 267: 1771–1779.

Courchamp, F., Clutton-Brock, T. and Grenfell, B. 1999. Inverse density dependence and the Allee effect. – Trends Ecol. Evol. 14: 405–410.

Drake, J. M., Shapiro, J. and Griffen, B. D. 2011. Experimental demonstration of a two-phase population extinction hazard. – J. R. Soc. Interface 8: 1472–1479.

Fay, M. P. and Shaw, P. A. 2010. Exact and asymptotic weighted logrank tests for interval censored data: the interval R package. – J. Stat. Softw. 36: jssv036i02.

Ferrer, M. and Donazar, J. A. 1996. Density-dependent fecundity by habitat heterogeneity in an increasing population of Spanish imperial eagles. – Ecology 77: 69–74.

Fountain, M. T. and Hopkin, S. P. 2005. *Folsomia candida* (Collembola): a 'standard' soil arthropod. – Annu. Rev. Entomol. 50: 201–222.

Gilpin, M. E. and Soulé, M. E. 1986. Minimum viable populations: processes of species extinction. – In: Soulé, M. E. (ed.), Conservation biology: the science of scarcity and diversity. Sin-auer Assoc, pp. 13–34.

Gosselin, F. and Lebreton, J.-D. 2000. Potential of branching pro-
cesses as a modeling tool for conservation biology. In Ferson,
S. and Burgman, M. (eds), Quantitative methods for conserva-
tion biology. – Springer, pp. 199–225.

Grambsch, P. M. and Therneau, T. M. 1994. Proportional hazards tests and diagnostics based on weighted residuals. – Biometrika 81: 515–526.

Griffen, B. D. and Drake, J. M. 2008. Effects of habitat quality and size on extinction in experimental populations. – Proc. R. Soc. B 275: 2251–2256.

Halley, J. M. and Iwasa, Y. 1998. Extinction rate of a population under both demographic and environmental stochasticity. – Theor. Popul. Biol. 53: 1–15.

Hartigan, J. A. and Hartigan, P. M. 1985. The dip test of unimodal-
ity. – Ann. Statist. 13: 70–84.

Hastings, A., Abbott, K. C., Cuddington, K., Francis, T., Gellner, G., Lai, Y. C., Morozov, A., Petrovskii, S., Scranton, K. and Zeeman, M. L. 2018. Transient phenomena in ecology. – Science 361: eaat6412.

Holm, S. 1979. A simple sequentially rejective multiple test procedure. – Scand. J. Stat. 6: 65–70.

Holyoak, M., Lawler, S. P. and Crowley, P. H. 2000. Predicting extinction: progress with an individual-based model of protozoan predators and prey. – Ecology 81: 3312–3329.

Hurvich, C. M. and Tsai, C.-L. 1989. Regression and time series model selection in small samples. – Biometrika 76: 297–307.

IPBES 2023. Thematic assessment report on invasive alien species and their control of the intergovernmental science-policy platform on biodiversity and ecosystem services. – IPBES Secretariat

Kramer, A. M. and Drake, J. M. 2010. Experimental demonstration of population extinction due to a predator-driven Allee effect. – J. Anim. Ecol. 79: 633–639.

Lande, R., Engen, S. and Sæther, B.-E. 2017. Evolution of stochastic demography with life history tradeoffs in density-dependent age-structured populations. – Proc. Natl Acad. Sci. USA 114: 11582–11590.

Le Bourlot, V., Tully, T. and Claessen, D. 2014. Interference versus exploitative competition in the regulation of size-structured populations. – Am. Nat. 184: 609–623.

Le Bourlot, V., Mallard, F., Claessen, D. and Tully, T. 2015. A simple graphical method for displaying structured population dynamics and STdiag, its implementation in an R package. – Soil Organ. 87: 203–213.
Lebreton, J.-D., Gosselin, F. and Niel, C. 2007. Extinction and viability of populations: paradigms and concepts of extinction models. – *Ecoscience* 14: 472–481.

Leslie, P. H. 1945. On the use of matrices in certain population mathematics. – *Biometrika* 33: 183–212.

Mack, R. N., Simberloff, D., Mark Lonsdale, W., Evans, H., Clout, M. and Bazzaz, F. A. 2000. Biotic invasions: causes, epidemiology, global consequences, and control. – *Ecol. Appl.* 10: 689–710.

Mallard, F., Farina, M. and Tully, T. 2015. Within-species variation in long-term trajectories of growth, fecundity and mortality in the Collembola *Folsomia candida*. – *J. Evol. Biol.* 28: 2275–2284.

Mallard, F., Le Bourlot, V., Le Coeur, C., Avnaim, M., Péronnet, R., Claessen, D. and Tully, T. 2020. From individuals to populations: how intraspecific competition shapes thermal reaction norms. – *Funct. Ecol.* 34: 669–683.

Moseby, K. E., Read, J. L., Paton, D. C., Copley, P., Hill, B. M. and Crisp, H. A. 2011. Predation determines the outcome of 10 reintroduction attempts in arid South Australia. – *Biol. Conserv.* 144: 2863–2872.

Pike, N., Tully, T., Haccou, P. and Ferrière, R. 2004. The effect of autocorrelation in environmental variability on the persistence of populations: an experimental test. – *Proc. R. Soc. B* 271: 2143–2148.

Primack, R. B. 2012. A primer of conservation biology, 5th edn. – *Sinauer Assoc.*

Ralls, K., Ballou, J. D., Dudash, M. R., Eldridge, M. D. B., Fenster, C. B., Lacy, R. C., Sunnucks, P. and Frankham, R. 2018. Call for a paradigm shift in the genetic management of fragmented populations. – *Conserv. Lett.* 11: e12412.

Sarrazin, F. and Legendre, S. 2000. Demographic approach to releasing adults versus young in reintroductions. – *Conserv. Biol.* 14: 488–500.

Sarrazin, F., Bagnolini, C., Pinna, J. L., Danchin, E. and Clobert, J. 1994. High survival estimates of griffon vultures (*Gyps fulvus fulvus*) in a reintroduced population. – *Auk* 111: 853–862.

Sarrazin, F., Bagnolin, C., Pinna, J. L. and Danchin, E. 1996. Breeding biology during establishment of a reintroduced griffon vulture *Gyps fulvus* population. – *Ibis* 138: 315–325.

Schreiber, S. J. and Ross, N. 2016. Individual-based integral projection models: the role of size-structure on extinction risk and establishment success. – *Methods Ecol. Evol.* 7: 867–874.

Schwartz, M. W., Iverson, L. R., Prasad, A. M., Matthews, S. N. and O'Connor, R. J. 2006. Predicting extinctions as a result of climate change. – *Ecology* 87: 1611–1615.

Shaffer, M. L. 1981. Minimum population sizes for species conservation. – *BioScience* 31: 131–134.

Short, J. 2009. The characteristics and success of vertebrate translocations within Australia. – *Wildlife Research and Management Pty Ltd*, and Australian Government Department of Agriculture, Fisheries and Forestry.

Sibly, R. M. and Hone, J. 2002. Population growth rate and its determinants: an overview. – *Philos. Trans. R. Soc. B* 357: 1153–1170.

Stacey, P. B. and Taper, M. 1992. Environmental variation and the persistence of small populations. – *Ecol. Appl.* 2: 18–29.

Sun, J. 1996. A non-parametric test for interval-censored failure time data with application to AIDS studies. – *Stat. Med.* 15: 1387–1395.

Therneau, T. M. and Grambsch, P. M. 2000. Modeling survival data: extending the Cox model. – Springer.

Tully, T. 2023. Diversity, plasticity and asynchrony of actuarial and reproductive senescence in the Collembola *Folsomia candida* (Willem, 1902). – *Front. Ecol. Evol.* 11: 1112045.

Tully, T. and Ferrière, R. 2008. Reproductive flexibility: genetic variation, genetic costs and long-term evolution in a Collembola. – *PLoS One* 3: e3207.

Tully, T. and Potapov, M. 2015. Intraspecific phenotypic variation and morphological divergence of strains of *Folsomia candida* (Willem) (Collembola: Isotomidae), the 'standard' test spring-tail. – *PLoS One* 10: e0136047.

Tully, T., D'Haese, C. A., Richard, M. and Ferrière, R. 2006. Two major evolutionary lineages revealed by molecular phylogeny in the parthenogenetic Collembola species *Folsomia candida*. – *Pedobiologia* 50: 95–104.

Van Dooren, T. J. M., Haccou, P., Hermus, G. and Tully, T. 2024. Extinction probabilities of small structured populations: adequate short-term model predictions in *Folsomia candida*. – *bioRxiv* 2024.08.26.609669.

Willem, V. 1902. Note préliminaire sur les Collemboles des grottes de Han et de Rochefort. – *Ann. Soc. Entomol. Belg.* 46: 275–283.

Wood, S. N. 2017. Generalized additive models: an introduction with R, 2nd edn. – Chapman and Hall/CRC.