The relative effect of past climate fluctuations and anthropogenic activities on current biome distribution is subject to increasing attention, notably in biodiversity hot spots. In Madagascar, where humans arrived in the last ~4 to 5,000 years, the exact causes of the demise of large vertebrates that cohabited with humans are yet unclear. The prevailing narrative holds that Madagascar was covered with forest before human arrival and that the expansion of grasslands was the result of human-driven deforestation. However, recent studies have shown that vegetation and fauna structure substantially fluctuated during the Holocene. Here, we study the Holocene history of habitat fragmentation in the north of Madagascar using a population genetics approach. To do so, we infer the demographic history of two northern Madagascar neighbouring, congeneric and critically endangered forest dwelling lemur species—Propithecus tattersalli and Propithecus perrieri—using population genetic analyses. Our results highlight the necessity to consider population structure and changes in connectivity in demographic history inferences. We show that both species underwent demographic fluctuations which most likely occurred after the mid-Holocene transition. While mid-Holocene climate change probably triggered major demographic changes in the two lemur species range and connectivity, human settlements that expanded over the last four millennia in northern Madagascar likely played a role in the loss and fragmentation of the forest cover.
Inferring the demographic history of species is crucial for understanding the evolutionary forces that shape genetic diversity (Keinan & Clark, 2012; Mitchell-Olds, Willis, & Goldstein, 2007). Genetic data can shed light on major (pre)historic events such as bottlenecks, expansions, migrations and admixtures that affected species at a regional or continental scales (e.g., Hewitt, 2004; Provan & Bennett, 2008). Additionally, inferring the demographic history can help to disentangle demographic from adaptive effects (Goldstein & Chikhi, 2002; Kelley, Madeoy, Calhoun, Swanson, & Akey, 2006; Nielsen, Hellmann, Hubisz, Bustamante, & Clark, 2007).
In Madagascar, as in other regions of the world, the relative importance of ancient human and climatic impacts on wildlife is still hotly debated (Faith, 2014; Godfrey & Irwin, 2007; Goodman & Jungers, 2014; Koch & Barnosky, 2006; Quéméré, Amelot, Pierson, Crouau-Roy, & Chikhi, 2012; Sandom, Faurby, Sandel, & Svenning, 2014; Stuart, 2015; Yoder et al., 2016). Human pauperization, economic activities and population growth are usually considered to be the main causes for the degradation of natural habitats in Madagascar. However, each region has its own history and human practices were not always the original drivers of forest loss and fragmentation and of species extinctions (Stuart, 2015). Reconstructing the demographic history of species from genetic data to tease anthropogenic and climatic factors apart is an attractive approach, which can complement historical, archaeological and palynological information (e.g., Agudo, Rico, Vilà, Hiraldo, & Donázar, 2010; Allentoft et al., 2014; Salmona et al., 2012). Furthermore, the demographic history of populations may provide information about species and habitat resilience to climate change and can therefore guide conservation management (Hoffmann et al., 2015; Shafer et al., 2015b).
The past 20 years has seen the tremendous development of methods for inferring the demographic history of populations. The first methods relied on the deviation of summary statistics from expected values under an equilibrium demographic model (e.g., Garza & Williamson, 2001; Luikart & Cornuet, 1998; Tajima, 1989). These methods were soon followed by likelihood or Bayesian approaches that allow estimating current and past effective population sizes, as well as dating demographic events (Beaumont, 1999, 2003; Storz & Beaumont, 2002). In recent years, more complex and realistic model-based approaches have emerged allowing the inference of several population size changes from sequence or microsatellites data (e.g., Heled & Drummond, 2008; Leblois et al., 2014; Nikolic & Chevalet, 2014; Wu & Drummond, 2011), and genomic data (e.g., Excoffier & Foll, 2011; Gutenkunst, Hernandez, Williamson, & Bustamante, 2009; Li & Durbin, 2011; Liu & Fu, 2015). Several of these methods are, however, limited by the difficulty of computing the likelihood function for large data sets or complex models, thus restricting their scope to simple evolutionary and molecular models which might not capture all relevant processes in complex demographic histories (Marjoram & Tavaré, 2006). For instance, inferences assuming oversimplistic demographic scenarios may lead selection signals to be confounded by demographic processes (Nielsen et al., 2005), and population structure to confound population size change inference (Beaumont, 2004; Chikhi, Sousa, Luisi, Goossens, & Beaumont, 2010; Heller, Chikhi, & Siegismund, 2013; Mazet, Rodríguez, & Chikhi, 2015; Peter, Wegmann, & Excoffier, 2010; Wakeley, 1999).
Modelling realistic habitat loss and fragmentation (HL&F) scenarios is particularly challenging because HL&F is a spatio-temporal process, which requires modelling changes both in population size and structure over time. The flexible approximate Bayesian computation (abc, Beaumont, Zhang, & Balding, 2002; Pritchard, Seielstad, Perez-Lezaun, & Feldman, 1999; Tavaré, Balding, Griffiths, & Donnelly, 1997) framework bypasses the difficulty of deriving likelihood functions from complex demographic models (Beaumont, 2010; Csilléry, François, & Blum, 2012; Sunnaker et al., 2013) and hence allows for population genetics inferences under such models. abc relies on simulated data from a set of models which are compared using the distance between simulated and observed data (usually using summary statistics). Finally, the posterior distribution of the model parameters can be approximated (Beaumont, 2010; Cornuet et al., 2008; Csilléry et al., 2012). It has been successfully applied to model the complex history of orangutans affected by episodes of HL&F (Nater et al., 2015), to infer the introduction history of macaques in Mauritius (Bonhomme, Blancher, Cuartero, Chikhi, & Crouau-Roy, 2008), to study the invasion histories of the bumblebee in New Zealand (Benazzo, Ghirotto, Vilaça, & Hoban, 2015), of rats in Madagascar (Brouat et al., 2014) and olive trees in Australia (Besnard et al., 2014b). However, in most population genetic inference studies including those using abc, population structure is often ignored, even when the habitat is clearly structured or fragmented.
HL&F has become a major concern worldwide due to its impact on biodiversity (Gibson et al., 2011; Laurance et al., 2012; Pimm & Raven, 2000). In Madagascar, one of the “hottest” biodiversity hot spots which harbours high species richness and endemism levels (Ganzhorn, Lowry, Schatz, & Sommer, 2001; Myers, Mittermeier, Mittermeier, Da Fonseca, & Kent, 2000), deforestation alone was estimated to have caused the loss of 9% of Malagasy plant and invertebrates species between 1950 and 2000 (Allnutt et al., 2008). Following Madagascar's colonization by humans, ~4 to 5,000 years BP (Crowley 2010; Dewar et al., 2013; Gommery et al., 2011), it has been suggested that a major part of its original forest cover was lost (Humbert, Darne, & Besairie, 1965) with only 10%–20% of Madagascar being forested today (Moat & Smith, 2007). The situation is alarming if we consider that more than 90% of the Malagasy species (including lemurs) live exclusively in forests and woodlands (Dufils, 2003; Goodman & Benstead, 2005). This narrative describing Madagascar as covered by woodlands when humans first arrived (Humbert, 1927; Perrier de La Bâthie, 1921) continues to be prevalent in management and other capacities, although it is not uncontroversial (Kull, 2000).
For instance, early descriptions (Gautier, 1902; Grandidier, 1898) and recent records (Burney et al., 2004; Gasse & Van Campo, 1998; Matsumoto & Burney, 1994; Virah-Sawmy, Willis, & Gillson, 2010) provide evidence that some regions of the island were covered by a mosaic of forests, scrublands and grasslands long before the first large human settlements (i.e., 1,000 year ago, Dewar & Wright, 1993) started to significantly impact the landscape (Burney, Robinson, & Burney, 2003; Burns et al., 2016; Gasse & Van Campo, 1998) and before the earliest evidence of human presence 4–5,000 yBP (Dewar et al., 2013; Gommery et al., 2011). Palaeontological records as well as genetic studies suggested that recent climatic events were responsible for decreases in wild mammal populations and habitat contraction in northern Madagascar (Jungers, Godfrey, Simons, & Chatrath, 1995; Quéméré et al., 2012; Rakotoarisoa, Raheriarisena, & Goodman, 2013; Simons et al., 1995). These findings contradicted the dominant narrative, but brought new exciting insights into the complex history of northern Madagascar lansdcape.
It is, however, undisputed that recent decades have brought significant human-mediated HL&F. Conservative estimates show that at least 52% of the forest cover was lost since the 1950s (Harper, Steininger, Tucker, Juhn, & Hawkins, 2007; ONE 2013; Schwitzer et al., 2014), from the use of fire, slash-and-burn cultivation (“tavy”), cattle raising, logging (Minten, Randrianarisoa, Randrianarison, & Food, 2003) and mining activities.
In this study, we investigate the HL&F history of two neighbouring and congeneric lemur species, Perrier's and Tattersall's sifaka (Propithecus perrieri and Propithecus tattersalli) to identify the historical drivers of wildlife dynamics (human vs. climatic drivers) in northern Madagascar. To reach these objectives, we use a comparative modelling approach of the recent demographic history of these two species. Our study capitalizes on four major advantages. We first benefit from a near-complete sampling of the distribution range of each species with detailed and reliable provenance. We further use various complementary modelling approaches as well as previous population genetic knowledge to build and test realistic demographic models for each species. Third, we use models incorporating population substructure to disentangle the effect of changes in population size from potential confounding effects. Finally, we compare alternative “temporally delineated” hypotheses to disentangle the potential climatic and anthropogenic effects on population decline and of HL&F.
2 MATERIAL AND METHODS
2.1 Study species
Sifakas (genus Propithecus) can be divided into two phylogenetic lineages: dry forest smaller sized sifakas to which P. tattersalli belongs and humid forest larger sized sifakas, to which P. perrieri belongs. Species from these two groups present wide, parapatric, and contiguous ranges along the western and eastern coast of Madagascar, except Perrier's and Tattersall sifakas that have a restricted distribution disjunct from the northern limit of their sister species (Fig. S1). These two tree-dwelling species are organized in matriarchal social groups of two to eight individuals (Lehman & Mayor, 2004; Meyers & Wright, 1993) and live in fragmented forests frequently connected by riparian corridors (Mittermeier et al., 2010; Quéméré et al., 2010a; Salmona et al., 2013; Figure 1). They typically disperse only short distances, but are also known to occasionally cross a large distance of open habitat (Mayor & Lehman, 1999; Meyers & Wright, 1993; Richard, Rakotomanga, & Schwartz, 1993). They show the smallest ranges and some of the lowest population sizes of any sifaka species. As a consequence, they are both considered Critically Endangered (CR; IUCN 2014), and Perrier's sifaka has even been listed in “The 25 most endangered primates of the world” on several occasions (Banks, Patel, Chikhi, & Salmona, In press, Banks, Patel, Chikhi, & Salmona, 2015). Subfossil data suggest that the two species may have had wider and perhaps even sympatric palaeodistribution (Godfrey, Jungers, Simons, Chatrath, & Rakotosamimanana, 1996; Jungers et al., 1995), illustrating the substantial influence of climate oscillation and refugia on today's biodiversity makeup (Wilmé, Goodman, & Ganzhorn, 2006; Wilmé et al., 2012).
2.2 Sample collection and DNA analysis
Faecal material from 244 P. tattersalli and 42 P. perrieri were collected during the dry season (April to October) from 2006 to 2013 over the entire range of the species and stored in dry condition with silica gel (Quéméré, Louis, Ribéron, Chikhi, & Crouau-Roy, 2010c; Figure 1). Individuals were genotyped with 13 and 24 microsatellites loci, respectively (Quéméré, Crouau-Roy, Rabarivola, Louis, & Chikhi, 2010b; Salmona et al., 2015). Field and laboratory procedures are described in previous studies (Quéméré et al., 2010b; Salmona, 2014; Salmona et al., 2015). Briefly, DNA was extracted using the 2-CTAB/PCI method (Vallet, Petit, Gatti, Levréro, & Ménard, 2008), and microsatellites markers were genotyped using a sequential repetition approach to ensure genotype accuracy as described in Quéméré et al. (2010c) and Salmona et al. (2015). Typical quality control applied to check for sample duplicate, null alleles, allele dropout and linkage disequilibrium is detailed in Quéméré et al. (2010b) and Salmona et al. (2015). All sample genotypes and geographic coordinates are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.8f45n.
2.3 Demographic history
The demographic history of both species was preliminarily investigated using two approaches, implemented in msvar1.3 (Storz & Beaumont, 2002) and vareff (Nikolic & Chevalet, 2014) that allow the detection, quantification and dating of changes in effective population size based on allelic frequency distributions. Both approaches make the assumption that population structure is negligible. From these “preliminary” analyses and from previous spatial analyses of Perrier and Tattersal's genetic diversity (Quéméré et al., 2010b; Salmona et al., 2015), we constructed more complex demographic scenarios incorporating population structure, habitat fragmentation and population size changes and compared them under an approximate Bayesian computation (abc) framework.
2.3.1 Generation time
There is no clear consensus regarding the generation time (GT) in sifakas. Sifakas first reproduce at ages as low as 3.5 years and have a lifespan of up to 32 years (recently reviewed in Zimmermann & Radespiel, 2015). Quéméré et al. (2012) used a value of 6 years in a previous study on P. tattersalli which corresponds to the median age at first reproduction for females in a related western sifaka species, P. verreauxi (Morris et al., 2011). However, they also considered other larger values including 17.5 years as estimated by Lawler et al. (2009) using demographic data from a long-term monitored P. verreauxi population. The true generation time is probably between these values. Here, to incorporate this source of uncertainty in the dating of demographic changes, we report estimates using generation time of 6 and 18 years that encompass the values used in previous work as well as the most recent GT estimations (Lawler et al., 2009; Morris et al., 2011) and the ones used by the IUCN.
2.3.2 Mutation rate
The estimation of the timing of demographic history events is affected by the assumed mutation rate. As there are no studies on microsatellites mutation rates in strepsirrhines, we used values between 10−4 and 10−3 that are widely assumed in demographic models (e.g., Goldstein, Linares, Cavalli-Sforza, & Feldman, 1995). This range is within the lower end of a range of pedigree-based estimates for autosomal microsatellites in humans (Ellegren, 2004) and within the higher range of humans–chimpanzees orthologous loci (Kelkar, Tyekucheva, Chiaromonte, & Makova, 2008). Our results—like most other demographic inferences using genetic data—depend on the assumption that the true mutation rate does not deviate dramatically from this commonly applied range.
The first method, developed by Storz and Beaumont (Storz & Beaumont, 2002) and implemented in msvar 1.3, assumes that the population underwent a single event of decline or growth and a strict stepwise-mutation model (SMM). It uses the information present in the full allelic distribution and a Bayesian coalescent-based MCMC approach to estimate the current and past effective population sizes N0 and N1 as well as the time T since the population change (in generations). We applied it to the Perrier's sifakas data so as to compare with the results of Quéméré et al. (2012). For Theta, we set a log mean of 3.5 with a standard deviation of 0.25, to favour mutation rate values between 10−4 and 10−3 (Storz & Beaumont, 2002). Wide “uninformative” priors and multiple runs with different starting points and different hyperprior parameters were used to avoid prior bias on posterior estimates (Table S1). At least four runs were performed for each sample with a total number of iterations always larger than 4 × 106 steps, discarding the first 10% to avoid influence in parameter estimation by starting conditions (burn-in period). The convergence of the four runs for each sample was checked both visually and using the Geweke convergence diagnostic (Geweke, 1992) implemented in the “coda” r (R Core Team 2014) package. The outputs of the runs were then merged to obtain robust estimates of the posterior distribution of the parameters.
The second method, implemented in the r package vareff, uses an approximate likelihood of the distribution of distance frequencies between alleles in a Monte Carlo Markov Chain framework (Nikolic & Chevalet, 2014). This approach offers several advantages over msvar: (i) it allows for several demographic changes, (ii) it implements the three most common microsatellites mutation models and (iii) it is much less computationally demanding than msvar 1.3. This gain in computational time enabled us to test several combinations of parameters such as the number of population size changes (JMAX parameter in vareff), mutation models (MODEL); the variance of the prior log-distribution of effective sizes (VARP1) and the maximal distance between alleles (DMAXPLUS). The final analyses were performed using each of the three mutation models (single step, two phase and geometric), a mutation rate of 5 × 10−3, allowing three population size changes (JMAX = 3). Additional parameters are detailed in supplementary material.
2.3.5 Potential effects of population structure on demographic inferences
msvar and VarrEff's models assume that samples are obtained from isolated populations, hence ignoring genetic substructure and migration from others populations. A growing number of studies showed that ignoring population substructure or using an inadequate sampling scheme may lead to spurious signatures of demographic change (Chikhi et al., 2010; Heller et al., 2013; Leblois, Estoup, & Streiff, 2006; Mazet, Rodriguez, Grusea, Boitard, & Chikhi, 2016; Mazet et al., 2015). We therefore performed several complementary analyses varying population sampling to test this potential bias. We considered either (i) pooled data (i.e., all samples from all possible subpopulations), (ii) population data (i.e., samples from each subpopulation or forest fragment analysed separately), (iii) random sampling (i.e., resampling a random subset of individuals from the data) and/or (iv) scattered data (one sample per sampling locality or forest fragment), when possible. Information on population genetic structure was taken from two previous studies that analysed spatial patterns of genetic differentiation (Quéméré et al., 2010b; Salmona et al., 2015; see supplementary material for further details).
2.3.6 Approximate Bayesian computation
To reconstruct the demographic history of Perrier's and Tattersall's sifakas, we also used abc approaches as implemented in abctoolbox version 1.1 (Wegmann, Leuenberger, Neuenschwander, & Excoffier, 2010) and in the r package “abc” (Csilléry et al., 2012). The principle behind abc is to compare data simulated under several alternatives scenarios to the real data, using (in general) summary statistics. Alternative scenarios can subsequently be compared and parameters of interest estimated from the most supported scenarios (Csilléry, Blum, Gaggiotti, & François, 2010). Within this framework, we simulated genetic data using the coalescent tool fastsimcoal version 1.1.2 (Excoffier & Foll, 2011). All simulations assumed a log-uniform mutation rate prior comprised between 10−3 and 10−4 as well as a set of specific parameters detailed in Tables S2 and S3 and explained below.
2.3.7 abc models
We first tested six simple scenarios assuming a single panmictic population, and differing from each other only by their history of population size change (Figure 2a). The first model (NULL, Figure 2a) posited constant population size and can be regarded as our null hypothesis model. The second model (1-SC) assumed a single population size change (expansion or decline). The third and fourth scenarios (1-BTL, 2-BTL) modelled, respectively, one and two bottlenecks without recovery. The fifth and sixth models (1-BTL − R, 2-BTL − R) are similar to the third and fourth, respectively, but incorporate population recovery after the first bottleneck (Figure 2a).
In addition, we also simulated two sets of scenarios assuming structured populations (Figure 2b) to (i) account for the P. tattersalli's population genetic structure revealed by Quéméré et al. (2010b, 2010c), (ii) test for potential effects of forest fragmentation on both sifaka species and finally (iii) consider confounding effects of structure and population size change (Chikhi et al., 2010; Heller et al., 2013; Peter et al., 2010). The first (see Figure 2b) and second sets of scenarios assuming structured populations (not shown in Figure 2b apart from model Str-NULL) differ in the state of the ancient population (panmictic vs. structured). The first structured scenario (Str-NULL) models demes with constant population size and connected by gene flow. Migration among demes is constant over time, mimicking long-term structure driven by natural barriers such as rivers, mountains or grasslands. The second scenario (Str-s-NULL) is similar to Str-NULL but aims to represent the creation of a barrier to gene flow without significant change of the population size (i.e., habitat fragmentation without loss scenario). This could represent (i) the disappearance of the main geneflow routes (opened habitat and corridors), or (ii) the enlargement of a river, while the core habitat (dense forest) is maintained. The last set of structured scenarios (Figure 2b) were based upon the Str-s-NULL and Str-NULL scenarios with the addition of one (Str-1-BTL) or two (Str-2-BTL) population decline(s) after fragmentation. The last scenario incorporates the idea that some subpopulations may have gone extinct in the past following a fragmentation event (Str-1-BTL + C). This scenario was built to test the potential extinction of various Perrier's and Tattersall's sifaka subpopulations outside their current distribution range (Dewar et al., 2013; Godfrey et al., 1996; Hawkins et al., 1990; Jungers et al., 1995 and references therein).
The results of the abc model choice procedure led us to identify scenarios that best explained the data for the two sifakas (see the model selection procedure below). For these scenarios, we tested two additional priors of timing (Ti) of decline (Tables S2 and S3), before and after the first documented human occurrence in Madagascar (T ~4000 years ago).
2.3.8 Summary statistics
We used arlsumstat version 188.8.131.52 (Excoffier & Lischer, 2010) to estimate summary statistics from the simulated and observed data sets (Table 1). These summary statistics were first chosen to capture as much information as possible from the microsatellites data sets about potential population size changes. Subsequently, the nine global summary statistics were filtered out on the basis of their linear relationship with parameters of interest and relative covariance; that is, when two statistics had a high covariance (>80%), we kept the one showing the best linear relationship with Ni, Ti, and Mi parameters. This was done as a preliminary filter to reduce the dimensionality of the summary statistics while retaining most of the information. Using this approach, the summary statistics sdK and sdR were therefore discarded for P. perrieri, while all statistics were kept for P. tattersalli. Structured and panmictic models were compared using the largest summary statistics set common to both data set (seven for P. perrieri and nine for P. tattersalli, Table 1). Tattersall sifaka structured models were compared using the largest possible number of summary statistics (up to 41) common to each pair of model. As the structured models had many more possible summary statistics, we furthermore performed a partial least square (PLS) regression to reduce the sets of summary statistics to a smaller number of independent components (see Supplementary material for further details).
Table 1. Summary statistics of the microsatellites data sets
All summary statistics were estimated with arlsumstat version 184.108.40.206 (Excoffier & Lischer, 2010).
Nind, number of diploid individuals; K, number of alleles; SD, standard deviation; H, observed heterozygosity; GW, Garza and Williamson (2001) index; R, allelic range.
Statistics calculated over the three clusters.
For P. perrieri, only the seven nonitalicized summary statistics were used for model choice and parameter inferences.
To assess model fit, we first calculated the marginal densities and the probability of the observed data with the generalized linear model (GLM) approach in abctoolbox (Leuenberger & Wegmann, 2010). The GLM was built from the 1,000 simulations (i.e., 0.2% of the 5 × 105) closest to the observed data. The p-value represents the proportion of the retained simulations showing a lower or equal likelihood under the inferred GLM as compared to the observed genetic data (Wegmann, Leuenberger, & Excoffier, 2009). Low p-values indicate that the observed data are unlikely to have been generated under the inferred GLM. Additionally, we estimated, from the model's marginal densities, the Bayes factor (BF), the ratio of the posterior densities of the two alternative hypotheses (i.e., scenario), over the ratio of the prior densities of the same alternative hypotheses. BF absolute values >3 were considered as significant evidence to reject the alternative hypothesis (Kass & Raftery, 1995). When two models showed BF absolute values <3, we kept the simplest model but considered both models‘ results for discussion. To confirm model choice, we also compared models with the “logistic” and “neuralnet” regression analysis and proportions of retained simulations ranging from 1 to 0.05% (i.e., 10,000–500 simulations/model) within the “abc” package in r (Csilléry et al., 2012).
2.3.10 abc validation and parameter estimation
Model selection and parameter estimation in an abc framework can suffer from the loss of information in the reduction in data to summary statistics (Csilléry et al., 2010). Therefore, we used a series of tests based on pseudo-observed data sets (pods) that allow to assess the accuracy of and validate the model selection and the parameter estimation procedures (see suppl. material for further details). To increase the accuracy of parameters, posterior distribution estimates from the best-fitting models (i.e., those with a significant Bayes factors as outlined above), we produced 2 × 106 simulated data sets under these models and estimated posteriors using the GLM approach with 0.1% (i.e., 2,000) simulations closest to the observed data based on exploratory analyses using a range from 1% to 0.05% simulations. To compare alternative temporally delineated hypothesis and identify the most likely time of demographic events within a Bayesian framework, we performed a BF analysis. We identified six time intervals corresponding to putative causes of sifaka historical demographic events in northern Madagascar. The BFs were computed for each of the six time intervals against all other periods taken together.
Similarly to the results of Quéméré et al. (2012) on P. tattersalli, we detected a clear bottleneck signal using the Storz and Beaumont (2002) method on P. perrieri. The posteriors of present and past effective population size log(N0) and log(N1) have distinct nonoverlapping distributions (Fig. S2) with respective median values of 2.24 (~170) and 4.43 (~27,000; Fig. S2a). All the posteriors were different from the priors and converged to similar distributions regardless of the priors used. The BF analysis favours a scenario with a ~21-fold population decrease. The posteriors of log(T), the time since population started to decrease, show a median around 4 (T = ~10,000 years BP; Table S4 and Figs S2b and S3) for a GT of 18 years (and T = ~3,300 years BP for a GT of 6 years) regardless of which prior distribution was used. The analyses carried out under various sampling schemes showed similar results (Table S5). Quéméré et al. (2012) similarly found a robust signal of P. tattersalli's population decline with medians posterior distribution of T ranging from ~7,000 yBP for a GT of 6 years to ~20,000 for GT values of 17 years. Altogether, msvar detected a bottleneck that started more than 3000 years ago in both species, but may have happened at different times for the two species.
Using the Nikolic and Chevalet (2014) method, we also detected a signal of a single bottleneck in both sifaka species (Figs S4 & S5), with log(N0) and log(N1) posteriors having distinct nonoverlapping distributions (Figs S4 & S5, bottom plots). While this approach allows for the detection of several changes in population size, we only detected one significant bottleneck (Figs S4 & S5). The method identifies a rather ancient population decline for Perrier's sifaka (~1,000 to 2,500 generations ago; ~6,000 to 45,000 yBP; lower range*low GT, larger range*large GT), and a more recent one for Tattersall's sifaka (~100 to 300 g. ago ~600 to 5,400 yBP; Figs S4 & S5 and Tables S6 and S7). The inferred decline is stronger for P. tattersalli than for P. perrieri, with log(N0) values around 3 for both, but log(N1) values, respectively, around 3.5–4 and 4–4.5, respectively (Figs S4 & S5 and Tables S6 and S7). The three mutation models tested show similar results with the exception of the timing of the demographic event for P. perrieri, which seems to be particularly influenced by the mutation model, hence leading to a large variance (between ~1,000 and 2,500 generation ago and ~6,000 to 45,000 yBP; Figs S4 and S5 and Tables S6 and S7). The analysis of random and/or scattered sampling schemes showed similar results, but wider range of posterior estimates for the time of decline (Figs S6–S8, Tables S6 and S7). The analysis of the three subpopulations of P. tattersalli identified by Quéméré et al. (2010b) shows variable results with the small population of the Antsaharaingy forest (North) showing very limited or no signal of population size change (Fig. S7). This result and the fact that Antsaharaingy forest (North) ancient population size shows no ovelap with those of the other populations suggest that Antsaharaingy forest may have been little connected with other forest for a long period of time (Fig. S7).
3.3.1 Model choice for P. perrieri
For P. perrieri, and when panmictic models were compared, we found strong support for a change in population size (1-SC vs. NULL, BF >10+105) with posterior estimates (low N0 and large N1) indicating a signal of population decrease (Fig. S9). The comparison with a scenario with one bottleneck (1-BTL) identified substantial support (1-BTL vs. 1-SC-Like, BF >19) and a good fit to the observed data (p-value = .99, Table S8). In addition, we tested the other scenarios including two sequential bottlenecks (2-BTL), with or without recovery after the first event (1-BTL+R, 2-BTL+R), and one bottleneck scenario with an ancient and recent time of bottleneck priors (1-BTL-O, 1-BTL-R, Figure 3a). All had lower support than 1-BTL under the GLM-BF analysis (Table S8; Figure 3a). When comparing models of population size decrease posterior or anterior to human arrival in Madagascar, we found slightly higher but not significant support for the recent bottleneck model (1-BTL-R vs. 1-BTL-O, BF = 1.95), with posterior distribution of the estimate of time T (for both models) skewed towards the values obtained in 1-BTL (Table S8; Fig. S10). In other words, the timing presented significant uncertainty.
We compared the HL&F or structured scenarios in a similar manner (Str-NULL, Str-1-BTL, etc. Figure 3b). Altogether, models showing fragmentation and population size changes showed greater support than models with constant-sized and structured population(s) (Str-s-1-BTL and Str-s-2-BTL vs. Str-s-NULL, Table S8; Figure 3b). The model with two successive bottlenecks showed slightly more support than the model with only one bottleneck (Str-s-2-BTL vs. Str-s-1-BTL, BF = 2.11) or the scenario modelling the loss of subpopulations (Str-s-2-BTL+C, Table S8; Figure 3b) but as above this was not significant. This suggests that models with population size changes and fragmentation are favoured, but that there is not enough information to infer the number of population size changes.
Models with a structured “ancestral” population showed better support than the equivalent scenarios where the ancestral population was panmictic (called hereafter “recently structured” models) or both ancestral and recent populations were panmictic (e.g., Str-NULL vs. Str-s-NULL, BF >10+9, Figure 3c). Similarly to previous comparisons, the scenario with structured but constant population size (Str-NULL) showed relatively low support when compared to models with one or several population size change (Str Constant vs. Non-constant, BF >1010). As there was little difference in support between the remaining structured models (Str-1-BTL, Str-2-BTL, Str-1-BTL+C, Table S8; Figure 3c), we selected the most parsimonious model with one bottleneck (Str-1-BTL) for further comparison. Finally, we compared the most supported model of each of the three categories, (i) panmictic, (ii) recently structured and (iii) structured. These last comparisons confirmed the support of the model Str-1-BTL (Figure 3d) for P. perrieri's population demographic history, a model with both recent and ancient structured population and with one population size decline.
3.3.2 Model choice for P. tattersalli
For P. tattersalli, we followed the same procedure to compare panmictic, recently structured and “ancestrally” structured models, with and without population size changes (Table S9; Fig. S11). Panmictic models (Fig. S11a) showed overall little support and a bad fit to the observed data (p-value <.01; Table S9; Figs S11a & S9b). As expected from the previous results of Quéméré et al. (2010b) who found that the population of P. tattersalli is genetically structured into subpopulations, the inclusion of population structure in our model simulations greatly increased their fit to the data (Table S9) with BF >1010 for all pair comparisons of similar structured and panmictic models. Models with an ancient panmictic population (partially structured models) showed overall good support for models with one or two population decline(s) (and with recent change in gene flow (Str-s-1BTL+GF; Table S9; Fig. S11b). The models with ancient structured population showed very similar model support and model choice than recently structured models (Table S9; Fig. S11c). Finally, from comparisons of the best models from the three sets (Fig. S11d), we found that there was no significant difference between the Str-s-1-BTL and the Str-1-BTL models (BF = 1.04) and kept both models for further parameter estimations.
3.3.3 Parameter estimation for P. perrieri and P. tattersalli
Under the abc framework, we estimated the parameters (Fig. S12) of the best models Str-1-BTL (M-35 in Tables 2 and S8) and Str-s-1-BTL (M-29) for P. perrieri as well as models Str-1-BTL (M-25 in Tables 2 and S8) and Str-s-1-BTL (M-22) for P. tattersalli (Table 3; Figure 4). Under these structured models, both species show small current and relatively large past deme size. Perrier's sifaka showed posterior values of current total population size of N ≈ 800, where N is the sum of all deme sizes. This total size is slightly smaller than for Tattersall's sifaka (N ≈1,250; Figure 4a,d), but both species showed large ancient population size values (N >30,000). Most estimates of the timing of the deme size change showed values coinciding with and/or following the most ancient evidence of human presence on the island as well as the start of a long-lasting drought period at the mid-Holocene, ~4 to 5,000 yBP (Table 2; Figures 4b,e and 6). The BF analysis (Figure 5) globally favoured the third scenario (decline occurring between 1,000 and 4,500 yBP) for both species. Tattersall's sifaka showed a more recent population decline than Perrier's sifaka (Tables 2 and 3; Figures 4b,e, 5 and 6) with mode values posterior to the first known dates of human traces in the region ~4 to 5,000 yBP (Tables 2 and 3; Figures 4b,e, 5 and 6). For models with two successive bottlenecks and models with the extinction of subpopulations (Figure 3), posterior estimates for the time of the second events, second bottleneck or reduced gene flow (T1; Figure 4), showed values in the recent past (T1 ≈ 300–400 y BP). Migration rate estimates were relatively high for P. perrieri with ~6.5 migrants between each pair of demes per generation, but they also exhibited a large variance. Tattersall's sifaka showed comparatively lower and narrower migration rates estimates with most posterior modes ~1 (Table 2; Figure 4c,f).
Table 2. Posterior parameter estimates from best ABC demographic models for P. perrieri
Nat (y.) GT18
Nat (y.) GT6
Values of median, mean, mode and highest posterior density (HPD) intervals of the five posterior distribution N0 (sum of current deme sizes), N1 (sum of ancient deme sizes before the demographic event), T time in logarithm 10 (Log10) of the number of generation, in generations (Nat(gen)), and in years (Nat(y.)) elapsed since demographic event for a generation time of 18 years (see main text for discussion) and of the number of migrant between demes per generation NM1 before and after the demographic event (NM) obtained from the best models selected with the ABC framework. Model numbers refer to Tables S2 and S8, as well as Figure 3. The scenario M-29 (Str-s-1-BTL) and M-35 (Str-1-BTL) both model a population size decrease from a panmictic (M-23) or structured ancient population (M-25) into a structured population since the decrease.
Table 3. Posterior parameters estimates from best ABC demographic models for P. tattersalli
Nat (y.) GT18
Nat (y.) GT6
Values of median, mean, mode and highest posterior density (HPD) intervals of the five posterior distribution N0 (sum of current deme sizes), N1 (sum of ancient deme sizes before the demographic event), T time in logarithm 10 (Log10) of the number of generation, in generations (Nat(gen)), and in years (Nat(y.)) elapsed since demographic event for a generation time of 18 years (see main text for discussion) and of the number of migrant between demes per generation NM1 before and after the demographic event (NM) obtained from the best models selected with the ABC framework. Model numbers refer to Tables S3 and S9, as well as Fig. S11. The scenarios M-23 (Str-s-1-BTL) and M-25 (Str-1-BTL) both model a population size decrease from a panmictic (M-23) or structured ancient population (M-25) into a structured population since the decrease.
The accuracy indicators for parameter moment estimates of P. perrieri (Csilléry et al., 2012; Leuenberger & Wegmann, 2010) showed that N0 and T2 are estimated with relatively good accuracy in contrast to N1 and NM which showed lower estimation accuracy (Fig. S13). For P. tattersalli, all parameter estimates showed high accuracy (i.e., low index values, Figure S14).
4.1 Inferring the past demography of sifakas and other lemurs in a complex historical context
With the exception of Olivieri, Sousa, Chikhi, and Radespiel (2008), who studied the demographic history of three mouse lemur species, previous demographic inferences in lemurs and other Malagasy species investigated only a single species and a restricted number of sampling locations (Lawler, 2008, 2011; Louis et al., 2005; Markolf, Roos, & Kappeler, 2008). Moreover, in most cases, population structure is ignored (Lawler, 2011; Markolf et al., 2008; Meyer et al., 2015) and even when it is accounted for, only a single inferential approach was used (Quéméré et al., 2012). Here, we used several methods and considered population structure and population size changes together in a combined framework. The abc models used were based on previously published population genetic studies on northern sifakas and samples from the entire species’ ranges (Bailey et al., 2016; Quéméré et al., 2010b, 2010c, 2012; Salmona et al., 2015). To provide a thorough and comparable set of results across the species, we analysed the data from the two species using the same methods. Despite some discrepancies, the three methods (msvar, vareff, and abc) showed a consistent and coherent bottleneck signal in both species data set when the abc model did not include population structure (and was hence comparable to the other methods).
4.2 On the importance of accounting for population structure in HL&F scenarios
However, several studies have shown that population structure and changes in connectivity can generate signals that methods such as msvar and vareff will interpret as population size changes. For example, the current population size inferred using msvar and vareff (which both assume panmixia) can be thought of as a “local” deme size, whereas the “ancient” population size pertains to a “meta” (unfragmented) population size. Further details about this reasoning can be found in Wakeley (1999), Nielsen and Beaumont (2004), Chikhi et al. (2010), Heller et al. (2013) and Mazet et al. (2015, 2016). The best abc models suggest that both species are likely to have been structured in the past and that this structure has shifted (forests appear to be more fragmented today than they were in the past) with a significant change in the size of the demes, and thus a change of the total population size. These “bottlenecks” cannot be dated precisely and may have happened at different times between the two species, likely around ~300 to 1,500 yBP for P. tattersalli and ~1,000 to 3,400 yBP for P. perrieri (Tables 2 and 3; Figures 4–6). For P. tattersalli, these dates are significantly more recent than those obtained by Quéméré et al. (2012) using msvar (>4,000 yBP). P. perrieri has not been investigated before, but the msvar dating obtained here was again older than under the abc models with structure. However, it is important to stress that the different methods are difficult to compare because the priors and mutation models are not identical. The comparison between the bottleneck dating in structured and nonstructured abc is more illustrative regarding this point. Here, we found that for P. perrieri the posterior is extremely wide and nearly identical to the prior (Table S11). For P. tattersalli, the model of population size change, in agreement with msvar, suggests an ancient event (between ~5,000 and ~15,000; Table S12), confirming that neglecting structure may in some cases bias the dating of an inferred bottleneck. Altogether, this works illustrates that caution should be taken when inferring and dating events, as model assumptions may lead to rather different conclusions. While we stress the importance of accounting for the structure, we should not discard methods such as msvar or vareff. We showed here that these two approaches can be used to detect the effect of structure by varying the sampling strategy as previously noted by several authors (Chikhi et al., 2010; Heller et al., 2013; Wakeley, 1999). For example, we observed that using different sampling strategy in vareff led to different estimates for the same ancient population size (e.g., Figs S5 & S8). This kind of behaviour may be caused by population structure (Nikolic & Chevalet, 2014).
4.3 On the respective role of past climate changes and human activities on sifaka populations
Our analyses allowed us to decipher different aspects of the demographic history of two sifaka species and shed light on past vegetation changes in northern Madagascar. Altogether, our results suggest that major changes in connectivity and population size may have happened more recently than previously believed. Despite the uncertainties in our bottleneck datings, we stress that they overlap (with or without structure) with the first documented human presence in the region (Dewar et al., 2013), but also with a worldwide major drought event termed the mid-Holocene boundary (Formal subdivision of the Holocene Series/Epoch: A Discussion Paper by a Working Group of INTIMATE, 2012; Figures 5 and 6). Given that both species are protected by a strong local taboo and are not hunted by local Sakalava population (Anania et al., 2017; Banks, Ellis, & Wright, 2007; Quéméré et al., 2010b), we can reasonably assume that the decline of sifaka populations in northern Madagascar is not due to a heavy direct human pressure. The abc analyses suggest that the two species have likely been structured prior to human arrival and that an increase in the fragmentation accompanied by a population size reduction probably occurred concomitantly with the first human settlements in the region.
Our data also suggest that at least P. perrieri's population was not able to recover after a period that coincides with the decline of a surprisingly large number of species in the western Indian Ocean (Table S10 and paragraphs below). Today, P. perrieri has a very small and fragmented population with very few individuals left in the wild (Banks, et al. In press, Banks et al., 2015) while our results suggest that 2–5000 years ago the species had a much larger population size and a wider distribution. Its habitat has been severely degraded in the last decades, and several studies suggest that P. perrieri has disappeared from several forests in which it was found within the last millenia. The actions of humans have likely played a critical role during this period. However, it may be worth noting that P. perrieri's decline appears more ancient than that of P. tattersalli under all the models and may therefore have been more influenced by climate change (Holocene droughts, Figures 5 and 6, see Supplementary material “Climate change “).
P. tattersalli's present-day situation is not as concerning as the total census population size is at least an order of magnitude higher (likely >18,000 individuals, Quéméré et al., 2010a). Its population is also divided into small forest fragments, but most of these fragments remain connected by a network of riparian corridors. Our results suggest that P. tattersalli's population decline was more recent than that of P. perrieri and that estimated by Quéméré et al. (2012). This collapse likely results from the combined effect of severe droughts occurring in the second half of the Holocene and human-induced changes in forest habitats (Figures 5 and 6, see Supplementary material “Human presence and impact” and “fire dynamics”).
4.4 Multidisciplinary perspectives
Drought events such as the mid-Holocene boundary likely led to forest fragmentation. Subsequently, the persistent aridity during the second half of the Holocene (Burney, 1993; Kiage & Liu, 2006), as well as more recent human-driven landscape modifications (at least since ~1,000 yBP on, Burns et al., 2016), has maintained open habitats impeding forest re-expansion. This scenario provides a coherent—if tentative—explanation to the range contraction and peculiar disjunct distribution of the two study species in the north of Madagascar, while the genus shows parapatric continuums in eastern and western regions (Supplementary material “Lemurs paleodistribution in northern Madagascar” and Fig. S1). It is also coherent with subfossil records (Figure 6 and Supplementary material “Lemurs paleodistribution in northern Madagascar”) and makes sense with regard to the increasing number of species showing approximately congruent genetic signals of demographic change in Madagascar and in the western Indian Ocean (Table S10).
MacPhee, Burney, and Wells (1985) and Matsumoto and Burney (1994) reported open grasslands vegetation in the dry western region well before any evidence of human settlements in that region. Furthermore, several studies documented the relative antiquity of grassland communities based on the large diversity of C4 grass lineages and the presence of plant and animal species endemic to Malagasy grassy biomes (Besnard et al., 2014a; Bond, Silander, Ranaivonasy, & Ratsirarson, 2008; Vorontsova et al., 2016; Willis, Gillson, & Virah-Sawmy, 2008). This evidence together with our abc results supports a scenario in which open habitat was present well before human arrival in northern Madagascar in agreement with the work of Quéméré et al. (2012). The relative contributions of human (i.e., “overkill scenario”) and natural environmental factors in the quaternary megafauna extinctions appear to vary strongly among regions of the world (Cooper et al., 2015; Grund, Surovell, & Lyons, 2012; Muldoon et al., 2012; Stuart, 2015) in relation to the rate of climatic changes (Nogués-Bravo, Ohlemüller, Batra, & Araújo, 2010; Prescott, Williams, Balmford, Green, & Manica, 2012) and the pattern of human settlement and activities (Koch & Barnosky, 2006; Stuart, 2015). In Madagascar, the situation is particularly interesting since humans were thought to have arrived only 2500 years ago, due to the lack of more ancient traces. This suggests that humans may have been present in Madagascar with limited impact on the endemic fauna or flora for a long period. The most notable environmental changes seem to have started during the shift towards agro-pastoralism around 1000 years ago. This would be >3,000 years after the first settlements (Figure 6). This agro-pastoralism shift represented a major change in fire regimes in Madagascar (Supplementary material “Fire dynamics”), but not necessarily (or not yet documented) in the North.
Our results also call for multidisciplinary collaborative research. Genetic data alone do not allow identifying the causes of the inferred histories. Sediment (e.g., Burney et al., 2003) and palynological cores (e.g., Virah-Sawmy et al., 2010) coupled with environmental ancient DNA analysis from sediment (sedaDNA; e.g., Parducci et al., 2017; Pansu et al., 2015) would be of great complementarity to unravel the landscape history of this region. To date, carbon dating of subfossils focuses mainly on extinct Madagascar fauna even though a substantial quantity of extant species subfossils are usually recovered in archaeological sites (Goodman & Jungers, 2014). Extending carbon dating to a larger range of taxa would provide a better understanding of Madagascar's history, clarify the palaeodistribution of extant species (e.g., Kistler et al., 2015) and may reveal unexpected patterns.
4.5 Genomic perspectives
The increasing availability of genomic data (including lemurs: mouse lemur, aye-aye and eulemur—Perry et al., 2012, 2013; Meyer et al., 2015; P. coquereli’s—GCA_000956105.1; and P. tattersalli—SRX701290-93) coupled with the recent development of methods to infer demographic history from genomic data (Li & Durbin, 2011; Liu & Fu, 2015; Mazet et al., 2015; Salmona, Heller, Lascoux & Shafer, 2017) promise to resolve demographic history of species with unprecedented resolution. Although their application is challenging for noninvasive DNA (that include exogenous plant, fungi and bacteria DNA; Lynn, Sechi, Chikhi, & Goossens, 2016), capture approaches (e.g., Fabbri et al., 2012; Perry, Marioni, Melsted, & Gilad, 2010) associated with reduction in genome complexity (e.g., RAD; Suchan et al., 2016; Ali et al., 2016; Hoffberg et al., 2016) promise to alleviate these challenges. The analysis of thousands of loci spread across the whole genome with abc - (Nater et al., 2015; Shafer, Gattepaille, Stewart, & Wolf, 2015a) or SFS-based approaches allow for complex model inferences (Excoffier, Dupanloup, Huerta-Sánchez, Sousa, & Foll, 2013; Gutenkunst et al., 2009) and may soon enable us to obtain a refined picture of the islands past by detecting older demographic events, and clarifying their timing. Finally, comparison with closely related species (i.e., Sgarlata et al., 2016) and sympatric forest dwelling species (Aleixo-Pais et al., 2017; Sgarlata et al., 2017), with contrasting landscape use (Knoop, Chikhi, & Salmona, 2017), may help to confirm the patterns identified for P. perrieri and P. tattersalli. In particular, the mouse lemur and sportive lemur species inhabiting the region (Microcebus tavaratra and Lepilemur milanoii and L. ankaranensis) have shorter generation time and their study may reveal northern Madagascar recent history (the past 5,000 years) in greater detail (Salmona 2015, Yoder et al., 2016).
We thank CAFF/CORE, the “Direction Générale de l'Environnement et des Forêts,” Madagascar National Park, the Fanamby NGO (including S. Rajaobelina, V. Rasoloarison, P. Ranarison and S. Wohlhauser), the “Direction Régionale de l'Environnement et des Forêts région DIANA, M. Banks for discussion and advice on Perrier's sifaka. The fieldwork was possible thanks to the continuous support of the “Département de Biologie Animale et Ecologie,” University of Mahajanga, the University of Antsiranana and to a large extent, thanks to the participation the Malagasy master students the field assistants and volunteers, of many great local guides and cooks which we warmly thanks for their help in the field and for sharing their incomparable expertise of the forest, misaotra anaero jiaby. Finally, we would like to thanks the Bioinformatic and Genomics Unit at the IGC, Sequencing Service for their collaboration and U. Radespiel, L. Wilmé, E Van-Campo and G. Besnard for comments and discussion on the manuscript. Financial support for this study was provided by the “Fundação para a Ciência e a Tecnologia” (grant number SFRH/BD/64875/2009 to J.S. and grant numbers Biodiversa/0003/2015, PTDC/BIA-BIC/4476/2012, PTDC/BIA-BEC/100176/2008 to L.C.), the GDRI Madagascar, the “Laboratoire d'Excellence” (LABEX) entitled TULIP (ANR-10-LABX-41), “Rufford Small Grant Foundation” (grant number 10941-1 to J.S.), the Instituto Gulbenkian de Ciência, the LIA BEEG-B (Laboratoire International Associé—Bioinformatics, Ecology, Evolution, Genomics and Behaviour) (CNRS), the European Science Foundation, ConGenomics Research networking programme, Grant Number 5094 to JS. Rasmus Heller was funded by research grants from the Danish National Research Council (DFF) and the Villum Foundation Young Investigator programme. This study was conducted in agreement with the laws of the countries of Portugal, France and Madagascar.
J.S., E.Q. and L.C. produced data, J.S., R.H., E.Q. & L.C. designed and discussed the experiments, J.S. performed the analyses, J.S., E.Q., R.H. and L.C. discussed the results, wrote and critically revised the paper.