Individual-based movement models reveals sex-biased effects of landscape fragmentation on animal movement

Animal movement and behavior may respond to habitat modification or fragmentation in non trivial ways, thereby strongly conditioning the fate of populations. This study aims to understand movement patterns of non-dispersal animals in both natural and altered landscapes, using the endangered terrestrial tortoise Testudo graeca as example. We used individual-based simulation models representing competing hypotheses on tortoise movement. Model parameterization and selection was based on radiotracking data and an inverse approach that is able to deal with observation uncertainty, individual variability, and process stochasticity. We find that land use intensification had a strong impact on the movement and behavior of non-dispersing individuals of T. graeca. In natural landscapes, males and females showed a similar movement and behavior profile with a strong home behavior component, and little individual variability. However, in altered landscapes, movement and behavior greatly varied among individuals, particularly in females, and males and females showed different movement patterns. Females showed a wide range of movement patterns, from strong home behavior to an unbounded movement. Our study shows that population or movement models that assume single behavioral states for animals inhabiting different landscape structures can be strongly misleading and, furthermore, that the impact of landscape modification on movement and behavioral patterns can be strongly sex-biased. Flexible, individual-based movement models coupled with inverse parameterization and model selection approaches proved useful in understanding the mechanisms controlling animal movement patterns.


INTRODUCTION
During the last decades there has been a shift in population ecology towards spatially-explicit approaches and animal movement has been recognized as a key ecological process (Tilman et al. 1997, Nathan et al. 2008).Introducing animal behavior into ecological studies of movement, as demanded by Lima and Zollner (1996), has improved the understanding of a wide array of issues such as spatial structure of populations and metapopulation dynamics (Revilla et al. 2004, Revilla andWiegand 2008), spread of exotic species and range expansion (Lensink 1997, Kramer-Schadt et al. 2004) or disease ecology (Kramer-Schadt et al. 2009).
From a conservation perspective, the key role of animal behavior and individual movement v www.esajournals.orgpatterns has been acknowledged when attempting to understand the effects of habitat alteration on population dynamics (Johnson et al. 1992, Darvill et al. 2010).The ability of animals to move allows exchange of individuals among habitat patches, (re)colonization of non occupied areas, avoidance of unsuitable areas, and maintenance of gene flow.Given the increasing concerns regarding habitat fragmentation and global change, all these aspects are key issues in conservation.There have been substantial advances in studying animal movement (e.g., Gautestad and Mysterud 2005, Getz and Saltz 2008, Nathan et al. 2008, Patterson et al. 2008, Schick et al. 2008, Smouse et al. 2010), and in the role of habitat heterogeneity on movement (Crist et al. 1992, Barraquand andBenhamou 2008), however, our understanding of animal movement in real human-altered landscapes is still limited although this is a key for conservation (Revilla et al. 2004).Observational studies developed in both natural and altered landscapes are needed to understand the spatial population processes as the ultimate consequence of the interaction of individual behavior with landscape structure and other constraints (Patterson et al. 2008).
Random walk models (Turchin 1998) are mechanistic approaches that have been frequently used to study animal movement and behavior (Kramer-Schadt et al. 2004, Bo ¨rger et al. 2008).The framework of random walk models allows estimation of movement metrics from movement data whereas simulations of mechanistic random walk models produce trajectories that can be compared to data (Patterson et al. 2008).This provides mean for identifying the underlying behavioral processes given the observed movement data (Moorcroft 2006, Bo ¨rger et al. 2008, Kie et al. 2010).Random walk models can also be expanded to include a habitat selection component (e.g., Morales et al. 2004, Revilla et al. 2004, Moorcroft and Barnett 2008, Forester et al. 2009).Crucially here, habitat selection processes are implemented at the same level as other behavioral processes (i.e., memory, home behavior) that eventually lead to the emergence of the home range (Cagnacci et al. 2010).Because different behavioral processes compete in their ability to explain observed movement patterns, the resulting habitat selection patterns are likely to be robust (Beyer et al. 2010, Cagnacci et al. 2010).This is a key issue when tackling the effects of habitat modification in the movement and behavior patterns of animals.
The general aim of our study is to understand movement patterns of adult (non-dispersing) individuals of the spur-thighed tortoise Testudo graeca in both undisturbed natural and altered landscapes and to infer the effects of habitat alteration or fragmentation on movement patterns.T. graeca is heavily threatened by habitat fragmentation (Anado ´n et al. 2006a) and it is ideally suited to pursue our aim because basic aspects of the use of habitat and space have been already addressed (e.g., Dı ´az-Paniagua et al. 1995, Pe ´rez et al. 2002, Anado ´n et al. 2005, 2006a, 2006b, 2007).Our data set comprises radiotracking data on 40 tortoises on five localities ranging from undisturbed to highly altered and fragmented landscapes.To reach our goal we developed a series of individual-based and spatially explicit movement models with increasing complexity to describe movements of T. graeca at an intraday temporal scale in a realistic way.These models are based on random walks and include our hypotheses on spatial autocorrelation in movement, habitat use, and homing behavior.These mechanisms are most commonly considered when characterizing individual animal movement (e.g., Bo ¨rger et al. 2008, Codling et al. 2008).The models are evaluated and fitted by comparing model predictions with observed movement data at a larger temporal scale (i.e., one year-cycle), both in natural and altered landscapes.To deal in model parameterization and selection simultaneously with observation uncertainty, individual variability, and process stochasticity we developed an inverse pattern-oriented approach to contrast the alternative movement models to our data (Wiegand et al. 2004, Hartig et al. 2011).

Study area and species
The spur-thighed tortoise T. graeca is a midsized terrestrial tortoise whose main population in Europe is located in Southeast Spain (Anado ´n et al. 2007).The species mainly inhabits semiarid shrublands and traditional agricultural landscapes and it is heavily threatened by habitat fragmentation and loss (IUCN 2002).In this sense, a good knowledge on the movement patterns of individuals in altered landscapes is the key for understanding the dynamics of their spatially structured populations (Revilla and Wiegand 2008).Current knowledge on the species shows that T. graeca individuals are non-territorial, with overlapping home ranges among individuals, which size ranges from 1 to 5 ha in natural areas (Anado ´n et al. 2006b).Radiotracking data together with long-term capture studies suggest that, as in other Testudinidae species (Freilich et al. 2000), individuals show strong site fidelity along many years (authors, unpublished data).

Radiotracking data and habitat mapping
Data for model development and fitting.-During2006 we radiotracked 29 tortoises on five different sites: Galera (1 males and 7 females), Sierrecica (3 males and 3 females), Palas (3 males and 3 females), Bas (3 males and 2 females) and Quiñoneros (3 males and 1 females).Whereas Galera comprises mostly natural habitats (.99%), the remaining localities represent landscapes heavily intensified due to new irrigated lands and high ways (intensive land uses from 24% to 44%; Appendix A: Table A1).All individuals were located once every two weeks during a one year cycle.
For characterization of intraday movement patterns we used radiotracking data from Anado ´n et al. (2006b).This database comprises 11 tortoises (6 females and 5 males) that were radiotracked in the Galera site (Appendix A: Fig. A1).Two weeks out every three weeks these individuals were located 5 times a day from sunset to sundown at equal periods of 2-3 hours (depending on the length of the day).All locations were registered with a GPS with an estimated error of location of 3 m.
Habitat maps.-Wecompiled habitat maps for the five sites considering four main habitats: areas heavily modified by humans (mainly new irrigated lands and main roads), dry crops (traditional agriculture lands) and two types of natural areas (Appendix A: Fig. A1).Because tortoises may show habitat preferences in relation to vegetation structure and relief (Anado ´n et al. 2005(Anado ´n et al. , 2006b)), we differentiated natural areas on relief from those on flat areas.This distinction describes main differences in vegetation structure in the study areas.Natural areas with no slope usually contain vegetation on early successional stages with very simplified structures originating from the abandonment of traditional dry crops, whereas areas on slope contain more complex vegetation structures with higher cover values (Anado ´n et al. 2005).

General modeling approach
In order to assess the potential differences in the movement and behavior of female and males tortoises in natural and fragmented landscapes we developed a series of nested movement models based on random walks.These models include our hypotheses on spatial autocorrelation in movement, habitat use, and home behavior.Each model run is a stochastic simulation of one trajectory, comprising individual movement steps during one year at an intra-day resolution and on a grid of 10 3 10 m cells (e.g., Revilla et al. 2004, Kramer-Schadt et al. 2004).
For parameterization of the competing models we employed a rejection filter approach similar to that used in Approximate Bayesian Computation (ABC) with a systematic parameter search through the parameter space (Wiegand et al. 2004, Csille ´ry et al. 2010, Hartig et al. 2011).Model predictions of different parameterizations (and models) were compared with observed movement data by comparing summary statistics derived from the raw movement data (observed and simulated) (Wiegand et al. 2003, 2004, Grimm et al. 2005).We accepted only model parameterizations that yield simultaneous match in several summary statistics.To this end we developed statistical acceptance criteria that considered measurement error, process stochasticity and individual variability.A detailed description and rationale of the modeling approach can be found in Appendix B.

Movement models
We hypothesized a series of eight alternative movement models that resulted from all combinations of four basic blocks: (1) a pure random walk (rw), (2) a correlated random walk (Crw) where the direction of consecutive movements steps was correlated, (3) a biased random walk (Brw) to describe home ranging behavior where movement is consistently biased by attraction to a focal point, (4) an habitat dependent random walk (Hrw) in which the animal enters different habitat types with different probabilities.The eight resulting candidate models were: rw, Crw, Brw, Hrw, CBrw, CHrw, BHrw and CBHrw.
Pure random walk (rw).-Thismodel simulates movement steps from a given cell to one of the eight neighboring cells in the simplest possible way: each of the eight neighboring cells had the same probability to be selected.To adapt our model to the radiotracking data structure and to give movement an explicit temporal dimension (Revilla et al. 2004) we used the (2-3 hour) interval between two consecutive radiotracking localizations as ''movement step'' and then simulated the cell-to-cell movement (i.e., ''cellsteps'') during one movement step.This structure required one parameter PMOV, the probability that the individual changed its location during one movement step, and the distribution DMOV(s), the probability to move s cell steps during one movement step.The PMOV and DMOV(s) were estimated, for each month and for males and females separately, directly from the detailed radiotracking data at the Galera site (Anado ´n et al. 2006b).All eight movement models were based on the algorithm described above, but in more complex models weights were used to describe how the probability to move to one of the eight adjacent cells (i.e., one cell step) may depend on autocorrelation in movement, habitat use, and home behavior.
Correlated random walk (Crw).-Tomodel a tendency of keeping the previous direction we introduced two parameters AU1 and AU2.AU1 describes autocorrelation between consecutive movement steps.It only affects the first cell step of each movement step, whereas the autocorre-lation in movement between consecutive cell steps (i.e., inside one movement step) is described by the parameter AU2.
Homing behavior (Brw).-Homingbehavior, i.e., the tendency of an individual to maintain a stable home range, was implemented as modified focalpoint attraction model (Bo ¨rger et al. 2008) governed by two parameters (dHB and rHB).The first parameter dHB controls the distancedependent strength of the homing behavior relative to the focal point attractor.The second parameter rHB describes a temporal delay in the start of the home behavior once a minimal distance from the attractor (i.e., dHB) has been reached.
Habitat dependent random walk (Hrw).-Inrelation to habitat quality, the eight neighboring cells received habitat-dependent weights h i ¼ Q i / P Q i , where Q i is the habitat quality of the cell i.For each of the four habitat types (i.e., intensive land use, traditional agriculture land, flat natural areas, and natural areas on slope) we introduced one parameter (H1W, H2W, H3W, and H4W, respectively) that gives the associated habitat quality.

Model parameterization and evaluation
In total we have 10 parameters that govern movement of male and female T. graeca in the different landscapes (Table 1).The values of the eight parameters that could not be directly parameterized from our data (all except PMOV and DMOV(s)) were estimated by means of an inverse approach based on rejection filters that identify from a large systematic sample of the parameter space those parameterizations that yield model outputs that agree simultaneously in several features (quantified by summary v www.esajournals.orgstatistics) with the observed data (for details see Appendix B).This required three steps: (1) the information contained in the raw data was condensed into summary statistics to capture different characteristic features of movement that are relevant for the scientific question asked (Wiegand et al. 2004, Grimm et al. 2005, Hartig et al. 2011, Martı ´nez et al. 2011), (2) a given parameterization of a model was accepted or rejected based on statistical criteria that evaluated the deviation between observed and simulated summary statistics considering observation uncertainty, individual variability, and process stochasticity (Wiegand et al. 2004), and (3) the candidate model(s) most likely given the data were selected (Wiegand et al. 2003).
Summary statistics.-Ourdata comprise the trajectory of 40 radiotracked individuals.We summarized the basic properties of each trajectory by 8 summary statistics: home range size (HRs), home range perimeter (HRp), distance to the starting point after one year cycle (Dy), distance of the furthest location to the starting point (Dmax), and proportion of locations in habitat type 1, 2, 3 and 4 (H1, H2, H3, H4, respectively).These descriptors of the trajectories cover the spatial pattern of the movement, home behavior, and the impact of habitat on movement; thus describing the characteristics of our data with respect to our objectives and hypothesis in a comprehensive way.
Assessing deviation of simulated data from observed data.-Inorder to compare model predictions with observations, we created for each parameterization and model 500 simulated data sets that followed the field sampling protocol of the observed data.The movement of each individual was simulated in the corresponding habitat map (site) starting from the coordinate where the individual was equipped with a transmitter.The model simulates the movement at an intraday resolution and, for every individual, the 500 simulated annual trajectories were sampled with the same protocol as used for field data (i.e., the same day and hour).
To take the different sources of uncertainty and stochasticity into account, statistical comparison between observed and simulated data was performed in a three-step procedure (Appendix B: Fig. B1) that accounts for (1) the stochasticity present in one simulation (one movement trajec-tory) due to the observation process (i.e., sampling protocol), (2) the stochasticity among simulations (i.e., process stochasticity in the model), and (3) the observed variability among individuals (see Appendix B for further details).The outcomes of the first two steps are two measures (m ij and v ij ) of the deviance between the summary statistics j of the observed and simulated trajectory of individual i.The third step determines thresholds tm ij and tv ij for acceptance of summary statistic j of individual trajectory i (i.e., m ij , tm ij and v ij , tv ij ).Basically we accept a parameterization if the deviance between simulated and observed trajectories is smaller than the deviance between observed trajectories.The result of this procedure is a binary matrix that tells us for a given model, parameterization, summary statistic j, and individual i whether the simulations successfully matched the observations or not, following a priori established statistical criteria (see Appendix B).
Model selection and analysis.-Toaddress our question regarding potential differences in movement at natural and disturbed landscapes we analyzed the alternative models separately for males and females and in natural and altered landscapes.A movement model is accepted (for each of these four cases) if we found at least one parameterization which successfully matched 95% of all summary statistics (i.e., n individuals 3 8 summary statistics 3 2 indices m ij and v ij ).The parsimony principle may be used to reject more complex models if several models are accepted.Studying in detail in which movement property (i.e., summary statistic) a given model failed provides important information on the processes governing tortoise movement.To assess where the different models failed we selected for each model the 5% best parameterizations, based on the number of successfully predicted patterns, and we calculated for each pattern the percentage of accepted parameterizations.
To quantify the similarity between the movement properties of males and females and in altered and natural landscapes, we calculated the overlap between the successful parameterizations for the model between males and females and in altered and natural landscapes (this is a type of cross validation).Finally, the inspection of the parameter values that successfully predicts the summary statistics allowed us to characterize and quantify movement properties of males and males in both natural and fragment landscapes in relation to single parameters and their interaction.

Empirical patterns
The probability of movement during one movement step ( p) assessed from field data ranged between 0 (no movement due to hibernation in December, January and February and due to aestivation in July) to 0.38/0.33 for males and females in April, where the peak of the spring activity period was reached (Fig. 1).
Considering the entire annual cycle, half of the movement steps were between 5 and 15 m (i.e., to the neighbored 10 3 10 m cell) for both males and females (Fig. 1).For females the maximum step length was ,100 m, whereas males made movement steps up to 330 m.Males also showed a peak of movements at mid-length values (35-45 m).

Observed patterns and individual variability
Inter-individual comparisons of the summary statistics of the observed data showed that female trajectories differ substantially between natural and altered landscapes.In altered areas, the interindividual variability in the summary statistics was much larger than that in natural areas (see Appendix C and Table C1 for further details).Thus, we needed to estimate the thresholds tm j and tv j (used for parameterization acceptance; see Methods and Appendices B and C) separately for individuals in natural and altered landscapes and for males and females (Appendix C: Table C2).Because individuals in natural landscapes behaved in a more homogeneous way than individual in altered landscapes, the thresholds for natural landscapes were more restrictive that for altered landscapes.Mean values of the summary statistics for the radiotracked males and females in natural and altered landscapes are given in Appendix C: Table C3.

Comparison among movement models
Natural landscapes.-Onlycandidate models that contained the home behavior component were able to match the summary statistics of the observed movement of females and males in natural landscapes: Brw (only females), CBrw, BHrw and the full model (CBHrw) (Table 2).As shown in Appendix D: Table D1, the four summary statistics of space use (HRa, HRp, Dmax and Dy) were for females in natural areas practically impossible to match if the home behavior component was not included in the model.For males in natural landscapes, these summary statistics (except HRp) were also difficult to match without home behavior, but to a lesser extent than in females.In this group, displacement after one year (summary statistic Dy) was difficult to match even in successful  D1).Interestingly, the habitat related summary statistics (i.e., H1, H2, H3, H4) were in natural landscapes in most cases successfully matched.For a given model, the percentage of successful parameterizations was slightly higher among females (2.5-5%) than males (,2%) (Table 2).
Altered landscapes.-Themovement of females was matched by all models containing a response to habitat (i.e., Hrw, CHrw, BHrw and the full model; Table 2) whereas the movement of males was matched only by the two models that included home behavior and a response to habitat (i.e., BHrw and the full model CBHrw; Table 2).The most restrictive summary statistic in altered landscapes was the proportion H1 of locations in intensive land uses (Appendix D: Table D1).We found that matching summary statistic H1 required a parameter value H1W ¼ 0 (i.e., no use of H1 at all).Thus, only models that contained the habitat component were able to match this summary statistic.For males, but not females, summary statistics HRa, Dmax and Dy were also restrictive and models including the home behavior component yielded in general a better match (Appendix D: Table D1).The number of successful parameterizations for females (12-14%) was much larger than for males (,1%; Table 2).
Differences between movement in altered and natural landscapes.-Weevaluated the ability of each model parameterization to match the summary statistics separately for females and males and for natural and altered landscapes.The interesting question that arises is if those four cases shared any accepted parameterization.Our analysis in altered landscapes showed that tortoises completely avoided intensive land uses (i.e., habitat type H1).Because only very few cells with intensive land uses are present in natural landscapes (,1%; Appendix A: Table A1), we obtained in natural landscapes also successful parameterizations with H1W 6 ¼ 0. We therefore used for the analysis of parameter overlap only parameterizations with H1W ¼ 0.
For the full model CBHrw, overlap analysis of the successful parameterizations (Table 3) showed that there was not a single parameterization shared by the four cases.Males and females in natural areas showed the highest similarity in model parameterizations (19% overlap) whereas movement properties of males and females in altered areas showed low overlap (1%).There was also a strong dissimilarity between the parameterizations found for natural vs. altered areas (0.07% overlap).Similarity in the movement properties of males between natural and altered landscapes was medium (11% overlap); that is, more similar than females in Notes: On every case (sex and landscape type) the total number of summary statistics is ¼ 8 3 number of individuals; np: number of summary statistics matched by the best parameterization; worst: number of summary statistics matched by the worst parameterization; n: number of accepted parameterizations (i.e., match in more than 95% of all indices m ij and v ij ); easy: number of summary statistics matched by all parameterizations; imp: number of summary statistics not matched by any parameterization.For a description of the models see Methods: Movement models.
v www.esajournals.orgnatural/altered landscapes, but more dissimilar than males and females in natural landscapes.Overlap analysis for the BHrw model, the second model that was able to produce match in all four cases, yielded similar overall results (not shown).

Parameter estimation
Natural landscapes.-Parameterscontrolling for home behavior (dHB and rHB) were the most relevant when simulating the movement of female and male tortoises (Fig. 2).For females, match of the summary statistics was achieved when the parameter dHB (that gives the distance at which attraction from the focal point starts) was smaller than 100 m, no matter which other movement processes acted (i.e., autocorrelation in movement or habitat) (Fig. 2).The delay in home behavior was in general smaller than 50 days.The same results were found for males, but here the parameter space was somewhat larger.Interestingly, the shape of the accepted dHB-rHB parameter space was basically the same for the four models Brw, CBrw, BHrw, and CBHrw (Fig. 2).It is more sparely filled for CBrw, BHrw, and CBHrw, probably because of the lower number of accepted parameterizations for these models.This result shows that tortoise show in natural landscapes a regular home range with excursion behavior.
Accepted parameterizations of CBrw models included only high values of intra-movement autocorrelation, indicating straight movements at a small scale (i.e., between consecutive 10 3 10 m cells).This pattern was clear in females but not in males due to the low number of successful parameterizations.However, this pattern disappeared in the full model, where all the range of parameters values achieved successful movement simulations (Appendix D: Fig. D1), thus indicating that the inclusion of autocorrelation in movement direction do not significantly alter the general behavior of the model.However, the CBrw model and the full model showed for both males and females a weak interaction between autocorrelation and home behavior: the higher the values of rHB and dHB the less autocorrelation is required in movement to match the observations (Appendix D: Fig. D2).
Models for males included mostly medium weights of natural habitat on flat areas (habitat 4) (Fig. 3 and Appendix D: Fig. D3).For females, the parameter space of accepted parameterizations of the full model (CBHrw) was not restricted by habitat weight parameters.
Altered landscapes.-Themost restrictive parameter for both males and females was the weight of intensive land uses.Only those parameterizations where HW1 was equal to zero (no use at all) matched the observed movement of individuals.For females the parameter space related with the other habitat types was not restricted, whereas successful parameterizations for males included more frequently very low weights for traditional land uses (habitat 2) (Fig. 3 and Appendix D: Fig. D3).
Females and males showed very different movement patterns in relation to home behavior (Fig. 2).In the BHrw model (i.e., excluding autocorrelation), only parameterizations with dHB , 125 m and rHB .20 days were selected, but it is not clear if this is simply a result of the low number of accepted parameterizations.However, accepted parameterizations for females, in both BHrw and the full model, only excluded parameterization with very low values of dHB and rHB (Fig. 2).For males and females, the area in the AU1-AU2 parameter space of the full model with accepted parameterizations was very similar to that in natural areas (Appendix D: Fig. D1).

DISCUSSION
A good knowledge of how the movement and behavior of individuals respond to different landscape structures is key for understanding Notes: Analysis based on the full model CBHrw.Absolute values in italics indicates number of successful parameterizations for males (M), females (F) and in natural (Nat) and altered landscapes (Alt), and their intersections (\).The % of the intersection is calculated as 2 3 A\B/(A þ B).
For parameterization in natural habitats only those parameterizations where H1W ¼ 0 were included (see Results for an explanation).
v www.esajournals.orgthe effects of habitat alteration and fragmentation on the species (Wang andGrimm 2007, Revilla andWiegand 2008).To the best of our knowledge, this study constitutes one of the first attempts to study the effects of habitat modification on the movement and behavioral patterns of animals by means of individual-based models.
Our study provides a clear picture on the movement patterns and behavior of non-dispersing individuals of T. graeca.Our main result is that land use intensification has a strong impact on the movement and behavior.While males and females showed in natural landscapes a similar movement and behavioral profile (characterized Movement in natural vs. altered landscapes Our results show that movement in natural landscapes can be described by random walk models that need to include a strong home behavior component.For females we accepted as simplest model the (homing) biased random walk (Brw), whereas males required slightly more complex models with autocorrelation or habitat bias (CBrw or BHrw).Both males and females presented a similar movement and behavior profile, although with slight variations in the parameters values.Females showed home behavior at somewhat shorter distances than males, and males showed habitat selection patterns, with a lesser preference for natural areas on flat relief than dry crops and natural areas on relief.The movement and behavior of individuals of T. graeca changed notably in altered landscapes with less natural landuses, larger patch size and the existence of intensive land uses (Appendix A: Fig. A1 and Table A1) These changes were threefold, related to (1) inter-individual variability, (2) dissimilarity between sexes and (3) movement properties.In relation to inter-individual variability, individuals moved and behaved in natural areas in a much more similar way than in altered landscapes, particularly females, that showed in altered landscapes great differences v www.esajournals.orgamong individuals.Regarding sex-related differences, in contrast to natural landscapes, movement patterns in altered landscapes were strikingly different among sexes.In relation to movement properties, both males and females completely avoided intensified patches (mainly new irrigated lands).In the case of males, they also avoided, but not completely, traditional dry crops.Interestingly, males maintained in altered landscapes largely their behavioral profile of natural landscapes, whereas females changed it.Movement and behavior of females in altered landscapes were in general only ruled by the avoidance of intensified patches, and the home behavior component varied among individuals, from strong home behavior, similar to individuals in natural areas, to an unbounded movement with no apparent home behavior component.
The differences in movement between males and females and altered vs. natural landscapes have a clear biological interpretation.The space used by non-dispersal individuals (i.e., the home range) has been classically defined as the spatial expression of behaviors that animals perform to survive and reproduce (Burt 1943), by providing the basic resources.In natural landscapes, males and females of T. graeca seem to meet their requirements in relatively small areas (1-3 ha; Appendix C: Table C3; Anado ´n et al. 2006b), probably because the environmental heterogeneity of these natural landscapes is organized at a fine grain due to relief and their history of traditional land uses (Anado ´n et al. 2005(Anado ´n et al. , 2006b)).The use of small areas yields to a strong home behavior component.The apparently unbounded movement of females, but not males, in altered landscapes suggests that in those landscapes basic conditions and resources needed by females are very sparsely distributed and require females to expand their home range considerably.The most evident resources needed exclusively by females are places for egg-laying.However, a visual inspection of the trajectories of females suggests that unbounded movements occur all along the annual cycle and not only during the egg-laying period (from April to June; Pe ´rez et al. 2002).Females are likely to have more specific nutrient needs than males for egg development (Henen 1997) and these resources could also be more patchy distributed in altered landscapes, that are organized at a coarser grain (larger patch size).This more extensive search for basic resources and conditions by females could explain that females use all habitat types (except intensive land uses) whereas movement of males was mainly confined to natural habitats.
On the other hand, as pointed out above, the role of memory of individuals has recently been advocated as key for understanding animal movements, particularly in non-dispersal individuals, since memory is a biological process that may lead to the appearance of home behavior (i.e., individuals returning to places already visited, Bo ¨rger et al. 2008, Van Moorter et al. 2009).In our work females showed much greater individual variability than males, a fact that could also suggest that memory may work differently for females and males.This hypothesis is plausible; for example, females fitness depends on spatially predictable resources (i.e., fixed in the space; e.g., places for egg lying, high quality feeding patches) and their movement could rely on spatially fine-grained cognitive spatial maps, thus being largely affected by alterations in the landscape.On the contrary, fitness of males depends on encountering the maximum number of females (Sutherland 1996), a more spatially unpredictable landscape attribute, for which memory (particularly finegrained) can be less useful and where other orientation mechanisms, not affected by land use changes, certainly participate (i.e., chemical cues; Poschadel et al. 2006).

Implications for population dynamics and conservation
Our results show that land use changes lead to strong alterations in the movement and behavioral patterns of non-dispersal tortoises.The effects of alteration of movement patterns on the population dynamics of non-dispersal animals are largely unknown and remain to be explored (but see Armsworth andRoughgarden 2005, Wang andGrimm 2007).Our study that identified different random walk models for different landscape structures and individual attributes (sex), reinforces a promising research line which allow detecting switches between different behavioral states (e.g., Morales et al. 2004, Delgado andPenteriani 2008).Our analysis indicates that movement models that employ a single behavioral state for animals inhabiting v www.esajournals.orgdifferent environmental situations (in our case landscape structures with varying degree of land use intensification) may be strongly misleading.This means also that land use changes may have a deeper impact on the species dynamics than expected because the effects of land use changes were strongly sex-biased in our case study, with stronger impact on females.Females of T. graeca travel much longer distances in altered landscapes than in natural areas and use traditional land uses (where human activities are frequent) more than males, and thus are likely to suffer higher mortality rates than males (i.e., Pe ´rez et al. 2004).Furthermore, energy employed for traveling large distances is likely to be at the expense of energy budgets otherwise dedicated to other lifehistory traits, such as reproduction or growth.Interestingly, the impact of land use changes on movement patterns could also have counterintuitive effects on population dynamics.In general it is assumed that land use intensification may reduce functional connectivity and constrains genetic flow compared to that in undisturbed areas (Young and Clarke 2000).In our case study the observed increase in movement distance in disturbed areas may increase functional connectivity and genetic flow between suitable patches at a landscape scale (provided the extant landscape structure allows for it) counteracting to some extent the negative effects of fragmentation.Further efforts should be dedicated to the development of spatially-explicit population models that can integrate demographics and dynamics with the differential movement and behavioral patterns among different landscapes structures, an issue yet barely explored in ecology and conservation (Johnson et al. 1992, Wang and Grimm 2007, Hawkes 2009).

Movement modeling issues
Our work proves that individual-based movement models coupled with inverse parameterization and model selection approaches can be useful in understanding the mechanisms controlling animal movement patterns.In particular, our results stress the importance of using several summary statistics when comparing simulated and real data, as pattern-oriented modeling advocates (Wiegand et al. 2003, 2004, Grimm et al. 2005, Hartig et al. 2011).This point has also been made recently by Dalziel et al. (2010) who outlined that summary statistics have the advantage of focusing on models' ability to predict emergent patterns, rather than focusing on fit and complexity.The summary statistics should correspond to the specific question asked and describe different emergent property of observed movement.For a detailed discussion on this issue see Wiegand et al. (2004) and Grimm et al. (2005).However, if important properties of the data are not used as summary statistics, the inverse approach may fail because the summary statistics do not contain enough information to constrain the value of model parameters (Wiegand et al. 2004).Using summary statistics instead of the (raw) data (see also Wood 2010 andMartı ´nez et al. 2011) may simplify the model fitting procedure but comes with the cost that information of the raw data may be lost (i.e., if the summary statistics are not sufficient; Csille ´ry et al. 2010, Hartig et al. 2011).Alternative approaches such as state-space models try to fit the position x of the animal at time t directly to the data to make inference on the movement properties (Patterson et al. 2008).Thus, the approach of using summary statistics shifts the problem from fitting high-dimension models to finding the ''right'' summary statistics.
Our approach of using rejection filters has advantages over approaches based on synthetic likelihood functions (e.g., Wood 2010, Martı ´nez et al. 2011).It is basically a statistical test to find out if we can distinguish the observed value of a summary statistic from the simulated value, given the uncertainty and stochasticity in the observed data.This avoids overfiting and more importantly, allows us to reject or accept a given parameterization.Model predictions can then be based on all accepted parameterizations of all models (because they were consistent with the data).The rejection filter approach offers the possibility to directly detect structural errors in the model that a single and relative goodness-offit measure glosses over.This allowed us to better understand why certain models failed (Appendix D: Table D1).However, the systematic parameter search required for the rejection filter approach may reach computational limits if the number of parameters becomes larger.In this case optimization methods are required which can be based on synthetic likelihood functions (e.g., Martı ´nez et al. 2011).

Model assumptions
We deliberately used very simple models to describe and analyze the tortoise movement.Our underlying strategy was to assess how well we can describe the observations with simple representations of reality.It was not clear a priori if our data were sufficient to sustain more complex and realistic models, given the large degree of stochasticity and relatively low sampling frequency in the observed trajectories.Our results showed that several of the candidate models were able to generate trajectory with properties that we can not distinguish statistically from the observed trajectories.It is thus unlikely that a substantial increase of complexity could be justified to fit the present data set.However, our approach tells us precisely where the bottlenecks are (Appendix D: Table D1) and his information can guide model improvement.
The models used here rely on several assumptions that could be improved in a number of ways.For example, the number of movements per day, assessed from PMOV and the length of the movement (DMOV) have been empirically obtained from individuals in natural areas (per month and sex) and applied to all individuals.We cannot exclude the possibility that some characteristics of these parameters, particularly for example the distribution of movement lengths may vary in altered landscapes (Moorcroft 2006, Moorcroft andBarnett 2008).Furthermore, different months and sex combinations may present different variability in their PMOV and DMOV values (e.g., Morales et al. 2004) which has not been taken into account in this work.In relation to the habitat, spatial resolution of the model (10 3 10 m) has been successful for describing main differences between natural vs. altered landscapes but may be too coarse for detecting finer habitat selection patterns, likely related to microhabitat conditions.
Our model relies on a home behavior component based on a modified focal point attraction approach that imposes a stable home range with a home range center defined a priori.However landscapes where home behavior was strong (i.e., natural landscapes) are quite homogeneous (Appendix A: Fig. A1), and therefore we expect that the exact location of the focal point will not condition year-movement patterns.On the other hand, in altered landscapes, where the exact location of the focal point can be important due to the habitat heterogeneity, we found that home behavior was much less important and will also not condition the year-movement patterns.Although methodologically useful, the phenomenological point attraction approach does not directly describe the true biological processes behind the emergence of home ranges (e.g., memory), that could shed light into a deeper biological understanding of animal movement (Bo ¨rger et al. 2008, Van Moorter et al. 2009).Up to now, memory effects have been successfully included in random walk models to describe animal movement in theoretical landscapes and more rarely on natural landscapes (Bo ¨rger et al. 2008) and considering only one basic resource.Introducing complex memory effects on real individuals and landscapes as a substitute of an imposed home behavior remains a challenge (but see Dalziel et al. 2008, Fryxell et al. 2008)., 18, 0, 3, 16, 11) and (1,1,16,16,17,3,3,3). v www.esajournals.org

APPENDIX B DETAILED DESCRIPTION AND RATIONALE OF THE MODELING APPROACH GENERAL MODELING APPROACH
In order to assess the potential differences in the movement and behavior of female and males tortoises in natural and fragmented landscapes we developed a series of nested movement models based on random walks.These models include our hypotheses on spatial autocorrelation in movement, habitat use, and home behavior.For parameterization of the competing models we employed a rejection filter approach with a systematic parameter search through the parameter space (Wiegand et al. 2004, Hartig et al. 2011) similar to rejection algorithms used in Approximate Bayesian Computation (Csille ´ry et al. 2010).Model predictions of different parameterizations (and models) were compared with observed movement data by comparing summary statistics derived from the raw movement data (observed and simulated), both in natural and altered landscapes (Wiegand et al. 2003, 2004, Grimm et al. 2005).A special technical challenge in quantitative comparison of models and data is the simultaneous inclusion of different sources of uncertainty and stochasticity such as measurement error and process stochasticity (Patterson et al. 2008).In our study we could cope with individual variability, measurement error and process stochasticity simultaneously, but on the cost of rather complex acceptance criteria of our rejection filter.
Each model run is a stochastic simulation of a trajectory, comprising individual movement steps during one year at an intra-day resolution and on a grid of 10 3 10 m cells (e.g., Kramer-Schadt et al. 2004, Revilla et al. 2004).This spatial resolution is fine enough to represent the relevant habitat types and biologically meaningful since previous works on T. graeca reports average daily movement distance of 50 m with a maximum daily distance of 1019 m (Dı ´az- Paniagua et al. 1995).The extent of the study landscapes was 270 3 270 cells (i.e., 729 ha).To compare models predictions with observations, we created for each parameterization and model 500 simulated data sets that followed the field sampling protocol of the observed data (Kramer-Schadt et al. 2004, Revilla et al. 2004).

MOVEMENT MODELS AND PARAMETERS
We hypothesized a series of eight alternative movement models that resulted from all combinations of four basic blocks: (1) a pure random walk (rw), (2) a correlated random walk (Crw) where the direction of consecutive movements steps was correlated, (3) a biased random walk (Brw) to describe home ranging behavior where movement is consistently biased by attraction to a focal point, (4) an habitat dependent random walk (Hrw) in which the animal enters different habitat types with different probabilities.The eight resulting candidate models were: rw, Crw, Brw, Hrw, CBrw, CHrw, BHrw and CBHrw.
Pure random walk (rw).-Thismodel simulates movement steps from a given cell to one of the eight neighboring cells in the simplest possible way: each of the eight neighboring cells had the same probability to be selected.To determine the total number of movement steps for a given individual and to give the movement an explicit temporal dimension (Revilla et al. 2004), we used a two step procedure.We first determined the probability (PMOV) that the individual changed its location within 2-3 hours (i.e., between 2 consecutive radiotracking localizations).This movement is called in the following ''movement step''.Because individuals where located 5 times per day at 2-3 h intervals, individuals may conduct between 0 and 4 movement steps per day.To simulate the detailed movement between two consecutive radiotracking localizations we used the empirical distribution DMOV(s) which gives the probability for a certain number of ''cell steps'' s (i.e., the number of cell-to-cell steps) the individual performed within a 2-3 hour period.The possible values range from 1 step to the maximal observed distance.All eight models were based on the algorithm described above, but in more complex models weights were used to describe how the probability to move to one of the eight adjacent cells may depend on autocorrelation in movement, habitat use, and home behavior.
The parameter PMOV and the distribution DMOV(s) were estimated directly from the detailed radiotracking data at the Galera site which represented natural areas (Anado ´n et al. 2006b).Because the movement pattern of tortoise shows a strong seasonality and may differ between males and females we estimated PMOV and DMOV(s) for each month and for males and females separately.PMOV(s) is given by the radiotracking protocol where individuals where located 5 times per day.We assumed that the tortoise have the same probability of moving ( p) in any of the four periods.Thus, the probability PMOV(s) of s movement steps in one day was ð 4 s Þp s ð1 À pÞ 4Às with p being a parameter to be determined from the data.Because of the spatial resolution of the model (10 m, one cell) the minimum distance considered was 5 m and the maximum distance matched the maximum distance registered in the field.To keep the model as simple as possible and because of lack of detailed movement data in altered landscapes we assumed in first approximation that PMOV and DMOV(s) were the same for natural areas and altered areas.Note that our model analysis procedure tests this assumption by comparing observed and simulated larger-scale properties of movement (i.e., one localization every 2 weeks).If this model assumption would not hold then simulations conducted in altered landscapes should have substantial difficulties to match the data.

Correlated random walk (Crw
).-To model a tendency of keeping the previous direction j, weights d i (i ¼ 1, . . ., 8) were calculated giving the probability to enter one of the eight neighboring cells.The index i ¼ 1 indicated the previous direction, whereas cells with index i ¼ 2, . . ., 8 were clockwise numbered starting from i ¼ 1.The weights were calculated for i ¼ 1, . . ., 5 as d i ¼ max [0, (1.5 þ (0.5 À i ) 3 AU)/(1.5À 0.5 3 AU)], and for the remaining weights we assumed isotropic movement, i.e., Finally the weights were normalized to sum up one.The parameter AU determines the degree of autocorrelation.For AU ¼ 0 all eight directions have the same weight (i.e., a pure random walk), and with increasing value of AU the weight of the previous direction (i.e., d 1 ) increases, and for AU ¼ 1, d 1 ¼ 1 (always following the same direction).
Because of the nested structure of two movement scales (cell step and movement step), we used two different parameters to describe autocorrelation in movement.We used the parameter AU1 to describe autocorrelation between consecutive movement steps.It only affects the first cell step of each movement step whereas the autocorrelation in movement between consecutive cell steps (i.e., inside one movement step) is described by the parameter AU2.
Home behavior (Brw).-Homebehavior, i.e., the tendency of an individual to maintain a stable home range, was implemented by means of a second set of weights b i (with i ¼ 1, . . ., 8 normalized to sum up one) following a modified focal-point attraction model (Bo ¨rger et al. 2008) and using two parameters (dHB and rHB).Here we arbitrary considered the first location of the individual as the focal-point attractor.As we argue in the discussion, due to the specific tortoise movement characteristics the exact position of the focal-point attractor is not likely to affect the outcome of our analyses.The first parameter dHB governs the home behavior: if the distance D from the current location of the individual to the focal-point attractor is smaller than parameter dHB all eight directions have the same weight (no tendency to move back), but for dHB , D , 2 3 dHB the bias toward the focal point increases proportionally with distance D up to distance 2 3 dHB where the bias is maximum (e.g., a vector of weights with b j ¼ 1 where j corresponds to the cell closest to the focal point and b i ¼ 0 for i 6 ¼ j ).
To increase the biological realism of this model we added a second parameter rHB describing a temporal delay in the start of the home behavior once the distance dHB had been reached.This parameter accounts for realistic biological situations where tortoises temporally leave their most regularly used areas to meet specific needs such as egg-laying in females or mate-searching in males (Dı ´az-Paniagua et al. 1995, Pe ´rez et al. 2002) Habitat dependent random walk (Hrw).-Inrelation to habitat quality, the eight neighboring cells received habitat-dependent weights h i ¼ Q i / P Q i , where Q i is the habitat quality of the cell i.For each of the four habitat types (i.e., intensive land use, traditional agriculture land, flat natural areas, and natural areas on slope) we introduced one parameter (H1W, H2W, H3W, and H4W, respectively) that gives the associated habitat quality.

MODEL PARAMETERIZATION AND EVALUATION
In total we have 10 parameters that govern movement of male and female T. graeca in the different landscapes (Table 1), two were related with the pure random walk, two with autocorrelation, two with home behavior and four with habitat quality.The values of the eight parameters that could not be directly parameterized from our data were estimated by means of an inverse pattern-oriented modeling approach (POM; Wiegand et al. 2004, Grimm et al. 2005, Hartig et al. 2011).POM is an inverse modeling technique (i.e., estimating model parameters indirectly from matching larger-scale observational data) based on rejection filters that identifies from a large systematic sample of the parameter space (see below, Determination of parameter ranges for the different models) those parameterizations that yield model outputs that agree simultaneously in several features (quantified by summary statistics) with the observed data.
Our procedure of model parameterization and model selection comprised three steps: (1) the information contained in the raw data was condensed into summary statistics to capture different characteristic features of movement that are relevant for the scientific question asked (Wiegand et al. 2004, Grimm et al. 2005, Hartig et 2011, Martı ´nez et al. 2011), (2) for a given parameterization and model the deviation between the values of the summary statistics of the simulated and observed data was quantified, considering stochasticity and uncertainty in the data and the model.In this step the decision was made if the given parameterization was rejected or considered to match the observations (Wiegand et al. 2004), and (3) the candidate model(s) most likely given the data were selected (Wiegand et al. 2003).

Summary statistics
Judging model fit is difficult in situations where the observed data have a strong stochastic component which prevents the real system to reproduce the exact course of the observed data if repeated (Wood 2010, Hartig et al. 2011; Fig. A2).Under such circumstances, which include radio-tracking data of tortoise, model fit can be judged using statistics that reflect what is dynamically important in the data, and to discard noise (Wood 2010).The observed data was therefore summarized into several summary statistics that capture the characteristic features of tortoise movement of interest.These summary statistics are called in POM ''patterns'' (e.g., Grimm et al. 2005, Hartig et al. 2011) and sometimes ''probes'' (e.g., Kendall et al. 1999, Dalziel et al. 2008, 2010, Wood 2010).
Our data comprise the trajectory of 40 radiotracked individuals.We summarized the basic properties of each trajectory by 8 summary statistics indexed with j ¼ 1, . . ., 8 (Table 2): home range size (HRs, j ¼ 1), home range perimeter (HRp, j ¼ 2), distance to the starting point after one year cycle (Dy, j ¼ 4), distance of the furthest location to the starting point (Dmax, j ¼ 3), and proportion of locations in habitat type 1, 2, 3 and 4 (H1, H2, H3, H4, respectively; j ¼ 5, . . ., 8).These descriptors of the trajectories cover the spatial pattern of the movement, home behavior, and the impact of habitat on movement; thus describing the characteristics of our data with respect to our objectives in a comprehensive way.Wiegand et al. (2004) investigated the impact of the selection of summary characteristics on the resulting model predictions in detail and found that many summary statistics may be redundant, and that rejection filters work generally well if summary characteristics are selected which cover different aspects (patters) of the observed system at different levels of organization and scales.Each summary statistic excludes model parameterizations that produce clearly unlikely model behavior, but when simultaneous match of several summary statistics is demanded the model behavior will be substantially constrained (Grimm et al. 2005).At a later step we explore this aspect in more detail.

Assessing deviation of simulated data from observed data
In order to compare model predictions with observations, we created for each parameterization and model 500 simulated data sets that followed the field sampling protocol of the observed data.The movement of each individual was simulated in the corresponding habitat map v www.esajournals.org(site) starting from the coordinate where the individual was equipped with a transmitter.The model simulates the movement at an intraday resolution and, for every individual, the 500 simulated annual trajectories were sampled with the same protocol as used for field data (i.e., the same day and hour).
The procedure to define match of observed and simulated summary statistics is complex (Fig. B1 and Table B1 ).This is grounded in the need to consider the variability among individuals, uncertainty in the observation process, and process stochasticity in the simulated data.For each individual and parameterization, the comparison between observed and simulated trajectory was performed in a three-step approach.In a first step we developed two indices of deviance between the summary statistics j of the simulated and observed trajectory that account for the uncertainty present in one trajectory of individual i due to the observation process (i.e., the sampling protocol).In the second step we considered the stochasticity among simulations (i.e., process stochasticity in the model).The outcomes of the first two steps are two measures (m ij and v ij ) of the deviance between the summary statistics j of the observed and simulated data of individual i.The third step accounts for the observation error and the variability among the observed trajectories of different individuals to define the criterion to reject or accept parameterizations.The third step determines thresholds tm ij and tv ij for acceptance of summary statistic j of individual trajectory i (i.e., m ij , tm ij and v ij , tv ij ).Basically we accept a parameterization if the indices m ij and v ij describing the joint variability due to observation error and model stochasticity range within the variability caused by observation error and individual variability.
Observation error in one observed or one simulated trajectory.-Thevalues of the summary statistics calculated from the trajectories of individual tortoises i are subject to uncertainty related with the sampling protocol (i.e., the observation process).Sampling protocol means here the points in time when the location of an individual was recorded.For example, the resulting home range size may vary depending on the specific time used to record the locations of the animal.For the same reason there is also variability in the resulting home range size depending on the selection of locations actually used for its calculation.To assess the variability of a summary statistics j for individual i in the observed data due to the observation process we re-sampled the observed data leaving out randomly 20% of the locations and calculated the summary statistics.This re-sampling procedure was repeated 200 times, yielding, for each summary statistic, and individual i a distribution Dij of the values of the summary statistic j (Fig. B1).The variability in the summary statistics j of one simulated trajectory was then estimated analogously for a given parameterization and individual i.For a given simulation s the trajectory was then sampled following the field protocol.For this data set we obtained for each summary statistic j a distribution D sim ijs of the simulated values after 200 re-sampling steps (Fig. B1).To compare the observed distribution Dij with the simulated distribution D sim ijs we calculated the differences between the mean Mð Dij Þ À MðD sim ijs Þand the differences between the standard deviations Vð Dij Þ À VðD sim ijs Þof the distributions of the observed and simulated trajectory.
Process error.-Toaccount for process stochasticity we simulated for each individual i and a given parameterization 500 trajectories.For each of the 500 simulated trajectories we repeated the analysis explained in the last section.This yields 500 differences Mð Dij Þ À MðD sim ijs Þ and Vð Dij Þ À V ðD sim ijs Þ (Fig. B1).Because we were here primarily interested in the average behavior of the model we calculated the average values of MðD sim ijs Þ and VðD sim ijs Þ and compared it to the observed values, i.e., Our method would also allow to calculate a additional indices that considers the variability of the distribution of the simulated MðD sim ijs Þ and VðD sim ijs Þ, but to not overload our already complex procedure of model data comparison we did not use this information.However note that we consider the analogous variability of the ob- served data by using the full distributions DM j and DV j for defining the rejection criteria (see step 3).This is in some way analogous to approaches that use synthetic likelihood functions as measure of deviance between summary statistics of the observed and the simulated data (e.g., Wood 2010, Martı ´nez et al. 2011).Wood (2010) uses in his synthetic likelihood function the vector of the observed summary statistics, and the vector of the mean values and the covariance matrix of the summary statistics of the simulated data.Martı ´nez et al. (2011) uses in her synthetic likelihood function the mean values of the simulated summary statistics, and the mean and standard deviation of the observed summary statistics.The latter is more appropriate if we cannot expect that all competing models would able to match the observations well (i.e., there could be a strong structural error in some models) and if we are primarily interested in the mean behavior of the model.
The differences m ij and v ij are thus complex Process error, 500 simulated trajectories Complex measures of the match between the mean values of the distribution of the summary statistics of the observed and simulated data that consider observation and process error for a given individual i and summary statistic j.
Complex measures of the match between the standard deviation of the distribution of the summary statistics of the observed and simulated data that consider observation and process error for a given individual i and summary statistic j.Distribution of m j kl for all pairs of individuals, describes noise in the observed data with respect to the mean of the summary statistic j (the inter-individual variability and the stochasticity due to the sampling protocol).DV j Distribution of v j kl for all pairs of individuals, describes noise in the observed data with respect to the standard deviation of the summary statistic j (the inter-individual variability and the stochasticity due to the sampling protocol).tm j Threshold of DM j for summary statistic that gives the 95% range of the variability in summary statistic j due to observation error and individual variability.tv j Threshold of DV j for summary statistic that gives the 95% range of the variability in summary statistic j due to observation error and individual variability.m ij , tm j

Variability among individuals
The value of summary statistic j for individual i cannot be distinguished from data, given the observation error, process error, and individual variability.v ij , tv j The value of summary statistic j for individual i cannot be distinguished from data, given the observation error, process error, and individual variability.
v www.esajournals.orgmeasures of the deviance between the summary statistics of the observed and simulated data that consider observation and process error for a given individual i and summary statistic j.In case of good match we expect a low difference m ij in the means and a low difference v ij in the standard deviations.Thus, for every parameterization we obtained a total of N 3 8 3 2 indices (differences), where N is the number of individuals, 8 refers to the eight summary statistics and 2 refers to the mean and standard deviation.
Variability among individuals.-Weused the differences m ij and v ij presented in the last section to define agreement between simulated and observed data.To this end we conducted a statistical test to find out if the simulated and observed differences m ij and v ij were statistically different, given the observed variability among individuals.We therefore assessed if a difference that captures observation error and process stochasticity (i.e., m ij or v ij ) is within the observed range of the analogous difference that captures observation error and individual variability.To achieve this we treated the data from different individuals in the same way as we treated different trajectories of the simulated data (Fig. B1).We thus calculated for each summary statistic j the resulting absolute difference in the means jjMð Dkj Þ À Mð Dlj Þjj and standard deviations jjVð Dkj Þ À Vð Dlj Þjj for all pairs (k, l ) of individuals.This procedure yields for each summary statistic j the two distributions DM j and DV j of the differences in the mean and in the standard deviation, given the observed interindividual variability and the stochasticity due to the sampling protocol.We used a one tailed test with a p , 0.05 to define match in the summary statistic.This was the case if the values of m ij and v ij were located within the left 95% tail of DM j and DV j .(Fig. B1).Thus, a simulation was accepted if the differences m ij and v ij between observed and simulated summary statistics (that consider observation error and process stochasticity) range within the observed 95% range of the variability due to observation error and individual variability (captured by thresholds tm ij and tv ij , respectively).
The result of our procedure to define match between simulated and observed data is a binary matrix that tells us for a given model, parameterization, and summary statistic j and individual i whether the simulations successfully matched the observations or not.We considered that a given parameterization successfully predicted the movement of tortoises if this parameterization successfully matched 95% of all n 3 8 3 2 indices m ij and v ij .

DETERMINATION OF PARAMETER RANGES FOR THE DIFFERENT MODELS
We divided the range of each movement parameter to be estimated into 21 values (the 0 value plus 20 intervals).Autocorrelation and habitat weight parameters (AU1, AU2, H1W, H2W, H3W and H4W) ranged between 0 and 1, dHB between 10 and 10,000 m, whereas rHB ranged between 1 and 300 days.Combinations of parameters for simulations were obtained by means of the Latin Hypercube Sampling (McKay et al. 1979) which provides a stratified sample of the parameter space and which is more efficient than random sampling.However, to avoid unrealistic movement and consider prior knowledge, some combinations of parameters were excluded.To avoid excessive wandering of tortoises inside the movement step (i.e., too many cell steps) that is unrealistic, and to obtain a direct correspondence between the parameter ''distance per movement step'' (which is measured as straight line between localizations) and the movement step in the model we only considered parameterizations with high autocorrelation of cell steps (i.e., AU2 .AU1).In relation to habitat, because of prior knowledge (authors, unpublished data) we only considered those parameterization where H1W , H2W, H3W and H4W, that, is, the habitat quality of intensive land uses is lower than the other land uses (dry crops and natural areas).
Given the above restrictions, we selected for the full model (that contained 8 unknown parameters) 80,000 parameterizations, representing approximately 0.002% of the all possible parameterizations (n ¼ 4.25e 9 ).For the pure random walk model with no unknown parameters we use one (empirical) parameterization, and for the other 6 models with 2 to 6 unknown parameters we created 500 parameterizations (representing from 100% to 0.002% of the all possible parameterizations).

MODEL SELECTION AND ANALYSIS
To address our question regarding potential differences in movement at natural and disturbed landscapes we analyzed the alternative models separately for males and females and in natural and altered landscapes.One advantage of our approach is that we can study in detail in which summary statistic, and therefore property of movement, a given model failed.This provides important information on the processes governing tortoise movement.Note that approaches which use one compound measure of deviance between model simulations and observed data (e.g., a likelihood function) do usually not study the behavior of single summary statistics.To assess where certain models fail we selected for each model the 5% best parameterizations, based on the number of successfully predicted patterns, and calculated for each pattern the % of accepted parameterizations.
To quantify the similarity between the movement properties of males and females and in altered and natural landscapes, as described by our movement parameters, we calculated the overlap between the successful parameterizations for the model between males and females and in altered and natural landscapes (this is a type of cross validation).Finally, the inspection of the parameter values that successfully predicts the summary statistics allowed us to understand the effect of single parameters (and their interaction) in the model.scapes.In altered areas, the inter-individual variability in the summary statistics was much larger than that in natural areas.
The thresholds tm j and tv j that define for each summary statistic j the left 95% range of the distributions DMj and DVj were estimated for individuals in natural and altered landscapes (males and females separately) (Table C2).They are used in the next paragraph as criterion to reject or accept a given parameterization.Because individuals in natural landscapes behaved in a more homogeneous way than individual in altered landscapes, the thresholds for natural landscapes were more restrictive that for altered landscapes.Mean values of the summary statistics for the radiotracked males and females in natural and altered landscapes are given in Table C3.
To explore potential effects of year and landscape we checked the eight summary statistics for differences between individuals in natural landscapes from different years.This analysis was only made for females, where both years were adequately represented.We did not find statistical differences in the DMj values, although there were some differences in the DVj values (Table C4).Thus, the data from individuals in natural landscapes contain a certain amount of internal heterogeneity due to inter-annual variability.However, Table C4 also shows that differences between altered vs. natural landscapes were much larger than those between individuals in natural landscapes sampled in different years.Notes: These acceptance thresholds describe the level of observation error and individual variability in the summary statistics.We accepted a summary statistic j calculated from movement data of individual i if the indices m ij and v ij that describe the difference between observed and simulated data, considering observation error and process stochasticity, were smaller than the corresponding thresholds tm j and tv j , respectively.For summary statistic 1, all acceptance thresholds values are 0 since none of the radiotracked individuals entered in intensive land uses.See Appendix B and Table B1 and Appendix C for details.For a description of summary statistics, see Table 2. Notes: Each candidate model is represented by the best 5% parameterizations of the model.The number in the table gives for each summary statistic j the percentage of indices m ij and v ij that are matched (note that we calculated for each individual i and summary statistic two indices m ij and v ij ).rw: pure random walk, C: autocorrelation in movement, B: homing behaviour, H: habitat.The symbols of the summary statistics: home range size (HRa), home range perimeter (HRp), maximal displacement distance (Dmax), displacement after one year (Dy), proportion of localizations within the four habitat types (H1-H4).% ¼ Percentage of accepted parameterization of the model (i.e., match in more than 95% of all indices m ij and v ij ).For a description of the models see Methods: Movement models and Appendix B.

Fig. 1 .
Fig. 1.Empirical estimates of movement probabilities and step length.(A) Change of the probability PMOV to change position during one movement step (i.e., a 2-3 hour period) during the annual cycle, separately for males (light grey) and females (dark grey).(B) Probabilities DMOV(s) to move s cell-steps for males (light grey) and females (dark grey), averaged over the entire annual cycle.Note that we assessed the distribution of DMOV(s) for modeling purposes for each month separately.

Fig. 2 .
Fig. 2. Distribution of accepted parameterizations in the parameter space spanned by the parameters Home behavior distance (dHB, Y-axis) and Home behavior delay (rHB, X-axis) for different movement models and for males and females and in natural and altered landscapes .The random walk models: Brw: homing-biased, CBrw: correlated homing-biased, BHrw: homing-and habitat-biased, CBHrw: correlated homing-and habitat-biased.
These results have a clear biological interpretation.Males cover the largest distances when searching for females (Dı ´az-Paniagua et al. 1995, Anado ´n et al. 2006b), whereas females show their largest movements in the annual cycle searching for egg-laying places (Dı ´az-Paniagua et al. 1995).Because egglaying places are ubiquitous in the natural study area (authors, personal observations), the movement of female tortoises far from their usual activity areas are rare (Anado ´n et al. 2006b).

Fig. 3 .
Fig. 3. Frequency of the different weights (H2W, H3W, H4W) for the habitat types 2 (traditional agriculture land), 3 (flat natural areas) and 4 (natural areas on slope) in accepted parameterizations of the most complex model, the correlated homing-and habitat-biased random walk model (CBHrw) separately for males and females and in natural and altered landscapes.The weights give the probability that an individual moves during one cell-step to a cell of a given habitat type, relative to that of the eight cells surrounding the current location.

Fig. A1 .
Fig. A1.Land use maps of the five study sites (from left to right and from top to bottom): Galera (natural landscape), Bas, Palas, Quiñoneros and Sierrecica (altered landscapes).Study localities are a grid of 270 3 270 cells with a cell size of 10 3 10 m.

Fig. B1 .
Fig. B1.Schematic representation of the procedure to assess match between simulated and observed trajectories.
m j kl ¼ jjMð Dkj Þ À Mð Dlj Þjj Difference between the mean values of the Dkj of individual k and Dlj of individual l. v j kl ¼ jjVð Dkj Þ À Vð Dlj Þjj Difference between the standard deviation of the Dkj of individual k and Dlj of individual l.DM j

Fig
Fig. D1.Intra-step movement autocorrelation (AU2, Y-axis) vs. Inter-step movement autocorrelation (AU1, Xaxis) values in successful parameterizations.For a description of the models see Methods: Movement models and Appendix B.

Fig
Fig. D2.Intra-step movement autocorrelation (AU2, Y-axis) vs. Home behavior distance (dHB, X-axis) values in successful parameterizations in the Full model.Right: Inter-step movement autocorrelation (AU1, Y-axis) vs. Home behavior delay (rHB, X-axis) values in successful parameterizations in the Full model.For a description of the models see Methods: Movement models and Appendix B.

Table 1 .
Parameters employed in the movement model of T. graeca.

Table 2 .
Evaluation of movement models for females and males in natural and altered landscapes.

Table 3 .
Overlap between movement properties for males and females and on natural and altered landscapes.

Table A1 .
Landscape composition (% of land uses) of the five study sites.

Table B1 .
Summary of the procedure to compare simulated and observed data.

Table C2 .
Acceptance thresholds tm j and tv j for the summary statistics.

Table D1 .
Assessment of the ability of the different candidate models to match the eight summary statistics.