Mitogenomic phylogeny and divergence time estimation of Artemia Leach, 1819 (Branchiopoda: Anostraca) with emphasis on parthenogenetic lineages

Abstract

The brine shrimp Artemia, a crustacean adapted to the extreme conditions of hypersaline environments, comprises nine regionally distributed sexual species scattered (island-like) over heterogeneous environments and asexual (parthenogenetic) lineages with different ploidies. Such sexual and asexual interaction within the genus raises questions regarding the origin and time of divergence of both sexual species and asexual lineages, including the persistence of the latter over time, a problem not yet clarified using single mitochondrial and nuclear markers. Based on the complete mitochondrial genome of all species and parthenogenetic lineages, this article first describes the mitogenomic characteristics (nucleotide compositions, genome mapping, codon usage, and tRNA secondary structure) of sexual species and asexual types and, secondly, it provides a comprehensive updated phylogenetic analysis. Molecular dating and geographical evidence suggest that the ancestral Artemia taxon originated in ca. 33.97 Mya during the Paleogene Period. The mitogenomic comparisons suggest that the common ancestor of diploid and triploid parthenogenetic lineages (ca. 0.07 Mya) originated from a historical ancestor (ca. 0.61 Mya) in the Late Pleistocene. Additionally, the common ancestor of tetraploid and pentaploid parthenogenetic lineages (ca. 0.05 Mya) diverged from a historical maternal ancestor with A. sinica (ca. 0.96 Mya) in the early Pleistocene. The parthenogenetic lineages do not share a direct ancestor with any sexual species. The Asian clade ancestor diverged more recently (ca. 14.27 Mya, Middle Miocene). The mitogenomic characteristics, maternal phylogenetic tree, and especially divergence time prove that A. monica and A. franciscana are two biological species.

Introduction

Artemia Leach, 1819 is a genus of halophilic zooplanktonic crustaceans commonly known as “brine shrimp”, which are distributed in isolated hypersaline habitats across the world [1] with an island biogeography dispersal model [2,3,4,5]. Artemia is highly plastic across morphometric and morphological traits in both natural and laboratory conditions [6,7,8,9,10,11] and tends to exhibit accelerated molecular evolution very likely due to the mutagenic effect associated with ionic strength variance [12] and other stringent ecological conditions [13], combined with natural population expansion and contraction cycles. Additionally, natural salt lake ionic compositions, temperatures, and salinity heterogeneity, further affected by climatic changes and human perturbations [14, 15], have left mitogenomic signatures in Artemia bisexual species [1], which are also reflected in private mitochondrial haplotypes [4, 5, 16,17,18,19,20,21,22,23]. There is no classic “identification key” for the genus Artemia and taxonomic delimitation of sexual species and parthenogenetic lineages has been limited to reproductive mode [4, 24], geographical distribution [25, 26], and molecular markers [4, 5, 18, 21, 27].

Artemia consists of nine sexual species and four obligate parthenogenetic lineages [4, 23], regionally distributed around the world (except Antarctica) in arid and semiarid areas. Three species in the New World: A. monica Verrill, 1869 (Mono Lake, USA), A. franciscana Kellogg, 1906, (North, Central and South America, but introduced for commercial aquaculture in Eurasia, Africa and Australia), and Artemia persimilis Piccinelli and Prosdocimi, 1968, (Chile and Argentina). Six species in the Old World: A. salina (Linnaeus, 1758) (Europe and Africa); A. urmiana Günther, 1899 (Lake Urmia, Iran and the Crimean Peninsula); A. sinica Cai, 1989 (China); A. amati Asem, Eimanifar, Hontoria, Rogers and Gajardo, 2023 (Kazakhstan); A. tibetiana Abatzopoulos, Zhang and Sorgeloos 1998; and A. sorgeloosi Asem, Eimanifar, Hontoria, Rogers and Gajardo, 2023, the last two from the Tibetan Plateau. Obligate parthenogenetic Artemia consists of four lineages with different ploidy levels (di-, tri-, tetra- and pentaploid), all occurring in the Old World [18, 19, 21, 25, 28,29,30,31,32,33,34,35,36,37,38] and Australia [39,40,41].

Artemia evolutionary relationships have been studied with different tools, with some concordances and differences, from electrophoresis-based enzyme assays [42], morphometrical characteristics [4, 25, 26, 30, 32, 43,44,45], sequences of partial mitochondrial genes (COI, 16S and 12S) [4, 5, 17,18,19, 21, 27, 34, 38, 45], as well as the complete sequence of the nuclear ITS1 gene [4, 21, 38, 46, 47], to nuclear simple sequence repeats (SSR) markers [4, 48, 49].

Triantaphyllidis et al. [43] argued that A. urmiana was significantly isolated from the other Asian groups based on morphometric characters. However, the use of mitochondrial markers demonstrates that Asian species have a common ancestor [4, 18, 38]. Baxevanis et al. [45] found that the Asian A. tibetiana and A. urmiana are genetically close taxa while morphometrically dissimilar. In contrast, South American A. persimilis and A. urmiana are significantly different genetically, while having similar morphometric patterns [45, 50]. Asem et al. [4] demonstrated no conformity between mitochondrial marker results and morphometric patterns. Examples of such inconsistencies follow.

The Mediterranean A. salina and A. franciscana group together based on mitochondrial 16S rDNA RFLP analyses [45]. In the mitochondrial COI-based study of Muñoz et al., A. salina and A. persimilis were placed in the same clade, and A. sinica as a sister clade of A. franciscana and other Asian Artemia, with A. franciscana in between [16]. Maniatsi et al. using nuclear ITS1 and mitochondrial COI, reported the American A. franciscana was sister to the Asian species and parthenogenetic lineages, and that A. persimilis was basal to both [18]. Using COI sequences, Eimanifar et al. found A. salina to be basal [21].

Overall, these studies show that Artemia maternal phylogenetic relationships are still unclear. Moreover, recent studies show that phylogenetic reconstructions based on short mitochondrial marker sequences do not necessarily reproduce the tree topologies and divergence estimates found using the complete mitogenome, which has increasingly become the marker of choice [51, 52].

Here, we use the complete mitochondrial genome (henceforth mitogenome) to refine Artemia phylogeny and divergence times and resolve the different results produced by the various studies on Artemia species and parthenogenetic lineages.

The mitogenome represents specific maternal origins and involves rapid evolutionary alterations without recombination [5, 53,54,55,56,57,58,59]. Additionally, mitochondria are responsible for ATP production, which is necessary to respond to critical environmental conditions such as high salinity and related hypoxia, demanding major oxygen consumption. The mitochondria regulate the expression of multiple genes in response to those and other environmental conditions [60]. The coordinated functioning of the whole mitochondrial complement is required for local adaptation. We used the Artemia mitogenome to (1) clarify Artemia phylogeny, (2) reconsider parthenogenetic lineage origins hypotheses, and (3) estimate divergence times and evolutionary ages of Artemia members. We also studied nucleotide compositions, genome mapping, codon usage, and tRNA secondary structure. A comprehensive updated phylogenetic study considering all species (some recently described) and parthenogenetic lineages is necessary to comprehend Artemia's evolutionary history, origin, and divergence.

Materials and methods

Sampling

Nine Artemia species and four parthenogenetic lineages were studied from topotype material, except A. salina, which has disappeared from its type locality (Lymington, England), and A. persimilis (Salinas Grandes de Hidalgo, Argentina), which was unavailable. We did not have access to the samples of “diploid parthenogenetic Artemia with A. urmiana type mitochondria (see Discussion). Further information of species and parthenogenetic lineages is summarized in Table 1.

Table 1 List of studied Artemia species and parthenogenetic lineages, their provenance and egg bank accession

Eggs from each locality were cultured following Hontoria and Amat [25], except A. monica, which was unhatchable; A. monica was sequenced directly from eggs (for more information, see DNA Extraction and Sequencing section). Only females were sequenced after reproduction mode confirmation by individual culture [4, 61], as in most Old-World populations, Artemia species and parthenogenetic lineages may coexist, and diploid parthenogens may produce rare males (see [62]). Ploidy levels were determined by karyotyping using cloned nauplii [38]. Taxonomic status of all sexual species were re-confirmed by applying NCBI’s BLAST online platform (https://blast.ncbi.nlm.nih.gov), using COI sequences datasets [63], to confirm Artemia bisexual species identities, as the exotic American A. franciscana has been widely introduced in many countries all over the world.

DNA extraction and sequencing

Total genomic DNA was extracted individually from two adult females of each species and each parthenogenetic lineage, using the Rapid Animal Genomic DNA Isolation Kit (Sangon Biotech Co., Ltd., Shanghai, China; NO. B518221). In the case of A. monica, two individual decapsulated eggs were used for DNA extraction. Total genomic DNA was amplified following MALBAC Single Cell Whole Genome Amplification Kit no. KT110700150 (Yikon Genomics Co., Ltd. Jiangsu, CN). The Microvolume Spectrophotometer (MaestroGen Inc., Hsinchu, Taiwan) was utilized to check the quantity of amplified DNA. From each individual 600 ng of the amplified DNA was pooled and employed to build a single paired end (2 × 150 bp) genomic library using the NEBNext® UltraTM II DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) [4, 64]. Next-generation sequencing (> 10 Gb) of the pooled DNA library was performed on one sequencing flow cell of a NovaSeq 6000 Illumina machine (Novogene Co., Tianjin, China).

Bioinformatics analysis

Quality control and sequence assembly

Adapter residues were removed from the sequencing data by Novogene Co., and only sequences consisting of both paired reads were used for further analyses. Quality control was performed with the software package FastQC [65], as in Asem et al. [61]. For mitogenome assembly of the Asian species and parthenogenetic lineages, the mitogenomes of both A. sinica (GenBank accession no. MK069595; [64]) and A. urmiana (GenBank accession no. MN240408; [61]) were used as reference sequences. The A. franciscana mitogenome (GenBank accession no. X69067; [66]) served as a reference sequence to assemble that of the A. salina and American species. Bowtie v2.2.9 software [67] and Geneious R9.1 software [68] were used for sequence mapping and reference-based assembly, respectively, with parameter settings as in Asem et al. [4]. To confirm mitogenome sequence validity (circular form), each mitogenome was assembled twice and considered in two different positions with a 5,000 bp difference.

Gene identification and annotation

The secondary structure and position of mitogenomic transfer RNA (tRNA) genes were determined with ARWEN online software using default parameters (http://130.235.244.92/ARWEN/). Nucleic acid folding and hybridization of the tRNA sequences were predicted using the mfold online platform (http://www.unafold.org/mfold/applications/rna-folding-form.php) [69]. Ribosomal RNA genes (rRNAs) and protein-coding genes (PCGs) were annotated, based on gene order in the reference mitogenomes, followed by NCBI's BLAST online platform (https://blast.ncbi.nlm.nih.gov), using default parameters. The BioEdit v.7.2 software [70] was used to help determine the orientation and position of rRNAs and PCGs, based on multiple sequence alignments between the reference mitogenomes and all examined sequences in this study. Finally, the ExPASy online program (https://web.expasy.org/translate/) was utilized (setting: “invertebrate mitochondrial” codes) to translate PCG sequences into amino acid sequences, to ensure each PCG can encode a functional protein.

The nucleotide composition (AT% and GC%) and codon usage were determined with DAMBE 7 [71]. The AT- and GC- skews were also assessed, as described in Perna and Kocher [72]. Relative synonymous codon usage (RSCU) patterns variation between Artemia species and parthenogenetic lineages were determined with principal component analysis (PCA) using SPSS v18 software.

MEGA X software [73] was employed to calculate the percentage of variable sites (VS), as well as the pairwise interspecific genetic distances (D) and nucleotide diversity (ND) for each rRNA, PCG, and for a concatenated sequence of two rRNAs and 13 PCGs from both Artemia species and parthenogenetic lineages, according to the uncorrected p-distance nucleotide model ([74], Meier, per. com). Heat map values were illustrated using the Plotly online software (https://chart-studio.plotly.com).

Phylogenetic analysis

Phylogenetic relationship analysis between Artemia members was performed on a concatenated dataset, including two rRNAs and 13 PCGs, using Maximum Likelihood (ML) and Bayesian Inference (BI) as implemented in RAxML-HPC BlackBox 8.2.12 and MrBayes on XSEDE 3.2.7a, respectively [75]. Both methods were conducted on the CIPRES Science Gateway online platform (https://www.phylo.org/portal2). Streptocephalus cafer (GenBank accession no. NC_046688; [76]) was chosen as an outgroup.

The best fitting nucleotide substitution model of DNA was calculated following MrModeltest 2.2 software [77]) The GTR + I + G was chosen as the best fit model for both methods (ML: bootstrap replicates: 1,000; BI: nst = 6, rates = invgamma; mcmcp ngen = 10,000,000 samplefreq = 100 nchains = 4, sump burnin = 25,000, sumt burnin = 25,000). FigTree v1.4.4 was used to visualize the phylogenetic trees [78]. For the ML bootstraps, values < 70 were regarded as low and ≥ 95 as high [79]. For the BI posterior probabilities, the values < 0.94 were considered low and ≥ 0.95 as high [80].

The SAW method (see [81]) was used to evaluate if “long branch length” was systematic error as “long branch attraction” (LBA) ([82, 83], Bergsten, per. com).

Divergence time estimation using BEAST

Manzi et al. [84] reported fossil A. salina from Cyprus without any evidence to its identity; it is only a part of a putative anostracan thorax. Artemia, Branchinectella media (Schmankewitsch, 1873) and Phallocryptus spinosa (Milne-Edwards, 1840) are hypersaline anostracans, and all are reported from the Mediterranean Basin, where they may coexist (see [85, 86]). There is no taxonomic evidence to accept Manzi’s finding is Artemia, let alone “Artemia salina”. Therefore, divergence time was calibrated at 145 Mya, based on the age of a fossil of Daphnia Müller, 1785 (Crustacea: Anomopoda) ([87], see also [22]). Bayesian tree reconstruction and divergence times were calculated using the software BEAST v1.10.1 [88]. We tested the “relaxed clock” and “strict clock”, with the following additional parameters: nucleotide substitution model = GTR with four rate categories; Gamma + Invariant Sites heterogeneity among species; Tree Prior = Yule process ([89, 90], see also [91]); Random starting tree; Ancestral State Reconstruction = Reconstruct states at ancestor and Tree Root [92]. XML files for all BEAST runs were created using BEAUti v1.10.1 [88].

Posterior probability distributions of parameters were calculated by Markov Chain Monte Carlo (MCMC) sampling [88]. All runs were combined after a 10% burn in using LogCombiner v1.10.1 [88]. TRACER v.1.7.2 [93] was used to verify the stationary distribution of acceptable MCMC steps mixing, and to ensure appropriate sampling of each parameter (i.e., effective sampling size (ESS) > 200) (https://beast.community/ess_tutorial, see also [22, 91, 94, 95]. TreeAnnotator v1.10.1 [88] was utilized to annotate the maximum clade credibility tree (see [22]). The phylogenetic tree and clade divergence times were visualized using FigTree v 1.4.0 [78].

The analysis was run using “strict clock” over 40 million generations (the first step that ESS > 200), taking samples every 1,000 generations using BEAST (for more information see Results and Discussion).

Results

Mitogenome organization and composition

The mitogenomes of the nine species and four parthenogenetic lineages represent typical circular DNA, including 22 tRNAs, 2 rRNAs, 13 PCGs and a noncoding control region exhibiting a total length ranging from 15,433 bp (SAL1, GenBank accession no.: OR423222) to 15,829 bp (FRA1, GenBank accession no.: OR423224), respectively. The gene arrangement is identical in all Artemia mitogenomes, where nine PCGs and 13 tRNAs are coded on the heavy (H-) strand while the other genes (both rRNAs, four PCGs, and nine tRNAs) are encoded on the light (L-) strand (Table S1A-V). The mitogenome size and nucleotide compositions PCGs + rRNAs sequences are given in Table 2.

Table 2 Detailed information of the mitogenome sequences from nine species and four parthenogenetic lineages of Artemia in this study (see Table 1)

Six start codons, including ATG, ATC, ATT, GTG, ATA and TTG, were recognized in Artemia mitogenomes. ATG codon was identified as common and rare start codons with a frequency of 48.52%. The protein-coding genes COX1, COX3, CYTB and ND1 were initiated by a common ATG codon in all species and parthenogenetic lineages. The rare TTG start codon was only identified in the ND5 and ND6 protein coding gene (Fig. 1, Tables S2, S3, S4). Two stop codons, TAA and TAG, and an incomplete T stop codon were identified with frequencies of 51.48%, 17.16%, and 31.36%, respectively. The TAA stop codon was identified in all ATP6 protein-coding genes and the incomplete T stop codon was observed in genes COX1, COX2, ND5 and ND4 of all Artemia mitogenomes (Fig. 1, Tables S2, S3, S4).

Fig. 1

Distribution of start and stop codon usages in the Artemia mitogenome. A frequency of start and stop codon usages; B) frequency of start codon usages in each PCG; C) frequency of stop codon usages in each PCG

Artemia mitogenomic tRNAs have the typical cloverleaf secondary structure, except tRNA-Ser_1, which lacks a D-arm structure (Fig. 2). The most and least conserved tRNA sequences belong to tRNA-Met (57 conserved sites vs. 65 total sites, 87.69%) and tRNA-Glu (38 conserved sites vs. 66 total sites, 57.58%), respectively. An abnormal structure was detected in the acceptor stem of tRNA-Ser₁ in some taxa (A. tibetiana, A. urmiana, A. amati, A. sorgeloosi, di- and triploid parthenogenetic lineages) due to a lack of linking between second pair nucleotides (G and A) using the ARWEN online software (Fig. 3A). The possibility of a secondary structure of tRNA-Ser₁ was reconsidered by using the mfold online platform, which demonstrated that the ARWEN predicted secondary structure could not be a valid folding. The acceptor stem secondary structure was rearranged and confirmed via the mfold online platform as shown in Fig. 3B. A similar condition was recognized in the tRNA-Leu₂ acceptor stem in the A. persimilis mitogenome (Fig. 4 A, B).

Fig. 2

The predicted secondary structure of tRNA genes in the Artemia mitogenomes of the 13 investigated taxa. Any alteration, deletion and insertion has been labeled (abbreviations listed in Table 1)

Fig. 3

The predicted secondary structure of the acceptor stem of tRNA-Ser₁ (A. tibetiana, A. urmiana, A. amati, A. sorgeloosi, di- and triploid parthenogenetic lineages). A) output of the ARWEN software; B) the revised structure (the red rectangle and green polygon show unlinked nucleotides and the position of the revised acceptor stem in the secondary structure, respectively)

Fig. 4

The predicted secondary structure of the acceptor stem of tRNA-Leu₂ of A. persimilis. A output of the ARWEN software, B) the revised structure (the red rectangle and green polygon show unlinked nucleotides and the position of the revised acceptor stem in the secondary structure, respectively)

The nucleotide compositions of the concatenated PCGs + rRNAs sequences are depicted in Fig. 5 (Supplementary material Table S5). GC% and AT% values demonstrate that Artemia members cluster in three groups. Group A members are American origin Artemia (A. monica, A. franciscana and A. persimilis) and some Asian origin Artemia (A. sinica, tetra- and pentaploid parthenogenetic lineages) (Fig. 5B). All Group B members are Artemia of Asian origin (A. urmiana, A. tibetiana, A. sorgeloosi, A. amati, di- and triploid parthenogenetic lineages) (Fig. 5C). The Mediterranean A. salina is in a distinct group alone (Group C) (Fig. 5A). The GC- and AT-skew values are negative in all PCGs + rRNAs sequences (Supplementary material Table S6). Although Artemia species and parthenogenetic members have a wide distribution based on GC- and AT-skews, A. persimilis and A. salina are notably more isolated than others (Fig. 5D).

Fig. 5

Distribution of Artemia taxa on scatter plots based on A; B; C) GC%- and AT% and D) GC- and AT-skew values of PCGs + rRNAs (abbreviations listed in Table 2)

Thirteen shared protein-coding genes from each individual (totaling 26 individuals, including two specimens each of 13 species and parthenogenetic lineages) were used to analyze the codon usage of Artemia mitogenomes. UUU (Phe, F) and AUU (Ile, I) were the most representative codons, with a total of 5,462 (ranging from 173 to 228) and 5,328 records (ranging from 135 to 223), respectively. AGG (Ser, S) was the least common with 158 records (ranging from 1 to 8). These most representative and least commonly used codons were noted in A. persimilis: UCU (Ser, S) with a usage count of 141 in both PER1 (RSCU: 3.00) and PER2 (RSCU:2.99) and AGG (Ser, S) with a usage count of 1 in PER1 (RSCU: 0.021) and 2 in PER2 (RSCU: 0.042) (Fig. 6). Supplementary material Tables S7–S9 provides details of RSCU listed for each individual of each species and parthenogenetic lineage. The PCA of RSCU is shown in Fig. 7. The first and second components represent 50.61% and 20.87% of the variation, respectively, with AUC (0.935), AUU (−0.935) and GUA (−0.920) contributing to the first component and CAC (0.919), CAU (−0.919) and GCG (−0.825) contributing to the second component. PCA showed four groups, with those Artemia of Asian origin clustered in two relatively homogenous groups (Groups B and D). American origin species clustered together (Group A), while the Mediterranean A. salina was significantly isolated as a distinct group (Group C).

Fig. 6

Heat-map values for “number of usages per codon”, based on 13 PCGs of the Artemia mitogenomes

Fig. 7

PCA plot based on RSCU value for 13 PCG of the Artemia mitogenomes (abbreviations listed in Table 2)

Intraspecific diversity and genetic distance

The percentage of variable sites (VS) and the nucleotide diversity (ND) results are summarized for each PCG, rRNA gene, and concatenated sequences of PCGs and rRNAs in Fig. 8. There is likely a positive correlation between the two parameters, as they seem to move together. The highest and lowest values of VS and ND belong to ATP8 / ND6 and COX1 / COX2 protein-coding genes, respectively (Fig. 8 and Table S10). The pairwise interspecific genetic distances based on 13 PCGs, two rRNAs and Artemia mitogenomes (PCG + rRNA genes) are shown in Fig. 9 (see also Table S11). The greatest distances were found in A. persimilis and A. salina for most genes, except for ATP8, ND4L and ND6 PCGs and 12S rRNA genes, where the highest values referred to the genetic distance between A. persimilis and A. franciscana/A. monica, di-/triploid parthenogenetic lineages, A. sorgeloosi and the pentaploid parthenogenetic lineage, respectively. The lowest values were shared between diploid and triploid lineages (ND2, COX2, ATP8, COX3, ND3, ND5, ND4L, ND6, 16S, 12S, PCG + rRNA genes), and between tetraploid and pentaploid lineages (COX1, ATP6, ND4, CYTB, ND1).

Fig. 8

Percentage of variable sites and nucleotide diversity for each PCG, rRNA and PCGs + rRNAs genes

Fig. 9

Heat-map values of interspecific genetic distances based on 13 PCGs and two rRNAs separately and concatenated sequence of PCG + rRNA genes (X and Y axes are referred to Artemia taxa, abbreviations listed in Table 1)

Phylogeny and origin

Figure 10A represents the phylogenetic relationships among all species and four parthenogenetic lineages using mitogenome sequences. Both reconstructed trees (ML and BI) show the same topology, with A. persimilis located as a basal clade. The phylogenetic trees contain a “long branch length” with A. salina (Fig. 10A). Following the SAW method (see Materials and Methods) the observed “long branch length” should be considered as a true organismal phylogeny, not a long branch attraction (LBA) as a form of systematic error (Figs. 10B-E), with A. salina sister to the North American Artemia, which includes separate, distinct, and well-supported clades for A. monica and A. franciscana. Asian Artemia is divided in two major clades. The first clade contains A. sinica, plus tetra- and pentaploid parthenogenetic lineages. The second clade consists of A. sorgeloosi, A. urmiana, A. amati, A. tibetiana, and di- and triploid parthenogenetic lineages. Artemia sinica and A. sorgeloosi are sisters to the Asian clades, respectively. It seems that the observed difference in the topology of the mitochondrial phylogeny trees is due to the neglect of the “long branch length” in previous studies (see introduction).

Fig. 10

Mitogenomic phylogeny of Artemia based on Bayesian inference (BI) and Maximum-Likelihood (ML). A) using all Artemia taxa B) following the SAW method, the A. persimilis clade is not considered C) following the SAW method, the North American clade not considered D) long branch length (A. salina clade) not considered (Bergsten, personal communication, 2023) E) Asian clade not considered (Bergsten, personal communication, 2023). Numbers above the branches in red show the length of major branches (left: BI and right: ML). The number behind major nodes denotes posterior probabilities. The Bayesian support (left) and ML bootstrap values (right) are shown for each major node. long branch length is indicated by a red branch (abbreviations listed in Table 1)

Divergence time

Ten setting analyses (10 to 100 million generations, taking samples every 1,000 generations) were performed following the “relaxed clock” method. In all ten analyses, the ESS for several parameters was lower than 200 (most of them were also lower than 100; see Materials and Methods). The topology of estimated trees using “relaxed clock” method were unreliable and displayed a false inference, so that A. sorgeloosi was located as a basal clade for American Artemia. Additionally, A. sinica and tetra-/pentaploid parthenogenetic lineages were located between them (e.g., see Fig. S1). This is in contrast with the results of phylogeny and genetic distance between Artemia members in the present study and former studies. Therefore the “strict clock” method was chosen to estimate divergence time (for more information see Discussion).

Following the result of “strict clock” method, the diversification of Artemia began ca. 33.97 Mya (30.74–37.27 Mya; node B) in the Paleogene Period (Early Oligocene–Late Eocene) following the divergence of the ancestral lineage of A. persimilis from the other taxa (Fig. 11). The second division occurred between ancestral lineages of Mediterranean A. salina and the other taxa ca. 26.32 Mya (23.57–28.96 Mya; node C) during the Paleogene Period (Early Miocene–Early Oligocene). The divergence between North American and Asian taxa took place in the Early Miocene, ca. 21.10 Mya (18.86–23.32 Mya; node D). Based on our calibration, diversification in common ancestors of North American (A. monica and A. franciscana) and Asian Artemia appears to have occurred ca.0.44 Mya (0.39–0.54 Mya) in the Late Pleistocene (node E) and ca. 14.27 Mya (12.73–15.96 Mya; node F) in the Middle Miocene, respectively. The second Asian clade formed ca. 4.85 Mya (4.23–5.49 Mya; node H), mainly in the Early Pliocene. Although Pleistocene diversification of the tetraploid parthenogenetic lineage was earlier than the pentaploid parthenogenetic lineage (Holocene) (0.01–0.15 Mya vs. 0.00–0.008 Mya), diversification of diploid and triploid parthenogenetic lineages occurred in the Holocene at the same time (0.0001–0.01 Mya, 0.0008–0.01 Mya, respectively). The data indicates that initial diversification to A. sinica, A. franciscana, A. salina, and di-/tri-/pentaploid parthenogens occurred during the Holocene, while divergence of other taxa (consisting of A. urmiana, A. tibetiana, A. Amati, A. sorgeloosi, A. persimilis, A. monica and tetraploid parthenogens) occurred during the Late Pleistocene.

Fig. 11

A chronogram for the genus Artemia using the mitogenome. The blue node bars indicate 95% posterior probability intervals. (Mya: million years ago; HPD: highest posterior density, ^a nsotes to lower and upper 95% HPD intervals, abbreviations listed in Table 1). Geologic time scale was determined following the U.S. Geological Survey (https://www.usgs.gov/)

Discussion

Mitogenome organization and composition

We sequenced and assembled the mitogenomes of all described Artemia species and parthenogenetic lineages (with their ploidies identified), assuming that the complete mitogenome will provide better phylogenetic resolution and divergence time relative to traditional single mitochondrial or nuclear markers previously considered [51, 52].

The Artemia mitochondrial gene order is uniform with the ancestral Pancrustacea model, including 22 tRNAs, 2 rRNAs and 13 PCGs [96]. However, a slight variation was observed in mitogenome length (15,433 vs. 15,829 bp), though concatenated sequences of PCGs and rRNAs exhibit a conserved size (12,398 vs. 12,426 bp). Previous studies showed that the control region, regulating mtDNA replication and transcription [97], is less conserved in Asian Artemia [4, 61]. Our results show that sequences of this region differ by 19.4% among Artemia species (1,467 vs. 1,820 bp). Thus, we conclude that the Artemia control region is not conserved. This mitochondrial variability in the cell “energy factory” may potentially influence the expression of key genes [4, 61, 98], and the metabolic paths generating energy for Artemia adaptation to the harsh hypersaline conditions [4]. The role of mitochondrial genes in the Artemia adaptation to the harsh conditions of hypersaline environments needs further experimental assessment.

The mitogenome nucleotide composition differs across taxa [72, 99,100,101,102], likely due to differential selective pressures of individuals adapted to ecologically heterogeneous and stringent environments, with habitat-specific mutation rates and DNA repair mechanisms (see [23, 101, 102]). The ecological barrier between A. monica and A. franciscana justifies their reproductive isolation, criteria used for considering them as independent species following the “Biological Species” concept. Our findings reveal clear mitogenome composition differentiation (AT-/GC- content, RSCU and structure of tRNA-Gly, tRNA-Arg and tRNA-Asn) of A. monica and A. franciscana.

Intraspecific diversity and genetic distance

Asem et al. [61] demonstrated that the ATP8 protein-coding gene, involved in ATP generation in the respiratory chain, exhibits a significant difference in nucleotide composition compared to other genes and ribosomal rRNA genes, and showed that the ATP8 gene has the highest percentage of variable sites and nucleotide diversity among PCGs and rRNAs genes among Artemia members. The ATP8 gene exhibits the highest genetic distance between A. salina and A. monica/A. franciscana (D = 0.451). In contrast, the COX1 protein-coding gene, also playing a role in the respiratory chain, specifically in complexes III and IV, exhibits the lowest intraspecific diversity and genetic distance among Artemia members and bears a similar structure between tetra- and pentaploid parthenogenetic lineages (D = 0). Therefore, ATP8 and COX1 genes are identified as non-conserved and conserved Artemia mitochondrial genes, respectively, recognizing, however, that they are part of a more complex structure, function and interaction network, including nuclear control [60]. This mitogenomic comparison among Artemia taxa provides a more realistic approach to the mitochondrial complexity and functionality than those based on single mitochondrial genes.

Phylogeny, origin and divergence time

The “relaxed clock” method allows each lineage and branch of a phylogenetic tree to have its own evolutionary rate [103, 104] and is used when the rate variations among lineages/clades are high [104]. However, “strict clock” is suitable for low-rate variation with shallow phylogeny [105].

According to the results of the “relaxed clock” method, the phylogenetic trees are unreliable, whereas A. sorgeloosi was located as a basal clade for American Artemia. Furthermore, A. sinica and tetra-/pentaploid parthenogenetic lineages were located. This problem is also observed in using “relaxed clock” to estimate divergence time in the previous study, where A. franciscana is located as a basal clade and A. salina is placed between Asian Artemia and A. franciscana (see [22]). Due to low-rate variation among North American and Asian Artemia, and the results of the estimated tree using a “relaxed clock” model, the “strict clock” was utilized as correct model for our data set [Brown per. com.; Yang per. com.].

Because Artemia species and lineages are regionally distributed and adapted to harsh and ecologically variable conditions, coupled with an island biogeography-type distribution of intraspecific diversity [3,4,5, 27], they represent unique conditions for studying evolutionary divergence. Such evolutionary change is likely to be accelerated or affected by climate change and anthropogenic impacts [106, 107]. Consequently, some species’ taxonomic status, particularly those in places difficult to access (Tibet, Kazakhstan, etc.) or less studied, is subject to debate [4, 5]. One problem is the uncoupling of evolutionary change at morphological and molecular levels, reflected by the occurrence of “sibling species”, e.g., morphologically similar but genetically divergent forms, which requires calibrating gene-level molecular markers (mitochondrial and nuclear) with different mutation rates (conserved, less conserved). In the late 1980s, Artemia had been confusingly referred to as a genus containing “sibling species” [108,109,110]. However, they are distinguishable morphologically using laboratory-cultured individuals [see 4]. Further confusion arose because most species were described typologically in the nineteenth and twentieth centuries, which led to the idea that Artemia was monotypic, consisting only of “A. salina” (see [111, 112]). This was exacerbated by toxicologists using “A. salina” as a test organism, without any information on the origins of their test subjects (see [2, 113]). More recently, the biological species concept (BSC) ([114], see also [115]) was applied to bisexual species of Artemia taxa reproductively isolated in nature [see 4]. However, the BSC has obvious limitations in Artemia, a genus with obligate parthenogenetic lineages [116]. Thus, our results provide a comprehensive mitogenomic phylogenetic comparison and divergence of Artemia, clustering the members into well-supported clades.

Artemia urmiana possesses significantly high mitochondrial and nuclear intraspecific genetic variation [4, 34]. Our results suggest that A. urmiana is the earliest established Artemia species and its intraspecific divergence dates back to ca. 0.44 Mya (0.36–0.54 Mya) in the Late Pleistocene. This explains the high genetic variation and population expansion of A. urmiana.

Sainz-Escudero et al. [117] argued “based on the occurrence of nuclear gene flow between the type locality of A. tibetiana and populations of A. urmiana” regarding the results of Maccari et al. [19], that A. tibetiana should be considered as a junior synonym of A. urmiana. Maccari et al. [19] analyzed two nuclear markers, ITS1 and a fragment of exon-7 of Na ⁺ /K ⁺ ATPase, but Sainz-Escudero et al. did not state which nuclear markers were affected by gene flow. Asem et al. [38] demonstrated that the fragment of exon-7 of Na ⁺ /K ⁺ ATPase is a highly conserved marker in Asian Artemia and cannot be considered a phylogenetic marker. Contrary to Sainz-Escudero et al.’s [117] claim, there are no shared haplotypes between ITS1 markers of topotype A. tibetiana (GenBank accession no. KF736290,91; [19]) and populations of A. urmiana (KF736249,52; [19]). Asem et al. [4] demonstrated that ITS1 is not an informative marker for Asian Artemia phylogeny. Additionally, nuclear gene pool isolation between A. tibetiana and A. urmiana has been previously documented by Nougué et al. [48] using nine SSR markers [see also 4,49]. Furthermore, morphological differentiation of A. urmiana in field and laboratory collections with a cercopod rudimentary/oligosetal pattern is morphological evidence to distinguish A. urmiana from other Asian species and parthenogenetic lineages (4, see also [116]). Although several experimental cross-studies showed weak or lack of mating isolating among Asian species, Asem et al. [4] argued that laboratory crossbreeding tests cannot prove fertility and/or infertility potential in nature. Following Sainz-Escudero et al. [117], Li et al. [118] treated A. urmiana and certain Tibetan populations (A. tibetiana and A. sorgeloosi) as the “Artemia urmiana species complex”, yet at the same time demonstrated that A. urmiana is completely isolated from Tibetan populations based on nuclear SSR markers. The species complex designation was explained as “a neutral term for a number of related taxonomic units, most commonly involving units in which the taxonomy is difficult or confusing” following Mayr and Ashlock [119]. A species complex is not a taxonomic rank or unit of organism classification recognized by the International Commission on Zoological Nomenclature (ICZN) (see [119]). Our mitogenome phylogenetic analysis demonstrates that A. urmiana and A. tibetiana each diverged from different ancestors in different geologic periods (ca. 1.32 Mya in the Early Pleistocene and ca. 0.37 Mya in the Late Pleistocene, respectively), by which they are linked via three noncommon ancestors. On the other hand, A. urmiana/A. sorgeloosi and A. tibetiana/A. sorgeloosi are connected via two and four noncommon ancestors, respectively. In conclusion, due to isolated SSR nuclear gene pools [4, 48, 49] and maternal gene pools [4, 5, 17,18,19, 21], and lack of common maternal ancestor, synonymisation of A. tibetiana or A. sorgeloosi with A. urmiana is not supported. Neither do they represent a “species complex”.

A recent Asian Artemia phylogenetic study inferred a lack of gene flow and probable reproductive isolation in nature due to the existence of private maternal haplotypes [4, 5]. Despite this, A. urmiana, A. tibetiana, and A. amati exhibited low genetic distances based on partial mitochondrial markers (COI, 16S and 12S) in comparison with other anostracan families (e.g. [120]). The low genetic distances within Artemia are likely related to its relatively recent divergence times compared to other anostracan genera (see [22, 121,122,123], also reviewed in 3,4). Our dating results show that the Asian Artemia ancestor diverged more recently (12.73–15.96 Mya, node F) than American and Mediterranean ancestors (30.74–37.27 Mya, node B; 23.57–28.96 Mya, node C; 18.86–23.32 Mya, node D). Reanalysis of intraspecific diversity using the mitogenome data refined the genetic differentiation resolution among A. urmiana, A. tibetiana and A. amati (A. urmiana/A. tibetiana: D = 0.026; A. urmiana/A. amati: D = 0.027; A. tibetiana/A. amati: D = 0.007). The low genetic distance between A. tibetiana and A. amati can also be attributed to the recent divergence of ca. 0.37 Mya (0.29–0.45 Mya, node M) in the Late Pleistocene (see below). The significant differentiation between A. tibetiana and A. amati nuclear genomes based on SSR markers is remarkable [4, 48, 49]. Maccari et al. [19], Sainz-Escudero et al. [117] and Li et al. [118] argued the taxonomic status of Asian Artemia using only close genetic distances, but ignored isolation based on morphological differentiation and molecular (mitochondrial and nuclear SSR) markers. Although “genetic distance” is an important molecular tool in intra- and interspecific studies, it alone cannot be used as evidence to determine taxonomic status within Artemia.

Asem et al. [4] demonstrated that A. tibetiana and A. sorgeloosi are monophyletic based on nuclear genome (SSR markers) but maternally polyphyletic. They also concluded that the maternal divergence time of the A. tibetiana clade was after A. sorgeloosi, suggesting that A. tibetiana may have originated from a past ancestral hybridization of a maternal ancestor of A. tibetiana with A. sorgeloosi or its ancestor [4]. Our phylogenetic results confirm that each taxon arose from a distinct maternal ancestor. We also show that the maternal divergence time of the A. tibetiana lineage originated in the Late Pleistocene (ca. 0.37 Mya, node M) after the A. sorgeloosi lineage in the Early Pliocene/Late Miocene (ca. 4.85 Mya, node H), while within A. sorgeloosi, diversification occurred later (ca. 0.08 Mya) than in A. tibetiana (ca. 0.14 Mya). Thus, A. sorgeloosi could not have a direct evolutionary function in the origination of A. tibetiana. Therefore, we have revised our previous hypothesis [see 4]: we think that A. tibetiana may have evolved from a past ancestral hybridization via a maternal ancestor of A. tibetiana with a paternal ancestor of A. sorgeloosi on the Tibet Plateau.

As previously mentioned, Artemia monica and A. franciscana are two independent species, ecologically and reproductively isolated (see [23]), which is confirmed by the mitogenomic phylogeny showing two well-supported clades (BI/ML = 1/100) and separated from a common ancestor ca. 0.44 Mya (0.36–0.54 Mya, node E) in the Late Pleistocene. According to the calibration performed, A. monica and A. franciscana would have diverged at different geologic times; the Late Pleistocene (0.13 Mya) and the Holocene (0.002 Mya), respectively. Therefore, it is necessary to highlight their specific status based on ecological and potential reproductive isolation (migrants are selected against, as demonstrated by Bowen [124]. Following our result, Artemia franciscana would have evolved more recently.

Artemia also consists of obligate parthenogenetic forms, independent of the Artemia species [116]. Although for several decades, parthenogenetic Artemia forms have been treated as A. parthenogenetica, Asem et al. [116] show it is an invalid binominal specific name, suggesting instead to consider them as “parthenogenetic lineage(s) of Artemia”. The evolutionary relationship and origin(s) of Artemia parthenogenetic lineages was explored using electrophoretic markers by Abreu-Grobois and Beardmore [114] in American and Mediterranean species, including parthenogenetic lineages with different ploidy levels. However, as Asian species were not included, all parthenogenetic lineages were clustered in one major clade with Mediterranean A. salina being the sister clade. They concluded that: 1) the diploid parthenogenetic lineage has arisen from Mediterranean A. salina; 2) the tetraploid parthenogenetic lineage diverged from diploid parthenogenetic lineages; 3) the pentaploid parthenogenetic lineage originated from tetraploid parthenogenetic lineages, and; 4) the triploid parthenogenetic lineage could have arisen independently from diploid parthenogenetic lineages. Later, Beardmore and Abreu-Grobois [42] included the Asian species A. urmiana and revised their previous hypothesis, stating that parthenogenetic lineages arose from A. urmiana. Both studies introduced parthenogenetic lineages as a monophyletic group. Nevertheless, a comprehensive study based on four ploidy levels of parthenogens and using three mitochondrial (maternal) markers (COI, 16S and 12S) revealed that parthenogenetic Artemia clustered in four distinct and well supported clades, each with different origins ([38], see also [18]). Our results also demonstrate that parthenogenetic Artemia lineages are a maternally polyphyletic group. Previously, it was suggested that diploid parthenogenetic lineages originated from A. amati ([17]; for more information see [49]) (A. amati was not described at that time [4]) and/or A. urmiana ([18, 38] see also [49]). Our mitogenome-based phylogenetic tree does not show a common ancestor between di-/triploid parthenogenetic lineages and A. urmiana. Thus, A. urmiana did not play a role in the origin of the diploid parthenogenetic lineage as previously suggested (see below). Our mitogenomic evidence shows that the common ancestor of A. amati/A. tibetiana (ca. 0.37 My, node M) and the common ancestor of diploid and triploid parthenogenetic lineages (ca. 0.07 Mya, node L) originated from a common historical ancestor (ca. 0.61 Mya, node K) in the Late Pleistocene. Thus, there is no direct maternal evolutionary link between A. amati and the diploid or triploid parthenogenetic lineages. We hypothesize that the origin of the diploid parthenogenetic lineage could be linked to mutational event(s) in the genome of the historical ancestor of diploid parthenogens. Regretfully, the historical correlate of the diploid parthenogenetic lineage is unclear.

Rode et al. [49] recognized two groups of diploid parthenogenetic Artemia (“diploid parthenogenetic Artemia with A. urmiana type mitochondria” and “diploid parthenogenetic Artemia with Artemia sp. Kazakhstan [later described as A. amati] type mitochondria”) using a partial sequence (< 700 bp) of COI haplotype network distribution. However, haplotype networks cannot reconstruct phylogenies (McFadden, per. com., 2022) or ancestral evolutionary relationships among taxa. Asem et al. [4] used three maternal markers (COI, 16S and 12S) to illustrate haplotype network distributions among Asian species. Using 16S sequences, A. tibetiana was located between A. urmiana and A. amati (Artemia sp. Kazakhstan in ([49], for more information see 4], but in contrast, using a 12S sequences haplotype network, A. amati was located between A. urmiana and A. tibetiana. Furthermore, A. tibetiana and A. amati have parallel positions in connection to A. urmiana, based on a COI haplotype network. The same problem can be found in multivariate analysis plots (PCoA, PCA, DA, etc.) using SSR/ISSR markers and morphological characters, where cluster position(s) can be changed with an increasing/decreasing number of studied groups. However, multivariant analysis plots cannot reconstruct ancestral evolutionary relationships and phylogenies. Both haplotype network distribution and multivariant analysis are useful bioinformatic tools to consider gene flow between taxa but are not reliable for phylogeny. Given that “diploid parthenogenetic Artemia with A. urmiana type mitochondria” [see [49]) was not analyzed in our study, its placement in the mitogenomic phylogeny is still open.

Generally, five hypotheses would explain the origin of the triploid parthenogenetic lineage: 1) fertilization of an unreduced ovum (secondary oocyte) from a diploid parthenogenetic Artemia by a sperm cell of A. urmiana ([18], see also [49]; 2) an unreduced A. urmiana ovum by sperm of a rare diploid parthenogen male [18]; 3) an unreduced A. urmiana ovum with A. urmiana sperm; 4) via polyspermy of a normal A. urmiana ovum, and; 5) via fertilization of two reduced nuclei (polar bodies) and a reduced ovum of a diploid parthenogen [38]. As discussed above, there is no common maternal ancestor between A. urmiana (and A. amati) and the triploid parthenogenetic lineage. Thus, an unreduced A. urmiana ovum (and an unreduced A. amati ovum) would not have had a maternal evolutionary function in the origin of the triploid parthenogenetic lineage. Following our mitogenomic phylogeny tree (see also [18]) and previous studies on SSR nuclear markers [49], triploid parthenogens could have originated through: 1) fusion of a normal ovum of diploid parthenogen with its two polar bodies 2) polyspermy of normal ovum of a diploid parthenogenetic Artemia via sperm of its own rare male(s), or; 3) fertilization of a diploid parthenogenetic Artemia unreduced ovum (secondary oocyte) with its rare male.

Consistent with previous studies [18, 38], our mitogenomic phylogeny demonstrates that A. sinica shares a historical maternal ancestor (ca. 0.96 Mya, in Early Pleistocene) with the common ancestor (ca. 0.05 Mya, Late Pleistocene) of tetraploid and pentaploid parthenogenetic lineages.

Rode et al. [49] proposed that tetraploid parthenogens could have originated by hybridization between A. sinica (maternal origin) and a rare male of a diploid parthenogenetic lineage (paternal origin). In contrast to Rode et al. [49], who referred to the historical hybridization event in East Siberia using a single sequence named “A. sinica (?)” (LC195586; [125]), we suggest that hybridization between A. sinica and a rare male of a diploid parthenogenetic Artemia should have occurred in East Asia. Biogeographically, A. sinica only occurs in East Asia, and there is no taxonomic evidence that LC195586 (TU13), obtained from a batch including two females from Lake Dus-Kholin in Siberia (see [125]), belongs to a bisexual population and/or A. sinica. Additionally, the sequence of TU13 most likely is a tetraploid (or pentaploid) parthenogenetic lineage, whereas the sequence quality and taxonomic status of specimens studied by Naganawa and Mura [125] are problematic and questionable (see [4, 34, 61, 63]). It is most unlikely that TU13 could be an introduced exotic specimen of A. sinica.

Asem et al. [38] suggested that the pentaploid parthenogenetic lineage may have evolved from an allopolyploid (hybridization) of A. sinica sperm and a tetraploid ovum or autopolyploid. The lack of gene exchange between the nuclear genes of A. sinica and pentaploid parthenogens (see [49]), cast doubts on the role of A. sinica in the origin of the pentaploid parthenogenetic lineage. Therefore, the autopolyploid or allopolyploid (historical hybridization with the paternal source of A. sorgeloosi and/or A. tibetiana) is most likely the origin of the pentaploid parthenogenetic lineage (see also [49, 4]).

Our mitogenomic study clarifies the phylogenetic tree topography controversies of previous studies. For example, the initial allozyme study [42] put the New World species together in a major clade, with A. salina and Asian species in another. Later, studies using the partial mitochondrial COI put either the New World A. persimilis [18] or the Mediterranean A. salina in a basal clade [21]. However, an abnormal phylogenetic tree using COI by Muñoz et al.’s [16] had the New World A. franciscana between A. sinica and other Asian members. Later trees based on the nuclear ITSI marker showed A. salina between the New World A. franciscana (A. monica not analyzed) and South American A. persimilis [18, 21, 38]. According to our results, the maternal ancestral clade of Artemia split into ancestral lineages of A. persimilis and A. salina in ca. 33.97 Mya in the Early Oligocene/Late Eocene, which should have originated in the Mediterranean area or South America. Although the Mediterranean area can be assumed as one of the possible geological origins of Artemia, the proposition of the Mediterranean area as the origin of Artemia divergence, as suggested by Beardmore and Abreu-Grobois [42], due to the extreme salinity rise there, now seems debatable. The Messinian salinity crisis in the Mediterranean basin occurred 6–5.3 Mya [126]. Therefore, this event cannot be referred to as the cause of Artemia diversification, which, according to our results, occurred ca. 33.97 Mya. Even so, the Messinian salinity crisis has probably played a role in the expansion of A. salina in the Mediterranean region.

Artemia diversification in the New World continued from the ancestral lineage of A. salina in the Mediterranean area via a common ancestor of A. monica and A. franciscana (ca. 21.10 Mya), whereas Asian Artemia diverged through a maternal ancestor lineage of A. sinica from East Asia (ca. 14.27 Mya). East Asia was the maternal origin of the tetraploid (ca. 0.1 Mya) and pentaploid parthenogenetic lineages (ca. 0.002 Mya). The Tibetan Plateau was the dispersal bridge for Artemia from East to West Asia [for the origin of A. urmiana; see [50, 127] via the maternal ancestor lineage of A. sorgeloosi in ca. 4.85 Mya. The Tibetan Plateau and Central Asia are recent scenarios of Artemia divergence (A. tibetiana and A. amati, respectively) through the maternal ancestor lineage of A. urmiana, which occurred ca. 1.3 Mya. Although the geologic origin(s) of diploid and triploid parthenogenetic lineages are unclear, their origin(s) should be in West and/or Central Asia.

Conclusion

The brine shrimp Artemia is a paradigmatic crustacean genus adapted to the variable hypersaline conditions, with sexual species and parthenogenetic lineages, whose evolutionary relationships and divergence times are updated by this study based on the whole mitochondrial genome. Relevant conclusions are:

• 1. Complete mitochondrial genomes (mitogenomes) can provide a clear phylogeny relationship and clarify common ancestor(s).

• 2. The Artemia diversification occurred from the ancestral lineages of A. salina in the Mediterranean area and A. persimilis in South America (ca. 33.97 Mya), followed by the ancestral lineage of “North America and Asia” (ca. 26.32 Mya). North America Artemia diverged in ca. 0.44 Mya and Asian Artemia in East Asia in ca. 14.27 Mya. East Asia was the maternal origin of the tetraploid (ca. 0.1 Mya) and pentaploid parthenogenetic lineages (ca. 0.002 Mya). Geological origins of diploid (ca. 0.0005 Mya) and triploid (ca. 0.007 Mus) lineages probably developed in West and/or Central Asia.

• 3. Asian Artemia originated from multiple ancestors that diverged in different geologic periods. Therefore, unlike previous studies, Asian species cannot be considered “species complex.”

• 4. The mitogenomic phylogenetic tree shows no direct association between maternal sources of A. amati and the diploid parthenogenetic lineage from Gahai Lake (“diploid parthenogenetic Artemia with Artemia sp. Kazakhstan type mitochondria” as previously reported [49]). They are, however, linked via three noncommon ancestors in our analysis. Therefore, maternal sources of diploid parthenogenetic Artemia cannot be traced back to the mitochondria of A. amati, which has its own independent mitogenome evolutionary origin. Moreover, our findings exhibit an identical evolutionary status for A. amati and A. tibetiana in relation to diploid parthenogenetic Artemia.

• 5. The taxonomic status of A. monica and A. franciscana has been controversial; however, in addition to the evidence of ecological and reproductive isolation, differentiation in tRNAs secondary structure, nucleotide compositions, RSCU, and especially intraspecific divergence time and geologic occurrence periods, A. monica and A. franciscana are demonstrably two distinct species (see also [23]).

• 6. Our case study shows that parthenogenetic Artemia is a maternally polyphyletic group, and that the tetraploid lineage divided in a different geologic time (Late Pleistocene vs. Holocene). Based on zoological species concepts and on Artemia polyphyletic evolutionary relationships, parthenogenetic lineages cannot be referred to as the binominal specific nomen “Artemia parthenogenetica”. The parthenogenetic lineages originated from several historical hybridization events, including crosses and backcrosses within Asian species.

• 7. Asian and North American Artemia are the most recently diversified taxa, and the value of genetic distance cannot be a reliable scale for species delimitation.