domain: spi phys nlin

text:
Laurette S Tuckerman

Focus on Fluids 1 Ricocheting inclined layer convection states

Keywords: Bénard convection, pattern formation, bifurcations

Inclining a fluid layer subjected to a temperature gradient introduces a profusion of fascinating patterns and regimes. Previous experimental and computational studies form the starting point for an extensive numerical bifurcation study by Reetz & Schneider (2020a); Reetz et al. (2020b). Intricate trajectories passing through multiple steady and periodic states organize the dynamics. The consequences for chaotic patterns in large geometries is discussed.

Introduction: inclined layer convection

The phenomenon of thermal convection is known to all: warm fluid rises and cool fluid falls. But the fluid cannot rise and fall in the same locations; the way in which the rising and falling locations are arranged in space provides the archetypal example of pattern formation. Reetz and coworkers (2020a;2020b) explore a generalization of Rayleigh-Benard convection in which the confining plates held at different temperatures are inclined at an angle γ with respect to gravity. From a physical point of view, the immediate consequence is a shear flow in the direction of inclination, upwards along the warmer plate and downwards along the cooler one. From a mathematical point of view, the inclination renders the fluid layer anisotropic, distinguishing the direction of inclination and that perpendicular to it. As the Rayleigh number is increased, inclined convection undergoes a first transition from a featureless base state B to straight rolls. For compressed CO 2 with a Prandtl number of 1.07, these rolls are buoyancy-driven, with axes parallel to the inclination (longitudinal rolls, LR seen in the title figure) for γ 78 • , and shear-driven, with axes perpendicular to the inclination (transverse rolls, T R seen in the title figure) for γ 78 • (e.g. [START_REF] Gershuni | Stability of plane-parallel convective motion with respect to spatial perturbations[END_REF][START_REF] Chen | Stability of free-convection flows of variable-viscosity fluids in vertical and inclined slots[END_REF], as shown in Figure 1. The limiting case in which the bounding plates are vertical (γ = 90 • ) and the imposed temperature gradient is horizontal is sometimes called natural convection.

The problem of inclined layer convection received a boost from experimental observations of exotic regimes such as crawling rolls and transverse bursts by Daniels and co-workers [START_REF] Daniels | Pattern Formation in Inclined Layer Convection[END_REF][START_REF] Daniels | Defect Turbulence in Inclined Layer Convection[END_REF][START_REF] Daniels | Localized Transverse Bursts in Inclined Layer Convection[END_REF]. These were followed by numerical and theoretical simulations of these patterns by timestepping [START_REF] Daniels | Competition and bistability of ordered undulations and undulation chaos in inclined layer convection[END_REF] and by Floquet and Galerkin analyses [START_REF] Subramanian | Spatio-temporal patterns in inclined layer convection[END_REF].

I B LR OW R OW R W R (a) (b) (c)
trajectory evolution from the unstable base state B at = 40 and 0.5). The trajectory visits unstable LR, followed by a transient visit of Subsequently, the trajectory undergoes a sequence of rapid oscillations ttracted to a heteroclinic cycle between equilibrium state OW R and ated equilibrium. The time series for 2000 < t < 10 000 in inset (a) creasing time spent near the equilibrium states. The phase portrait of cycle is magnified in inset (b). Since the energy transfer rates do not try related states, the heteroclinic cycle appears as a homoclinic cycle. r a projection that distinguishes the two equilibrium states. e ⌧ xy OW R (Figure 7). Here, ⌧ xy = ⌧ (0.25, 0.25) is a shift operator transflow structure by half a pattern wavelength in the direction of the domain hich the wavy convection rolls of OW R are aligned. OW R and ⌧ xy OW R l periods of the wavy pattern along the domain diagonal and both are transformations of S owr = h⇡ xyz , ⌧(0.5, 0.5)i. Without imposing the sym-OW R and ⌧ xy OW R have each a single purely real unstable eigenmode, vector e u and ⌧ xy e u , respectively, that breaks the ⌧ (0.5, 0.5)-symmetry one spatial period along the domain diagonal (Figure 7c). Perturbing itiates the temporal transition OW R ! ⌧ xy OW R. Perturbing ⌧ xy OW R ates the temporal transition ⌧ xy OW R ! OW R. Together, the two dyions form a heteroclinic cycle. The two involved unstable eigenmodes eserve the symmetries ⇡ xyz and ⇡ xyz ⌧ (0.5, 0.5), respectively. Thus, they o sine-and cosine-eigenmodes in that they first, are orthogonal such er product is he u , ⌧ xy e u i = 0 and second, have a reflection symmetry di↵erent reflection points, namely ⇡ xyz and ⇡ xyz ⌧ (0.5, 0.5) (Figure 7c). ⇡ xyz -symmetry, the unstable eigenmode ⌧ xy e u is disallowed, ⌧ xy OW R ically stable, and the associated symmetry subspace ⌃ 1 contains only the ction OW R ! ⌧ xy OW R. When imposing ⇡ xyz ⌧ (0.5, 0.5)-symmetry, the ode e u is disallowed, OW R becomes dynamically stable, and the associubspace ⌃ ⌧ contains only the dynamical connection ⌧ xy OW R ! OW R. roclinic cycle satisfies all three conditions for the existence of a robust e between two symmetry related equilibrium states (Krupa 1997): saddle and ⌧ xy OW R is an attractor (or sink) in a symmetry subspace entire state space of the [2 x , y ]-periodic domain.

saddle-attractor connection OW R ! ⌧ xy OW R in ⌃ 1 .

Overview: a numerical bifurcation study

Reetz and coworkers (2020a; 2020b) have sought the bifurcation-theoretic origin of these exotic patterns. Their motivation is in part because of the features that inclined layer connection shares with wall-bounded shear flows such as plane Couette flow and pipe flow. These shear flows were not known to have any non-trivial solutions until their resistance to bifurcation-theoretic approaches was breached by [START_REF] Nagata | Three-dimensional finite-amplitude solutions in plane Couette flow: bifurcation from infinity[END_REF] and [START_REF] Clever | Three-dimensional convection in a horizontal fluid layer subjected to a constant shear[END_REF], who continued solutions from Taylor-Couette flow and Rayleigh-Bénard convection to plane Couette flow. Shortly thereafter, traveling waves in pipe flow were computed by [START_REF] Faisst | Traveling waves in pipe flow[END_REF] and Wedin & Kerswell (2004). Since then, dozens of solutions have been computed (e.g., [START_REF] Gibson | Visualizing the geometry of state space in plane Couette flow[END_REF]), all unstable.

These multiple solutions are thought to be of more than zoological interest, since [START_REF] Cvitanović | Periodic orbit expansions for classical smooth flows[END_REF]; [START_REF] Kawahara | The Significance of Simple Invariant Solutions in Turbulent Flows[END_REF] have proposed that turbulence in shear flows could be viewed as a collection of trajectories ricocheting between these multiple solutions via the connections between them. Indeed, for pipe flow, [START_REF] Hof | Experimental observation of nonlinear traveling waves in turbulent pipe flow[END_REF] compared exact solutions with snapshots during the evolution of the flow.

Reetz & Schneider (2020a) focus on several key states and dynamical regimes. One example is a sequence of bifurcations at γ = 40 • from the base state to longitudinal rolls to wavy rolls to oblique wavy rolls. The oblique wavy rolls participate in a robust heteroclinic cycle, shown in figure 2. Reetz et al. (2020b) expand on these results by surveying the parameter space of Rayleigh number and inclination angle to construct complete bifurcation diagrams and interpreting the transitions in the context of the large-aspect-ratio experiments. They have computed eight different invariant solutions all at Rayleigh number 21 266 and angle γ = 90 • , where the longitudinal rolls (LR) do not exist: the transverse rolls (TR) already mentioned, and various periodically modulated versions of the longitudinal rolls. Carrying out direct numerical simulations at the same parameter values, they extracted snapshots that resembled these states, shown in figure 2.

Future: from bifurcation diagrams to turbulence?

There is no doubt that patterns and temporal behavior are controlled and explained by the plethora of underlying dynamical objects -fixed points and traveling waves, periodic orbits and heteroclinic cycles. The ability to compute such objects for the full three-dimensional Navier-Stokes and Boussinesq equations has resulted from advances in several fields: first, the spectacular growth of dynamical systems theory and symmetry (e.g., the monographs by [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory: Volume II[END_REF][START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF] following the discovery of deterministic chaos; second, the increasing power of computation; and third, the incorporation of matrix-free methods for linear algebra [START_REF] Dijkstra | Numerical Bifurcation Methods and their Application to Fluid Dynamics: Analysis beyond Simulation[END_REF] into algorithms for tracking dynamical objects. See, e.g., [START_REF] Marques | Centrifugal effects in rotating convection: axisymmetric states and three-dimensional instabilities[END_REF]; Borońska & Tuckerman (2010a,b) for computational bifurcation studies of convection. However, the application of the theories of [START_REF] Cvitanović | Periodic orbit expansions for classical smooth flows[END_REF]; [START_REF] Kawahara | The Significance of Simple Invariant Solutions in Turbulent Flows[END_REF] to a turbulent hydrodynamic state would require the computation of an even larger number of dynamical objects, on a scale that is not yet possible. The challenge is to bridge the gap between bifurcation theory and the large-scale statistical phenomenon of turbulence.

Figure 1 .

 1 Figure 1. Left: evolution from base state B to longitudinal rolls LR to wavy rolls W R (inset (c)), followed by a sequence of widening oscillations leading to a heteroclinic cycle between the oblique wavy roll state OW R and its image under a shift τxy in x and y (inset (b)). The trajectory is projected onto coordinates (D/I, I), with D the viscous dissipation and I the energy input. Right: temperature isosurfaces of W R and OW R. Adapted from Reetz & Schneider (2020a).

Figure 2 .

 2 Figure 2. (a)-(h) Eight invariant states in vertical layer convection (γ = 90 • ), all at Rayleigh number 21 266, computed by continuation. (i)-(l) Four snapshots from DNS of turbulent regime at the same parameters. Adapted from Reetz et al. (2020b).

Conflicts of interest

The author reports no conflict of interest.


Jelmer Poelstra

Jordi Salmona

George P Tiley

Dominik Schüssler

Marina B Blanco

Jean B Andriambeloson

Olivier Bouchez

Ryan Campbell

Paul D Etter

Paul A Hohenlohe

Dominik Schüßler

Kelsie E Hunnicutt

Amaia Iribar²

Eric A Johnson

Peter M Kappeler

Peter A Larsen

Sophie Manzi

José M Ralison

Blanchard Randrianambinina

Rodin M Rasoloarison

David W Rasolofoson

Amanda R Stahlke

David Weisrock

Rachel C Williams

Lounès Chikhi

Edward E Louis

Ute 14

Radespiel

Anne D Yoder
email: anne.yoder@duke.edu

Ute Radespiel
email: ute.radespiel@tiho-hannover.de

Cryptic Patterns of Speciation in Cryptic Primates: Microendemic Mouse Lemurs and the Multispecies Coalescent

Keywords: effective population size, cryptic species, multispecies coalescent, species delimitation, speciation, microendemism

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Mouse lemurs (Microcebus spp.) are small, nocturnal primates that are widespread in the forests of Madagascar [START_REF] Mittermeier | Lemurs of Madagascar, 3 rd Edition[END_REF], one of the world's most biodiverse environments [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF][START_REF] Goodman | Updated estimates of biotic diversity and endemism for Madagascar[END_REF][START_REF] Estrada | Impending extinction crisis of the world's primates: Why primates matter[END_REF]. Mouse lemur diversity was long overlooked [START_REF] Zimmermann | Species concepts, diversity, and evolution in primates: Lessons to be learned from mouse lemurs[END_REF] until the introduction of genetic analyses made it feasible to identify diverging lineages comprised of lemurs with similar phenotypes and ecological niches. This genetic perspective has led to the description of many new species, with 24 species recognized at present. In one such study, [START_REF] Radespiel | Exceptional diversity of mouse lemurs (Microcebus spp.) in the Makira region with the description of one new species[END_REF] surveyed the forests of the Makira region (Fig. 1) in northeastern Madagascar and found evidence for three divergent mitochondrial lineages occurring in sympatry. One of these was identified as M. mittermeieri [START_REF] Louis | Revision of the mouse lemurs (Microcebus) of eastern Madagascar[END_REF], while the second was newly described as M. macarthurii. A third lineage, provisionally called M. sp. #3, was hypothesized to represent a new species closely related to M. macarthurii but was not formally named because the data were limited to mtDNA sequence data from one individual. Furthermore, two other species occur in the region, M. lehilahytsara (Roos and Kappeler in [START_REF] Kappeler | Morphology, behaviour and molecular evolution of giant mouse lemurs (Mirza spp.) Gray, 1870, with description of a new species[END_REF]) at higher elevations, and M. simmonsi [START_REF] Louis | Revision of the mouse lemurs (Microcebus) of eastern Madagascar[END_REF] in lowland forests in the south (Fig. 1).

Given that many previous taxonomic descriptions of mouse lemurs have relied strongly, if not entirely, on mtDNA sequence divergence, there has been criticism that mouse lemurs (and CRYPTIC SPECIATION IN MOUSE LEMURS 8 lemurs more generally) may have been oversplit [START_REF] Tattersall | Madagascar's lemurs: Cryptic diversity or taxonomic inflation?[END_REF][START_REF] Markolf | On species delimitation: Yet another lemur species or just genetic variation?[END_REF]). Species delimitation using only mtDNA is now widely regarded as problematic, given that the mitochondrial genome represents a single non-recombining locus whose gene tree may not represent the underlying species tree (e.g., [START_REF] Pamilo | Relationships between gene trees and species trees[END_REF][START_REF] Maddison | Gene trees in species trees[END_REF]. Mitochondria are also maternally inherited and therefore susceptible to effects of male-biased dispersal (e.g., [START_REF] Dávalos | Sex-biased dispersal produces high error rates in mitochondrial distance-based and tree-based species delimitation[END_REF], which is prevalent in mouse lemurs (reviewed in [START_REF] Radespiel | Can behavioral ecology help to understand the divergent geographic range sizes of mouse lemurs?[END_REF].

Moreover, previous attempts to resolve mouse lemur relationships using nuclear sequences have been complicated by high gene tree discordance, consistent with strong incomplete lineage sorting (e.g., [START_REF] Heckman | Multiple nuclear loci reveal patterns of incomplete lineage sorting and complex species history within western mouse lemurs (Microcebus)[END_REF][START_REF] Weisrock | Delimiting species without nuclear monophyly in Madagascar's mouse lemurs[END_REF]. These issues can be overcome with genomic approaches, which provide power for simultaneously resolving phylogenetic relationships and estimating demographic parameters such as divergence times, effective population sizes, and rates of gene floweven among closely related species (e.g., [START_REF] Palkopoulou | A comprehensive genomic history of extinct and living elephants[END_REF][START_REF] Pedersen | A southern African origin and cryptic structure in the highly mobile plains zebra[END_REF].

Given that cryptic species are by definition difficult to identify based on phenotypic characters [START_REF] Bickford | Cryptic species as a window on diversity and conservation[END_REF], recently developed methods for genomic species delimitation have advanced our ability to recognize and quantify their species diversity. In the past decade, both theory and methods for species delimitation have seen substantial progress, especially those which leverage the multispecies coalescent (MSC) model [START_REF] Pamilo | Relationships between gene trees and species trees[END_REF][START_REF] Rannala | Bayes estimation of species divergence times and ancestral population sizes using DNA sequences from multiple loci[END_REF]. MSC-based species delimitation methods have been increasingly applied to genomic data (e.g. [START_REF] Carstens | Species delimitation using a combined coalescent and informationtheoretic approach: an example from North American Myotis bats[END_REF][START_REF] Yang | Bayesian species delimitation using multilocus sequence data[END_REF][START_REF] Grummer | Species delimitation using Bayes factors: simulations and application to the Sceloporus scalaris species group (Squamata: Phrynosomatidae)[END_REF][START_REF] Dincă | The conundrum of species delimitation: a genomic perspective on a mitogenetically super-variable butterfly[END_REF][START_REF] Hundsdoerfer | Genome-wide SNP data reveal an overestimation of species diversity in a group of hawkmoths[END_REF]), though they have also been considered controversial [START_REF] Edwards | Species detection and individual assignment in species delimitation: can integrative data increase efficacy?[END_REF][START_REF] Sukumaran | Multispecies coalescent delimits structure, not species[END_REF][START_REF] Barley | Impact of model violations on the inference of species boundaries under the multispecies coalescent[END_REF].

The controversy largely relates to the idea that strong population structure can be mistaken for species boundaries, which may lead to oversplitting [START_REF] Jackson | Species delimitation with gene flow[END_REF][START_REF] Sukumaran | Multispecies coalescent delimits structure, not species[END_REF][START_REF] Luo | Comparison of methods for molecular species delimitation across a range of speciation scenarios[END_REF][START_REF] Leaché | The spectre of too many species[END_REF][START_REF] Chambers | The multispecies coalescent over-splits species in the case of geographically widespread taxa[END_REF]). To overcome this potential weakness, [START_REF] Jackson | Species delimitation with gene flow[END_REF] proposed a heuristic criterion, the genealogical divergence index (gdi), with [START_REF] Leaché | The spectre of too many species[END_REF] further suggesting that gdi helps to differentiate between population structure and species-level divergence. In parallel, sophisticated statistical approaches have been developed that can detect the presence and magnitude of gene flow during or after speciation [START_REF] Gronau | Bayesian inference of ancient human demography from individual genome sequences[END_REF][START_REF] Dalquen | Maximum likelihood implementation of an isolation-with-migration model for three species[END_REF][START_REF] Wen | Coestimating reticulate phylogenies and gene trees from multilocus sequence data[END_REF]. Taken together, these analytical developments are crucial to our ability to recognize the patterns that characterize the speciation process, despite the challenge of identifying species without universally agreed upon criteria [START_REF] De Queiroz | Species concepts and species delimitation[END_REF].

In this study, we use a structured framework starting with phylogenetic placement of lineages and culminating with the MSC to delimit species, estimate divergence times, identify post-divergence gene flow, and to estimate both current and ancestral effective population sizes (Fig. S1). We take advantage of increased geographic, population-level, and genomic sampling to comparatively examine speciation dynamics for two pairs of closely related lineages in the CRYPTIC SPECIATION IN MOUSE LEMURS 10 region (described below as Clades I and II) and perform MSC species delimitation methods with Restriction-site Associated DNA sequencing (RADseq) data to infer divergence times, effective population sizes, and rates of gene flow between these lineages. We also provide a novel wholegenome assembly for the previously undescribed lineage and compare inferences of effective population size (N e ) through time from whole-genome versus RADseq data. We find notably different species delimitation results for the lineages in the two mouse lemur clades and believe that the comprehensive analytic framework here used can be applied more generally to allow investigators to test hypotheses of population-versus species-level differentiation.

MATERIALS AND METHODS

Summary of Analyses

We generated RADseq data for 63 individuals from 6 lineages, of which 48 were from the two focal clades and passed quality control. First, we used maximum likelihood approaches to infer relationships among lineages and to provide a framework for subsequent species delimitation analyses (Fig. 2A andC). To delimit species, we performed clustering (Fig. 2B) and PCA analyses (Fig. 3A-C), as well as formal MSC species delimitation analyses using SNAPP and BPP. We also used the recently developed genealogical divergence index gdi based on BPP parameter estimates (Fig. 3D) and performed an isolation-by-distance analysis (Fig. 4). To determine to what extent ongoing and ancestral gene flow may have contributed to current patterns of divergence, we used G-PhoCS and D-statistics (Fig. 5). Finally, we generated and assembled whole-genome sequencing data for a single M. sp. #3 individual. The whole-genome sequences of M. sp. #3 reported here for the first time and M. mittermeieri from a previous study [START_REF] Hunnicutt | Comparative genomic analysis of the pheromone receptor class 1 family (V1R) reveals extreme complexity in mouse lemurs (Genus, Microcebus) and a chromosomal hotspot across mammals[END_REF] were used to infer N e though time with Multiple Sequentially Markovian Coalescent (MSMC) analysis and compared with estimates from G-PhoCS (Fig. 6). Below, we describe the methods in more detail, while further details can be found in the Supplementary Material.

Study Sites and Sampling

Microcebus samples were obtained by taking ~2 mm² ear biopsies of captured (and thereafter released) individuals between 2008 and 2017 at seven humid evergreen forest sites (50-979 m a.s.l.) in the Analanjirofo and Sava regions of northeastern Madagascar (Fig. 1; Table S1). Additional samples were used from Riamalandy, Zahamena National Park (NP), Betampona Strict Nature Reserve (SNR) and Tampolo [START_REF] Louis | Revision of the mouse lemurs (Microcebus) of eastern Madagascar[END_REF][START_REF] Weisrock | Delimiting species without nuclear monophyly in Madagascar's mouse lemurs[END_REF] Louis and Lei 2016) (Fig. 1). With this sampling strategy, we expected to include all mouse lemur species thought to occur in the region (from north to south): M. mittermeieri, M. macarthurii, M. sp. #3, M. lehilahytsara, and M. simmonsi (Fig. 1). Microcebus murinus, which occurs in western and southeastern Madagascar, was used as an outgroup.

Sequencing Data, Genotyping and Genome Assembly

We generated RADseq libraries using the SbfI restriction enzyme, following three protocols (Supplementary Methods, Table S1). Sequences were aligned to the M. sp. #3 nuclear genome generated by this study, and to the published M. murinus mitochondrial genome [START_REF] Lecompte | Complete mitochondrial genome of the gray mouse lemur, Microcebus murinus (Primates, Cheirogaleidae). Mitochondrial DNA A DNA Mapp[END_REF]. We used two genotyping approaches to ensure robustness of our results.

First, we estimated genotype likelihoods (GL) with ANGSD v0.92 [START_REF] Nielsen | SNP calling, genotype calling, and sample allele frequency estimation from new-generation sequencing data[END_REF][START_REF] Korneliussen | ANGSD: Analysis of Next Generation Sequencing Data[END_REF], which retains information about uncertainty in base calls, thereby alleviating some issues commonly associated with RADseq data such as unevenness in sequencing depth and allele dropout [START_REF] Lozier | Revisiting comparisons of genetic diversity in stable and declining species: assessing genome-wide polymorphism in North American bumble bees using RAD sequencing[END_REF][START_REF] Pedersen | A southern African origin and cryptic structure in the highly mobile plains zebra[END_REF][START_REF] Warmuth | Genotype-free estimation of allele frequencies reduces bias and improves demographic inference from RADSeq data[END_REF]. Second, we called genotypes with GATK v4.0.7.0 [START_REF] Depristo | A framework for variation discovery and genotyping using next-generation DNA sequencing data[END_REF], and filtered GATK genotypes following the "FS6" filter of O' Leary et al. (2018;their Table 2). We furthermore used three mtDNA fragments [Cytochrome Oxidase II (COII), Cytochrome B (cytB), and d-loop] that were amplified and Sanger sequenced.

The genome of the M. sp. #3 individual sampled in Mananara-Nord NP (Table S3) was sequenced with a single 500bp insert library on a single lane of an Illumina HiSeq 3000 with paired-end 150bp reads. We used MaSuRCA v3.2.2 [START_REF] Zimin | The MaSuRCA genome assembler[END_REF] for contig assembly and SSPACE [START_REF] Boetzer | Scaffolding pre-assembled contigs using SSPACE[END_REF] for scaffolding. Scaffolds potentially containing mitochondrial or X-chromosome sequence data were removed in downstream analyses.

Phylogenetic Analyses

We used three phylogenetic approaches to infer relationships among lineages: (1) maximum likelihood using RAxML v8.2.11 [START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF], (2) SVDquartets, an MSC method that uses phylogenetic invariants, implemented in PAUP v4a163 [START_REF] Chifman | Quartet inference from SNP data under the coalescent model[END_REF], and (3) SNAPP, a full-likelihood MSC method for biallelic data that does not require joint gene tree estimation (v1.3.0; [START_REF] Bryant | Inferring species trees directly from biallelic genetic markers: bypassing gene trees in a full coalescent analysis[END_REF]. Analyses with RAxML and SVDquartets used all available individuals, whereas SNAPP analyses were performed with subsets of 12 and 22 individuals for computational feasibility (see Supplementary Methods).

Species Delimitation

Clustering approaches and summary statistics. -Clustering analyses were performed using corresponding methods based on ANGSD genotype likelihoods [clustering in NgsAdmix v32 [START_REF] Skotte | Estimating individual admixture proportions from next generation sequencing data[END_REF] and PCA in ngsTools va4d338d [START_REF] Fumagalli | ngsTools: methods for population genetics analyses from next-generation sequencing data[END_REF]] and on GATKcalled genotypes [clustering in ADMIXTURE v1.3.0 [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF]) and PCA using the glPca() function in adegenet v2.1.1 [START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF]]. These analyses were run for Clade I and Clade II together and separately.

MSC-based approaches. -We used SNAPP to test if the two lineages each in Clade I and II could be delimited using Bayes factors [START_REF] Leaché | Species delimitation using genome-wide SNP data[END_REF], interpreting 2ln Bayes factors greater than six as strong evidence for a given model [START_REF] Kass | Bayes Factors[END_REF]. We also applied gdi. -Coalescent node heights (τ) and ancestral effective population sizes (θ) estimated by BPP were used to compute the genealogical divergence index (gdi; [START_REF] Jackson | Species delimitation with gene flow[END_REF][START_REF] Leaché | The spectre of too many species[END_REF] for the lineages in Clade I and Clade II. We calculated gdi as in [START_REF] Leaché | The spectre of too many species[END_REF], using their equation 7 (gdi = 1-e -2τ/θ ), where 2τ/θ represents the population divergence time between two taxa in coalescent units. θ is taken from one of the two taxa and therefore, gdi was calculated twice for each species pair, alternating the focal taxon. We computed gdi using τ and θ parameter estimates for each posterior BPP sample to incorporate uncertainty in the estimates. [START_REF] Jackson | Species delimitation with gene flow[END_REF] suggested the following interpretation of gdi values: the taxon pair (a) is unambiguously a single species for gdi < 0.2, (b) is unambiguously two separate species for gdi > 0.7, and (c) falls in an ambiguous zone for 0.7 > gdi < 0.2.

Isolation-by-distance. -We tested for isolation-by-distance using the VCF file produced by GATK with the gl.ibd() function in the R package dartR 1.1.11 [START_REF] Gruber | dartr: An r package to facilitate analysis of SNP data generated from reduced representation genome sequencing[END_REF].

Inference of gene flow and divergence times

G-PhoCS v1.3 [START_REF] Gronau | Bayesian inference of ancient human demography from individual genome sequences[END_REF]), a Bayesian MSC approach that allows for the estimation of periods of gene flow (i.e. "migration bands"), was used to jointly infer divergence times, population sizes, and rates of gene flow between specific lineages. Based on the results of exploratory models that each contained a single "migration band" between two lineages, we ran a final model with a migration bands allowing gene flow from mittermeieri to lehilahytsara and from macarthurii to M. sp. #3. Given the observed mitonuclear discordance between M. sp. #3 and M. macarthurii (see Results), we investigated gene flow between them in more detail by running G-PhoCS using a dataset with only M. sp. #3, M. macarthurii, and M. lehilahytsara individuals, wherein M. sp. #3 was divided into the two populations detected using clustering approaches.

The D-statistic and related formal statistics for admixture use phylogenetic invariants to infer post-divergence gene flow between non-sister populations or taxa. We used the qpDstat tool in admixtools v4.1 [START_REF] Patterson | Ancient admixture in human history[END_REF] to compute four-taxon D-statistics for all possible configurations in which gene flow could be tested between non-sister lineages among the five ingroup lineages. We additionally tested for gene flow between M. macarthurii and M. sp. #3 by separately treating (1) the two distinct M. sp. #3 populations detected by clustering approaches, and (2) M. macarthurii individuals with and without "M. sp. #3-type" mtDNA (see Results). In all tests, M. murinus was used as P4 (outgroup).

Effective Population Size Through Time

Studies have shown that population structure can generate spurious signals of population size change [START_REF] Beaumont | Recent developments in genetic data analysis: what can they tell us about human demographic history?[END_REF][START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size change[END_REF][START_REF] Heller | The confounding effect of population structure on Bayesian skyline plot inferences of demographic history[END_REF]. For example, sequentially Markovian coalescent approaches such as MSMC (Schiffels and Durbin 2014) actually estimate the inverse instantaneous coalescence rate, which is only equivalent to an effective size in panmictic models [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution--lessons for ancestral population size inference?[END_REF][START_REF] Rodríguez | The IICR and the non-stationary structured coalescent: towards demographic inference with arbitrary changes in population structure[END_REF]. We therefore inferred and compared population size histories using two methods. We estimated N e over time with MSMC for two species, using the whole-genome data of M. sp. #3 and M. mittermeieri [START_REF] Hunnicutt | Comparative genomic analysis of the pheromone receptor class 1 family (V1R) reveals extreme complexity in mouse lemurs (Genus, Microcebus) and a chromosomal hotspot across mammals[END_REF]) mapped to the chromosome-level genome assembly of M. murinus [START_REF] Larsen | Hybrid de novo genome assembly and centromere characterization of the gray mouse lemur (Microcebus murinus)[END_REF]).

These estimates were compared to inferred changes in N e over time based on θ estimates from G-PhoCS for each predefined extant or ancestral population. Although G-PhoCS was not expressly developed to estimate change in N e over time, this allowed us to explicitly examine broad demographic trends, even with small population-level sampling.

Mutation Rate and Generation Time

We used empirical estimates of mutation rate and generation time to convert coalescent units from BPP, G-PhoCS and MSMC analyses into absolute times and population sizes. We incorporated uncertainty by drawing from mutation rate and generation time distributions for each sampled generation of the MCMC chains in BPP and G-PhoCS (MSMC parameter estimates were converted using the point estimates). For the mutation rate, we used a gamma distribution based on the mean (1.236 x 10 -8 ) and variance (0.107 x 10 -8 ) of seven pedigree-based mutation rate estimates for primates (see [START_REF] Campbell | Pedigree-based measurement of the de novo mutation rate in the gray mouse lemur reveals a high mutation rate, few mutations in CpG sites, and a weak sex bias[END_REF], Table S1). For the generation time, we used a lognormal distribution with a mean of ln(3.5) and standard deviation of ln(1.16) based on estimates of 4.5 years calculated from survival data [START_REF] Zohdy | Teeth, sex, and testosterone: aging in the world's smallest primate[END_REF][START_REF] Yoder | Geogenetic patterns in mouse lemurs (genus Microcebus) reveal the ghosts of Madagascar's forests past[END_REF]) from M. rufus, and 2.5 years from average parent age based on capture-mark-recapture and parentage data in the wild [START_REF] Radespiel | An empirical estimate of the generation time of mouse lemurs[END_REF] for M. murinus.

RESULTS

RADseq Data and Whole-Genome Assembly

We used three library generation protocols, two sequencing lengths, and a combination of single and paired-end sequencing, yielding data for all 63 individuals in the study and demonstrating the utility of cross-laboratory RAD sequencing, as previously shown in other taxa (e.g., [START_REF] Gonen | Exploring the utility of cross-laboratory RAD-sequencing datasets for phylogenetic analysis[END_REF]. From more than 447 million raw reads (Table S1), over 394 million passed quality filters, with approximately 182 million successfully aligning to the M. sp. #3 reference genome. We obtained an average of 120,000 loci per individual with coverage ranging from ~1 to ~22x (Table S1).

We assembled approximately 2.5 Gb of nuclear genome sequence data for M. sp. #3 with a contig N50 around 36 Kb (Table S3). While the final assembly was fragmented, as expected for a single Illumina library genome, only 6.4% of mammalian BUSCOs (Benchmarking Universal Single-Copy Orthologs, [START_REF] Simão | BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs[END_REF] were found to be missing. The genome sequence and associated gene annotations can be accessed through NCBI (Bioproject PRJNA512515).

Phylogenetic Relationships

RAxML and SVDquartets recovered well-supported nDNA clades for M. simmonsi, M. macarthurii, and M. sp. #3, the latter two as sister taxa with 100% bootstrap support (Fig. 2; Fig. S2). SNAPP also supported M. sp. #3 as sister taxon to M. macarthurii (referred to as Clade I) and placed M. lehilahytsara as sister taxon to M. mittermeieri (referred to as Clade II) (Fig. S2). However, M. lehilahytsara was not monophyletic in RAxML analyses of nDNA (Fig. 2C) or mtDNA (Fig. 2A), and a SVDquartets analysis of nDNA placed one M. lehilahytsara individual from Ambavala as sister to all other M. lehilahytsara and M. mittermeieri, and only weakly supported a monophyletic M. mittermeieri (Fig. S2A).

Although mtDNA analyses placed several individuals from Anjiahely in a well-supported clade with M. sp. #3, individuals from Ambavala (see Fig. 1), Mananara-Nord NP, and Antanambe (Fig. 2A; see lower gray box), nuclear RADseq data placed them unambiguously within the M. macarthurii clade (Fig. 2B,C). This suggests that individuals from Anjiahely are in fact M. macarthurii, but carry two divergent mtDNA lineages, and that true M. sp. #3 are only found between Ambavala and Antanambe (Fig. 1). The cause of this mitonuclear discordance for macarthurii in Anjiahely was investigated further (see the section "Interspecific Gene Flow").

Species Delimitation

Genetic structure. -A PCA with both pairs of sister lineages (Clade I: M. macarthurii and M. sp. #3; Clade II: M. mittermeieri and M. lehilahytsara) distinguished the two clades along PC1, and distinguished M. macarthurii and M. sp. #3 along PC2 (Fig. 3B). When restricting clustering analyses to Clade I, K = 2 was the best-supported number of clusters with both approaches, distinguishing M. macarthurii and M. sp. #3 (Fig. S5; Fig. S7B). At K = 3, M.

sp. #3 was divided into two clusters with individuals from Mananara-Nord NP and Antanambe separated from Ambavala individuals (Fig. S7B, Fig. S10). A separate PCA analysis for Clade I also distinguished these two groups along PC2 (Fig. 3C), which we hereafter refer to as "southern M. sp. #3" (Mananara-Nord NP and Antanambe are south of the Mananara river)

and "northern M. sp. #3" (Ambavala is north of the river, and 24.0 km from Mananara-Nord NP and 35.2 km from Antanambe; Fig. 1), respectively. When restricting clustering analyses to Clade II, ADMIXTURE and ngsAdmix suggested optimal values of 1 and 2, respectively; at K=2, M. mittermeieri and M. lehilahytsara were largely but not entirely separated by both approaches (Fig. S5, Fig. S7C, Fig. S11). A PCA distinguished M. mittermeieri and M. lehilahytsara along PC1 but with little separation (Fig. 3D).

SNAPP and BPP.

-SNAPP Bayes factors strongly favored splitting Clade I into two species (2lnBF = 17,304 and 34,326 for two different datasets, Table S6), as well as splitting Clade II, although with a smaller difference in marginal likelihood scores (2lnBF = 1,828 and CRYPTIC SPECIATION IN MOUSE LEMURS 20 993). All putative species assignments were recovered by the guided delimitation analysis with BPP (Fig. S12).

Genealogical divergence index (gdi).

-For the Clade I sister pair, gdi was 0.727 (95% HPD: 0.718-0.737) from the perspective of M. macarthurii (i.e. above the upper threshold for species delimitation), and 0.500 (0.488-0.511) from the perspective of M. sp. #3 (i.e. in the upper ambiguous zone for species delimitation; Fig. 3A). In contrast, gdi values for the Clade II putative species pair were much lower and even below the lower threshold for species delimitation: 0.080 (0.074-0.086) from the perspective of M. lehilahytsara, and 0.193 (0.187-0.201) from the perspective of M. mittermeieri (Fig. 3A).

Isolation-by-distance (IBD). -While comparisons within and between lineages appeared

to follow a single isolation-by-distance pattern for M. mittermeieri and M. lehilahytsara (Clade II, r=0.693, p=0.002, Fig. 4B), comparisons within versus between lineages differed strongly for M. macarthurii and M. sp. #3 (Clade I, Fig. 4A). Specifically, genetic distances between M. macarthurii and M. sp. #3 were much larger than within lineages and were also much larger than between M. mittermeieri and M. lehilahytsara, despite similar geographic distances. Clade I, from M. sp. #3 to macarthurii [2Nm = 0.08 (95% HPD: 0.07-0.09), migrants per generation: 0.10% (0.05-0.15%)] (Fig. 5C). G-PhoCS also inferred low levels of gene flow between the two clades, most likely between ancestral populations, but the timing and direction of gene flow could not be determined (Supplementary Results; Fig. S13), and D-statistics testing for gene flow between the clades were not significant (Fig. S14).

We further investigated gene flow between M. sp. #3 and M. macarthurii by taking the strong population structure within M. sp. #3 into account. D-statistics suggested that northern M. sp. #3 and M. macarthurii with "M. sp. #3-type" mtDNA share a slight excess of derived alleles in relation to southern M. sp. #3, significantly deviating from 0, which indicates gene flow (Fig. S15A). Using a G-PhoCS model with separate northern and southern groups for M. sp. #3, we found that (1) gene flow with M. macarthurii took place before and after the onset of divergence between northern and southern M. sp. #3, (2) gene flow between extant lineages occurred or occurs only between northern (and not southern) M. sp. #3 and M. macarthurii, and

(3) gene flow is asymmetric, predominantly into M. macarthurii (Fig. S15B).

Divergence Times

We estimated divergence times under the MSC model using BPP and G-PhoCS both with and without interspecific gene flow (Fig. 5; Fig. S17). Results were similar across these approaches, with the exception of divergence times between sister lineages in G-PhoCS models CRYPTIC SPECIATION IN MOUSE LEMURS 22 with versus without gene flow (Fig. 5). Specifically, the divergence time between M. sp. #3 and M. macarthurii (Clade I) without gene flow was estimated at 115 kya (95% HPD range: 52-190 ka across G-PhoCS and BPP models) (Fig. 5; Fig. S17), but at 193 kya (95% HPD: 89-318 ka) when incorporating gene flow (Fig. 5C-D). In Clade II, this difference in estimated divergence times was considerably larger: under an isolation model it was estimated to be 103 kya (95% HPD: 49-171 ka; Fig. 5) and as much as 520 kya (95% HPD: 249-871 ka) when modeled with gene flow (Fig. 5). Deeper nodes were not as strongly affected: divergence time between Clades I and II was estimated at 687 kya (95% HPD: 337-1126 ka) across G-PhoCS and BPP isolation models, and at 796 kya (95% HPD: 360-1311 ka) in a G-PhoCS model with gene flow (Fig. 5D).

Effective Population Sizes

We found large differences in N e among lineages, with considerably larger N e for the lineages in Clade II, M. lehilahytsara (mean estimate and 95% HPD range across the BPP and G-PhoCS models with and without interspecific gene flow: 159 k; 58-265 k) and M. mittermeieri (78 k; 36-140 k), than the lineages in Clade I, M. sp. #3 (24 k; 12-38 k) and M. macarthurii (12 k; 5-19 k) (Fig. 5A-C). Wide HPD intervals for M. mittermeieri and lehilahytsara are due to differences between models with and without gene flow. Using the G-PhoCS model focused on Clade I, fairly similar effective population sizes were estimated separately for northern (47 k; 17-78 k), southern (23 k; 12-37 k), and ancestral (33 k; 17-53 k) M. sp. #3 lineages (Fig. S13).

Using the whole-genome data for one individual of M. sp. #3 (from the southern group) and for M. mittermeieri, a comparison of MSMC analyses and G-PhoCS models with and without gene flow (Fig. 6) showed highly similar and markedly declining estimates of population sizes towards the present for M. sp. #3 (Fig. 6A). Estimates for M. mittermeieri were more variable across analyses but none showed a consistent decline towards the present (Fig. 6B).

DISCUSSION

We used an MSC-based framework for genomic species delimitation and identified rapid and recent diversification of mouse lemurs in a relatively small area in northeastern Madagascar.

The same region was previously identified to harbor high levels of lemur microendemism and to be vulnerable to the effects of climate change [START_REF] Brown | Shifting ranges and conservation challenges for lemurs in the face of climate change[END_REF] and anthropogenic habitat alteration [START_REF] Schüßler | Identification of crucial stepping stone habitats for biodiversity conservation in northeastern Madagascar using remote sensing and comparative predictive modeling[END_REF], marking it as a region of conservation concern. Species-level divergence was strongly supported for M. sp. #3 and its sister species M. macarthurii (Clade I, Fig. 2), but not for the pair of M. mittermeieri and M. lehilahytsara (Clade II, Fig. 2), overturning our a priori expectation that the latter were distinct species [START_REF] Olave | Upstream analyses create problems with DNA-based species delimitation[END_REF]. We inferred that the focal species all diverged from their common ancestors within the past million years and documented two cases of sympatric occurrence, each with one representative from Clade I and one from Clade II. The combined findings of recent divergence and sympatric overlap suggest that reproductive isolation can evolve rapidly in mouse lemurs.

Support for Separate Sister Species Differs Sharply Between the Two Clades

Evidence for distinguishing M. sp. #3 and M. macarthurii as separate species was strong and consistent across analyses. They were reciprocally monophyletic across all phylogenetic analyses of RADseq data (Fig. 2C; Fig. S2; Fig. S3), separated unambiguously in clustering and PCA analyses (Fig. 2B; Fig. 3BC; Fig. S4; Fig. S6-S10), were strongly supported as separate lineages using SNAPP Bayes factors (Table 1) and BPP (Fig. S12), and passed the heuristic species delimitation criterion of gdi (Fig. 3A). A comparison of genetic and geographic distances moreover showed a clear distinction between intra-and interspecific genetic distances (Fig. 4). Finally, gene flow between these two lineages was estimated to have occurred at very low levels (G-PhoCS migration band = 0.08; Fig. 5c).

In contrast, separate species status of M. lehilahytsara and M. mittermeieri (Clade II) was not supported by comprehensive analysis of the genomic data. These species were paraphyletic in RAxML and SVDquartets analyses (Fig. 2A,C Although the Bayes factor support from SNAPP was strong by standard guidelines [START_REF] Kass | Bayes Factors[END_REF], the evidence was much weaker relative to species in Clade I and decreased when more individuals were included (Table S6). It is unsurprising that Bayes factors will support splitting lineages with genetic structure [START_REF] Sukumaran | Multispecies coalescent delimits structure, not species[END_REF][START_REF] Leaché | The spectre of too many species[END_REF] even with low levels of gene flow [START_REF] Barley | Impact of model violations on the inference of species boundaries under the multispecies coalescent[END_REF]. Therefore, standard guidelines for interpreting Bayes factors may be of limited value for delimiting species, as informed by the lack of monophyly, high levels of inferred gene flow, and failure of additional delimitation tests observed here. Guided delimitation also separated M. lehilahytsara and M. mittermeieri (Fig. S11), but similar criticisms of oversplitting (e.g. [START_REF] Barley | Impact of model violations on the inference of species boundaries under the multispecies coalescent[END_REF]) lead us to not interpret MSC delimitation results as evidence of species status. Most strikingly, reciprocal gdi statistics for Clade II were <0.2, thus falling in the range suggested to unambiguously indicate a single species [START_REF] Jackson | Species delimitation with gene flow[END_REF][START_REF] Leaché | The spectre of too many species[END_REF]; Fig. 3A). Finally, comparing genetic and geographic distances within Clade II showed that a single isolation-by-distance pattern fits both intra-and interspecific comparisons (Fig. 4). While the range of M. lehilahytsara expands considerably further south than the populations examined here, our results strongly suggest that M. mittermeieri and M. lehilahytsara are best considered a single species. Sampling gaps are expected to cause false positive species delimitations rather than false negatives [START_REF] Barley | Impact of model violations on the inference of species boundaries under the multispecies coalescent[END_REF][START_REF] Chambers | The multispecies coalescent over-splits species in the case of geographically widespread taxa[END_REF][START_REF] Mason | Coalescent-based species delimitation is sensitive to geographic sampling and isolation by distance[END_REF], thus indicating that additional sampling of M. lehilahytsara populations farther south should not affect our conclusions regarding genetic continuity. As such, and given that the original description of M. lehilahytsara precedes that of M. mittermeieri, we recommend that M. mittermeieri be synonymized as M. lehilahytsara.

Mitonuclear Discordance and Gene Flow

Mitonuclear discordance was observed for a subset of M. macarthurii individuals from Anjiahely. These individuals carried mtDNA similar to that of M. sp. #3 (see [START_REF] Radespiel | Exceptional diversity of mouse lemurs (Microcebus spp.) in the Makira region with the description of one new species[END_REF]) but had nDNA indistinguishable from sympatric M. macarthurii. Although genealogical discordance could be due to incomplete lineage sorting (e.g., [START_REF] Heckman | Multiple nuclear loci reveal patterns of incomplete lineage sorting and complex species history within western mouse lemurs (Microcebus)[END_REF][START_REF] Weisrock | Delimiting species without nuclear monophyly in Madagascar's mouse lemurs[END_REF], mitochondrial introgression is supported by D-statistics (Fig. S15) and the inferred low levels of gene flow from the northern M. sp. #3 population into M. macarthurii by G-PhoCS (Fig. S13). Besides a possible case in [START_REF] Sgarlata | Genetic and morphological diversity of mouse lemurs (Microcebus spp.) in northern Madagascar: The discovery of a putative new species?[END_REF], mitochondrial introgression has not previously been reported in mouse lemurs. We conclude that, somewhat curiously, the discovery of a divergent mtDNA lineage at Anjiahely [START_REF] Radespiel | Exceptional diversity of mouse lemurs (Microcebus spp.) in the Makira region with the description of one new species[END_REF]) which prompted the current work was apparently the result of mtDNA introgression from an undescribed species into its sister species.

Population Size and Species Delimitation

The comparison of effective population sizes in Clades I and II revealed marked differences, which can affect species delimitation tests such as gdi [START_REF] Leaché | The spectre of too many species[END_REF]. The gdi is calculated using population sizes and divergence times estimated under models with no gene flow, and since divergence time estimates in these models were highly similar in both clades (Fig. 5), differences in effective population sizes also appear to play a role in the stark difference in gdi. Indeed, gdi aims to quantify the probability that two sequences from the focal taxon coalesce more recently than the divergence time between the taxa, and larger effective population sizes result in slower sorting of ancestral polymorphisms [START_REF] Maddison | Gene trees in species trees[END_REF].

Assessing "progress" in speciation by quantifying rates of neutral coalescence, however, implies that the magnitude of genetic drift is a good predictor of species limits. At least when considering reproductive isolation (i.e., biological species), this can be problematic, given that the role of drift in speciation is generally thought to be small [START_REF] Rice | Laboratory experiments on speciation: what have we learned in 40 years?[END_REF]; Czekanski-Moir and Rundell 2019; but see [START_REF] Uyeda | Drift promotes speciation by sexual selection[END_REF]. Therefore, additional measures of divergence should be taken into account, including those that do not depend on population size, such as rates of gene flow and divergence time itself [START_REF] Yang | Bayesian species delimitation using multilocus sequence data[END_REF][START_REF] Leaché | The spectre of too many species[END_REF].

Sympatric Occurrence and the Tempo of Speciation in Mouse Lemurs

Sympatric Microcebus species were found at two study sites, with a representative of each of the two focal clades in Anjiahely (M. macarthurii and M. mittermeieri) as well as in Ambavala (M. sp. #3 and M. lehilahytsara; Fig. 1). These cases of sympatric occurrence, with no evidence for recent admixture, imply that the two clades are reproductively isolated. Though our methods cannot address the mechanisms underlying reproductive isolation, possible barriers include male advertisement calls, which tend to differ strongly among species [START_REF] Braune | Acoustic divergence in the communication of cryptic species of nocturnal primates (Microcebus ssp.)[END_REF], and timing of reproduction, which has previously been found to differ among sympatric mouse lemur species [START_REF] Schmelting | Reproduction of two sympatric mouse lemur species (Microcebus murinus and M. ravelobensis) in northwest Madagascar: first results of a long term study[END_REF][START_REF] Evasoa | Variation in reproduction of the smallest-bodied primate radiation, the mouse lemurs (Microcebus spp.): A synopsis[END_REF] including the focal species (Schüßler et al., in revision). Only six other cases of sympatry among mouse lemur species are known, five of which include M. murinus as one of the co-occurring species [START_REF] Radespiel | Can behavioral ecology help to understand the divergent geographic range sizes of mouse lemurs?[END_REF][START_REF] Sgarlata | Genetic and morphological diversity of mouse lemurs (Microcebus spp.) in northern Madagascar: The discovery of a putative new species?[END_REF].

Given that the sympatrically occurring species were estimated to have had a common ancestor as recently as ~700-800 kya (i.e., the divergence time between Clade I and Clade II, see Fig. 5), this suggests rapid evolution of reproductive isolation and a short time to sympatry among mouse lemurs. By comparison, [START_REF] Pigot | Dispersal and the transition to sympatry in vertebrates[END_REF] estimated that after 5 Ma of divergence, only 21-23% of primate species have attained sympatry. In fact, the one sympatric pair within their dataset of 74 sister species pairs younger than 2.5 million years consisted of Galago gallarum and G. senegalensis [START_REF] Pigot | Dispersal and the transition to sympatry in vertebrates[END_REF], which are also Strepsirrhini. Moreover, [START_REF] Curnoe | Timing and tempo of primate speciation[END_REF] compiled data for naturally hybridizing primate species, and found the median estimated divergence time to be 2.9 Ma. More broadly, primate speciation rates do not appear to be lower than those for other mammals or even vertebrates [START_REF] Curnoe | Timing and tempo of primate speciation[END_REF][START_REF] Upham | Inferring the mammal tree: Species-level sets of phylogenies for questions in ecology, evolution, and conservation[END_REF]. It should be noted, however, that the temporal estimates reported in our study are based on MSC analyses using mutation rates estimated from pedigree studies, whereas dates for other primate clades were largely calculated from fossil-calibrated relaxed clock methods.

Complexities of Divergence Time Estimates

CRYPTIC SPECIATION IN MOUSE LEMURS 29

There are two noteworthy discrepancies in divergence time estimates highlighted by this study. First, the age estimate between the M. mittermeieri and M. lehilahytsara lineages increased from approximately 100 kya (Fig. 5b) to more than 500 kya (Fig. 5c) when the MSC model allowed for gene flow. The substantial effect of incorporating or disregarding gene flow on divergence time estimation has been previously noted [START_REF] Leaché | Species delimitation using genome-wide SNP data[END_REF][START_REF] Tseng | Influence of gene flow on divergence dating -implications for the speciation history of Takydromus grass lizards[END_REF] and we here reiterate its significance. Second, the coalescent-based estimates of divergence times presented here differ drastically from estimates based on fossil-calibrated relaxed-clock methods. In the present study, we estimated the mean age of the most recent common ancestor (MRCA) of mouse lemurs to be under 1.5 Ma, with the highest upper bound of 95% HPDs across models at 2.40 Ma. This age estimate is in stark contrast to previous fossilcalibrated estimates of 8 -10 Ma [START_REF] Yang | Comparison of likelihood and Bayesian methods for estimating divergence times using multiple gene Loci and calibration points, with application to a radiation of cutelooking mouse lemur species[END_REF][START_REF] Dos Reis | Using phylogenomic data to explore the effects of relaxed clocks and calibration strategies on divergence time estimation: primates as a test case[END_REF].

Several factors likely contribute to this large difference. First, the MSC estimate uses a de novo mutation rate sampled from a distribution based on available pedigree-based mutation rates in primates, including mouse lemurs [START_REF] Campbell | Pedigree-based measurement of the de novo mutation rate in the gray mouse lemur reveals a high mutation rate, few mutations in CpG sites, and a weak sex bias[END_REF]. This rate is nearly two-fold higher than the estimated substitution rate for M. murinus [START_REF] Dos Reis | Using phylogenomic data to explore the effects of relaxed clocks and calibration strategies on divergence time estimation: primates as a test case[END_REF]. Second, converting coalescent units to absolute time also requires a generation time estimate. We attempted to account for uncertainty in generation time by similarly drawing from a distribution based on empirical parent age estimates [START_REF] Zohdy | Teeth, sex, and testosterone: aging in the world's smallest primate[END_REF][START_REF] Radespiel | An empirical estimate of the generation time of mouse lemurs[END_REF] in mouse lemurs.

Thus, either overestimation of the mutation rate and/or underestimation of the generation time would lead to divergence time estimates that are too recent. However, theoretical considerations suggest that instead, mouse lemur divergence time estimates from fossil-calibrated clock models are too old.

When incomplete lineage sorting is common, clock models that assume a single topology underlies all loci can overestimate species divergences compared to MSC estimates that allow gene trees to vary [START_REF] Stange | Bayesian divergence-time estimation with genome-wide single-nucleotide polymorphism data of sea catfishes (Ariidae) supports Miocene closure of the Panamanian Isthmus[END_REF][START_REF] Feng | Estimating uncertainty in divergence times among three-spined stickleback clades using the multispecies coalescent[END_REF]. This is likely to apply to mouse lemurs given that high levels of incomplete lineage sorting have been previously documented [START_REF] Heckman | Multiple nuclear loci reveal patterns of incomplete lineage sorting and complex species history within western mouse lemurs (Microcebus)[END_REF][START_REF] Weisrock | Delimiting species without nuclear monophyly in Madagascar's mouse lemurs[END_REF]Hotaling et al., 2016). Moreover, due to the absence of a post-K-Pg terrestrial fossil record for Madagascar, clock-model estimates of divergence times in mouse lemurs have relied on fossil calibrations from the distantly related African sister lineage of lemurs, the Lorisiformes [START_REF] Seiffert | Fossil evidence for an ancient divergence of lorises and galagos[END_REF], as well as from anthropoid primates and other mammals. This scenario -estimation of divergence times for younger, internal nodes with calibrations placed on much older nodes -should lead to overestimation of divergence times [START_REF] Angelis | The impact of ancestral population size and incomplete lineage sorting on Bayesian estimation of species divergence times[END_REF]. Therefore, it is likely that divergence times between mouse lemur species have been overestimated by previous studies with fossil-calibrated clock models (e.g. [START_REF] Yang | Comparison of likelihood and Bayesian methods for estimating divergence times using multiple gene Loci and calibration points, with application to a radiation of cutelooking mouse lemur species[END_REF][START_REF] Dos Reis | Using phylogenomic data to explore the effects of relaxed clocks and calibration strategies on divergence time estimation: primates as a test case[END_REF], and we suggest that the mutation rate-calibrated MSC divergence times presented here are more accurate.

Our estimates of divergence times imply that the entire mouse lemur radiation originated in the Pleistocene, in turn suggesting that Pleistocene climatic oscillations represent a likely factor leading to geographic isolation and subsequent genetic divergence. Periods of drought during glacial maxima are hypothesized to have caused dramatic contraction of forest habitats [START_REF] Burney | Environmental change, extinction and human activity: evidence from caves in NW Madagascar[END_REF][START_REF] Gasse | Late Quaternary environmental changes from a pollen and diatom record in the southern tropics (Lake Tritrivakely, Madagascar)[END_REF][START_REF] Wilmé | Biogeographic evolution of Madagascar's microendemic biota[END_REF][START_REF] Kiage | Late Quaternary paleoenvironmental changes in East Africa: a review of multiproxy evidence from palynology, lake sediments, and associated records[END_REF] and to isolation of previously connected populations. Notably, the patterns of differentiation observed in this study are consistent with the predictions of [START_REF] Wilmé | Biogeographic evolution of Madagascar's microendemic biota[END_REF] wherein Madagascar's river drainage systems created high-elevation retreat-dispersal corridors during periods of climatic oscillation. That is, whereas the lineages in Clade I (highly differentiated and low N e ) appear to occur only in lowland forests, those in the Clade II (poorly differentiated and high N e ) occur at both higher and lower elevations (Schüßler et al., in revision). Moreover, the Mananara river runs between the fairly distinct northern and southern populations of M. sp. #3, further emphasizing the potential of large rivers to act as phylogeographic barriers in lemurs [START_REF] Martin | Adaptive radiation and behaviour of Malagasy lemurs[END_REF][START_REF] Pastorini | A molecular approach to comparative phylogeography of extant Malagasy lemurs[END_REF][START_REF] Goodman | Biogeography of lemurs in the humid forests of Madagascar: The role of elevational distribution and rivers[END_REF][START_REF] Olivieri | The ever-increasing diversity in mouse lemurs: Three new species in north and northwestern Madagascar[END_REF].

Population Size Dynamics

A long-term decline in population size was inferred for the lineage leading to M. sp. #3.

While changes in inferred N e may be confounded by changes in population structureespecially for single-population sequential Markovian coalescent (PSMC/MSMC) models that do not explicitly consider population subdivision [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution--lessons for ancestral population size inference?[END_REF][START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF])we recovered CRYPTIC SPECIATION IN MOUSE LEMURS 32 similar results with both MSMC and G-PhoCS analyses (Fig. 6A). This congruence is especially persuasive given the underlying differences between the G-PhoCS and MSMC models and their input data. Moreover, Markovian coalescent approaches are robust to genome assembly quality [START_REF] Patton | Contemporary demographic reconstruction methods are robust to genome assembly quality: A case study in Tasmanian Devils[END_REF], yielding further confidence in the results. The inferred decline and population subdivision of M. sp. #3 was initiated long before anthropogenic land use, supporting the emerging consensus that human colonization in Madagascar alone does not explain the occurrence of open habitats and isolated forest fragments [START_REF] Quéméré | Genetic data suggest a natural prehuman origin of open habitats in northern Madagascar and question the deforestation narrative in this region[END_REF][START_REF] Vorontsova | Madagascar's grasses and grasslands: anthropogenic or natural?[END_REF][START_REF] Yoder | Geogenetic patterns in mouse lemurs (genus Microcebus) reveal the ghosts of Madagascar's forests past[END_REF][START_REF] Salmona | Climate change and human colonization triggered habitat loss and fragmentation in Madagascar[END_REF][START_REF] Salmona | Late Miocene origin and recent population collapse of the Malagasy savanna olive tree (Noronhia lowryi)[END_REF][START_REF] Hackel | Grass diversification in Madagascar: In situ radiation of two large C3 shade clades and support for a Miocene to Pliocene origin of C4 grassy biomes[END_REF]. Conversely, results for the M. mittermeieri lineage do not indicate a declining N e through time (Fig. 6B). This latter result may well be a simple corollary of the evidence described above, that this lineage is part of a single species complex represented by Clade II and thus occurs at both higher and lower elevations in northeastern Madagascar.

Conclusions

We have shown that substantial mouse lemur diversity exists within a 130-km-wide stretch in northeastern Madagascar, including two instances of sympatric occurrence between representatives of two closely related clades. Within one of these clades, our comprehensive approach that uses a variety of genomic analyses indicates that the undescribed lineage M. sp. #3 represents a distinct species, while the two named species in the other clade, M. mittermeieri and M. lehilahytsara, are better considered a single, widespread species with significant population structure. Given that the original description of M. lehilahytsara precedes that of M. mittermeieri, primate taxonomists should synonymize the two as M. lehilahytsara.

The divergence times calculated here using pedigree-based mutation rate estimates with the MSC are much younger than those of previous studies that used external fossil-based calibrations with relaxed-clock methods. The younger dates suggest rapid evolution of reproductive isolation in mouse lemurs as well as a Pleistocene origin of the radiation, likely following population isolation due to climatic oscillations. This departure from previous hypotheses of mouse lemur antiquity emphasizes the need for future studies focused on resolving discrepancies in divergence time estimates, both in mouse lemurs and in other recently evolved organismal groups for which such comparisons have yet to be made.

  guided species delimitation analyses with BPP (Yang and Rannala 2010; Rannala and Yang CRYPTIC SPECIATION IN MOUSE LEMURS 14 2013) using full-length fasta files for a subset of individuals based on the species tree estimated by SVDquartets and SNAPP.

  Interspecific gene flow. -G-PhoCS inferred high levels of gene flow in Clade II, from M. mittermeieri to M. lehilahytsara [population migration rate (2Nm) = 1.59 (95% HPD: 1.50-1.68), migrants per generation: 0.18% (95% HPD: 0.09-0.27%)], and much lower levels of gene flow in CRYPTIC SPECIATION IN MOUSE LEMURS 21

  ; Fig. S2) and not as clearly separated in clustering and PCA analyses (Fig. 2B; Fig. 3BD; Fig. S5, S7-S9; Fig. S10).

  Biology at The University of Western Australia to AM, the Fundação para a Ciência e a Tecnologia, Portugal (PTDC/BIA-BEC/100176/2008, PTDC/BIA-BIC/4476/2012, and SFRH/BD/64875/2009), the Groupement de Recherche International (GDRI) Biodiversité et développement durable -Madagascar, the Laboratoire d'Excellence (LABEX) TULIP (ANR-10-LABX-41) and CEBA (ANR-10-LABX-25-01, the Instituto Gulbenkian de Ciência, Portugal to LC and JS, the ERA-NET BiodivERsA project: INFRAGECO (Inference, Fragmentation, Genomics, and Conservation, ANR-16-EBI3-0014 & FCT-Biodiversa/0003/2015) the LIA BEEG-B (Laboratoire International Associé -Bioinformatics, Ecology, Evolution, Genomics and Behaviour, CNRS) to LC and JS. Further financial support came from the Institute of Zoology, University of Veterinary Medicine Hannover and UR acknowledges the long-term support of the late Elke Zimmermann for her research activities on Madagascar. The genomic data were generated with funds from NSF DEB-1354610 to ADY and DWW and from the EDB Lab to JS. ADY also gratefully acknowledges support from the John Simon Guggenheim Memorial Foundation, the Alexander von Humboldt Foundation, and the Duke Tropical Conservation Initiative. EELJ would like to acknowledge support from the Ahmanson Foundation for the data generation. This work was performed in collaboration with the GeT core facility, Toulouse, France (http://get.genotoul.fr), and was supported by France Génomique National infrastructure, funded as part of "Investissement d'avenir" program managed by Agence Nationale pour la Recherche (contract ANR-10-INBS-09). JS, UR & LC also gratefully CRYPTIC SPECIATION IN MOUSE LEMURS 36 acknowledge support from the Get-Plage sequencing and Genotoul bioinformatics 614 (BioinfoGenotoul) platforms Toulouse Midi-Pyrenees. This is DLC publication #1459.

ACKNOWLEDGMENTS

We thank the editors and three anonymous reviewers for their constructive input which substantially improved our manuscript. The study was conducted under the research permit No. 197/17/MEEF/SG/DGF/DSAP/SCB.Re (DS), 072/15/MEEMF/SG/DGF/DCB.SAP/SCB (MBB), 137/13/MEF/SG/DGF/DCB.SAP/SCB (DWR), 175/14/MEF/SG/DGF/DCB.SAP/SCB (A . Miller), kindly issued by the directeur du système des aires protégées, Antananarivo and the regional authorities (Direction Régional de l'Environnement, de l'Ecologie et de Forêts). We are endebted to J.H. Ratsimbazafy, N.V. Andriaholinirina, C. Misandeau, B. Le Pors and S.

Rasoloharijaona, for their help with administrative tasks, to A. Miller for sharing samples, and to G. Besnard for facilitating this study. We thank our field assistants (T. Ralantoharijaona, I.

Sitrakarivo, C. Hanitriniaina and T. Ralantoharijaona), the Wildlife Conservation Society

Madagascar and the ADAFAM (Association Des Amis de la Forêt d'Ambodiriana-Manompana, C. Misandeau in particular) for their valuable help during sample collection. We warmly thank the many local guides and cooks for sharing their incomparable expertise and help in the field, misaotra anareo jiaby.

Funding was granted by the Bauer Foundation and the Zempelin Foundation of the "Deutsches Stiftungszentrum" under grant no. T237/22985/2012/kg and T0214/32083/2018/sm to DS, Duke Lemur Center/SAVA Conservation research funds to MBB, the School of Animal

Statement

-Raw sequencing data is available at NCBI BioProject PRJNA560399.

-All other data is available from the Dryad Digital Repository at https://datadryad.org/stash/share/-q8Ez_lXqi7LiDsBHJTAgM1JpdPOCGep0nU5gZ_EH9k