Brine shrimps adrift: historical species turnover in Western Mediterranean Artemia (Anostraca)

Brine shrimps (Artemia) have undergone geographic range and demographic expansions as a result of their interaction with humans since the beginning of salt harvesting. This interaction has favoured the expansion of some species but compromising the survival of others. Mediterranean native populations of Artemia salina from coastal salterns and lagoons are facing the presence and expansion of the introduced and invasive American species Artemia monica (= A. franciscana). However, this species could not be the only threat. Parthenogenetic populations of the Asian species A. urmiana and A. sinica are widespread along the Mediterranean and other areas of the world. In this work, with the use of large cox1 and mitogenomic datasets, phylogenetic and phylogeographic inferences, and a time calibrated tree, we confirmed the Asian origin and recent arrival of the current Western Mediterranean parthenogenetic populations of Artemia. In addition, the replacement of Iberian populations of A. salina by Asiatic parthenogenetic populations lead us to recognize parthenogens as invasive. Current salterns development and commercial importance of Artemia make human-mediated introduction probable. These results demonstrate again the impact that changing human interests have on population expansion or decline of species adapted to anthropogenic habitats. Artemia salina decline makes urgent the implementation of conservation measures such as its use in fish farming and salt production or its inoculation in inland salterns.


Introduction
Wild species able to colonize man-modified habitats are often capable to expand vastly their geographic ranges. The newly created habitats or transport provided by humans and lack of competitors may result in demographic bursts that significantly contribute to geographic expansions of their ranges (Senar et al. 2019;Tollenaere et al. 2010;Lewis et al. 2019). However, human needs and production systems change rapidly through time, and what was a favourable system for a certain species might change drastically, either favouring a different species, or limiting dramatically their available habitat (Bøhn et al. 2008;Dafni et al. 2010). As a consequence of these changes, populations of the once widespread, highly successful species, might suffer dramatic declines or become extirpated, and if the original natural habitat would be no longer available, the entire species might face extinction.
One of these taxa whose interactions with humans have produced local demographic bursts and vast geographic range expansions are brine shrimps, small crustaceans of the genus Artemia Leach, 1819 (Branchiopoda: Anostraca). The currently recognized five species of Artemia (Sainz-Escudero et al. 2021) are inhabitant of worldwide continental aquatic saline ecosystems (Van Stappen 2002). Their original habitat, inland saline lakes and coastal lagoons are limited, and their original distribution was probably restricted to certain geographic areas. This is still the case for the American Southern Cone species, Artemia persimilis Piccinelli andProsdocimi, 1968 (Amat et al. 1994). However, Neolithic originated salt harvesting for human needs (Weller and Dumitroaia 2005;Fíguls et al. 2007;Manrique 2011) was based for centuries in the creation of salterns, that accumulate salty water from wells or sources associated to subterranean diapirs or sea water, providing a completely new and suitable habitat for brine shrimps to be colonized (Martínez-Abraín and Jiménez 2015). The species of Artemia expanded their geographic ranges favoured by passive dispersal of their resistance eggs mediated by migratory birds (Green et al. 2005;Muñoz et al. 2013) or by anthropogenic movements motivated by salterns development or commercialization for fish farming (Vikas et al. 2012;Sorgeloos et al. 2001). However, human commercial interests have also come along with disturbances in brine shrimp species distribution, even compromising the survival of some. The large demand of Artemia for fish farming triggered the introduction and invasion of the North American cultivated species Artemia monica Verrill, 1869 (= A. franciscana Kellogg, 1906) into coastal salterns all over the world (Triantaphyllidis et al. 1994;Amat et al. 2007;Mura et al. 2006;Ruebhart et al. 2008;Scalone and Rabet 2013;Saji et al. 2019). This species is displacing the Mediterranean native species Artemia salina (Linnaeus, 1758) (Oscoz et al. 2010;Horváth et al. 2018) possibly due to its high adaptive potential and physiological plasticity that enhance its invasion range capacity (Dlugosch and Parker 2008).
Contrary to this idea, we hypothesize that the historically expanded populations of the Mediterranean species A. salina (Muñoz et al. 2014) are being affected not only by the human-induced introduction and invasion of A. monica but also by those of the Asian parthenogenetic populations of A. urmiana and A. sinica, which likely would have arrived recently. These invasions, probably associated to salt farming interests, seem to be displacing the native populations of A. salina, which might lead the species at the brink of extinction in some areas. However, there is very little previous evidence about the introduced and invasive character of parthenogens in the Mediterranean region. To test for these hypotheses, we need to check the historical and current presence of A. salina in Western Mediterranean salterns, confirm the recent arrival of the Asian parthenogenetic populations, and discuss the potential damage they are causing on A.
salina. To do this, (1) we surveyed Iberian and some other Western Mediterranean inland and coastal salterns to document the current presence of Artemia and confirmed their identity by studying their morphology and sequencing a fragment of the mitochondrial cytochrome c oxidase I gene (cox1). Historical data on the presence of Artemia in Iberian salterns contained in previous bibliography was used to discuss the invasive character of parthenogens. (2) We identified the approximate geographical origin of the parthenogenetic A. urmiana and A. sinica populations in the Western Mediterranean region by using our cox1 data, together with all the previously available cox1 dataset of Artemia. Finally, (3) we used next generation sequencing (NGS) to generate the first complete mitogenomes for parthenogentic individuals of A. urmiana and A. sinica, and phylogenetically compare them with already published bisexual Artemia mitogenomes. We inferred divergence times estimation on our mitogenomic tree in order to find out the approximate time of appearance of parthenogenesis across the evolutionary history of the genus, and estimate their possible timing of expansion into the Western Mediterranean region.
In this work, we use Western Mediterranean brine shrimps (mainly from Iberian Peninsula) as a model to exemplify the human dependent fate to which human-favoured species are subjected following changes during the economic and productive procedures.

Sampling and sequencing
We collected around five to ten adult individuals of Artemia from artificial inland and coastal salterns from a total of 30 localities, of which 25 correspond to Iberian Peninsula and 5 to other Western Mediterranean and North Africa locations. Information about sampling localities and their geographical coordinates is included in Table 1. Specimens were captured with an aquarium net, photographed ( Fig. 1), georeferenced, and preserved in absolute ethanol and stored at − 20 °C at the Museo Nacional de Ciencias Naturales (MNCN-CSIC) from Madrid (Spain).
Total genomic DNA (samples deposited at the MNCN DNA-tissues Collection) was isolated according to the protocols described by Hwang et al. (2019). A partial fragment of the cytochrome c oxidase subunit I (cox1) was amplified via the polymerase chain reaction (PCR) employing the universal primers LCO1490 and HCO2198 (Folmer et al. 1994. PCR reactions were performed in a total volume of 25 µL that contained 2µL of extracted DNA, 1 µL of dNTP (10 mM), 2.5µL of reaction buffer 10x (Tris-HCl, pH 8.3, Biotools), 0.8 of µL MgCl 2 (50 mM), 0.5 µL of each primer (10 µM), 17.3 µL of distilled water and 0.4 µL of Taq DNA polymerase (Biotools, 5 ud/ µL). PCR reactions consisted of 1 cycle of 5 min at 95 °C for initial denaturation, 40 cycles of 45 s at 42 °C and 1 min at 72 °C, and a final extension of 10 min at 72 °C. Amplified PCR products were visualized by electrophoresis in a 1.5% agarose gel. Samples   Additionally, to our data, all the available Artemia cox1 sequences available in GenBank (Valsala et al. 2005;Hou et al. 2006;Tizol-Correa et al. 2009;Muñoz et al. 2008Muñoz et al. , 2010Muñoz et al. , 2013Maniatsi et al. 2009Maniatsi et al. , 2011Maccari et al. 2013b;Eimanifar and Wink 2013;Eimanifar et al. 2014Eimanifar et al. , 2015Eimanifar et al. , 2016Asem et al. 2016Asem et al. , 2019Asem et al. , 2020Naganawa and Mura 2017;Horváth et al. 2018) and one of Branchinecta ferox used as outgroup (LT821334 [Rodríguez-Flores et al. 2017]) were retrieved in order to build a dataset represented by 1505 sequences, that allowed us to depict the structuring of the genus through the Neighbour Joining analysis and to perform phylogeographic analyses. Some dissimilar sequences that featured stop codons when traduced to amynoacids were removed from the analyses due to the existence of pseudogenes according to Rode et al. (2021). In parallel, four parthenogenetic individuals (one corresponding to A. urmiana lineage and three to A. sinica) from four localities in the Iberian Peninsula were chosen to be high-throughput sequenced (at AllGenetics, A Coruña. Spain; Table 1). The DNA extraction protocol, library preparation, and sequencing process, were the same as those explained in Sainz-Escudero et al. (2021). Genome assembly of the parthenogenetic individuals of A. urmiana and A. sinica was carried out using the cox1 sequence as reference. Finally, annotation was performed using MITOS2 (Bernt et al. 2013), checking manually the start and stop codons of all coding genes. Mitogenomes of parthenogenetic Artemia urmiana and A. sinica were deposited in GenBank. GenBank accession numbers of the mitogenomes are indicated in Table 1.

Mitochondrial (cox1) analyses and phylogeography
In order to visualize the general structure inside the genus Artemia, the cox1 dataset and a sequence of   and Swofford 2003) ( Fig. 2A). A bootstrap (BS) analysis (1000 replicates) was used to assess node support (implemented in PAUP*). Part of this dataset was also used to perform phylogeographic analyses. We made two cox1 haplotype networks with sequences related to A. urmiana and A. sinica clades (Fig. 3) in order to represent the geographical distribution of the allele diversity and find the position of the diploid and tetraploid parthenogenetic populations from Iberian salterns inside these networks. We first used DNA Sequences Polymorphism 6.12.01 (Rozas et al. 2017) to generate the collapsed matrix of unique alleles, and then, networks were constructed using the TCS algorithm applied through Population Analysis with Reticulate Trees (PopART) 1.7 software (Leigh and Bryant 2015) to shape the relationships between the population individuals. Sequences coming from GenBank which had more than the 50% of missing data were removed from the collapsed matrix. Information about the sequence-haplotype correspondence and their bibliographic source is included in Tables 2 and 3 (A. urmiana and A. sinica, respectively).

Mitogenomic phylogeny and divergence time estimates
We reconstructed a new phylogenetic hypothesis based on complete mitochondrial genomes, to sort out the main Artemia lineages, and to identify a possible time of appearance of the parthenogenetic populations (Figs. 2,4). The data matrix of Artemia mitochondrial genomes was composed by those available at the GenBank database (see Sainz-Escudero et al. 2021) and the new four genomes from the Asian parthenogenetic individuals of A. urmiana and A. sinica. The complete matrix includes 16 terminal taxa.
Phylogenetic reconstruction was performed using a Bayesian Inference approach implemented in MrBayes version 3.2.6 (Ronquist et al. 2012) and divergence time estimation was carried out in BEAST 1.7 (Drummond et al. 2012). These analyses were performed following the protocol and considerations concerning molecular clock calibration and priors are specified in Sainz-Escudero et al. (2021).

Cytochrome c oxidase I allele diversity and phylogeography
Although the Neighbour Joining analysis was not aimed to clarify the phylogenetic relations between Artemia groups, it helped us to visualize sequence similarities. Accordingly, ten more or less differentiated mitochondrial groupings exist all over the world ( Fig. 2A). This exploratory method supports the monophyly of all recognized species (BS values over 95%). The bisexual species A. monica shows a large genetic diversity with many geographical or ecologically isolated populations (Browne and Bowen 1991;Muñoz et al. 2013). Artemia urmiana is represented by three main clades constituted by bisexual populations from Western Asia, Tibet, Kazakhstan, and diploid/ triploid parthenogenetic populations from Western Asia, Africa, Australia, Europe, Madagascar and Russia. Artemia sinica is formed by two reciprocally Fig. 3 Fig. 2A).
The phylogeographic analyses of the mitochondrial cox1 for the Eastern Asian Clade (A. sinica) resulted in a total of 29 haplotypes (Fig. 3A, Table 2).

Phylogeny of Artemia
The topology of the mitogenomic Bayesian analysis was totally congruent with the topology of the ultrametric tree obtained with BEAST (Fig. 2B). All nodes are supported with a posterior probability of 1 (PP = 1). Phylogenetic relationships within the main Artemia lineages concur with those previously proposed by Sainz-Escudero et al. (2021). Parthenogenetic samples were recovered in two non-sister clades. One parthenogenetic sample from Murcia (Spain) with unknown ploidy (mtArt_1), falls within the A. urmiana clade closely related with the sample of Urmia Lake, Iran (Fig. 2B). The samples of unknown ploidy from Santarem (Portugal) (mtArt_2), Guadalajara (Spain) (mtArt_3) and Zaragoza (Spain) (mtArt_4) conform a monophyletic group closely related with A. sinica (Fig. 2B). The divergence time estimates for the main clades concur with that previously proposed by Sainz-Escudero et al. (2021) (Fig. 4). The time to the most recent common ancestor of the clades of bisexual A.
urmiana and its parthenogentic variants is placed during the Pleistocene (Mean 1.25 Ma, . Separation between the A. sinica lineage and its parthenogenetic variants occurred about  more than a half million years later (Mean 0.68 Ma, 95% HPD 1.59-0.95 Ma) (Fig. 4).

Discussion
Our study supports the multiple independent origin of parthenogenesis in Artemia (Muñoz et al. 2010;Maniatsi et al. 2011;Asem et al. 2016) and the geologically recent origin of parthenogenetic populations within A. urmiana and A. sinica lineages (Eimanifar et al. 2015;Sainz-Escudero et al. 2021;Rode et al. 2021). In the Iberian Peninsula, parthenogenetic populations were noted since the end of the twentieth century (Amat 1979(Amat , 1980Triantaphyllidis et al. 1998). Later, a few populations from coastal and inland salterns were identified as part of the Asian lineages because of their phylogenetic relation with bisexual A. urmiana and A. sinica (Muñoz et al. 2010;Maniatsi et al. 2011;Maccari et al. 2013a). Our phylogeographic analyses corroborate the Asian origin of Iberian parthenogenetic samples. Cox1 haplotypes shown by parthenogenetic Iberian populations are shared with those of parthenogens from Asian localities, or are differentiated from them by no more than three point mutations, indicating shallow and recent divergences. Parthenogenetic populations of the A. urmiana lineage inhabit Iberian inland and coastal salterns, although they are more frequent in coastal ones (Amat 1980;Amat et al. 1995;Muñoz et al. 2010;Maccari et al. 2013a;this study). By contrast, parthenogenetic populations of the A. sinica lineage are mainly found in inland salterns (Amat 1980, this study) with only two coastal locations previously registered (Maniatsi et al. 2011;Maccari et al. 2013a).
Causes for the arrival of parthenogens to the Iberian Peninsula might be diverse. Bird-mediated dispersal seems likely for the colonization of large coastal salterns currently dominated by A. urmiana parthenogens (Persoone and Sorgeloos 1980;Green et al. 2005;Sánchez et al. 2012;Muñoz et al. 2013Muñoz et al. , 2014, but not for inland ones, predominantly occupied by A. sinica parthenogens, because of their smaller size, and generally less suitable conditions for bird feeding or nesting. The rise of salterns development makes anthropogenic transport also likely Vol.: (0123456789) (Léger et al. 1986;Sorgeloos et al. 2001;Dhont and Sorgeloos 2002;Van Stappen et al. 2007;Muñoz et al. 2008). The deliberate inoculation of Artemia by man in salterns have been a common practice because of the benefits that they generate for sea-salt production (Davis 1974;Persoone and Sorgeloos 1980;Dhont and Sorgeloos 2002), as it has been already reported in Australia (McMaster et al. 2007). In the Iberian case, the presence of A. urmiana in only a few inland salterns could be explained by current anthropogenic dispersal from coastal salterns to these precise inland salterns, most of them recently remodelled and conditioned for salt production (e.g. Oro and Imón salterns, in Navarra and Guadalajara respectively) or where the species has been possibly introduced for aquaristic purposes (La Seca, Madrid).
But, is the presence of parthenogenetic strains disturbing the stability of Mediterranean native A. salina populations? Historical records indicate that the native A. salina was present in many Iberian coastal and inland salterns (Amat 1980;Amat et al. 1995;Triantaphyllidis et al. 1998). Previous data also indicate that A. salina and Asian parthenogens were able to coexist in some places (Amat 1980;Amat et al. 1995;Barata et al. 1995;Van Stappen 2002). However, two salterns sampled in this study, Poza de la Sal (Burgos) (Fig. 5C) (Maccari et al. 2013a, b; this study) and probably will have a similar fate. The decline of A. salina likely reflects the higher colonization success of polyploid parthenogens (Browne and MacDonald 1982;Browne et al. 1988;Zhang and King 1993;but see Browne and Wanigasekera 2000). In fact, the colonization ability of parthenogens may even overcome that of the invasive American species in continental environments, because so far, only one population of this species has been found in Iberian inland salterns (Gerri de la Sal-Lleida, Amat et al. 2007;Muñoz et al. 2014 sub A. franciascana). These data provide an evidence of the displacement effect and competitive potential of parthenogens against A. salina populations and thus, confirming their invasive character.
However, the decline of A. salina in Iberia (Muñoz et al. 2008; this study) is not only a consequence of the parthenogens invasion, but also of a generalized abandonment and decay of inland salterns (Amat et al. 2007). Many inland salterns where the presence of Artemia was previously reported (Triantaphyllidis et al. 1998;Muñoz et al. 2008) are no longer in use and their tanks and ponds abandoned and dried-out (Armallá and Rienda in Guadalajara, Arcos de las Salinas in Teruel, Peralta de la Sal in Huesca, and a long etc.) (Fig. 5A, B). New surveys are still needed in order to check for the continuity of the remaining populations in the near future (Fig. 5D). The evident decline of A. salina, together with the relative facility to reverse the situation in inland habitats, makes conservation actions urgent.

Conclusions
Our phylogenetic and phylogeographic analyses support a recent origin of parthenogenesis within the lineages of A. sinica and A. urmiana and also the recent colonization of the Western Mediterranean by Asian parthenogenetic populations. In Iberia, Asian parthenogens have replaced previously established populations of A. salina, which together with a dramatic loss of habitat due to the abandonment and decay of inland salterns, is driving A. salina to a critical situation, worsened by the invasion of coastal salt lakes by the American species, A. monica (Amat et al. 1995(Amat et al. , 2007Muñoz et al. 2008). In this situation, feasible conservation measures are necessary. Development of A. salina hatcheries would be desirable to avoid the import and spread of the American and Asian species. The restoration of Iberian historical inland salterns (Neolithic, Roman, Medieval, or pre-Industrial Revolution) and the inoculation of A. salina native cysts would be useful and effective steps to maintain an inland network of populations less susceptible to further replacement.
The case of Artemia invasions shows that the distribution range and population fate of the species of Artemia currently depends on human economical activities, not only because of the strong commercial displacement to which they are subjected, but because of their almost complete affiliation and success in artificial salterns (Carscadden et al. 2020). If human activities have been responsible to drive a species from centuries of historical population success to a current dramatic population decline, as it is the case of Mediterranean A. salina, we should be also responsible to revert its current pathway to extinction.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This research was funded by the project-grant IND2018/AMB9692 ("Doctorado Industrial" program of Comunidad Autónoma de Madrid) to MG-P. LS-E is supported by the "Doctorado Industrial" grant (IND2018/AMB9692), through the Fundación Global Nature, with support of MITECO (Ministerio para la Transición Ecológica y el Reto Demográfico, Spain). EKL-E is supported by a doctoral scholarship from CONACyT-Mexico (330519/472100). PCR-F is supported by a E.O Wilson postdoctoral fellowship from the MCZ (Museum of Comparative Zoology) at Harvard University.

Availability of data and material
The datasets generated and analysed during the current study are available for their visualization and download at the GenBank repository. Gen-Bank accession numbers are included in the manuscript. DNA samples have been deposited at the Crustacean Collection of the MNCN (Museo Nacional de Ciencias Naturales) in Madrid, Spain. MNCN codes are included in the manuscript.

Declarations
Conflict of interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Consent for publication
All authors agreed to publish the manuscript.

Consent to participate
All authors agreed to participate in the study and its publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included