Problems seeded in the past: lagged effects of historical land-use changes can cause an extinction debt in long-lived species due to movement limitation

Land-use change is one of the main threats to biodiversity on the global scale. Legacy effects of historical land-use changes may affect population dynamics of long-lived species, but they are difficult to evaluate through observational studies alone. We present here an interdisciplinary modelling approach as an alternative to address this problem in landscape ecology. Assess effects of agricultural abandonment and anthropisation on the population dynamics of long-lived species. Specifically, we evaluated: (a) how changes in movement patterns caused by land-use change might impact population dynamics; (b) time-lag responses of demographic variables in relation to land-use changes. We applied an individual-based and spatial-explicit simulation model of the spur-tighed tortoise (Testudo graeca), an endangered species, to sequences of real-world landscape changes representing agricultural abandonment and anthropisation at the local scale. We analysed different demographic variables and compared an “impact scenario” (i.e., historical land-use changes) with a “control scenario” (no land-use changes). While agricultural abandonment did not lead to relevant changes in demographic variables, anthropisation negatively affected the reproductive rate, population density and the extinction probability with time-lag responses of 20, 30 and 130 years, respectively, and caused an extinction debt of 22%. We provide an understanding of how changes in animal movement driven by land-use changes can translate into lagged impacts on demography and, ultimately, on population viability. Implementation of proactive mitigation management are needed to promote landscape connectivity, especially for long-lived species for which first signatures of an extinction debt may arise only after decades.


Abstract
Context Land-use change is one of the main threats to biodiversity on the global scale. Legacy effects of historical land-use changes may affect population dynamics of long-lived species, but they are difficult to evaluate through observational studies alone. We present here an interdisciplinary modelling approach as an alternative to address this problem in landscape ecology.
Objectives Assess effects of agricultural abandonment and anthropisation on the population dynamics of long-lived species. Specifically, we evaluated: (a) how changes in movement patterns caused by land-use change might impact population dynamics; (b) time-lag responses of demographic variables in relation to land-use changes.
Methods We applied an individual-based and spatial-explicit simulation model of the spur-tighed tortoise (Testudo graeca), an endangered species, to sequences of real-world landscape changes representing agricultural abandonment and anthropisation at the local scale. We analysed different demographic

Introduction
Land-use change is one of the main threats to biodiversity worldwide, and together with climate change a main driver of global change (Sala et al. 2000;Daskalova et al. 2020). Land-use change does not follow a single trajectory, but many and, possibly, divergent processes of change operate with different spatial-temporal dynamics depending on the history of the landscape, socio-economic conditions and ecological context (Foley et al. 2005). On the one hand, traditional agricultural practices in rural landscapes are abandoned (hereafter agricultural abandonment) and lead to the recovery of natural or seminatural landscapes (Moreira and Russo 2007;Queiroz et al. 2014;Vadell et al. 2016). On the other hand, anthropisation processes lead to deforestation, agricultural intensification, development of linear infrastructures, tourist and industrial developments as well as urban areas (Tzanopoulos and Vogiatzakis 2011;Laurance et al. 2014;Martínez-Fernández et al. 2015).
The consequences of these land-use changes for biodiversity are complex (Daskalova et al. 2020). Agricultural abandonment, which has accounted for about 13% of global land-use changes since 1992 (Haščič and Mackie 2018), leads to the regeneration of forests and other natural habitats through passive management approaches (Navarro and Pereira 2012), and to the recovery of ecological functions (Moreira and Russo 2007;Navarro and Pereira 2012;Queiroz et al. 2014). In contrast, the effect of intensified land use, which has accounted for approximately 41% of global land-use changes since 1992 (Haščič and Mackie 2018), impoverishes biological communities (Sala et al. 2000) and brings about loss of ecological functions and ecosystem services (Foley et al. 2005;Maes et al. 2012). In the Mediterranean Basin, where land uses are distributed in a fine-grained manner, current landscapes are strongly polarised with abandoned areas often directly adjacent to intensified areas (Tzanopoulos and Vogiatzakis 2011).
An important task for the conservation of landscapes with such mosaic structures is to better understand how animal species may respond to abandonment and intensification. Animals typically alter their movement behaviour in response to landuse changes because they need to adapt to changes in resource availability and landscape structure (Lange et al. 2013). Interestingly, this response is highly idiosyncratic. Species may show less movement if they are reluctant to cross patch boundaries (Kuefler et al. 2010), such as infrastructures blocking movement, which can lead to isolation (Soanes et al. 2013;Lister et al. 2015). In the other extreme, species may show increased movement if they must forage across many distant patches for resources that are left after habitat loss (Andreassen and Ims 1998;Saunders 1980). These contrasting responses of movement to land-use change may lead to strong differences in the resulting population dynamics. However, this critical aspect of land-use change has hardly been explored. Most studies that investigate how land-use changes impact population dynamics are based on effects mediated by direct (i.e., imposed) changes in reproductive and survival rates, whereas the impact of disrupting movement patterns is often neglected (Doherty and Driscoll 2017). The present study, therefore, makes inferences about the effects of landuse change on different aspects of population dynamics and viability by accounting for observed changes in movement behaviour with land use (Anadón et al. 2012). It does so by considering the dynamics of two contrasting land-use changes that dominate many landscapes: abandonment and anthropisation processes.
Whereas the dynamics of some species might rapidly respond to landscape changes (e.g. Alcocer-Rodríguez et al. 2020), this response may often show time-lags, especially in long-lived species (Morris et al. 2008;Metzger et al. 2009). Time-lag response can generate two types of land-use legacies (Lira et al. 2012): species credit and extinction debt, which refer to species that may become more or less abundant following the improvement or worsening of habitat quality, respectively (Hanski 2000;Daskalova et al. 2020). Land-use legacies have mainly been analysed in the context of extinction debt responses in communities (Figueiredo et al. 2019), and very few studies have analysed the extinction debt of single species (see Pandit et al. 2017;Sergio et al. 2019;Alcocer-Rodríguez et al. 2020) and mostly in a metapopulation context (Schtickzelle et al. 2005).
The potentially lagged response to land-use changes makes observational studies unfeasible for long-living species (Lindenmayer et al. 2012;Ewers et al. 2013). An alternative to observational studies is use of the interdisciplinary approach of computer simulation modelling (Stillman et al. 2015). For example, individual-based simulation models (IBM; Wiegand et al. 2004;Grimm and Railsback 2005;Grimm et al. 2006;Jiménez-Franco et al. 2020a) can integrate the knowledge and field data that are available on the spatially explicit movement of individuals and demographic processes into a common framework. IBMs are effective tools in assessing emerging spatial effects (DeAngelis and Grimm 2014), such as that of mate-finding Allee effects (Jiménez-Franco et al. 2020a), while the effect of low survival or reproduction rates on population dynamics can be assessed well with age-structured matrix models (Caswell 2001;Graciá et al. 2020a). The IBMs developed with specific management goals in mind usually require a high degree of detail, which is essential for making predictions, but IBMs are ideally suited to obtain a more principle understanding of the complex phenomena emerging in landscape ecology, such as lagged demographic responses as a consequence of altered movement behaviour in dynamic landscapes (DeAngelis and Grimm 2014). Here the structural (spatial) realism of IBMs allows the essential spatial relations that are important for the question under study to be included, while hypothetical simulation experiments may manipulate the underlying landscape, initial population densities or other important aspects for the question. By taking the latter approach, we can assess the land-use legacies of agricultural abandonment and anthropisation by comparing the development of a simulated population in landscape change scenarios to that of the initial landscape without changes.
The aim of this study is to better understand impacts of land-use changes (agricultural abandonment and anthropisation) on the population dynamics of longlived species inhabiting fine-grained cultural landscapes. Specifically, we used an individual-based modelling approach to quantify: (a) how changes in movement patterns caused by land-use change might translate into impacts on different aspects of population dynamics; (b) to what extent the response of different aspects of species' demography are lagged in relation to such land-use changes. To address our general objective, we used the spatially-explicit model STEPLAND (Graciá et al. 2020a;Jiménez-Franco et al. 2020a) that is based on data from intense field studies in southeast Spain of the spatial and demographic processes of the spur-thighed tortoise (Testudo graeca). This species is ideally suited to address our general objective because it is a long-lived species inhabiting an area subject to both, agricultural abandonment and anthropisation, it is threatened by landuse change (Anadón et al. 2007;Graciá et al. 2020b), it shows strong movement responses to land-use change (Anadón et al. 2012), and its spatial and demographic processes are well studied (e.g., Pérez et al. 2002;Rodríguez-Caro et al. 2013. We hypothesise that changes in movement patterns due to land-use change might interact with a matefinding Allee effect (Graciá et al. 2020a;Jiménez-Franco et al. 2020a) that could result in either a negative impact on reproduction due to fewer mating opportunities in fragmented landscapes, or in a positive effect due to rising reproduction rates if increased movement distances lead to more mating opportunities.

Model species: the spur-thighed tortoise in SE Spain
The spur-thighed tortoise T. graeca is a medium-sized long-lived tortoise with slow population dynamics, steady population growth rates (Rodríguez-Caro et al. 2017, delayed maturation (9-12 years; Rodríguez-Caro et al. 2013), and a clutch size of 1-7 eggs. It is endangered (Vulnerable; IUCN 1996), with main threats in Spain being habitat loss and fragmentation caused by anthropogenic activities (e.g. irrigated lands, greenhouses and tourist urban development; Anadón et al. 2007;Graciá et al. 2020b). In Spain T. graeca inhabits heterogeneous cultural landscapes in the coastal mountains of the most arid region of SW Europe, formed by a shrubland matrix that includes other uses, mainly irrigated and non-irrigated crops (Anadón et al. 2006b). Since the 1950s, this area has undergone intense rural depopulation (De Aranzabal et al. 2008), and since the 1980s the region has been subject to agricultural intensification, the development of urban and touristic areas and infrastructure development (Martínez-Fernández et al. 2015). During the last 60 years these processes have created a polarised landscape where natural and anthropisated areas cooccur, often in close proximity, in an area traditionally inhabited by this tortoise.

General rationale
We used the spatially explicit individual-based model STEPLAND ( Fig. 1), which synthesises data of intense field studies from almost two decades to simulate the species demography, habitat selection and movement of T. graeca in real landscapes. We simulated the population dynamics of T. graeca from the 1950s to the present-day based on time-series data of real landscapes that represent historical land-use changes. We also projected the population trend into the future and recorded changes in several emerging demographic characteristics, including the reproductive rate (RR; i.e., mean number of offspring of reproductive females), population density (number of individuals per hectare), and extinction probabilities after 200 years (P ext ). To check the legacy effects of historical land-use changes on population dynamics, we also quantified the extinction debt.
Habitat preference and changed movement patterns impact population dynamics in STEPLAND through a mate-finding Allee effect (Jiménez-Franco et al. 2020a) in two main ways. According to Anadón et al. (2012), tortoises, especially females, respond to a high proportion of intensified areas by moving longer distances and by shifting from a narrow home range to a more unbounded movement pattern, which may lead to increased mating opportunities. Secondly, tortoises very much avoid areas with intensified use (i.e., main roads, dense human infrastructures, intensive agriculture or pine forests) (Anadón et al. 2012), which act as barriers that decrease mating opportunities. However, our extensive empirical work found no evidence for variation in survival (based on age-structure analyses) and clutching between natural and altered landscapes (own data), probably because tortoises never used areas of intensified use (see Table C3 in Anadón et al. 2012). For this reason, and as our goal was to assess the potential impact of land-use changes on population dynamics based on the observed movement patterns, we used the same survival rates and reproduction parameters (i.e., number of clutches, number of eggs per clutch) for the tortoises inhabiting natural and altered landscapes, but different movement parameterizations. Decreasing survival and reproduction will lead to a lower population growth rate, which can be assessed using Leslie matrix models ( The STEPLAND model ( Fig. 1) was developed to integrate demographic processes, i.e., reproduction, mortality and ageing (Graciá et al. 2020a;Jiménez-Franco et al. 2020a), with a previously developed individual-based model of tortoise movement (Anadón et al. 2012). For model description, we followed the Overview, Design concepts and Details protocol (ODD) proposed by Grimm et al. (2010). We summarised the model in the paragraphs below, while a full ODD version is provided in Appendix S1. The model was implemented in Python 2.7, and its code, parameterisation, and the main result files are available in the Figshare repository (Jiménez-Franco et al. 2021). started with an initial population distributed over the given landscape (Los Mayorales or Los Estrechos) and was governed by a specific model parameterisation. After parameter uploading, tortoises were subject to the three basic biological processes of movement, reproduction and mortality (see dashed arrows; details in Appendix S1). Landscape characteristics influence tortoise movement and mating (see solid arrows): individual tortoises obtained by birth one movement parameterisation, that was randomly assigned from a pool of 5627 parameterisations, depending on their sex and the landscape type. The movement parameterisations governs the three subprocesses homing, habitat dependency and autocorrelation (Anadón et al. 2012). Reproduction (specifically mating) was conditioned by a maximal encounter distance between males and females and landscape barriers (Graciá et al. 2020a;Jiménez-Franco et al. 2020a). Clutching and mortality, were parameterised based on empirical data and literature (see Appendix S3). The model output comprised the detailed information about the simulated tortoise populations every 10 years that were used to estimate derived population variables such as density, reproductive rates, or extinction probabilities STEPLAND contains two types of entities: landscape and tortoises. Landscapes are composed of a grid of 10 m 9 10 m cells. Each grid cell is characterised by its position (x and y coordinates), and also by its assignment to one of the five habitat quality parameters in STEPLAND. We reclassified the original land-use map that comprised eight land-use cover categories to five habitat types that influence movement (see more details of landscape reclassification in Appendix S2). Tortoises have the attributes sex, age and their location over time, and the model is composed of specific submodels for movement, reproduction and mortality (Fig. 1). The overall time step of the model is 1 year, annual mortality occurs at the end of the year before the age of the surviving tortoises is updated, and newborn tortoise are added at the beginning of the year. However, the movement submodel operates at a finer intraday scale with sexspecific movement probabilities depending on the month.
The movement submodel reflects the empirical finding that the activity and movement patterns of T. graeca individuals strongly differed between ''natural'' and ''altered'' landscapes (for a definition of landscape types see the following subsection) as found by Anadón et al. (2012). Movement tracks of tortoise consisted in the model of a sequence of steps to neighbour 10 m 9 10 m cells, with probabilities to enter one of the eight neighbouring cells, or to stay in the current cell, depending on three sets of weights that described autocorrelation in movement, home behaviour, and habitat dependence. Each day was divided into four activity intervals, and a probability PMOV decided if a tortoise moved during that interval, and a probability distribution DMOV decided how many steps the individual moved. PMOV and DMOV were determined from the radio-tracking data (see Fig. S4 in Appendix S1), and PMOV followed a seasonal pattern driven by hibernation and aestivation periods, and differed between males and females. To reflect the individual variability in the movement patterns, Anadón et al. (2012) determined in a simulation study 5627 individual movement profiles, which included four types of profiles (female-natural, female-altered, male-natural, and male-altered). They consisted of different combinations of the nine movement parameters that produced trajectories that were in agreement with that observed by radio-tracking data. Following previous studies (Graciá et al. 2020a;Jiménez-Franco et al. 2020a), we randomly assigned one individual movement profile to a tortoise at birth and for all other tortoises we assigned a new movement profile at the end of each year (corresponding to its sex and landscape type). Further details are provided in the ODD protocol in Appendix S1 (Movement section).
The implementation of the demographic processes of reproduction and mortality reflect the species biology with most demographic parameters determined directly from our data (see Tables S1 and S3 in Appendix S3). Reproduction consists of the subprocesses mating, sperm storage and clutching. Mating probabilities in T. graeca emerge from the interaction of movement with landscape structure and are strongly conditioned by landscape connectivity (Fig. 1). STEP-LAND assumed that mature females will mate: (i) if during the peak mating season (end of April) at least one reproductive male is located within an encounter distance of 500 m (as determined by Graciá et al. 2020a); and (ii) if females and males were not separated by barriers such as main roads, dense human infrastructures or intensive agriculture (see landscape barriers in Fig. S2, Appendix S2). Sperm storage of T. graeca females can allow them to reproduce for the next 3 years after mating (Cutuli et al. 2013;Jiménez-Franco et al. 2020a). Adult females harboring viable sperm have the opportunity to lay between one to three clutches per year (with 1-6 eggs per clutching), as observed in real populations ). The eggs were placed at the location of the female in spring and early summer when clutching occurs (from the end of April to the end of July), and the newborns emerge from this location.
Annual survival rates are applied at the end of the year and vary among age classes, including newborns (representing hatching success and survivals of individuals under 1 year), immature individuals (aged 1-3), subadults (aged 4-6) and adults (aged C 7). We did not include a carrying capacity or negative density dependence in demographic rates because the population sizes occurring during model simulations ranged between 0 and 0.4 ind/ha, which correspond to ''very low'' and ''low'' densities, considering that the highest density observed in the study area was 13 ind/ha (Anadón et al. 2009).
Landscape scenarios: historical habitat transformations in the Mediterranean ecosystem The two selected study areas are located in the Almenara mountain range at the core of the species' range distribution area of T. graeca in SE Spain (Anadón et al. 2006a(Anadón et al. , 2007 (Fig. 2a) and faithfully represent the history of land-use changes of semiarid Mediterranean landscapes in the past century and, to a certain extent, the history of many rural landscapes in developed countries. The two modelling areas covered 5 9 5 km (Fig. 2b, c). The historical dynamics of land uses were characterised by mapping land-use covers in 1956, 1987 and 2010. One of the study areas, Los Mayorales, has a representative history of the agricultural abandonment with secondary succession of increased Mediterranean scrub. This landscape is composed of a mosaic of natural areas (Mediterranean Fig. 2 General framework to study the effects of landscape evolution on Testudo graeca populations by using the individual-based model STEPLAND. a Range of distribution of the species in the study area (SE Spain, in grey) and location of the two study areas of 5 9 5 km. b Overview of the spatiotemporal land-use scenarios. Black arrows denote the temporal period from 1956 up to 200 simulation years; red arrows depict the two landscape change times, 1987 and 2010, respectively, and the final simulation period (in the year 2156). The two study areas have land-use scenarios: ''without change'' and ''with change''. The ''without change'' land-use scenarios assume that no changes occurred at all throughout the simulation period (control scenarios). In the mid-twentieth century, the two areas comprised a mosaic of natural vegetation and traditional rainfall agriculture. The ''with change'' land-use scenarios assume real changes for each landscape areas undergoing a land-use polarisation process: on the left, study area Los Mayorales has undergone the agricultural abandonment process; on the right, study area Los Estrechos has undergone an anthropisation process (i.e. increased intensive agriculture, linear infrastructures, urban areas). c Percentage of land-use cover of the landscape scenarios at three historic times for the two study areas of 5 9 5 km scrub and pine forest) and traditional agricultural areas (dry-land crops such as cereals, almonds and olives). The other study area, Los Estrechos, represents the anthropisation process, that is, an increase in agricultural intensification, urban areas and infrastructure development with highways (which are no permeable for T. graeca) and roads (which have low permeability).
We built six landscape maps (2 study areas 9 3 periods) with eight land-use cover categories (Fig. 2b,  c), which were transformed to five habitat types that influence movement (see more details of the mapping process in Appendix S2). We selected these habitat types because they influence T. graeca movement patterns (Anadón et al. 2012; see the ODD for details). Each habitat type is associated with one parameter (i.e., H1W, H2W, …, H5W) that gives the weight in the habitat-biased random walk. The habitat types include (i) with null habitat quality (H5W): urban areas and highways (which are never used); (ii) with very low habitat quality (H1W): intensive agriculture (irrigated crops, citrics and greenhouses), pine forest areas (Mediterranean pine forest of Pinus halepensis) and roads (paved roads except highways); (iii) with good habitat quality (H4W, H3W and H2W, respectively): natural areas on slope (Mediterranean scrub on slope), flat natural areas (flat Mediterranean scrub) and traditional agricultural lands (see more details in Appendix S2).
As the movement parameterizations differ between natural and altered landscapes, we need to assign each of the six landscape maps either to natural or to altered, based on their habitat composition (Anadón et al. 2012). The 2010 landscape in the study area Los Estrechos was classified as altered, whereas all other five landscapes were classified as natural. A landscape with more than 23.8% of cover types (i) and (ii) defined above was defined as altered and natural otherwise.

Design of simulation experiments
We used initial population sizes of 250 individuals in a 5 9 5 km area that represented low T. graeca population densities in natural Spanish populations (0.1 tortoises/ha; Anadón et al. 2009). We selected this density because previous work showed that it is related to an extinction threshold (by assuming sperm storage durability of 3 years): if the initial density is below, spatial effects lead to a mate-finding Allee effect that can cause, especially in altered landscapes, growth rates below 1, but otherwise growth rates exceed 1 (Jiménez-Franco et al. 2020a). Thus, for higher initial densities we would most likely not observe responses to altered movement patterns.
As in previous parameterisations (Graciá et al. 2020a;Jiménez-Franco et al. 2020a), the initial age structure was the same for all the simulations. It was based on a stable age distribution predicted by an agestage structured deterministic matrix population model using the same parameterisation of survival and reproduction parameters as in the simulation model. For this purpose, we employed the popbio package (Stubben and Milligan 2007) in the statistical R programme, version 3.5.1 (R Core Team 2021; R scripts are shown in Appendix S4). To avoid individuals' unrealistic initial distribution (e.g., too many individuals placed in low suitability areas), we allowed the initial population to adjust to habitat types. We accomplished this by initially distributing tortoises randomly across all areas of good habitat quality (i.e., habitat types H2, H3, and H4). We then simulated the movement of these individuals (but without demography; i.e., reproduction, mortality or ageing) over 50 years to let them adjust to the initial landscape map. We then took these individual locations as the starting position in the following simulations for the two selected areas with land-use scenarios ''without''/''with'' change since 1956 (Fig. 2).
We assessed the effect of both, the agricultural abandonment process and the anthropisation process on the reproductive rate (RR), population density (ind/ ha) and extinction probability (P ext ) of T. graeca over 200 years. To this end, we conducted ''impact'' scenarios where we simulated the population dynamics in the two study areas: Los Mayorales (representing agricultural abandonment) and Los Estrechos (representing anthropisation). We also conducted ''control'' simulations by assuming that the landscape had not changed since 1956 (Fig. 2). To account for the stochastic processes, especially for estimating extinction probabilities, we repeated each scenario 256 times. Therefore, the total number of independent model simulations was 1024 (2 study areas 9 2 landuse scenarios 9 256 replicates).

Data analysis
The simulations in each landscape scenario were analysed to calculate three demographic variables, which were averaged over the 256 replicate trajectories every 10 years. First, the reproductive rate RR at time step t was the mean number of offspring of reproductive females, calculated by dividing the number of individuals born in time step t by the number of reproductive females (those aged 10 years or older) that were present in time step t. Only the nonextinct trajectories in each time step were used to calculate RR. As a second demographic variable we used the population density at time step t. Finally, we calculated the extinction rate P ext for each year t as the proportion of extinct replicates (i.e., having a population size N t B 1).
To quantify the total legacy effects of the historical land-use changes at the end of the simulation period, we calculated the differences in population densities and P ext -values between the corresponding control vs. impact scenarios. The extinction debt for a given scenario was calculated as the difference between the population viability after the habitat change and the end of the simulation period (Kuussaari et al. 2009). The time-lag for RR and population density was quantified by considering the initial period of the landuse change (in 1987) until the period when the confidence intervals of the control/impact scenarios did not overlap. The time-lag for population viability was quantified by considering a difference in the control/impact scenarios values equal or higher than 0.05.

Results
The two historical study areas showed initially (i.e., in 1956) similar characteristics, being typical mosaic structures of cultural and historical Mediterranean landscapes composed mainly of Mediterranean scrub and traditional agricultural land with a 72:28 proportion (Los Mayorales) and a 66:32 proportion (Los Estrechos; Fig. 2c), respectively. Other land uses (urban, intensive agriculture, linear infrastructures) were practically non-existent (0.5% and 1.9%, respectively; Fig. 2c).
The results of the simulations in the control scenarios with the land-use cover of 1956 showed for both study areas similar mean reproductive rates over time, with similar values at the end of the simulation period: 1.61 (1.52-1.70, CI) and 1.50 (1.41-1.59, CI) for Los Mayorales and Los Estrechos, respectively (Fig. 3a, b). Population densities increased over time in both study areas (Fig. 3c, d; see lines with triangle symbols) and reached at the end of the simulation period in 2156 a mean of 0.34 ind/ha (Los Mayorales) and 0.35 ind/ha (Los Estrechos) (Table S1), with an extinction probability of 0.07 for both study areas at the end of the simulation period (Fig. 3e, f; Table S1).
Landscape transformations over time  caused polarised landscapes (Fig. 2b). Los Mayorales underwent a crop-scrub transition with 20% traditional agricultural land loss and a 17.8% increase in natural areas ( Fig. 2c; see Appendix S2 for details). In Los Estrechos, changes were more pronounced with decreases of 19.3% in traditional agricultural land and 17.5% in natural areas, but increases of 31% for intensive agriculture, 8.3% for urban areas, and 1.6% for roads and highways (Fig. 2c, see Appendix S2 for details).
The results of the simulations of the ''impact'' scenarios showed for both study areas clear differences in the demographic variables over time, with opposite tendencies for reproductive rates (Fig. 3a, b) and population densities (Fig. 3c, d; see the lines with circle symbols). In Los Mayorales, which underwent agricultural abandonment, tortoise density increased over time (Fig. 3c) with a mean density of 0.31 ind/ha at the end of the simulation period (Table S1). On the contrary in Los Estrechos, which underwent an anthropisation process, tortoise density decreased over time (Fig. 3d) reaching a mean density of 0.04 ind/ha at the end of the simulation period (Table S1). The time-lags in response to historical land-use change for Los Estrechos differed among demographic variables ( Fig. 3; see red arrows), however, effects of the two land-use changes were difficult to separate. Changes in reproductive rate were immediately visible with a time-lag of 20 years and occurred already before the second land-use change that amplified the effect (Fig. 3b). Population density responded slower, larger effects were not imminent before 2016 and resulted in a time-lag response of approximately 30 years (Fig. 3d). Finally, effects on the extinction risk appeared only around year 2116, 130 years after the beginning of land-use change (Fig. 3f).  Fig. 2). a, b Reproductive rate (RR, mean number of offspring per reproductive female ± confidence intervals). Note that the RR at time 0 is the deterministic RR det = 1.7. Horizontal lines represent the critical RR of T. graeca (1.35) that would result from a critical growth rate of 1. c, d Mean population density (individuals per hectare, mean ± confidence intervals), with 0.1 at time 0 for all the scenarios. e, f Population viability (1 -proportion of extinct replicates). The arrows in panels e and f depict the extinction debt (calculated as the difference between the population viability after the habitat change and the end of the simulation period in the scenario ''with change''). Each simulation was repeated 256 times. Triangles (filled triangle) denote the ''without change'' land-use scenarios (i.e., control scenarios, with no change in land-use cover during the simulation period); circles (filled circle) depict the ''with change'' land-use scenarios, i.e. impact scenarios, with a change in land-use cover during the simulation period (agricultural abandonment process and anthropisation process for each study area, respectively). Vertical dashed black lines represent the two land-use change periods of the ''with change '' scenarios in 1987 and 2010 In the Los Mayorales study area that represents agricultural abandonment, land-use change did not cause differences in tortoise density (Figs. 3, 4a; Table S1) and extinction rates (Fig. 4b). In contrast, in the Los Estrechos study area the demographic variables differed strongly between the control and impact scenario (Fig. 4), with both reproductive rates (Fig. 3b) and individual density (Fig. 3d) being considerably lowered. This resulted at the end of the simulation period in a decrease in population density of 0.31 ind/ha (Table S1), and an increased in the extinction probability by 314% (Table S1). Finally, the land-use legacy effects resulted in Los Estrechos in an extinction debt of 22% additional extinctions at the end of the simulation period (Fig. 3f, red arrow), but only a negligible extinction debt resulted in Los Mayorales (0.09; Fig. 3e, red arrow).

Discussion
This study demonstrates the applicability of individual-based simulation models (IBMs) to provide suitable solutions to practical problems in landscape ecology, such as assessing emergent spatial and demographic phenomena that are difficult to capture otherwise. We show how an IBM can be used as effective tool to assess the lagged effects of historical land-use changes on species' population dynamics, an endeavour which would not be possible with empirical field monitoring (Lindenmayer et al. 2012;Stillman et al. 2015;Jiménez-Franco et al. 2020b). Our interdisciplinary approach allowed us to analyse and better understand the complex ways of how land-use change can impact population viability and create an extinction debt, as mediated by changes in movement patterns, and as a result of reproduction and population dynamics, even without assuming direct land-use driven changes in reproduction and survival rates.
We generally found a strong impact of the anthropisation process on the target species' population dynamics. Remarkably, this impact showed a delay of several decades, which indicates a marked extinction debt in these landscapes (Krauss et al. 2010). In contrast, agricultural abandonment did not lead to changes in population dynamics. However, our results cannot be taken literally as an inference of what is happening in the real landscapes, for example because the initial population densities for the year 1956, which critically influence the occurrence of the Allee effect (Graciá et al. 2020a;Jiménez-Franco et al. 2020a), are unknown. However, as aimed, our study provided a general understanding of effects in the population dynamics of species inhabiting fine-grain mosaic areas that can emerge under ongoing land-use polarisation (Navarro and Pereira 2012), which are Fig. 4 Effects of land-use changes on the demographic dynamics of Testudo graeca in the two landscapes zones of 5 9 5 km during the simulated period (200 year). The relation between the values of the demographic parameters for the landuse scenarios ''without change'' (x-axis) and ''with change'' (yaxis) are displayed (see details of these process in Fig. 2): a Population density (mean ± confidence intervals), b Extinction probability. Dots depict the values of demographic parameters during the simulation period (from 0 to 200 years, every 10 years); stars (*) depict the initial period and arrows represent trends over time. Colours represent study areas: Los Mayorales (green), Los Estrechos (blue). Diagonal dashed lines comprise values where the demographic parameters of the ''with change'' land-use scenarios do not differ from the parameter values of the ''without change'' land-use scenarios. Note that axes are on a logarithmic scale particularly helpful for the management of species with limited movement ability such as T. graeca.

Legacy effects of anthropisation
In our model, landscape anthropisation had a negative effect on the population dynamics of T. graeca in the three studied demographic parameters (i.e. reproductive rate, population density, population viability). Thus, increased movement distances due to land-use change could not compensate for the negative effects of habitat loss and fragmentation. This type of altered landscape likely reduces mate-finding between reproductive individuals (Gascoigne et al. 2009) and lowers the reproductive rate of females (RR) every year. Consequently, population density may decline below the thresholds that generate a mate-finding Allee effect, which impacts population viability (Graciá et al. 2020a;Jiménez-Franco et al. 2020a).
Remarkably, our results showed that anthropisation effects on the population viability of T. graeca appear only at large temporal scales. In the study area Los Estrechos (with the anthropisation process starting 1987 and intensifying in 2010), the negative effects on population density appeared 30 years after the occurrence of the disturbances. Our extinction debt measure (difference in the extinction probability at the beginning of the land-use change and the end of the simulation period) showed an additional 22% of extinctions at the end of the period (130 years after the land-use change occurred, about six T. graeca generations). Species like T. graeca that show longterm responses to habitat changes typically have ''slow life histories'' (Jeschke and Kokko 2009), and long-term legacy effects are expected to be commoner in long-lived species (Metzger et al. 2009). Ignoring time-lags in population and metapopulation dynamics will lead to an underestimated number of effectively endangered species (Hanski and Ovaskainen 2002).
Our results are important for the management of T. graeca populations because they show that marked anthropisation processes can indeed have detrimental effects on population viability, but with considerable time-lags (Fig. 3). It should be stressed that our simulations did not consider differences in mortality and clutching between natural and altered habitats (due to lack of empirical evidences of our T. graeca populations up to date), which would most likely lead to even stronger negative effects on population viability. Moreover, we also did not consider possible rises in mortality rates due to habitat destruction for 1987 and 2010 (e.g. direct or traffic mortality; Shine et al. 2004;Van Der Ree et al. 2011;Florencio et al. 2014;Marsh et al. 2017), or loss of individuals' fitness during this process (e.g. population stress). Nevertheless, these additive and/or synergic effects on the long term T. graeca demography could be relevant to be evaluated in future studies.

Legacy effects of agricultural abandonment
For agricultural abandonment, our models indicated only a slightly reduced reproductive performance due to changes in movement patterns, but it was not relevant as the confidence interval between the control/impact scenarios overlapped. This effect showed a negligible time-lag in population viability. It reveals how even minor impacts on demographic parameters, in this case reproductive rates, can show up decades after habitat modification. However, these changes in reproduction did not promote more marked long-term changes in density or extinction probability as reproductive rates were higher than the critical threshold for RR (Fig. 3a). Although agricultural abandonment in our simulations brought about an increase in the population densities of T. graeca, we observed a similar population increase in the control simulations. This was probably related to our study species' habitat preferences for cultural landscapes composed of habitat with scrubs and traditional agricultural areas (Rodríguez-Caro et al. 2017). Previous studies support this pattern by identifying a unimodal response (i.e. with an optimum at intermediate values) of T. graeca occurrence in relation to scrubland (Anadón et al. 2006b).
In order to evaluate the effect of agricultural abandonment on species and communities in general, we need to consider that the passive recovery of natural vegetation on abandoned agricultural landscapes is a complex process (Navarro and Pereira 2012), which often occurs towards the intermediate successional stages that diverge from the original biodiversity on a local scale (Baeten et al. 2010;Queiroz et al. 2014). Thus depending on the species' (or communities') biological characteristics, using real-world landscape scenarios of agricultural abandonment may not be sufficient, but would require a more detailed description of intrinsic abandonment processes, which can take distinct regeneration pathways depending on differences in adjacent land-use covers (De Aranzabal et al. 2008;Martínez-Fernández et al. 2015), with intermediate natural succession stages to primitive ecosystems. Additionally in these successions, natural perturbations can occur, such as invasive species or fires (Campos et al. 2020).
On the nature of time-lagged responses Our modelling approach also allowed us to do more than simply estimating the lagged impacts of a perturbation on species dynamics, but to gain insights into two barely explored aspects in demography. We used three population variables to describe three different demography levels that are sequentially time-related (i.e. movement impacts reproductive rates which, in turn, influence population density which, in turn, impacts extinction probability). This demographic framework allowed us to understand how a time-lag response propagated across population dynamics, especially in the study area suffering anthropisation process. That is, the time-lag cumulates by being short for reproduction, increases for density and is at its longest for population viability (Fig. 3). In a similar vein, our results also suggested propagation thresholds in the study area Los Mayorales. This was noted for the minor effects on reproduction, which did not result in changes of density or population viability (Fig. 3). To the best of our knowledge, how population time-lagged responses propagate across demographical parameters has virtually not been explored at all (Benton et al. 2006). Our study offers the first insights of time-lagged responses that are propagated across demographical parameters considering dynamic landscapes in general, (i.e. we have not explored differences between the two specific land-use change periods), but further research is necessary in this barely explored aspect of demography.

Conclusion
Our study contributes to biodiversity conservation by anticipating time-lagged demographic responses populations to land-use changes. Our results highlight the need to control the expansion of intensified land uses in fine-grained mosaic landscapes, which are widely increasingly worldwide. In this line, we conclude that measures to promote landscape connectivity, such as wildlife corridors and wildlife passages, are necessary to maintain natural habitats and, therefore, species' ecological functions (i.e. movement, reproduction and, ultimately, their population viability). We also conclude that although abandonment had in our example case no impact on the population dynamics of our study species, it is highly relevant to study agricultural abandonment in the context of habitat rewilding and natural perturbations like wildfires.
We exemplified how an interdisciplinary approach that integrates field studies, movement ecology and individual-based simulation modelling can be applied to solve a specific problem in landscape ecology (Teckentrup et al. 2019;Rocha et al. 2021), and to better understand and quantify potential impacts of land-use changes on long-term population dynamics (Rustigian et al. 2003;Alderman et al. 2005;Schumaker et al. 2014;Synes et al. 2016). Such an understanding is needed for providing suitable solutions in proactive landscape planning and population management. We expect that our approach can also be applied to different species inhabiting dynamic landscapes, especially for long-lived species with limited movement ability, and to be integrated into management plans (Stillman et al. 2015). This calls for a new research directions based on the development of more general tools that adapt to different species, such as new platforms for individual-based simulations to study species' responses to environmental changes Malchow et al. 2021).
Regional Development Fund, MINECO/FEDER, UE), as well as the Biodiversity Foundation and the OECC of the Ministry for Ecological Transition and the Demographic Challenge of the Government of Spain (Project CORREDOR, call 2020 for the Evaluation of Spanish Terrestrial Biodiversity). M.V.J.F. was supported by a postdoctoral grant cofunded by the Regional Valencian Government and the European Social Fund (APOSTD/2018/043), and a ''Juan de la Cierva-Incorporación'' research contract from the Spanish Ministry of Economy and Competitiveness (Reference IJC2019-039145-I). R.C.R.C. was supported by the Regional Valencian Government and the European Social Fund with a postdoctoral grant (APOST/2020/090). JDA is currently supported by a Ramón y Cajal contract (RYC-2017-22783) co-funded by the Spanish Ministry of Science, the Agencia Estatal de Investigación and the European Social Fund. Consent to participate Not applicable.

Consent for publication Not applicable.
Ethical approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.