High genetic variation and phylogeographic relations among Palearctic fairy shrimp populations reflect persistence in multiple southern refugia during Pleistocene ice ages and postglacial colonisation

1WasserCluster Lunz, Lunz am See, Austria 2Department of Limnology and Bio‐Oceanography, University of Vienna, Vienna, Austria 3Animal Ecology, Global Change and Sustainable Development, KU Leuven, Leuven, Belgium 4Sorbonne Université, Muséum national d’Histoire naturelle, UCN, UA, CNRS, IRD, Biologie des organismes et écosystèmes aquatiques, BOREA, UMR 7208, CP26 75231, 43 rue Cuvier, 75005, Paris Cedex 05, France 5Department of Genetics and Biosystematics, Faculty of Biology, University of Gdańsk, Gdańsk, Poland 6Community Ecology Laboratory, Department of Biology, Vrije Universiteit Brussel (VUB), Brussels, Belgium 7Water Research Group, Unit for Environmental Sciences and Management, North‐West University, Potchefstroom, South Africa 8Centre for Environmental Management, University of the Free State, Bloemfontein, South Africa


| INTRODUC TI ON
Although temporary ponds are common aquatic habitats across many regions, they are increasingly threatened by human activities including urbanisation, draining, and intensification of agriculture (Silva, Phillips, Jones, Eldridge, & O'Hara, 2007;Van den Broeck, Waterkeyn, Rhazi, & Brendonck, 2015a; Van den Broeck, Waterkeyn, Rhazi, Grillas, & Brendonck, 2015b). Due to their typically small size and the fact that their filling and drying depends on rainfall and temperature, they are also particularly vulnerable to climate change (Moss, 2012;Stoks, Geerts, & Meester, 2014;Tuytens, Vanschoenwinkel, Waterkeyn, & Brendonck, 2014). However, these systems have a high ecological importance as feeding grounds for migratory birds, stepping-stones for dispersal of aquatic organisms, and habitats of a specialised aquatic fauna and flora with high degrees of endemicity (Williams, 2006).
Large branchiopod crustaceans (Crustacea, Branchiopoda; group including the fairy shrimps) are an iconic group of temporary pond inhabitants. They typically grow and mature fast as an adaptation to the short growing seasons, determined by the time-constrained wet phase of the pond. In addition, they bridge dry periods through the production of drought-resistant dormant stages (Dumont & Negrea, 2002). Dormant stages also serve as propagules for spatial dispersal via wind, flowing water, or through animal vectors (Bilton, Freeland, & Okamura, 2001;Pinceel, Brendonck, & Vanschoenwinkel, 2016).
The fairy shrimp Branchipus schaefferi Fischer 1934 occurs in temporary freshwaters across Europe, Northern Africa, and Asia (Al-Sayed & Zainal, 2005;Brtek & Thiéry, 1995). Given its distribution and phenology, B. schaefferi is considered to be a warm water species (Mura, 1999;Vanschoenwinkel, Brendonck, Pinceel, Dupriez, & Waterkeyn, 2013). In most of Europe, the species usually occurs from late spring to early autumn (Eder, Hödl, & Gottwald, 1997;Petrov & Cvetković, 1997). In warmer regions in Northern Africa, the Mediterranean, and Asia, populations have been reported throughout the year (Marrone & Mura, 2006). While the ecology, taxonomy and range of occurrence of the species have been addressed to some extent (Brtek & Thiéry, 1995;Gandolfi, Rossi, & Zarattini, 2015;Vanschoenwinkel et al., 2013), a range-wide molecular phylogeography is lacking. Given that a number of closely related species still needs to be validated (Belk & Brtek, 1995;Gandolfi et al., 2015), a full-scale genetic study would provide essential complementary information to resolve taxonomic relationships and point to meaningful taxonomic units.
Genetic data can provide highly valuable information to study the history and diversity of a species. It can, for instance, be used to detect historical gene flow among populations and to reconstruct dispersal events more precisely than with only traditional methods based on morphological features (Freeland, Kirk, & Peterson, 2012).
Conservation of the full adaptive potential of a species should have priority over simple species conservation (Moritz, 1994;Ryder, 1986;Waples, 1995). Large branchiopods are known for high levels of cryptic genetic diversity among individuals that look morphologically similar (Aguilar et al., 2017;Pinceel et al., 2013aPinceel et al., ,2013bSchwentner et al., 2013). Such individuals could differ extensively in physiology and may represent distinct evolutionary significant units (ESUs) for conservation (Pinceel et al., 2013b). Finally, phylogeographic studies may improve our understanding of the effect of past climate events, which, in turn, may serve to forecast consequences of future environmental changes (Pinceel et al., 2013a,b).
The phylogeography of a number of large branchiopod species in Europe and North Africa has been reconstructed and many studies show limited genetic divergence among populations, especially in more northern regions of mainland Europe (Kappas et al., 2017;Reniers, Vanschoenwinkel, Rabet, & Brendonck, 2013;Vanschoenwinkel et al., 2012). This has been explained as a consequence of relatively recent range expansion from a small number of refugia after the Pleistocene glacials. Two studies have been undertaken to investigate specific aspects of the phylogeny of B. schaefferi, one based on allozymes and CO1 (Gandolfi et al., 2015) and 4. Overall, the studied B. schaefferi dataset comprises high levels of genetic differentiation, without a clear morphological signal. Phylogenetic searches and pairwise genetic distances suggest that the studied lineages belong to a complex of four morphologically cryptic species. Since these four evolutionary old clades persist (±2 million years), despite overlapping geographic ranges and since they span a variety of ecological conditions, they should be considered as separate evolutionary significant units for conservation.

K E Y W O R D S
dispersal, freshwater, habitat destruction, molecular clock, temporary ponds another on 18S (Mioduchowska et al., 2018). These studies were, however, restricted to 11 populations in Italy, Spain, and Morocco (Gandolfi et al., 2015) and 11 populations in Poland, Italy, and Algeria (Mioduchowska et al., 2018).
Here, we conduct a large-scale phylogeographic study of the fairy shrimp species B. schaefferi across its range of occurrence. For this, we study the mitochondrial CO1 and nuclear ITS1 gene regions of individuals from 79 populations from wide areas in Europe and northern Africa and a single population in the Middle East.
First of all, we perform phylogenetic searches and use sequence divergence based methods to verify if molecular data support the species status of the studied specimens, which were all identified as B. schaefferi based on morphological traits. Given the extensive geographic and ecological range of occurrence of the species, we expect high levels of genetic differentiation among certain populations. Second, we use genetic divergence data among genetically distinct groups and standard molecular clocks, to assess the likelihood of different historic scenarios as explanation for the current distribution of genetic lineages. Given the fact that B. schaefferi is mostly successful under relatively high temperatures, the Pleistocene ice ages would have driven B. schaefferi to extinction in Northern regions. Therefore, we hypothesise low levels of genetic diversity in Northern regions compared to high levels of diversity around glacial refugia. Finally, based on the level of genetic differentiation between identified haplotype groups we aim to delineate the ESUs important for conservation of the adaptive potential within B. schaefferi.  Table 1) were obtained after hatching field-collected sediment with B. schaefferi egg banks in the laboratory.

| DNA extraction, polymerase chain reaction, DNA purification and sequencing
The molecular laboratory procedures to acquire the DNA sequences for the targeted genes were performed in two laboratories separately, at the Department of Genetics and Biosystematics, University of Gdansk in Poland (45 specimens from 10 Polish populations) and at the Laboratory for Animal Ecology, Global Change and Sustainable Development at KU Leuven in Belgium (all other specimens; see Supporting Information for both protocols).

| Phylogenetic and phylogeographical reconstructions
All generated B. schaefferi CO1 sequences were assembled and visually checked for quality using SeqScape v2.5. Consensus sequences were edited in BioEdit Sequence Alignment Editor (Hall, 1999). All sequences that contained insertions and/or deletions (15 in total) were removed from the CO1 alignments to avoid the risk of co-amplified nuclear mitochondrial pseudogenes interfering with the analyses (Song et al., 2008). The newly generated sequences of B. schaefferi, together with existing B. schaefferi sequences from GenBank (Gandolfi et al., 2015), one sequence of Branchipus blanchardi Daday 1908 (KP702861.1) and one outgroup taxon (CO1: Branchipodopsis drakensbergensis GU139737.1 and ITS1: Branchipodopsis wolfi MN325155), were aligned with the CLUSTALW multiple alignment tool in BioEdit. All sequences were uploaded to GenBank (for accession codes see Table 1). The most probable evolutionary model for both markers was determined in PhyML (Lefort, Longueville, & Gascuel, 2017) based on both the Bayesian information criterion and Akaike information criterion (AIC). For CO1, the AIC selected for a general time reversible model (GTR) with discrete γ model (+G; γ = 1.83) with invariable sites (I = 0.57) which was used to assemble the Bayesian inference (BI) and maximum likelihood (ML) tree. To assemble neighbour joining (NJ) trees, we used a Tamura Nei evolutionary model (TN93; Tamura & Nei, 1993) with a discrete γ distribution, which was the best scored available model for the NJ method. For ITS1, the Bayesian information criterion selected a Kimura 2-parametric model (K2P; Kimura, 1980), which was used for constructing the ML and NJ tree. The GTR with invariable sites was selected as most suitable evolutionary model by the AIC and used for assembling the BI tree since K2P models are not embedded within MrBayes. Substitution saturation was tested in DAMBE v. 7.0.28 (Xia & Kumar, 2018). The index of substitution saturation was significantly smaller than the critical index of substitution saturation, indicating little saturation (Xia & Lemey, 2009;Xia, Xie, Salemi, Chen, & Wang, 2003) for both markers. The haplotype number was determined based on calculated pairwise distances in MEGA X (Kumar, Stecher, Li, Knyaz, & Tamura, 2018).
The consensus phylogeny was constructed based on CO1 sequences by comparing phylogenetic trees obtained with four different methods of inference: NJ, ML, maximum parsimony (MP), and BI. ML analyses were performed in MEGA X and PhyML (Guindon et al., 2010) according to the GTR + G + I evolutionary model for the CO1 and K2P model for ITS1 with 1,000 bootstrap replicates.
The MP analyses for CO1 were performed in PAUP* v4.0 (Swofford, 2001)   TA B L E 1 (Continued) of 90% (<10% alignment gaps, missing data, and ambiguous bases were allowed at any position). Bayesian inference was performed in MrBayes (Huelsenbeck & Ronquist, 2001;Ronquist & Huelsenbeck, 2003;Ronquist et al., 2012). Markov Chain Monte Carlo analysis ran for 2 × 10 6 generations with a sampling frequency of 1,000 generations and a 25% burn-in. Pairwise genetic distances (based on K2P model) between all generated sequences and the mean genetic distances within and among the main groups in the phylogeny of B. schaefferi were calculated in MEGA X (Kumar, Stecher, & Tamura, 2016) with partial deletion of 90% (373 positions in the final data set).
To estimate the approximate timing of divergence among phylogenetic groups we applied various molecular clocks to the minimum and maximum of the calculated pairwise distances. For CO1, molecular clocks of 1.4 and 2.6% divergence per million years were applied as in Reniers et al. (2013). Eventually, the time range of the group split was estimated between the highest pairwise distance (D/2) divided by the lowest rate of evolution, for the most distant time scenario, and the lowest pairwise distance divided by the highest rate of evolution for the most recent likelihood of the events.
We used the automatic barcoding gap discovery method (ABGD) (Puillandre, Lambert, Brouillet, & Achaz, 2012)  can be considered as separate species when the mean distance between them is at least four times larger than the mean divergence within the lineages (Birky et al., 2005). The ABGD method separates species based on the barcode gap that occurs when the divergence between individuals of the same species is lower that the divergence between the individuals of the different species (Puillandre et al., 2012). To run the ABGD method, we used the online version (http:// wwwabi.snv.jussi eu.fr/publi c/abgd/abgdw eb.html) following the default settings.

| Genetic diversity
We identified a total of 31 unique CO1 haplotypes among the stud-

| Phylogenetic analyses based on CO1
The four different methods of phylogenetic inference (ML, MP, NJ, and BI) produced trees with a highly similar topology for the studied   (Table 2; for TN93 + G between group distances see Table S2). Mean between group distances of B. blanchardi and four B. schaefferi clades were overall higher than distances between the B. schaefferi clades (Table 2). Both ABGD (prior maximal distance p = 0.035; Table S1) and the 4x-rule suggest that the clades should be considered as different species.

| Phylogenetic analyses based on ITS1
All phylogenetic search methods divided the studied haplotypes in two groups. The first group corresponded to the clade B recognised for CO1 (three haplotypes), while the second group included all other groups (A, C, and D, overall 10 haplotypes; Figure S1). Mean F I G U R E 1 Consensus phylogenetic tree for Branchipus schaefferi, based on the mitochondrial CO1 gene fragment (maximum likelihood-ML, maximum parsimony-MP, neighbour joining-NJ and Bayesian inference-BI). The ML tree was used as a template. The supporting values of four evolution reconstruction methods are included close to the nodes (ML/MP/NJ/BI). The unsupported groupings are indicated with '-'. Codes within the first pair of brackets indicate the codes of sequenced specimens and numbers in the second pair of brackets specify the number of specimens from the same region that belong to the same haplotype. The groups (clades) identified by the phylogenetic search methods are indicated with the same colour-coding as in Figure  All phylogenetic searches grouped lineages from Tunisia together as a sub-group with two haplotypes. This is consistent with the reconstructions based on CO1, with the exception that the haplotype from Malta is not separated from the other lineages. It should be noted that, based on the ITS1 region and the partial CO1 fragment, the studied specimen from Bahrain (R1) belongs to clade D. The table contains minimum, maximum, and mean genetic distances between groups (in %) and an assessment of the timing of divergence among groups (millions of years ago, mya). The number of base substitutions per site, averaged over all sequence pairs between groups, represents the mean distances. Analyses were conducted using the Kimura 2-parameter model (Kimura, 1980). Fewer than 10% alignment gaps, missing data, and ambiguous bases were allowed at any position. Sequences were 373 base pairs in length in the final alignment. Molecular clocks ranged between 1.4% and 2.6% of substitution per my (cf. Reniers et al., 2013).

| Pleistocene glaciations shaped the evolutionary history of B. schaefferi
Temperate Europe was characterised by extreme climatic fluctuations throughout the Pleistocene (2.6 million years ago [mya]-11,000 years ago). Periods of glaciation and extending ice cover were alternated with milder interglacial periods (Paillard, 1998). The Arctic ice cover was formed 2.4 mya and until 0.9 mya the ice coverage advanced and retreated in cycles of approximately 41,000 years (Hewitt, 2000). Later on, 0.9 mya-present, the glacial/interglacial cycles became more severe and typically lasted around 100,000 years (Paillard, 1998). Most species from temperate regions were affected by prolonged ice cover, which resulted in local extinctions, genetic bottlenecks, and range shifts. In contrast, warmer interludes were typically associated with demographic and range expansions (Hewitt, 2000). Even based on the least conservative molecular clocks available for CO1, the four major B. schaefferi clades identified in our study diverged before or during the first Pleistocene ice ages, around 6.8-1.7 mya. This suggests that at least four separate B. schaefferi refugia may have existed during the last glacial periods.
An alternative-and non-mutually exclusive-explanation could be that the clades were reproductively isolated prior to the ice ages.
It seems reasonable to assume that the clade A and C lineages would have been least affected by the consequences of glaciation.
Representatives from these clades all originate from localities in the Mediterranean and in northern Africa, areas far less affected by climate change throughout the Pleistocene than the more northern regions where B. schaefferi occurs today (Hewitt, 2000).  (Hewitt, 2000), including fairy shrimps (Muñoz et al., 2008;Reniers et al., 2013).
Branchipus schaefferi and Chirocephalus diaphanus are both widespread species in Europe. However, they differ in some ecological traits. While C. diaphanus is typically a cold-tolerant species and appears in fall and early spring, B. schaefferi is a thermophilic species, generally present during late spring and summer in Europe (e.g. Petrov & Cvetković, 1997). Phylogenetic reconstructions revealed genetic differentiation between C. diaphanus populations from eastern and western Europe (Reniers et al., 2013). Our reconstructions for B. schaefferi are largely consistent with this observation.
However, a number of shared haplotypes do occur among eastern and western regions (e.g. specimens from Serbia and Hungary present in clade D and the specimen from Spain in clade B). This is suggestive of recent long-distance dispersal events via vectors such as migrating water birds (Brochet et al., 2009;Green et al., 2005) and motorised vehicles (Waterkeyn et al., 2010).  (Brendonck, Pinceel, & Ortells, 2017;Vanschoenwinkel et al., 2013). These eggs act as propagules for passive dispersal and can easily be ingested by water birds or sporadically stick to their feathers (Sánchez, Hortas, Figuerola, & Green, 2012). Results from field studies demonstrate that eggs of the fairy shrimps can be dispersed across long distances in such a way (Brochet et al., 2009;Green et al., 2005;Lovas-Kiss et al., 2018;Rogers, 2014). It is also likely that longdistance dispersal of B. schaefferi is mediated by migratory birds. For instance, haplotype links between populations from Algeria and those in Belgium and Hungary in clade D match the yearly migration routes of some wader birds (Svensson, 2009 Combined, our results suggest that the studied B. schaefferi lineages probably represent four morphologically cryptic species, which correspond to the four clades in the reconstructed phylogeny. Therefore, the species should be subject to a taxonomic revision, based on combined information from morphological, genetic and ecological analyses.

| Indications for long distance dispersal
Clades A and C are especially vulnerable due to their restricted distribution in the Mediterranean region where their habitats are disappearing at an alarming rate (2015a; Van den Broeck, Waterkeyn, Rhazi, Grillas, et al., 2015b;Zacharias & Zamparas, 2010). Considering the current threats to B. schaefferi across its range of occurrence, we would for now at least like to promote the recognition of the four clades within the phylogeny as separate ESUs for directing conservation efforts.

| CON CLUS IONS
Overall, our results illustrate the importance of assessing the phylogeography of a species for the development of conservation strategies, especially for morphologically cryptic taxa with high genetic diversity. Temporary ponds contain many rare and specialist species, such as the studied freshwater crustacean. However, across the studied regions temporary ponds are also essential as sources of food, water, and breeding grounds to many other organisms including threatened birds and amphibians. Therefore, their protection should be considered a priority in nature management plans.

ACK N OWLED G M ENTS
The authors would like to thank Marko Šćiban, Christopher Rogers,

CO N FLI C T O F I NTE R E S T
Authors declare no conflict of interest.

DATA AVA I L A B I L I T Y
The DNA sequence data that support the findings of this study are openly available in GenBank at https ://www.ncbi.nlm.nih.gov/genba nk/, accession numbers are listed in Table 1.