# Species-Specific Traits Shape Genetic Diversity During an Expansion–Contraction Cycle and Bias Demographic History Reconstruction

Ravi Vishwakarma, Gabriele Maria Sgarlata, David Soriano-paños, Rita

Rasteiro, Tiago Maié, Tiago Paixão, Rémi Tournebize, Lounès Chikhi

► To cite this version:

Ravi Vishwakarma, Gabriele Maria Sgarlata, David Soriano-paños, Rita Rasteiro, Tiago Maié, et al. "Species-Specific Traits Shape Genetic Diversity During an Expansion–Contraction Cycle and Bias Demographic History Reconstruction." Molecular Ecology, 2024, 34 (1), pp.e17597. 10.1111/mec.17597. hal-05390609

HAL Id: hal-05390609

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

Submitted on 1 Dec 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
Species-specific traits shape genetic diversity during
expansion-contraction and bias demographic history
reconstruction

Ravi Vishwakarma¹,*, Gabriele Maria Sgarlata¹,²,*†, David Soriano-Paños¹,³,*†, Rita Rasteiro¹,⁴, Tiago Maié¹,⁵, Tiago Paixão¹, Rémi Tournebize⁶,⁷,†, and Lounès Chikhi¹,⁶,⁸

¹Instituto Gulbenkian de Ciência, Oeiras, Portugal.

²Department of Evolution and Ecology and Center for Population Biology, University of California, Davis, California, USA.

³Universitat Rovira i Virgili.

⁴MRC Integrative Epidemiology Unit, University of Bristol, Bristol, United Kingdom.

⁵Institute for Computational Genomics, RWTH Aachen University, Germany.

⁶Centre de Recherche sur la Biodiversité et l'Environnement, CNRS, IRD, UPS, Université de Toulouse
Midi-Pyrénées, Toulouse, France.

⁷DIADE, IRD, Université de Montpellier, Montpellier, France.

⁸Centre for Ecology, Evolution and Environmental Changes (cE3c), Faculdade de Ciências da Universidade de Lisboa, Campo Grande, Lisboa, Portugal.

*shared first author

†Current affiliation: Department of Evolution and Ecology and Center for Population Biology, University of California (Gabriele Maria Sgarlata); Universitat Rovira i Virgili (David Soriano-Paños); DIADE, IRD, Université de Montpellier (Rémi Tournebize)

November 5, 2024

**Corresponding authors:**

Ravi Vishwakarma: rvishwakarma@igc.gulbenkian.pt,

Lounès Chikhi: lounes.chikhi@univ-tlse3.fr / chikhi@igc.gulbenkian.pt

David Soriano-Paños: sorianopanos@gmail.com

**Running Title:** Genetics of expansion and contraction

**Keywords:** demographic history, dispersal, glacial cycle, life-history traits, refugium, spatial pro-
cesses
Abstract

Species ranges are dynamic, experiencing expansions, contractions, or shifts in response to habitat changes driven by extrinsic factors such as climate change or human activities. While existing research examines the genetic impacts of spatial processes, few studies integrate species-specific traits to analyze how habitat changes affect co-existing species. In this study, we address this gap by investigating how genetic diversity patterns vary among species with different traits (such as generation length, population density and dispersal) experiencing similar habitat changes. Using spatial simulations and a simpler panmictic population model, we investigate the temporal genetic diversity in refugium populations undergoing range expansion, followed by stationary and contraction period.

By varying contraction speed and species traits, we identified three distinct temporal dynamics of genetic diversity during population contraction: i) a decrease in genetic diversity, ii) an initial increase followed by a decrease, and iii) a continuous increase throughout the contraction. We show that genetic diversity trajectories during population decline can be predicted by comparing sampled population diversity to equilibrium values expected under expanded and contracted habitat ranges.

Our study also challenges the belief that high genetic diversity in a refugium population is due to a recent and rapid habitat loss. Instead, we found contrasting effects of contraction speed on genetic diversity depending on the interaction between species-specific traits and the dynamics of habitat change. Finally, using simulated genetic data, we found that demographic histories inferred from effective population size estimates may vary across species, even when they experience similar habitat changes.
# Introduction

The current spatial distribution of species results from complex interactions between extrinsic (e.g., habitat and climate change) and intrinsic (e.g., dispersal) factors over timescales ranging from decades to millennia. Major climatic changes are thought to lead to significant changes in species distribution, leaving putative genetic signatures in the genomes of populations (Hewitt, 1996, 2004). Recently, human activities have also been altering species distributions (Finn et al., 2023), leading to range contractions (Laliberte and Ripple, 2004; Wolf and Ripple, 2017; Ceballos et al., 2017), expansions (Pacifici et al., 2020; Beyer and Manica, 2020) or fragmentation (Cushman, 2006; Hadad et al., 2015). Understanding the genetic consequences of changing habitats, population sizes (Nei et al., 1975) and connectivity (Alcala and Vuilleumier, 2014) is crucial, given the central role of genetic diversity in determining a species capacity to adapt to new environmental conditions and the negative effects of inbreeding in small populations (Frankham, 2005). Investigations into the genetic effects of population size changes have a long history. For instance, both theoretical (Nei et al., 1975; Maruyama and Fuerst, 1985; Allendorf, 1986) and empirical (Cabe, 1998; Dures et al., 2019; Leigh et al., 2019) studies found that population declines reduce genetic diversity in natural populations. However, some empirical studies have also observed little to no loss (Welch et al., 2012; Maebe et al., 2016a; Crates et al., 2019) or even an increase in the genetic diversity of species experiencing population declines (Townsend et al., 2023).

The apparent disconnection between genetic diversity and contemporary population decline in a number of studies presents a significant hurdle, particularly when trying to incorporate genetic diversity into policy and conservation practices (Hoban et al., 2021; Laikre, 2010). One reason is that it remains challenging to separate and identify the contribution of major evolutionary factors such as drift, migration, mutation, selection and life-history traits in shaping the genetic diversity of a species (De Kort et al., 2021). To investigate the evolution of genetic diversity amidst the intricate interplay of demographic and ecological factors over space and time, spatially explicit models have been particularly useful. They not only provide a better representation of natural populations but also complement, validate, or extend analytical results that are often limited to much simpler models (Edmonds et al., 2004; Klopfstein et al., 2006; Schneider et al., 2010; Rasteiro et al., 2012; Szép et al., 2022). While spatial simulations may be more difficult to interpret due to the potentially large number of parameters, they have yielded new insights into the genetics of spatial expansions (Klopfstein et al., 2006; Edmonds et al., 2004; Schneider et al., 2010) or contractions (Leblois et al., 2006; Pflüger et al., 2019; Rogan et al., 2023; Arenas et al., 2012; Sgarlata et al., 2022a) and have significantly contributed to our understanding of how spatial processes can influence the evolution and distribution of genetic diversity. However, with few exceptions (Rasteiro et al., 2012; Garnier and Lafontaine, 2022), they generally do not incorporate life-history traits into their analyses. Consequently, the extent to which genetic diversity patterns vary across different species undergoing the same habitat changes remains to be elucidated.

To address this gap, we simulated scenarios inspired by climatic oscillations of the Quaternary. During these periods, habitats and geographical ranges of species expanded from refugial areas and subsequently contracted into refugia over the glacial-interglacial periods (Hewitt, 1996; Schmitt, 2007). Similar refugia probably played a vital role in allowing species to survive unfavorable periods and to recolonize the landscape when conditions improved (Sommer and Nadachowski, 2006; Hewitt, 2004). We used spatial simulations and analytical derivations under a panmictic population model to investigate the temporal dynamics of genetic diversity within refugium populations. Specifically, we investigated how these dynamics were affected by a range expansion followed by a contraction, considering various contraction speeds and species-specific traits. We built upon the simulation study of Arenas et al. (2012), by considering comparable habitat changes. We first focused on the effect of generation length variation and its impact on genetic diversity patterns across different scenarios of habitat changes. We then simulated habitat changes using more realistic species-level traits (such as population density, dispersal distance, and generation length) based on reported allometric relationships of mammalian species. Allometric relationships mathematically describe how biological characteristics scale with body mass or other measurable traits in organisms (Peters
and Peters, 1986; Blueweiss et al., 1978). Our simulations allowed us to (a) describe the temporal dynamics of genetic diversity in contracting populations, (b) identify different genetic diversity trajectories as a function of generation length, population density, and dispersal, and (c) highlight contrasting effects of contraction speeds on the patterns of genetic diversity at the end of contraction. Moreover, using the simulated microsatellite genetic data, we applied the VarEff method (Nikolic and Chevalet, 2014) to estimate effective population sizes across various contraction speeds and species-specific traits. We found that demographic histories inferred using these estimates can vary across species, even when they experience similar habitat changes.

## Materials and Methods

### Simulation tool

All the simulations were performed using the in-house SINS (Simulating Individuals in Space) program. An earlier version of SINS has been used in (Rasteiro et al., 2012) whereas the latest version has been used in (Sgarlata et al., 2022a,b). SINS allows the simulation of the joint spatio-temporal evolution of both genetic and demographic data over non-overlapping discrete generations. Specifically, SINS is a population-centered but individual-based software that uses a forward-time approach to simulate diploid individuals (males and females) over a two-dimensional grid of demes, similar to a 2D stepping stone lattice (Kimura and Weiss, 1964), and along general principles inspired by SPLATCHE (Currat et al., 2004). SINS and SPLATCHE (and the subsequent SPLATCHE versions 2, 3) (Currat et al., 2004; Ray et al., 2010; Currat et al., 2019) share several features but also differ in some properties as noted in (Rasteiro et al., 2012; Sgarlata et al., 2022a,b). For instance, SINS can simulate several types of molecular markers at the same time (mtDNA, Y chromosome) and can account for variance in reproductive success, unlike SPLATCHE or its later versions. More details on the comparison between the SINS and SPLATCHE implementations can be found in (Rasteiro et al., 2012; Sgarlata et al., 2022a,b).

In SINS, each deme is characterized by its carrying capacity ($K$) and friction ($F$) values. The population size in each deme evolves using the logistic growth formula (Smith and Slatkin, 1973), corrected to account for the fact that the growth is limited by the number of reproductive females (Rasteiro et al., 2012). Therefore, the expected number of diploid individuals in generation $t+1$ is given by

$$E(N_{t+1}) = 2N_{f,t} \frac{1+r}{1+r(\frac{2N_{f,t}}{K})} \quad (1)$$

where $N_{f,t}$ is the number of females at generation $t$, $r$ is the intrinsic growth rate per generation and $K$ is the deme carrying capacity. Note that deme population size $N_{t+1}$ is not computed deterministically as the number of individuals in generation $t+1$ is drawn from a Poisson distribution with mean $E(N_{t+1})$. We consider random mating within demes, meaning that one reproductive individual from each sex is randomly chosen to generate new offspring, which inherit one allele from each parent for each locus. At each generation, mating occurs iteratively until $N_{t+1}$ is reached.

Each focal deme can then exchange migrants at a certain rate, with up to four neighboring demes, depending upon the geographical location of the focal deme within the available habitat boundaries. Such migrations can only occur towards demes with suitable habitat (i.e. any demes having a non-null carrying capacity, $K>0$, and friction rate strictly lower than one, $F<1$).

The number of individuals migrating from the focal deme at each generation is drawn from a Poisson distribution with mean $M = \frac{N_t m n_d}{4}$, where $m$ is the migration rate per generation and $n_d$ is the number of available neighboring demes. The probability $P_{i,j}$ that one migrant from the focal
deme i moves to a neighboring deme j depends on the friction values of all neighboring demes: 144

$$P_{i,j} = \frac{1 - F_j}{\sum_{k \in \Gamma(i)} (1 - F_k)} \qquad (2)$$

where $\Gamma(i)$ denotes the set of demes in the neighborhood of deme i in the lattice. The number of migrants moving to each neighboring deme from a focal deme is sampled from a multinomial distribution with probabilities given by Eq. 2. While several types of molecular markers can be simulated, for the current study, we exclusively used independent microsatellite loci, mutating according to the Stepwise Mutation Model (Ohta and Kimura, 1973). 145 146 147 148 149

## Methodological details

We simulated a hypothetical landscape of 18 × 18 demes with two refugia of 2 × 2 demes each, inspired by Arenas et al. (2012) (for detailed comparison, see Supplementary S1, S2, Figure 3, Figure S2-S4). In Figure 1A, we illustrate the underlying spatial structure of our stepping-stone model, while Figure 1B depicts the expansion-contraction scenarios that were simulated onto the landscape. For clarity, we consider that the landscape is oriented north-south, and the two refugia (Figure 1A) were located in the southeastern and southwestern corners (hereafter referred to as SE refugium and SW refugium, respectively). The simulation began with an ancestral population at the top-left deme of the 2 × 2 SW refugium (Figure 1A), with one allele per microsatellite (considering 50 microsatellites in total, see below). The habitat then underwent three phases as described in Figure 1B:(i) expansion, where habitat grew northward from the SW refugium; (ii) stationary, when the habitat has filled the entire landscape; and (iii) contraction, during which the habitats reduced step by step until only the two refugia remained. During the expansion phase, the species colonizes new areas as a function of its intrinsic growth rate, migration rate and generation time. All demes with habitat were assumed to have the same carrying capacity (K) throughout the landscape. 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164

At the start of expansion, the mean number of alleles (hereafter MNA) was thus 1 and the expected heterozygosity (hereafter $H_e$) was 0. The size of the ancestral population was set to K, which was allowed to vary according to the simulated scenarios (Table S1). The total population size across the landscape reached the total number of demes (18 × 18 + 8) times the carrying capacity K when the entire landscape was initially colonized. After the contraction phase, it decreased to $(2 \times 4) \cdot K$ (where 4 corresponds to the number of demes in each of the SW and SE refugia, Figure S3). Altogether, this corresponds to an increase of more than (and a decrease of around) two orders of magnitude of the total population size. Once the whole landscape was colonized, the overall population size ranged from 16,600 when K = 50 to 332,000 when K = 1000 individuals. We used a growth rate of r = 0.8 per generation across all simulations. We simulated 50 independent microsatellite loci per individual, with a stepwise mutation rate of $10^{-3}$ per locus per generation. This rate aligns with estimates reported in the literature for several mammalian species (Schlötterer et al., 1998). 165 166 167 168 169 170 171 172 173 174 175 176 177

While the number of generations serves as a fundamental unit of time in population genetics, habitat changes are often described in years as they more accurately reflect the temporal scale of extrinsic environmental change. This choice of describing changes in habitats in years allowed us to examine the effects of the same habitat change for multiple co-existing species with varying species-specific traits. The total duration of the simulation (in years) was the same across range contraction scenarios and was set to ≈ 14,150 years. This duration was inspired by glacial-interglacial cycles, with interglacial periods typically lasting around 10,000-15,000 years (Berger and Loutre, 2010). The number of generations spanned by habitat changes thus differed between species, depending on their respective generation length (henceforth denoted by GL). 178 179 180 181 182 183 184 185 186

Habitats expanded over the entire landscape over 360 years, following a series of 18 expansion steps (Figure 1B, green period). For each northward expansion step, occurring every 20 years, we created a new horizontal line of 18 demes available for colonization. Each deme was assumed to cover an area of 100 × 100 km², resulting in a landscape that extended over 1,800 km from the southernmost 187 188 189 190
to the northernmost demes. The habitat expansion speed was set at approximately 5 km/year, which is higher but comparable to the estimated rate of spread observed in European pine, hazel, and alder (estimated to be 1.5 to 2 km/year according to Hewitt (1996)). Following the expansion phase, the habitat remained stable for 2,500 years before habitat contraction began. Throughout the stationary period, different species may colonize the landscape at varying rates, depending on their respective migration and growth rates. Except for one specific parameter combination ($GL = 10$ years, $K = 1,000$, $m = 0.003$), species with other parameter combinations colonized the entire landscape before habitat contraction began, as indicated in Table S1. To simulate habitat contraction, we set a horizontal line of 18 demes from the northern edge of the landscape as uninhabitable by setting $K = 0$ and $F = 1$. We define the onset of contraction as the time when the first line of habitats is emptied at time = 2,860 years (i.e 360 + 2500 years). This process was then repeated at fixed time intervals, considering three durations for each contraction step: 20, 240 or 640 years. This resulted in a total duration of range contraction (hereafter $t_{cont}$) of 340, 4,080, and 10,880 years, respectively (Figure 1B). These durations of contractions correspond to contraction speeds of approximately 5.3, 0.44, and 0.16 km/year. The longer the duration of the total contraction phase, the slower the contraction speed.

## Sampling design and statistics summarizing genetic diversity

To investigate the temporal dynamics of genetic diversity, we conducted 15 independent simulations for each of the scenarios depicted in Figure 1B. Specifically, we sampled 40 diploid individuals from the SW refugium, with 10 individuals selected from each of the four demes, across 60 time points. The timing of sampling (in years) remained the same for all species across scenarios with different contraction speeds. To assess the molecular diversity at the end of contraction with different speeds, we followed two sampling schemes: (i) we sampled 40 individuals from the SW refugium, selecting 10 individuals from each of its four demes, as done for the time series; and (ii) we sampled 20 individuals randomly across four demes each in both the SW and SE refugia.

We computed two summary statistics, averaged over 50 microsatellites: mean number of alleles (MNA), expected heterozygosity ($H_e$). We performed the computations of summary statistics using R software version 4.2.3 (Team et al., 2016). We calculated $H_e$ using hierfstat R package (v0.5.11) (Goudet, 2005), while MNA was computed using the adegenet R package (v2.1.10) (Jombart, 2008).

## Parameter selection for species: Allometric relationships

We employed allometric relationships and published data to generate biologically relevant parameter values for population density ($D$), mean dispersal distance ($\sigma$), and generation length ($GL$). Allometric relationships relate life-history traits with easily measurable traits, e.g. body mass (BM) (Peters and Peters, 1986; Blueweiss et al., 1978). Except for $GL$, population density, and dispersal distance can be described in the form of

$$Y = a(BM)^b \qquad (3)$$

where $Y$ represents the biological characteristic of interest, $BM$ represents the body mass of the species, and $a$ and $b$ are empirically derived constants.

Allometric relationships show that body mass negatively correlates with population density and positively with dispersal distance. We used the allometric relation reported for population density and body mass (Damuth, 1981):

$$D = 68.72(BM)^{-0.75}$$

where $D$ is population density in units of individual per km$^2$ and $BM$ is body mass in grams. For dispersal distance, we used a dataset linking empirical estimates of mean dispersal distance ($\sigma$)
to body mass across 164 mammal species (Santini et al., 2013). Since our simulations assume no bias in migration rate between sexes, we only used data points where dispersal distance was not sex-specific, or where the ratio in sex-specific dispersal (max(male, female)/min(male, female)) was less than 1.2. The resulting data points were analyzed using linear regression on a logarithmic scale obtaining the following relation

$$ \sigma = 0.05(BM)^{1.26} \qquad (4) $$

where σ is dispersal distance in km per generation and BM is body mass in grams. The data points between BM and σ are represented in Figure S6.

Given that no reported allometric relationship exists between GL and BM, we manually sampled nine combinations of body mass and generation length values from a public dataset on mammalian generation lengths (Pacifi et al., 2013). We considered species with body masses ranging from 0.1 kg to 100 kg because Equation 3 is only valid within this range (Silva and Downing, 1995). By visually inspecting the data represented in Figure S6A, we can distinguish two different trends:

* Species that vary in GL (1-20 years) but show little variation in BM (1-5 kg).

* Species that vary in BM (between 5-100 kg) but show little variation in GL (occurring mostly in the range of 4-10 years).

The first trend corresponds to species differing in generation lengths but which have similar population density and dispersal distances, while the second trend refers to species with similar generation lengths but varying in population density and dispersal distances. We selected three generation lengths (2, 5, 10 years) and five different body mass points (Table S1). For the chosen body masses, we used allometric relationships to compute corresponding population densities and dispersal distances.

## Conversion of dispersal distances to migration rates in a stepping-stone model

Dispersal distances of species are measured in continuous space and can be related to migration rates between demes in the stepping stone model via the following equation (Baird and Santos, 2010; Al-Asadi et al., 2019)

$$ \sigma^2 = mX^2 $$

where σ is the root mean squared axial parent-offspring distance (dispersal distance) per generation, m is the probability that an individual from any focal deme moves to neighboring demes (migration rate) per generation and X denotes the edge length of a deme (here set as 100 km). Thus, given the dispersal distance σ, we have $m = (\sigma/X)^2$.

## Panmictic model with time-varying population size

Although the former stepping-stone framework allowed us to vary the parameters related to species-level traits, we kept certain parameters constant due to the computational cost associated with spatial simulations, for instance, the duration of the stationary phases. To better understand the results obtained from spatial simulations, we proposed a simplified panmictic population model (Nei et al., 1975) which assumes similar changes in population sizes to the ones in the spatial simulation. This panmictic model allowed us to derive the expected heterozygosity ($H_e$) analytically. Forward in time, we started with a population of $N_i$ diploid individuals that were all homozygous (i.e. $H_e = 0$ or homozygosity, $F = 1$). This population immediately underwent an instantaneous population
expansion to $N_{exp}$ individuals. The population then remained stable for $t_{stat}$ generations till the start of the contraction phase. During contraction, population size reduced linearly over a duration of $t_{cont}$ generations, eventually reaching $N_i$ individuals (schematized in Figure 5A). The population size changes can be summarised as:

$$N(t) = \begin{cases} N_i & \text{if } t = 0 \\ N_{exp} & \text{if } 0 < t \le t_{stat} \\ N_{exp} - \left( \frac{N_{exp}-N_i}{t_{cont}} \right) (t - t_{stat}) & \text{if } t_{stat} < t \le t_{stat} + t_{cont} \\ N_i & \text{if } t > t_{stat} + t_{cont} \end{cases} \quad (5)$$

As previously mentioned, $N_i$ refers to initial population size (i.e., size of refuge), while $N_{exp}$ represents the expanded population size (i.e., size of the landscape). The term $t_{stat}$ indicates the duration (in generations) during which population remains at $N_{exp}$, and $t_{cont}$ denotes the period (in generations) over which $N_{exp}$ decreases linearly to $N_i$.

To find the analytical expression for the expected heterozygosity ($H_e$), we utilized homozygosity $F = 1 - H_e$. The recurrence formula for the homozygosity $F$, defined as the probability that two alleles, randomly chosen without replacement, are identical by descent at generation $t$, is given by:

$$F_t = (1-u)^2 \left\{ \frac{1}{2N_{t-1}} + \left(1 - \frac{1}{2N_{t-1}}\right) F_{t-1} \right\} \quad (6)$$

where $N_{t-1}$ represents the size of the population at time $(t-1)$ generation and $u$ is the mutation rate per generation.

At equilibrium, $F_t = F_{t-1}$, the expected homozygosity for a population of $N$ individuals, denoted by $F_{eq}(N)$, reads:

$$F_{eq}(N) = \frac{(1-u)^2}{(1-u)^2 + 2N(1-(1-u)^2)} \quad (7)$$

The above equation reduces to the known form $F_{eq}(N) = \frac{1}{1+4Nu}$ for $u \ll 1$. Approximating $F(t) - F(t-1) \approx \frac{dF}{dt}$, equation (6) reads:

$$\frac{dF}{dt} = (1-u)^2 F \left[ 1 - \frac{1}{2N} - \frac{1}{(1-u)^2} \right] + \frac{(1-u)^2}{2N} \quad (8)$$

The former differential equation can be solved analytically to derive an expression of $F$ over the stationary and contraction phases.

For the stationary period ($t > 0$ to $t \le t_{stat}$ generations), integrating equation (8) with $N = N_{exp}$ and initial condition $F(t = 0) = 1$, we get (for detailed steps, refer to Supplementary S3):

$$F_{stat}(t) = F_{eq}(N_{exp}) + [1 - F_{eq}(N_{exp})] e^{\left[ (1-u)^2 \left( 1 - \frac{1}{2N_{exp}} - \frac{1}{(1-u)^2} \right) \right] t} \quad (9)$$

where $F_{stat}(t)$ represents the homozygosity during time interval $0 < t \le t_{stat}$ generations and $F_{eq}(N_{exp})$ represents the equilibrium homozygosity for $N_{exp}$ individuals.

Similarly, we can obtain the expression of the homozygosity at any time $\Delta t_{cont} = t - t_{stat}$ during the contraction phase (denoted by $F_{cont}$). Specifically, setting $N = N_{exp} - \left( \frac{N_{exp}-N_i}{t_{cont}} \right) (t - t_{stat})$ and $F_{cont}(0) = F_{stat}(t_{stat})$ as the initial condition, we obtain (for detailed steps, refer to Supplementary S3):
$$F_{cont}(\Delta t_{cont}) = F_{stat}(t_{stat})e^{(a-1)\Delta t_{cont}} \left(1 - \frac{c}{b}\Delta t_{cont}\right)^{a/c} + e^{(a-1)(c\Delta t_{cont}-b)/c} \left(\frac{a}{c}\right) \left(\frac{(a-1)(c\Delta t_{cont}-b)}{c}\right)^{a/c} \\ \left\{ \Gamma\left(\frac{-a}{c}, \frac{(a-1)(c\Delta t_{cont}-b)}{c}\right) - \Gamma\left(\frac{-a}{c}, \frac{(1-a)b}{c}\right) \right\}, \quad (10)$$

where $\Gamma(s, x)$ denotes the incomplete Gamma function. Note that, we defined $\Delta t_{cont} = t - t_{stat}$, $a = (1 - u)^2$, $b = 2N_{exp}$, $c = 2\left(\frac{N_{exp}-N_i}{t_{cont}}\right)$ to ease notation.

Once contraction ends, population size remains constant ($N = N_i$). Therefore, integrating equation 8 with $N = N_i$ and initial condition $F_{cont}(t_{stat} + t_{cont})$ (homozygosity at the end of contraction), we get:

$$F_{cont-refuge}(\Delta t_{cont-refuge}) = F_{eq}(N_i) + [F_{cont}(t_{stat} + t_{cont}) - F_{eq}(N_i)] e^{(1-u)^2 \left(1 - \frac{1}{2N_i} - \frac{1}{(1-u)^2}\right) \Delta t_{cont-refuge}}, \quad (11)$$

where we have defined $\Delta t_{cont-refuge} = t - t_{stat} - t_{cont}$.

From the homozygosity values, we can express the temporal evolution of $H_e$ as:

$$H_e(t) = \begin{cases} 0 & \text{if } t = 0 \\ 1 - F_{stat}(t) & \text{if } 0 < t \le t_{stat} \\ 1 - F_{cont}(t - t_{stat}) & \text{if } t_{stat} < t \le t_{stat} + t_{cont} \\ 1 - F_{cont-refuge}(t - t_{stat} - t_{cont}) & \text{if } t > t_{stat} + t_{cont} \end{cases} \quad (12)$$

The expression of $H_e$ above is given in units of generation as time units. To account for variation in generation lengths among species, time should be rescaled as $\tilde{t} = \frac{t}{GL}$, where GL represents the generation length of the species.

## Results

### Temporal dynamics of genetic diversity in the refugium across different speeds of contraction

In Figure 1A, we illustrate the underlying spatial structure of our stepping-stone model, while Figure 1B depicts the expansion-contraction scenarios that were simulated onto the landscape. For most parameter combinations (eight out of nine) used here, the different species colonized the whole habitat before the contraction started. They thus experienced a stationary phase, where population size remained constant, before decreasing when habitat contracted. The three scenarios that we considered differed in the duration of the contraction phase ($t_{cont}$) and are represented in Figure 1B: i) $t_{cont}$ = 340 years, ii) $t_{cont}$ = 4,080 years, or iii) $t_{cont}$ = 10,880 years. Once the contraction ended, the different species survived in the two SW and SE refugia until the simulation was stopped at time 14,150 years.

Figure 2 depicts the evolution of $H_e$ in the SW refugium for two species with generation length $GL = 0.5$ (Figure 2A) and $GL = 5$ (Figure 2B) years respectively, having the same carrying capacity ($K = 50$) and migration rates ($m = 0.2$ per generation). For both species, $H_e$ increases during the expansion and stationary phases, though at a significantly faster pace for the species with shorter generation length. Then, $H_e$ trajectories between the two species depart significantly during the contraction phase.
The species with shorter *GL* exhibits persistent *H*<sub>e</sub> levels for the two long-duration contractions, plateauing before decreasing only towards the end of the contraction (Figure 2A). In contrast, species with a longer *GL* continue accumulating genetic diversity during the contraction phase, regardless of the contraction duration (Figure 2B). For instance, for *t*<sub>cont</sub> = 4,080 years (Figure 2B, blue line), the genetic diversity of the species with *GL* = 5 years keeps increasing for over 3,000 years (between time interval 2,860 – 6,000 years) after habitat had started contracting (Figure 2B, black line). Overall, the species with a shorter *GL* shows a quicker response (in both increase and decrease of its genetic diversity) to habitat changes throughout the simulation process. Regarding the impact of contraction speed on genetic diversity, we observed that for both species, at any time point following the start of contraction, faster contraction scenarios consistently led to lower genetic diversity than slower contraction scenarios. This is expected, since at any time point during the contraction phase, faster contractions imply smaller habitat ranges and population sizes than slower contractions.

We also assessed the impact of contraction speeds on the genetic diversity at the end of contraction, following the approach in Arenas et al. (2012). Figure 3 illustrates $H_e$ in the SW refugium population at four different time points (10, 60, 160, and 310 years), after the end of the contraction period for the three contraction speeds. For shorter GL (0.5 years, Figure 3A), the genetic diversity at the end of the fast contraction ($t_{cont}$ = 340 years, red bars) is consistently greater than for the slowest contraction speed ($t_{cont}$ = 10,880 years, green bars), regardless of the sampling times following the end of contraction. This observation aligns with the findings reported in Arenas et al. (2012), where faster contractions were interpreted to preserve higher levels of initial genetic diversity than slower contractions. While the difference is large shortly after the contraction ends (10 years post end of the contraction), it diminishes over time because all scenarios eventually reach the equilibrium of the refuge after contraction ends, regardless of contraction speed. However, for a species with longer GL (5 years, Figure 3B), the opposite pattern is observed, with a higher genetic diversity at the end of the slow contraction ($t_{cont}$ = 10,880 years, green bars) compared to the end of the fast contraction ($t_{cont}$ = 340 years, red bars). These results are thus in contrast with those of Figure 3A and remained consistent across different combinations of carrying capacities, migration rates, and sampling schemes (see Supplementary Table S3 and S4). Although we extensively described $H_e$, we also found consistent results on temporal dynamics and effect of contraction speed on genetic diversity when considering MNA (Figure S4, Supplementary Table S3, S4). Moreover, there appears to be a relationship between the temporal dynamics of genetic diversity during contraction phase and effect of contraction speed on genetic diversity, which is clarified in the discussion section. These results suggest that depending on their traits, some species may still accumulate genetic diversity during slow contraction periods (see next section).

These simulations show that variation in generation lengths alone can strongly impact the patterns of genetic diversity sampled in the refuge population during and at the end of different contraction speeds. In Figure 2, we observed that such effect seems to be related to the amount of genetic diversity accumulated in the SW refugium before the start of contraction (time = 2,860 years, black vertical line), with respect to the expected mutation-drift-migration equilibrium values (horizontal lines in Figure 2). These expected mutation-drift-migration equilibrium values were calculated using simulations, with details provided in Supplementary Table S2. Specifically, for the parameter combinations considered, these equilibria correspond to $H_e^{refuge} \simeq 0.36$ when the species is isolated in the SW refugium. Additionally, $H_e^{landscape} \simeq 0.77$ represents the value of $H_e$ measured in the SW refugium population when the entire landscape has been colonized for extended periods. In the above simulations, at the beginning of contraction (time = 2,860 years), the species with $GL = 0.5$ years have $H_e \simeq 0.73$ (close to $H_e^{landscape}$) whereas the species with $GL = 5$ years has $H_e \simeq 0.56$ (above $H_e^{refuge}$, below $H_e^{landscape}$).

The results obtained in this first set of simulations assumed that the two simulated species differed only in their generation length, all other demographic parameters remained the same (per generation), namely carrying capacity (K), growth rate (r), and migration rate (m). Yet, real species vary not only in their generation lengths but also in other species-specific traits such as population density and dispersal, and may thus exhibit patterns that differ from those reported in the present section. In the next section, we investigate the dynamics of genetic diversity for various combination
## Temporal genetic diversity with allometric relations across different speeds

We investigated the $H_e$ trajectories in our spatial framework considering three variable species-specific traits (generation length, dispersal distance, and population density) which are likely to influence local and global genetic diversity. Using empirical datasets, as detailed in the Methods section, we identified demographic parameters inspired by real species allometric relationships (see Supplementary Table S1 for parameter values). We found that for all considered combinations of traits values in Supplementary Table S1, the genetic diversity at the start of habitat contraction was always lower than their associated $H_e^{refuge}$ (Supplementary Table S1, S2), the $H_e$ value at equilibrium when the species is isolated in the SW refugium (Figure 4). These patterns thus differ from those described in the previous section, where $H_e$ before the start of contraction was greater than $H_e^{refuge}$ (Figure 2).

In Figure 4, we represented the results for two species with different generation lengths, carrying capacities, and migration rates. In Figure 4B, we considered a species with $GL = 5$ years and a body mass $BM \approx 76$ kg (leading to $K = 150$ and $m = 0.059$, which aligns with species like wild goats or bush pigs). We found that the amount of $H_e$ at the start of contraction (Figure 4B, time = 2,860 years, black vertical line, $H_e \approx 0.53$) is close to $H_e^{refuge}$ (horizontal dashed line). Then, under the slow contraction (green curve), $H_e$ continues to increase for 8,000 years during the contraction phase (between the time interval 4,000 – 12,000 years), before decreasing. During the contraction phase under the fast contraction (Figure 4B, red line), $H_e$ in the refugium increases monotonically towards $H_e^{refuge}$ and then plateaus after the end of the contraction. Surprisingly, despite a significant loss in population size ($\approx 98\%$ of the habitat), the trajectory of $H_e$ for fast contraction never decreases. In contrast, $H_e$ for slower contractions shows an initial increase above $H_e^{refuge}$, before decreasing during the contraction phase. Thus, for more realistic combinations of biological parameters, different contraction speeds also lead to distinct temporal trajectories of genetic diversity during the contraction phase, as observed in the previous section (Figure 2).

In Figure 4C, we considered a species with $GL = 10$ years and $BM \approx 6$ kg (leading to $K = 1,000$ and $m = 0.003$, characteristic of species like lion-tailed macaque). In this case, the amount of genetic diversity at the beginning of contraction (time = 2,860 years) is significantly lower than $H_e^{refuge}$. The duration of the contraction phase was such that the genetic diversity remained below $H_e^{refuge}$, regardless of the contraction speeds considered. We observed little to no differences in $H_e$ across contraction speeds (Figure 4C). This similarity in trajectories arises because genetic diversity in the refugium population is increasing towards $H_e^{refuge}$ and keeps on accumulating genetic diversity within the refugium, despite the fact that the species as a whole loses a large amount of population size across the landscape. The habitat changes outside of the refugium areas therefore have a negligible effect on the genetic diversity observed in the refugium.

## Panmictic model of range contraction: analytic solutions

The computational demands of spatial simulations constrained the amount of variable parameters that could be tested, thereby limiting the generality of conclusions that could be made regarding the genetic signatures of habitat change. To overcome this limitation and generalize our findings for a wider range of scenarios and species-level traits, we considered a simplified panmictic model where the temporal evolution of $H_e$ could be solved analytically (cf. Materials and Methods section). The panmictic model represents a single population experiencing similar but simplified dynamics of habitat change (here modeled as effective population size variations, Figure 5A). Specifically, we assumed an initial population of $N_i$ genetically homogeneous individuals ($H_e = 0$). This population underwent an instantaneous expansion to $N_{exp}$ individuals and remained stable for $t_{stat}$ generations. Afterward, the population size contracted back to $N_i$, linearly over $t_{cont}$ generations.

Figure 5B illustrates distinct temporal dynamics of genetic diversity during the contraction phase
428

for a stationary duration of 1,000 generations and a contraction phase of 800 generations. Varying $N_i$ and $N_{exp}$ influenced the mutation-drift equilibrium values ($H_{eq}^{N_i} = H_e^{refuge}$ and $H_{eq}^{N_{exp}} = H_e^{landscape}$), leading to different temporal dynamics of genetic diversity as the ones observed in the refugium population from the previous spatial simulations. In particular, our panmictic model could reproduce the same temporal dynamics of genetic diversity as the ones observed in the SW population from the previous spatial framework.

429

430

431

432

433

434

Lastly, we analyzed the impact of varying generation length on patterns of genetic diversity at the end of contraction by varying the stationary ($t_{stat}$) and contraction ($t_{cont}$) durations in our panmictic model, using equation 12 and assuming $N_i = 100$, $N_{exp} = 10^4$ and $u = 0.001$. We represent in Figure 5C the variations of $H_e$ at the end of contraction, as a function of the contraction and stationary times. Our results are qualitatively similar to those observed for the spatial simulations. For instance, considering a stationary phase of 1,250 years, this duration would translate into $t_{stat} = 2,500$ ($t_{stat} = 250$) generations for a species with $GL = 0.5$ ($GL = 5$) years before habitat contraction starts. These two species will experience different trajectories of $H_e$, qualitatively similar to those represented in the first two panels of Figure 5B.

435

436

437

438

439

440

441

442

443

Figure 5D shows how $H_e$ values differ at the end of contraction between both species, in a similar way comparable to Figure 3. For $t_{stat} = 2,500$ generations (short $GL = 0.5$ years, black horizontal line in Figure 5C and black curve in Figure 5D), there is enough time for $H_e$ to reach the mutation-drift equilibrium ($H_e^{N_{exp}}$) before the contraction starts. Consequently, genetic diversity at the end of contraction decreases as contraction duration increases. Thus, $H_e$ at the end of faster contractions is higher than at the end of slower ones. This finding aligns with the results obtained from spatial simulations conducted with the species having $GL = 0.5$ year (see Figure 3A) and with (Arenas et al., 2012). Conversely, for the species with long $GL = 5$ years, the population is far from the equilibrium value $H_e^{N_{exp}}$ at the start of habitat contraction. Therefore, for $t_{stat} = 250$ generations (grey horizontal line in Figure 5C and grey curve in Figure 5D), we observe that $H_e$ at the end of contraction increases as contraction time increases (slower contractions maintain more diversity) till $t_{cont} \approx 1,000$ generations. This is in agreement with the spatial simulation for $GL = 5$ years (Figure 3B) where we found higher genetic diversity at the end of slow contraction compared to fast. However, the panmictic model predicts that $H_e$ starts to decrease as the contraction time is further increased (Figure 5D, grey curve, when $t_{cont} > 1,000$ generations), resulting in an overall non-monotonous trajectory for $H_e$. The two types of relationships we observe for $H_e$ at the end of contraction and contraction speed are distinguished by the amount of genetic diversity accumulated before the start of the contraction phase. In the current case, the boundary between both behaviors corresponds to $t_{stat} \approx 1,220$ generations. In general, for $t_{stat} > 1,220$, the dependence is monotonous, whereas for $t_{stat} < 1,220$ generations the dependence is non-monotonous. For a stationary time of 1,220 generations, the panmictic population reaches $H_e = 0.896$, which is approximately $\approx 90\%$ of the equilibrium value of genetic diversity of the expanded population size ($N_{exp}$).

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

The non-monotonic trend between $H_e$ at the end of contraction and $t_{cont}$ identified under the panmictic model was not observed in the spatial simulations due to the limited duration of the range contractions that were tested. To overcome this limitation, we conducted additional spatial simulations, considering longer contraction times. Figure 5E confirms that $H_e$ eventually decreases for longer contraction times, giving rise to the non-monotonic trend predicted by the panmictic model. These results are further supported when analyzing the dependence of MNA at the end of the contraction, as represented in Supplementary Figure S7.

467

468

469

470

471

472

## Discussion

473

474

Natural populations undergo changes in population size and geographical distribution across various spatial and temporal scales. As a result, many species are likely to deviate from mutation-migration-drift equilibrium due to many (non-exclusive) factors, including fluctuations in habitat

475

476
477

478

479

480

481

482

483

484

485

486

487

488

489

490

availability, variation in species-specific traits, and dispersal patterns, as a consequence of geographical or environmental barriers. Here, we explored how species-specific traits can greatly influence our expectations of the genetic consequences of habitat changes among co-occurring species. We investigated the temporal dynamics of genetic diversity (as measured through the expected heterozygosity, $H_e$, and the mean number of alleles, MNA) in a refugium population subjected to a range expansion and subsequent contraction. We used spatially explicit individual-based simulations and analytical results from a single panmictic population model. Our analyses revealed that generation length significantly impacted the dynamics of genetic diversity. We explored various combination of parameters such as generation length, population density and dispersal, and identified three distinct temporal dynamics of $H_e$ and MNA in the refugium population during population contraction phase. These dynamics were: (i) an initial increase followed by a continuous decrease; or (ii) an increase that continues for a very long period after habitat started contracting, followed by a delayed decrease; or (iii) a continuous increase during the whole expansion and contraction period, as if habitat range had never changed beyond the refugium.

491

492

493

494

495

496

497

498

499

500

501

We also observed two contrasting patterns in how contraction speed affected the level of genetic diversity sampled at the end of contraction phase. The first pattern, consistent with findings of Arenas et al. (2012) and shown in Figure 3A, indicates that quick contractions maintain higher levels of diversity compared to slower contractions. The second pattern reveals a non-linear relationship between contraction speed and genetic diversity at the end of contraction. Specifically, genetic diversity initially increases as contraction speed decreases, up to a certain threshold speed. However, beyond this threshold, any further reduction in contraction speed leads to a decrease in genetic diversity (Figures 5D, E). In other words, fast and slow contraction speeds maintain less genetic diversity (at the end of the respective contraction phases) than intermediate speeds. Under such scenarios, observing high diversity in a refugium, after habitat contraction, does not necessarily reflect a recent habitat loss.

502

503

504

505

506

507

508

509

510

Altogether, our spatial and panmictic analysis allowed us to identify different temporal trajectories between population size, genetic diversity, and contemporary habitat extent, suggesting that accurately inferring habitat extent or population size from genetic data can be more challenging. More specifically, our results suggest that the increasing use of genetic data to reconstruct the demographic history of species (interpreted as changes in effective population size, $N_e$) may not be straightforwardly connected to ancient habitats and population sizes. Our work also emphasizes the importance and relevance of demographic inferences based on simulations, to better account for the impact of species-level traits on patterns of genetic evolution. In the next sections, we try to discuss and clarify some consequences of our work.

511

## Temporal trajectories of genetic diversity in contracting populations

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

Population genetics theory states that a decrease in effective population size ($N_e$) leads to a reduction in genetic diversity in a panmictic model ($H_e$ and MNA). Under simplifying assumptions, $N_e$ is also thought to be positively correlated with census population size ($N_c$), which is itself correlated to the geographical range, and possibly to habitat extent (if no major barriers to gene flow are present in the range). Thus, a decline in $N_c$ or habitat range is expected to lead to a reduction in genetic diversity, and to a decrease in $N_e$. These predictions find solid support in numerous empirical studies that have reported a loss of genetic diversity in species undergoing habitat loss (Nyström et al., 2006; Dures et al., 2019; Taron et al., 2021; Gauthier et al., 2020). However, several recent studies have observed no apparent reduction in genetic diversity despite drastic decreases in population sizes (Maebe et al., 2016b; Le Gouar et al., 2009; Fourcade et al., 2020; Carvalho et al., 2019), and in some cases, even an increase in genetic diversity over time has been documented (Townsend et al., 2023). These observations have been attributed to either the immigration of individuals from neighboring populations compensating for the diversity loss due to population decline or to species-specific traits, such as generation time, which could delay the genetic response (Goossens et al., 2005; Hailer et al., 2006). Our results show that $H_e$ trajectories exhibiting opposite trends with respect to population size and habitat trajectories are expected and may happen under reasonable conditions. For instance,
528

529

530

531

532

533

534

535

536

we found that populations that deviate from mutation-drift-migration equilibrium at the start of a
contraction phase can present these unexpected temporal dynamics of increasing genetic diversity
over hundreds or thousands of years in the refugium population when habitats and population range
are contracting. Importantly, we found that genetic diversity trajectories during population decline
can be predicted by comparing the genetic diversity of the population (e.g. $H_e$) with the genetic
diversity expected under expanded and contracted habitat ranges (i.e. $H_e^{\text{landscape}}$ and $H_e^{\text{refuge}}$,
respectively). Specifically, the temporal trajectory of $H_e$ in the refugium population depends on
which of the following three situations is observed at the start of contraction: (i) $H_e \approx H_e^{\text{landscape}}$
(ii) $H_e^{\text{landscape}} > H_e > H_e^{\text{refuge}}$ (iii) $H_e < H_e^{\text{refuge}}$.

537

538

539

540

541

542

543

544

545

In the first regime ($H_e \approx H_e^{\text{landscape}}$, as in Figures 2A, 5B1), the whole population is at its mutation-drift-migration equilibrium before contraction starts. Any decline in population size thus results in a decrease in genetic diversity eventually reaching the equilibrium of the refugium population size ($H_e^{\text{refuge}}$). If the population decline occurs at different speeds, genetic diversity at the end of faster contraction tends to be higher than that at the end of slower contractions. This occurs because, when contractions occur quickly, a substantial amount of genetic diversity is retained at the end of contraction, even though this diversity is lost rapidly afterwards. On the contrary, during slower contractions, the refugium population is subjected to more genetic drift, leading to a greater and gradual loss of genetic diversity by the end of contraction.

546

547

548

549

550

551

552

553

554

For the second regime ($H_e^{\text{landscape}} > H_e > H_e^{\text{refuge}}$, when contraction starts), as in Figure 2B), we find that during contractions, genetic diversity can keep increasing for rather long periods during contraction phase, before decreasing. A larger difference between $H_e^{\text{landscape}}$ and $H_e$, at the beginning of the contraction phase, is associated with a longer period during which genetic diversity increases before starting to decrease. Interestingly, with this temporal dynamics during contraction, fast range contractions do not preserve higher levels of genetic diversity. Instead, intermediate contraction speeds preserve a larger amount of genetic diversity compared to fast and slow contraction speeds (see Figure 5D, E, Figure S7). This means that, under this regime, low levels of genetic diversity in a refugium at present could be either due to a recent or ancient contraction.

555

556

557

558

559

560

Lastly, in the third regime ($H_e < H_e^{\text{refuge}}$, when the contraction starts), $H_e$ will increase until it eventually reaches $H_e^{\text{refuge}}$ and, if by that time the contraction is finished, $H_e$ will stay at this level until further habitat change occurs. Under such a scenario, we will thus observe an increase of diversity in refugium populations through time despite a major decrease in census population size and habitat across the landscape. However, if the contraction is slow, such that at some time point during contraction $H_e > H_e^{\text{refuge}}$, then the temporal dynamics will be similar to the second regime.

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

The identification of genetic diversity trajectories during population contraction is crucial for conservation efforts. Our results show that these trajectories could be identified by calculation of three quantities, $H_e$, $H_e^{\text{refuge}}$, and $H_e^{\text{landscape}}$, and knowing fundamental species-level traits. Identifying temporal dynamics also contribute to the concept of "genetic extinction debt", highlighting time lags where the effects of habitat disturbances on genetic diversity may not be immediately apparent, but instead manifest over time (Epps and Keyghobadi, 2015; Gargiulo et al., 2024). Moreover, temporal dynamics of genetic diversity during contraction phase determines how different contraction speeds impact genetic diversity at the end of contraction. We observed these patterns in both spatially structured (using simulation) and simple panmictic populations (analytically), though they differ in several important ways, as discussed below. Our findings suggest that the widespread belief that fast (or recent) contractions better preserve initial genetic diversity than slow (long-term decline) (Arenas et al., 2012) may not always hold true. For instance, Wang et al. (2022) concluded that high genetic diversity of Yellow-breasted Bunting was due to a recent rapid population decline. However, our results indicate that this expectation holds true if the populations were at equilibrium before the contraction occured. Furthermore, co-distributed species living in the same habitats are often expected to exhibit similar demographic responses to habitat changes. Consequently, it is anticipated that co-distributed species will show congruence in their trajectories of genetic diversity and, as a result, in their inferred effective population size ($N_e$). Our results suggest that even closely related species undergoing similar habitat changes may exhibit divergent temporal dynamics of genetic diversity due to variations in the initial levels of genetic diversity before the habitat changes. As a
581
582

consequence, this questions our ability to link $N_e$ trajectories inferred from genetic data to ancient habitat changes, as we will discuss in the next section.

583

## Comparison of panmictic vs. structured populations

584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600

The field of population genetics has historically relied on simple models to capture the properties of natural populations. The panmictic model, widely used for its simplicity and mathematical tractability, has proven valuable in numerous scenarios (Tang et al., 2024). We found that it was useful here too, since it allowed us to explore a large spectrum of parameters, identifying $H_e$ trajectories that would have otherwise been difficult and computationally intensive to explore with spatial simulations (e.g., non-monotonic relationship between contraction duration and $H_e$ estimated at the end of the contraction). In this study, the panmictic and spatially structured population approaches yielded consistent qualitative results regarding the temporal dynamics of $H_e$ and the influence of contraction speed on $H_e$ at the end of contraction. However, natural populations are generally spatially structured and typically require a larger number of generations to reach mutation-drift-migration equilibrium compared to panmictic populations (Leblois et al., 2006; Epps and Keyghobadi, 2015). As a result, spatial and panmictic populations may exhibit different temporal dynamics of genetic diversity for given population size changes (Figure S5). Our panmictic model was both useful and limited since it could not include the effect of migration, even though the latter plays a major role in the temporal decoupling of the different trajectories ($N_c$, $H_e$, $N_e$). The differences between panmictic and spatial models become particularly relevant for the demographic inference, which we discuss in next section.

601
602

## Coupling between genetic diversity, species-specific traits, habitat change and demographic history

603
604
605
606
607
608
609
610
611
612
613
614

While species genetic diversity and spatial distribution are influenced by habitat changes over time (Hewitt, 2004; Svenning et al., 2015), we rarely have a direct measure of these changes. In the past few decades, population geneticists have increasingly relied on genetic data and advanced statistical methods to reconstruct the demographic history of populations or species (Salmona et al., 2019; Marchi et al., 2021). Coalescent theory provides a solid theoretical framework to link genomic diversity with past demography, leading to development of numerous tools for demographic inference (Spence et al., 2018; Nadachowska-Brzyska et al., 2022). Many of these tools assume that genetic data are sampled from a neutrally evolving panmictic population. This simplification makes demographic inference possible since only one parameter is allowed to change ($N_e$). Consequently, the "demographic history" is often implicitly or explicitly interpreted in terms of changes in $N_e$, which are then related to putative changes in climate or habitat (Nadachowska-Brzyska et al., 2015; Natesh et al., 2020; Teixeira et al., 2021).

615
616
617
618
619
620
621
622
623
624
625
626
627
628

However, it is increasingly evident that interpreting the past demography of a species using reconstructed $N_e$ trajectories under a simplistic panmictic model can lead to misinterpretations and biases (Loog, 2021). This problem is particularly pronounced when the populations have evolved under more complex scenarios, such as interconnected populations (or structured population) (Wakeley, 1999; Beaumont, 2004; Chikhi et al., 2010; Mazet et al., 2016) or experiencing selective pressures (Schrider et al., 2016; Johri et al., 2021; Boitard et al., 2022). An increasing number of studies have shown that statistical methods that infer $N_e$ while ignoring population structure may infer, quantify, and date with great precision changes in $N_e$ that never occurred (Wakeley, 1999; Chikhi et al., 2010; Heller et al., 2013; Chikhi et al., 2018). This is because structured population models have genealogical properties that can be interpreted erroneously under panmictic models (Herbots, 1994; Notohara, 1990; Wakeley, 1999; Mazet et al., 2016). It is important to note that the idea that inferred $N_e$ trajectories do not necessarily reflect changes in census population size (and by extension, in habitat range) has been around for some time (Leonardi et al., 2021; Miller et al., 2021).
629

630

631

632

633

634

635

636

637

638

Irrespective of the way in which $N_e$ is interpreted, co-distributed species living in the same habitat are often expected to exhibit similar demographic responses to climatic changes. As a result, it is anticipated that co-distributed species will show congruence in their genetic diversity trajectories and, consequently, in their inferred $N_e$. However, empirical studies suggest that such congruence is not necessarily observed. For example, studies inferring $N_e$ across various taxa have reported discordant or asynchronous demographic responses to Pleistocene climate changes (Burbrink et al., 2016; Kuhn et al., 2022; Walton et al., 2021; Bai et al., 2018). Additionally, even closely related species sharing similar ecological or life-history traits may exhibit incongruent demographic histories in similar environment. This has been observed across diverse taxonomic groups such as in trees (Bai et al., 2018), hymenopteran parasitoids (Walton et al., 2021), herpetofauna (Kuhn et al., 2022).

639

640

641

642

643

644

645

646

647

In our simulated scenarios, we tested whether coalescent-based approaches to infer $N_e$, such as the method of Nikolic and Chevalet (2014) implemented in the R package VarEff would help us better elucidate the relationships between $N_c$, $H_e$, and inferred $N_e$ trajectories (for details, see Supplementary S6). This method assumes a panmictic population model, utilizes microsatellites to infer changes in $N_e$ (Figures S11, S12), and has been applied to various endangered species (Salmona et al., 2017; Srinivas and Jhala, 2023). Our goal was to assess if simulated scenarios, involving species with varying traits undergoing similar changes in habitat, would exhibit comparable inferred changes in $N_e$. Additionally, we explored the ability to distinguish fast and slow contractions using inferred $N_e$ trajectories.

648

649

650

651

652

653

654

655

656

657

658

659

660

661

662

663

664

We conducted the analysis using genetic data sampled at the end of simulation, which corresponds to 14,000 years forward in time. We examined various contraction duration and species-specific parameters (for details refer to Supplementary S6). In scenarios with fast contractions, populations were smaller (confined to refuge area) for extended periods in the past before sampling. Conversely, for slower contractions, populations decreased gradually over time before their sampling. For each scenario, we plotted the temporal trajectories of $N_c$, $H_e$, and inferred $N_e$. Our results suggest that the inferred $N_e$ trajectories are not generally consistent with the simulated $N_c$ changes, and can have variable similarity with $H_e$ trajectories (Figures S9, S10). Fast and slow contraction scenarios showed distinguishable inferred $N_e$ trajectories (Figures S9, S10, S13). For instance, for three scenarios, the slow contraction exhibits higher inferred $N_e$ compared to the fast contraction scenarios in the past (Figures S9b, c, d). However, we note that in Figure S9a, the inferred $N_e$ of fast contraction was surprisingly higher than for the slower contraction in the past. This is surprising, as fast contraction experiences a greater loss in genetic diversity while population size remained limited to the SW refugium for a longer duration till the present. We further provide a detailed discussion in Supplementary S6 (Figures S14, S15, S16). In particular, it appears that this output of the VarEff method, which is used to represent changes in $N_e$ should be interpreted with care as the posterior distribution may be rather wide (Supplementary Figure S16).

665

666

667

668

669

670

671

672

673

674

675

676

677

678

We stress that these results are admittedly limited to a few scenarios. It is possible that the $N_e$ trajectories could represent the "demographic history" of species under other scenarios, yet to be identified. However, in our scenarios, these $N_e$ trajectories should be taken with a grain of salt, as they provide little information about past $N_c$ or $H_e$ trajectories or habitat extent. The observed discordance in inferred $N_e$ across multi-species studies in similar environments should not surprising. Co-distributed species likely differ in their ecology and species-specific traits relevant to genetics, and applying the same demogenomic inference methodology is unlikely to provide the same trajectories across species, even if they underwent similar climatic and habitat changes. These inconsistencies are further supported by empirical studies (Bai et al., 2018; Walton et al., 2021; Kuhn et al., 2022). For example, a study by Miller et al. (2021) used a Bayesian Skyline Plot (BSP) to infer $N_e$ and employed Species Distribution Modelling (SDM) to reconstruct the geographical range of 102 Holarctic bird species. The study revealed no correlation between the magnitude of inferred changes in $N_e$ and range size inferred using species distribution modeling and past climate reconstructions.

679

680

681

Our work suggests that the trajectories of $N_c$, habitat size, $H_e$, and inferred $N_e$ cannot be expected to be coupled temporally in a simple way. Instead, they represent different conceptual objects and $N_e$ trajectories should be seen as estimates of complex summary statistics, rather than
682
683
684
685
686

objects that can be interpreted directly (Chikhi et al., 2018, 2024). Hence, the reconstruction of demographic history should be interpreted with extreme caution and be validated with realistic simulations and external evidence from other disciplines, when available. Spatially explicit models remain essential for assessing the genetic effects of habitat loss and fragmentation, over both short and long time scales.

687

## Perspectives

688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712

Our results suggest that more work is needed to improve our understanding of spatially struc-
tured models undergoing cycles of expansion and contraction. We found connections between some
species-specific traits, $H_e$, and $N_c$ dynamics, and identified potentially serious issues in the use of
$N_e$ as a proxy for any of these. We showed how comparing current $H_e$ to $H_e^{refuge}$ and $H_e^{landscape}$
can be used to predict the future trajectory of $H_e$ with important implications for the conserva-
tion of endangered species currently experiencing habitat contraction (within or outside refugia).
Although we simulated different species by taking into account empirical allometric relationships,
we still made many simplifying assumptions. For instance, we considered homogeneous carrying
capacities, migration, mutation, and growth rates. We also assumed random mating within demes
for all species across the landscape and considered co-distributed species as ecologically independent,
with no consideration for species interactions. However, deviations from these factors could impact
the time to mutation-migration-drift equilibrium and would affect the temporal scale of the genetic
responses we identified here (Epps and Keyghobadi, 2015). Nevertheless, testing different mutation
models gave qualitatively similar effect of contraction speed on $H_e$ at the end of contraction phase
(Supplementary Figure S8). Additionally, we simulated population contraction by gradually remov-
ing a horizontal layer of demes from north to south. While species might preferentially migrate
towards refugia in response to contraction, Arenas et al. (2012) found that qualitative effects of
contraction speed on genetic diversity remained similar between isotropic and anisotropic migration
towards refugia. Among many other factors, social structure is one that we consider to be relatively
understudied within spatial models (Chesser, 1991a,b; Sugg et al., 1996; Parreira and Chikhi, 2015).
We also assumed that the connectivity between the refugium and the neighboring populations re-
mained similar throughout the simulation. Also, the species we simulated have non-overlapping
generations but theoretical studies show that overlapping generations increase the time to reach
equilibrium (Lloyd et al., 2013). Incorporating non-overlapping generations could increase the lag
effect identified here, disconnecting, even more, the temporal trajectories of habitat, $N_c$ and $H_e$.

713
714
715
716
717
718
719
720
721
722
723
724
725
726
727

Regarding our simulated scenarios, we considered specific initial conditions where all individuals
were homozygous and we focused on a single cycle of range expansion and contraction, although
most species likely experienced multiple glacial cycles during the Holocene and Pleistocene. Our
results indicate that the latter two limitations should not impact the main conclusions of our work,
because we can still predict the behavior of the trajectories as a function of the amount of genetic
diversity accumulated before a new cycle of expansion and contraction, forward in time. Yet, an
interesting perspective would be to investigate the impact of multiple cycles while accounting for
multiple refugia scenarios. Indeed, independent refugia can maintain divergent alleles or haplotypes
for longer periods of cyclic expansions and contractions, which can complexify the trajectories of $H_e$
and inferred $N_e$ when these alleles mix during population reconnections. A natural extension of our
work will be to explore the temporal dynamics of genetic diversity in samples from other regions
within the landscape as the dynamics of their genetic diversity could be different from the ones
observed within the present-day refugium. Future research should also investigate more complex
patterns of range contraction (Rogan et al., 2023) and behaviors in continuous space (Leblois et al.,
2006; Bradburd and Ralph, 2019).
Acknowledgements
728

We thank all the members of the Population and Conservation Genetic group at Instituto Gulbenkian de Ciência (IGC) for their support and useful discussions on this topic. We would also like to acknowledge the Bioinformatics Unit and the Informatics Team of the IGC for their help and support with computational resources. We acknowledge the financial support received from the Fundação para a Ciência e Tecnologia through a PhD fellowship awarded to RV (reference number: UI/BD/154372/2022), a grant awarded to GMS (PD/BD/114343/2016), and LC (PTDC-BIA-EVL/30815/2017), as well as the financial support provided by the Calouste Gulbenkian Foundation through the PONTE program to DSP. LC's research was also supported by the 2015-2016 Bio-divERsA COFUND call for research proposals, with national funders ANR (ANR-16-EBI3-0014) and the Fundação para a Ciência e Tecnologia (reference: Biodiversa/0003/2015 and PT-DLR 01LC1617A). LC additionally received support from the DevOCGen project, funded by the Oc-citanie Regional Council's "Key challenges BiodivOc" program. This work was also supported by the LABEX entitled TULIP (ANR-10-558 LABX-41 and ANR-11-IDEX-0002-02) as well as the IRP BEEG-B (International Research Project Bioinformatics, Ecology, Evolution, Genomics and Be-haviour). We acknowledge an Investissement d'Avenir grant of the Agence Nationale de la Recherche (CEBA: ANR-10-LABX-25-01). RV would also like to acknowledge Victor Mello for support with Figure 6 and for valuable discussions, and Pragya Kaul for her help in reviewing the manuscript. He also extends his thanks to Dhiraj Kumar, Navneet Kumar Karn and Dr. Rahul Kumar (Assistant Professor at Sheodeni Sao College, Bihar) for reviewing Hindi translation of the abstract.
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747

References
748

Al-Asadi, H., Petkova, D., Stephens, M., and Novembre, J. (2019). Estimating recent migration and population-size surfaces. *PLoS genetics*, 15(1):e1007908.

Alcala, N. and Vuilleumier, S. (2014). Turnover and accumulation of genetic diversity across large time-scale cycles of isolation and connection of populations. *Proceedings of the Royal Society B: Biological Sciences*, 281(1794):20141369.

Allendorf, F. W. (1986). Genetic drift and the loss of alleles versus heterozygosity. *Zoo biology*, 5(2):181-190.

Arenas, M., Ray, N., Currat, M., and Excoffier, L. (2012). Consequences of range contractions and range shifts on molecular diversity. *Molecular biology and evolution*, 29(1):207-218.

Bai, W.-N., Yan, P.-C., Zhang, B.-W., Woeste, K. E., Lin, K., and Zhang, D.-Y. (2018). Demographically idiosyncratic responses to climate change and rapid pleistocene diversification of the walnut genus juglans (juglandaceae) revealed by whole-genome sequences. *New Phytologist*, 217(4):1726-1736.

Baird, S. J. and Santos, F. (2010). Monte carlo integration over stepping stone models for spatial genetic inference using approximate bayesian computation. *Molecular ecology resources*, 10(5):873-885.

Beaumont, M. A. (2004). Recent developments in genetic data analysis: what can they tell us about human demographic history? *Heredity*, 92(5):365-379.

Berger, A. and Loutre, M.-F. (2010). Modeling the 100-kyr glacial-interglacial cycles. *Global and Planetary Change*, 72(4):275-281.

Beyer, R. M. and Manica, A. (2020). Historical and projected future range sizes of the world's mammals, birds, and amphibians. *Nature Communications*, 11(1):5633.

Blueweiss, L., Fox, H., Kudzma, V., Nakashima, D., Peters, R., and Sams, S. (1978). Relationships between body size and some life history parameters. *Oecologia*, pages 257-272.
Boitard, S., Arredondo, A., Chikhi, L., and Mazet, O. (2022). Heterogeneity in effective size across the genome: effects on the inverse instantaneous coalescence rate (iicr) and implications for demographic inference under linked selection. *Genetics*, 220(3):iyac008.

Bradburd, G. S. and Ralph, P. L. (2019). Spatial population genetics: it's about time. *Annual Review of Ecology, Evolution, and Systematics*, 50:427-449.

Burbrink, F. T., Chan, Y. L., Myers, E. A., Ruane, S., Smith, B. T., and Hickerson, M. J. (2016). Asynchronous demographic responses to pleistocene climate change in eastern nearctic vertebrates. *Ecology Letters*, 19(12):1457-1467.

Cabe, P. R. (1998). The effects of founding bottlenecks on genetic variation in the european starling (Sturnus vulgaris) in north america. *Heredity*, 80(4):519-525.

Carvalho, C. S., Lanes, É. C., Silva, A. R., Caldeira, C. F., Carvalho-Filho, N., Gastauer, M., Imperatriz-Fonseca, V. L., Nascimento Júnior, W., Oliveira, G., Siqueira, J. O., et al. (2019). Habitat loss does not always entail negative genetic consequences. *Frontiers in Genetics*, 10:1011.

Ceballos, G., Ehrlich, P. R., and Dirzo, R. (2017). Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines. *Proceedings of the national academy of sciences*, 114(30):E6089-E6096.

Chesser, R. K. (1991a). Gene Diversity and Female Philopatry. *Genetics*, 127:437-447.

Chesser, R. K. (1991b). Influence of gene flow and breeding tactics on gene diversity within populations. *Genetics*, 129(2):573-583.

Chikhi, L., Rodríguez, W., Grusea, S., Santos, P., Boitard, S., and Mazet, O. (2018). The iicr (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice. *Heredity*, 120(1):13-24.

Chikhi, L., Rodriguez, W., Paris, C., Ha-Shan, M., Jouniaux, A., Arredondo, A., Nous, C., Grusea, S., Corujo, J., Lourenco, I., et al. (2024). Extending the iicr to multiple genomes and identification of limitations of some demographic inferential methods. *bioRxiv*, pages 2024-08.

Chikhi, L., Sousa, V. C., Luisi, P., Goossens, B., and Beaumont, M. A. (2010). The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes. *Genetics*, 186(3):983-995.

Crates, R., Olah, G., Adamski, M., Aitken, N., Banks, S., Ingwersen, D., Ranjard, L., Rayner, L., Stojanovic, D., Suchan, T., et al. (2019). Genomic impact of severe population decline in a nomadic songbird. *PLoS One*, 14(10):e0223953.

Currat, M., Arenas, M., Quilodran, C. S., Excoffier, L., and Ray, N. (2019). Splatche3: simulation of serial genetic data under spatially explicit evolutionary scenarios including long-distance dispersal. *Bioinformatics*, 35(21):4480-4483.

Currat, M., Ray, N., and Excoffier, L. (2004). Splatche: a program to simulate genetic diversity taking into account environmental heterogeneity. *Molecular Ecology Notes*, 4(1):139-142.

Cushman, S. A. (2006). Effects of habitat loss and fragmentation on amphibians: a review and prospectus. *Biological conservation*, 128(2):231-240.

Damuth, J. (1981). Population density and body size in mammals. *Nature*, 290:699-700.

De Kort, H., Prunier, J. G., Ducatez, S., Honnay, O., Baguette, M., Stevens, V. M., and Blanchet, S. (2021). Life history, climate and biogeography interactively affect worldwide genetic diversity of plant and animal populations. *Nature communications*, 12(1):516.
815

816

817

Dures, S. G., Carbone, C., Loveridge, A. J., Maude, G., Midlane, N., Aschenborn, O., and Gottelli, D. (2019). A century of decline: Loss of genetic diversity in a southern african lion-conservation stronghold. *Diversity and distributions*, 25(6):870–879.

818

819

Edmonds, C. A., Lillie, A. S., and Cavalli-Sforza, L. L. (2004). Mutations arising in the wave front of an expanding population. *Proceedings of the National Academy of Sciences*, 101(4):975–979.

820

821

Epps, C. W. and Keyghobadi, N. (2015). Landscape genetics in a changing world: disentangling historical and contemporary influences and inferring change. *Molecular ecology*, 24(24):6021–6040.

822

823

Finn, C., Grattarola, F., and Pincheira-Donoso, D. (2023). More losers than winners: investigating anthropocene defaunation through the diversity of population trends. *Biological Reviews*.

824

825

826

Fourcade, Y., Richardson, D. S., and Secondi, J. (2020). No evidence for a loss of genetic diversity despite a strong decline in size of a european population of the corncrake crex crex. *Bird conservation international*, 30(2):260–266.

827

Frankham, R. (2005). Genetics and extinction. *Biological conservation*, 126(2):131–140.

828

829

Gargiulo, R., Budde, K. B., and Heuertz, M. (2024). Mind the lag: understanding delayed genetic erosion.

830

831

Garnier, J. and Lafontaine, P. (2022). Life history traits and dispersal shape neutral genetic diversity in metapopulations. *Journal of Mathematical Biology*, 84(6):45.

832

833

834

Gauthier, J., Pajkovic, M., Neuenschwander, S., Kaila, L., Schmid, S., Orlando, L., and Alvarez, N. (2020). Museomics identifies genetic erosion in two butterfly species across the 20th century in finland. *Molecular Ecology Resources*, 20(5):1191–1205.

835

836

837

Gonçalves, M., Siegismund, H. R., van Vuuren, B. J., Ferrand, N., and Godinho, R. (2021). Evolutionary history of the roan antelope across its african range. *Journal of Biogeography*, 48(11):2812–2827.

838

839

840

841

Goossens, B., Chikhi, L., Jalil, M., Ancrenaz, M., Lackman-Ancerenaz, I., Mohamed, M., Andau, P., and Bruford, M. W. (2005). Patterns of genetic diversity and migration in increasingly fragmented and declining orang-utan (pongo pygmaeus) populations from sabah, Malaysia. *Molecular Ecology*, 14(2):441–456.

842

843

Goudet, J. (2005). Hierfstat, a package for r to compute and test hierarchical f-statistics. *Molecular ecology notes*, 5(1):184–186.

844

845

Grant, W. S. and Cheng, W. (2012). Incorporating deep and shallow components of genetic structure into the management of alaskan red king crab. *Evolutionary Applications*, 5(8):820–837.

846

847

848

Grant, W. S., Liu, M., Gao, T., and Yanagimoto, T. (2012). Limits of bayesian skyline plot analysis of mtDNA sequences to infer historical demographics in pacific herring (and other species). *Molecular phylogenetics and evolution*, 65(1):203–212.

849

850

851

Haddad, N. M., Brudvig, L. A., Clobert, J., Davies, K. F., Gonzalez, A., Holt, R. D., Lovejoy, T. E., Sexton, J. O., Austin, M. P., Collins, C. D., et al. (2015). Habitat fragmentation and its lasting impact on earth's ecosystems. *Science advances*, 1(2):e1500052.

852

853

854

855

Hailer, F., Helander, B., Folkestad, A. O., Ganusevich, S. A., Garstad, S., Hauff, P., Koren, C., Nygård, T., Volke, V., Vila, C., et al. (2006). Bottlenecked but long-lived: high genetic diversity retained in white-tailed eagles upon recovery from population decline. *Biology letters*, 2(2):316–319.

856

857

Heller, R., Chikhi, L., and Siegismund, H. R. (2013). The confounding effect of population structure on bayesian skyline plot inferences of demographic history. *PloS one*, 8(5):e62992.
858
859
Herbots, H. M. J. D. (1994). *Stochastic models in population genetics: genealogy and genetic differentiation in structured populations*. PhD thesis, Queen Mary University of London.

860
861
Hewitt, G. M. (1996). Some genetic consequences of ice ages, and their role in divergence and speciation. *Biological journal of the Linnean Society*, 58(3):247-276.

862
863
Hewitt, G. M. (2004). Genetic consequences of climatic oscillations in the quaternary. *Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences*, 359(1442):183-195.

864
865
866
867
Hoban, S., Campbell, C. D., da Silva, J. M., Ekblom, R., Funk, W. C., Garner, B. A., Godoy, J. A., Kershaw, F., MacDonald, A. J., Mergeay, J., et al. (2021). Genetic diversity is considered important but interpreted narrowly in country reports to the convention on biological diversity: Current actions and indicators are insufficient. *Biological Conservation*, 261:109233.

868
869
Hudson, R. R. (1983). Properties of a neutral allele model with intragenic recombination. *Theoretical population biology*, 23(2):183-201.

870
871
872
Johri, P., Riall, K., Becher, H., Excoffier, L., Charlesworth, B., and Jensen, J. D. (2021). The impact of purifying and background selection on the inference of population history: problems and prospects. *Molecular biology and evolution*, 38(7):2986-3003.

873
874
Jombart, T. (2008). adegenet: a r package for the multivariate analysis of genetic markers. *Bioinformatics*, 24(11):1403-1405.

875
876
Kimura, M. and Weiss, G. H. (1964). The stepping stone model of population structure and the decrease of genetic correlation with distance. *Genetics*, 49(4):561.

877
Kingman, J. F. C. (1982). The coalescent. *Stochastic processes and their applications*, 13(3):235-248.

878
879
Klopfstein, S., Currat, M., and Excoffier, L. (2006). The fate of mutations surfing on the wave of a range expansion. *Molecular biology and evolution*, 23(3):482-490.

880
881
882
Kuhn, A., Gehara, M., Andrianarimalala, M. S., Rabibisoa, N., Randriamahatantsoa, B., Overcast, I., Raxworthy, C. J., Ruane, S., and Burbrink, F. T. (2022). Drivers of unique and asynchronous population dynamics in malagasy herpetofauna. *Journal of Biogeography*, 49(4):600-616.

883
884
Laikre, L. (2010). Genetic diversity is overlooked in international conservation policy implementation. *Conservation Genetics*, 11:349-354.

885
886
Laliberte, A. S. and Ripple, W. J. (2004). Range contractions of north american carnivores and ungulates. *BioScience*, 54(2):123-138.

887
888
889
Le Gouar, P. J., Vallet, D., David, L., Bermejo, M., Gatti, S., Levréro, F., Petit, E. J., and Ménard, N. (2009). How ebola impacts genetics of western lowland gorilla populations. *PLoS One*, 4(12):e8375.

890
891
Leblois, R., Estoup, A., and Streiff, R. (2006). Genetics of recent habitat contraction and reduction in population size: does isolation by distance matter? *Molecular Ecology*, 15(12):3601-3615.

892
893
894
Leigh, D. M., Hendry, A. P., Vázquez-Domínguez, E., and Friesen, V. L. (2019). Estimated six per cent loss of genetic variation in wild populations since the industrial revolution. *Evolutionary applications*, 12(8):1505-1512.

895
896
Leonardi, M., Barbujani, G., and Manica, A. (2021). *Genetic demography: What does it mean and how to interpret it, with a case study on the Neolithic transition*, page 91-100. Kerns Verlag.

897
LLC, W. A. (2009). Wolfram | alpha.

898
899
Lloyd, M. W., Campbell, L., and Neel, M. C. (2013). The power to detect recent fragmentation events using genetic differentiation methods. *PLoS One*, 8(5):e63981.
900
901

Loog, L. (2021). Sometimes hidden but always there: the assumptions underlying genetic inference of demographic histories. *Philosophical Transactions of the Royal Society B*, 376(1816):20190719.

902
903
904

Lourenço, A., Álvarez, D., Wang, I. J., and Velo-Antón, G. (2017). Trapped within the city: Integrating demography, time since isolation and population-specific traits to assess the genetic effects of urbanization. *Molecular Ecology*, 26(6):1498-1514.

905
906
907

Maebe, K., Meeus, I., Vray, S., Claeys, T., Dekoninck, W., Boevé, J.-L., Rasmont, P., and Smagghe, G. (2016a). A century of temporal stability of genetic diversity in wild bumblebees. *Scientific Reports*, 6(1):38289.

908
909
910

Maebe, K., Meeus, I., Vray, S., Claeys, T., Dekoninck, W., Boevé, J.-L., Rasmont, P., and Smagghe, G. (2016b). A century of temporal stability of genetic diversity in wild bumblebees. *Scientific Reports*, 6(1):38289.

911
912

Marchi, N., Schlichta, F., and Excoffier, L. (2021). Demographic inference. *Current Biology*, 31(6):R276-R279.

913
914
915

Maruyama, T. and Fuerst, P. A. (1985). Population bottlenecks and nonequilibrium models in population genetics. ii. number of alleles in a small population that was formed by a recent bottleneck. *Genetics*, 111(3):675-689.

916
917
918

Mazet, O., Rodríguez, W., Grusea, S., Boitard, S., and Chikhi, L. (2016). On the importance of being structured: instantaneous coalescence rates and human evolution—lessons for ancestral population size inference? *Heredity*, 116(4):362-371.

919
920
921

Miller, E. F., Green, R. E., Balmford, A., Maisano Delser, P., Beyer, R., Somveille, M., Leonardi, M., Amos, W., and Manica, A. (2021). Bayesian skyline plots disagree with range size changes based on species distribution models for holarctic birds. *Molecular Ecology*, 30(16):3993-4004.

922
923

Nadachowska-Brzyska, K., Konczal, M., and Babik, W. (2022). Navigating the temporal continuum of effective population size. *Methods in Ecology and Evolution*, 13(1):22-41.

924
925
926

Nadachowska-Brzyska, K., Li, C., Smeds, L., Zhang, G., and Ellegren, H. (2015). Temporal dynamics of avian populations during pleistocene revealed by whole-genome sequences. *Current Biology*, 25(10):1375-1380.

927
928
929
930

Natesh, M., Vinay, K., Ghosh, S., Jayapal, R., Mukherjee, S., Vijay, N., and Robin, V. (2020). Contrasting trends of population size change for two eurasian owlet species—athene brama and glaucidium radiatum from south asia over the late quaternary. *Frontiers in Ecology and Evolution*, 8:608339.

931
932

Nei, M., Maruyama, T., and Chakraborty, R. (1975). The bottleneck effect and genetic variability in populations. *Evolution*, pages 1-10.

933
934

Nikolic, N. and Chevalet, C. (2014). Detecting past changes of effective population size. *Evolutionary applications*, 7(6):663-681.

935
936

Notohara, M. (1990). The coalescent and the genealogical process in geographically structured population. *Journal of mathematical biology*, 29:59-75.

937
938

Nyström, V., Angerbjörn, A., and Dalén, L. (2006). Genetic consequences of a demographic bottleneck in the scandinavian arctic fox. *Oikos*, 114(1):84-94.

939
940

Ohta, T. and Kimura, M. (1973). A model of mutation appropriate to estimate the number of electrophoretically detectable alleles in a finite population. *Genetics Research*, 22(2):201-204.

941
942
943

Pacifici, M., Rondinini, C., Rhodes, J. R., Burbidge, A. A., Cristiano, A., Watson, J. E., Woinarski, J. C., and Di Marco, M. (2020). Global correlates of range contractions and expansions in terrestrial mammals. *Nature Communications*, 11(1):2840.
944

945

Pacifici, M., Santini, L., Di Marco, M., Baisero, D., Francucci, L., Marasini, G. G., Visconti, P., and Rondinini, C. (2013). Generation length for mammals. *Nature Conservation*, 5:89–94.

946

947

948

Parreira, B. R. and Chikhi, L. (2015). On some genetic consequences of social structure, mating systems, dispersal, and sampling. *Proceedings of the National Academy of Sciences*, 112(26):E3318–E3326.

949

950

Peters, R. H. and Peters, R. H. (1986). *The ecological implications of body size*, volume 2. Cambridge university press.

951

952

Pflüger, F. J., Signer, J., and Balkenhol, N. (2019). Habitat loss causes non-linear genetic erosion in specialist species. *Global Ecology and Conservation*, 17:e00507.

953

954

955

Postaire, B. D., Devloo-Delva, F., Brunnschweiler, J. M., Charvet, P., Chen, X., Cliff, G., Daly, R., Drymon, J. M., Espinoza, M., Fernando, D., et al. (2023). Global genetic diversity and historical demography of the bull shark. *Journal of Biogeography*.

956

957

958

Rasteiro, R., Bouttier, P.-A., Sousa, V. C., and Chikhi, L. (2012). Investigating sex-biased migration during the neolithic transition in europe, using an explicit spatial simulation framework. *Proceedings of the Royal Society B: Biological Sciences*, 279(1737):2409–2416.

959

960

961

Ray, N., Currat, M., Foll, M., and Excoffier, L. (2010). Splatche2: a spatially explicit simulation framework for complex demography, genetic admixture and recombination. *Bioinformatics*, 26(23):2993–2994.

962

963

964

Rogan, J. E., Parker, M. R., Hancock, Z. B., Earl, A. D., Buchholtz, E. K., Chyn, K., Martina, J., and Fitzgerald, L. A. (2023). Genetic and demographic consequences of range contraction patterns during biological annihilation. *Scientific Reports*, 13(1):1691.

965

966

Salmona, J., Heller, R., Lascoux, M., and Shafer, A. (2019). Inferring demographic history using genomic data. *Population genomics: Concepts, approaches and applications*, pages 511–537.

967

968

Salmona, J., Heller, R., Quéméré, E., and Chikhi, L. (2017). Climate change and human colonization triggered habitat loss and fragmentation in madagascar. *Molecular Ecology*, 26(19):5203–5222.

969

970

Santini, L., Di Marco, M., Visconti, P., Baisero, D., Boitani, L., and Rondinini, C. (2013). Ecological correlates of dispersal distance in terrestrial mammals. *Hystrix*, 24(2).

971

972

973

Schlötterer, C., Ritter, R., Harr, B., and Brem, G. (1998). High mutation rate of a long microsatellite allele in drosophila melanogaster provides evidence for allele-specific mutation rates. *Molecular Biology and Evolution*, 15(10):1269–1274.

974

975

Schmitt, T. (2007). Molecular biogeography of europe: Pleistocene cycles and postglacial trends. *Frontiers in zoology*, 4:1–13.

976

977

Schneider, N., Chikhi, L., Currat, M., and Radespiel, U. (2010). Signals of recent spatial expansions in the grey mouse lemur (microcebus murinus). *BMC Evolutionary Biology*, 10(1):1–17.

978

979

Schrider, D. R., Shanku, A. G., and Kern, A. D. (2016). Effects of linked selective sweeps on demographic inference and model selection. *Genetics*, 204(3):1207–1223.

980

981

Sgarlata, G. M., Maié, T., de Zoeten, T., Rasteiro, R., and Chikhi, L. (2022a). The effect of habitat loss and fragmentation on isolation-by-distance and time. *bioRxiv*, pages 2022–10.

982

983

Sgarlata, G. M., Maié, T., de Zoeten, T., Rasteiro, R., and Chikhi, L. (2022b). On the genetic consequences of habitat contraction: edge effects and habitat loss. *bioRxiv*, pages 2022–10.

984

985

Silva, M. and Downing, J. A. (1995). The allometric scaling of density and body mass: a nonlinear relationship for terrestrial mammals. *The American Naturalist*, 145(5):704–727.

986

Smith, J. M. and Slatkin, M. (1973). The stability of predator-prey systems. *Ecology*, 54(2):384–391.
987

988

Sommer, R. S. and Nadachowski, A. (2006). Glacial refugia of mammals in europe: evidence from fossil records. *Mammal Review*, 36(4):251-265.

989

990

Spence, J. P., Steinrücken, M., Terhorst, J., and Song, Y. S. (2018). Inference of population history using coalescent hmms: review and outlook. *Current opinion in genetics & development*, 53:70-76.

991

992

Srinivas, Y. and Jhala, Y. (2023). Genetic diversity, structure, and demographic histories of unique and ancient wolf lineages in india. *Conservation Genetics*, pages 1-16.

993

994

995

Storz, J. F. and Beaumont, M. A. (2002). Testing for genetic evidence of population expansion and contraction: an empirical analysis of microsatellite dna variation using a hierarchical bayesian model. *Evolution*, 56(1):154-166.

996

997

Sugg, D. W., Chesser, R. K., Stephen Dobson, F., and Hoogland, J. L. (1996). Population genetics meets behavioral ecology. *Trends in Ecology & Evolution*, 11(8):338-342.

998

999

1000

Svenning, J.-C., Eiserhardt, W. L., Normand, S., Ordonez, A., and Sandel, B. (2015). The influence of paleoclimate on present-day patterns in biodiversity and ecosystems. *Annual Review of Ecology, Evolution, and Systematics*, 46:551-572.

1001

1002

1003

Szép, E., Trubenová, B., and Csilléry, K. (2022). Using gridcoal to assess whether standard population genetic theory holds in the presence of spatio-temporal heterogeneity in population size. *Molecular Ecology Resources*, 22(8):2941-2955.

1004

1005

Tang, Q., Fung, T., Hart, D. E. T., and Rheindt, F. E. (2024). Rate and extent of genetic diversity loss under non-equilibrium scenarios of habitat loss. *Biological Conservation*, 289:110381.

1006

1007

1008

Taron, U. H., Salado, I., Escobar-Rodríguez, M., Westbury, M. V., Butschkau, S., Paijmans, J. L., VonHoldt, B. M., Hofreiter, M., and Leonard, J. A. (2021). A sliver of the past: The decimation of the genetic diversity of the mexican wolf. *Molecular ecology*, 30(23):6340-6354.

1009

1010

Team, R. C. et al. (2016). R: A language and environment for statistical computing. r foundation for statistical computing, vienna, austria. http://www.R-project.org/.

1011

1012

1013

Teixeira, H., Montade, V., Salmona, J., Metzger, J., Bremond, L., Kasper, T., Daut, G., Rouland, S., Ranarilalatiana, S., Rakotodravony, R., et al. (2021). Past environmental changes affected lemur population dynamics prior to human impact in madagascar. *Communications biology*, 4(1):1084.

1014

1015

1016

Townsend, A. K., Jones, M. L., Chen, N., Chivily, C., McAndrews, C., Clark, A. B., McGowan, K. J., and Eimes, J. (2023). Increased genetic diversity and immigration after west nile virus emergence in american crows: No evidence for a genetic bottleneck. *Molecular Ecology*.

1017

Wakeley, J. (1999). Nonequilibrium migration in human history. *Genetics*, 153(4):1863-1871.

1018

1019

Walton, W., Stone, G. N., and Lohse, K. (2021). Discordant pleistocene population size histories in a guild of hymenopteran parasitoids. *Molecular Ecology*, 30(18):4538-4550.

1020

1021

Wang, P., Hou, R., Wu, Y., Zhang, Z., Que, P., and Chen, P. (2022). Genomic status of yellow-breasted bunting following recent rapid population decline. *IsScience*, 25(7):104501.

1022

1023

1024

Welch, A. J., Wiley, A. E., James, H. F., Ostrom, P. H., Stafford Jr, T. W., and Fleischer, R. C. (2012). Ancient dna reveals genetic stability despite demographic decline: 3,000 years of population history in the endemic hawaiian petrel. *Molecular biology and evolution*, 29(12):3729-3740.

1025

1026

Wolf, C. and Ripple, W. J. (2017). Range contractions of the world's large carnivores. *Royal Society open science*, 4(7):170052.
1027 Conflict of Interest

1028 The authors declare no conflict of interest.

1029 Data and Code Availability

1030 The data that support the findings of this study are available in the supplementary material of
1031 this article. All the code used for data analysis and simulations are publicly available on a GitHub
1032 repository: https://github.com/ravi-15ms/Genetics-of-expansion-and-contraction.git

1033 Author Contribution

1034 Conceptualization: GMS, LC, RV; Methodology: LC, RV; Investigation: DSP, GMS, RV; Formal
1035 Analysis: DSP, GMS, LC, RV, RT; Contributed materials/analysis tools: DSP, GMS, LC, RV, RT,
1036 RR, TM, TP; Writing-Original Draft: DSP, GMS, LC, RT, RV; Writing-Review and Editing: DSP,
1037 GMS, LC, RV, RT, RR, TM, TP; Authors listed in alphabetical order.
<figure><img src="image_2.png" /><figcaption>Figure 1: **Simulated landscape and habitat change scenarios:** Panel A) The simulated 2D landscape consists of 18 × 18 demes with two refuge areas (SW and SE refugium) of 2 × 2 demes. The red-grey deme (top-left deme in SW refugium) corresponds to the ancestral population from which the spatial expansion starts. Panel B) Schematic representation of habitat changes over time, assuming an initial phase of expansion (green), after which habitat remains stationary (orange) before contraction (red). We considered three different scenarios characterized by the same expansion and stationary periods but different contraction durations ($t_{cont}$, y-axis). After the end of the contraction, the population was allowed to survive in the refuge areas. Individuals were sampled only from the refuge areas (see text for details). Altogether, the population starts with K individuals in the SW refugium, grows to a maximum possible of 324 demes × K, and decreases to two refugium populations of 4K each.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2: Dynamics of genetic diversity in the Southwest (SW) refugium population:**
Temporal evolution of $H_e$ estimated for 40 individuals sampled in the SW refugium population at discrete times. Two generation lengths are considered: $GL = 0.5$ years (Panel A) and $GL = 5$ years (Panel B). In both panels, each point represents the average over 15 simulations, while the error bars represent the standard deviations. The discrete time points are connected by dashed lines and their symbols and colors represent the duration of the contraction phase for three scenarios, corresponding to $t_{cont} = 340$ years (red circles), $t_{cont} = 4,080$ years (blue triangles) and $t_{cont} = 10,880$ years (green squares) respectively. The dashed horizontal line (- -) represents the value of $H_e$ at mutation-migration-drift equilibrium when only the SW refugium is occupied ($H_e^{refuge}$), while the dotted line (· · ·) represents the equilibrium value when the entire landscape is occupied ($H_e^{landscape}$). The four vertical lines, from left to right, represent the following events: i) start of the contraction phase (black, time = 2,860 years for all scenarios) and end of the contraction phase for (ii) $t_{cont} = 340$ years (red, time = 3,200 years) iii) $t_{cont} = 4,080$ years (blue, time = 6,940 years) and iv) $t_{cont} = 10,880$ years (green, time = 13,740 years). The other parameter values for both plots are ($K, m$)=(50, 0.2). Panel B shows that when contraction is slow (here 10,880 years), $H_e$ can keep increasing for more than 7,000 years after contraction started and keep high and stable levels of genetic diversity that are much higher than $H_e^{refuge}$ for more than 10,000 years.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: *H*<sub>e</sub> in the SW refugium population at the end of contraction: Two generation lengths are considered: *GL* = 0.5 years (Panel A) and *GL* = 5 years (Panel B). In both panels, colors represent the duration of the contraction phase for three scenarios, corresponding to *t*<sub>cont</sub> = 340 years (red), *t*<sub>cont</sub> = 4,080 years (blue), and *t*<sub>cont</sub> = 10,880 years (green) respectively. The height of the histogram bars represents the means and the error bars represent the standard deviations measured across 15 simulations. The other parameter values for both plots are (*K*,*m*)=(50,0.2). Panel A suggests that quick contraction "maintains" more genetic diversity than slow contraction. Panel B suggests the opposite with a significant difference even when populations are sampled 310 years after the end of contraction.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4: Species-specific traits and dynamics of H<sub>e</sub> in SW refugium:** Panel A: Generation length as a function of body mass across different orders of mammalian species (data from (Pacifici et al., 2013)). Panels B-C: Temporal dynamics of H<sub>e</sub> in the SW refugium population for 40 individuals sampled at discrete times. Panel B represents a species with (GL = 5 years, K = 150 and m = 0.059, orange square in panel A) whereas Panel C represents a species with (GL = 10 years, K = 1,000 and m = 0.003, violet triangle in panel A). Colors, symbols, and lines in panels B and C encode the same information as in Figure 2.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Genetic diversity in the panmictic model and spatial simulation validation:

Panel A: Schematic representation of population size changes in our panmictic single population model. Panel B: Time evolution of $H_e$ in a panmictic population under size changes represented in panel A. Three distinct temporal dynamics of $H_e$ are observed by changing the initial ($N_0$) and expanded population ($N_{exp}$) sizes, keeping other parameters constant. The $N_0$ represents the refugium population size and $N_{exp}$ represents the population size when the whole habitat is colonized. Specifically, we set ($N_0, N_{exp}$) = (50, 250) (Panel B1), ($N_0, N_{exp}$) = (100, 10000) (Panel B2) and ($N_0, N_{exp}$) = (4000, 12000) (Panel B3), three scenarios corresponding to a five, hundred and three-fold increase in population size, respectively. In the three subpanels (B1, B2, and B3), the two vertical lines, from left to right, represent the following events, respectively: i) the start of the contraction phase (black, time = 1,000 years) and ii) the end of contraction (red, time = 1,800 years). Panel C: $H_e$ estimated at the end of the contraction phase, as a function of the duration of the contraction (x-axis, $t_{cont}$) and stationary (y-axis, $t_{stat}$) phases. The two horizontal dashed lines represent two transects for two $t_{stat}$ values, $t_{stat}$ = 2,500 generations (black) and $t_{stat}$ = 250 generations (grey). Panel D: $H_e$ estimated at the end of contraction as a function of contraction durations, for two stationary times, $t_{stat}$ = 2,500 generations (black squares) and $t_{stat}$ = 250 generations (grey circles), represented by the dashed horizontal lines in Panel C. Panels C and D were obtained for ($N_0, N_{exp}$) = (100, 10000), i.e. a scenario with two orders of magnitude change in population size. Panel E: Spatial simulation analysis: $H_e$ estimated at the end of contraction for ($t_{stat}$ = 2,500 years) as a function of contraction duration. The lines and symbols correspond to two generation lengths, $GL$ = 0.5 years (black, 5,000 generations in stationary phase) and $GL$ = 5 years (grey, 500 generations in stationary phase). Each point represents the mean and the error bars represent the standard deviations measured across 15 simulations. The other parameters for the spatial simulations are set to ($K, m$)=(100, 0.2).</figcaption></figure>
<figure><img src="image_7.png" /></figure>

**Figure 6: Summary of main results:** This figure summarizes the main findings of our study, illustrating the trajectory of genetic diversity during habitat contraction in a scenario (left panel) and comparing diversity levels at the end of the contraction period across scenarios with varying contraction durations (right panel). The top section of the left panel shows a schematic representation of the demographic scenarios simulated (see Figure 1 for further details). Three possible trajectories of genetic diversity are shown, representing the diversity within the refugium after contraction begins. These trajectories depend on the genetic diversity level reached before the contraction starts, in relation to two equilibrium diversity levels: diversity expected when the entire landscape is inhabited ($H_e^{landscape}$) and diversity expected when only the refugium is occupied ($H_e^{refuge}$). The right panel highlights the effect of contraction duration (or speed) on genetic diversity at the end of the contraction phase across scenarios. In some scenarios, fast contractions retain more genetic diversity than slow contractions, while in others, the opposite effect is observed. Between the left and right panels, we illustrate the individual effects of species-specific traits (population density, dispersal, and generation length) on the genetic diversity patterns shown in the left and right panels. For instance, everything else being equal, species with long generation lengths will tend to increase their diversity for long periods after the habitat contraction has started (left panel, bottom curve) as they require more time to accumulate genetic diversity and reach the equilibrium within the refugium. In contrast, species with low population densities reach mutation-migration-drift equilibrium quicker than species with high densities. Note that this representation does not illustrate the correlations among traits but rather how each variable influences the genetic diversity patterns, with other parameters being constant. We stress that these traits are often correlated (see Materials and Methods), making predictions difficult in real species. Overall, these findings highlight some general trends while also emphasizing the complex interplay between species-specific traits and habitat dynamics in shaping genetic diversity through time.
# Supplementary Materials

## Species-specific traits shape genetic diversity during expansion-contraction and bias demographic history reconstruction

Ravi Vishwakarma, Gabriele Maria Sgarlata, David Soriano-Paños, Rita Rasteiro, Tiago Maié, Tiago Paixão, Rémi Tournebize, Lounès Chikhi

# Contents

<table><tr><td>S1</td><td>Comparison of simulated habitat changes with Arenas et al. (2012)</td><td>2</td></tr><tr><td>S2</td><td>Patterns of F<sub>st</sub> within SW and between SW and SE refugium</td><td>2</td></tr><tr><td>S3</td><td>Expression for heterozygosity in the panmictic population under size changes</td><td>5</td></tr><tr><td>S4</td><td>Tables</td><td>7</td></tr><tr><td>S5</td><td>Figures</td><td>10</td></tr><tr><td>S6</td><td>Demographic inference analysis using VarEff</td><td>15</td></tr></table>
## S1 Comparison of simulated habitat changes with Arenas et al. (2012)

The spatial simulated scenarios carried out in the present study are inspired by and adapted from the simulations carried out by Arenas et al. (2012). They simulated scenarios of expansion and contraction using a spatial simulation tool (SPLATCHE Ray et al. (2010)) that inspired the development of SINS (Rasteiro et al., 2012). These similar frameworks allows for meaningful validation and comparisons between their study and ours. Our scenarios are similar to Arenas et al. (2012) when we set the generation length as 0.5 years and when carrying capacities are 50, 100, or 200, but can differ for other generation lengths and carrying capacities. Using a generation length of 0.5 years in our simulations, the time of expansion and stationary phase were 720 and 5,000 generations respectively, whereas in Arenas et al. (2012), these were 550 and 5,000 generations. The discrepancy stems from the fact that Arenas et al. (2012) simulated 51 demes per size of their landscape, while we simulated 18 for computational optimization. However, we kept the speed ratios for expansion and contraction roughly the same. For example, in Arenas et al. (2012), the ratio of slowest to fastest contraction was 50, whereas in our case it was 32. We found that, despite using a forward simulation approach as opposed to a coalescent approach in SPLATCHE and a smaller world, the simulation framework we used was capable of reproducing the results of Arenas et al. (2012) when replicating their sampling design (i.e. when sampling after the end each of contraction phase) (Figure 3A).

## S2 Patterns of $F_{st}$ within SW and between SW and SE refugium

To calculate $F_{st}$ within the SW refugium, we employed the method outlined in the Materials and Methods section, which included sampling 10 individuals from each of the four demes. Similarly, to assess $F_{st}$ between the SW and SE refugia, we sampled 20 individuals from each of the SW and SE refugia. Note that Data points do not start at time = 0 because the SE refugium is empty at that time. Figure S1 represents the evolution of $F_{st}$ between SW and SE refugium for two species that differ in GL but have the same carrying capacity ($K = 50$) and migration rate ($m = 0.2$). Figure S1 A refers to a species with $GL = 0.5$ years, and Figure S1 B refers to another species with $GL = 5$ years. For both species, $F_{st}$ between refuge areas are similar at the start of the contraction phase, $t = 2,860$ years. After range contraction starts, $F_{st}$ increases positively as a function of contraction speed. This is because the fastest the contraction, the more time the SW and SE refugia will remain isolated, thus leading to greater genetic differentiation between the two refugia. We also observe that $F_{st}$ increases steeply after the end of each contraction phase, since the SW and SE refugia get completely disconnected at the last step of the contraction process.

When comparing $F_{st}$ between refuge areas at the end of contraction (Figure S2), following the approach in Arenas et al. (2012), we observe that for both generation lengths, $F_{st}$ at the end of slow contraction was higher compared to fast contraction speed. For $F_{st}$ within SW refugium at the end of contraction, the values did not differ significantly across different generations lengths (Table S3). These results, regarding the effect of contraction speed on $F_{st}$ within and between refugia, align with the findings reported by Arenas et al. (2012).
<figure><img src="image_8.png" /><figcaption>Figure S1: *F*<sub>st</sub> dynamics between SW and SE refugium: Temporal evolution of F<sub>st</sub> between SW and SE refugium at discrete times. Two generation lengths are considered: GL = 0.5 years (Panel A) and GL = 5 years (Panel B). In both panels, each discrete point represents the mean and the error bars represent the standard deviations measured across 15 simulations. The discrete time points are connected by dashed lines and their symbols and colors represent the duration of the contraction phase for three scenarios, corresponding to t<sub>cont</sub> = 340 years (red circles), t<sub>cont</sub> = 4080 years (blue triangles) and t<sub>cont</sub> = 10880 years (green squares) respectively. The four vertical lines, from left to right, represents the following events respectively: i) start of the contraction phase (black, time = 2860 years) and end of the contraction phases for (ii) t<sub>cont</sub> = 340 years (red, time = 3200 years) iii) t<sub>cont</sub> = 4080 years (blue, time = 6940 years) and iv) t<sub>cont</sub> = 10880 years (green, time = 13740 years). The other parameter values for both plots are (K,m)=(50,0.2).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure S2: *F*<sub>st</sub> between SW and SE refugium population after the end of the contraction period: *F*<sub>st</sub> between SW and SE refugium for different sampling times after the end of contraction. Two generation lengths are considered: *GL* = 0.5 years (Panel A) and *GL* = 5 years (Panel B). In both panel, colors represent the duration of the contraction phase for three scenarios, corresponding to *t*<sub>cont</sub> = 340 years (red), *t*<sub>cont</sub> = 4080 years (blue) and *t*<sub>cont</sub> = 10880 years (green) respectively. The bars represent the mean and the error bars represent the standard deviations measured across 15 simulations. The other parameter values for both plots are (K,m)=(50,0.2).</figcaption></figure>
## S3 Expression for heterozygosity in the panmictic population under size changes

The time evolution of the population size in our panmictic model is given by:

$$N(t) = \begin{cases} N_i & \text{if } t = 0 \\ N_{exp} & \text{if } 0 < t \le t_{stat} \\ N_{exp} - \left( \frac{N_{exp}-N_i}{t_{cont}} \right) (t-t_{stat}) & \text{if } t_{stat} < t \le t_{stat} + t_{cont} \\ N_i & \text{if } t > t_{stat} + t_{cont} \end{cases}$$

where $N_i$ is the initial population size, $N_{exp}$ is the final population size, $t_{stat}$ is the duration of stationary time and $t_{cont}$ is the duration in contraction time. Times are in generations.

Recalling Equation 8 in the Materials and Methods section, the differential equation for the homozygosity $F$ is given by:

$$\frac{dF}{dt} = (1-u)^2 F \left[ 1 - \frac{1}{2N} - \frac{1}{(1-u)^2} \right] + \frac{(1-u)^2}{2N},$$

We start by assuming that a homozygous population ($F(t=0)=1$) of $N_i$ undergoes an instantaneous expansion to reach $N_{exp}$ individuals. Afterwards, the population size remains constant to $N_{exp}$ individuals during a stationary period of $t_{stat}$ generations. To find the evolution of the homozygosity in the stationary phase, hereinafter referred to as $F_{stat}$, we solve the following differential equation:

$$\frac{dF_{stat}}{dt} = (1-u)^2 F_{stat} \left[ 1 - \frac{1}{2N_{exp}} - \frac{1}{(1-u)^2} \right] + \frac{(1-u)^2}{2N_{exp}},$$

integrating using initial condition $F(t=0)=1$, gives

$$F_{stat}(t) = F_{N_{exp}} + (1 - F_{N_{exp}}) e^{\left[ (1-u)^2 \left( 1 - \frac{1}{2N_{exp}} - \frac{1}{(1-u)^2} \right) \right] t}.$$

After the stationary phase, the population starts contracting for a duration of $t_{cont}$. We modeled the population size contraction as a linear reduction in population size from $N_{exp}$ to $N_i$ individuals. Mathematically,

$$N(t) = N_{exp} - \frac{(N_{exp} - N_i)}{t_{cont}} (t - t_{stat})$$

Substituting $N(t)$ in equation 8, we get the differential equation governing the time evolution of the heterozygosity in the contraction phase, $F_{cont}$:

$$\frac{dF_{cont}}{dt} = \left\{ (1-u)^2 F_{cont} \left[ 1 - \frac{1}{2 \left( N_{exp} - \frac{(N_{exp}-N_i)}{t_{cont}} (t-t_{stat}) \right)} - \frac{1}{(1-u)^2} \right] + \frac{(1-u)^2}{2 \left( N_{exp} - \frac{(N_{exp}-N_i)}{t_{cont}} (t-t_{stat}) \right)} \right\},$$

To ease the notation, we define $\Delta t_{cont} = t - t_{stat}$, $a = (1-u)^2$, $b = 2N_{exp}$, $c = 2 \left( \frac{N_{exp}-N_i}{t_{cont}} \right)$.

The former equation transforms to:
$$ \frac{dF_{cont}}{d\Delta t_{cont}} = aF_{cont}\left[1 - \frac{1}{b - c\Delta t_{cont}} - \frac{1}{a}\right] + \frac{a}{b - c\Delta t_{cont}}, $$

We utilize WolframAlpha LLC (2009) to integrate the above equation, setting initial condition as $k_1 = F_{stat}(t = t_{stat})$. The specific command we used to solve this equation is:

$$ \text{solve} \quad y' = y(t) a \left(1 - \frac{1}{b - ct} - \frac{1}{a}\right) + \frac{a}{b - ct}, y(0) = k_1, $$

Following the LLC (2009) integration, the homozygosity $F_{cont}(\Delta t_{cont})$ during contraction phase takes the form:

$$ F_{cont}(\Delta t_{cont}) = F_{stat}(t_{stat})e^{(a-1)\Delta t_{cont}} \left(1 - \frac{c}{b}\Delta t_{cont}\right)^{a/c} + e^{(a-1)(c\Delta t_{cont} - b)/c} \left(\frac{a}{c}\right) \left(\frac{(a-1)(c\Delta t_{cont} - b)}{c}\right)^{a/c} \left\{ \Gamma\left(\frac{-a}{c}, \frac{(a-1)(c\Delta t_{cont} - b)}{c}\right) - \Gamma\left(\frac{-a}{c}, \frac{(1-a)b}{c}\right) \right\}, $$

After the end of contraction, the population size remains constant ($N_i$ individuals). Therefore, the homozygosity after the end of contraction, denoted by $F_{cont-refuge}$, is governed by the following equation:

$$ \frac{dF_{cont-refuge}}{d\Delta t_{cont-refuge}} = (1-u)^2 F_{cont-refuge} \left[1 - \frac{1}{2N_i} - \frac{1}{(1-u)^2}\right] + \frac{(1-u)^2}{2N_i}, $$

where $\Delta t_{cont-refuge} = t - t_{stat} - t_{cont}$. This equation is similar to the one during the stationary phase but has the initial condition $F_{cont-refuge}(0) = F_{cont}(t_{cont})$. The expression of homozygosity after the end of contraction thus reads:

$$ F_{cont-refuge}(\Delta t_{cont-refuge}) = F_{N_i} + (F_{cont}(t_{cont}) - F_{N_i}) e^{\left[(1-u)^2 \left(1 - \frac{1}{2N_{exp}} - \frac{1}{(1-u)^2}\right)\right]} \Delta t_{cont-refuge}. $$
## S4 Tables

Abbreviations:

* GL: generation length (years)

* K: carrying capacity of a deme

* σ: dispersal distance (km/generation)

* m: migration rate (per generation)

* MNA: mean number of alleles within SW refugium

* $H_e$: expected heterozygosity within SW refugium

* $F_{st}$: genetic differentiation

<table><caption>Table S1: **Parameters values:** density and σ deduced using known allometric relationship (cf. Materials and Methods). $H_e$ at time = 2,860 years, is the amount of genetic diversity accumulated when contraction starts, values (mean ± standard deviation) estimated over 15 independent simulations.</caption><thead><tr><th>GL (years)</th><th>K (ind/deme)</th><th>density (ind/km²)</th><th>body mass (kg)</th><th>σ (km/gen)</th><th>m (in stepping stone model)</th><th>$H_e$ at time = 2,860 years</th></tr></thead><tbody><tr><td>2</td><td>650</td><td>0.065</td><td>10.81</td><td>8.15</td><td>0.007</td><td>0.712 ± 0.008</td></tr><tr><td>2</td><td>1000</td><td>0.1</td><td>6.09</td><td>5.91</td><td>0.003</td><td>0.724 ± 0.008</td></tr><tr><td>5</td><td>150</td><td>0.015</td><td>76.34</td><td>24.30</td><td>0.059</td><td>0.530 ± 0.013</td></tr><tr><td>5</td><td>250</td><td>0.025</td><td>38.65</td><td>16.61</td><td>0.028</td><td>0.534 ± 0.018</td></tr><tr><td>5</td><td>400</td><td>0.04</td><td>20.65</td><td>11.70</td><td>0.014</td><td>0.541 ± 0.024</td></tr><tr><td>5</td><td>1000</td><td>0.1</td><td>6.09</td><td>5.91</td><td>0.003</td><td>0.567 ± 0.001</td></tr><tr><td>10</td><td>150</td><td>0.015</td><td>76.34</td><td>24.30</td><td>0.059</td><td>0.374 ± 0.018</td></tr><tr><td>10</td><td>250</td><td>0.025</td><td>38.65</td><td>16.61</td><td>0.028</td><td>0.375 ± 0.016</td></tr><tr><td>10</td><td>1000</td><td>0.1</td><td>6.09</td><td>5.91</td><td>0.003</td><td>0.379 ± 0.012</td></tr></tbody></table>

<table><caption>Table S2: **Mutation-Drift-Migration equilibrium:** $H_e^{refuge}$ and $H_e^{landscape}$, estimated over 15 independent simulations. Equilibrium was estimated by running simulations for 99,000 generations, sampling every 1,000 generation, and calculating the mean ± standard deviation based on the last 10 time points.</caption><thead><tr><th>K (ind/deme)</th><th>m (in stepping stone model)</th><th>$H_e^{refuge}$</th><th>$H_e^{landscape}$</th></tr></thead><tbody><tr><td>50</td><td>0.1</td><td>0.375 ± 0.005</td><td>0.716 ± 0.004</td></tr><tr><td>50</td><td>0.2</td><td>0.377 ± 0.004</td><td>0.784 ± 0.002</td></tr><tr><td>100</td><td>0.1</td><td>0.510 ± 0.003</td><td>0.829 ± 0.003</td></tr><tr><td>100</td><td>0.2</td><td>0.510 ± 0.006</td><td>0.869 ± 0.001</td></tr><tr><td>150</td><td>0.059</td><td>0.589 ± 0.003</td><td>0.841 ± 0.003</td></tr><tr><td>250</td><td>0.028</td><td>0.672 ± 0.003</td><td>0.853 ± 0.001</td></tr><tr><td>400</td><td>0.014</td><td>0.738 ± 0.002</td><td>0.867 ± 0.002</td></tr><tr><td>650</td><td>0.007</td><td>0.788 ± 0.002</td><td>0.877 ± 0.001</td></tr><tr><td>1000</td><td>0.003</td><td>0.836 ± 0.001</td><td>0.892 ± 0.001</td></tr></tbody></table>
<table><caption>Table S3: Molecular diversity measures at the end of contraction for Southwest (SW) *refugium*. Values are given as mean ± standard deviation, estimated over 15 independent simulations.</caption><thead><tr><th>GL (years)</th><th>m (/gen)</th><th>K</th><th>contraction duration (years)</th><th>time post contraction (years)</th><th>sampling time (years)</th><th>MNA</th><th>H<sub>e</sub></th><th>F<sub>st</sub> (within SW)</th></tr></thead><tbody><tr><td rowspan="10">0.5</td><td rowspan="10">0.1</td><td rowspan="10">100</td><td rowspan="3">340</td><td>10</td><td>3210</td><td>7.08 ± 0.06</td><td>0.727 ± 0.005</td><td>0.047 ± 0.002</td></tr><tr><td>60</td><td>3260</td><td>5.96 ± 0.06</td><td>0.674 ± 0.005</td><td>0.044 ± 0.002</td></tr><tr><td>160</td><td>3360</td><td>4.96 ± 0.06</td><td>0.601 ± 0.006</td><td>0.046 ± 0.002</td></tr><tr><td rowspan="3">4080</td><td>10</td><td>6950</td><td>6.29 ± 0.07</td><td>0.676 ± 0.005</td><td>0.045 ± 0.002</td></tr><tr><td>60</td><td>7000</td><td>5.53 ± 0.05</td><td>0.639 ± 0.006</td><td>0.042 ± 0.002</td></tr><tr><td>160</td><td>7100</td><td>4.81 ± 0.03</td><td>0.580 ± 0.006</td><td>0.041 ± 0.002</td></tr><tr><td rowspan="3">10880</td><td>10</td><td>13750</td><td>5.94 ± 0.07</td><td>0.663 ± 0.004</td><td>0.046 ± 0.002</td></tr><tr><td>60</td><td>13800</td><td>5.25 ± 0.05</td><td>0.624 ± 0.004</td><td>0.045 ± 0.002</td></tr><tr><td>160</td><td>13900</td><td>4.67 ± 0.04</td><td>0.581 ± 0.004</td><td>0.047 ± 0.002</td></tr><tr><td rowspan="10">5</td><td rowspan="10">0.1</td><td rowspan="10">100</td><td rowspan="3">340</td><td>10</td><td>3210</td><td>4.06 ± 0.04</td><td>0.521 ± 0.006</td><td>0.043 ± 0.002</td></tr><tr><td>60</td><td>3260</td><td>3.92 ± 0.03</td><td>0.519 ± 0.006</td><td>0.043 ± 0.002</td></tr><tr><td>160</td><td>3360</td><td>3.86 ± 0.04</td><td>0.513 ± 0.007</td><td>0.048 ± 0.002</td></tr><tr><td rowspan="3">4080</td><td>10</td><td>6950</td><td>4.89 ± 0.04</td><td>0.636 ± 0.003</td><td>0.041 ± 0.002</td></tr><tr><td>60</td><td>7000</td><td>4.80 ± 0.04</td><td>0.629 ± 0.004</td><td>0.045 ± 0.002</td></tr><tr><td>160</td><td>7100</td><td>4.64 ± 0.04</td><td>0.619 ± 0.004</td><td>0.046 ± 0.002</td></tr><tr><td rowspan="3">10880</td><td>10</td><td>13750</td><td>5.64 ± 0.07</td><td>0.677 ± 0.004</td><td>0.042 ± 0.002</td></tr><tr><td>60</td><td>13800</td><td>5.49 ± 0.07</td><td>0.671 ± 0.004</td><td>0.045 ± 0.002</td></tr><tr><td>160</td><td>13900</td><td>5.39 ± 0.06</td><td>0.661 ± 0.003</td><td>0.049 ± 0.002</td></tr><tr><td rowspan="10">0.5</td><td rowspan="10">0.2</td><td rowspan="10">100</td><td rowspan="3">340</td><td>10</td><td>3210</td><td>7.84 ± 0.07</td><td>0.776 ± 0.003</td><td>0.020 ± 0.001</td></tr><tr><td>60</td><td>3260</td><td>6.34 ± 0.07</td><td>0.707 ± 0.004</td><td>0.023 ± 0.001</td></tr><tr><td>160</td><td>3360</td><td>5.15 ± 0.05</td><td>0.617 ± 0.007</td><td>0.026 ± 0.002</td></tr><tr><td rowspan="3">4080</td><td>10</td><td>6950</td><td>6.86 ± 0.06</td><td>0.721 ± 0.003</td><td>0.023 ± 0.001</td></tr><tr><td>60</td><td>7000</td><td>5.80 ± 0.05</td><td>0.667 ± 0.005</td><td>0.022 ± 0.001</td></tr><tr><td>160</td><td>7100</td><td>4.89 ± 0.06</td><td>0.600 ± 0.005</td><td>0.020 ± 0.002</td></tr><tr><td rowspan="3">10880</td><td>10</td><td>13750</td><td>6.45 ± 0.05</td><td>0.700 ± 0.005</td><td>0.024 ± 0.002</td></tr><tr><td>60</td><td>13800</td><td>5.47 ± 0.06</td><td>0.651 ± 0.005</td><td>0.025 ± 0.001</td></tr><tr><td>160</td><td>13900</td><td>4.75 ± 0.04</td><td>0.587 ± 0.005</td><td>0.020 ± 0.001</td></tr><tr><td rowspan="10">5</td><td rowspan="10">0.2</td><td rowspan="10">100</td><td rowspan="3">340</td><td>10</td><td>3210</td><td>4.21 ± 0.03</td><td>0.549 ± 0.004</td><td>0.021 ± 0.002</td></tr><tr><td>60</td><td>3260</td><td>4.12 ± 0.05</td><td>0.544 ± 0.004</td><td>0.023 ± 0.002</td></tr><tr><td>160</td><td>3360</td><td>3.98 ± 0.03</td><td>0.542 ± 0.004</td><td>0.023 ± 0.002</td></tr><tr><td rowspan="3">4080</td><td>10</td><td>6950</td><td>5.19 ± 0.02</td><td>0.664 ± 0.003</td><td>0.020 ± 0.001</td></tr><tr><td>60</td><td>7000</td><td>5.05 ± 0.02</td><td>0.659 ± 0.003</td><td>0.022 ± 0.001</td></tr><tr><td>160</td><td>7100</td><td>4.87 ± 0.03</td><td>0.646 ± 0.004</td><td>0.024 ± 0.001</td></tr><tr><td rowspan="3">10880</td><td>10</td><td>13750</td><td>6.03 ± 0.05</td><td>0.711 ± 0.003</td><td>0.021 ± 0.001</td></tr><tr><td>60</td><td>13800</td><td>5.94 ± 0.05</td><td>0.702 ± 0.003</td><td>0.023 ± 0.002</td></tr><tr><td>160</td><td>13900</td><td>5.64 ± 0.06</td><td>0.690 ± 0.003</td><td>0.025 ± 0.001</td></tr></tbody></table>
<table><caption>Table S4: Molecular diversity measures at the end of contraction when sampling individuals from both Southwest (SW) and Southeast (SE) refugia. Values are given as mean ± standard deviation, estimated over 15 independent simulations.</caption><thead><tr><th>GL (years)</th><th>m (/gen)</th><th>K</th><th>contraction duration (years)</th><th>time post contraction (years)</th><th>sampling time (years)</th><th>MNA</th><th>H<sub>e</sub></th><th>F<sub>st</sub> (between SW and SE)</th></tr></thead><tbody><tr><td rowspan="10">0.5</td><td rowspan="10">0.1</td><td rowspan="10">100</td><td rowspan="3">340</td><td>10</td><td>3210</td><td>8.92±0.07</td><td>0.760±0.002</td><td>0.142±0.003</td></tr><tr><td>60</td><td>3260</td><td>8.02±0.06</td><td>0.697±0.002</td><td>0.216±0.004</td></tr><tr><td>160</td><td>3360</td><td>7.03±0.06</td><td>0.624±0.005</td><td>0.299±0.005</td></tr><tr><td rowspan="3">4080</td><td>10</td><td>6950</td><td>8.76±0.08</td><td>0.701±0.003</td><td>0.227±0.005</td></tr><tr><td>60</td><td>7000</td><td>8.03±0.05</td><td>0.667±0.003</td><td>0.267±0.004</td></tr><tr><td>160</td><td>7100</td><td>7.24±0.05</td><td>0.603±0.004</td><td>0.339±0.005</td></tr><tr><td rowspan="3">10880</td><td>10</td><td>13750</td><td>8.46±0.08</td><td>0.689±0.003</td><td>0.239±0.004</td></tr><tr><td>60</td><td>13800</td><td>7.77±0.07</td><td>0.649±0.003</td><td>0.283±0.005</td></tr><tr><td>160</td><td>13900</td><td>7.11±0.06</td><td>0.601±0.003</td><td>0.339±0.004</td></tr><tr><td rowspan="10">5</td><td rowspan="10">0.1</td><td rowspan="10">100</td><td rowspan="3">340</td><td>10</td><td>3210</td><td>4.38±0.03</td><td>0.530±0.005</td><td>0.105±0.003</td></tr><tr><td>60</td><td>3260</td><td>4.35±0.03</td><td>0.530±0.004</td><td>0.114±0.003</td></tr><tr><td>160</td><td>3360</td><td>4.28±0.03</td><td>0.528±0.004</td><td>0.124±0.004</td></tr><tr><td rowspan="3">4080</td><td>10</td><td>6950</td><td>5.53±0.02</td><td>0.658±0.003</td><td>0.121±0.004</td></tr><tr><td>60</td><td>7000</td><td>5.49±0.03</td><td>0.654±0.004</td><td>0.124±0.004</td></tr><tr><td>160</td><td>7100</td><td>5.41±0.03</td><td>0.641±0.003</td><td>0.143±0.006</td></tr><tr><td rowspan="3">10880</td><td>10</td><td>13750</td><td>6.66±0.06</td><td>0.697±0.003</td><td>0.154±0.004</td></tr><tr><td>60</td><td>13800</td><td>6.57±0.05</td><td>0.691±0.003</td><td>0.165±0.004</td></tr><tr><td>160</td><td>13900</td><td>6.45±0.06</td><td>0.679±0.002</td><td>0.178±0.003</td></tr><tr><td rowspan="10">0.5</td><td rowspan="10">0.2</td><td rowspan="10">100</td><td rowspan="3">340</td><td>10</td><td>3210</td><td>9.30±0.06</td><td>0.787±0.002</td><td>0.106±0.002</td></tr><tr><td>60</td><td>3260</td><td>8.13±0.06</td><td>0.717±0.002</td><td>0.184±0.004</td></tr><tr><td>160</td><td>3360</td><td>7.14±0.06</td><td>0.630±0.003</td><td>0.289±0.005</td></tr><tr><td rowspan="3">4080</td><td>10</td><td>6950</td><td>9.14±0.05</td><td>0.734±0.001</td><td>0.178±0.003</td></tr><tr><td>60</td><td>7000</td><td>8.19±0.09</td><td>0.680±0.004</td><td>0.241±0.004</td></tr><tr><td>160</td><td>7100</td><td>7.23±0.07</td><td>0.611±0.005</td><td>0.321±0.004</td></tr><tr><td rowspan="3">10880</td><td>10</td><td>13750</td><td>8.52±0.05</td><td>0.711±0.003</td><td>0.194±0.005</td></tr><tr><td>60</td><td>13800</td><td>7.80±0.07</td><td>0.665±0.004</td><td>0.241±0.004</td></tr><tr><td>160</td><td>13900</td><td>6.97±0.06</td><td>0.601±0.003</td><td>0.322±0.004</td></tr><tr><td rowspan="10">5</td><td rowspan="10">0.1</td><td rowspan="10">100</td><td rowspan="3">340</td><td>10</td><td>3210</td><td>4.56±0.03</td><td>0.559±0.003</td><td>0.070±0.003</td></tr><tr><td>60</td><td>3260</td><td>4.47±0.03</td><td>0.559±0.004</td><td>0.080±0.003</td></tr><tr><td>160</td><td>3360</td><td>4.45±0.03</td><td>0.552±0.003</td><td>0.095±0.003</td></tr><tr><td rowspan="3">4080</td><td>10</td><td>6950</td><td>5.71±0.03</td><td>0.676±0.002</td><td>0.090±0.003</td></tr><tr><td>60</td><td>7000</td><td>5.60±0.02</td><td>0.667±0.002</td><td>0.102±0.003</td></tr><tr><td>160</td><td>7100</td><td>5.51±0.04</td><td>0.658±0.003</td><td>0.119±0.003</td></tr><tr><td rowspan="3">10880</td><td>10</td><td>13750</td><td>6.92±0.05</td><td>0.722±0.003</td><td>0.119±0.004</td></tr><tr><td>60</td><td>13800</td><td>6.81±0.04</td><td>0.716±0.002</td><td>0.129±0.003</td></tr><tr><td>160</td><td>13900</td><td>6.67±0.05</td><td>0.702±0.002</td><td>0.147±0.003</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Figure S3: Dynamics of the census population size ($N_c$) in the whole landscape: Two generation lengths are considered: $GL = 0.5$ years (Panel A) and $GL = 5$ years (Panel B). In both panels, each discrete point represents the mean and the error bars represent the standard deviations measured across 15 simulations. The discrete time points are connected by dashed lines and their symbols and colors represent the duration of the contraction phase for three scenarios, corresponding to $t_{cont} = 340$ years (red circles), $t_{cont} = 4,080$ years (blue triangles) and $t_{cont} = 10,880$ years (green squares) respectively. The horizontal black line represents the maximum census population size for the carrying capacity in the landscape (including the two refugia). The four vertical lines, from left to right, represent the following events respectively: i) start of the contraction phase (black, time = 2,860 years) and end of the contraction phases for (ii) $t_{cont} = 340$ years (red, time = 3,200 years) iii) $t_{cont} = 4,080$ years (blue, time = 6,940 years) and iv) $t_{cont} = 10,880$ years (green, time = 13,740 years). The other parameter values for both plots are ($K,m$)=(50,0.2). These figures show that $N_c$ follows habitat change with a slight delay during the expansion phase, with a shorter delay when the generation time is shorter. For instance, when the expansion phase is finished (end of the green area) $N_c$ is close to its maximum for $GL = 0.5$ years (Panel A) whereas it only reaches about a third of its maximum size for $GL = 5$ years (Panel B).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure S4: Dynamics of the mean number of alleles (MNA) in the SW refugium:** Two generation lengths are considered: GL = 0.5 years (Panel A) and GL = 5 years (Panel B). In both panels, each discrete point represents the mean and the error bars represent the standard deviations measured across 15 simulations. The discrete time points are connected by dashed lines and their symbols and colors represent the duration of the contraction phase for three scenarios, corresponding to t<sub>cont</sub> = 340 years (red circles), t<sub>cont</sub> = 4,080 years (blue triangles) and t<sub>cont</sub> = 10,880 years (green squares) respectively. The dotted (dashed) horizontal line represents the value of MNA at mutation-migration-drift equilibrium when the refugium (entire landscape) is occupied *MNA*<sup>refuge</sup> (*MNA*<sup>landscape</sup>). Both equilibria are computed by keeping stable conditions across 99,000 generations. The four vertical lines, from left to right, represent the following events respectively: i) start of the contraction phase (black, time = 2,860 years) and end of the contraction phases for (ii) t<sub>cont</sub> = 340 years (red, time = 3,200 years) iii) t<sub>cont</sub> = 4,080 years (blue, time = 6,940 years) and iv) t<sub>cont</sub> = 10,880 years (green, time = 13,740 years). The other parameter values for both plots are (K,m)=(50,0.2). The curves show trajectories that are very similar to those of *H*<sub>e</sub>. In addition they resemble the change in *N*<sub>c</sub> shown in figure S3 but are shifted in time in such a way that the MNA continues increasing when *N*<sub>c</sub> started decreasing hundreds if not thousands of years before.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

(a) $t_{cont} = 340$ years

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

(b) $t_{cont} = 10,880$ years

Figure S5: Dynamics of the $N_c$, $H_e$ and MNA in spatial landscape and equivalent single panmictic population model: Two contraction durations are considered: $t_{cont} = 340$ years (Panel a) and $t_{cont} = 10,880$ years (Panel b). In each panel, three subpanels represent $N_c$ (A), $H_e$ (B), and MNA (C). In each subpanel, each discrete point represents the mean and the error bars represent the standard deviations measured across 15 simulations. The three vertical lines, from left to right, represent the following events respectively: i) start of the contraction phase (black, time = 2,860 years) and end of the contraction phases for (ii) $t_{cont} = 340$ years (red, time = 3,200 years) iii) $t_{cont} = 10,880$ years (green, time = 13,740 years). Dashed lines connect discrete time points, and symbols colors denote different population parameters for both the spatial stepping stone model and the equivalent panmictic model. Note that in the spatial landscape, genetic diversity was estimated based on 40 individuals sampled from the SW refugium (10 individuals/deme), whereas, in the single panmictic model, an equal number of individuals were sampled from the panmictic population.
<figure><img src="image_14.png" /><figcaption>**Figure S6: Life-history traits data points:** Relation between generation length (*GL*) and body mass (*BM*) (Panel A) and between dispersal distance (*DD*) and *BM* (Panel B) for mammalian species. In Panel A, the color code represents empirical data points from various orders of mammalian species, as reported in Pacifici et al. (2013). Black solid circles represent the parameter combinations selected for the spatial simulations (Table S1). In panel B, grey circles represent the empirical data points from Santini et al. (2013), while black circles represent the parameter combinations selected for the spatial simulations considering the allometric relationship between *DD* and *BM* (Table S1). The red triangles denote the carrying capacities and migration rates parameter values used in the simulations that did not specifically consider allometric relationships.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure S7: Contraction speed and genetic diversity at the end of contraction in the spatial simulations:** Relation between the contraction speed and the genetic diversity at the end of contraction in SW refugium at discrete times. Panel A shows expected heterozygosity (*H*<sub>e</sub>), while Panel B displays the *MNA*. Three generation lengths are considered: *GL* = 0.5 years (black circles), *GL* = 2 years (light-grey diamonds) and *GL* = 5 years (dark-grey triangles). The expansion and stationary phase were 360 and 2,500 years respectively, similar to Figure 1 B. The dotted (dashed) horizontal line represents the value of genetic diversity at mutation-migration-drift equilibrium when the refugium (entire landscape) is occupied. For this plot, the parameter values were (*K*,*m*)=(100,0.2).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure S8: Effect of contraction speed on genetic diversity at the end of contraction using *msprime* simulations: The genetic diversity at the end of contraction in a panmictic population as a function of contraction duration (x-axis, *t*<sub>cont</sub>) for three different mutation models, represented in three panels. Panel A illustrates Infinite Alleles Model (similar to Figure 5D), while Panel B displays Stepwise Mutation Model and Panel C represents Two Phase Model. In each panel, we plot two stationary times *t*<sub>stat</sub> = 250 generations (black) and *t*<sub>stat</sub> = 2500 (blue). Each data point represents mean across 100 replicates, whereas associated error bar indicate standard deviation. For this plot, the parameter values were (*N*<sub>0</sub>, *N*<sub>exp</sub>, mutation rate)=(400, 12000, 0.003).</figcaption></figure>
# S6 Demographic inference analysis using VarEff

## Methodology

Our simulations showed that during the process of habitat expansion and contraction, the census size ($N_c$) expanded and contracted by around two orders of magnitude, whereas genetic diversity ($H_e$ or $MNA$) followed a dynamics that could significantly differ from that of $N_c$. Here, we used the simulated microsatellite data and applied the method of Nikolic and Chevalet (2014) implemented in the VarEff R package. This method uses the coalescent theory to approximate the likelihood of the distribution of allele frequencies of microsatellites conditional on the mutation process and demographic model, assuming panmixia and population size changes. We used this method because it is less computationally demanding than other methods such as MsVar 1.3 (Storz and Beaumont, 2002), and because it allows for several $N_e$ changes. We ran VarEff for each of 15 simulations across various scenarios changing in contraction speed and life-history traits.

We wrote our own R script to convert the simulated genetic data from a genind object (in adegenet (Jombart, 2008), see Methods below) to a VarEff input format. The VarEff method requires two initial parameter estimates: $DMAX$ (the maximum distance between alleles) and $NBAR$ (global prior mean of effective size $N_e$). To estimate these parameters, we ran the $Theta()$ function in VarEff to obtain the distribution of allele frequencies (to estimate $DMAX$) and mean $N_e$. Each run of the $Theta()$ function yielded a $DMAX$ value and range of prior values for $N_e$. We took the mean of the prior $N_e$ range for each simulation. Thus, for 15 simulations, we had 15 values of $DMAX$ and 15 prior $N_e$. We then selected the median of $DMAX$ and prior $N_e$ values as representative for that particular scenario and executed the VarEff function with the following parameters:

<table><tr><td>NBLOC</td><td>50</td></tr><tr><td>JMAX</td><td>2</td></tr><tr><td>MODEL</td><td>S</td></tr><tr><td>MUTAT</td><td>0.001</td></tr><tr><td>VARP1*</td><td>3</td></tr><tr><td>VARP2*</td><td>3</td></tr><tr><td>RHOCORN*</td><td>0</td></tr><tr><td>NumberBatch*</td><td>10,000</td></tr><tr><td>LengthBatch*</td><td>10</td></tr><tr><td>SpaceBatch*</td><td>10</td></tr><tr><td>Diagonale*</td><td>0.5</td></tr><tr><td>DMAX</td><td>median over 15 simulations</td></tr><tr><td>NBAR</td><td>(global prior mean of N<sub>e</sub>) median over the 15 simulations</td></tr><tr><td>GBAR</td><td>(number of generations since the assumed origin of the population) depending upon the number of generations</td></tr></table>

The parameters marked with * are set based on the authors' recommendations. We set the number of times the effective size changes ($JMAX$) equal to 2 for all scenarios, regardless of contraction speed. This choice aligns with our simulation, which involved a population undergoing expansion, stationarity, and contraction, resulting in a simplified history with three effective population size estimates. Testing other values for $JMAX$ ($JMAX = 3,4$) did not significantly alter the shape of the inferred effective population sizes. This observation is consistent with previous studies, where variations in the $JMAX$ parameter had only a slight impact on the curve's steepness (Gonçalves et al., 2021; Lourenço et al., 2017). The parameters that differed across various scenarios were $NBAR$ (global prior mean of effective size $N_e$), $GBAR$ (number of generations since the assumed origin of the population), and $DMAX$ (the maximal distance between alleles).

For a given combination of life-history traits, we ran VarEff for each of the 15 simulations under the scenarios with contraction duration of 340 (fast contraction) and 10,880 years (slow contraction). In all instances, genetic data were sampled at 14,000 years, irrespective of contraction time (unless
126

127

128

129

130

131

132

133

134

otherwise stated). We focused on the median of posterior distribution since it was found to be the
most robust estimator (Nikolic and Chevalet, 2014). Consequently, we obtained 15 inferred curves
of effective population size change for each scenario. In each of the inferred $N_e$ plots, we selected
the median of inferred $N_e$ across 15 simulation for each time point, alongside the interquartile range
(25% to 75%, Figure S9). The inferred curves of effective population size ($N_e$) were plotted along
with the known dynamics of census population size change ($N_c$) and expected heterozygosity ($H_e$).
Each point in $N_c$ and $H_e$ represents the average value over 15 simulations. It is important to note
that, in contrast to other figures in this study, the time axis is reversed (ancient times to the right
and present time to the left) in the following figures and presented in years before the present.
<figure><img src="image_17.png" /><figcaption>Figure S9: Dynamics of the $N_c$, $H_e$ and inferred $N_e$ in the spatial landscape with sampling in SW refugium: Four different combinations of life-history parameters were considered (a) $GL = 0.5$ years, $K = 100$, $m = 0.1$ (b) $GL = 5$ years, $K = 100$, $m = 0.1$ (c) $GL = 5$ years, $K = 150$, $m = 0.059$ (d) $GL = 10$ years, $K = 1000$, $m = 0.003$. Each panel consists of three subpanels representing $N_c$ (top), $H_e$ (middle), and inferred $N_e$ (bottom). In $N_c$ panel, the three vertical lines from right to left, represent the following events: i) start of the contraction phase (black) and end of the contraction phase for (ii) $t_{cont} = 340$ years (red) iii) $t_{cont} = 10,880$ years (green). Additionally, the black horizontal line in the $N_c$ panel represents the maximum census size when the landscape is fully occupied. In $N_c$ and $H_e$ subpanels, each point represents the average value over 15 simulations. In the inferred $N_e$ subpanel, each point represents the median of inferred $N_e$ across 15 simulation, alongside the interquartile range (25% to 75%) across 15 simulations. The population was sampled at time = 14,000 years (as in the spatial case) after the start of the simulation. The time axis is reversed and presented in years before the sampling.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure S10: **Dynamics of the N_c, H_e and inferred N_e in equivalent panmictic model:** Three different combinations of life-history parameters were considered (a) GL = 0.5 years, K = 50, (b) GL = 5 years, K = 150 (c) GL = 10 years, K = 1000. Each panel consists of three subpanels representing N_c (top), H_e (middle), and inferred N_e (bottom). In N_c and H_e subpanels, each point represents the average value over 15 simulations. In the inferred N_e subpanel, each point represents the median of inferred N_e across 15 simulation, alongside the interquartile range (25% to 75%) across 15 simulations. The population was sampled at time = 14,000 years (as in the spatial case) after the start of the simulation. The time axis is reversed and presented in years before the sampling. For more details, cf. caption of S9.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure S11: Dynamics of the N_c, H_e and inferred N_e in spatial landscape with just SW refugium occupied (no expansion or contraction): Four different combinations of K and m values were considered. The GL considered here was 1 year/gen. The population was sampled at time = 14,000 generations after the start of the simulation. The time axis is reversed and presented in years before the sampling. For more details, cf. caption of S9.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure S12: Dynamics of the N_c, H_e and inferred N_e in fully occupied spatial landscape: Four different combinations of K and m values were considered. The GL considered here was 1 year/gen. The population was sampled at time = 14,000 generations after the start of the simulation. The time axis is reversed and presented in years before the sampling. For more details, cf. caption of S9.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure S13: **Dynamics of the N_c, H_e and inferred N_e in spatial landscape for GL = 0.5, K = 100, m = 0.1 from genetic data sampled at different time points:** Four different sampling time points were considered, a) 3,560 years, b) 8,700 year, c) 11,220 year and d) 14,000 year. The figure shows the trajectories for fast (in red, tcont = 340 years) and slow (in green, tcont = 10,880 years) contraction. For more details, cf. caption of S9. The inferred N_e for fast contraction reveals a bottleneck event across all sampling times, whereas slow contraction shows an initial demographic expansion which changes into a demographic decline at the latest sampling time. The time axis is reversed and presented in years before the sampling.</figcaption></figure>
<figure><img src="image_22.png" /></figure>

Figure S14: **Distribution of the frequencies of allele distances:** The two panels representing frequencies of allele distances from genetic data sampled at 14,000 years (forward in time) for two parameter combinations, namely a) GL = 0.5 years, K = 100, m = 0.1 and b) GL = 5 years, K = 100, m = 0.1 (from Figure S9a,b). Each panel displays two colored lines representing the fast (red, $t_{cont}$ = 340 years) and slow (green, $t_{cont}$ = 10880 years) contractions, respectively. Each point indicates the mean frequency across simulations, and the shaded region represents the standard deviation.

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

Figure S15: **Distribution of the frequencies of allele distances at different sampling times:** The two panels representing frequencies of allele distances for a) fast ($t_{cont}$ = 340 years) and b) slow ($t_{cont}$ = 10880 years) contraction scenarios at two different time points for GL = 0.5 years, K = 100, m = 0.1. In each panel, two colored lines represent the distribution of allele distances at two sampling times at 3, 560 years (light red, early sampling) and 14,000 (blue, late sampling), respectively. Each point indicates the mean frequency across simulations, and the shaded region represents the standard deviation.
<figure><img src="image_24.png" /></figure>

(a) $t_{cont} = 340$ yr, GL = 0.5 yr, K = 100, m = 0.1 (b) $t_{cont} = 10,880$ yr, GL = 0.5 yr, K = 100, m = 0.1

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

(c) $t_{cont} = 340$ yr, GL = 5 yr, K = 100, m = 0.1 (d) $t_{cont} = 10,880$ yr, GL = 5 yr, K = 100, m = 0.1

Figure S16: **Posterior Distribution of inferred $N_e$ with 90% confidence interval:** Two different values of GL while maintaining K = 100 and m = 0.1 were considered for both fast and slow contraction. Each panel represents the posterior distribution of inferred $N_e$ for three different simulations. Each point represents the median of the posterior distribution whereas the shaded region represents a 90% confidence interval of posterior distribution.
## Results

Figure S9 illustrates the parallel dynamics of the census population size ($N_c$), expected heterozygosity ($H_e$), and inferred $N_e$ across four different sets of life-history traits parameters. Overall, we observe that the demographic history ($N_e$) inferred by the coalescent-based approach VarEff is not coherent with the simulated $N_c$, but sometimes correlates with the $H_e$ dynamics (Figure S9b, c and d). For instance, the slow contraction scenario for $GL = 5$ years, and the fast and slow contraction scenarios for $GL = 10$ years show similarity in $H_e$ and $N_e$ trajectories, though timings of size changes differ. Specifically, for $GL = 5$ years (Figure S9b, c), slow contraction shows the ancient increase in $N_e$, capturing ancient $H_e$ trajectory. In contrast, the inferred $N_e$ for fast contraction appears constant over time, with VarEff being unable to detect the ancient population size increase. In this particular case, the constant inferred $N_e$ could be attributed to the rapid reduction of the population to refugia, where it remained unchanged for a long time. For $GL = 10$ years (Figure S9d), both fast and slow contractions reveal a steady increase in inferred $N_e$ over time, aligning with $H_e$ trajectories. Two rates of contraction can be distinguished from the inferred $N_e$, with slow contraction showing higher inferred $N_e$ compared to the fast contraction. However, for $GL = 0.5$ years (Figure S9a), the fast contraction surprisingly exhibits higher past $N_e$ than slow contraction, detecting a reduction in $N_e$ around the same time for both contraction speeds.

We further explored how the inferred $N_e$ varies with sampling time (Figure S13) for the scenario in Figure S9a, and found that average $N_e$ decreases with sampling times for slow contraction ($N_e \approx 30,000$ at 3,560 years versus $N_e \approx 7,500$ at 14,000 years) while increases for fast contraction ($N_e \approx 5,000$ at 3,560 years versus $N_e \approx 12,500$ at 14,000 years). While the inferred $N_e$ for fast contraction consistently shows a bottleneck event in the recent past across sampling times, slow contraction reveals a bottleneck only at a later sampling time (14,000 years).

To validate the inferred $N_e$ for different contraction speeds, we conducted the same VarEff analyses on various scenarios, including an equivalent panmictic population model (undergoing similar population size changes as in the spatial case, Figure S10), scenario where only the SW refugium was occupied (Figure S11) and a completely occupied landscape (Figure S12). In these scenarios, we started with the same initial condition (i.e. $MNA = 1$ and $H_e = 0$) and sampled genetic data at 14,000 years. In the equivalent panmictic population model, a consistent pattern is observed between two contraction durations, revealing that fast contractions exhibit lower inferred $N_e$ compared to the slow contraction in the past, across different combinations of life-history traits parameters. The fast contraction scenarios for $GL$ of 0.5 and 5 years, the inferred $N_e$ remains constant as the population quickly reduced to refugia and remained unchanged for long time. Conversely, the slow contraction scenarios exhibit variability for these two generation lengths. In the case of $GL = 0.5$ years, the inferred $N_e$ captures the decline in $H_e$ but misses the ancient increase. On the other hand, for $GL = 5$ years, the inferred $N_e$ shows both expansion and (slight) decline, though timings differ compared to $H_e$ and $N_c$. In the case of $GL = 10$ years (Figure S10c), the inferred $N_e$ follows the $H_e$, again differing in timings of increase and decrease.

In the spatial scenario depicted in Figure S12, where only the SW refugium is occupied for the entire simulation, the inferred $N_e$ consistently exhibits a pattern of constant values in the past, aligning with the trajectories of $N_c$ and $H_e$. Moreover, the inferred $N_e$ values are close to the actual $N_c$. In the fully occupied landscape (Figure S12), we observe different patterns across different combinations of $K$ and $m$ values. In Figure S12b, c, d, the inferred $N_e$ trajectories show an ancient expansion (aligning with $H_e$ trajectory). The recent inferred $N_e$ coincided closely with $N_c$ for $K = 100, m = 0.1$, whereas not in the other cases. For the particular case of $K = 50, m = 0.1$, we observe constant inferred $N_e$ in the past with a recent bottleneck. This pattern of recent bottleneck despite a constant population size (spatially structured) could be related to the issue of interpreting $N_e$ as population size change when fitting an inference method assuming panmixia to a structured population model (Wakeley, 1999; Beaumont, 2004; Chikhi et al., 2010).
Discussion

Our results reveal distinct results regarding the correspondence between the inferred effective pop-
ulation size ($N_e$) and the simulated census size ($N_c$) depending on the simulation scenario. For
the simplest cases (for instance, panmictic population and constant population sizes in structured
populations), the inferred $N_e$ trajectory correlated with simulated $N_c$ and $H_e$ trajectories (Figure
S11, S10). However, the simulated scenarios of interest involving spatial expansion and contraction,
led to infer $N_e$ trajectories that were decoupled with the simulated $N_c$ history, while it could un-
der some scenarios show broad similarities with $H_e$ trajectories (Figure S9). Interestingly, fast and
slow contraction simulated scenarios were distinguishable from the inferred $N_e$ trajectories for some
scenarios (Figure S9). Surprisingly, in one case (Figure S9a), we found that the inferred $N_e$ of fast
contraction was higher than for slow contraction, contrary to expectations given the greater loss of
genetic diversity experienced by fast contraction in the recent past. In this particular scenario, the
distribution of allele distances showed very slight differences between fast and slow contraction when
sampling at time = 14,000 years (Figure S14a).

We also observed that, for fast contractions (Figure S9a), the median $N_e$ across simulations
in the past tended to increase with increasing sampling time backward in time (Figure S13, red).
Conversely, slow contraction exhibited a transition from an inferred demographic expansion (when
sampling at 3,560 years) to an inferred bottleneck (when sampling at 14,000 years) with increasing
sampling time (Figure S13, green). The fast contraction scenario showed a wider distribution of
allele distances (Figure S15a) whereas slow contraction shows a consistent cumulative decrease in the
frequency of divergent alleles with sampling time, with a significant increase in identical alleles (allele
distance = 0 in Figure S15b for 14,000 year). These divergent alleles may contribute to the higher
inferred $N_e$, as these divergent alleles would take a relatively long time before coalescing, having a
lower coalescence rate, which is inversely proportional to $N_e$ (Kingman, 1982; Hudson, 1983). On
the other hand, the posterior distribution of inferred $N_e$ along with a 90% confidence interval for
sampling at 14,000 years shows a wide confidence interval in the case of both $GL = 0.5$ and $GL = 5$
years (Figure S16). Across simulations, the medians of the posterior were significantly different for
$GL = 0.5$ years, while they were more similar for $GL = 5$ years. Given these observations, explaining
the higher inferred $N_e$ for fast contraction compared to slow contraction would require more detailed
analysis.

Additionally, across the simulated scenarios, we observe that the inferred $N_e$ reflected relatively
recent part of the history rather than the full genetic history (i.e., $H_e$ trajectory). This becomes
evident when comparing the demographic history inferred at early versus late sampling times. For
instance, in a slow contraction scenario, inferred $N_e$ identifies a demographic expansion from genetic
data sampled at 3,560 year, which changes to a constant $N_e$ at 11,220 year and subsequently to
a recent bottleneck at 14,000 year (Figure S13). This suggests that information on more ancient
demographic events is progressively lost with increasing sampling time. Similar issues have been
previously identified with other coalescence-based approaches applied to mtDNA (Grant et al., 2012;
Grant and Cheng, 2012) and genomic data (Postaire et al., 2023), showing that ancient population
history tends to get hidden by more recent demographic events.

Our analyses, though limited, indicate that even in simple scenarios of population expansion
and subsequent contraction, trajectories of inferred $N_e$ should be interpreted with a grain of salt,
as they may provide little information about $N_c$ and $H_e$ changes through time. Life-history traits
play a significant role in deviations from genetic equilibrium as well as in generating different levels
of population structure across multiple species, which could lead to discordant $N_e$ even though
they underwent similar changes in habitats. Thus more work is required to understand the relation
between inferred $N_e$, $H_e$, and $N_c$ in cases where the populations deviate from genetic equilibrium.