Reanalysis and Revision of the Complete Mitochondrial Genome of Artemia urmiana Günther, 1899 (Crustacea: Anostraca)

In the previously published mitochondrial genome sequence of Artemia urmiana (NC_021382 [JQ975176]), the taxonomic status of the examined Artemia had not been determined, due to parthenogenetic populations coexisting with A. urmiana in Urmia Lake. Additionally, NC_021382 [JQ975176] has been obtained with pooled cysts of Artemia (0.25 g cysts consists of 20,000–25,000 cysts), not a single specimen. With regard to coexisting populations in Urmia Lake, and intraand inter-specific variations in the pooled samples, NC_021382 [JQ975176] cannot be recommended as a valid sequence and any attempt to attribute it to A. urmiana or a parthenogenetic population is unreasonable. With the aid of next-generation sequencing methods, we characterized and assembled a complete mitochondrial genome of A. urmiana with defined taxonomic status. Our results reveal that in the previously published mitogenome (NC_021382 [JQ975176]), tRNA-Phe has been erroneously attributed to the heavy strand but it is encoded in the light strand. There was a major problem in the position of the ND5. It was extended over the tRNA-Phe, which is biologically incorrect. We have also identified a partial nucleotide sequence of 311 bp that was probably erroneously duplicated in the assembly of the control region of NC_021382 [JQ975176], which enlarges the control region length by 16%. This partial sequence could not be recognized in our assembled mitogenome as well as in 48 further examined specimens of A. urmiana. Although, only COX1 and 16S genes have been widely used for phylogenetic studies in Artemia, our findings reveal substantial differences in the nucleotide composition of some other genes (including ATP8, ATP6, ND3, ND6, ND1 and COX3) among Artemia species. It is suggested that these markers should be included in future phylogenetic


Introduction
Mitogenomes have been extensively utilized in phylogeny and population genetic studies because mitochondrial genes share some particular characteristics, including maternal origin, rapid evolutionary rates and lack of recombination [1,2]. Sequences of mitochondrial genes have been widely applied as informative molecular markers; therefore, they have been widely used in the molecular phylogenetic studies [3][4][5][6]. However, compared to partial mitochondrial sequences, such as COX1, 16S, 12S, etc., complete mitochondrial genome sequences could provide more of a higher resolution and provide better understanding about evolutionary relationships and structures of animals [4,7].
The number of complete mitochondrial sequences from metazoa has been rising rapidly in the past decade. With the rapid advancement in genomics techniques, such as next-generation sequencing (NGS) technology, many genomes and mitogenomes have been characterized [6,8,9]. The current methodology provides high quality results and has allowed us to accurately determine the gene order and organization in many taxa and to compare the outcomes with conventional sequencing methods [10][11][12][13]. Complete mitochondrial genomes have proved useful in solving evolutionary and biosystematics questions over a broad taxonomic range [6,14].
Overall, the phylogenetic relationship within Artemia is presently still ambiguous [22]. Obligate asexuality is one of the critical challenges in biosystematics of the genus Artemia. Accordingly, asexual forms of Artemia do not form defined species but parthenogenetic population(s) [18,23]. A comprehensive phylogenetic study using mitochondrial and nuclear markers provided evidence that the parthenogenetic Artemia taxa form a polyphyletic group. The diploids and triploids are maternally related to A. urmiana, whereas tetrapod and pentaploid lineages share a common ancestor with Artemia sinica [19]. Although the mitochondrial marker Cytochrome Oxidase Subunit I (COX1) has confirmed that parthenogenetic populations were divided into two polyphyletic groups, the nuclear marker ITS1 suggests a common haplotype with the participation of all ploidy degrees [19,24]. Additionally, Eimanifar et al. [24] showed that parthenogenetic populations share a major haplotype with A. urmiana and A. tibetiana based on the ITS1 nuclear maker.
Maniatsi et al. [25] and Eimanifar et al. [24] have proposed different phylogenetic positions for A. salina and A. persimilis, based on nucleotide sequences of COX1. However, results of morphological and genetic investigations were contradictory [22]. To date, the taxonomy and biosystematics of Artemia are still controversial, especially concerning the Asian species [21]. Mitogenomic information could provide a better reconstruction of the maternal evolutionary mechanism and phylogenetic status of Artemia. However, only four mitochondrial genomes of Artemia species have been published to date.
The previously published mitogenome analysis, performed by Zhang et al. [26] for two Artemia species, Artemia urmiana and Artemia tibetiana, detected an unprecedented unusual nucleotide sequence in the mitochondrial genome of A. urmiana from Urmia Lake, Iran. The length of the complete mitochondrial genome of A. urmiana (GenBank accession number NC_021382 [JQ975176]), is 15,945 bp; its control region would be significantly longer than other species. In addition, the taxonomic status of the assembled mitogenome had not been identified. For this examined mitogenome, Zhang et al. have used Artemia cysts stored at the Laboratory of Aquaculture & Artemia Reference Center (ARC, Gent University), which are labeled as ARC 1227 [26]. This sample has been collected in 1991 and it has been only registered as "from Urmia Lake" in ARC documents (Christ Mahieu, personal communication, ARC, Gent University).
Artemia urmiana has coexisted with a low percentage of parthenogenetic populations in the main body of Urmia Lake in western Iran, whereas coastal areas and neighboring lagoons of Urmia Lake hosted mostly parthenogenetic populations. Barigozzi and Baratelli [27] documented that all samples collected from Urmia Lake in 1987 were parthenogenetic and contained di-, tetra-, pentaploid individuals. However, Azari Takami [28] reported that A. urmiana and parthenogenetic populations coexisted in the lake. Asem et al. [29] identified two partly isolated parthenogenetic populations from the main body of Urmia Lake and lagoons neighboring the lake. Generally, some of the adjacent lagoons are connected with Urmia Lake when the lake level raises annually during rainy seasons in spring and autumn [29], which increases the probability of collecting parthenogenetic specimens along the shoreline of the lake (Atashbar, personal communication, Urmia University). On the other hand, genetic variation between parthenogenetic populations and A. urmiana is quite low [19,21]. Additionally, a current study based on barcoding with the mitochondrial COX1 marker found that parthenogenetic populations in some localities share same haplotypes with A. urmiana [30]. Based on this fact, the use of COX1 sequencing alone to distinguish A. urmiana from di-and triploids is questionable. However, these populations need to be further examined with special emphasis on the status of reproductive mode (bisexual or parthenogenesis) [31], morphological traits [29,32,33] and a Single Nucleotide Polymorphism (SNP) in the Na + /K + ATPase α-1 subunit [34] to distinguish and identify A. urmiana and parthenogenetic populations. The taxonomic status of Artemia from Urmia Lake should be identified before genome sequencing, see [21,35], which has not been considered by Zhang et al. [26]. There is a critical problem in Zhang et al. [26] mitogenome sequence that it has been procured by pooled cysts of Artemia (0.25 g cyst consists of 20,000-25,000 cysts) from Urmia Lake. It is a pivotal puzzle how differences among individuals have been ignored in assembling process. In this condition, the NC_021382 [JQ975176] mitogenome cannot be a valid sequence and any attempt to attribute it to A. urmiana or a parthenogenetic population is unreasonable.
Through the present study, we confirmed the taxonomic status of A. urmiana sampled from Urmia Lake. The purpose of the current research was to re-sequence the complete mitochondrial genome of A. urmiana using Illumina Hi-Seq X next generation sequencing technology with a high-coverage. For that reason, we first aimed to fully re-annotate and characterize the complete mitochondrial genome of A. urmiana and to compare it with previously published assembly by Zhang et al. [26]. Due to lack of consistency regarding the common molecular markers to study Artemia phylogeny, our results will contribute to understanding interspecific variation among different mitochondrial genes, and to identify mitochondrial markers with high variability, and the evolutionary origin of Asian Artemia.

Artemia Sample, DNA Extraction and Sequencing
The cysts of Artemia urmiana were collected from Urmia Lake, Iran (Kholman-khaneh Port; 45 • 29 E, 37 • 64 N) in 2004. The cysts were transferred to the laboratory and stored at the optimum storage condition until analysis. The cysts were processed and incubated at the standard hatching conditions according to procedures suggested by Sorgeloos et al. [36]. Specimens were identified as belonging to the A. urmiana, as confirmed by the reproductive mode (bisexual or parthenogenetic) [19] and a SNP in the Na + /K + ATPase α-1 subunit [34]. The reared specimen showed the rudimentary furca characterized by oligosetae, which is a main morphological character for A. urmiana [32], while parthenogenetic populations of Urmia Lake have a furca with two lobes and polysetae [29][30][31][32][33]. In order to distinguish it from other bisexual species of Artemia (especially invasive American A. franciscana), Diversity 2021, 13, 14 4 of 17 the taxonomic status of the chosen specimen was ultimately re-confirmed by applying the standard nucleotide BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) using COX1 sequences based on Eimanifar et al. [24] datasets.
Total genomic DNA was extracted from a female using a Rapid Animal Genomic DNA Isolation Kit (Sangon Biotech Co., Ltd., Shanghai, China; NO. B518221), following the manufacturer's instructions and procedures performed by Asem et al. [8]. The quality of extracted DNA was checked on a 1.5% agarose gel and then quantified using a Microvolume Spectrophotometer (MaestroGen Inc., Hsinchu City, Taiwan). An amount of 600 ng of total DNA was used to construct the genomic library with paired-end (2 × 150 bp) using NEBNext ® Ultra TM II DNA library preparation kit, followed by next-generation sequencing (10 Gb) approach. The sequencing was performed by a single constructed library, pooled on Illumina HiSeq X-ten sequencing flow cell (Novogene Co., Tianjin, China).

Sequence Quality Control, Assembly and Annotation
Adapter residues were removed before further analyses. The Illumina sequence reads were qualitatively checked using FastQC [37]. Assembly and coverage metrics are summarized in Figure S1 and Table S1.

Gene Identification and Annotation
The position and secondary structure of the tRNA genes were determined by AR-WEN (http://130.235.46.10/ARWEN/) online software. Protein-coding genes (PCGs) and ribosomal RNA genes (rRNAs) were annotated based on the gene order on the reference mitochondrial map using BLAST analysis (https://blast.ncbi.nlm.nih.gov). The position and orientation of the PCGs and rRNAs were recognized by the analysis of multiple sequence alignments to the reference mitogenome using BioEdit program [39]. In addition, all PCGs were translated into amino acids by the ExPASy online program (https://web.expasy.org/translate/), and sequences were examined to ensure that each could encode a functional protein.

Bioinformatics and Phylogenetic Analysis
GenomeVx online tool was utilized to draw the circular map of the mitochondrial genome of A. urmiana [40]. The complete mitochondrial sequences of A. sinica (MK069595), A. tibetiana (NC_021383) and A. franciscana (X69067) were downloaded from GenBank. The nucleotide composition and codon usage were calculated with DAMBE 6 [41], and AT-and GC-skew were also calculated using the following formulas: AT-skew = [A-T]/[A+T] and GC-skew = [G-C]/[G+C] [42]. The secondary structure of tRNAs were visualized using ARWEN online software. Sequences were aligned using MEGA X selected with MUSCLE default setting [43]. Pair-wise distances were computed for 13 PCGs and 2 rRNAs using the uncorrected p-distance nucleotide model as implemented in MEGA X [43]. Heat-maps of Euclidean distance among AT-and GC-skew were generated using the Heatmapper online tool [44].
Phylogenetic analyses of the concatenated datasets, including 13 PCGs and two rRNAs, were carried out using two different tree-building methods, Maximum Likelihood (ML) and Bayesian Inference (BI). The ML analysis was performed using MEGA X [43]. The BI was implemented in MrBayes 3.2.2 on XSEDE [45]. For ML and BI the best-fitting nucleotide substitution model was calculated based on the results of the MrModelltest 2.2 [46] test and HKY+G was chosen as the best-fit model (ML: bootstrap replicates: 1000, maximum parsimony analyses were run using TNT (Nearest-Neighbor-Interchange); BI: mcmcp ngen = 10,000,000, samplefreq = 100, nchains = 4, sump burnin = 25,000). The trees were visualized using FIGTREE v1.4.0 [47]. For the ML bootstraps, the values <70 were considered as low, 70-94 as middle, and ≥95 as high [48]. For the BI posterior probabilities, the values <0.94 were regarded as low, and ≥0.95 as high [49].

Complementary Experiments
To confirm the presence or absence of the 311 bp partial nucleotide in the mitochondrial sequence of A. urmiana, 48 further specimens have been examined following designed primers by Zhang et al. [26] to amplify this partial DNA region.
In order to consider sequence quality of the previously published mitogenome by Zhang et al. [26], the standard nucleotide BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was performed to check similarity of 16S rRNA sequence of NC_021382 [JQ975176] among A. urmiana and parthenogenetic populations (di-and triploids).

Complementary Experiments
To confirm the presence or absence of the 311 bp partial nucleotide in the mitochondrial sequence of A. urmiana, 48 further specimens have been examined following designed primers by Zhang et al. [26] to amplify this partial DNA region.
In order to consider sequence quality of the previously published mitogenome by Zhang et al. [26], the standard nucleotide BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was performed to check similarity of 16S rRNA sequence of NC_021382 [JQ975176] among A. urmiana and parthenogenetic populations (di-and triploids).

Transfer RNAs
As expected, the mitogenome of A. urmiana included all of the 22 tRNAs that are commonly found in metazoan mtDNA. Except for tRNA-Ser 1 , which shows a D-loop structure, all tRNAs possess the typical cloverleaf structure. The highest and lowest values of AT percentage in tRNAs were detected in tRNA-Glu (80.3%) and tRNA-Ser 1 (47.7%), respectively. The shortest and longest tRNAs were tRNA-Ala, tRNA-Ley (59 bp), and tRNA-Ser 2 (67 bp), respectively ( Figure 3).

Transfer RNAs
As expected, the mitogenome of A. urmiana included all of the 22 tRNAs that are commonly found in metazoan mtDNA. Except for tRNA-Ser1, which shows a D-loop structure, all tRNAs possess the typical cloverleaf structure. The highest and lowest values of AT percentage in tRNAs were detected in tRNA-Glu (80.3%) and tRNA-Ser1 (47.7%), respectively. The shortest and longest tRNAs were tRNA-Ala, tRNA-Ley (59 bp), and tRNA-Ser2 (67 bp), respectively ( Figure 3).
The relative synonymous codon usage (RSCU) results are summarized in Table S4 and  Table S4). (ND3 and ND6) and incomplete codons T (COX1, COX2 and ND5, and ND4) (Table S2). In A. urmiana, the AT-and GC-skews were negative in all encoded genes, except in ND5 and ND1, showed a positive GC-skew value ( Figure 2 and Table S3). The relative synonymous codon usage (RSCU) results are summarized in Table S4 and

Ribosomal RNAs
The 16S ribosomal RNA and 12S ribosomal RNA are positioned between tRNA-Leu2 and control region, where they are separated by tRNA-Val. Both ribosomal RNA genes are encoded on the light strand from 12,099 to 13,255 (1157 bp) and from 13,317 to 14,027 (711 bp), respectively ( Figure 1 and Table S2). The 16S ribosomal gene has a rather higher A+T content than 12S (63.96% vs. 60.90%). Both 16S and 12S ribosomal RNAs show positive values for AT-and GC-skew ( Figure 2 and Table S3).

Correction of the Previously Published A. urmiana Mitogenomic Structure
The current assembled mitochondrial genome of A. urmiana (MN240408) has revealed a shorter length than one that has been deposited in GenBank under NC_021382 [JQ975176] accession number (15,699 bp vs. 15,945 bp). There were high sequence similarities between both mitochondrial genomes (98%) except for a length difference of 311 bp difference. The 16S rRNA sequence of Zhang et al. [26] represented high similarity with parthenogenetic populations (di-and triploids) compared to A. urmiana (96.76-96.21% vs. 95.63-96.07%).
The length difference concerned the sequence of the control region (1672 bp vs. 1932 bp). Our analysis showed that the partial DNA sequence, including 12S ribosomal RNA, and the control region (311 bp length), extending from 13,884 bp to 14,194 bp, consisting of 130 bp from 12S rRNA and 181 bp from the control gene, were duplicated in the control region of the previously published GenBank sequence. The corresponding duplication was from 14,196 bp to 14,506 bp ( Figure 6). This part could not be identified in the current study. In addition, the 311 bp duplication (14,196 bp to 14,506 bp) could not be amplified in 48 specimens of A. urmiana using designed primers by Zhang et al. [26].

Ribosomal RNAs
The 16S ribosomal RNA and 12S ribosomal RNA are positioned between tRNA-Leu 2 and control region, where they are separated by tRNA-Val. Both ribosomal RNA genes are encoded on the light strand from 12,099 to 13,255 (1157 bp) and from 13,317 to 14,027 (711 bp), respectively ( Figure 1 and Table S2). The 16S ribosomal gene has a rather higher A+T content than 12S (63.96% vs. 60.90%). Both 16S and 12S ribosomal RNAs show positive values for AT-and GC-skew ( Figure 2 and Table S3).

Correction of the Previously Published A. urmiana Mitogenomic Structure
The current assembled mitochondrial genome of A. urmiana (MN240408) has revealed a shorter length than one that has been deposited in GenBank under NC_021382 [JQ975176] accession number (15,699 bp vs. 15,945 bp). There were high sequence similarities between both mitochondrial genomes (98%) except for a length difference of 311 bp difference. The 16S rRNA sequence of Zhang et al. [26] represented high similarity with parthenogenetic populations (di-and triploids) compared to A. urmiana (96.76-96.21% vs. 95.63-96.07%).
The length difference concerned the sequence of the control region (1672 bp vs. 1932 bp). Our analysis showed that the partial DNA sequence, including 12S ribosomal RNA, and the control region (311 bp length), extending from 13,884 bp to 14,194 bp, consisting of 130 bp from 12S rRNA and 181 bp from the control gene, were duplicated in the control region of the previously published GenBank sequence. The corresponding duplication was from 14,196 bp to 14,506 bp ( Figure 6). This part could not be identified in the current study. In addition, the 311 bp duplication (14,196 bp to 14,506 bp) could not be amplified in 48 specimens of A. urmiana using designed primers by Zhang et al. [26].
In the GenBank NC_021382 [JQ975176], tRNA-Phe has been mistakenly reported in the heavy strand but it is positioned in the light strand (see Table S2 and Figure 1). This mistake also was observed in annotation of mitogenome of A. franciscana (X69067) and A. tibetiana (JQ975177; NC_021383). Additionally, there was a major problem in position of ND5, which included the tRNA-Phe and would be biologically incorrect. Furthermore, there were several differences in 15 tRNA sequences between NC_021382 [JQ975176] and the current assembly (Figure 7). Our sequence has revealed GTG as a start codon for ATP6 whilst ATG was reported in NC_021382 [JQ975176]. TAA and TAG have been incorrectly reported as stop codons where these were belonging to COX1 and COX2 genes by Zhang et al. [26], while in NC_021382 [JQ975176] and MN240408 assemblies, the incomplete codons T have been revealed in COX1 and COX2 (Table S2). In the GenBank NC_021382 [JQ975176], tRNA-Phe has been mistakenly reported in the heavy strand but it is positioned in the light strand (see Table S2 and Figure 1). This mistake also was observed in annotation of mitogenome of A. franciscana (X69067) and A. tibetiana (JQ975177; NC_021383). Additionally, there was a major problem in position of ND5, which included the tRNA-Phe and would be biologically incorrect. Furthermore, there were several differences in 15 tRNA sequences between NC_021382 [JQ975176] and the current assembly (Figure 7). Our sequence has revealed GTG as a start codon for ATP6 whilst ATG was reported in NC_021382 [JQ975176]. TAA and TAG have been incorrectly reported as stop codons where these were belonging to COX1 and COX2 genes by Zhang et al. [26], while in NC_021382 [JQ975176] and MN240408 assemblies, the incomplete codons T have been revealed in COX1 and COX2 (Table S2).

Figure 6.
The partial map (a) and sequence alignment (b) of NC_021382 by Zhang et al. [26] showing duplicated part in control region. (NC_021382 has been corrected by NCBI staff. The reference sequence is identical to JQ975176).   (up) and previous mitochondrial assembly NC_021382 by Zhang et al. [26] (down). Nucleotides in red color display differences between two sequences (NC_021382 has been corrected by NCBI staff. The reference sequence is identical to JQ975176).

Artemia Mitogenome Variation
The complete mitochondrial sequences of four Artemia species had 10,709 conserved and 5047 variable sites, of which 1226 sites were parsimony informative; 3780 sites were singletons; 335 sites were gap. Among the mitogenomes sequenced so far, A. tibetiana (TIB) and A. franciscana (FRA) had the longest (15,826 bp and 15,822 bp) and A. sinica (SIN) and A. urmiana (URM) had the shortest length (15, 689 bp and 15,699 bp), respectively. Even so, regardless of the CR region, no significant difference could be observed among the concatenated tRNAs, PCGs and rRNAs sequences of the four recognized species of Artemia (URM: 14,027 bp, TIB: 14,014 bp, SIN: 14,027 bp, FRA: 14,000 bp).
The AT% content of tRNAs among Artemia mitogenomes is summarized in Table 1. Ten tRNAs reveal indicative differences of AT% (>5%), in which the highest dissimilarities were found between tRNA-Ser 1 (URM/TIB-FRA: 8.8%) and tRNA-His (URM/TIB-SIN: 8.5%), respectively. Nucleotide composition analysis indicates that the AT content of whole mitogenomes ranged from 62.54% (URM) to 64.51% (SIN). The CYTB gene exhibits the lowest AT content in all mitogenomes and varies between 58.18% (TIB) and 60.38% (FRA). The highest AT content was determined in ND3 gene of A. sinica (70.83%). In addition, the ND4L gene displayed a high AT content in other mitogenomes, ranging from 66.28% (TIB) to 69.77% (FRA) (Figure 8 and Table S5).
Nucleotide composition analysis indicates that the AT content of whole mitogenomes ranged from 62.54% (URM) to 64.51% (SIN). The CYTB gene exhibits the lowest AT content in all mitogenomes and varies between 58.18% (TIB) and 60.38% (FRA). The highest AT content was determined in ND3 gene of A. sinica (70.83%). In addition, the ND4L gene displayed a high AT content in other mitogenomes, ranging from 66.28% (TIB) to 69.77% (FRA) (Figure 8 and Table S5).  Detailed information of AT-and GC-skews in each individual protein-condoning gene and the ribosomal ones are provided in Table S3 and Figure 2. It was determined that AT-and GC-skews were positive for both ribosomal RNA genes. Generally, nucleotide skewness was positive with some exceptions. With the exception of FRA, the GC-skews of ND5 were positive for all species, with TIB (0.087) presenting the highest value. In COX1 and ND5, the GC-skew of FRA was reported zero with equal content of G and C. There were two obvious heterogeneity in the value of GC-skew in COX3 and ND1, where TIB (0.013)/FRA (0.096) and URM (0.035)/SIN (0.003) exhibited positive skewness, respectively. The distributional pattern of Artemia species based on AT-and GC-skew values of CmtG and PCG+rRNA sequences were shown on scatter plots (Figure 9). The results have documented the apparent dissimilarity of A. franciscana from other taxa. Detailed information of AT-and GC-skews in each individual protein-condoning gene and the ribosomal ones are provided in Table S3 and Figure 2. It was determined that AT-and GC-skews were positive for both ribosomal RNA genes. Generally, nucleotide skewness was positive with some exceptions. With the exception of FRA, the GC-skews of ND5 were positive for all species, with TIB (0.087) presenting the highest value. In COX1 and ND5, the GC-skew of FRA was reported zero with equal content of G and C. There were two obvious heterogeneity in the value of GC-skew in COX3 and ND1, where TIB (0.013)/FRA (0.096) and URM (0.035)/SIN (0.003) exhibited positive skewness, respectively. The distributional pattern of Artemia species based on AT-and GC-skew values of CmtG and PCG+rRNA sequences were shown on scatter plots (Figure 9). The results have documented the apparent dissimilarity of A. franciscana from other taxa.  Table S4 and Figures 4 and 5.

Genetic Distance and Phylogenetic Analysis
The genetic distances among Artemia mitogenomes (PCG+rRNA genes) are summarized in Table Altogether, the most frequently used codons were UUU and AUU. AGG codon was observed with the lowest frequency (URM: 0.2%; TIB: 0.17%; SIN: 0.17%; FRA: 0.2). The summary of RSCU analysis and percentage of amino acid in each species are documented in Table S4 and Figures 4 and 5.

Genetic Distance and Phylogenetic Analysis
The genetic distances among Artemia mitogenomes (PCG+rRNA genes) are summarized in Table 2 in which the highest and lowest distances were determined between A. franciscana-A. tibetiana/A. urmiana (0.214/0.123) and A. tibetiana-A. urmiana (0.090). Phylogenetic analyses (ML and BI) using concatenated PCGs and rRNAs mitogenome sequences have generated a concordant trees topology. According to the phylogenetic analysis, Artemia franciscana was placed as a clade, sister to the Asian spp. (Figure 10).

Discussion
In the present study, the complete mitochondrial genome of a taxonomically confirmed A. urmiana specimen from Urmia Lake was re-sequenced and compared with previous mitochondrial assembly (NC_021382 [JQ975176]) published by Zhang et al. [26]. We detected a difference in the sequence length of NC_021382 [JQ975176]: A partial sequence of 311 bp length (consisting of the 130 bp from 12S gene and 181 bp from the CR) is duplicated in the control region. This part covers more than 16% of the whole control region of NC_021382 [JQ975176] and would make it significantly longer than in other Artemia mitogenomes. In addition, we discovered several mistakes in characterization of mitogenome annotation. The absence of 311 bp duplicated partial in 48 specimens documents that this partial cannot be referred to an intraspecific variation in A. urmiana population.  [50,51]. Our finding shows that the control region of the Artemia mitogenome can differ by 10% between species (A. urmiana: 1672 bp vs. A. franciscana: 1822 bp). Therefore, it is suggested that the control region is less conserved in the genus Artemia. The control region has a particular function in initiation and regulation of DNA transcription and replication [52]. Although mutations in this part may not

Discussion
In the present study, the complete mitochondrial genome of a taxonomically confirmed A. urmiana specimen from Urmia Lake was re-sequenced and compared with previous mitochondrial assembly (NC_021382 [JQ975176]) published by Zhang et al. [26]. We detected a difference in the sequence length of NC_021382 [JQ975176]: A partial sequence of 311 bp length (consisting of the 130 bp from 12S gene and 181 bp from the CR) is duplicated in the control region. This part covers more than 16% of the whole control region of NC_021382 [JQ975176] and would make it significantly longer than in other Artemia mitogenomes. In addition, we discovered several mistakes in characterization of mitogenome annotation. The absence of 311 bp duplicated partial in 48 specimens documents that this partial cannot be referred to an intraspecific variation in A. urmiana population. Although  finding shows that the control region of the Artemia mitogenome can differ by 10% between species (A. urmiana: 1672 bp vs. A. franciscana: 1822 bp). Therefore, it is suggested that the control region is less conserved in the genus Artemia. The control region has a particular function in initiation and regulation of DNA transcription and replication [52]. Although mutations in this part may not produce faulty mitogenome products, its variability may change the rate of expression of mitochondrial genes [26,53] The mitochondrial gene contents of Artemia are uniform with the ancestral Pancrustacea model, exhibiting 22 tRNAs, 13 PCGs and two rRNAs genes (1). Contrary to the dissimilar gene arrangements in mitochondrial genomes among species of decapods [54], as well as the branchiopod Daphnia [55], the gene order recognized for Asian Artemia species was identical to the previously sequenced American A. franciscana. This finding suggests that the gene arrangements in the mitochondrial genome of Artemia from the Old World to New World were highly conserved (see [56]).
As with other branchiopods [55], ATG was found as the start codon in most proteincoding genes in A. urmiana and two other species, A. tibetiana (n = 8) and A. franciscana (n = 5), have displayed the highest and lowest abundances, respectively [26,57]. Most PCGs terminate with a TAA codon in Artemia. We recognize that the COX1, COX2, ND5 and ND4 genes terminate with an incomplete T stop codon (T-or TA-) in Asian species. Additionally, COX3 terminates with an incomplete stop codon in A. sinica. However, no incomplete stop codon was found in the mitogenome of A. franciscana [57].
In this study, the AT content of the A. urmiana mitogenome was calculated, and subsequently compared with the other species. Similar to decapod mitogenomes that are generally AT-rich [54], our results show that the AT content was substantial in Artemia. Though complete mitogenomes of Artemia species had identical AT content ( Figure 8 and Table S5), there considerable differences in AT content (>5%) were found in ATP8, ATP6, ND3 and ND6, and almost half of the tRNAs (Table 1) among species. Additionally, the ATand CG-skew of the mitogenome provide useful information respecting the strand-specific nucleotide frequency bias [42,58,59]. COX3 and ND1 genes exhibited more asymmetries of AT-and GC-skews among the available Artemia mitogenomes.
The mitochondrial markers, COX1 and 16S, have been successfully used in phylogeny of branchiopods [60][61][62][63][64][65]. To date only COX1 has been utilized for phylogenetic studies on Artemia, nevertheless mitogenomic results demonstrated significant difference in the nucleotide composition of ATP8, ATP6, ND3, ND6, ND1 and COX3. Due to the phylogenetic problems in the genus Artemia [21], other mitochondrial markers with high nucleotide variation should be considered in future investigations of intra-and interspecific variation of the genus Artemia. Further studies are required to characterize remaining complete Artemia mitogenomes of bisexual and parthenogenic lineages. This would clarify their phylogenetic relationships and pattern of diversification, as well as suggest new mitochondrial markers for evolutionary and population genetics studies.
Due to phylogenetic problems of the genus Artemia, there is an ongoing project related to sequencing of the complete mitochondrial genomes of all bisexual species and parthenogenetic Artemia with different ploidy levels by Hainan Tropical Ocean University (China) and the result of this study would clarify the phylogenetic relationships among Artemia taxa.