Extinction of canid populations by inbreeding depression under stochastic environments in Southwestern Goiás State : A simulation study

A frequently addressed question in conservation biology is what is the chance of survival for a population for a given number of years under certain conditions of habitat loss and human activities. This can be estimated through an integrated analysis of genetic, demographic and landscape processes, which allows the prediction of more realistic and precise models of population persistence. In this study, we modeled extinction in stochastic environments under inbreeding depression for two canid species, the maned wolf (Chrysocyon brachiurus) and the crab-eating fox (Cerdocyon thous), in southwest Goiás State. Genetic parameters were obtained from six microsattelite loci (Short Tandem Repeats STR), which allowed estimates of inbreeding levels and of the effective population size under a stepwise mutation model based on heterozygosis. The simulations included twelve alternative scenarios with varying rates of habitat loss, magnitude of population fluctuation and initial inbreeding levels. ANOVA analyses of the simulation results showed that times to extinction were better explained by demographic parameters. Times to extinction ranged from 352 to 844, in the worst and best scenario, respectively, for the large-bodied maned wolf. For the small-bodied crab-eating fox, these same estimates were 422 and 974 years. Simulations results are within the expectation based on knowledge about species’ life history, genetics and demography. They suggest that populations can persist through a reasonable time (i.e., more than 200 years) even under the worst demographic scenario. Our analyses are a starting point for a more focused evaluation of persistence in these populations. Our results can be used in future research aiming at obtaining better estimates of parameters that may, in turn, be used to achieve more appropriate and realist population viability models at a regional scale.


Introduction
Conserving biological diversity often requires researchers and environmental planners to face questions regarding the persistence of local populations.One of them is what is the chance of population persistence for a given number of years under certain conditions of habitat and human activity.This basic question motivated the creation and quick development of the research program in Population Viability Analysis (PVA) (Soulé, 1987;Morris and Doak, 2002).PVA is usually based on stochastic population models and is frequently used as a practical exercise in applied population ecology.It is based on real life history parameters of the target species or populations and on actual environmental features of the species habitat (Frankham et al., 2002;Morris and Doak, 2002).Thus, PVA allows an integrated analysis of genetic, demographic and landscape parameters that can provide an effective attempt to understand extinction risks and biodiversity loss (Pray et al., 1994;Britten, 1996;Henle et al., 2004).Genetic factors have always been recognized as an important component in evaluating population persistence (O'Brien et al., 1983;Lynch and Lande, 1993, Soulé, 1987, Lande, 1993, 1999), but the development of easier ways to assess population variability and structure was necessary before they became really useful for practical purposes (Frankham et al., 2002).Most of the research in the area focuses on the endogamy levels in natural or managed populations, on how these levels reflect in inbreeding depression and which components of fitness they affect (Barret and Charlesworth, 1991;Charlesworth and Charlesworth, 1987;Beisinger, 2000).These empirical estimates are usually the basis of population genetics models, which can then be combined with demographic and life history data to generate future scenarios at a landscape level using PVA.
The southwestern region of Goiás State, Brazil, dominated by the Cerrado biome, started to be agriculturally explored in the late 1930s, with the expansion of the agricultural activities from the neighboring states of São Paulo, Minas Gerais and Paraná.In the 1950-1960s, the implementation of governmental programs led to rapidly increasing agricultural and agro-industrial activities in the region, generating an inevitable high level of land conversion from natural into human-modified landscapes (Klink and Moreira, 2002;Klink and Machado, 2005;Marris 2005).The Parque Nacional das Emas (PNE), created in the core of this region in the decade of 1960, is currently the largest conservation unit of integral protection in the Cerrado biome, with a relatively high population density of many large-bodied mammals (Rodrigues et al., 2002).Although large-bodied mammals perceive the habitat at course grain and would disperse in these human-modified landscapes, these new landscapes pose difficulties for effective support of local populations and at least in part must limit their dispersion (e.g., Gehring and Swihart, 2003;Blanco et al., 2005) In this study, we simulated extinction by inbreeding depression under stochastic environments using Tanaka's (2000a) model to evaluate the time to extinction T e in two populations of Canids, the maned-wolf (Chrysocyon brachiurus) and the crab-eating fox (Cerdocyon thous), in southwestern Goiás State.We first conducted a population genetic analysis of these two species based on microsatellite (STR) data of specimens captured in the Parque Nacional das Emas and surroundings (Rodrigues, 2005).We used the estimated heterozygosis (He) to assess regional effective population size under a stepwise mutation model and to obtain the inbreeding coefficients.Other parameters of the simulations were obtained using a macroecological approach based on the overall knowledge of the region and species' life history.In the simulations, we assumed that combined effects of restricted dispersion, habitat loss (i.e., reduction in the carrying capacity of the region), and environmental stochasticity can drive regional populations to extinction.

Methodology
The genetic model Our population dynamics simulations were based on Tanaka's model (2000a) of extinction by inbreeding depression in stochastic environments (see also Tanaka, 1997Tanaka, , 1998Tanaka, , 2000b)).It assumes that inbreeding depression involves recessive deleterious genes, distributed among n diallelic autossomal loci and with similar mutation rates, selection coefficients and dominance.The model also assumes that effects of linkage disequilibrium and epistatic interactions are negligible.The genotypes AA, Aa and aa in each locus have mean fitness 1, 1 and 1 -s, where s is the selection coefficient.
Tanaka's model ( 2000a) is based on three dynamic recurrence equations for the mean allele frequency q, the inbreeding coefficient F and the population size N.These equations are: (1) In equation ( 1), μ is the mutation rate and γ is a random normal variate with mean 0 and standard deviation equal to (q (1 -q)/ 2nN) 1/2 , where n is the number of loci.Equation 3shows that population size N grows according to a discrete logistic model in time t, with carrying capacity K and fluctuating around the average N with a standard deviation ε.The Malthusian parameter (intrinsic growth rate) is given by r max , and thus the discrete increase in N by generation is λ, given by: However, population is also being disturbed by external events, so that we assume here that the carrying capacity K varies in time and is being reduced at a rate h per time t (in years, since this is the scale of λ), so that: The original model by Tanaka (2000a) was based on an exponential decrease, but we assumed a linear and constant decrease through time, generating a sharper decrease in K.Because of the relatively small intervals of N (see below), this change will make little difference in relation to the original model.Besides, population size is reduced by the factor δ, which is a function of allele frequencies q and by how deleterious genes affect the population growth rate.Thus, population growth is related to the genetic load L and is given as λδ, where:

Ecological context and model parametrization
We used the model described above to evaluate the population dynamics of the maned wolf and the crab-eating fox in the Parque Nacional das Emas and surroundings.The expected heterozygosis for each population, based on the genetic variation of six microsatellite (STR) loci (Rodrigues, 2005), was 0.533 and 0.676 for the maned wolf and the crab-eating fox, respectively.These values are within the expected range based on estimates performed with relatively close species worldwide (e.g., Wayne, 1996).We used a stepwise-mutation model (Nei and Kumar, 2000) to estimate effective population size with the following equation: where He is the expected heterosigosity, Ne is the effective population size and u is the mutation rate, equal to 10 -4 .We can thus approximate the effective population size of the maned wolf and the crab-eating fox in the region (not necessarily within the PNE) as Ne ≈ 500 and Ne ≈ 1000, respectively.Multiplying this number by a conservative factor of 2 (Morris and Doak, 2002) gives a reasonable initial estimate of population size N in the region, in the absence of better empirical estimates.Standard deviations in population size N were set to 5 or 10 individuals, so that this parameter will cause small effects at high population densities, but large effects at small densities, promoting an extinction vortex.
We simulated population dynamic processes in a region surrounding a relatively large conservation unit (the PNE), so we ran the simulations with two relatively low rates of decrease in carrying capacity K, equal to 0.005 and 0.01 individuals * year -1 .This also indicates that even if the PNE area is maintained for a long time, the carrying capacity of the population will still be reduced due to indirect border effects, such as the reduction of prey populations, pollution or loss of specific suitable habitats within the park.
There are no available growth rates data for the two species analyzed, and they were thus estimated with the following empirical equation (Alroy, 2001): where logM is the natural logarithm of the body mass of the species, equal to 24 kg for the maned wolf and to 6 kg for the crab-eating fox (Einsenberg and Redford, 1999).Thus, the maximum intrinsic growth rates were equal to 0.10 and 0.17, for the maned wolf and the crab-eating fox, respectively.
A large variation in the inbreeding coefficients F estimated for six STR loci was observed in both species.F ranged from -0.394 to 0.197 (average equal to -0.09) for the maned wolf and from -0.055 to 0.205 (average equal to 0.032) for the crab-eating fox.We thus ran the simulations with three low initial F-values, equal to 0.01, 0.05 and 0.10.Other parameters of genetic load were equal to those used by Tanaka (2000a), so that μ = 10 -4 , n = 15000 loci and the selection coefficient s equals 1.
Simulations were run for each species using 12 combinations of three varying parameters: 'low' and 'high' loss rates in carrying capacity K (0.005 and 0.01); 'low', 'medium' and 'high' inbreeding coefficients F (0.01, 0.05 and 0.10) and 'low' and 'high' standard deviations in population size (5 and 10 individuals) (Table 1).For each combination of parameters, 500 simulations were performed and the time to extinction T e was obtained.To evaluate how the combination of parameters used affect T e , a three-way factorial analysis of variance (ANOVA) (Sokal and Rohlf, 1995) was used.

Results
The average times to extinction T e for both species were affected by the three parameters analyzed and their interactions, according to the ANOVA, but T e was specially affected by K and SDN (Table 2).As expected, the worst scenario (i.e., shorter T e ) occurs for both species when populations fluctuate more, have highest inbreeding values and highest rates of decrease in carrying capacity (Table 1, scenario 12).On the other hand, the best scenario occurred in the opposite conditions (Table 1, scenario 1).For similar combinations of parameters, the decrease in the maned wolf population is always faster than that of the wild dog population, as expected by its 'slow' life history parameters (smaller initial population size N and lower growth rate r).

Simulation of extinction in canids from the Parque Nacional das Emas 123
Table 1 -Twelve simulation scenarios based on the combination of three variable parameters: Carrying capacity (K), standard deviation of population size (SDN) and inbreeding coefficient (F).For the maned wolf, T e values ranged from 291 to 413 years, in the worst scenario, and from 756 to 932 years, in the best scenario.The most important factors affecting T e in this species were K and SDN, and their interaction (Table 2), and the full ANOVA model explained 98% of the total variance in T e .

Scenarios
Similar results were found for the crab-eating fox.The T e values ranged from 359 to 486 years, in the worst scenario, and from 849 to 1099 years, in the best scenario.The most important factors affecting T e in this species were also K and SDN, and their interaction, but there was also a significant component for the interaction between F and K (Table 2).The ANOVA model also explained 98% of the total variance in T e .
Table 3 shows the detailed descriptive statistics for Te in both species and one simulation example in the worst and best scenarios for each species are shown in Figures 1 and  2. The population size decreases, due to both the reduction in carrying capacity and the effects of inbreeding depression.The average allele frequency q decreases quickly, whereas the F-values increase up to 1.0 (see also Tanaka, 2000a).

Discussion
The most important factors affecting the time to extinction in both species analyzed were those related to the demographic process, i.e., the relative magnitude of population fluctuations and the reduction in the carrying capacity.These factors have always been considered as the most important in driving populations to extinction (Frankham et al., 2002).The advantage of using Tanaka's model (2000a) is that genetic factors were also incorporated, contributing to reduction in population size, since the genetic load directly affects fitness.The absence of statistical significance of the levels of F-values in the model does not mean that these factors are not important, but just that the variation interval used was not large enough to generate variation in final T e , specially considering the stochastic variation around N. Neverthless, these values are probably close to real, considering the results obtained for other canid species in similar conditions (see Wayne, 1996;DeMattos et al., 2004) and our own analysis with six STR loci.
It is important to consider that some studies showed that large-bodied mammal populations can have relatively low genetic loads due to past bottleneck events that filter deleterious variability, thus being less affected by inbreeding (Vilà et al., 2002;Merola, 1994;O'Brien et al., 1983;Sage and Wolff, 1986;Avise, 2004).
The comparative analyses of the model parameters for both populations are in agreement with expectations based on demographic parameters and life history differ-124 Rodrigues and Diniz-Filho  ences among species (e.g.Wayne, 1996).The maned wolf is a large-bodied canid with a large home range and a relatively low growth rate.Although it can disperse more effectively and is probably less sensitive to habitat loss and fragmentation, in some circumstances this absence of a strong population structure could increase the possibility of regional extinction (Henle et al., 2004).On the other hand, in small-bodied species, such as the crab-eating fox, growth rates are usually higher, which generates more population fluctuations.When associated with higher habitat specialization (not the case with this species), this can generate more local extinction, which is usually counteracted by recolonization processes.
It is important to stress that the simulations that we performed should be viewed as mere initial evaluations of extinction risks and population viabilities.There are several reasons for that.First, there is no information on detailed genetic load for these populations and how this can be linked with basic population dynamics parameters, such as mortality rate and fecundity.Secondly, we derived the initial population size N based on the expected heterozygosis estimated from a few microsatellite loci and the relation between effective and observed population sizes is unknown.The estimated population sizes provide approximations of populations in equilibrium in an evolutionary context.Recent events of population reduction and mortality, such as those caused by recent habitat loss and hunting in the last 50 years of intense human occupation in the region, could have occurred without changing the estimates based on He because their effects on genetic data are usually delayed.Finally, our model assumes that only habitat loss (generating reduction in K) affects genetic parameters, whereas in fact more complex extinction patterns, such as critical thresholds, can appear because of the interaction between habitat loss and habitat fragmentation.
Our model is based on the persistence of a relatively large population at a regional scale, which is not unrealistic for the species studied considering their geographic range (see Eisenberg and Redford, 1999).We could also consider more drastic scenarios, specifically for southwestern Goiás State, by, for example, setting that core remaining populations of both species are concentrated in PNE.Rough estimates of abundance by people performing field work in PNE suggest that population sizes for both species are at least one order of magnitude lower than the one we used in our simulations (Leandro Silveira, personal communication).Starting with a population of 100 individuals of maned wolves in PNE, for example, the Te would be 126 years in average in the worst scenario and 364 years in the best scenario.
Our simulations showed that the populations of maned wolf and crab-eating fox in southwestern Goiás can persist for a reasonable time (i.e., for over 300 years) even in the worst possible scenario.These are initial theoretical projections and further detailed analyses of these populations should be conducted.Nevertheless, our analysis is a starting point for a more focused evaluation of persistence in these populations and can drive future research towards obtaining better estimates of parameters that can, in turn, be used to achieve more appropriate and realistic models.

Figure 1 -
Figure 1 -Results of one simulation for the crab-eating fox population of the Parque Nacional das Emas.The variation in F-values (A) and in population size (B) through time in the two extreme scenarios are shown.The solid line shows the trajectory of the parameters in the best scenario and the dashed line indicates the trajectories in the worst scenario.

Table 2 -
F-values and their significance levels, from the factorial ANOVA applied to understand the components affecting time to extinction (Te) in the two species analyzed, including variation in carrying capacity (K), in standard deviation of population size (SDN) and in the inbreeding coefficient (F) and their first and second order interactions.

Table 3 -
Time to extinction for the maned wolf (Cb) and the crab-eating fox (Ct) in 12 simulation scenarios (see Table1).