Using ISSR Genomic Fingerprinting to Study the Genetic Di ﬀ erentiation of Artemia Leach, 1819 (Crustacea: Anostraca) from Iran and Neighbor Regions with the Focus on the Invasive American Artemia franciscana

: Due to the rapid developments in the aquaculture industry, Artemia franciscana , originally an American species, has been introduced to Eurasia, Africa and Australia. In the present study, we used a partial sequence of the mitochondrial DNA Cytochrome Oxidase subunit I ( mt- DNA COI ) gene and genomic ﬁngerprinting by Inter-Simple Sequence Repeats (ISSRs) to determine the genetic variability and population structure of Artemia populations (indigenous and introduced) from 14 di ﬀ erent geographical locations in Western Asia. Based on the haplotype spanning network, Artemia urmiana has exhibited higher genetic variation than native parthenogenetic populations. Although urmiana represented a completely private haplotype distribution, no apparent genetic structure was recognized among the native parthenogenetic and invasive A. franciscana populations. Our ISSR ﬁndings have documented that despite that invasive populations have lower variation than the source population in Great Salt Lake (Utah, USA), they have signiﬁcantly revealed higher genetic variability compared to the native populations in Western Asia. According to the ISSR results, the native populations were not fully di ﬀ erentiated by the PCoA analysis, but the exotic A. franciscana populations were geographically divided into four genetic groups. We believe that during the colonization, invasive populations have experienced substantial genetic divergences, under new ecological conditions in the non-indigenous regions.


Introduction
The brine shrimp Artemia is a unique zooplankton that has a limited number of species, distributed globally, except in Antarctica [1]. This tiny crustacean has potentially adapted to live in extreme environmental conditions, such as hypersaline environments [2,3].
Artemia has been mainly used as live food in larviculture and fishery industries, especially in Asia [4]. Artemia has been used to improve the quality of sodium chloride production in solar salt-fields [5,6]. It was also introduced as a model organism in many bioscience studies, including

Sample Collection and DNA Extraction
Artemia cyst specimens were collected from 14 geographical localities across Iran and neighbor countries ( Figure 1). All studied populations had been confirmed to be bisexual or parthenogenetic according to Asem et al. [3]. The sample localities with their geographical coordinates, abbreviations and IPMB voucher are summarized in Table 1. Total genomic DNA was extracted from part of the antenna of adult males and females using the Chelex ® 100 Resin method (6%, Bio-Rad Laboratories, Hercules, CA, USA) [16]. All extracted DNA was stored at −80 • C for subsequent genetic characterization.  Table 1).

Population Identification and Phylogenetic Analyses
A partial sequence of the mitochondrial marker cytochrome c oxidase subunit I (COI) was utilized to identify the taxonomical status of the studied populations using phylogenetic analyses as implemented in the MEGA X program (Temple University, Philadelphia, USA) [2,17,18]. To identify the taxonomical status of the studied populations, the COI reference sequences from the recognized bisexual species and parthenogenetic populations were downloaded from GenBank (Table 2). Sequences were aligned using MEGA X with default settings [31]. Phylogenetic trees were reconstructed based on the Maximum Likelihood approach included in the MEGA X program. To reveal the genealogical and geographical relationships, a median haplotype network was established, following the median-joining algorithm in the Network program ver. 5.0.1.1 (Universität Hamburg, Hamburg, Germany) [32].

Genomic Fingerprinting by ISSR-PCR
Genomic variability was examined by inter simple sequence repeat ISSR-PCR using the same DNA template used for phylogenetic analyses. Initially, 15 ISSR primers were analyzed to distinguish the intra-and inter-specific genetic variability within and among 83 randomly selected individuals, belonging to 14 geographically different localities of Artemia. Out of 15 tested ISSR primers, five were selected because of unambiguous banding patterns of the PCR products (Table 3). . PCR amplifications were executed in a thermal cycler based on the following conditions: 94 • C denaturation for 1 min, 35 cycles of 46-54 • C annealing for 50 s and 72 • C extension for 2 min. The final cycle was continued for 7-min at 72 • C. Final PCR products were mixed with 8 µl bromophenol blue and run on high-resolution denaturing polyacrylamide gels 6% (0.2 mm) for 3 h at 65 W (size 45 × 30 cm) including 1× TBE buffer. The gels were dried and exposed for 2 days to X-ray hyperfilm (Kodak, Taufkirchen, Germany) and subsequently developed. Finally, the autoradiograms were scanned to identify the polymorphic bands [25].

ISSR Statistics
The quality and quantity of ISSR bands were inspected visually. Ambiguous and smeared bands were excluded from the analysis, and only unequivocally reproducible bands were scored for each individual as present (1) or absent (0). The binary data matrix (presence = 1; absence = 0) was formulated in MS Excel v.2016 and used for subsequent genetic analyses.
ISSR data were analyzed via the Bayesian model-based clustering algorithm as implemented in the STRUCTURE v. 2.3 program (University of Oxford, Oxford, United Kingdom) [34,35]. We analyzed the genetic structure among populations by assigning individuals into potential numbers of clusters (K = 1 − 10). ISSR genotypes were processed with a period of burn-in 50,000 and 100,000 MCMC repetitions [34]. The CLUMPAK online program was employed to identify the pattern of clustering modes and packaging population structure [35]. The online programs, CLUMPAK [36] and STRUCTURE HARVESTER [37] were implemented to assess and visualize the most appropriate number of K by calculating the likelihood of the posterior probability [38].
The binary data matrix was employed to calculate the genetic diversity parameters of each population using GenAlex ver. 6.5 (Australian National University, Acton, Australia) [39]. The population genetic parameters were as follows: Na (number of different alleles), Ne (number of effective alleles), I (Shannon's information index), He (expected heterozygosity), uHe (unbiased expected heterozygosity), PPL (percentage of polymorphic loci), NB (number of bands) and NPB (number of private bands) and pairwise population matrix of Nei genetic distance [28,29].
Intra-and inter-specific molecular variations and genetic relationships among populations were implemented by Principal Coordinate Analysis (PCoA) and Analysis of Molecular Variance (AMOVA) as utilized by GenAlex ver. 6.5, respectively [39].
To better understand the population genetic variations, ISSR analyses were performed on three platforms separately, as follows: whole populations, native populations and invasive American A. franciscana.

Phylogenetic Analyses and Haplotype Distribution
Our phylogenetic analyses provide evidence that the studied bisexual specimens from three localities of Iran, Nough Catchment (NOG), Mahshar port (MSH) and Maharlu Lake (MAHR), and a locality from Iraq, Garmat Ali (GAA), clustered in the clade of A. franciscana ( Figure 2). In addition, all parthenogenetic populations clustered in two separated clades that shared a common ancestor with urmiana. Although most of the parthenogenetic populations were located in clade P1, the majority of CAM specimens (eight out of nine). The clade P1 contained two sub-clades, consisting of diploids and triploid populations.  Table 1). Figure 3 represents the haplotype spanning network of COI among native A. urmiana and parthenogenetic populations. Results demonstrated that A. urmiana has wider genetic variation compared to parthenogenetic ones. Genetic differentiation and a close relationship of Camalti Lake (CAM) population (Turkey) with A. urmiana was clearly revealed in the tree. While COI sequences of nine parthenogenetic populations were distributed in five Haplotypes (H2-6), A. urmiana showed private haplotypes without a shared haplotype in other populations. The COI sequences of invasive A. franciscana were grouped into six distinct haplotypes. No population with private haplotypes was observed ( Figure 4). The numbers of individuals and population composition were calculated for each haplotype (Appendix A, Tables A1 and A2).  Table 1).  Table 1).

ISSR Profiling
We have shown previously that ISSR can detect substantial genetic variation in the genus Artemia [28,29]. ISSR fingerprints can differ within and among populations. Altogether, 152 observed and unambiguously identified ISSR bands were analyzed. The total number of polymorphic loci showed a mean value of 25.42 ± 2.28%; the lowest number was seen in Camalti Lake (CAM) (7.24%) and the highest in Nough Catchment (NOG) (38.82%). Generally, the lowest values of genetic indices belonged to the parthenogenetic CAM, while the highest values were shared between the native parthenogenetic population from Eastern lagoon around Urmia Lake (LAGE) (Ne = 1.256 ± 0.031, I = 0.208 ± 0.024, He = 0.143 ± 0.017, uHe = 0.171 ± 0.020) and invasive A. franciscana from NOG (Na = 1.013 ± 0.071, I = 0.208 ± 0.023) ( Table 4). In summary, 134 and 126 distinguished ISSR bands were examined for ten native and four American A. franciscana populations, respectively (Tables 5 and 6). According to the results of the Nei genetic matrix, parthenogenetic CAM from Turkey and invasive GAA from Iraq showed the high genetic distance with parthenogenetic and invasive populations, respectively (Tables 7 and 8).       Based on the total results of ISSR, the AMOVA analysis documented that most of the genetic variations were attributed among native and invasive populations more than within populations (69% vs. 31%). There was no indicative genetic variability observed in the midst of among-and within variation (55% vs. 45%) in native populations. In contrast, the high differentiation was represented within populations (71% vs. 29%) of non-indigenous A. franciscana in Asia (Table 9, Figure 5A-C).  Bayesian clustering analysis using STRUCTURE was performed to investigate the genetic patterns of the studied populations. The optimum K was obtained at K = 2 for the whole 14 populations and ten native populations, and K = 9 for four invasive A. franciscana, respectively. Figure 6 showed the clustering of genetic structures, where the first highest posterior probability (K) was represented by different colors for each population. With regard to the genetic patterns of ISSR, native and exotic populations could be completely divided into two groups ( Figure 6A). The results of the analysis for native populations documented that parthenogenetic CAM is a distinct population with a differing clustering pattern. The STRUCTURE analysis could not fully distinguish A. urmiana and parthenogenetic populations ( Figure 6B). The high value of optimum K (K = 9) was prominent in the clustering analysis for A. franciscana populations, which generally revealed a complex pattern ( Figure 6C). Proportions of genetic clusters (percentage) for each locality were summarized in Figure 7 (see also Tables A3-A5).  Table 1).  Table 1).
The first and second PCoA coordinates explained 64.66% and 7.00% of the variance, respectively (overall, 71.66% of total variation). The results showed all populations clustered into three groups, where invasive A. franciscana has been significantly separated from the native populations based on the first coordinate. Ten native populations (including A. urmiana and parthenogenetics) were divided into two distinct groups, G2 and G3. G2 included all of the CAM population and a single individual of Eastern lagoon around Urmia Lake (LAGE). Bisexual A. urmiana and other parthenogenetics were placed in G3 (Figure 8). The results of the separate analyses of PCoA for native and invasive populations are shown in Figures 9 and 10, which include overall 42.55% and 41.55% of the total variation, respectively. Although native populations produced almost the same result with the "whole populations" analysis ( Figure 9), the separated analyses of PCoA for invasive A. franciscana could separate all four populations in isolated groups ( Figure 10).  Table 1).  Table 1).  Table 1).
Similar to the phylogeny based on COI sequences, clustering analysis of ISSR sequences by STRUCTURE has revealed a distinguished structure for non-indigenous A. franciscana populations. In this analysis, all native populations (bisexual A. urmiana and parthenogenetic ones) have displayed almost similar patterns, but a separated analysis for native populations has revealed a non-identical structure for the CAM population. Additionally, an inconsistent pattern of the Eastern lagoon around Urmia Lake (LAGE) was also observed. A separated analysis for invasive populations could not reveal different patterns among examined populations by STRUCTURE ( Figure 6A-C). On the other hand, PCoA has divided A. franciscana from native populations. Although PCoAs could not branch off natives by localities, the separated PCoA has divided invasive populations based on localities in the four groups. Contrary to COI haplotype distribution, the ISSR marker was unable to reveal a private pattern for bisexual A. urmiana in comparison to native parthenogenetic populations.

Discussion
The present study was performed to compare the population structure and genetic differentiation of native and invasive Artemia populations. The mitochondrial COI gene has been established as a useful molecular marker to determine the intra-and inter-specific evolutionary associations [2,3,17,18]. Asem et al. [2] have documented that di-and triploid parthenogenetic brine shrimps are maternally related to A. urmiana, while tetra-and pentaploid lineages shared a common maternal ancestor with A. sinica. Based on the mitochondrial COI dataset, all examined parthenogenetic individuals have grouped in a close evolutionary relationship with A. urmiana (Figure 2). Our results have demonstrated that they should include di-and/or triploids. This observation has also been confirmed by the phylogenetic tree. Although individuals of eight Artemia populations have been exactly located in sub-clades of diand/or triploid, eight out of nine specimens of CAM from Turkey have been placed in a particular clade (P2), in a close connection with A. urmiana ( Figure 2). This finding has also been confirmed by haplotype distribution (Figure 3). Previously, Sayg [40] has confirmed that triploid and pentaploid parthenogenetic populations coexisted in Camalti Lake (CAM). Our observation has confirmed that Camalti Lake (CAM) populations had a parthenogenetic reproductive mode and those examined specimens should be considered as a diploid and/or triploid population (see [3]). Despite the fact that A. urmiana shared a common ancestor with di-and triploids, it has presented an unexpectedly high level of haplotype diversity of the COI marker. These results have also been documented by Eimanifar and Wink [28]. The high level of haplotype variation might be attributed to the evolutionary life history of A. urmiana [30] and/or its large population size [28]. Urmia Lake has undergone considerable changes in environmental conditions, such as salinity and temperature [41,42], which could have influenced genetic variation and population size during evolution.
Although COI sequences should reveal a phylogeographic structure in the closely related species [43], our results could not determine a level of geographical differentiation among native as well as among invasive populations. Contrary to the results of the mitochondrial marker, genomic fingerprinting ISSR could not reveal a significant high level of genetic variation in A. urmiana. This result might be due to differences in the rate of variation of mitochondrial and nuclear genes and potential hybridization events in the past.
The ISSR method had been utilized to study ten diploid parthenogenetic Artemia populations from China by Hou et al. [27]. Their genetic variation was significantly higher than our examined parthenogenetic populations. For example, the percentage of polymorphic loci ranged from 54.12-87.06% vs. 8.21-39.55%. Western Asia is the origin of di-and triploid parthenogenetic populations [3], so high genetic diversity of Chinese populations could be the result of their adaptation during colonization through biological dispersal to the new environments in East Asia under historical evolutionary progress.
In general, it was reported that the invasive populations have lower genetic variation in the new environments as compared with their origin populations [44]. The colonized population of A. franciscana in Vietnam has a lower genetic diversity than its native source population from San Francisco Bay (SFB) [22]. In contrast, Eimanifar et al. [2] have documented that Asian invasive A. franciscana populations had higher genetic variation than the American Great Salt Lake (GSL) population and native Asian species. Similar results have been reported for some invasive populations from Mediterranean regions [23,24]. A recent study assumed that invasive populations of A. franciscana show a wide degree of genetic differentiation in Australia [17].
Overall, our ISSR results have documented that invasive A. franciscana populations had distinctly higher genetic variation than Western Asian native parthenogenetic populations. On the other hand, native A. franciscana from Great Salt Lake (GSL) have represented higher variation than examined invasive populations in this study, as the percentage of polymorphic loci differed from 67-81% vs. 30.95-46.83% (see [29]). Additionally, all four invasive A. franciscana populations clearly revealed a different genetic structure. Observation of low genetic diversity in native populations might be attributed to the effect of asexual reproduction in parthenogenetic populations and/or critical climatic conditions in West Asia, especially Urmia Lake in the last two decades (see [3,30]). We believe that interactions between different ecological conditions in the new environments and the high potential of physiological plasticity and genetic adaptation of A. franciscana could exert different evolutionary pathways during the introduction of exotic populations, which would have ultimately caused intra-specific variations and genetic divergence in the examined invasive populations.
In conclusion, it is expected that the non-indigenous species should have a lower genetic variation than their source populations [25,44,45]. However, non-indigenous A. franciscana populations gave opposite results in comparison with native populations from GSL and SFB. Since there is neither a taxonomical identification key nor morphological identifications to distinguish bisexual species and parthenogenetic populations [11,46], it would not be possible to identify the exotic population at the earliest time of invasion. Therefore, there is a lack of information to regularly determine the colonization progress and evolutionary development of A. franciscana in the new habitats. We assume that differences in the genetic variation of non-indigenous populations could be due to the study on different invasion periods consisting of i) introduced, ii) establishing/colonizing, iii) established/colonized populations.
Author Contributions: Research design, material preparation, and data collection were performed by A.E. Data analysis was carried out by A.A. The first draft of the manuscript was written by A.A. and A.E. P.-Z.W., W.L. and M.W. reviewed the draft. All authors have read and agreed to the published version of the manuscript.
Funding: This study was carried out at IPMB, Department of Biology, University Heidelberg and A/10/97179. Amin Eimanifar was supported by a Ph.D. fellowship from the Deutscher Akademischer Austauschdienst (DAAD, German Academic Exchange Service).