Does Litter Size Variation Affect Models of Terrestrial Carnivore Extinction Risk and Management?

Background Individual variation in both survival and reproduction has the potential to influence extinction risk. Especially for rare or threatened species, reliable population models should adequately incorporate demographic uncertainty. Here, we focus on an important form of demographic stochasticity: variation in litter sizes. We use terrestrial carnivores as an example taxon, as they are frequently threatened or of economic importance. Since data on intraspecific litter size variation are often sparse, it is unclear what probability distribution should be used to describe the pattern of litter size variation for multiparous carnivores. Methodology/Principal Findings We used litter size data on 32 terrestrial carnivore species to test the fit of 12 probability distributions. The influence of these distributions on quasi-extinction probabilities and the probability of successful disease control was then examined for three canid species – the island fox Urocyon littoralis, the red fox Vulpes vulpes, and the African wild dog Lycaon pictus. Best fitting probability distributions differed among the carnivores examined. However, the discretised normal distribution provided the best fit for the majority of species, because variation among litter-sizes was often small. Importantly, however, the outcomes of demographic models were generally robust to the distribution used. Conclusion/Significance These results provide reassurance for those using demographic modelling for the management of less studied carnivores in which litter size variation is estimated using data from species with similar reproductive attributes.


Introduction
Demographic variation, resulting from extrinsic and intrinsic sources, fundamentally affects population dynamics and is particularly important when assessing extinction risk for threatened species [1,2]. Predictions of population dynamics depend on the ability to attribute sources of stochasticity accurately in population models [3,4]. Of particular importance is the distinction between demographic stochasticity and demographic heterogeneity. Demographic stochasticity is the random fate of an individual arising from a chance event drawn from a specified uniform vital rate, whereas demographic heterogeneity is the variation in the underlying parameter value arising from withinpopulation variability in individual condition [5]. Both types of demographic variation make important contributions to a populations' total demographic variance [3]. Indeed, accounting for demographic stochasticity in fecundity can lead to increased predictions of extinction risk; for example, overall demographic variance is increased when this parameter is Poisson-distributed [5]. Here, we focus on stochasticity in demographic fates, which can easily be accounted for by drawing rates from appropriate probability distributions [6,7].
Mean litter (or clutch) size has long been the focus of evolutionary and population biologists concerned with causes of interspecific variation [8][9][10][11], correlations with environmental gradients [9,[12][13][14] and optimality in this trait [15][16][17][18]. However, intra-population variation in litter size has been largely overlooked (but see [19]). Limited knowledge of the underlying measures of empirical litter size distributions, such as the degree of dispersion, hinders the accurate representation of the stochasticity of this parameter in population models. Demographic stochasticity in offspring number is most commonly modelled with Poisson or normal distributions [6,20,21], although there is little theoretical justification for these choices [19]. Furthermore, many demographic modelling programmes (e.g. RAMAS [7] and VORTEX [21]) have limited provision for specifying distributions. Unlike survival, which is a Bernoulli process [20], choosing a distribution to describe variation in litter sizes in multiparous species can be complex because the biology of reproduction differs substantially among species and is ultimately limited by physiological capacity. Standard probability distributions might lack the flexibility required to account for litter size variation in many species.
In population modelling, the influence of distribution choice has only been considered previously for demographic parameters other than litter size, with a focus on environmental stochasticity. Studies that modelled environmental stochasticity found that population growth rate (l) estimates were underestimated as a result of inaccurately defined, symmetrical survival distributions [22] and large differences in l estimates were found when drawing recruitment rates from different distributions [23]. Yet, the shape of the distribution may also be important for populations that are susceptible to fluctuations in vital rates as a result of demographic stochasticity, such as small populations. Failing to account for demographic stochasticity in litter size may lead to inaccurate predictions of extinction risk [19]. In this context, it is useful to establish whether failing to incorporate an appropriate theoretical distribution for litter size, to describe demographic stochasticity, could lead to erroneous estimates of model outputs.
Here, we examine the fit of specified candidate probability distributions to empirical data on terrestrial carnivore litter size frequencies. The Carnivora exhibit some of the most diverse life history traits of all mammalian orders, as reflected in their broad range of litter sizes [24]. While many carnivores are at increasing risk of extinction [25], others are predators of economic importance or important hosts of zoonotic and wildlife diseases such as rabies [26]. Although data collection is often challenging [27], both categories of carnivore are frequently the subject of population models (e.g. [28][29][30]). Given the importance of carnivore management and the sparseness of much of the data used to model carnivore demography, it is useful to establish whether the choice of distribution used to model demographic stochasticity in litter sizes affects the inferences drawn from models of carnivore population dynamics. To illustrate the applied importance of using appropriate distributions, three previously published population models are replicated to determine the consequences of mis-specifying litter size distributions for inferences regarding extinction probabilities or disease dynamics.

Probability distribution fitting
Litter size frequency data were collated for 32 terrestrial multiparous carnivore species, from 63 published studies of 73 wild populations, to reflect the diversity of life history within the order. Each species has a single annual breeding attempt. None of the studies included litters of zero; modelling litter size inherently assumes that an individual has bred. If studies presented data for multiple conspecific populations or for multiple methods of litter size determination, these were analysed as discrete datasets. For 15 species, data were obtained for between two and ten populations. For three species, data from multiple methods of litter size determination (e.g. placental scars and direct counts) were available. We thus also considered whether there was strong support for genuine underlying difference in litter size distributions between conspecific populations or between data determined by different methods (again, for a given population) (see Appendix S1 for details of the analyses).
Twelve probability distributions were selected based on a review of previous studies. Specifically, four discrete distributions were chosen: the Poisson distribution [6]; the generalised Poisson, which has a wide-ranging suitability for describing litter size frequencies [19]; the binomial distribution, previously fitted successfully to carnivore litter data [19]; and the negative binomial, widely used to describe ecological processes (e.g. [31]). For each discrete distribution, both a 'right shifted' and 'zerotruncated' form were fitted (Appendix S2), to exclude litter sizes of zero. For zero-truncation, the probability mass function was scaled by the exclusion of predicted zeros. Shifting involved moving the entire distribution one interval to the right. Three continuous probability distributions were chosen: the normal and lognormal distributions are both widely used [6], although logtransformation is not recommended for count data [32]; and the stretched beta (two and three parameter forms), as proposed by Morris and Doak [6]. Appendix S2 provides details of how these continuous distributions were converted into discrete forms.
Maximum-likelihood parameters, denotedĥ h, were estimated using the 'optim' function in R 2.14.0 (R Development Core Team 2011). Here, the multinomial log-likelihood defined by h and given all the data is: where N is the total number of litters observed, N i is the number of litters observed of size i, P i is the predicted litter size probability determined by a given distribution (Appendix S2), x max is the maximum litter size, and C(x) is the complete gamma function. The fits for each probability distribution were compared using Akaike's Information Criterion (AIC); all distributions having a DAIC#6 of the best fitting distribution (i.e. lowest AIC) were considered to have some support [33]. To check that our bestfitting models were consistent with the data, and because of the small sample sizes of the predicted frequencies, we performed goodness-of-fit tests using Fisher's Exact Test. Variance-mean ratios [34] were determined to measure the dispersion of the empirical and fitted distributions.

Carnivore population models
Published stochastic population models for three management scenarios were used to illustrate the broader applied significance of this study. The Canidae were chosen because they provide the widest range of litter sizes within the Carnivora [24]. Models were chosen to depict a range of conservation and management scenarios that could be replicated from published data; the intention was to identify whether the choice of distribution used to represent litter sizes influences predicted model outcomes. By ''outcomes'', we refer to a major emergent parameter from the models, on which further inference would be based (see below). The emergent parameter of interest varied because the three models were created for different applications. Using the parameters that were estimated by maximum likelihood as described above, 10,000 stochastic replicates of the models were simulated drawing litter sizes from each of the 12 probability distributions. This enabled calculation of 95% confidence intervals around mean outcome values. For each case study, disparities were determined between the outcome values of the 12 model versions. This allowed us to evaluate the effect on each model of employing different litter size distributions, in relation to the degree of empirical support for those distributions. See Appendix S3 for full descriptions of each case study model and Table 1 for the initial parameter values for the three models.
First, we investigated the island fox Urocyon littoralis, which reached near extinction on Santa Catalina Island due to an outbreak of canine distemper virus [35]. We conducted a densitydependent population viability analysis (PVA) for two subpopulations, based on Kohlmann et al. [28]; the outcome of interest was the probability of quasi-extinction, defined in this model as the probability of the population declining to 50 individuals, due to a disease epidemic. Second, we investigated the red fox Vulpes vulpes, a locally abundant carnivore that is the focus of much attention due to its economic importance as a predator and role in the spread of rabies [36]. A density-dependent model simulating control after a rabies outbreak [29] was replicated to illustrate, as the outcome of interest, the probability of successful disease control. Finally, we investigated the African wild dog Lycaon pictus, which is restricted throughout much of its range and susceptible to several diseases, including rabies [37]. A density-dependent PVA [30] for small wild dog populations was reproduced to determine quasi-extinction probabilities (the outcome variable), defined here as the probability of only one sex remaining. Following Vial et al. [37], we also investigated the effects of including a component Allee effect (a positive relationship between population size and a measurable component of fitness [38]) with respect to recruitment. Here, rather than reducing pup mortality, individual litter size was assumed to be an increasing function of group size, sensu Vial et al. [37]. These three investigations illustrate canids with small, medium, and large mean litter sizes, respectively (Table S1). All modelling and analyses were conducted in R 2.14.0 (R Development Core Team 2011).

Carnivore population models
The demographic modelling showed that the distribution chosen to represent litter size uncertainty in the three canid models has limited impacts, regardless of the fit of the distributions. PVA models for island foxes showed that estimating extinction probability was largely unaffected by the choice of distribution, with less than 1% difference in quasi-extinction probabilities between models that used the best and worst fitting litter size distributions ( Fig. 2A, B). Similarly, regardless of whether the litter size distributions used in the model provided a good fit to empirical litter size data, there was only a 2% difference in the probability of successful disease control in the rabies model for red foxes (Fig. 2C, D). Likewise, quasi-extinction probabilities for The number of datasets tested for each species (denominator, see Table S1 for details) and the number of datasets that were adequately fitted by a given distribution (numerator, see Table S2 for   African wild dogs showed only a 1% difference among models that employed different litter size distributions (Fig. 2E, F). When litter size was reduced as a function of group size, to simulate an Allee effect, the influence of the distributions was slightly greater (Fig. 2G, H), with an increase of approximately 4% between quasiextinction probabilities for the best and worst-fitting distributions. Even in this case, only models employing the worst-fitting distributions differed substantially in their predictions from those of models employing other distributions. The variation in the skew and variance of the fitted distributions (Fig. 2) may be attributed to process and sampling error in the data, as well as properties of the distributions such as the tendency to favour overdispersion, e.g. the negative binomial. However, for all parsimonious distributions, these measures were generally consistent with the empirical distributions for all models except island foxes (Fig. 2). In this latter case, the variation in agreement between distributions with DAIC #6 and the empirical properties ( Fig. 2A, B) is probably due to the small sample size increasing the uncertainty of the observed parameter estimates, translating into the selection of multiple distributions. Despite the widely varying variance, the resultant model outcomes were in general unaltered by the choice of distribution. Coefficients of variation (CV) were small for all model outcomes (Table S4), with the greatest variation in the African wild dog model with an Allee effect; the best-fitting distribution (CV = 0.712) was 1.07 times more variable than for the worst fitting model (CV = 0.668).

Discussion
Multiple distributions were shown to be consistent with the data for describing litter size frequencies for a range of carnivore species. However, the outcomes of demographic models appear robust to the choice of litter size distribution. These findings are discussed in light of the biological implications of litter size distribution choice and the applied importance of incorporating suitable probability distributions in demographic models.

Model selection for describing litter size variation
Unlike many biological parameters, offspring number is often underdispersed [39,40] and positively skewed [41,42]. Litter size frequencies are best fitted by probability distributions able to describe the biological constraints on the upper limit of offspring production. While the Poisson distribution is most commonly used for fitting count data in general, it does not allow for underdispersion. In contrast, the generalised Poisson separates the variance from the mean [19], allowing greater flexibility, but at the cost of additional parameters. Of the continuous functions, the discretised normal distribution is the most flexible and is suitable for data characterised by low variance.
In a recent model of vertebrate reproductive success, the zerotruncated generalised Poisson was consistently the best-fitting of several parametric distributions fitted to litter size [19]. However, that study only included one carnivore population, (lion, Panthera leo), which was fitted solely by the zero-truncated-binomial. In our study, that distribution performed less well, perhaps because more competitive functions were considered (including shifted discrete distributions and discretised continuous distributions) that were not assessed in the earlier study [19]. The better fit of shifted forms over zero-truncation suggests that further work is needed to determine whether there is an underlying probabilistic mechanism in the distribution of litter size.
The lack of evidence for intraspecific variation in underlying litter size distributions (Appendix S1) could indicate that biological limitations on reproduction allow for little intraspecific variation in this trait. The known biases associated with litter size determination methods for red foxes [43,44] probably explain the observed differences in litter size distributions (Appendix S1), although the results of the management scenarios analysed in this study (see next section) suggest that this finding is unlikely to be of consequence for future modelling efforts. Given the pooling of litter size datasets in this study over multiple years, due to insufficient data, the results must be interpreted with caution in light of potential temporal variation.
These analyses assumed that individuals had the same underlying expected reproductive capacity. However, demographic heterogeneity in offspring production is influenced by many factors, including female age, body condition or social status [45,46], as well as maternal versus offspring trade-offs in reproductive success [47,48]. The methods in these analyses could be incorporated into population models that address such intrinsic individual variation, as well as those modelling environmental stochasticity.

Applied importance of litter size distributions
Despite interspecific variability in the consistency of distributions to describe litter size data, we have shown that model outcomes of applied management scenarios, e.g. extinction risk, may be robust to the distribution chosen to represent litter sizes. The lack of any apparent effect of litter size distribution choice in carnivore models might be because mammalian litter sizes are generally small due to physiological limitations. Underdispersion will promote sampling of offspring closer around the mean; therefore, sampling variation will only weakly impact model outcomes. There are indications that the distribution choice could be important in limited circumstances. In the case of African wild dog populations, the example presented here illustrates how modelling a component Allee effect in reproduction using an illfitting, underdispersed distribution can result in an overestimation of extinction risk (see Fig. 2E-H).
Further work is required to determine the potential influence of temporal variation in the underlying litter size distribution on predictions of extinction risk. This is particularly important given that temporal or environmental variability means that combining data over time will inflate estimates of litter size variation, leading to erroneous predictions of extinction risk. In spite of these concerns, the lack of available data meant that pooling data was necessary for our purposes; consequently, our results are indicative only of how mis-specified distributions could affect model predictions. As in [19], we stress that determining appropriate distributions is a step towards a more mechanistic understanding of litter size variability that could provide insight into a species' response to selective pressures or management actions.
That litter size distributions have limited effects on the outcomes of management models may also reflect the relative contributions of life history traits to population growth. For long-lived species such as carnivores [49], the elasticity of adult survival typically contributes more to population growth than fecundity. Indeed, variance in demographic parameters with low elasticities will have little effect on the variance of the population growth rate, due to the near linear relationship between population growth and vital rates [50]. Notably, for all three canid populations in the models presented here, the elasticity of survivorship is as high or higher than fecundity [28,30,51], which is consistent with the limited impact of litter size variation observed in the case studies.
Although this study focused on the Carnivora, our findings should apply to taxa with multiparous females, including other mammals, birds and lizards. While it is hard to determine the exact ecological and physiological mechanisms generating a litter size distribution, insight into the drivers of these empirical distributions could aid our understanding of the adaptation of reproductive strategies to extrinsic and intrinsic population pressures. Recent work demonstrating that female red foxes exhibit sex-biased investment in offspring as a function of body mass and population density suggests that altering litter size composition rather than litter size could be an alternative mechanism for increasing fitness [52]. Ultimately however, applied models for carnivores appear to be robust to choice of litter size distribution, which has positive implications for modelling species with limited data.    Appendix S1 Testing for intraspecific variation in litter size distributions, using the red fox Vulpes vulpes as an example.

(DOC)
Appendix S2 Functional forms for the 12 probability distributions fitted to empirical litter size frequency data.

(DOC)
Appendix S3 Model descriptions for the three canid management scenarios used to illustrate the consequence of using different distributions to model litter size. (DOC)