Past environmental changes affected lemur population dynamics prior to human impact in Madagascar

Quaternary climatic changes have been invoked as important drivers of species diversification worldwide. However, the impact of such changes on vegetation and animal population dynamics in tropical regions remains debated. To overcome this uncertainty, we integrated high-resolution paleoenvironmental reconstructions from a sedimentary record covering the past 25,000 years with demographic inferences of a forest-dwelling primate species (Microcebus arnholdi), in northern Madagascar. Result comparisons suggest that climate changes through the African Humid Period (15.2 – 5.5 kyr) strongly affected the demographic dynamics of M. arnholdi. We further inferred a population decline in the last millennium which was likely shaped by the combination of climatic and anthropogenic impacts. Our findings demonstrate that population fluctuations in Malagasy wildlife were substantial prior to a significant human impact. This provides a critical knowledge of climatically driven, environmental and ecological changes in the past, which is essential to better understand the dynamics and resilience of current biodiversity.

P ast climate changes have often been invoked as important drivers of population demographic dynamics, species distribution, and genomic diversity patterns, as species directly respond to environmental changes and to changes in resource availability 1 . Mainly studied in temperate regions [2][3][4][5] , such dynamics remain much less explored and more debated in the tropics 6,7 , which are characterized by a very high biodiversity richness. Among tropical regions, Madagascar, with its long history of isolation, high rates of endemism 8,9 , and no significant human impact until the late Holocene 10 , is a key model area to understand how natural climate changes impacted evolutionary trajectories. Past climate fluctuations have been proposed as one of the main drivers of species diversification on the island, for example in lemurs, amphibians, or reptiles [11][12][13][14] . Mouse lemurs (Microcebus spp.) are the world's smallest extant primates and form a speciose genus among lemurs with 25 currently recognized species 15 . They represent a suitable model system to investigate the impact of past environmental changes on demographic dynamics because (i) they are forest-dependent and should therefore be strongly influenced by vegetation changes; (ii) they have a comparably short generation time and a high reproductive rate, thereby quickly accumulating genetic signatures of past demographic events; and (iii) populations are still large enough and may have preserved enough genetic diversity to accurately reconstruct their demographic history 11,16,17 . A wide range of molecular approaches have been recently developed to infer species demographic history using genomic data [18][19][20][21][22][23] . However, interpretations are challenging, given the overall lack of historic population data and the scarcity of paleoenvironmental reconstructions. In particular, the evaluation of realistic demographic scenarios (e.g., combining population size and population structure) requires information about past habitat fluctuations that can only be provided by local paleoenvironmental reconstructions.
To investigate how past environmental changes shaped current biodiversity in a tropical rainforest ecosystem, we combined highresolution paleoenvironmental reconstructions and demographic modeling of mouse lemur genomics data generated for the same study area, Montagne d'Ambre National Park (NP, Fig. 1a). Located in northern Madagascar, this NP fulfills two essential prerequisites: (i) mouse lemurs are known to occur in the evergreen humid forest of this NP 24 and (ii) this NP hosts volcanic (maar) lakes. Such maar lakes are well known as ideally suited archives as they continuously accumulate sediments, which serve as natural recorders of local and regional environmental and climate dynamics. Since the maar lakes in Montagne d'Ambre are at least of Neogene or Quaternary age, their deposits are assumed to reflect changes in the climatic and environmental conditions during critical periods such as the Pleistocene-Holocene transition. This transition led to abrupt climate and vegetation changes in mainland Africa through the so-called African Humid Period (AHP) 25,26 . Using genome-wide data, four complementary demographic approaches (Stairway Plot, PSMC, IICR simulations, and fastsimcoal2) were used to infer past population size and connectivity changes of Microcebus arnholdi in two parts of the NP, about 18 km from each other. Specifically, mouse lemur dynamics were studied in the northern part of the NP inside the evergreen humid forest (Mahasarika, 1073 m asl, Fig. 1a) and in the southern part of the NP along the forest-grassland ecotone (Fantany, 848 m asl, Fig. 1a). The inferred demographic dynamics were then compared and integrated with paleoenvironmental reconstructions that were derived from a continuous sediment record from Lac Maudit, located at the center of the NP in the evergreen humid forest (1250 m asl Fig. 1a) and less than 6 km away from either mouse lemur site. Such a dual paleoecological-genomic approach that compares genomic inferences with predictions derived from local paleoenvironmental datasets is unique for tropical regions. Our result comparisons demonstrate that the demographic dynamics of M. arnholdi were exclusively shaped by natural climate changes until the late Holocene, whereas the dynamics during the last millennium were likely shaped by both natural climate and anthropogenic changes. Altogether, our study illuminates how natural climate fluctuations and anthropogenic activities impacted Malagasy wildlife across time and provides a critical understanding of biological responses of forest-dwelling species to ongoing and accelerating climate change.

Results and discussion
Five periods of distinct environmental conditions during the past 25 kyr. A robust chronology was established for the sedimentary record of Lac Maudit by using a set of 19 14 C ages with an age-modeling approach (Supplementary section 1.2). For the environmental reconstruction, a multiproxy approach was applied, which combined palynological (and charcoal counting), granulometric and inorganic geochemical (X-ray fluorescence core scan, XRF) analyses (see Paleoenvironmental analyses in Methods). Based on pollen zones, granulometric properties as well as a Principal Component Analysis (PCA) of the XRF data, our results unveil five distinct periods of different climatic conditions associated with different plant assemblages ( Fig. 2 Supplementary  Table S2. b Assignment of 38 individuals to two genetic clusters (K = 2) using genotype likelihoods. Each individual is represented by a single horizontal bar and each color represents a distinct genetic cluster. Samples are sorted by sampling site. The animals sampled in the northern and southern site were assigned to two distinct genetic clusters with low levels of admixture between the two sites.
precipitation decline after the AHP, and a fifth period (>0.9 kyr) related to the increase of fire activity.
Prior to 15.2 kyr (period 1; 25-15.2 kyr), pollen assemblages show that montane vegetation (25-70%; Fig. 2c) and evergreen humid forest (10-40%; Fig. 2d) were growing together in Montagne d'Ambre. Today, the NP (>800 m asl) harbors only evergreen humid forest up to the mountaintop (1475 m asl). The montane vegetation that grows elsewhere in northern Madagascar above 1800 m asl under lower temperature and precipitation is absent in Montagne d'Ambre ( Fig. 1a and see Environmental settings in Methods). Therefore, the occurrence of montane vegetation in Montagne d'Ambre prior to 15.2 kyr suggests drier and cooler climate conditions than in modern times, similar to what was described for tropical East Africa 27,28 . Dry conditions at the study site are also evidenced by low sediment-accumulation rates of ca. 0.07 mm a −1 revealed by the age model (Fig. S2). Furthermore, texture and grain size show a very compact sediment with silty-to-clayey deposits (mean grain size of 7.5 µm) related to reduced precipitation and runoff amounts 29 ( Fig. 2b and Supplementary section 1.1). The PCA carried out with XRF data resulted in three main axes, with Axis 1 representing 61.5%, Axis 2 24.4%, and Axis 3 9.4% of the total data variance (Supplementary section 1.3). The first axis shows negative loadings for all measured elements, which is interpreted as high minerogenic input and low organic matter content consistent with dry conditions during the first period (Fig. S3). For the second axis, positive (negative) values mainly reflect increase (decrease) of iron (Fe), titanium (Ti), and nickel (Ni) that correspond to an enhanced (reduced) chemical weathering related to more (less) humid conditions and associated forest soil formation ( Fig. 2a and Supplementary section 1.3). Increased siliciclastic input (high loadings for silicon (Si), calcium (Ca), and strontium (Sr); Fig. S3) into the lake prior to 15.2 ka, reflected by low values in Axis 2 ( Fig. 2a), points to dry environmental conditions favoring physical weathering of bare, exposed materials (bare rock) 30 . Due to the dry environmental conditions, a low soil-formation rate and development of montane vegetation was favored, associated with enhanced physical weathering. During period 1, we therefore conclude that montane vegetation was present in the upper parts of Montagne d'Ambre at lower altitudes than is observed today in other regions of northern Madagascar. The evergreen humid forest also shifted to lower elevations in comparison with modern times and/or was restricted to sheltered habitats with humid local conditions in Montagne d'Ambre (e.g., valleys, rivers). However, a wide lowland expansion of this forest would be unexpected due to the rather dry conditions recorded between 25 and 15 kyr in tropical East Africa 26 , which were accompanied by a savanna expansion in African lowlands 31 .
During the second period (15.2-11.8 kyr), an increase in sediment accumulation rate by one order of magnitude (mean ca. 0.6 mm a −1 ; Fig. S2) along with larger, but variable mean particle sizes (6-23 µm; Fig. 2b) is pointing to higher transport energies and generally more erosive processes in the catchment. A sudden sediment color change from gray to dark-brown likely indicates higher amounts of organic matter and thus a stronger bioproductivity (Supplementary section 1.1 and Fig. S1). Increasing values in Axis 2 of the PCA show a trend toward higher portions of Fe, Ti, and Ni within the deposits ( Fig. 2a and Supplementary section 1.3), which is interpreted as enhanced chemical weathering and deep soil formation. Contemporaneously, a rapid development of evergreen humid forest taxa (>50%; Fig. 2d) at the expense of montane vegetation (<25%; Fig. 2c) is observed. From these proxies, we infer a spectacular increase in moisture availability and thus enhanced precipitation in Montagne d' Ambre. The major increase of evergreen humid forest also indicates an upward shift of the altitudinal vegetation zones and suggests a temperature increase in addition to precipitation increase. During this period, a short increase of fire events (14.5-13.5 kyr; Fig. 2f) is indicated by charcoal influx. However, the low values of charcoal influx (less than three particles per cubic-centimeter sample) do not support local fire occurrences at the catchment of the lake. Instead, this may rather suggest a development of vegetation prone to fire in the region,  supported by an increase of Myrica pollen (Fig. S4). A likely increase of thunderstorms with lightning frequency and the development of suitable combustible biomass may best explain the detected episodes of wildland fires. Starting at 15.2 kyr, period 2 matches perfectly with the onset of the AHP documented in Madagascar's central highlands 32 and in East Africa 26 . This shows that climate changes were synchronous across this entire region. Maximum proportions of evergreen humid forest (>60%; Fig. 2d) start at 11.8 kyr and indicate a third period (11.8-5.5 kyr) at the transition to the Holocene. From 10 to 5.5 kyr, even enhanced chemical weathering is evidenced by maximum values in Axis 2 of the PCA performed on XRF data (Fig. 2a) and rather fine-grained, dominantly silty and clayey deposits, comparable with those from period one (mean grain size ca. 7 µm), are observed (Fig. 2b). This likely points to maximum rates of precipitation within the entire record and to a quite high lake level. The latter is further supported by low abundance of aquatic plants (<5%; Fig. S4), which are generally related to shallow lake conditions. We conclude that during this period, evergreen humid forest grew on the entire mountaintop as it is observed today. In addition, the maximum rates of precipitation recorded during this period imply a lower altitudinal limit and thereby wider distribution of this forest type than today.
During the fourth period (5.5-0.9 kyr), a coarsening in the particle size (>10 µm; Fig. 2b) as well as higher abundance of aquatic plants (10-30%; Fig. S4) suggest a lake level drop, likely indicating generally drier conditions associated with reduced precipitation. A contemporaneous weakening in chemical weathering derived from lower values in Axis 2 (PCA of XRF data; Fig. 2a) is also consistent with a reduced moisture availability. These results denote the termination of the AHP and match with paleoecological data from East Africa where the termination of the AHP starts around 5 kyr 26 . In Montagne d'Ambre, this precipitation decrease may have led to a contraction of the evergreen humid forest by moving upward its lower altitudinal limit as also observed in East Africa 33 . However, with its central and higher location in Montagne d'Ambre, the lake site was probably not able to record this vegetation shift. This certainly explains why, except for a few changes in forest composition after 5.5 kyr, the evergreen humid forest surrounding the lake has not been markedly impacted by this reduction of humidity ( Fig. 2d and Supplementary section 1.4).
The last~0.9 kyr, corresponding to the fifth period, are characterized by reduced evergreen humid forest taxa (ca. 40%; Fig. 2d), although still present at Lac Maudit. This period is marked by a tremendous increase of charcoal particles reflecting an important increase in fire frequency (Fig. 2f). Combined with the high amounts of Poaceae pollen (>40%; Fig. 2e While human arrival during the early or mid-Holocene is still hotly debated for Madagascar 34,35 , it has been shown that humans impacted northwestern Madagascar, only a couple of hundred kilometers south of Montagne d'Ambre, significantly through a rapid transition to a grass-dominated landscape with an abrupt increase in fire activity at~1.2 kyr 36,37 . Furthermore, large settlements developed from 1 kyr along the coast such as at Lakaton'i Anja located next to Montagne d'Ambre 38,39 . In congruence with these findings, our results suggest an important human impact on the northern landscapes at least over the last millennium. Late quaternary lemur population dynamics. The four complementary demographic approaches implemented in this study confirm that M. arnholdi underwent substantial demographic changes during the past 25 kyr. The Stairway Plot 21 and the Pairwise Sequentially Markovian Coalescent (PSMC) 19 methods were consistent in inferring the occurrence of an apparent population decrease in recent times (periods 5 and 4, <5.5 kyr) in both Mahasarika and Fantany (Fig. 3a,c). This population decrease was preceded by a period of a large population size, itself preceded by a more ancient population bottleneck. While the two methods identified a similar overall trend, the temporal dynamics were different among them. The ancient population bottleneck inferred by the Stairway Plot (Fig. 3c) dated >25 kyr, whereas the PSMC method (Fig. 3a) suggested lowest population sizes at the beginning of period 1 (25-15.2 kyr). Although the methods differ in the duration of the subsequent population recovery, both methods agreed that mouse lemurs had a relatively large population size during the maximum humid conditions of the AHP (i.e., period 3, 11.8-5.5 kyr) at both sampling sites. Differences on the temporal dynamics among the two methods are not surprising, as it has been shown that PSMC performs best for ancient times, while the Stairway Plot tends to perform best toward the present [40][41][42][43] . However, the interpretations of the previous demographic approaches rely on the assumption that population structure can be neglected and that mating is random. In an attempt to overcome this limitation of panmictic-based interpretations of the PSMC 44-46 , we simulated the Inverse Instantaneous Coalescent Rate (IICR) 47 of the southern site (Fig. 3b) by changing connectivity parameters, but assuming a constant population size. The model parameters that generated the best fit between the IICR in such a structured model and the dynamics inferred by the PSMC method suggest an n-island model of migration (29 demes interconnected by symmetrical gene flow) with four changes in connectivity. Those changes in connectivity dated around 65, 40, 10-15 and 1-2 kyr (vertical lines in Fig. 3b). More specifically, the IICR simulations suggest a period of high connectivity that started long before period 1 (~40 kyr), which was followed by a period of reduced connectivity during the AHP (~13 kyr) (see Supplementary section 2.1.9 for details about connectivity inferences). Such a scenario is unlikely, since the paleoenvironmental records revealed relatively unfavorable conditions (=lower extension of evergreen humid forest) during period 1 (25-15.2 kyr) and increasingly favorable conditions during the AHP (= maximum extension of evergreen humid forest). We therefore conclude that population size changes explain the PSMC dynamics during this time period better than changes in population connectivity (but assuming no change in population size).
To disentangle the confounding effects of population-size changes and structure, and allow for both factors to potentially explain our genomic data, we built and compared the congruence of 13 alternative demographic models of increasing complexity (M1-M13) to the RADseq data using a composite likelihood framework implemented in fastsimcoal2 20 . The models incorporated various combinations of changes in population size, population structure (with models involving population splits), and changes in connectivity among the sampling sites (Supplementary Figs. S8 and S9). The likelihood comparison of all 13 demographic models revealed that models including population structure had a better fit than those assuming panmixia, but also that recently structured models (i.e., models where a panmictic population become recently structured) with population size changes showed the highest congruence to the observed data (Supplementary results 2.2.3 and Table S6). The best fitting model (M7) suggested the existence of a large ancient population that only became structured (i.e., split in several demes) and suffered a first reduction on population size around 5 kyr. This event was followed by a more recent population decline at~1 kyr (Fig. 3d, Table 1). However, the high estimated levels of gene flow among the populations after they became structured (2NM0~7.23, where 2NM0 = average number of haploid immigrants entering the population per generation), suggest that the northern and southern site were probably never completely isolated from each other. This last approach not only allowed us to identify the most likely scenario incorporating both changes in population size and structure, but also demonstrates that both types of events are necessary to explain the data.
Late quaternary paleoenvironmental shifts and first human impact on endemic fauna. Population genomics and paleoenvironmental reconstructions were generated from the same study area to investigate how past climatic changes impacted the population dynamics of a forest-dwelling species. One of the most striking results is the significant impact of the AHP on the paleoecology and on the demographic history of M. arnholdi. In particular, the PSMC inferred a population decline preceding the AHP (contemporaneous to the Last Glacial Maximum; period 1; 25-15.2 kyr), after which the abrupt precipitation increase of moisture conditions related to the AHP (period 2 and 3; it predicted a lower connectivity during the AHP, which is unlikely, given the large expansion of evergreen humid forest during this time period. In particular, during period 3 (11.8-5.5 kyr), climate conditions supported an even wider extent of evergreen humid forest than in present times, which likely allowed mouse lemur populations to expand and grow further. Although we cannot determine the exact distribution limits of this forest during the AHP, specifically in lowland areas, the current occurrence of M. arnholdi in isolated humid-habitat patches outside the present limits of Montagne d'Ambre 24 suggests a temporarily broad expansion of evergreen humid forest in the past. Connectivity between lemur populations should have been high during such periods of increased humidity, as it could have been the case during the AHP. This hypothesis is further supported by the discovery of subfossil lemurs, currently restricted to humid forests in eastern Madagascar (i.e., Propithecus diadema, Indri indri, and Hapalemur simus) 48-50 , in caves located in the currently dry forests of Ankarana, close to our study area (~25 km south of Fantany). Conversely, it is plausible that Montagne d'Ambre acted as a refugium for rainforestadapted species during dry periods and facilitated recolonizations during subsequent humid periods 8,51 . Our interdisciplinary study revealed a strong mid-Holocene population decline of M. arnholdi, which occurred under precipitation decrease related to the termination of the AHP and thereby long before important human impact started in Madagascar 52 . Such a mid-Holocene population decline was also reported in multiple Malagasy taxa, including lemurs 12 , frogs 53 , endemic Eagle species 54 , and in the Malagasy giant rodents 55 . The first population decline was followed by a second decline in M. arnholdi at the beginning of the last millennium, which again parallels previously documented population bottlenecks for other lemur species 17,56,57 , all attributed to human-induced changes in forest habitats. Based on local paleoenvironmental reconstructions, we demonstrate that the latest population decline of M. arnholdi indeed coincided with a large increase of fire activity accompanied by a distinct ecosystem shift to a grassland-dominated landscape in the lowland areas surrounding Montagne d'Ambre. These changes most likely resulted from a regional increase in human activities and a shift to agro-pastoralism practices 58 . Although natural climatic changes might have also been relevant for ecosystem shifts in lowland areas during the past millennium, our results imply that human impact was certainly one of the dominant drivers in these latest and strong ecosystem changes in northern Madagascar.

Conclusion
Our dual paleoecological-genomic approach demonstrates signatures of both past natural climate change and strong recent anthropogenic impact in northern Madagascar. While demographic dynamics of M. arnholdi were exclusively shaped by natural climate changes until the late Holocene, the dynamics during the last millennium were also shaped by anthropogenic impact. In particular, we argue that the population decline of M. arnholdi at 5 kyr was climatically driven by the termination of the AHP, but that human occupation combined with the increasing aridification after the AHP impacted this species even further. Therefore, the termination of the AHP had severe impacts also in northern Madagascar and potentially on large parts of the island, leading to substantial population fluctuations and declines in Malagasy wildlife prior to human impact 34,35 . This interdisciplinary study reveals that understanding extant biodiversity requires a critical knowledge of natural, climatically driven environmental and ecological changes during the past and particularly beyond the late Holocene when significant human impact started on the island. Additionally, this study illustrates that it is extremely valuable to use multidisciplinary datasets to validate demographic approaches and explore both population size and connectivity changes across different time scales.

Methods
Environmental settings. Montagne d'Ambre NP comprises an isolated volcanic massif in northern Madagascar that stretches 30 km from north to south and 15 km from east to west with an altitude range of 400-1475 m asl (Fig. 1a). The location close to the coast (Indian Ocean/Mozambique Channel) with frequent trade winds triggers high orographic rainfalls. This process, frequently observed in tropical mountains, is mainly explained by the altitude increase that cools the air causing frequent fogs and high orographic precipitation 59   Maximum-likelihood (ML) estimates were obtained from the run with the highest composite likelihood. The respective 95% confidence intervals (CI) were generated by blockbootstrapping. All population size estimates are given in number of haploid copies. Time changes were scaled considering GT = 2.5 years. N0POP = effective population size for each population after the most recent population decline; N1POP = effective population size for each population after the older population decline; NANC = ancestral population size; 2NM0 = average number of haploid immigrants entering the population per generation. T1 = time of the most recent population decline; T2 = time of the older population decline. See Fig. 3d for a schematic illustration of M7. In bold: ML estimates that lie within the 95% confidence intervals. Note that the NANC estimate should be treated with caution as it falls outside the 95% confidence interval range. and combined to a composite core. Using the XRF data, a Principal Component Analysis (PCA) was carried out, using the selected elements (Si, K, Ca, Ti, Mn, Fe, and Ni) as input variables. Prior to this multivariate statistic, the normalized elements were transformed using standard "z-transformation". Grain-size distribution measurements were realized on LM1A, while paleoecological analyses, including pollen and charcoal counting, were realized on LM1B. For the grain size analysis 1 cm 3 subsamples at an interval of 5 cm were taken and measured after destroying carbonates and organic matter according to standard protocols with HCl (15, 30% p.A.) resp. H 2 O 2 (10, 32% p.A.). Measurements were carried out with a laser diffraction particle-size analyzer (LS 13320 Beckman Coulter) in seven cycles of 60 s each. The first reproducible signal was considered as reliable and final distribution data were calculated using the Fraunhofer optical model. Statistical parameters were derived from the raw data using a modified version of the MS-Excel macro Gradistat (V.4.5 s) 64 . For pollen extraction,~0.5 cm 3 subsamples at an interval varying between 8 and 48 cm on LM1B were prepared following a standard chemical protocol. A minimum sum of 300 terrestrial pollen grains was counted for each subsample using a light microscope at 400x magnification. Pollen and fernspore percentages were calculated on the terrestrial pollen sum. For charcoalparticle extraction, 1 cm 3 of sediment was sampled every cm along the core LM1B and soaked in a 3% NaP 2 O 4 solution plus bleach for several hours to deflocculate sediments and oxidize organic matter. Samples were sieved through a 160 μm mesh and charcoal particles were counted using a stereomicroscope at x40 magnification coupled to a digital camera. Influx of charcoal particles was calculated using the age-depth model (number of particles/cm 2 /year). For a reliable and robust chronology, 19 radiocarbon measurements were obtained by Accelerator Mass Spectrometry radiocarbon analyses from bulk organic matter of 1 cm 3 sediment samples (Supplementary Table S1). All subsamples for radiocarbon dating were measured at Poznan Radiocarbon Laboratory (Poland) and are calibrated to calendar kiloyears before present (kyr) using the southern hemisphere terrestrial calibration curve SHCal13 65 Fig. 1a). The first site lies at 1073 m asl, while the second (approx. 18 km apart from Mahasarika) is located at slightly lower altitude at the ecotone between forest and grassland at 848 m asl. Field work took place between August and October 2017 (dry season). Mouse lemurs were captured by hand near the trees that they used during nighttime activity 66  DNA extraction, RADseq, and whole-genome sequencing. Total genomic DNA was extracted from ear biopsies using the DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer's protocol with few modifications 67 . Two genomic datasets were generated to reconstruct the demographic history of M. arnholdi: a restriction site-associated DNA sequencing (RADseq) dataset and whole-genome sequences. Library preparation and RAD sequencing of all 46 samples were conducted at the GeT-PlaGe platform (Toulouse, France). In parallel, two mouse lemur samples (one from Mahasarika and one from Fantany) were selected for whole-genome sequencing (Supplementary Table S2). Both samples were derived from females in order to avoid a bias due to potential sex differences in dispersal behavior, which is common in other mouse lemur species 68,69 . The raw genomic data were trimmed and aligned to the Microcebus murinus reference genome 70 . All analyses using the RADseq dataset were performed based on genotype likelihoods [71][72][73] , considering only individuals with a mean sequencing coverage >4X [73][74][75] . ANGSD 73 was used to call genotype likelihoods for sites present in at least 75% of the individuals, and with a base quality above 20, a mapping quality above 30, and with a minimum minor allele frequency above 0.05. For the analyses with the whole-genome sequences, we applied a missing data filter of 30% and controlled for read depth and quality of the variant sites by considering a base quality above 20, a mapping quality above 30, a minimum read depth of three, and a maximum read depth of 100 (see Supplementary sections 2.1.2−2.1.4 for further methodological details and 2.2.1 and Table S3 for details about the genomic datasets).
Relatedness analyses and detection of population structure. In order to overcome a potential demographic history reconstruction bias 44,45,76 , the relatedness between the 46 sampled mouse lemurs was investigated first. The relatedness between two individuals is usually described by the concept of identity-by-descent (IBD). The software NGSrelate 74 implemented in ANGSD framework 73 was used to calculate the IBD coefficients (k 0 , k 1 , and k 2 ) between all pairs of individuals of our dataset (see Supplementary section 2.1.5 for details). Of each dyad of closely related individuals (e.g., parent-offspring or full siblings), only one partner was retained in the dataset. Subsequently, NGSadmix 72 was used to investigate population structure of 38 non-related individuals (12 from Mahasarika and 26 from Fantany). NGSadmix uses a maximum-likelihood clustering method to assign individuals to a predefined number of populations (K) based on genotype likelihoods 72 . A total of 10 independent runs of NGSadmix were performed to assign all individuals to two (K = 2) or three (K = 3) ancestral populations, assuming that the existence of more than three populations is not realistic at our local scale. The most suitable value of K was determined with Clumpak 77 , following the Evanno method 78 . The population structure analyses support the existence of two genetic clusters with low levels of admixture between the two sites ( Fig. 1b; Supplementary section 2.2.2 and Fig. S7). Consequently, the two sampling sites were considered as two independent demes in all downstream analyses.
Generation time. Demographic analyses are known to be affected by generation time (GT) 79 . Different generation time values ranging from 1 to 4.5 years can be found in the literature on mouse lemurs 11,17,80,81 . A GT = 1.0 was widely used in early demographic studies 17,80 and corresponds to the age at first reproduction. A longer generation time (3.0-4.5 years) has been estimated based on survival data for a population of M. rufus in eastern Madagascar 11 . Recently, an intermediate GT value of 2.5 years was estimated for a free-living M. murinus population in northwestern Madagascar based on parentage data 81 , where the authors defined GT as the average age of parents. A GT = 2.5 seems biologically appropriate for mouse lemurs, because it considered the entire reproductive career of individuals and took into account age-and sex-related differences in reproductive performance 81 . This estimate was used to generate all genomic results presented in the main text. However, in order to evaluate the consequences of choosing different generation times, all demographic analyses were repeated using the three aforementioned GTs (1.0, 2.5, and 4.5 years) and the results are presented in the Supplementary section 2.2.4.
PSMC analyses and IICR. PSMC 19 analyses were performed on two whole-genome autosomal sequences considering a minimum read depth per site of 3 (-d3) and a maximum read depth per site of 100 (-D100). The upper limit of the TMRCA was set to 5 (-t5) and the initial θ/ρ value to 5 (-r5). The N e was inferred using a total of 100 bootstrap replicates. PSMC plots were scaled using generation times of 1.0, 2.5, and 4.5 years per generation (see Supplementary section 2.2.4 and Fig. S15), and a mutation rate of 1.2 × 10 −811,82,83 . This mutation rate is the most accurate estimate available for mouse lemurs and was calculated from average pedigree-based estimates of seven primate species, including that of one mouse lemur species 83,84 . To evaluate the impact of the relatively low genome-wide coverage in our PSMC inferences, we also ran simulations for variant coverage divergence based on minimum read depth option per site (-d) between 1 and 9 for both study sites. The results are presented in the Supplementary section 2.2.5 and Fig. S18. Recent work has shown that PSMC dynamics can also result from changes in connectivity and not only from population size change, if populations are structured (e.g., n-island model, Supplementary section 2.1.6) 44,47 . It was shown that the PSMC method actually estimates a time-and sample-dependent function called the Inverse Instantaneous Coalescent Rate (IICR) whose interpretation is model-dependent (i.e., population panmixia vs population structure) 44 . The IICR 44,47 can be predicted under a very large family of structured models using simulations to identify complex models that can replicate the observed PSMC curve. In order to investigate if the demographic trajectories revealed by the PSMC method could be explained exclusively by changes in population connectivity, the IICR of Fantany was computed considering an n-island model of migration with a constant population size. The model simulations were performed by varying the number of connectivity changes between two and four, the migration rate between 0 and 100 and considering the number of islands up to 100. Genomic data was simulated using the ms software 85 , and the IICR was computed and plotted with a python script available at: https://github.com/willyrv/IICREstimator 47 , assuming a mutation rate of 1.2 × 10 −8 and the three generation times (see Supplementary section 2.2.4 and Fig.  S16). The IICR obtained under the best structure model was finally compared to the one inferred by the PSMC method.
Stairway Plot. The complete RADseq dataset of non-related animals (N = 38 individuals) was used for estimating a Site Frequency Spectrum (SFS; i.e., distribution of the allele frequencies of a given set of single nucleotide polymorphisms in a population) 16 . As a first step, ANGSD 73 was used to calculate the site allele frequency likelihoods based on individual genotype likelihoods. The realSFS tool implemented in ANGSD 73 was then used to estimate the 1d-SFS for each sampling site separately. Given that there is no suitable outgroup available to determine the ancestral state of each allele, we considered the minor allele-frequency spectrum (i.e., folded SFS) 18 for all SFS-related analyses. The two-folded 1d-SFS were finally used as input data for Stairway Plot v2.0 21 . Inferences were made based on 200 SFS for each sampling site, as suggested by the authors. The remaining 199 SFS were generated by bootstrapping, using the script provided by the software 21 . The plots were scaled using all three generation times (see Supplementary section 2.2.4 and Fig. S14) and 1.2 × 10 −8 as mutation rate 11,82 . The analyses were performed considering the entire dataset for each forest site, and also an equal sample size for Mahasarika and Fantany (N = 12). For this analysis, the individuals from Fantany were randomly selected using R (R Development Core Team 2005).
Demographic modeling with fastsimcoal2. A total of 13 potential and plausible demographic scenarios were designed for M. arnholdi and their composite likelihoods were compared following a hierarchical approach. First, five simple models, including a constant population size model, a single-population-size change model (expansion or bottleneck), an ancient or recent population structure model, and a connectivity change model, were explored (M1-M5; Supplementary sections 2.2.3, Fig. S8). Second, more complex demographic models were formulated based on the results of the previous step (Supplementary Table S6). Specifically, we compared four recently structured models (i.e., models where a panmictic population become recently structured) with the inclusion of one or two population bottlenecks after the populations become structured, or an ancient population expansion preceding that event (M6-M9; Supplementary Fig. S9a). Additionally, we also compared four models considering changes in connectivity among populations (with or without size change), where the changes in migration rate occurred in parallel to the population structuring or afterward (M10-M13; Supplementary Fig. S9b). See Supplementary section 2.2.3 and Table S6 for details of model selection. The analyses were performed with the RADseq dataset but, for computational reasons, only 20 individuals (i.e., 10 individuals from each sampling site) were selected for the 2d-SFS estimation. As a first step, ANGSD 73 was used to calculate the allelefrequency likelihoods for each sampling site, and the realSFS tool implemented in the same program combined together with custom scripts was used to estimate the two-dimensional minor allele frequency spectrum (i.e., folded 2d-SFS 18 ). It is important to emphasize that (i) Mahasarika is about 18 km apart from Fantany, (ii) mouse lemurs are known to disperse over very small distances (e.g., below 1 km for M. murinus) [86][87][88][89][90] , and (iii) our genetic structure analyses revealed only low levels of admixture among our study sites (Fig. 1b). For these reasons, an additional "ghost population" was included in our demographic models to control for the effect of an unsampled subpopulation on the estimation of migration rates between Mahasarika and Fantany (i.e., by representing a potential source of admixture between the northern and southern sites) 91,92 . The inclusion of such a "ghost population" is recommended when gene flow among the sampled populations is limited 92 , and it was shown to improve the fit of the demographic models in human populations 92 . Fastsimcoal2 v.2.6 20 was subsequently used to compute expected SFSs from simulations under the 13 demographic scenarios (see Supplementary Table S4 for details about model parameters). The best demographic model was selected according to three criteria. The Akaike Information Criteria (AIC) 93 was initially used to rank all 13 demographic models, and the models with the lowest ΔAIC (i.e., the models that minimized the distance between observed and expected SFS) were selected as the best ranked demographic models (Supplementary Table S6). Second, the fitting of the 2d-SFS obtained for the parameters that maximized the likelihood under model simulations were visually compared with the 2d-SFS estimated from the RADseq dataset 94 (see Supplementary Fig. S11-S13). Third, a block-bootstrap approach was used to calculate 95% confidence intervals (CI) for the parameters estimated under the three best ranked demographic models 20 (see Table 1 and  Supplementary Table S5 Impact of protein-coding sites in the demographic modeling of M. arnholdi. We reran the PSMC, Stairway Plot, and the best fastsimcoal2 demographic model without protein-coding sites to exclude possible confounding effects of sites under selection. All the analyses were performed using the same command options aforementioned and considering 2.5 years as generation time. The results are discussed in the Supplementary section 2.2.6 and shown in Figs. S19-S21.
Statistics and reproducibility. The Methods section and supplementary materials provide all details for statistical analyses. Specific details are also provided in figures when appropriate. All experiments and analyses performed in this paper are entirely reproducible.