Large effective size as determinant of population persistence in Anostraca (Crustacea: Branchiopoda)

The fairy shrimp Branchinectella media, because of its passive dispersal capacity and scarce and irregularly distributed habitats (temporary saline aquatic systems), is an intriguing organism from a population genomics and conservation perspective. Stochasticity of dispersal events and the irregular distribution of its habitat might lead to low levels of population connectivity and genetic diversity, and consequently, populations with limited persistence through time. Indeed, by using genomic data (SNPs), we found a strong genetic structure among some of the geographically isolated Iberian populations of B. media. Interestingly, we also obtained high estimates of effective population sizes. Lack of suitable habitat between populations (absence of a “stepping stone” network) and strong genetic differentiation suggest limited dispersal success in B. media. However, the high effective population sizes observed ensure persistence of B. media populations against genetic stochasticity (genetic drift). These results indicate that rescue-effect might not be essential for population persistence if they maintain high effective population sizes able to hold adequate levels of genetic diversity. Should high population sizes be reported in other low dispersing Anostraca, one might be optimistic with regard to their conservation status and fate, provided that their natural habitats remain undisturbed.


Introduction
Genetic diversity is a key factor for species persistence, as it allows individuals to face selection and adapt to environmental changes (Ceballos et al. 2017;Frankham 1996).Genetic variability in a population may decrease through time as a result of non-random mating, genetic drift and/or natural selection (Newman and Pilson 1997; but see Reed 2007).
The intensity of genetic drift increases in small and isolated populations (mutations usually have a negligible effect in allele frequencies at short timescale).Indeed, loss of genetic diversity in small and isolated populations may be enhanced by subsequent inbreeding (Crnokrak and Roff 1999;Hedrick and Kalinowski 2000;Jangjoo et al. 2016;Star and Spencer 2013;Wright 1931).Species with fragmented distributions, either due to natural processes like long-distance colonisations (Green et al. 2005;Recuerda et al. 2021;Sánchez-Vialas et al. 2021) and geologic/climatic events (García-Marín et al. 1999;Saunders et al. 1991), or due to human induced habitat destruction (Keller and Largiadèr 2003;Templeton et al. 1990), are therefore, exposed to genetic diversity loss to a larger extent because of connectivity hampering.In contrast, genetic drift will be weak in large populations (with numerous reproductive individuals) because they are expected to harbour more genetic diversity and transfer it through generations to a larger extent.In addition, a larger number of breeders will lead to more mutation events (Méndez et al. 2011;Rich et al. 1979).The number of breeding individuals in a population represent the effective population size (N e ) (Wright 1938).The most frequent definition of N e is the size of an ideal population that would experience the level of genetic drift or inbreeding observed in the real population.Other calculation procedures gave rise to eigenvalue, mutation and coalescence effective sizes (Caballero 2020).Therefore, isolated populations with large contemporary N e maintain a high level of genetic diversity in spite of limited gene flow (Palstra and Ruzzante 2008), and will retain their evolutionary potential (Franklin and Frankham 1998).
Dispersal and connectivity between populations in species with geographically disjunct habitats may be further limited when they disperse passively, depending thus on external abiotic and biotic agents like meteorological phenomena or animal vectors (Beladjal and Mertens 2009;Coughlan et al. 2017;Green et al. 2005;Vanschoenwinkel et al. 2008).This may lead to species with very strong genetic structure due to geographic isolation (Makino and Tanabe 2009;Mills et al. 2007;Muñoz et al. 2008;Rodríguez-Flores et al. 2019;Xu et al. 2009) and unlikely rescue effect (Vilà et al. 2003).In the absence of migration, the persistence of populations would thus rely directly on a large N e .The estimation of N e is particularly important to evaluate the conservation status of isolated populations, being one of the main objectives in conservation genetics/genomics (Benestan et al. 2016;Carreras et al. 2017;Hohenlohe et al. 2021;Plough 2017).
Fairy shrimps (Branchiopoda: Anostraca) are invertebrates that disperse passively across geographically fragmented populations (García- de-Lomas et al. 2015).They are small continental aquatic crustaceans that inhabit temporary ecosystems.In Anostraca, connectivity among populations depends on both habitat availability (Sainz-Escudero et al. 2022) and stochasticity of dispersal events (Kappas et al. 2017;Rodríguez-Flores et al. 2017), often mediated or directed by vertebrates (Rogers 2014a;Sainz-Escudero et al. 2022).Species whose habitats are more continuously distributed have better chances to maintain a connected metapopulation system using a stepping-stone network, keeping higher levels of genetic exchange (Bellin et al. 2020;García-de-Lomas et al. 2017;Hanski and Simberloff 1997).An increase of connecting success might occur when anostracan eggs pass through the digestive tract of dispersing vertebrates (Rogers 2014a).In contrast, less common and distant habitats may result in high inter-population genetic differentiation and strong phylogeographic structure due to the limited connectivity and gene flow, and, consequently, in threatened species (Ketmaier et al. 2008;Lukić et al. 2019;Rodríguez-Flores et al. 2019; but see Kappas et al. 2017).One species likely subjected to these conditions is Branchinectella media (Schmankewitsch 1873), which inhabits temporary saline lakes, ponds and puddles present in steppes scattered along the Palaearctic region (Alonso 2010;Sainz-Escudero et al. 2019) (Fig. 1).In the Iberian Peninsula, this kind of habitat is patchily distributed in separated complexes along the eastern and central areas of Spain (Alcorlo et al. 1997;Alcorlo and Baltanás 1999;García et al. 1997;Pons et al. 2018;Rodríguez-Flores et al. 2016;Sala et al. 2005).These saline lakes are protected under the EU Habitats Directive (Council Directive 92/43/ EEC of 21 May 1992 on the conservation of natural habitats and of wild fauna and flora).
Since fairy shrimp conservation, and B. media in particular, depends on knowledge of population demography (i.e., which lakes present long term viable N e ), in this work we used single nucleotide polymorphisms (SNPs) generated by MobiSeq (Rey-Iglesia et al. 2019) to genotype 103 individuals of B. media sampled at five Iberian saline ponds/lakes, and develop a population and conservation genomics study.Our specific aims were (i) to investigate the geographic distribution and levels of genetic diversity in the Iberian populations of B. media and evaluate the existence of isolation, (ii) to estimate contemporary N e of the resulting genetic units and (iii) to assess the viability and relevance of N e for B. media populations persistence and the effectiveness of this evolutionary strategy in Anostraca species with limited dispersal capacities.

Illumina library preparation and sequencing
High quality genomic DNA was extracted from thoracopods of the 120 individuals of B. media using a Zymo Research Quick-DNA Microprep Plus Kit (Zymo Research Corporation, Portugal) following the manufacturer's protocol.DNA quantity was assessed by using a Qubit Fluorometer (Thermo Fisher Scientific®) and 4 individuals (1 from Altillo and 3 from Alcahozo) did not achieve the minimum DNA amount required for sequencing.All laboratory processing from the low coverage whole genome sequencing of one individual and candidate region selection of a transposable element (TE) to libraries sequencing was carried out by All Genetics & Biology SL (www.allge netics.eu) following the reduced representation library protocol of MobiSeq (Rey-Iglesia et al. 2019).Once the TE was chosen, the conserved region of the TE was used to design the primer MobiBranchi_578R (5′GAA GCT CAT TCT CGT AGG C3′), which was selected for the MobiSeq experiment.Libraries and amplification of selected genomic regions were performed.The resulting libraries were quantified and sequenced in a fraction of a NovaSeq PE150 flow cell.

Raw data processing, SNP genotyping and filtering
Bioinformatic processing from raw data curation to SNP calling was also carried out by AllGenetics & Biology SL.Illumina paired-end raw data (between 1,358,794 and 5,123,706 paired-end reads per individual/library) consisted of forward (R1) and reverse (R2) reads.Quality of raw fastq files was checked in FastQC (v.0.11.5) (Andrews 2010) and summarized in MultiQC (Ewels et al. 2016).The clumplify.sh script from the BBmap package 38.90 (Bushnell B.https:// sourc eforge.net/ proje cts/ bbmap/) was used to remove exact forward and reverse duplicated reads.Trimmomatic 0.39 (Bolger et al. 2014) was employed to remove the adapters (ILLUMINACLIP option) and low-quality regions (AVGQUAL:26).Cutadapt 2.9 (Martin 2011) was used to remove reads that did not contain the primer sequences (allowing up to 1 mismatches).De novo assembly of loci (reference catalogue construction), read mapping and SNP calling was performed through dDocent pipeline (Puritz et al. 2014).Previous difficulties in the primer selection, and the high levels of obtained PCR duplicates indicated high levels of nucleotide polymorphisms.This led to the use of dDocent script, which takes full advantage of both paired-end reads and performs a quality trimming instead of the quality filtering carried out by STACKS (Catchen et al. 2011(Catchen et al. , 2013)).
The decomposed VCF (variant calling format) file obtained was filtered using VCFtools v4.1 (Danecek et al. 2011).Indels (insertions and deletions) were removed and bi-allelic SNPs were kept, as well as loci with a missing data proportion equal or less than 15% and with a minimum allelic frequency (MAF) of 0.02.There were 13 individuals with more than 75% of missing data for all the loci, so they were discarded keeping a total of 103 individuals (15 from Longar, 12 from Altillo, 23 from Alcahozo, 23 from Tirez and 30 from Piñol).In order to detect loci under balancing or purifying selection (outlier loci) and remove them from the neutral dataset, we used the Bayesian-based software BayeScan v2.1 (Foll 2012;Foll and Gaggiotti 2008).We ran the program with the default 5000 iterations (sample size), with a thinning interval of 200, 20 pilot runs of 10,000 iterations each, and a burn-in of 100,000.Prior odds for the model were set at 100.After the convergence checking, outlier loci were selected through a false discovery rate of 0.1, obtaining 13,689 neutral SNPs and 98 outliers.In order to remove library and sequencing artefacts or allelic dropout (Waples 2015), we removed loci departing from Hardy-Weinberg Equilibrium (HWE) following the "out all" approach, which avoids the removal of informative loci for population structure through the choice of a locus as departed from HWE when it did in all populations (Pearman et al. 2021).We used the function gl.hwe.pop of the R package "HardyWeinberg" 1.7.5 (Graffelman 2015), and 151 loci were detected at disequilibrium in each of the five populations.Finally, our dataset consisting of 13,538 loci grouped in 1771 contigs, was thinned to keep one locus per contig (O'Leary et al. 2018) to avoid possible linked loci due to proximity, resulting in 1771 SNPs as final dataset.This dataset was used to describe the relatedness degree among individuals, population genetic diversity parameters, population genetic structure and N e .

Genetic diversity and relatedness
The usual descriptors of genetic diversity parameters were estimated for each population (Table 1): allelic richness (A R ), observed and expected heterozygosities (H o and H e ) and inbreeding coefficient (F IS ).These parameters were estimated with the function basic.stats of the R package "hierfstat" 0.5-11 (Goudet 2005).We also reported the population average of the inbreeding coefficient estimates (F L ) calculated for each individual using the two likelihood methods implemented in COANCESTRY 1.0.1.10(Wang 2011).Both likelihood methods produced the same results.Relatedness coefficients and 95% confidence intervals were calculated per pair of individuals using both the dyadic (Milligan 2003) and triadic (Wang 2007) likelihood methods, assuming inbreeding and default parameters as implemented in COANCESTRY.Pairwise relatedness results were represented as heatmaps through the corrplot function of the R package "corrplot" 0.92 (Wei et al. 2021).Both methods produced similar results, so only triadic results are shown.

Population genetic structure
Population structure was inferred through three different methods.We used PGDSPIDER (Lischer and Excoffier 2012) to convert the VCF file to the different formats required by each approach.
We used the Bayesian clustering method carried out by the software STRU CTU RE 2.3.4 (Pritchard et al. 2000).We ran a preliminary analysis to infer the lambda value with K = 1, and once lambda was fixed, the analysis was configuring assuming admixture and correlated allele frequencies (Falush et al. 2003) between populations.The analysis was run five times per K value, each one including a burnin of 50,000 iterations and 100,000 iterations after the burnin.The first analysis included all populations and K values ranged from 1 to 5.Separately, we ran a specific analysis of La Mancha populations to improve the structure resolution in this region and K values ranged from 1 to 4. The STRU CTU RE plots were generated using CLUMPAK (Kopelman et al. 2015) and the number of clusters was explored by the rate of change in the log probability of data [Pr(X|K]) between successive K values (∆K) (Evanno et al. 2005).
As an alternative to the Bayesian clustering method, we used the multivariate analyses of Principal Components (PCA) (Hotelling 1933) and Discriminant Analysis of Principal Components (DAPC) (Jombart et al. 2010) in order to summarize the overall variability between individuals and to visualize more clearly the relationships between different clusters, respectively.PCA was performed with the function glPCA implemented in the R package "adegenet" (Jombart 2008).For DAPC analysis (retaining 50 principal components and 1 discriminant function) we used the find.clusters and dapc functions in the R package "adegenet" to identify the number of clusters or genetic groups (K's) under a Bayesian Criterion Inference (BIC) and to describe the relation between and within these clusters (Jombart et al. 2010), respectively.The first DAPC analysis included all populations.An additional analysis included only La Mancha populations to improve the genetic structure resolution of this group of lakes.
Differentiation between populations was estimated by pairwise F ST (Weir and Cockerham 1984).Point estimates and their associated p-values were calculated using stamppFst function in the R package "StAMPP" 1.6.3(Pembleton et al. 2013).

Effective population size estimation
Contemporary N e was estimated for each genetic unit using a bias-corrected version of the linkage-disequilibrium (LD) single-sample method by (Waples 2006;Waples and Do 2008) implemented in NeESTIMATOR 2.01 (Do et al. 2014).Although missing data in our final dataset is not high (less than 15% per site), NeEstimator v2 accounts adequately for the lack of missing data (Peel et al. 2013).Due to the close proximity between Longar and Altillo (0.7 km apart) and the total absence of genetic differentiation among them, they were considered as a single population to avoid the effect of low population sizes in estimations and obtain the "metapopulation" N e (Waples and Do 2010).We used PGD-SPIDER 2.1.1.5(Lischer and Excoffier 2012) to convert the VCF file into the "genepop" file required for the software.The analysis was set to allow for random mating due to the polygyny mating behaviour of anostracans (Belk 1991) and to calculate confidence intervals (lower and upper bounds).
Due to the absence of a reference genome, we could not ensure the absence of physically linked loci in our dataset (LD not caused by genetic drift but by physical proximity within chromosomes), which would downwardly bias N e estimates (Waples et al. 2016).Downwardly biased N e (^N e ) was estimated with two different datasets.We used the unpruned dataset of 13,538 SNPs (no LD filter applied), and also the thinned dataset (1 SNP per contig) in which possible physically linked loci are minimized to a larger extent to reduce bias.To estimate N e (unbiased), we took into account the number of chromosomes as suggested by Waples et al. (2016): ^Ne /N e = 0.098 + 0.219 * ln (n° chromosomes).Number of chromosomes was set to 24 due to B. media belonging to Chirocephalidae (Kořínková and Godyn 2011).Parametric and non-parametric (jackknife-Waples and Do 2008) methods were used to calculate confidence intervals.

Genetic diversity and relatedness
Allelic richness was similar in all populations (ranging from 1.363 to 1.383) (Table 1).Expected (0.135-0.139) and observed (0.044-0.067) heterozygosities were lower in Longar, Altillo and Alcahozo and higher in Tirez and Piñol.Hardy-Weinberg deviations were obtained for all the populations, being expected heterozygosity higher than the observed one (Table 1).The inbreeding coefficients were fairly high, showing Piñol the lowest values (F IS = 0.379, F L = 0.441; Table 1).
Values of genetic relatedness between pairs of individuals within and between populations were obtained by the two likelihood methods implemented in COANCESTRY were similar.Global genetic relatedness values between pairs of individuals ranged from 0 to 0.25 (Fig. 3).Relationships in the whole sample included parent-offspring, full siblings, half siblings, avuncular, grandparent-grandchild, double first cousins, first cousins, in addition to unrelated individuals.However, no twins were present in the sample (relatedness value equal or higher than 0.5) (Wang 2011).Close relations were very frequent within Piñol with relatedness values between 0.125 and 0.25, meaning that most of pairs of this population showed half sibs, avuncular, grandparent-grandchild and double-first cousin relations (relatedness = 0.125-0.25),and only two pairs of full siblings (relatedness ≥ 0.25) (Fig. 3).These close relationships also appeared within La Mancha populations.About half of the pairs of individuals within each population show relatedness values between 0.06 and 0.25, meaning the existence of half sibs, avuncular, grandparent-grandchild and double-first cousins while the other half are unrelated.Pairs of individuals between

Genetic structure
STRU CTU RE revealed strong genetic differentiation between Piñol and La Mancha complex (Fig. 4).When all populations were included in the analysis, the optimal number of clusters was K = 2.A hierarchical reanalysis removing Piñol from the dataset revealed substantial substructure within La Mancha, i.e., between Tirez and the remaining sites, resulting in an optimal number of clusters of K = 2.Most of the genetic variability was described in the first two principal components of the PCA analysis.Population scores from the first and second principal components of the PCA were plotted on two axes (PC1 and PC2), which cumulatively explained 8.67% of the total genetic variation (PC1 = 7.20%, PC2 = 1.47%) (Fig. 5A).The first component separated two clusters represented by La Mancha complex and Piñol population (Fig. 5A).The second component showed a slightly separated cluster represented by Tirez population (Fig. 5A).Intra-population genetic variability was high in all populations (Fig. 5A).DAPC analysis identified the lowest BIC value at K = 2 (463.092) in the dataset including all populations and at K = 1 (327.129) in that including only La Mancha populations, so only the result of the first dataset is shown (Fig. 5B).DAPC grouped individuals into two well differentiated genetic clusters corresponding to La Mancha complex and Piñol population (Fig. 5B).
The highest significant F ST values were obtained between Piñol and La Mancha complex (0.120-138).Within La Mancha the only differentiated population is Tirez, as the p-values for pairwise comparisons between Longar, Altillo and Alcahozo were ≥ 0.18 (Table 2).

Effective population size
Both the thinned and unpruned datasets yielded downwardly biased estimations of contemporary N e (^N e ) due to an excess of LD (specially the unpruned one).Biased N e obtained with the thinned dataset ranged from 490.3 (Longar + Altillo) to 4405.7 (Alcahozo) (Table 3).The lower bound of parametric and non-parametric confidence intervals showed the same trend, but the higher bounds were infinite for Longar + Altillo, Tirez and Alcahozo.The unpruned dataset produced finite intervals for the four sites when calculated by parametric method, however, higher bounds were infinite with jackknife method.Results indicated a significantly lower ^Ne for Longar + Altillo (280.1), while the other three populations had approximately twice this size.We applied the correction suggested by Waples et al. (2016) obtaining a bias of 0.794.Corrected N e became larger in all populations after bias adjustment.In the two cases, Longar + Altillo has the lowest N e value, approximately half the size of the remaining populations.

Population genomics and biology of Anostraca
We found that the genetic diversity (allelic richness and heterozygosity) of Iberian B. media populations is generally low with respect to other non-Anostraca branchiopods for mitochondrial, nuclear or microsatellite markers (Cesari et al. 2007;Horn et al. 2014).However, studies of Branchiopoda or other phylogenetically close crustaceans, in which genetic diversity indices are estimated based on SNPs are still scarce (e.g., Muñoz et al. 2016;Schwentner et al. 2018).Accordingly, obtained inbreeding levels are higher (Rogers 2015;Rowe et al. 2007;Velonà et al. 2009), which may be the reason for the resulted deviations from HWE.These comparatively (with Branchiopoda) low levels of genetic variability within anostracan populations may be a consequence of partial hatching during life cycles (Wahlund effect) (Brendonck et al. 2000a, b;Vanschoenwinkel et al. 2010;Wahlund 1928), abortive hatching and/or habitat monopolization by highly successful individuals following genetic drift or founder effects (Rogers 2015).In the case of B. media we cannot attribute low genetic diversity to limited dispersal and reduced gene flow because the lowest estimates do not concur with the most isolated populations (Ortego et al. 2010).Anyway, the values observed so far in Anostraca are not particularly low if compared to those of many other passively dispersed organisms (Montero-Pau et al. 2017;Unal and Bucklin 2010).In fact, genetic diversity estimates in B. media concur with those obtained so far in other Anostraca species (Bohonak 1998;Muñoz et al. 2009).High levels of inbreeding, and thus, low levels of heterozygosity in Anostraca might not be rare (Davies et al. 1997).Generation overlapping (cohabitation with relatives) and partial eggs hatching might lead to a reduction of heterozygosity in populations.Another possibility could be that inbreeding estimates were biased by the existence of loci with non-Mendelian inheritance contained in our dataset.However, some filters such as thinning (linkage disequilibrium) and keeping bi-allelic SNPs (eliminating loci with multiple alleles) decrease their occurrence probability.
Level of genetic differentiation (as measured by F ST ) between our populations is intermediate if compared to previous analyses on Anostraca (Boileau et al. 1992;Brendonck et al. 2000a, b;Ketmaier et al. 2008Ketmaier et al. , 2012)).Nevertheless, the three approaches performed in this study concurred that the strongest genetic differentiation occurs between Piñol (part of the Ebro Basin population complex) and La Mancha population complex.The strong genetic structure and large geographic distance between Piñol and La Mancha point to two independent and historically isolated genetic clusters.We can rule out such a pattern to be an artefact caused by the higher relatedness found in Piñol, as sequence data from the cox1 mitochondrial marker show a 4% of genetic distance (p) between Piñol and La Mancha complex and 0% among La Mancha lakes (unpublished data).These estimates indicate that obtained differentiation between these two subunits has been caused by a significant isolation through time, and not by the high inbreeding of Piñol.Piñol lake is in the Ebro 1 3  Optimal number of clusters is identified by the lowest Bayesian Information Criterion (BIC).PCA eigenvalues represent the number of principal components retained in the data transformation step.

La Mancha Piñol
Basin while La Mancha complex is in the Tajo Basin.These two basins, although of ancient origins (Tertiary Period), are currently covered by Quaternary sediments (Benito-Calvo et al. 2009).The time since the origin of the lakes seem to be sufficiently long to separate their populations as genetic units.The distinct climatic and morphogenetic nature of these two areas (northern and central of Iberian Peninsula) could have been supported the differentiation of the units (Rogers 2014b).With respect to La Mancha complex, all the methods agree that populations of Longar, Altillo and Alcahozo conform a single subunit which also correspond to the same geomorphological structure.Conversely, the different genetic composition of Tirez revealed by STRU CTU RE (Fig. 4), the partially non-overlapped distribution of individuals of this population displayed in the PCA (Fig. 5A in red) and its intermediate F ST values compared to the other La Mancha populations, may indicate this lake as a different genetic subunit.This population is settled in an artificially cut off adjacent section of the lake, which together with its smaller size may influence the frequency of visits by the dispersal agents, likely large birds.Consequently, dispersal events may become rarer (Boronat-Chirivella 2003;Rogers 2014a).This could imply that the population inhabiting the contiguous unsampled natural lake may likely share the genetic structure displayed by the other studied lakes (Longar, Altillo and Alcahozo) but additional analyses are needed to test for this hypothesis.Either way, this artificial population does conform a slightly different genetic unit, and therefore must be considered for conservation management.Therefore, Iberian B. media populations conform a minimum of three genetic subunits corresponding to (1) Longar, Altillo and Alcahozo, (2) Tirez (artificial pond) and (3) Piñol.The existence of a certain level of kinship between Piñol and La Mancha populations could be indicating some incidental dispersal.These results show that geographical distance affects B. media dispersal capacity at large scale (~ 300 km), but not at small or intermediate ones (~ 50 km).Longar, Altillo and Alcahozo are practically identical in genetic composition, which implies that dispersal and consequently, gene flow among these populations are very high.
Connectivity among these separated lakes seem to be consequence of passive vectors, birds like flamingos or ducks (Amat et al. 2005;Figuerola et al. 2003;Green et al. 2005).
On the other hand, the apparent inconsistent pattern of high relatedness and low inbreeding estimates in Piñol (Fig. 3), may have several explanations.The irregular filling of this lake might account for this seemingly contradictory pattern.The marginal areas of Piñol do not fill every year, so some eggs in the dry area will not hatch (diapause).When the filling is complete and the diapause eggs hatch, there will be co-occurrence of different generations (Rogers 2014a, b).Therefore, it may be more likely to have sampled related individuals at Piñol after complete filling increased the degree of relatedness.In addition, inconsistent values of relatedness and inbreeding (F IS, F L ) are feasible to be found due to the different genetic data used for their estimation (allele frequencies for inbreeding-Wright 1951-and proportions of modes of identity by descent for relatedness -Wang 2011).Observed differences in inbreeding between Piñol and La Mancha populations (these last with a high inbreeding level) might be a consequence of the high concentration of additional salt lakes across surrounding area of Piñol (e.g., at least eleven known populations in a 25 km radius, Alcorlo and Baltanás 1999).These unstudied populations might be a constant source of new immigrants to Piñol contributing to decrease F IS values.
Contemporary N e values in B. media are high compared with estimates in other crustacean lineages, such as freshwater (Ada 2021) and marine decapods (Anagnostou and Schubart 2017;Cabezas et al. 2012;Heras et al. 2019), amphipods (Gergs et al. 2019) and isopods (Guzik et al. 2012).The highest N e estimates were obtained with the thinned dataset, which seems to be reasonable because of the excess of removed LD (Waples and Do 2010).Conversely, the non-LD pruned dataset, which contains the whole LD, physical linkage and that originated by genetic drift, produces lower N e estimates.In both cases, the excess of LD is balanced with Waples et al. (2016) correction, which increases N e estimates.Highest N e concur with the largest sized lakes (Tirez, Alcahozo and Piñol).The populations of Longar and Altillo are located in small to medium-size ponds near a large lake, so a lower N e is expected in relation to the possible smaller population size (Nunney and Elam 1994).However, these values are still high, which suggests that habitat size does not considerably influence the N e of Anostraca populations.Large population sizes have been previously registered in small habitats like rock pools or those from wheel tracks (Brendonck et al. 2000a, b;Timms 2006;Vanschoenwinkel et al. 2013).Comparing with standardized N e thresholds, a minimum value of 50-70 is needed to minimize inbreeding depression, and 500 to maintain an adequate evolutionary potential (Franklin 1980;Caballero et al. 2017).Accordingly, the large observed N e is buffering the effect of relatedness and genetic drift on genetic diversity levels (Giesel 1971;Keller and Waller 2002).In addition, high N e estimates suggest a balanced sex ratio and reproductive success within all populations (Caballero 1994;Frankham 1995;Nunney and Elam 1994).These results proved invaluable for a first assessment of the conservation status of B. media, as well as to infer that this species likely maintains large N e to compensate limited gene flow and subsequent inbreeding effects and maintain adequate levels of genetic diversity.However, we must highlight that the estimation of N e in anostracans following LD method infringes the assumption of non-overlapping generations that can be underestimating the results (Waples et al. 2014), which is confirmed with the high level of relatedness obtained in all the study sites.This means that our estimates of contemporary N e should be taken with caution.Future studies should confirm our estimates by using the temporal approach, proved to be more robust to scenarios with overlapping generations (reviewed by Wang et al. 2016).In our case, this issue is not so worrying since despite having overlapping generations, our N e estimates are still high, giving an optimistic approximation about the conservation status of Iberian Branchinectella media.
Population extirpation in metapopulational systems may come from natural causes such as high levels of inbreeding or limited reproduction effectivity, which lead to the decrease of genetic diversity.The extent of the effect of these factors within populations depends on population size.Small populations usually tend to disappear randomly and faster than larger populations (Vrijenhoek 1994).Most species of Anostraca frequently occupy small, human-modified ponds (Horváth and Vad 2015;Sainz-Escudero et al. 2021), usually forming metapopulation systems, and therefore likely subjected to high population turnover.However, the case of Iberian B. media, in which habitat restrictions (saline ponds) determine extremely limited availability of possible population settling, population turnover or recolonization is highly unlikely.This fact together with the limited connectivity (low levels of genetic diversity) between the three studied genetic units would make short term local extinction a probable output for B. media populations.Only large population sizes can assure long-term population survival under such circumstances, as it is precisely the case observed in our sampling.According to our data, B. media diversified within Iberia at least from the late Pliocene (unpublished data), and therefore it was established in the territory at least since then.In addition, an historical record of B. media in Piñol (Brehm and Margalef 1948) also assures the continuity of the population from, at least, fifty years ago.Following these considerations, we put forward that other species of Anostraca whose populations are likely subjected to genetic isolation, such as those of Branchipodopsis in southern Africa (Brendonck et al. 2000a, b), may be relying on their large N e to persist through time despite genetic diversity losses due to stochasticity.Indeed, there are other invertebrates adapted to temporary aquatic ecosystems which also display cases of strong genetic differentiation such as Cladocera, Copepoda and Ostracoda (Haag et al. 2006;Rossi et al. 1994;Scheihing et al. 2011).In order to recognize threatened units of species living in temporary aquatic ecosystems and with limited dispersal capacity, conservation studies should include the estimation of effective sizes, accompanied, if possible, by an approximate observation of the total size of the population as an additional support (Wang & Rogers 2018).

Conservation implications
At present, Iberian lineages of Branchinectella media would not require specific nor urgent conservation actions.Large effective population sizes in B. media populations, together with current local strict protection of their inhabited and potential salt lakes guarantee its short-term persistence.However, our optimistic N e estimates should not preclude this species from being monitored in the future.Continental aquatic ecosystems are currently threatened by direct and indirect anthropogenic disturbances (Dudgeon et al. 2006;Reid et al. 2019).Among them, temporary ecosystems are especially vulnerable (Zacharias and Zamparas 2010).Climatic change and the resulting unpredictable rainfall distribution threaten the fill of endorheic basins, halting the formation of the necessary aquatic habitats for development of anostracans (Asem et al. 2019;Rautio et al. 2011) and other invertebrates (Williams 1997).Similarly, changes in light and temperature alter hatching success in Anostraca (Tladi et al. 2020).Besides that, agricultural and farming development causes aquifer collapse due to water overuse (Löffler 1993) and water pollution deriving of the use of pesticides, chemical fertilizers (Migliore et al. 1993;Sánchez-Bayo 2006) and veterinary pharmaceuticals (Udebuani et al. 2021).In addition, modifications of the temporary nature of lakes with artificial water inputs (Brock et al. 1999) impede fairy shrimps and other organisms adapted to habitat filling-drying conditions to complete their life cycle.The persistence of temporary habitats is thus compromised by anthropogenic disturbances, which may condition the survival of B. media despite the evolutionary efforts (high N e ) nature is making for its persistence.

Fig. 1 Fig. 2
Fig. 1 Branchinectella media general habitus and typical habitat of the species.A Adult male from Salada de Piñol (Zaragoza, Spain) photographed ex situ.Note the large size of male eyes, a characteristic trait of Branchinectella.B Adult female from the same locality.C Two young females swimming in shallow pools of Salada de Piñol (Zaragoza, Spain) photographed in situ.D Laguna de Alcahozo (Ciu-

Fig. 3
Fig.3Heatmap of pairwise relatedness between specimens of B. media as revealed by the triadic likelihood method implemented in COANCESTRY.Individuals appear in the same order in the horizon-

Fig. 4
Fig. 4 Genetic structure inferred by STRU CTU RE for 103 individuals of B. media based on 1771 SNPs (thinned dataset).Each vertical bar represents and individual coloured according to that individual

Fig. 5
Fig. 5 Genetic structure represented by multivariate analyses of A PCA and B DAPC for all B. media populations.Principal Component Analysis represent the distribution of the sample genetic variation in the two first axes.Discriminant Analysis of Principal Components

Table 1
Genetic diversity parameters for each of the analysed populations of B. media n number of individuals, A r allelic richness, H o observed heterozygosity, H e expected heterozygosity, F IS inbreeding coefficient, F L averaged likelihood-estimated inbreeding

Table 2
Pairwise genetic differentiation between B. media populations calculated by F ST parameter and corresponding p-values above and below de diagonal, respectively

Table 3
Contemporary biased effective population size (^N e ) estimates produced by N E ESTIMATOR (LD approach) and bias-corrected (N e ) following Waples et al. (2016) adjustment Thinned dataset contains 1 SNP/contig (1771 SNPs final dataset) to reduce LD because of proximity between loci in chromosomes.Dataset without LD pruning contain an excess of LD