Introduction
Species are cornerstone in evolution and are used as study units in both evolutionary and conservation biology (Butlin et al. 2012, Wake 2006). Though divergence with gene-flow has been theorized and observed (review in Pinho & Hey 2010), the model of allopatric speciation predominates by far in the literature (Stroud & Losos 2016, Butlin et al. 2012). In this model, a physical barrier to gene flow catalyzes genetic differentiation between populations, through selection and/or genetic drift, eventually followed by other pre- or post-zygotic barriers (Coyne & Orr 2004). In practice however, the mechanisms that impede gene flow and promote differentiation are multifactorial and still poorly understood (Butlin et al. 2012, Price 2008). In particular, geographic barriers alone may not explain the differentiation of populations in highly dispersive species, e.g. marine birds (e.g. Genovart et al. 2007), birds of prey (e.g. Doyle et al. 2016), mammals (Hassanin et al. 2018) or plants (Sanz et al. 2014). Seabirds are a case in point: their wide geographic distribution and dispersal ability should theoretically maintain high levels of gene flow, but many seabirds show surprisingly strong geographic population structure, a pattern attributed to their high degree of philopatry (Abbott & Double 2003, Smith et al. 2007a).
Present or historic landmasses were identified as the most important barriers to gene flow in seabirds (Friesen 2015, Friesen et al. 2007a) though they cannot explain every differentiation patterns, and other factors such as philopatry, isolation-by-distance, or segregation between breeding and non-breeding zones may play a role. So far, however, no study has contrasted all these processes simultaneously. Moreover, most previous studies were based on maternally-inherited mitochondrial markers (Friesen 2015), while gene flow is expected to vary between sexes given stronger male philopatry in seabirds and other birds in general (Greenwood 1980). Sex-specific patterns of divergence are expected, and indeed mito-nuclear discordance has been detected in seabirds, with more genetic structure in mt than in nuDNA (Burg & Croxall 2001, Deane 2013, Gangloff et al. 2013, Silva et al. 2015, Welch et al. 2011, Friesen et al. 2006, Welch Fleicher 2012, Genovart 2012, Sonsthagen et al. 2016; but see Pons et al. 2014). Such discordance was suggested to result from incomplete lineage sorting in nuDNA due to a higher effective sample size than mtDNA (McKay & Zink 2010), but other mechanisms were proposed such as adaptive introgression of mtDNA, demographic disparities and sex-biased asymmetries (review in Toews & Brelsford 2012). Disentangling these latter processes can be difficult (McKay & Zink 2010), thus a combined use of mitochondrial and biparentally-inherited nuclear markers, as well as coalescent theory-based methods using the lineage species concept coupled with population genetics within a phylogeographic context, are advocated (Hudson & Coyne 2002).
Shearwaters (order Procellariiformes) are long-lived birds showing slow demographic rates (Warham 1996). They breed in large colonies on remote oceanic islands, are pelagic (González-Solís et al. 2007) and highly philopatric (Brooke 2004). Geographic structuring of populations can be strong (Genovart et al. 2012, Gómez-Díaz et al. 2009), and the systematics of some taxa is still highly controversial, especially for the little and Audubon’s Shearwaters Puffinus assimilis-lherminieri complex (Austin et al. 2004). This widespread small-sized, black-and-white shearwater breeds from equatorial to sub-arctic seas (see map in Fig. 1). Previous studies, based only on the mitochondrial gene cytb , indicated that population genetic clusters matched with geographic distribution (Austin et al. 2004, Kawakami et al. 2018). Using this single marker, three distinct lineages were recognized in the North Atlantic: lherminieri in the Caribbean and off Brazil, baroli in the Azores, Canaries and Madeira, and boydi in Cape Verde (Fig. 1). These taxa are characterised by non-overlapping breeding and non-breeding distributions at sea (Ramos et al. in review), and are thus geographically and genetically (for at least one marker) separated. Still, they are morphologically and ecologically highly similar (Precheur et al. 2016, Calabrese et al. in review), a pattern typical of the first stages of the speciation process, the so-called grey-zone (De Queiroz 2007). Unsurprisingly, their taxonomic ranking has been hotly debated, e.g. baroli being considered as belonging to assimilis (Shirihai et al. 1995), lherminieri (Austin et al. 2004) or a species of its own (Sangster et al. 2005). Lineages belonging to this complex are also found in the Indian and Pacific Ocean, with populations breeding in the Seychelles (nicolae), Réunion (bailloni) and many islands in the Pacific Ocean (dichrous). Breeding populations of bailloni are characterised by different breeding phenologies on the northern and southern parts of Réunion (Bretagnolle & Attie 1996), potentially impacting genetic structuration (Friesen et al. 2007b). Indian Ocean birds were alternatively considered as a P. lherminieri subspecies (Warham 1990) or subspecies of a bailloni pan-tropical taxon (Austin et al. 2004). The exact taxonomic status of these six lineages is thus largely unresolved.
Here we consider these six lineages, covering Atlantic and the Indo-Pacific branches of the Puffinus assimilis-lhermineri complex. We sampled birds from all but one known breeding localities in the North Atlantic (four Caribbean breeding sites for lherminieri, two main breeding sites in Cape Verde for boydi, one breeding site in Azores, Madeira, Selvagens Canary Islands for baroli) and in the northern Indian Ocean (northern and southern populations of bailloni in Réunion and Seychelles for nicolae; Fig. 1). Only a single breeding locality was unsampled, the Fernando de Noronha archipelago, off Brazil (Fig. 1). We analysed three mitochondrial and six nuclear markers to delineate genetic units and investigated patterns of genetic differentiation and divergence among populations. As conflicting geographic patterns between mitochondrial and nuclear markers (hereon referred to as mito-nuclear discordance; Toews & Brelsford 2012) was systematically found in petrels so far (references above), we inferred population structure among females and males independently to test for sex-biased dispersal on nuDNA. We then used multispecies coalescent inference and an ABC framework to investigate evolutionary scenarios of breeding site colonisation over the last million years, trying to disentangle contrasted processes shaping evolutionary history and the contemporary population structure of this complex, such as landmass presence, isolation-by-distance, oceanographic conditions (sea surface temperature) and climatic oscillations.
Material and Methods