Limited long-distance dispersal success in a Western European fairy shrimp evidenced by nuclear and mitochondrial lineage structuring.

Anostraca are known by their ability for long-distance dispersal, but the existence in several species of deep, geographically structured mtDNA lineages suggests their populations are subjected to allopatric differentiation, isolation, and prevalence of local scale dispersion. Tanymastix stagnalis is one of the most widespread species of Anostraca and previous studies revealed an unclear geographical pattern of mtDNA genetic diversity. Here, we analyze populations from the Iberian and Italian Peninsulas, Central Europe, and Scandinavia, with the aim to characterize the patterns of genetic diversity in a spatio-temporal framework using mtDNA and nuclear markers to test gene flow among close populations. For these aims we built a time-calibrated phylogeny and carried out Bayesian phylogeographic analyses using a continuous diffusion model. Our results indicated that T. stagnalis presents a deeply structured genetic diversity, including 7 ancient lineages, some of them even predating the Pleistocene. The Iberian Peninsula harbors high diversity of lineages, with strong isolation and recent absence of gene flow between populations. Dispersal at local scale seems to be the prevailing dispersal mode of T. stagnalis, which exhibits a pattern of isolation-by-distance in the Iberian Peninsula. We remark the vulnerability of most of these lineages, given the limited known geographic distribution of some of them, and the high risk of losing important evolutionary potential for the species.

The study of phylogeographic patterns of fresh water Anostraca has increased in the last decades (Ketmaier et al. 2005, 2008, 2012; Aguilar 2011; Reniers et al. 2013; Zarattini et al. 2013; Lindholm et al. 2016; Aguilar et al. 2017; Kappas et al. 2017;Rodríguez-Flores et al. 2017; Mioduchowska et al. 2018; Lukić et al. 2019) challenging our vision on how they disperse and colonize new environments and also improving our knowledge of the historical events affecting them. Although having potential for long-distance dispersal helped by their resistance eggs (Muñoz et al. 2014; Rodríguez-Flores et al. 2017; and references therein), the existence of deep, geographically structured mtDNA lineages indicates allopatric fragmentation and isolation and prevalence of local scale dispersion over long-distance colonization (but see Kappas et al. 2017; Rodríguez-Flores et al. 2017). On the other hand, these markedly structured phylogeographic patterns have been explained by founder effect and local adaptation followed by rapid growth of founding populations (Rogers 2015). The difficulties in the establishment of migrants due to demographic impediments (Monopolization Hypothesis) (De Meester et al. 2002) reinforce these patterns of genetic structure. Additionally, historical events such as Quaternary climate fluctuations may have left their mark in the genetic structure across current patterns of distribution, as has been observed for many Holarctic organisms (Verovnik et al. 2005; Recuero and García-París 2011; Reniers et al. 2013; Teixeira et al. 2018), including the fairy shrimp Tanymastis stagnalisLinnaeus, 1758 (Arukwe and Langeland 2013).

Tanymastix stagnalis is one of the most widespread species of Anostraca. It ranges from the British Islands eastwards to the Black Sea (Ukraine) and latitudinally from North Africa (Algeria) to Scandinavia (Sweden and Norway) (Young 1976; Vekhov 1991; Brtek and Thiéry 1995; Petkovski 1995; Samraoui and Dumont 2002; Sіdorovsкіy 2012; Arukwe and Langeland 2013; Zavadil et al. 2013). The species is typically found in granite or sandstone rock pools, dune or sand slacks, and road ditches, with oligotrophic water, poor mineralization, and high dissolved oxygen (Brtek and Thiéry 1995; Boix 2002; Olmo et al. 2015). Patterns of spatial distribution of related taxa (e.g., Chirocephalus diaphanusPrévost, 1803) are quite different from those shown by Tanymastix. In fact, T. stagnalis presents a patchy distribution, hypothetically as a consequence of its narrower ecological requirements. For instance, marked temperature fluctuations along its biological cycle can limit the viability of T. stagnalis populations in Mediterranean areas (Mura and Zarattini 2000; Zarattini et al. 2017). Specifically, in the Iberian Peninsula the species has been recorded from central and southwestern areas with some isolated populations in the northwestern and northeastern corners and a single coastal population in Valencia (Alonso 1985, 1996; Boix 2002; Olmo et al. 2015; Rodríguez-Flores et al. 2016; Sala et al. 2017). Even though about 10% of all known Anostraca are included in the IUCN Red List (García-de-Lomas et al. 2015), widely distributed species such as T. stagnalis are not generally considered as threatened (but see, e.g., Eder and Hödl 2002), even though T. stagnalis presents a highly fragmented distribution, with isolated populations and documented cases of local extinctions (Eder and Hödl 2002; Ketmaier et al. 2005; Arukwe and Langeland 2013; Olmo et al. 2015).

Ketmaier et al. (2005) studied the phylogeographic patterns of T. stagnalis in Western Europe and obtained an unclear geographical pattern of mtDNA genetic diversity suggesting that there is no gene flow among the analyzed populations. After that, Arukwe and Langeland (2013), using partial sequences of mitochondrial cytochrome b, found that Scandinavian populations were more closely related to Italian populations than to the single Spanish population studied. They concluded that the species likely survived during the last glaciation isolated in glacial refugia in Scandinavia. Here, we complemented these previous studies obtaining new mtDNA sequences from additional populations from the Iberian Peninsula in an attempt to characterize patterns of genetic diversity of the species in a spatio-temporal framework. Also, we test hypotheses on the effect of long-distance dispersal on phylogeographic patterns and the relationship between isolation by distance (IBD) and geographic structure across T. stagnalis lineages. We analyze geographically close Iberian populations of T. stagnalis, which theoretically should present high levels of gene flow between them, by using a nuclear marker (ITS2), already tested in phylogeographic studies of Anostraca (Rodríguez-Flores et al. 2017).

Material and Methods

Iberian populations of T. stagnalis were sampled from all habitat types previously described. Samples were obtained from 16 different locations in the Iberian Peninsula (Table 1). Specimens were collected using a small hand net, preserved in absolute ethanol, and stored at −20°C at the Museo Nacional de Ciencias Naturales (MNCN-CSIC).

Table 1.

Samples employed for this study, including MNCN DNA codes, locality, habitat typology, altitude, geographic coordinates, and GenBank accession numbers

MNCN: ADN	Locality	Habitat typology	Altitude (m)	Latitude	Longitude	cox1	ITS2	cox1 lineage	Mito-nuclear concordance
103441	Portugal: Faro: Lagõa da Nave A	Temporary small lake	155	37°13′07.0″N	8°02′53.9″W	MN190212	—	G	—
103442	Portugal: Faro: Lagõa da Nave A	Temporary small lake	155	37°13′07.0″N	8°02′53.9″W	MN190213	—	G	—
103448	Portugal: Faro: Lagõa da Nave A	Temporary small lake	155	37°13′07.0″N	8°02′53.9″W	MN190219	—	G	—
103449	Portugal: Faro: Lagõa da Nave A	Temporary small lake	155	37°13′07.0″N	8°02′53.9″W	MN190220	—	G	—
103443	Portugal: Faro: Lagõa da Nave B	Temporary small lake	155	37°13′09.1″N	8°02′58.5″W	MN190214	—	G	—
103444	Portugal: Faro: Lagõa da Nave B	Temporary small lake	155	37°13′09.1″N	8°02′58.5″W	MN190215	—	G	—
103420	Spain: Ávila: El Losar del Barco	Granitic rock pools	1,006	40°24′41.1″N	5°32′19.3″W	MN190191	MN190224	C	Yes
103421	Spain: Ávila: El Losar del Barco	Granitic rock pools	1,006	40°24′56.1″N	5°32′10.4″W	MN190192	MN190225	C	Yes
103422	Spain: Cáceres: Membrío	Road tracks in pasture, granite sands	345	39°33′44.6″N	7°04′32.1″W	MN190193	MN190226	D	Yes
103423	Spain: Cáceres: Membrío	Road tracks in pasture, granite sands	345	39°33′44.6″N	7°04′32.1″W	MN190194	—	D	—
103418	Spain: Cádiz: Puerto de Santa María	Dune slacks	5	36°32′39.6″N	6°12′18.3″W	MN190187		A	Yes
103417	Spain: Cádiz: Puerto de Santa María	Dune slacks	5	36°32′39.6″N	6°12′18.3″W	MN190188	—	A	—
103416	Spain: Cádiz: Puerto de Santa María	Dune slacks	5	36°32′39.6″N	6°12′18.3″W	MN190189	MN190222	A	—
103424	Spain: Ciudad Real: Las Virtudes	Flooded arable field, granite sands	876	38°34′04.2″N	3°24′40.7″W	MN190195	—	E	—
103426	Spain: Córdoba: Fuente Obejuna	Flooded arable field, granite sands	537	38°16′52.6″N	5°25′02.8″W	MN190197	MN190228	E	No
103427	Spain: Córdoba: Fuente Obejuna	Flooded arable field, granite sands	537	38°16′52.6″N	5°25′02.8″W	MN190198	—	E	—
103428	Spain: Girona: Sant Climent Sescebes	Road tracks	152	42°23′09.9″N	2°57′33.8″E	MN190199	—	E	—
103429	Spain: Girona: Sant Climent Sescebes	Road tracks	152	42°23′09.9″N	2°57′33.8″E	MN190200	—	E	—
103430	Spain: Huelva: Mazagón	Sand slacks, granite sands	52	37°10′09.0″N	6°47′26.5″W	MN190201	MN190229	A	Yes
103447	Spain: Huelva: Moguer	Sand slacks, granite sands	53	37°12′05.4″N	6°45′48.7″W	MN190218	—	A	—
103431	Spain: Madrid: La Cabrera	Granitic rock pools	1,072	40°51′35.9″N	3°37′39.1″W	MN190202	MN190230	C	Yes
103432	Spain: Madrid: La Cabrera	Granitic rock pools	1,072	40°51′26.0″N	3°37′38.7″W	MN190203	—	C	—
103419	Spain: Madrid: La Cabrera	Granitic rock pools	1,072	41°51′33.1″N	3°37′35.6″W	MN190186	—	C	Yes
103415	Spain: Madrid: La Cabrera	Granitic rock pools	1,072	41°51′33.1″N	3°37′35.6″W	MN190190	MN190223	C	—
103438	Spain: Madrid: Valdemaqueda	Granitic rock pools	1,005	40°29′39.3″N	4°18′59.4″W	MN190209	—	C	—
103439	Spain: Madrid: Valdemaqueda	Granitic rock pools	1,005	40°29′39.3″N	4°18′59.4″W	MN190210	—	C	—
103440	Spain: Madrid: Valdemaqueda	Granitic rock pools	1,005	40°29′39.3″N	4°18′59.4″W	MN190211	—	C	—
103433	Spain: Sevilla: Coria del Río	Sand slacks	26	37°14′32.2″N	5°59′43.5″W	MN190204	MN190231	A	Yes
103434	Spain: Sevilla: Coria del Río	Sand slacks	26	37°14′32.2″N	5°59′43.5″W	MN190205	MN190232	A	Yes
103435	Spain: Toledo: Calera y Chozas	Road ditches, granite sands	409	39°55′39.8″N	5°01′12.9″W	MN190206	MN190233	B	Yes
103436	Spain: Toledo: Calera y Chozas	Road ditches, granite sands	409	39°55′39.8″N	5°01′12.9″W	MN190207	MN190234	B	Yes
103425	Spain: Toledo: Talavera de la Reina	Flooded arable field, granite sands	403	39°56′53.1″N	4°56′40.9″W	MN190196	MN190227	B	Yes
103437	Spain: Toledo: Talavera de la Reina	Flooded arable field	404	39°56′53.1″N	4°56′45.0″W	MN190208	MN190235	B	Yes
103445	Spain: Valladolid: Llano de Olmedo	Granite sand slacks	762	41°16′06.7″N	4°36′14.5″W	MN190216	MN190236	C	Yes
103446	Spain: Valladolid: Llano de Olmedo	Granite sand slacks	762	41°16′06.7″N	4°36′14.5″W	MN190217	MN190237	C	No
	France: Corse-du-Sud: Bonifacio			41°22′49″N	9°12′43″E	AY555262	—	F	—
	France: Île-de-France: Fontainebleau			49°00′51″N	0°05′34″W	AY555254	—	E	—
	Germany: Brandenburg: Rühstädt			52°58′37″N	9°56′17″E	AY555253	—	E	—
	Italy: Toscana: Capraia Island			43°02′17″N	9°49′35″E	AY555265	—	F	—
	Italy: Umbria: Forca Canapine			42°45′07″N	13°12′01″E	AY555261	—	E	—
	Italy: Valle d′Aosta: Gran Paradiso			45°40′22″N	7°20′26″E	AY555256*	—	C	—
	Italy: Abruzzo: Lago dell′Orso			42°24′21″N	13°14′47″E	AY555260	—	E	—
	Italy: Abruzzo: Piani di Fugno			42°22′09″N	13°27′35″E	AY555258	—	E	—
	Italy: Abruzzo: Piano di Raccollo			42°23′37″N	13°39′27″E	AY555259	—	E	—
	Italy: Sardegna: Romazzino			41°05′22″N	9°33′47″E	AY555263	—	F	—
	Italy: Sardegna: Carloforte			39°09′08″N	8°13′57″E	AY555264	–	F	—
	Italy: Emilia-Romagna: Val Cedra			44°22′09″N	10°05′28″E	AY555257	—	E	—
	Spain: Madrid: Sierra Guadarrama			40°51′37″N	3°41′04″W	AY555255*	—	E	—
	Norway: Sør-Trøndelag: Trollheimen			62°11′15″N	9°55′17″E	AY555252	—	E	—

Em dashes (—) mean no data available. Locality of GenBank sequences with * seems to be mislabeled and were treated accordingly.

Tissue for the molecular analyses was obtained from the first thoracopods. DNA was extracted according to the protocols described by Vörös et al. (2016). The barcoding fragment of the mitochondrial Cytochrome c oxidase subunit I gene (cox1) was amplified by polymerase chain reaction (PCR) employing the primers LCO1490 (Folmer et al. 1994) and COI-H (Machordom et al. 2003). The nuclear ITS2 fragment was also amplified for some of the populations, using the pair of primers CAS28sB1d/CAS5p8sFt (Ji et al. 2003). PCR reactions were performed in a total volume of 25 µL, including one unit of DNA polymerase (Biotools, 5 U/µL), 2.5 mM of each primer, 0.4 mM of dNTPs, 1.5 mM of MgCl2, and 67 mM of a reaction buffer (Tris–HCl, pH = 8.3, Biotools). PCR reactions typically consisted of 35 cycles with a denaturing temperature of 94°C (30 s), annealing at 42–45°C (45 s), and extension at 72°C (90 s). PCR products were visualized in a 1.5% agarose gel, purified with sodium acetate precipitation, and sequenced at Macrogen Inc. (Macrogen Spain) using a AppliedBiosystem 3730XL DNA analyzer.

Sequences were edited and assembled using Sequencher v.5.4.1, then aligned with MAFFT (Katoh and Toh 2008), and the final alignments were corrected manually in Mesquite v 3.51 (Maddison and Maddison 2018). To identify possible pseudogenes, cox1 sequences were translated to look for internal stop codons, and nucleotide alignments were revised to discard the existence of indels. Available sequences of cox1 of T. stagnalis were obtained from GenBank (GenBank accession numbers: AY555238–AY555265) (Ketmaier et al. 2005) and included in the analyses. Table 1 shows sampling localities, and accession numbers of the sequences of cox1 and ITS2 generated in this study. A sequence of Tanymastix affinis from Ouezanne (Morocco) was used as outgroup (GenBank accession numbers: MN190221, MN190238).

Divergences between lineages and within lineages were calculated using uncorrected p-distances with MEGA7 (Kumar et al. 2016). Haplotypic, gene and nucleotidic diversity were calculated with DnaSP v. 5.10 (Rozas et al. 2003). We used Mantel’s tests (Mantel 1967) to test potential existence of “IBD” processes among populations from central and southern Spain, considering both mtDNA and nuclear data. Mantel’s tests were performed with GENALEX v.6 (Peakall and Smouse 2006), with 999 permutations to estimate the 95% upper tail probability of the matrix correlation coefficients.

Phylogenetic relationships were estimated using coalescence Bayesian inference with BEAST v1.10.2 (Drummond et al. 2012) and maximum likelihood (ML). ML analyses were performed using RAxMLGUI v.1.5 platform (Silvestro and Michalak 2012), calculating bootstrap support values with 1000 pseudo replicates using the GTR substitution model.

Coalescence Bayesian phylogenetic inference and time divergence estimates for mitochondrial and nuclear lineages were calculated using BEAST v1.10.2. There is not fossil available for contemporary Anostraca to calibrate the molecular clock, so we used a rate of substitution which provides a tentative estimate of divergence times (Reniers et al. 2013; Eimanifar et al. 2015; Lindholm et al. 2016). First, substitution model that fits best our data was estimated using JmodelTest (Posada 2008). Second, we ran and compared 2 analyses in BEAST in order to select the best clock model. We compared 2 clock models: a relaxed uncorrelated lognormal clock model and a strict clock model. We estimated the marginal likelihood of each clock model with a stepping stone and path sampling analysis in BEAST (Baele et al. 2016) and selected the best model with a Bayes factor test as was performed in Rodríguez-Flores et al. (2019). We implemented a GTR substitution model; a Bayesian skyline tree prior, which considers past population dynamics based on the coalescent model (Drummond et al. 2005) and a molecular clock rate with an uniform distribution from 0.013 to 0.007. The analyses were run for 200 million generations sampling every 20,000, and repeated independently to check consistency of results. We assessed the convergence of effective sampled size (ESS) in Tracer v1.7 (Rambaut et al. 2018). Once analyses reached convergence and ESS higher than 200, the first 25% of the initial trees were discarded as burn-in. Support of the nodes was estimated with posterior probabilities (pP), considering high support values of 0.95 or higher. In this study, clades with moderate pP (pP = 0.8) are only considered for discussion if they have bootstrap values higher than 70.

Continuous diffusion analysis was performed with BEAST 1.10.2, using cox1 dataset and the same configuration as for the time-calibrated phylogenetic reconstruction (i.e., GTR as substitution model, strict clock with a molecular clock rate with uniform distribution from 0.013 to 0.007, and a Bayesian Skyline tree prior). A location trait was generated using geographical coordinates of each sample and a Cauchy RRW model was used in order to infer geographical diffusion processes through time across branches of the inferred tree (Lemey et al. 2010). The analysis was run for 200 million generations and repeated independently to check consistency of results. Convergence was assessed using Tracer. A maximum clade credibility tree was reconstructed with TreeAnnotator applying a 25% burn-in and used to generate the reconstruction of the diffusion process using the module “Continuous Tree” in Spread 1.0.7 (Bielejec et al. 2011).

Results

The alignments included a total of 43 individuals for cox1 (658 bp), including new data as well as some sequences from GenBank, and 16 individuals for the ITS2 (713 bp) (Table 1). The results of the Bayes factor test supported the strict clock model (1) vs. relaxed lognormal model (2): logmarginal likelihood under path sampling model (1) = −2901.837 and (2) =−2902.940; logBF (1) = 0 and (2) =−2.206; p (1) = 0.901 and (2) = 0.099.

Six main mitochondrial lineages of T. stagnalis are recovered as monophyletic with high pP and bootstrap support values (Figure 1) although relationships among them are not always fully resolved. Lineage A includes populations from southern Iberian Peninsula (El Puerto de Santa María, Cádiz; Mazagón, Huelva; Coria del Río, Sevilla); lineage B includes populations from the central area of the Tagus River basin (Toledo); lineage C groups populations from the Sistema Central Mountain Range and the Northern Iberian Plateau (Sierra de Guadarrama in Madrid, Ávila, and Valladolid); lineage D includes a single population from Western Iberian Peninsula (Membrío, Cáceres); lineage E is formed by 4 subclades: E1 formed by a single, isolated population from the southernmost edge of the Iberian Southern Plateau (Sierra Morena, Ciudad Real); E2, including a population from Fuente Obejuna, Córdoba; E3, with all populations from Central Europe (France and Germany) and northern Italy nested together with the population from northeastern Iberia (San Climent Sescebes, Girona), and E4 incorporating Central Italian and Scandinavian haplotypes; lineage F incorporates all populations from the Tyrrhenian Sea islands (Corsica, Sardinia, Capraia); and lineage G is formed by populations from southern Portugal. Lineages D and C are recovered as sister groups with high pP (pP = 0.96), but this relationship lacks bootstrap support. Lineage F appears as sister to a clade including lineages A, B, C, D, and E, with moderate support (pP = 0.92, not supported by ML analysis). However, the placement of lineage G as the sister group of a clade including all other lineages is highly supported. ITS2 tree topology (including only Iberian populations) was congruent with cox1 tree topology, with the following exceptions: Llano de Olmedo (Valladolid) (an17239) which is clustered together with lineage B instead of being within lineage C (pP = 0.9), and Fuente Obejuna (Córdoba), which is clustered with lineage A instead of lineage E (pP = 1) (Figures 1 and 2). ML analyses recovered poor bootstrap support for most lineages (Figure 2), with the exception of lineage C from Central Spain.

Figure 1.

mtDNA chronogram obtained from BEAST, including posterior probabilities (values under 0.85 are not displayed)/ML bootstrap support values, high posterior density (HPD95%) intervals of the time to most recent common ancestor (TMRCA), and map showing geographic distribution of the mtDNA lineages. Hyphens (-) indicate ML boostrapt support <65.

Figure 2.

nDNA tree (ITS2) obtained from BEAST, including posterior probabilities (values under 0.85 are not displayed) and lineages. Asterisks indicate nodes supported by ML analysis. Pictures represent Tanymastix stagnalis (male and female) habitus and typical habitat for each lineage.

Genetic divergences (uncorrected p-distances) ranged from 2.22% to 7.95% in cox1 and from 1.12% to 1.86% in the ITS2 among the analyzed lineages. Distances within lineages ranged from 0 to 0.83% and from 0 to 2.39% in ITS2 and cox1, respectively (Table 2). Mantel test’s for “IBD” indicated significative correlation between geographical and genetic distances among populations from central and southern Spain, both for cox1 (P = 0.001) and ITS2 (P = 0.034) data (Figure 3). However, the signal of IBD is much weaker in the ITS2 dataset, probably because of its much lower substitution rate. Haplotypic, gene and nucleotidic diversity, number of specimens per lineage, and mean intralineage P-distances are summarized in the Supplementary Material.

Table 2.

Genetic distances (uncorrected P-distance) in average percentages between the recovered lineages of Tanymastix stagnalis for the cox1 (values under the diagonal) and ITS2 (values above the diagonal)

LINEAGE	A	B	C	D	E1	E2	E3	E4	F	G
A	1.07/0.54	1.12	1.64	1.72	—	—	—	—	—	—
B	4.36	0.41/0.57	1.72	1.86	—	—	—	—	—	—
C	4.22	3.34	0.58/0.83	1.57	—	—	—	—	—	—
D	5.19	4.41	3.60	0/0	—	—	—	—	—	—
E1	4.88	4.18	3.73	4.86	0/—	—	—	—	—	—
E2	4.28	3.34	3.08	4.56	2.74	0.91/—	—	—	—	—
E3	4.86	4.54	4.19	5.06	3.19	2.86	2.39/—	—	—	—
E4	4.17	3.77	3.98	4.82	3.01	2.22	2.81	0.58/—	—	—
F	5.13	4.62	4.76	5.00	5.00	4.74	5.03	4.49	1.42/—	—
G	7.95	6.65	6.63	7.17	7.00	7.41	7.37	7.13	6.78	0.51/—

Numbers in bold show cox1 and ITS2 intra-lineage distances. Em dashes (–) mean no data.

Figure 3.

Isolation-by-distance plot of uncorrected P-distances versus geographic distance (km) with all pairwise combinations within the Tanymastix stagnalis populations from Central to Southern Spain.

According to our estimates, the time to most recent ancestor (TMRCA) of the mtDNA lineages of T. stagnalis is placed at 4.38 Ma [high posterior density (HPD95%): 2.95–6.43 Ma] when the splitting of lineage G from the others occurred. Lineage F split from the others at 2.93 Ma (HPD95% 2.00–4.27). TMRCA of the rest of lineages is set at 2.49 Ma (HPD95%: 1.73–3.62 Ma). TMRCA for lineages C and D is 1.69 Ma (HPD95%: 1–2.62 Ma). The lineage E1 split from the others within the clade E at 1.64 Ma (HPD95%: 1.07–2.46 Ma); populations from Scandinavia are estimated to split from those from Central Italy at 0.55 Ma (HPD95%: 0.21–0.95 Ma). Within clade F, the lineage from Capraia split from the rest of the Tyrrhenian islands at 1.56 Ma (HPD95%: 0.8–2.55 Ma). TMRCA estimates for the rest of clades are more recent: (A) 0.72 Ma (HPD95%: 0.44–1.49 Ma), (C) 0.57 Ma (HPD95%: 0.26–1 Ma), (G) 0.51 Ma (HPD95%: 0.18–0.99 Ma), (B) 0.32 Ma (HPD95%: 0.08–0.66 Ma), (D) 0.02 Ma (HPD95%: 0–0.13 Ma) (Figure 1).

The continuous diffusion analyses suggest ancestral areas for T. stagnalis at around 4.5 Ma located in the southwestern and central regions of the Iberian Peninsula, with subsequent expansions at the beginning of the Pleistocene in central Spain, and also to the Tyrrhenian Islands and the Italic Peninsula. From this later area, the species would have expanded into different areas in northeastern Iberia and Central Europe, in a progressive range expansion. According to the analyzed data, colonization of Scandinavia would be the result of a northern expansion of the central Italian clade in very recent times (Figure 4).

Figure 4.

Continuous diffusion reconstruction of the fairy shrimp T. stagnalis, showing ancestral distribution and range expansion of the obtained mitochondrial lineages. (A) 4.5 Ma (million years old), (B) 1.4 Ma; and (C) 0.09 Ma.

Discussion

Phylogeographic history

Our study unveils a high previously unnoticed genetic diversification in T. stagnalis. Cox1 data exhibit a clear, geographically structured genetic diversity formed by 6 main lineages, some of them presenting also some degree of substructure. All of these divergent lineages, except the Tyrrhenian clade are represented within the Iberian Peninsula, some of them restricted to small geographic areas. This structure was not detected previously since studies on the phylogeography of T. stagnalis focused on populations from the Tyrrhenian islands, the Italic Peninsula, and central and northern Europe, including a single population from the Iberian Peninsula (Ketmaier et al. 2005; Arukwe and Langeland 2013).

The mtDNA phylogeographic structure is mirrored by the nuclear ITS2 data, which allows to reject the existence of effective gene flow among the characterized lineages from central and southern Iberia. Only populations from Llano de Olmedo, Valladolid; and Fuente Obejuna, Córdoba, presented ITS2 alleles that were not grouped together in a single lineage, which could be a consequence of gene flow or introgression (Perea et al. 2016), but also of the multicopy nature of this locus (Yao et al. 2010).

Most studies on the phylogeographic patterns of Anostraca are based in mtDNA sequence data, repeatedly showing deep mtDNA lineages (Ketmaier et al. 2005, 2008, 2012; Aguilar 2011; Reniers et al. 2013; Eimanifar et al. 2015; Lindholm et al. 2016; Kappas et al. 2017). However, existence of gene flow and introgression among populations from different lineages, for instance after range expansions leading to secondary contact (Recuero et al. 2014), cannot be tested in absence of rapid nuclear marker data. Our results from the ITS2 marker demonstrated recent isolation and absence of gene flow between Iberian populations of Branchinecta feroxMilne-Edwards, 1840 (Rodríguez-Flores et al. 2017). Similarly, the high degree of genetic structure found in T. stagnalis indicates a history of strong geographic isolation of populations maintained through time. There are no obvious geographic barriers to dispersal that could be preventing gene flow among lineages. Alternatively, lineage extinction and independent evolutionary histories in glacial refugia or sanctuaries during the Pleistocene could have shaped these patterns of diversity (Gómez and Lundt 2007; Recuero and García-París 2011). The estimated TMRCA of the main T. stagnalis lineages suggests ancient processes even predating the Pleistocene glacial episodes, with the earliest differentiations occurring during Mid to Late Pliocene and therefore with a long evolutionary history subjected to lineage replacement or extinction. Main T. stagnalis lineages and sub-lineages persisted mostly in the Iberian Peninsula, but also the Italic Peninsula and the Tyrrhenian islands, probably during the last 2.4 Ma in a context of climatic oscillations. This scenario usually leads to phylogeographic patterns known as “refugia within refugia” (Gómez and Lunt 2007), in which successive extinctions and range contractions promote isolation and allopatric differentiation (García-París et al. 2003; Ribera and Vogler 2004; Teixeira et al. 2018), contrasting with the patterns of low genetic diversity with no phylogeographic structure due to bottlenecks and rapid range expansion, mostly observed outside Pleistocene refugia and sanctuaries (Hewitt 2001; Recuero and García-París 2011; Kappas et al. 2017).

The Iberian Peninsula may have acted as a center of lineage accumulation, harboring populations from northern areas, that were progressively colder at the end of the Pliocene. The Peninsula acted then, either as a sanctuary, incorporating the northern genetic diversity to the diversity of the native Iberian populations (Recuero and García-París 2011), or as a refugium with multiple refugia (sensuGómez and Lunt 2007) with each lineage relatively isolated initially. Successive glacial cycles could have either promoted isolation or alternatively the intermixing of lineages (Baselga et al. 2011). According to our phylogeographic data, lineages within the Iberian Peninsula remained isolated from each other for a long period (Figure 1), supporting thus the idea of strong geographic isolation during the Pleistocene cycles. The high diversity of T. stagnalis lineages and their geographic isolation in the area studied could be consequence of 2 alternative scenarios: (1) isolation maintained by different ecological niches or (2) isolation maintained by limited dispersal abilities (genetic isolation in geographic areas with independence of the ecological niche). We could choose the first hypothesis over the second if the niche occupation would have phylogenetic signal, but in T. stagnalis there is no evidence of monophyly for habitat preferences, as has been observed in other branchiopods that lack genetic segregation depending on habitat type (Ventura et al. 2014). For instance, in spite of a clear preference of the species for acid and granitic soils, populations inhabiting strictly in granitic rock pools appeared in different lineages (Sistema Central and Scandinavia) (Arukwe and Langeland 2013), as happens with occupation of dune or sand slacks. Moreover, the antiquity of the lineages suggests long persistence of the populations in each area, with high degree of local isolation, and high genetic diversity.

Local versus long-distance dispersal patterns and considerations for conservation

Tanymastix stagnalis presents a patchy population distribution with lineages present in small, defined regions, with scattered records and discontinuous distribution across its global range (Brtek and Thiéry 1995). These observations are coincident with the genetic patterns of diversity shown in our study, including the positive correlation between geographic and genetic distances. According to these data, it seems that local dispersal is the prevailing mode of colonization for T. stagnalis, and that long-distance dispersal events are uncommon, if not absent at least in the Iberian Peninsula.

Other species of Anostraca present patterns of dispersal at a local scale coincident with isolation-by-distance models (Hulsmans et al. 2007). However, according to our continuous diffusion analyses, we cannot discard that the colonization of regions outside the Iberian and Italic peninsulas could have been the result of long-distance colonization events. Alternatively, current distribution patterns could be explained by a continuous range expansion in times of more favorable ecological conditions, followed by extinction in large regions currently not occupied by the species. For instance, our data suggest an event (likely in the last 500,000 years) of long-distance colonization from central Italy to Scandinavia. This reconstruction is compatible with the scenario of coastal glacial refugia in Scandinavia during the Last Glacial Maximum proposed by Arukwe and Langeland (2013), but to test this hypothesis it would be necessary to include populations from the Balkans and Eastern Europe, since we cannot exclude a more recent colonization of Scandinavia from potential eastern refugia, and not necessarily by a long-distance dispersal event but as a continuous range extension from yet uncharacterized refugia.

Tanymastix stagnalis constitutes a steno-topic organism, with strict ecological requirements as its preference for oligotrophic and cold waters, being especially sensible to sudden changes and high temperatures (Olmo et al. 2015). Tanymastix stagnalis presents particular resting eggs that harbor high embryo volume which would permit them to survive long drought periods (Mura and Zarattini 2000; Thiéry et al. 2007; Zarattini et al. 2017). Mura and Zarattini (2000) indicated that T. stagnalis biological cycle is adapted to longer hydroperiods than other species; this may explain the presence of the species in high rainfall areas like Ireland or northwestern Iberian Peninsula, where T. stagnalis is the only anostracan present. However, some populations occur in ephemeral ponds with very short hydroperiods (for instance, rock pools in La Cabrera, Madrid) (personal observation), which could represent local adaptations.

Local-scale dispersal observed among Iberian populations, steno-topic ecological requirements, together with a strong phylogeographic pattern and scattered isolated populations distributed across the western Palearctic (Brtek and Thiéry 1995), indicate a complex evolutionary history. A likely scenario to explain the observed phylogeographic pattern and current distribution is that the species had its optimum during glaciations, when colder temperatures allowed for a larger range with a continuous distribution. Extant populations would be located currently in interglacial Holocene refugia, as has been observed in other organisms favoring cold temperatures (Teixeira et al. 2018; Sánchez‐Montes et al. 2019), having suffered recurrent periods of population extinction and range contraction. This would favor the hypothesis of continuous range expansion followed by extinction as the cause of the current pattern of isolated populations along the distribution range of T. stagnalis over the alternative hypothesis of long-distance colonization. The incorporation of genetic data from eastern Palearctic populations, and comparative phylogeographic studies including co-distributed anostracan, would help to clarify the evolutionary history, colonization patterns, and dispersal mechanisms of the species, and to support the hypotheses here presented.

The eggs of Anostraca are dispersed passively through different vectors like wind, birds, flying insects, mammals, amphibians, other crustaceans, and humans (Hulsmans et al. 2007; Beladjal and Mertens 2009; Muñoz et al. 2014; Rogers 2014), being in some cases an effective way to colonize habitats far away from the original population (Muñoz et al. 2014; Rodríguez-Flores et al. 2017). However, our study adds to many other that revealed a strong decoupling between dispersal potential and complex genetic structure (Martínez‐Solano et al. 2006; Korn et al. 2010; Recuero and García-París 2011; Marrone et al. 2013), leading to reject the Everything is Everywhere Hypothesis (Incagnone et al. 2015).

The observed strong isolation of T. stagnalis lineages, even at regional scale, indicate a strong vulnerability to risk factors, both at local and global scales (habitat deterioration, global warming) (Arukwe and Langeland 2013; Olmo et al. 2015). The effect of global warming might be particularly strong for Iberian populations, where highest diversity is accumulated. In this area, droughts are predicted to be more frequent and generalized (Grillakis 2019), menacing the viability of T. stagnalis populations. Indeed, those populations represent long-term reservoirs of species genetic diversity and, consequently, sources of evolutionary and speciation potential (Hampe and Petit 2005; Zavadil et al. 2013; Mee et al. 2015). Given the patchy distribution pattern of T. stagnalis, with several lineages known from few, isolated populations, and taking into account the species reduced dispersal ability, we consider that the species is threatened of losing a large part of its evolutionary legacy and potential.

Limitations of this work derive from the relatively small sample sizes used and from the absence of material for study from the easternmost areas of the distribution of the species, precisely where it seems to be in close parapatry with Tanymastix motasiOrghidan, 1945 (Petkovski 1995). Time divergence estimations should be taken cautiously because we have used a general sequence divergence rate for Anostraca in absence of fossil data (Lindholm et al. 2016; Lukić et al. 2019), and also, they might be affected by sampling bias and density (Wang et al. 2017). Similarly, some phylogeographic patterns, particularly those concerning the relationships of the Scandinavian populations, might be veiled because of the lack of eastern European specimens. Additional sampling at a local scale, and an increased number of specimens per population, should be analyzed to identify secondary contact zones and possible introgression across lineages. However, our work includes a large number of previously unknown populations, some of them representing completely new lineages (Toledo, Cáceres, Ciudad Real…); in this respect, historical absence of adequate geographic data, obscures the true distribution of T. stagnalis, which currently presents enormous, real or not, distribution gaps (most of France south to Fontainebleau, for example; Rabet 1994). Accordingly, ecological studies to identify niche clade preferences are still premature.

Click here for additional data file.