\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,romanian,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Coalescence times, life history traits and conservation concerns: an
example from four coastal shark species from the Indo-Pacific}



\author[1]{Pierre Lesturgie}%
\author[2]{Serge Planes}%
\author[1]{Stefano Mona}%
\affil[1]{Museum National d'Histoire Naturelle}%
\affil[2]{Université de Perpignan}%


\vspace{-1em}



  
  \date{March 14, 2025}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Dispersal abilities play a crucial role in shaping the extent of
population genetic structure, with more mobile species being panmictic
over large geographic ranges and less mobile ones organized in
meta-populations exchanging migrants to different degrees. In turn,
population structure directly influences the coalescence pattern of the
sampled lineages, but the consequences on the estimated variation of the
effective population size (\emph{Ne}) over time obtained by means of
\emph{unstructured} demographic models remain poorly understood.
However, this knowledge is crucial for biologically interpreting the
observed Ne trajectory and further devising conservation strategies in
endangered species. Here we investigated the demographic history of four
shark species (\emph{Carharhinus melanopterus, Carharhinus limbatus,
Carharhinus amblyrhynchos, Galeocerdo cuvier}) with different degrees of
endangered status and life history traits related to dispersal
distributed in the Indo-Pacific and sampled off New Caledonia. We
compared several evolutionary scenarios representing both
\emph{structured} (meta-population) and~\emph{unstructured} models and
then inferred the~\emph{Ne} variation through time. By performing
extensive coalescent simulations, we provided a general framework
relating the underlying population structure and the observed~\emph{Ne}
dynamics. On this basis, we concluded that the recent decline observed
in three out of the four considered species when assuming
\emph{unstructured} demographic models can be explained by the presence
of population structure. Furthermore, we also demonstrated the limits of
the inferences based on the sole site frequency spectrum and warn that
statistics based on linkage disequilibrium will be needed to exclude
recent demographic events affecting meta-populations.%
\end{abstract}%



\sloppy


\textbf{Introduction}

Reconstructing the evolutionary history of a species is a challenging
exercise only partially eased by the growing size of genetic data
available. Indeed, larger amounts of data will indeed provide more
precision but not more accuracy if the model(s) chosen to infer
demographic parameters is distant from the true one. Species are dynamic
entities whose geographic range has often changed in time through range
expansions, contractions and shifts (Arenas, Ray, Currat, \& Excoffier,
2012; Excoffier, Foll, \& Petit, 2009; Mona, Ray, Arenas, \& Excoffier,
2014). As a consequence, many species are most likely organized in
meta-populations (i.e. groups of demes or sub-populations exchanging
migrants to some extent), even though the more vagile ones might be
panmictic at a large scale (Corrigan et al., 2018; Karl et al., 2010).
Neglecting the meta-population structure (i.e., performing demographic
inferences under \emph{unstructured} models) may lead to spurious
inference of population size change (Chikhi, Sousa, Luisi, Goossens, \&
Beaumont, 2010; Maisano Delser et al., 2019, 2016; Olivier Mazet,
Rodr\selectlanguage{ngerman}íguez, \& Chikhi, 2015), which is particularly worrisome for species
of conservation concern. Unfortunately, the link between the inferred
temporal trajectory of \emph{Ne} and the real demographic history of the
meta-population remains largely under explored. However, the role of
connectivity, particularly the number of migrants \emph{Nm} exchanged
each generation and the migration matrix, has been put forward as a key
actor in shaping the gene genealogy of lineages sampled from a deme
belonging to a meta-population (Chikhi et al., 2010; Mona et al., 2014;
Ray, Currat, \& Excoffier, 2003; Städler, Haubold, Merino, Stephan, \&
Pfaffelhuber, 2009).

Understanding the relations between meta-population structure, the
inferred \emph{Ne} variation under \emph{unstructured} models, and
species dispersal abilities, is crucial to correctly interpret the
pattern of genetic variability and to establish conservation priorities.
To search for general rules describing such relations, we followed an
inductive approach investigating species: i) with large distribution
(which in principle should guarantee an organization in
meta-populations); ii) with different life history traits (LHT) related
to dispersal; iii) of conservation concerns. In this spirit, we selected
for our genomic study four shark species (\emph{Carcharhinus
amblyrhynchos} , \emph{Carcharhinus limbatus} , \emph{Carcharhinus
melanopterus,} and \emph{Galeocerdo cuvier} ) from New Caledonia. These
species have a large and overlapping distribution in the Indo-Pacific
(https://sharksrays.org/) and they differ for LHT features such as size
(which is positively correlated with the capacity for long distance
swimming and oceanic migration (Parsons, 1990)), residency pattern, and
long-distance dispersal ability as measured by tagging data (Table S1).
Moreover,~the IUCN red list reported that~the black-tip shark (\emph{C.
limbatus})~and~the tiger shark\emph{~}(\emph{G. cuvier})~are Near
Threatened (with a decreasing trend in the tiger shark),~the black-tip
reef shark (\emph{C. melanopterus})~is Vulnerable~with decreasing trend,
and~the grey reef shark (\emph{C. amblyrhynchos})\emph{~}is Endangered
with decreasing trend as well\emph{.~}We first compared several
population genetics models by means of coalescent simulations coupled
with an \emph{approximate Bayesian computation} framework (Bertorelle,
Benazzo, \& Mona, 2010) to detect whether panmixia or a meta-population
model best describe the genomic variation of each species. Then, we
inferred the demographic parameters under the most likely model and
applied the \emph{stairwayplot} , which assumes a panmictic unstructured
population (Liu \& Fu, 2015), to detect the\emph{Ne} variation through
time in each species. We finally run extensive coalescent simulations
under the tested meta-population models with parameters compatible to
those observed in real data. The simulated datasets were in turn
analysed with the \emph{stairwayplot} to: i) help interpreting the
observed data in the four shark species; ii) providing general
coalescent arguments relating the demographic history of a
meta-population and the reconstructed variation in \emph{Ne} through
time by means of \emph{unstructured} models.

\textbf{Material \& Methods}

\textbf{Sampling}

Eight specimens of tiger shark (\emph{G. cuvier} ), 13 black tip shark
(\emph{C. limbatus} ), and 12 grey reef shark (\emph{C. amblyrhynchos} )
were collected off New Caledonia. Total genomic DNA was extracted from
muscle tissue or fin clips, and preserved in 96\% ethanol using QIAGEN
DNeasy Blood and Tissue purification kit (Qiagen, Hilden, Germany)
according to the manufacturer's protocols. Double-digest
restriction-associated DNA (ddRAD) libraries were prepared following
Peterson et al. (2012) using EcoRI and MspI restriction enzymes and a
400-bp size selection. The genomic libraries obtained were sequenced
with a HiSeq 2500 Illumina sequencer (single-end, 125~bp). Exon capture
data of eight black tip reef shark (\emph{Carcharhnius melanopterus} )
from New Caledonia (Maisano Delser et al. 2019) were included in this
study for comparative purposes.

\emph{\textbf{De novo}} \textbf{assembly and data filtering (dd-RADseq
samples)}

Raw reads were first demultiplexed and quality filtered through
the~\emph{process\_radtags.pl} pipeline in stacks v.2.5 (Rochette,
Rivera-Colón, \& Catchen, 2019). In the absence of a reference genome
for any of the three species, RAD-seq loci were~\emph{de novo} assembled
independently in each species under the~\emph{denovo\_map.pl} pipeline
in stacks. We used the following assembly parameters:\emph{m} =3
(minimum read depth to create a stack),~\emph{M}=4 (number of mismatches
allowed between loci within individuals), and~\emph{n=4}(number of
mismatches allowed between loci within catalogue). We found an average
coverage per species of \textasciitilde{}10x (see results). A consensus
on the threshold below which SNP calling may be considered unreliable is
still lacking. However, genotype free estimation of allele frequency is
generally recommended with low to medium coverage (Korneliussen,
Albrechtsen, \& Nielsen, 2014). This approach, implemented in the
software angsd v.0.923 (Korneliussen et al., 2014), has been rarely
applied to Rad-seq data (however, see Warmuth and Ellegren (2019) for an
exception) and, to our knowledge, never to Rad-seq data from non-model
organisms, probably due to need of a reference sequence for the software
to work. Here, we followed the approach of Heller et al. (2021) and
Khimoun et al. (2020) by creating an artificial reference sequence.
First, we used the~\emph{population} script in stacks to assemble loci
present in at least 80\% of the individuals (using the flag~\emph{r=0.8}
); then, we concatenated the consensus sequences of the retrieved loci
spaced by a stretch of 120 N (unknown) characters (the same length of
the Rad-loci) to facilitate the subsequent mapping. Raw reads were then
mapped back to the novel reference sequence by means of
the~\emph{bwa-mem} algorithm with default parameters (Li \& Durbin,
2009). Using custom bash scripts coupled with angsd, we applied a number
of filters to the aligned data and eliminated:~\emph{i)} sites with
coverage \textless{}3 (-\emph{minIndDepth} =\emph{3} flag),\emph{ii)}
bad quality bases and poorly aligning reads (-\emph{minQ} and
-\emph{minMapQ} and~\emph{-C} flags with default
values);~\emph{iii)}poor quality sites based on the per base alignment
quality (-\emph{baq=1} flag); :~\emph{iv)} SNPs in the last 5~bp of each
locus;\emph{v)} SNPs heterozygote in at least 80\% of
individuals;~\emph{vi)}loci with more than 5 SNPs that could potentially
be paralogous;\emph{vii} ) sites with missing data by setting
the~\emph{-minInd} flag to the total number of individuals retained in
each species. The filtered dataset was then used to generate a site
allele frequency likelihood file, with the genotype likelihoods computed
with the SAMtools method (\emph{-GL=1} flag), further optimised to
compute a folded~\emph{site frequency spectrum} (SFS) with no missing
data for downstream analyses. An alternative (and simpler) approach
would have been to augment~\emph{m} to achieve an higher coverage
(Paris, Stevens, \& Catchen, 2017). However, beside the considerable
loss in the number of assembled loci (and hence of retrieved SNPs), we
found by extensive simulation of~\emph{in silico} Rad experiments that
selecting high coverage loci biases the SFS towards low frequency
variants (Mona, Bertorelle, Benazzo, in preparation). The SFS
for~\emph{C. melanopterus}was estimated directly from the high coverage
exon-capture dataset of Maisano Delser et al. (2019).

\textbf{Genetic diversity and demographic inferences}

Nucleotide diversity (\selectlanguage{romanian}\emph{Ɵ\textsubscript{π}}\selectlanguage{greek}
\selectlanguage{english}),\selectlanguage{romanian}\emph{Ɵ\textsubscript{w}} (Watterson's theta, based on segregating
sites (Watterson, 1975)) and Tajima's \emph{D (TD,} (Tajima, 1989)) were
computed from the SFS for each species with custom scripts. Significance
of \emph{TD} was evaluated after 1,000 coalescent simulations of a
constant population model with scaled size \selectlanguage{greek}\emph{Ɵ\textsubscript{π}} .
\selectlanguage{english}To test whether sampled demes are isolated or belong to a structured
meta-population and to eventually estimate connectivity, we devised
three alternative evolutionary models for each species (Figure 1) within
an \emph{approximate bayesian computation} framework (ABC) framework.
Model NS (non-structured) defined an isolated population characterized
by a modern effective population size (\emph{Nmod} ) switching
instantaneously into an ancestral population size (\emph{Nanc} )
at\emph{Tc} generations before present. Model FIM specifies a
non-equilibrium finite island model defined by 100 demes
exchanging\emph{Nm} migrants each generation under a symmetric migration
matrix. The array of demes is instantaneously colonized \emph{Tcol}
generations before present from a population with an ancestral size
(\emph{Nanc} ). Model SST is similar to FIM but demes exchanging
migrants only with their four neighbours (or less, if they are at the
border of the array), in a steppingstone fashion. We performed 50,000
coalescent simulations from prior distributions using fastsimcoal
v.2.6.0.3 (Excoffier, Dupanloup, Huerta-S\selectlanguage{ngerman}ánchez, Sousa, \& Foll,
2013)\emph{,}reproducing the exact number of individuals and loci for
each species (Table 1). We first performed model selection through the
random forest (RF) classification method implemented in the abcRF R
package (Pudlo et al., 2016). We then performed 50,000 additional
simulations under the most supported model in order to estimate
demographic parameters with the abcRF regression method (Raynal et al.,
2019). Both model selection and parameter estimation were computed with
the following set of summary statistics: the SFS,
\selectlanguage{romanian}\emph{Ɵ\textsubscript{π}}\selectlanguage{greek} ,\selectlanguage{english}\selectlanguage{romanian}\emph{Ɵ\textsubscript{w}} and \emph{TD.} The
first two axes of a Linear Discriminate Analysis performed on the
previous statistics were also included for model selection in order to
increase the accuracy of the estimates (Pudlo et al., 2016). Even though
\selectlanguage{greek}\emph{Ɵ\textsubscript{π}} ,\selectlanguage{english}\selectlanguage{romanian}\emph{Ɵ\textsubscript{w}} and \emph{TD} are
function of the SFS, they convey additional information by the
non-linear feature of the functions. Information redundancy among the
considered summary statistics is accounted for by the RF algorithm.
Model selection and parameter estimation were~run~twice~on each set of
simulations~to check the consistency of the analyses, and cross
validation (or confusion matrix for the model selection) was performed
on the first of the two runs. The number of trees in each RF algorithm
was chosen by monitoring the evolution of the out-of-bag error (Pudlo et
al., 2016).

We investigated the variation in the effective population size
(\emph{Ne} ) through time by running the composite likelihood approach
implemented in the stairwayplot v.0.2 software (Liu \& Fu, 2015). We set
the generation time to seven years for \emph{C. melanopterus} (Maisano
Delser et al., 2016) and to 10 years for the other species (Cort\selectlanguage{ngerman}és,
2002; Pirog et al., 2019) for all demographic inferences. We applied a
mutation rate per generation per site of 8.4×10\textsuperscript{-9} to
the exon capture data of \emph{C. melanopterus} (Maisano Delser et al.
2016) and of 1.93×10\textsuperscript{-8} to the RADseq data for the
remaining three species. This mutation rate was determined by scaling
genetic diversity between ddRAD (obtained under the same protocol of
this study) and Exon Capture data from 12 \emph{C. melanopterus}
individuals from Moorea, French Polynesia (Supplementary Material).

\textbf{Simulation study}

We ran coalescent simulations under FIM, SST and their modified version
FIM-CH and SST-CH, where the~\emph{Nm}parameter is changed
at~\emph{T}\textsubscript{\emph{CH}}~generations B.P. (Figure 1), to
first inspect the shape of the SFS and to further uncover the variation
of~\emph{Ne}~over time assuming a panmictic population by means of
the~\emph{stairwayplot}. We investigated in total 288 demographic
scenarios under the four meta-population models (Tables 2, S2, S3, S4,
S5 and S6).~Similarly to the analyses performed on the real data, all
scenarios were represented by~\emph{d}=100~demes exchanging migrants.~We
sampled 10 diploid individuals either from a randomly selected deme in
the case of FIM/FIM-CH (since all demes have the same
coalescencehistory) or from the central deme of the array in the case of
SST/SST-CH (to avoid border effects).~Deme size was fixed
to~\emph{N}\textsubscript{\emph{DEME}}~=5000 with~\emph{m}~varying
accordingly to obtain a long-term~\emph{Nm}~of 1, 5, 10, and 15 in order
to encompass the range of the estimated values (see
results).~\emph{T}\textsubscript{\emph{COL}}~was fixed to 5,000, 15,000
and 50,000 generations B.P. or to \selectlanguage{english}[?] (i.e., equilibrium model), and the
ancestral effective size was fixed
to~\emph{N}\textsubscript{\emph{ANC}}~= 50,000.~\emph{~}Change of
connectivity occurred at~\emph{T}\textsubscript{\emph{CH}}=10 or 50
generations B.P., to mimic human induced effects due to overfishing
and/or habitat modifications (i.e., climate changes). Looking forward in
time, we modelled the change in connectivity by instantaneously
decreasing~\emph{m~}or~\emph{N}\textsubscript{\emph{DEME}}{}by a factor
10 or 100 with respect to the long-term~\emph{Nm~}(Tables S3, S4, S5 and
S6). For each combination of parameters, we performed 100 coalescent
simulations of 50,000 Rad-like loci of 115 bp. Mutation rate per site
per generations was set to 1.93x10\textsuperscript{-8}~and the
generation time to 10 years. We computed for each scenario (averaged
over the 100 replicates): a) summary statistics
(\selectlanguage{greek}θ\selectlanguage{english}\selectlanguage{greek}\textsubscript{\emph{π}}, θ\selectlanguage{english}\textsubscript{\emph{w}}, and~\emph{TD});
b) the normalised SFS as in~(Lapierre, Lambert, \& Achaz, 2017); c)
the~\emph{stairwayplot}, to reconstruct the apparent variation
of~\emph{Ne}~through time.~~We note that the number of diploid
individuals and simulated loci were chosen to be consistent with our
data (preliminary analyses conducted on a subsample of 5,000 loci
produced consistent results).~\emph{}

\textbf{Results}

Summary statistics (number of assembled loci, SNPs, genetic diversity
and Tajima's~\emph{D}) are presented in Table 1. Mean coverage (and
standard deviation) per sample was 9.02 (\selectlanguage{ngerman}± 2.62), 7.93 (± 0.48), 8.39 (±
0.81) for~\emph{G. cuvier, C. limbatus}~and~\emph{C.
amblyrhynchos~}respectively.

We compared the models NS, FIM, and SST (Figure 1) in the four species
by means of an ABC-RF algorithm and estimated demographic parameters for
the most supported model. After checking for the evolution of the
out-of-bag error of the RF, model selection and parameter estimation
were computed using respectively 500 and 1,000 trees in each species. We
found that NS had the higher posterior probability (\emph{p}=0.84)
for~\emph{G. cuvier}~(Tables 1 and S7). In contrast, demographic
histories~of the~three other species~were best described by SST, with a
posterior probability ranging from 0.53 to 0.88 (Tables 1~and S7). The
estimated median number of migrants per generation~\emph{Nm}~was 1.8
(95\% CI: 0.7-3.0) for~\emph{C. melanopterus}, 6.6 (95\% CI: 1.5-15.4)
for~\emph{C. limbatus,}~and 11.5 (95\% CI: 3.0-22.0) for~\emph{C.
amblyrhynchos~}(Figure 2, Table 1). The posterior distribution
of~\emph{Nm}~strongly differed from the prior distribution and showed a
clear unimodal peak with small credible intervals, and low mean square
error (SME) and mean root square error (SMRE) in all three species
(Figure 2, Table S8, suggesting that these estimates are highly
reliable. Conversely,
both~\emph{T}\textsubscript{\emph{COL}}~and~\emph{N}\textsubscript{\emph{ANC}}~had
larger SME and SMRE errors in all species (Table S8), but it was only
in~\emph{C. melanopterus~}where posterior and prior distribution could
not be distinguished (Figure 2).~\emph{T}\textsubscript{\emph{COL}}~has
a clear unimodal distribution in~\emph{C. amblyrhynchos}~but a more
disperse one (and with wider credible intervals) in~\emph{C.
limbatus~}(Figure 2, Table 1).~

The~\emph{stairwayplot}~showed a nearly similar dynamic for~\emph{C.
amblyrhynchos}~and~\emph{C. limbatus}, characterized by a strong
ancestral expansion (Figure 3). When approaching T=0, both species
underwent a bottleneck but of distinct strength. This is consistent with
the shape of the normalized SFS, which clearly shows a stronger deficit
in low frequency variants for~\emph{C. limbatus}~compared to~\emph{C.
amblyrhynchos}~(Figure 3).~Similarly to~\emph{C. limbatus},~\emph{C.
melanopterus}~experienced a recent 10-fold population collapse around
20,000 years B.P. starting from a long term constant~\emph{Ne}.
However,~\emph{C. melanopterus}showed no signature of ancestral
expansion,~consistent to the results obtained by~Maisano Delser et al.
(2019)~using~\emph{abc-skyline}~method. Finally,~\emph{G.
cuvier}~displayed an ancestral expansion around 100,000 years B.P.
with~\emph{Ne}~reaching \textasciitilde{}12,000 before dropping to
\textasciitilde{}3000 at T \textasciitilde{}1,600 years B.P. Remarkably,
the ancestral expansion retrieved by the~\emph{stairwayplot}(Figure 3)
for both~\emph{C. amblyrhynchos}~and~\emph{C. limbatus~}overlap with the
posterior distribution of~\emph{T}\textsubscript{\emph{COL}}~estimated
by the SST model (Table 1). This analogy holds too for~\emph{C.
melanopterus,~}where\emph{~T}\textsubscript{\emph{COL}}~could not be
properly estimated under the structured model (we obtained a flat
posterior distribution, Figure 2) and there was no signature of
ancestral expansion in the~\emph{stairwayplot}~(Figure 3).~~\emph{~}

The first set of coalescent simulations was run under FIM and SST only
(Table 2 and S2 to check if simulated data could reproduce the pattern
of genetic variability (both~q~estimators and~\emph{TD}) observed
for~\emph{C. melanopterus},~\emph{C.limbatus}, and~\emph{C.
amblyrynchos}. The simulated~q~values (excluding the equilibrium model)
ranged between 0.001 and 0.003 per site, in line with the observed
values (Table 1 and 2).~\emph{TD}~follows a U-shaped distribution for
each~\emph{Nm}~value as a function
of~\emph{T}\textsubscript{\emph{COL}}, being more positive at
recent~\emph{T}\textsubscript{\emph{COL}}~and at equilibrium and less
positive (or negative for higher~\emph{Nm}) at intermediate values.
Therefore, species demography with~\emph{Nm}~\textasciitilde{}10 (and
higher) and~\emph{T}\textsubscript{\emph{COL}}~within 15k and 50k
generations B.P. will have negative~\emph{TD}~values. In contrast,
species with lower~\emph{Nm}~and very recent or very
ancient~\emph{T}\textsubscript{\emph{COL}}~will have positive~\emph{TD}.
This matches strikingly the~\emph{TD}~observed for the three shark
species and their estimated demographic parameters under SST (Table 1).
We plot the normalized SFS and the~\emph{stairwayplot}~for all scenarios
presented in Table 2 (Figures 4, 5, S1, S2 and S3). First, we note that
none of our scenarios, even those at equilibrium and with no variation
in~\emph{Nm~}through time, showed a normalized SFS compatible with a
constant size population (Figures 4, 5, S1, S2 and S3). The normalized
SFS and the reconstructed~\emph{stairwayplot}~depend generally on the
interaction
between~\emph{Nm}~and~\emph{T}\textsubscript{\emph{COL}}~with a dynamic
strikingly similar to~\emph{TD}~(which is indeed a summary of the SFS).
For~\emph{Nm=1}~we observed the signature of a recent decrease
in~\emph{Ne~}for all scenarios and independently
of~\emph{T}\textsubscript{\emph{COL}}~(Figure 4). The normalized SFS
showed consistently a strong deficit of low frequency variants, typical
of a demographic bottleneck and in agreement with the
positive~\emph{TD}~(Figure 4 and Table 1). Furthermore,
the~\emph{stairwayplot}~could never detect the ancestral expansion for
any~\emph{T}\textsubscript{\emph{COL}}. For growing~\emph{Nm,}~the
interplay with~\emph{T}\textsubscript{\emph{COL}}~becomes more complex.
A general result is that, once again, all scenarios were characterized
by a recent decrease of~\emph{Ne~}when looking at
the~\emph{stairwayplot}~and a deficit of singletons compared to the
other low frequency classes when looking at the normalized SFS (Figures
5, S1 and S2). However, a strong signature of ancestral expansion
appeared for~\emph{Nm}~\textgreater{}10
and~\emph{T}\textsubscript{\emph{COL}}~between 15k and 50k generations
B.P., mirroring the results of~\emph{TD}~for which most of these
scenarios displayed a negative value. Remarkably,
the~\emph{stairwayplot}~retrieved the ancestral expansion only slightly
overestimating the simulated~\emph{T}\textsubscript{\emph{COL}}~(Figures
5, S1 and S2).~~Similar results were obtained for FIM (Figures S4 and
S5).~

We compared SST vs SST-CH model~(Figure 1)~by means of the same ABC-RF
model selection framework previously adopted.~The two models cannot be
clearly distinguished in any of the three structured species
since:~\emph{i}) they showed similar posterior probability
(\textasciitilde{}0.50);~\emph{ii}) the prior error rates are large
\textasciitilde{}0.40 (Table S9);~\emph{iii}) posterior distributions
of~\emph{Nm}before and after~\emph{T}\textsubscript{\emph{CH}}~are wide
and largely overlapping (Table S10);~\emph{iv}) the normalized SFS
closest to the observed data retrieved under the two models are very
similar (Figure S6). We ran a second set of coalescent simulations
focusing on the consequences of a recent change in connectivity on the
observed SFS and the reconstructed~\emph{stairwayplot~}(Table S3 and
S4). The decrease in connectivity was simulated by reducing
either~\emph{m}~(the migration rate per generation)
or~\emph{N}\textsubscript{\emph{DEME}}~(the effective population size of
each deme). As expected, we found a signature of recent population
decline in all simulated scenarios, with its intensity only slightly
affected by the change in~\emph{Nm}~(Figures 6, S7, S8 and S9). However,
the drop in~\emph{N}\textsubscript{\emph{DEME}}~(Figures 6 and S9) had
larger effect compared to the drop in~\emph{m}~(Figures S7 and S8) on
both the normalized SFS and the expansion time estimated by
the~\emph{stairwayplot}. In scenarios with 100x reduction
in~\emph{N}\textsubscript{\emph{DEME}}, the~\emph{stairwayplot}~could
not retrieve the ancestral expansion even for large
long-term~\emph{Nm}~(Figure 6). FIM-CH models displayed a behaviour
similar to SST-CH models but more pronounced (Figures S10, S11, S12 and
S13, Table S5 and S6). While at~\emph{T}\textsubscript{\emph{CH}}~= 10 a
decrease in~\emph{Nm}~slightly affected the SFS and the
reconstructed~\emph{stairwayplot}, the consequence of the change in
connectivity are more substantial at~\emph{T}\textsubscript{\emph{CH}}~=
50, with a stronger deficit in singletons and a more pronounced recent
decline in~\emph{Ne}~particularly in scenarios with a 100-fold reduction
of~\emph{N}\textsubscript{\emph{DEME}}~(Figures S11 and S13).

\textbf{Discussion}

\textbf{\emph{Life history traits and demographic history of the four
shark species}}

Discriminating whether the most appropriate model to reconstruct the
demographic history of a species
is~\emph{structured}~or~\emph{unstructured}~should be the first step in
empirical population genetics investigations, particularly when
targeting species of conservation concerns. Even when an extensive
spatial sampling is lacking,~an ABC model selection approach can
actually distinguish whether the sampled deme belongs or not to a
meta-population (similarly to previous studies(Maisano Delser et al.,
2019; Peter, Wegmann, \& Excoffier, 2010)).~Among the four species
considered here, the tiger shark is the only panmictic. The three other
species conversely are best described by the SST model, i.e., the
sampled populations belong to a meta-population exchanging migrants
following a stepping stone matrix. Our results reflect the tight link
between the level of meta-population structure (or its absence) and life
history traits. The panmictic~\emph{G. cuvier}unsurprisingly can
accomplish transoceanic movements and has the largest body size among
the sharks here considered (Table S1). In the three other species, the
estimated number of migrants (\emph{Nm}) remarkably follows the increase
of movement range (Table 1 and S1) and it is consistent with their
behaviour and habitat use. Indeed,~\emph{C. melanopterus}, a strongly
lagoon dependent species, displays the lowest level of connectivity
among the studied species (Tables 1 and S1). These results bring
meaningful hints about the influence of life history traits on
population structure in sharks, but more studies addressing this topic
will be needed to accurately detect which traits best predict its
extent.~

\textbf{\emph{Gene genealogies in the four shark species and simulated
scenarios}}

While it may seem counterintuitive to apply~\emph{unstructured~}models
to demes belonging to a meta-population, we further investigated the
demographic history of the four species by means of
the~\emph{stairwayplot}. When enough data is available,
non-parametric~\emph{unstructured}~models (such as the PSMC~(Li \&
Durbin, 2011), the extended Bayesian skyline plot~(Heled \& Drummond,
2008)~and the~\emph{stairwayplot}~among others) provide a careful
description of the distribution of coalescence times of the gene
genealogy, which ultimately depends from the ``true'' demographic
history (whether it is known or not) of the sampled lineages. If
panmixia is the most likely scenario, the distribution of coalescence
times is directly related to the variation of~\emph{Ne}~through time and
can therefore have a direct biological interpretation. This is the case
for~\emph{G. cuvier~}(Table 1),
whose~reconstructed~\emph{stairwayplot}~suggests that this species
experienced a mild ancestral expansion and a recent
\textasciitilde{}4-fold bottleneck around 2,000 years B.P. (consistent
with the results of~Pirog et al. (2019), Figure 3). Conversely, signals
detected by~the~\emph{stairwayplot}~in the remaining three species,
better described by the SST model (Table 1), cannot be directly
interpreted as changes in~\emph{Ne~}over time.~In this light,~we
ran~coalescence simulations to provide helpful and general insights into
the understanding of the relation between the inferences performed
under~\emph{unstructured}~and~\emph{structured}~models.~

We first focus on scenarios simulated under the SST, with parameters
close to those estimated in real data. The first and most striking
result is that we systematically observed a recent bottleneck under all
simulated scenarios (Table 2, Figures 4, 5, S1, S2 and S3). This result
could seem at a first glance surprising and due to an artefact. However,
this is not the case, as: i) the signal does not depend on the
inferential algorithm chosen to analyse the data (i.e.,
the~\emph{stairwayplot}), since the normalized spectra showed a deficit
in singletons compared to the other low frequency classes (Figures 4, 5,
S1 and S2), which is typical of a recent population decline; ii) it is
consistent~with the distribution of the Inverse Instantaneous
Coalescence~Rate (IICR) computed in one diploid individual, which shows
a signature of decline under similar meta-population models~(Chikhi et
al., 2018; Mazet, Rodríguez, Grusea, Boitard, \& Chikhi, 2016; Rodríguez
et al., 2018).~The results of our simulations are consistent with the
recent bottleneck observed in the three shark species (Figure 3), with
its intensity inversely correlated to the
estimated~\emph{Nm~}(i.e.,\emph{~}stronger for~\emph{C.
melanopterus}~and~\emph{C. limbatus}~than for~\emph{C. amblyrhynchos}).
In our SST model there is an instantaneous colonization of the array of
demes at~\emph{T}\textsubscript{\emph{COL}}, which corresponds also to a
demographic expansion (i.e., the total number of individuals in the
array of deme is larger than those in the ancestral deme). However, this
signature is detected only for~\emph{Nm}\selectlanguage{english}[?]5
when~\emph{T}\textsubscript{\emph{COL}}~is neither too recent nor too
old (at equilibrium) (Figures 5, S1, S2 and S3). In these scenarios, the
beginning of the expansion retrieved by the~\emph{stairwayplot}broadly
corresponds to the simulated~\emph{T}\textsubscript{\emph{COL}}. This
again corroborates the results obtained for the three shark species,
since the two species with higher~\emph{Nm}~displayed indeed an
ancestral expansion in the~\emph{stairwayplot}~with a timing consistent
with the estimated~\emph{T}\textsubscript{\emph{COL}}~(Table 1, Figures
2 and 3). Similarly, it explains why we could not retrieve the ancestral
expansion for~\emph{C. melanopterus}~nor
estimate~\emph{T}\textsubscript{\emph{COL}}~under the SST model: this
appears to be a property of the coalescence pattern and it is not
related to the amount of data available (see below).~

\textbf{\emph{Coalescence phases in structured models}}

It is now straightforward to frame all these findings under the
coalescence~perspective. The coalescence~history of the lineages sampled
from a single deme in an SST (or FIM) model can be separated for
simplicity into three phases: the~\emph{scattering},
the~\emph{collecting}~and the~\emph{ancestral}~phase~(Figure 7). Going
backward in time, lineages will coalesce in the sampled deme with a rate
according to
both~\emph{Nm}~and~\emph{N}\textsubscript{\emph{DEME}}~until all
lineages either have coalesced or migrated to another deme. This is
the~\emph{scattering~}phase described in the seminal works
of\emph{~~}(Wakeley, 1998, 1999). The~\emph{scattering~}phase was
considered instantaneous for mathematical tractability, with its outcome
dependent on~\emph{Nm}~only, but later works could disentangle the
effect of~\emph{N}\textsubscript{\emph{DEME}}~and~\emph{m}~on the shape
of the gene genealogy~(Mona, 2017). The~\emph{collecting~}phase starts
when the lineages which did not coalesce have migrated to other demes of
the array: they will then coalesce according to a Kingman process with a
rate scaled by~\emph{Nm}~and the number of demes~\emph{d}~of the
array~(Wakeley, 1999)~(Figure 7).~Finally, all surviving lineages (in
non-equilibrium model) will reach the ancestral deme
at~\emph{T}\textsubscript{\emph{COL}}, where they will coalesce at a
rate depending only on the~\emph{N}\textsubscript{\emph{ANC}}~parameter
(Figure 7).~The interplay between the demographic parameters
(\emph{N}\textsubscript{\emph{DEME,~}}\emph{Nm},~\emph{N}\textsubscript{\emph{ANC}},\emph{~d})
and the historical events
(\emph{T}\textsubscript{\emph{COL}}~and\emph{~T}\textsubscript{\emph{CH}})
determines the length of each coalescence phase and the resulting shape
of the gene genealogy of the sampled lineages (Figure 7).{}

In species with low~\emph{Nm}, the rate of coalescence~during
the~\emph{scattering}~phase is very fast since lineages have low
probability of emigrating from the sampled deme and high probability of
coalescence due to the small~\emph{N}. Once all the lineages are
dispersed in the array of demes, there will be two possible outcomes: i)
in equilibrium model, we shift to the~\emph{collecting~}phase, where the
rate of coalescence~drops since lineages will hardly fall in the same
deme again; ii) in non-equilibrium model, with the parameters we have
simulated here, there will be very few (if any) coalescence~events
during the collecting phase and the transition will be directly from
the~\emph{scattering~}to the~\emph{ancestral}~phase. Both
the~\emph{collecting~}and the~\emph{ancestral~}phases have a rate of
coalescence~lower than the~\emph{scattering~}phase, which determines the
observed recent drop in~\emph{Ne~}for all simulated scenarios.
Remarkably, the decline in~\emph{Ne}~is much stronger in equilibrium
model, since the rate of coalescence~is much lower in
the~\emph{collecting}~than in the~\emph{ancestral~}phase (Figures 4, 5,
S1, S2 and S3). Low~\emph{Nm}~species will therefore have only two
coalescence phases, the~\emph{scattering}~and either
the~\emph{collecting~}(in equilibrium model) or the~\emph{ancestral~}(in
non-equilibrium model) which is why the signature of the ancestral
expansion is lost.~

For growing~\emph{Nm}, in equilibrium model there will be again only two
coalescence~phases, namely the~\emph{scattering}~and~\emph{collecting},
with the latter having a lower rate of coalescence~than the former
independently of the simulated parameters. This is why we observed
always a strong bottleneck consistent with the distribution of the IICR
statistics in any equilibrium model~(Chikhi et al., 2018; Mazet et al.,
2015; Rodriguez et al., 2018).~In non-equilibrium model, there will be
two different situations: a)~\emph{T}\textsubscript{\emph{COL}}~(in
generations) is of the same order of the deme
size~\emph{N}\textsubscript{\emph{DEME}}. In this setting, going
backward in time few lineages would have escaped the sampled demes
before~\emph{T}\textsubscript{\emph{COL}}. This corresponds to a shift
in the coalescence~rate directly from the~\emph{scattering~}to
the~\emph{ancestral}~phase, resulting in a bottleneck of lower intensity
compared to an equilibrium model (Figures 4, 5, S1 and S2), for the same
reasons as above; b)\emph{~T}\textsubscript{\emph{COL}}~(in generations)
is larger than~\emph{N}\textsubscript{\emph{DEME}}. In this setting,
some coalescence events may occur during the~\emph{collecting~}phase, at
a rate much slower than the two other phases. This determines the hump
observed in the~\emph{stairwayplot}~(Figures 4, 5, S1 and S2) and
explains why in this window of parameters it is also possible to
correctly estimate~\emph{T}\textsubscript{\emph{COL}}~using our ABC
framework.~Further simulations under the FIM model confirmed those
patterns even though the ancestral expansion could be detected for lower
long-term~\emph{Nm~}than the corresponding SST scenario\emph{~}(Figure
S4). This is probably due to a higher apparent connectivity underlined
the by FIM, where lineages can move more freely during the collecting
phase in comparison to SST where migrants only come from the closest
neighbours. If many coalescence~events occur during
the~\emph{collecting~}phase, the change in coalescence~rate will affect
the resulting gene genealogy and it will be detected by
the~\emph{stairwayplot~}(or any other~\emph{unstructured}~method based
on coalescent~theory).{}

\textbf{\emph{Changes in connectivity~~}\emph{{}}}

Using coalescence arguments, we clarified why simple meta-population
models with constant connectivity generate a gene genealogy harbouring a
signature of a recent decline for any parameters' combination. The
signature of bottleneck detected by the~\emph{stairwayplot}~in the three
shark species best described by SST can be therefore interpreted as a
consequence of the underlying structure. However, connectivity likely
changes through time. For instance, human activities have likely
impacted the evolutionary history of a large number of species either by
decreasing their effective population size and/or by fragmenting their
habitat (i.e., reducing migration rates between demes). This intuitively
should exacerbate the signature of population decline in the resulting
gene genealogy. However, it remains to be shown whether this signature
is qualitatively and quantitatively distinguishable from models with
constant connectivity. This is a question of fundamental importance to
understand whether it is possible to detect recent bottleneck in
structured populations. To this end, we further investigated by
coalescent simulations the expected gene genealogy in SST-CH (and
FIM-CH) models with a change in connectivity 10 or 50 generations B.P.,
which matches the beginning of extensive anthropogenic influence on
biodiversity considering our species' generation time~(Ceballos et al.,
2015). The resulting gene genealogies were poorly affected by the recent
drop in connectivity, with both the normalized SFS and the
inferred~\emph{Ne}~dynamic following the same trajectory of the
corresponding scenario with the same
long-term~\emph{Nm}~and~\emph{T}\textsubscript{\emph{COL}}~(Figures 6,
S7, S8, S9, S10, S11, S12 and S13). We noticed the drop
in~\emph{N}\textsubscript{\emph{DEME}}~(Figures~6, S9, S12 and S13) had
stronger influence than the drop in~\emph{m~}(Figures~S7, S8, S10 and
S11), consistent with previous finding showing that the distribution of
coalescence~events depends not only by the~\emph{Nm}~compound parameter
but also by their individuals values~(Mona, 2017).~This can be explained
once again in the light of the length of the coalescence phases (Figure
7). Reducing~\emph{N}\textsubscript{\emph{DEME}}~will increase
exponentially the number of coalescence events, drastically shortening
the~\emph{scattering}~phase and the number of surviving lineages.
Reducing~\emph{m}~will only linearly reduce the probability of
migrations outside the deme, marginally affecting the length of
the~\emph{scattering~}phase and the number of surviving lineages
compared to constant~\emph{Nm}~scenarios. This is why a 100-fold
reduction in~\emph{N}\textsubscript{\emph{DEME}}~significantly reduces
the number of lineages entering in the~\emph{collecting~}phase, almost
hiding the ancestral expansion in high long-term~\emph{Nm}~scenarios
(Figures 6, S9, S12 and S13), while a 100-fold reduction in~\emph{m}~is
barely detectable (Figures S7, S8, S10 and S11). Similarly, the recent
reduction in
either~\emph{N}\textsubscript{\emph{DEME}}\emph{~}or\emph{~m}~cannot be
detected for lower long-term~\emph{Nm~}scenarios, where
the~\emph{collecting}~phase is already missing.~This explains why the
general pattern is strikingly similar between SST-CH and SST
simulations, which implies~that the simulated change in connectivity is
too recent to significantly alter the pattern of coalescence~events and
that a recent drop~can be hardly detected on the basis of the SFS only.
Our empirical data are consistent with these findings: when we
compared~SST vs. SST-CH models in the three shark species using the ABC
framework, we failed to clearly distinguish the two models (Tables S9
and S10, Figure S6).~This seems to be a paradox: we observed a recent
bottleneck in species of conservation concern
using~\emph{unstructured~}model, but we cannot exclude that this is just
the consequence of population structure.~{}

\textbf{\emph{Practical recommendations and conservation concerns}}{}

This study highlight once more the importance~to explicitly test for
meta-population structure before interpreting~the demographic signals
detected by~\emph{unstructured}~models, similarly to what advocated
previously by~(Maisano Delser et al., 2019; Rodriguez et al., 2018). If
the meta-population structure hypothesis is rejected, the variation
of~\emph{Ne~}through time can be directly interpreted as the demographic
history of the population under investigation, such as the case of tiger
shark. Otherwise, this variation is still related to demographic events,
but it has to be explained in the light of population structure and its
consequence on the rate of coalescence~events. We showed by coalescent
simulations how to interpret such variation: the recent bottleneck
detected by the~\emph{stairwayplot~}in demes belonging to a
meta-population is a consequence of the coalescence process. In other
words,~any inferential method implementing an\emph{~unstructured}~model
will detect such decline (if enough data is available) since it is a
property of the gene genealogy. Importantly, the gene genealogy is only
slightly affected by recent changes in connectivity if the time of this
change in generations is of the same order of the size of the deme.~

Our study underscore a key issue in conservation genetics as a~recent
decline inferred by an~\emph{unstructured}~model can be mis-interpreted
as a consequence of recent anthropic pressures~(Ceballos et al.,
2015)~when it actually results from meta-population structure. This is
all the more alarming~since the majority of species is likely organised
in meta-populations across their range rather than panmictic at a large
scale. We therefore stress the necessity for an educated choice of tools
to correctly uncover the recent trend of a species and design proper
conservation programs. For instance, detecting a recent bottleneck in
meta-populations will require~summary statistics measuring the linkage
disequilibrium~(Boitard, Rodriguez, Jay, Mona, \& Austerlitz, 2016;
Kerdoncuff, Lambert, \& Achaz, 2020)~and/or the inferential framework
based on the IICR~(Chikhi et al., 2018; Rodriguez et al., 2018)~coupled
with whole genome data.~On a positive note, we showed that the
colonization time of the array of demes can be estimated to some extent
(and under some combinations of parameters)
by~\emph{unstructured~}models. We believe that this is particularly
important because it has been shown that the simple instantaneous
colonization process we used here behaves similarly to a spatial
explicit range expansion~(Hamilton, Stoneking, \& Excoffier, 2005; Mona,
2017), which is certainly a more realistic model but more difficult to
investigate. We are aware that the meta-population models here tested
are simple and the parameters chosen are specific of the three shark
species we focused on. Nevertheless, the time-scale separation of the
coalescence~process is general, and it allows explaining intuitively any
structured models. The four shark species here used as an example has
the merit to cover a large spectrum of LHT and consequently a large
spectrum of demographic scenarios, going from a highly structured to a
panmictic population: this has strong implications on the distribution
of coalescence~times and therefore on the interpretation of the observed
data.

\textbf{\emph{Conclusion}}

In this study we found that population structure, independently from the
degree of connectivity between demes and the migration matrix relating
them, intrinsically determines a variation in the rate of
coalescence~events through time. We showed that the intensity and the
direction(s) of such variation related to the demographic parameters of
the meta-population in a predictable way. Our results highlight the
importance of detecting population structure (which depends on LHT among
other factors) before performing any demographic inferences but, at the
same time, they reveal the utility of~\emph{unstructured}~models to
describe the shape of the gene genealogy, which is the final product of
the evolutionary history of a species. A combination of structured
and~\emph{unstructured~}models (better if non-parametric) is therefore
the key to best characterize the evolutionary history of a species.~We
call for a change in perspective when investigating the demographic
history of a species: the focus should be put in the reconstruction of
the variation of both~\emph{N}~and~\emph{m}~through time, which requires
certainly new methodological development and probably more data.~

\textbf{Acknowledgement}

We are grateful to the Genotoul bioinformatics platform Toulouse
Midi-Pyrenees (Bioinfo Genotoul; http://bioinfo.genotoul.fr/) for
providing computing resources. We are indebted to Oscar Lao for fruitful
discussions and careful reading of the manuscript. This work was
supported by two ATM grants (2016 and 2017) from the Museum National
d'Histoire Naturelle to S.M.

\textbf{References}

Arenas, M., Ray, N., Currat, M., \& Excoffier, L. (2012). Consequences
of range contractions and range shifts on molecular
diversity.~\emph{Molecular Biology and Evolution},~\emph{29}(1),
207--218. https://doi.org/10.1093/molbev/msr187

Bertorelle, G., Benazzo, A., \& Mona, S. (2010). ABC as a flexible
framework to estimate demography over space and time: Some cons, many
pros.~\emph{Molecular Ecology},~\emph{19}(13), 2609--2625.
https://doi.org/10.1111/j.1365-294X.2010.04690.x

Boitard, S., Rodriguez, W., Jay, F., Mona, S., \& Austerlitz, F. (2016).
Inferring Population Size History from Large Samples of Genome-Wide
Molecular Data - An Approximate Bayesian Computation
Approach.~\emph{PLoS Genetics},~\emph{12}(3), e1005877.
https://doi.org/10.1371/journal.pgen.1005877

Ceballos, G., Ehrlich, P. R., Barnosky, A. D., Garcia, A., Pringle, R.
M., \& Palmer, T. M. (2015). Accelerated modern human-induced species
losses: Entering the sixth mass extinction.~\emph{Science
Advances},~\emph{1}(5), 9--13. https://doi.org/10.1126/sciadv.1400253

Chikhi, L., Rodriguez, W., Grusea, S., Santos, P., Boitard, S., \&
Mazet, O. (2018).~The IICR (inverse instantaneous coalescence rate) as a
summary of genomic diversity: Insights into demographic inference and
model choice.~\emph{Heredity},~\emph{120}(1), 13--24.
https://doi.org/10.1038/s41437-017-0005-6

Chikhi, L., Sousa, V. C., Luisi, P., Goossens, B., \& Beaumont, M. A.
(2010).~The Confounding Effects of Population Structure, Genetic
Diversity and the Sampling Scheme on the Detection and Quantification of
Population Size Changes.~\emph{Genetics},~\emph{186}(3), 983--995.
https://doi.org/10.1534/genetics.110.118661

Corrigan, S., Lowther, A. D., Beheregaray, L. B., Bruce, B. D., Cliff,
G., Duffy, C. A., \ldots{} Rogers, P. J. (2018). Population Connectivity
of the Highly Migratory Shortfin Mako (\emph{Isurus
oxyrinchus~}Rafinesque 1810) and Implications for Management in the
Southern Hemisphere.~\emph{Frontiers in Ecology and
Evolution},~\emph{6}(NOV), 1--15.
https://doi.org/10.3389/fevo.2018.00187

Cortes, E. (2002). Incorporating uncertainty into demographic modeling:
Application to shark populations and their
conservation.~\emph{Conservation Biology},~\emph{16}(4), 1048--1062.
https://doi.org/10.1046/j.1523-1739.2002.00423.x

Excoffier, L., Dupanloup, I., Huerta-Sanchez, E., Sousa, V. C., \& Foll,
M. (2013). Robust Demographic Inference from Genomic and SNP
Data.~\emph{PLoS Genetics},~\emph{9}(10), e1003905.
https://doi.org/10.1371/journal.pgen.1003905

Excoffier, L., Foll, M., \& Petit, R. J. (2009). Genetic Consequences of
Range Expansions.~\emph{Annual Review of Ecology, Evolution, and
Systematics},~\emph{40}(1), 481--501.
https://doi.org/10.1146/annurev.ecolsys.39.110707.173414

Hamilton, G., Stoneking, M., \& Excoffier, L. (2005). Molecular analysis
reveals tighter social regulation of immigration in patrilocal
populations than in matrilocal populations.~\emph{Proceedings of the
National Academy of Sciences of the United States of
America},~\emph{102}(21), 7476--7480.
https://doi.org/10.1073/pnas.0409253102

Heled, J., \& Drummond, A. J. (2008). Bayesian inference of population
size history from multiple loci.~\emph{BMC Evolutionary
Biology},~\emph{8}(1), 289. https://doi.org/10.1186/1471-2148-8-289

Heller, R., Nursyifa, C., Garcia\selectlanguage{english}-Erill, G., Salmona, J., Chikhi, L.,
Meisner, J., \ldots{} Albrechtsen, A. (2021). A reference\selectlanguage{english}-free approach
to analyse RADseq data using standard next generation sequencing
toolkits.~\emph{Molecular Ecology Resources},~\emph{21}(4), 1085--1097.
https://doi.org/10.1111/1755-0998.13324

Karl, S. A., Motta, P. J., Stewart, B. S., Wilson, S. G., Bowen, B. W.,
Castro, A. L. F., \ldots{} Hueter, R. E. (2010). Population genetic
structure of Earth's largest fish, the whale shark (\emph{Rhincodon
typus}).~\emph{Molecular Ecology},~\emph{16}(24), 5183--5192.
https://doi.org/10.1111/j.1365-294x.2007.03597.x

Kerdoncuff, E., Lambert, A., \& Achaz, G. (2020). Testing for population
decline using maximal linkage disequilibrium blocks.~\emph{Theoretical
Population Biology},~\emph{134}, 171--181.
https://doi.org/10.1016/j.tpb.2020.03.004

Khimoun, A., Doums, C., Molet, M., Kaufmann, B., Peronnet, R., Eyer, P.
A., \& Mona, S. (2020). Urbanization without isolation: The absence of
genetic structure among cities and forests in the tiny acorn ant
Temnothorax nylanderi.~\emph{Biology Letters},~\emph{16}(1).
https://doi.org/10.1098/rsbl.2019.0741

Korneliussen, T. S., Albrechtsen, A., \& Nielsen, R. (2014). ANGSD:
Analysis of Next Generation Sequencing Data.~\emph{BMC
Bioinformatics},~\emph{15}(1), 1--13.
https://doi.org/10.1186/s12859-014-0356-4

Lapierre, M., Lambert, A., \& Achaz, G. (2017). Accuracy of Demographic
Inferences from the Site Frequency Spectrum: The Case of the Yoruba
Population.~\emph{Genetics},~\emph{206}(1), 439--449.
https://doi.org/10.1534/genetics.116.192708

Li, H., \& Durbin, R. (2009). Fast and accurate short read alignment
with Burrows-Wheeler transform.~\emph{Bioinformatics},\emph{25}(14),
1754--1760. https://doi.org/10.1093/bioinformatics/btp324

Li, H., \& Durbin, R. (2011). Inference of human population history from
individual whole-genome sequences.~\emph{Nature},~\emph{475}(7357),
493--496. https://doi.org/10.1038/nature10231

Liu, X., \& Fu, Y.-X. (2015).~Exploring population size changes using
SNP frequency spectra.~\emph{Nature Genetics},~\emph{47}(5), 555--559.
https://doi.org/10.1038/ng.3254

Maisano Delser, P., Corrigan, S., Duckett, D., Suwalski, A., Veuille,
M., Planes, S., \ldots{} Mona, S. (2019). Demographic inferences after a
range expansion can be biased: the test case of the blacktip reef shark
(\emph{Carcharhinus melanopterus}).~\emph{Heredity},~\emph{122}(6),
759--769. https://doi.org/10.1038/s41437-018-0164-0

Maisano Delser, P., Corrigan, S., Hale, M., Li, C., Veuille, M., Planes,
S., \ldots{} Mona, S. (2016). Population genomics of~\emph{C.
melanopterus}~using target gene capture data: Demographic inferences and
conservation perspectives.~\emph{Scientific Reports},~\emph{6}(April),
1--12. https://doi.org/10.1038/srep33753

Mazet, O., Rodriguez, W., Grusea, S., Boitard, S., \& Chikhi, L. (2016).
On the importance of being structured: instantaneous coalescence rates
and human evolution---lessons for ancestral population size
inference?~\emph{Heredity},~\emph{116}(4), 362--371.
https://doi.org/10.1038/hdy.2015.104

Mazet, Olivier, Rodriguez, W., \& Chikhi, L. (2015). Demographic
inference using genetic data from a single individual: Separating
population size variation from population structure.~\emph{Theoretical
Population Biology},~\emph{104}, 46--58.
https://doi.org/10.1016/j.tpb.2015.06.003

Mona, S. (2017). On the role played by the carrying capacity and the
ancestral population size during a range
expansion.~\emph{Heredity},~\emph{118}(2), 143--153.
https://doi.org/10.1038/hdy.2016.73

Mona, S., Ray, N., Arenas, M., \& Excoffier, L. (2014). Genetic
consequences of habitat fragmentation during a range
expansion.~\emph{Heredity},~\emph{112}(3), 291--299.
https://doi.org/10.1038/hdy.2013.105

Paris, J. R., Stevens, J. R., \& Catchen, J. M. (2017). Lost in
parameter space: a road map for stacks.~\emph{Methods in Ecology and
Evolution},~\emph{8}(10), 1360--1373.
https://doi.org/10.1111/2041-210X.12775

Parsons, G. R. (1990). Metabolism and swimming efficiency of the
bonnethead shark~\emph{Sphyrna tiburo}.~\emph{Marine
Biology},~\emph{104}(3), 363--367. https://doi.org/10.1007/BF01314338

Peter, B. M., Wegmann, D., \& Excoffier, L. (2010). Distinguishing
between population bottleneck and population subdivision by a Bayesian
model choice procedure.~\emph{Molecular Ecology},~\emph{19}(21),
4648--4660. https://doi.org/10.1111/j.1365-294X.2010.04783.x

Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., \& Hoekstra,
H. E. (2012). Double Digest RADseq: An Inexpensive Method for De Novo
SNP Discovery and Genotyping in Model and Non-Model Species.~\emph{PLoS
ONE},~\emph{7}(5), e37135. https://doi.org/10.1371/journal.pone.0037135

Pirog, A., Jaquemet, S., Ravigne, V., Cliff, G., Clua, E., Holmes, B.
J., \ldots{}~Magalon, H. (2019). Genetic population structure and
demography of an apex predator, the tiger shark~\emph{Galeocerdo
cuvier}.~\emph{Ecology and Evolution},~\emph{9}(10), 5551--5571.
https://doi.org/10.1002/ece3.5111

Pudlo, P., Marin, J.-M. M., Estoup, A., Cornuet, J.-M. M., Gautier, M.,
\& Robert, C. P. (2016). Reliable ABC model choice via random
forests.~\emph{Bioinformatics},~\emph{32}(6), 859--866.
https://doi.org/10.1093/bioinformatics/btv684

Ray, N., Currat, M., \& Excoffier, L. (2003). Intra-deme molecular
diversity in spatially expanding populations.~\emph{Molecular Biology
and Evolution},~\emph{20}(1), 76--86.
https://doi.org/10.1093/molbev/msg009

Raynal, L., Marin, J. M., Pudlo, P., Ribatet, M., Robert, C. P., \&
Estoup, A. (2019). ABC random forests for Bayesian parameter
inference.~\emph{Bioinformatics},~\emph{35}(10), 1720--1728.
https://doi.org/10.1093/bioinformatics/bty867

Rochette, N. C., Rivera-Colon, A. G., \& Catchen, J. M. (2019).~Stacks
2: Analytical methods for paired-end sequencing improve RADseq-based
population genomics.~\emph{Molecular Ecology},~\emph{28}(21),
4737--4754. https://doi.org/10.1111/mec.15253

Rodriguez, W., Mazet, O., Grusea, S., Arredondo, A., Corujo, J. M.,
Boitard, S., \& Chikhi, L. (2018). The IICR and the non-stationary
structured coalescent: towards demographic inference with arbitrary
changes in population structure.~\emph{Heredity},~\emph{121}(6),
663--678. https://doi.org/10.1038/s41437-018-0148-0

Stadler, T., Haubold, B., Merino, C., Stephan, W., \& Pfaffelhuber, P.
(2009). The impact of sampling schemes on the site frequency spectrum in
nonequilibrium subdivided populations.~\emph{Genetics},~\emph{182}(1),
205--216. https://doi.org/10.1534/genetics.108.094904

Tajima, F. (1989). Statistical method for testing the neutral mutation
hypothesis by DNA polymorphism.~\emph{Genetics},~\emph{123}(3),
585--595. https://doi.org/10.1093/genetics/123.3.585

Wakeley, J. (1998). Segregating Sites in Wright's Island
Model.~\emph{Theoretical Population Biology},~\emph{53}, 166--174.

Wakeley, J. (1999). Nonequilibrium migration in human
history.~\emph{Genetics},~\emph{153}(4), 1863--1871.

Warmuth, V. M., \& Ellegren, H. (2019). Genotype-free estimation of
allele frequencies reduces bias and improves demographic inference from
RADSeq data.~\emph{Molecular Ecology Resources},~\emph{19}(3), 586--596.
https://doi.org/10.1111/1755-0998.12990

Watterson, G. A. A. (1975). On the number of segregating sites in
genetical models without recombination.~\emph{Theoretical Population
Biology},~\emph{7}(2), 256--276.
https://doi.org/10.1016/0040-5809(75)90020-9{}

\textbf{Data availability statement}

Fastq sequence files, SFS and scripts are available from the Dryad
Digital Repository:~\url{https://doi.org/10.5061/dryad.b8gtht7d1}.~

\textbf{Authors contribution}

S.M. and P.L. conceived the project. S.P. provided reagents and samples.
S.M. and P.L analysed the data and wrote the manuscript with input from
S.P.

\textbf{Table 1.} Summary statistics and ABC estimation. Number of loci
and SNPs after filtering, mean pairwise difference (\selectlanguage{greek}θ\selectlanguage{english}\selectlanguage{greek}\textsubscript{π}\selectlanguage{english}),
Watterson theta (\selectlanguage{greek}θ\selectlanguage{english}\textsubscript{w}), Tajima's D (TD), posterior
probability of the most supported model and its parameters (median value
and 95\% credible interval in parentheses).\selectlanguage{english}
\begin{longtable}[]{@{}llllllllll@{}}
\toprule
& N\selectlanguage{ngerman}° Loci & N\selectlanguage{ngerman}° SNP & \selectlanguage{greek}θ\selectlanguage{english}\selectlanguage{greek}\textsubscript{π} & \selectlanguage{greek}θ\selectlanguage{english}\textsubscript{w} & TD &
Model (probability)\textsuperscript{+} & \emph{Nm} & \emph{Tcol}
\textsuperscript{++} & \emph{Nanc}\tabularnewline
\midrule
\endhead
\emph{G. cuvier}\selectlanguage{ngerman}\textsuperscript{§} (N=8) & 117976 & 25785 & 0.00057 &
0.00051 & -0.03 & NS (0.84) & - & - & -\tabularnewline
\emph{C. amblyrhynchos} (N=12) & 69490 & 68355 & 0.00216 & 0.00229 &
-0.23\textsuperscript{*} & SST (0.85) & 11.5 (3.0-22.0) & 20456
(12567-75649) & 40961 (1315-49276)\tabularnewline
\emph{C. limbatus} (N=13) & 60812 & 43449 & 0.00180 & 0.00166 &
0.43\textsuperscript{*} & SST (0.55) & 6.6 (1.5-15.4) & 50198
(475-245440) & 25521 (1913-52820)\tabularnewline
\emph{C. melanopterus}\selectlanguage{ngerman}\textsuperscript{¶} (N=8) & 926 & 784 & 0.00040 &
0.00030 & 0.691\textsuperscript{*} & SST (0.89) & 1.8 (0.7-3.0) & 91719
(5000-291341) & 34607 (2760-95380)\tabularnewline
& & & & & & \emph{Priors} \textsuperscript{a} & U: 0.001 - 100 & U: 1 -
300000 & U: 100 -- 100000.\tabularnewline
\bottomrule
\end{longtable}

\textsuperscript{*} Tajima's \emph{D} values are significant (\emph{p}
\textless{}0.001).

\textsuperscript{+} Most supported model and its posterior probability.

\textsuperscript{++} \emph{Tcol} is expressed in generations.

\selectlanguage{ngerman}\textsuperscript{§} \emph{G. cuvier} is best represented by the NS
model: its demography is depicted through the
\emph{stairwayplot}algorithm (see discussion).

\textsuperscript{~¶} Data from Maisano Delser et al. (2019).

\textsuperscript{a} Uniform prior distribution. The prior distribution
of \emph{Nm} is the product of two uniforms (one for \emph{N} and one
for \emph{m} ).

\textbf{Table 2.} Coalescent simulations of 50,000 Rad-loci under SST
model, with mutation rate fixed to 1.93*10\textsuperscript{-8} per site
per generation and \emph{Nanc} fixed to 50,000. Mean pairwise difference
(\selectlanguage{greek}θ\selectlanguage{english}\selectlanguage{greek}\textsubscript{π}\selectlanguage{english}), Watterson theta (\selectlanguage{greek}θ\selectlanguage{english}\textsubscript{w}), Tajima's D
(TD), and number of segregating sites (S) are averaged over 100
replicates.\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
Nm & Tcol & \selectlanguage{greek}\textbf{θ\textsubscript{π}} \selectlanguage{english}\textsuperscript{++} &
\selectlanguage{greek}\textbf{θ\textsubscript{s}}\selectlanguage{english} \textsuperscript{++} & \textbf{TD} &
\textbf{S}\tabularnewline
\midrule
\endhead
1 & 5000 & 0.0013 & 0.0011 & 0.531 & 23599\tabularnewline
& 15000 & 0.0013 & 0.0012 & 0.405 & 24094\tabularnewline
& 50000 & 0.0017 & 0.0016 & 0.406 & 32201\tabularnewline
& [?]\textsuperscript{+} & 0.0161 & 0.0139 & 0.669 & 283564\tabularnewline
5 & 5000 & 0.0017 & 0.0016 & 0.361 & 32443\tabularnewline
& 15000 & 0.0019 & 0.0018 & 0.191 & 37712\tabularnewline
& 50000 & 0.0028 & 0.0028 & 0.035 & 56474\tabularnewline
& [?] & 0.0177 & 0.0150 & 0.749 & 306786\tabularnewline
10 & 5000 & 0.0019 & 0.0018 & 0.180 & 36561\tabularnewline
& 15000 & 0.0021 & 0.0022 & -0.087 & 44380\tabularnewline
& 50000 & 0.0031 & 0.0034 & -0.364 & 69436\tabularnewline
& [?] & 0.0180 & 0.0158 & 0.585 & 321619\tabularnewline
15 & 5000 & 0.0019 & 0.0019 & 0.048 & 38919\tabularnewline
& 15000 & 0.0022 & 0.0024 & -0.274 & 48479\tabularnewline
& 50000 & 0.0032 & 0.0038 & -0.608 & 77391\tabularnewline
& [?] & 0.0181 & 0.0163 & 0.465 & 331816\tabularnewline
\bottomrule
\end{longtable}

\textsuperscript{+} Equilibrium model obtained by simulating\emph{Tcol}
=[?].

\textsuperscript{++} Theta values are expressed per site per generation.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure1-models/Figure1-models}
\caption{{Evolutionary scenarios considered in this study (to both infer
parameters in real data under an ABC framework and to perform coalescent
simulations). SST (FIM) model is a simplified version of SST-CH (FIM-CH)
in which connectivity~\emph{Nm}~is constant
after~\emph{T}\textsubscript{\emph{COL}}. Details on each parameter are
presented in the main text.
{\label{986351}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure2-abc/Figure2-abc}
\caption{{Posterior distribution of the number of migrants per
generation~\emph{Nm~}(panel A) and of the colonisation time of the array
of deme~\emph{T\textsubscript{COL}~}(panel B) estimated under the
stepping stone model (SST) for \emph{Carcharhinus amblyrhynchos} (red),
\emph{Carcharhinus limbatus} (green) and \emph{Carcharhinus
melanopterus} (blue). ~
{\label{714459}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure3-sp-species/Figure3-sp-species}
\caption{{Panel A: variation of the effective population size (\emph{Ne}) through
time and its 75\% confidence interval estimated by
the~\emph{stairwayplot}. Panel B: normalized SFS computed as in
(Lapierre et al., 2017).~\emph{Carcharhinus amblyrhynchos} is
represented in red,~\emph{Carcharhinus limbatus} in
green,~\emph{Carcharhinus melanopterus} in blue, and~\emph{Galeocerdo
cuvier~}in purple.
{\label{261878}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure4-SST-nm1/Figure4-SST-nm1}
\caption{{\emph{stairwayplot} (maximum likelihood \emph{Ne} and 75\% confidence
interval) (panel A) and normalized SFS (panel B) computed in simulated
non-equilibrium SST scenarios with \emph{Nm}=1, averaged over 100
replicates. Colonisation time of the array of deme\emph{Tcol} occurred
5,000 (red), 15,000 (blue), and 50,000 (green) generations B.P.,
visually represented by the vertical dashed lines in panel A. The
normalized SFS expected under a constant size non-structured model (NS
constant size) is also shown (grey dashed line in panel B).
{\label{871080}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure5-SST-nm10/Figure5-SST-nm10}
\caption{{\emph{stairwayplot} (maximum likelihood~\emph{Ne} and 75\% confidence
interval) (panel A) and normalized SFS (panel B) computed in simulated
non-equilibrium SST scenarios with~\emph{Nm}=10, averaged over 100
replicates. Colonisation time of the array of deme\emph{Tcol} occurred
5,000 (red), 15,000 (blue), and 50,000 (green) generations B.P.,
visually represented by the vertical dashed lines in panel A. The
normalized SFS expected under a constant size non-structured model (NS
constant size) is also shown (grey dashed line in panel B).
{\label{106315}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure6-SST-ch-N-10/Figure6-SST-ch-N-10}
\caption{{\emph{stairwayplot} (maximum likelihood \emph{Ne}) (panel A) and
normalized SFS (panel B) computed in simulated non-equilibrium SST
scenarios with \emph{Tcol}=15,000 generations B.P. and an instantaneous
decrease of the deme size (\emph{N}) forward in time at \emph{Tch}=10
generations B.P. Colours represent the long-term connectivity values:
\emph{Nm}=1 (blue), \emph{Nm}=5 (green),\emph{Nm}=10 (red), \emph{Nm}=15
(black). Line style represents the 10-fold (small dashes) or 100-fold
(dots) reduction of \emph{N}, or constant \emph{Nm} (continuous line).
The vertical grey dashed line in panel A represents the simulated
colonisation time \emph{Tcol}.
{\label{740803}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure7-coal-rates/Figure7-coal-rates}
\caption{{Schematic diagram representing the different coalescence phases in the
history of lineages sampled from a deme belonging to a non-equilibrium
meta-population. Each phase and related parameters are represented by a
colour. Parameters influencing the coalescence rate in each phase are:
the effective size of the deme (\emph{N}\textsubscript{\emph{DEME}}) and
the migration rate (\emph{m}) for the~\emph{scattering}~phase; the
number of migrants exchanged per generation (\emph{Nm}) and the number
of demes (\emph{d}) for the~\emph{collecting}~phase; and the ancestral
effective size (\emph{N}\textsubscript{\emph{ANC}}) for
the~\emph{ancestral}~phase.~
{\label{513108}}%
}}
\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

