The impact of one-decade ecological disturbance on genetic changes: a study on the brine shrimp Artemia urmiana from Urmia Lake, Iran

Urmia Lake, the largest natural habitat of the brine shrimp Artemia urmiana, has progressively desiccated over the last two decades, resulting in a loss of 80% of its surface area and producing thousands of hectares of arid salty land. This ecological crisis has seriously affected the lake’s native biodiversity. Artemia urmiana has lost more than 90% of its population during the decade from 1994 (rainy period) to 2004 (drought period) due to salinity increasing to saturation levels (∼300 g/l). We studied the influence of this ecological crisis on the genetic diversity of A. urmiana in Urmia Lake, based on one cyst collections in 1994 and 2004. AMOVA analysis on ISSR data demonstrated a 21% genetic variation and there was a 5.5% reduction of polymorphic loci between samples. PCoA showed that 77.42% and 68.75% of specimens clustered separately in 1994 and 2004, respectively. Our analyses of four marker genes revealed different genetic diversity patterns with a decrease of diversity at ITS1 and an increase for Na+/K+ ATPase. There was no notable difference in genetic variation detected for COI and 16S genes between the two periods. However, they represented distinctly different haplotypes. ITS1 and COI followed a population expansion model, whereas Na+/K+ ATPase and 16S were under demographic equilibrium without selective pressure in the 1994 samples. Neutrality tests confirmed the excess of rare historical and recent mutations present in COI and ITS1 in both samples. It is evident that a short-term ecological disturbance has impacted the genetic diversity and structure of A. urmiana.


INTRODUCTION
Urmia Lake (37 • 42 N, 45 • 19 E) is a landlocked thalassohaline lake with oligotrophic characteristics located in Northwest Iran. Its historical water surface area has ranged from 4,750 to 6,100 km 2 with the average and greatest recorded depths being 6 and 16 m, respectively (Azari Takami, 1993;Van Stappen, Fayazi & Sorgeloos, 2001). It is among the largest hypersaline lakes in the world, like Great Salt Lake, USA, which has an average surface area from 4,400 to 8,500 km 2 (Abatzopoulos et al., 2006;USGS, 2013), and it is inhabited by the brine shrimp Artemia urmiana Günther, 1899. Contrary to widespread opinions, Urmia Lake and its adjacent wetlands are the habitat of various organisms. Based on its unique biodiversity, environmental gradients, socioeconomic importance and existence of indigenous communities, Urmia Lake has been registered as a protected area since 1967 and as a national park since 1975. Because of its importance for migratory birds, it was also registered in the Ramsar Convention on Wetlands as a wetland of international importance in 1975 and considered as one of the 59 biosphere reserves by UNESCO in 1976(Eimanifar & Mohebbi, 2007Asem et al., 2014;Asem, Eimanifar & Sun, 2016;Asem, Eimanifar & Wink, 2016).
Historical records document that Urmia Lake has grappled with drought crises over centuries to the extent that locals were able to walk across the lake (about 20 km) via a paved road (Morier, 1818;Curzon, 1892). Additionally, it was reported that due to the lack of freshwater and food, herbivorous animals, inhabiting the islands within the lake, deserted the islands by swimming and migrated into the surrounding mountains (Binder, 1887).
Urmia Lake had the highest water-level elevation in 1994-1996 (1277.8 m a.s.l.) over the past six decades (from 1955 to 2015). Based on the first resource assessment of Artemia cysts and biomass in 1994-1995, Artemia urmiana cyst production in the upper 50 cm of the lake's water column ranged from 4,200 to 4,500 tonnes/year (dry weight) (Sorgeloos, 1997). Analysis estimated that the cyst concentration was 399 cysts/l in that period (Asem, Mohebbi & Ahmadi, 2012). The water level of the lake fell below the ''minimum ecological water level '' after 2001(1274Abbaspour & Nazaridoust, 2007) (Fig. 1). Later estimates of Artemia cyst production declined to 27 and 25 cysts/l in 2003 and 2004, respectively, when salinity increased to saturated levels (∼300 g/l) (Ahmadi, 2005). In the following years cyst production dropped from an estimated 11 cysts/l in 2005 to 3 cysts/l in 2007(Ahmadi, 2007. The lake lost most of its area after 2007 and no further assessments were performed, but some estimations indicated that the cyst concentration decreased to below 1 cyst/l (Asem, Mohebbi & Ahmadi, 2012). No live Artemia were observed in the main body of the lake during the summer of 2016, but did occur in the surrounding lagoons and estuaries (IARC, 2016). Now Urmia Lake has become an ecological disaster, receiving international attention. Iran's Department of the Environment and the United Nations Development Programme (UNDP) ratified a project to save the lake and the surrounding wetlands (UNDP, 2014), at an estimated restoration cost of $1.3 billion (Hecht, 2014). Critical environmental conditions can affect biodiversity and species distribution (Menendez et al., 2006;Barrett & Schluter, 2008;Jump, Marchant & Penuelas, 2009;Berkhout et al., 2014). Genetic diversity plays a decisive role in evolutionary history and future evolutionary directions of taxa (May, 1994;Forest et al., 2007;Jump, Marchant & Penuelas, 2009). While there have been many studies that focused on the effect of environmental changes on biodiversity, few studies have focused on intraspecific genetic variation (Pauls et al., 2013). Li et al. (1999), Li et al. (2000) and Li et al. (2001) documented significant genetic differentiation of wild emmer wheat (Triticum dicoccoides) in response to ecological change. A similar pattern was observed in slender oat (Avena barbata), followed by environmental variations (Jump, Marchant & Penuelas, 2009). A comparable pattern of genetic variation was observed in European white birch populations (Betula pendula), which showed different genotypes in warm and cool years (Kelly et al., 2003). An additional example was the observed significant negative correlation between Gly-3 allele frequency and increasing summer precipitation in pinon pine (Pinus edulis) (Mitton & Duran, 2004).
Most studies on rapid evolutionary responses focus on morphological, physiological and nutritional variation. There are few studies that consider altered environmental conditions, especially short-term crises, on genetic variability (Thompson, 2009). Based on cyst collections and assessments in 1994 and 2004, Artemia has lost more than 90% of its reproductive potential and population size, and reproduction has stopped in the main body of Urmia Lake (Asem, Mohebbi & Ahmadi, 2012). We hypothesized that the tremendous changes in the environmental conditions of Urmia Lake and the reduction of the Artemia population size may have affected the genetic diversity. And, if the short-term ecological disturbance has influence on genetic variation, this change would be able to document in the intraspecific genomic dissimilarity. The specific objective of this study was to determine the independent impressionability of mitochondrial and nuclear genes through ecological crises.

Sampling strategy and DNA extraction
To assess intraspecific variation and population structure, we examined quiescent Artemia embryos collected in 1994 and in 2004 from Urmia Lake (Kholman-khaneh station; 45 • 29 E, 37 • 64 N). The samples were obtained from the upper 50 cm of the water column (see Sorgeloos, 1997). Because the bisexual A. urmiana coexists with a low ratio of a parthenogenetic population (Azari Takami, 1989;Agh et al., 2007) every specimen studied was first identified as belonging to the bisexual A. urmiana using SNP polymorphism in the Na + /K + ATPase α-1 subunit (Manaffar et al., 2011), and re-certified by the phylogenetic analyses using the COI mitochondrial marker.
Total DNA of each specimen decapsulated embryo (number of specimens for each experiment below) were extracted following the Chelex R 100 Resin method (Bio-Rad Laboratories, Hercules, CA, USA). The embryos were crushed via a sterilized pipette tip, incubated for 2.5-3 h at 60 • C (tubes were shaken by vortex every 30 min) and eventually for 10 min at 80 • C. The tubes were centrifuged at 10,000 rpm for 1 min and the supernatant phase was directly used in the PCR reaction (Montero-Pau, Gómez & Muñoz, 2008;Eimanifar & Wink, 2013;Asem, Eimanifar & Sun, 2016;Asem, Eimanifar & Wink, 2016). The extracted DNA was stored at −80 • C for further genetic analyses.

ISSR statistics
The binary matrix (1 = presence; 0 = absence of a band) was determined for each year and population genetic information was computed separately. Genetic relationships among ISSR genotypes were established by principal coordinate analysis (PCoA) using GenAlex version 6.5 (Peakall & Smouse, 2012). The partition of genetic variation within and between 1994 and 2004 was determined using the Analysis of Molecular Variance implemented in GenAlex ver. 6.5 program (Peakall & Smouse, 2012).
A partial fragment of the nuclear gene, Na + /K + ATPase α-1 subunit, was amplified using the primers of Manaffar et al. (2011). PCR amplification was performed under the following conditions: 94 • C for 2 min, 32 cycles of 94 • C for 25 s and 56 • C for 25 s and 72 • C for 1 min, and final extension with 72 • C for 3 min.
A fragment of the nuclear DNA containing a partial sequence of the 18S ribosomal RNA (18S), the complete sequence of internal transcribed spacer 1 (ITS1) and a partial sequence of the 5.8S ribosomal RNA (5.8S) genes, was PCR-amplified using the primers 18d-5 /R58 (Baxevanis, Kappas & Abatzopoulos, 2006). The thermal cycler PCR conditions were as follows: 4 min at 93 • C, 32 cycles of 40 s at 93 • C, 40 s at 62 • C, 1 min at 72 • C, and a final extension of 5 min at 72 • C.
Amplification of a partial fragment of the mitochondrial cytochrome oxidase subunit 1 (COI ) gene was performed using the invertebrate universal primers LCOI 490/HC02198 (Folmer et al., 1994). PCR amplification was carried out using the following program: a cycle of 3 min at 95 • C, followed by 35 cycles of one min at 95 • C, one min at 40 • C and one and half min at 72 • C, with a final step of 7 min at 72 • C.
The fragment of 16S ribosomal RNA (16S) was amplified using the primers 16S-SP/12S-SP (Bossier et al., 2004). PCR amplification was carried out under the following conditions: 1 cycle of 94 • C for 2 min, 34 cycles of 1 min 15 s at 94 • C, 45 s at 52 • C, 2 min at 72 • C and a final extension cycle of 72 • C for 4 min.
Our DNA dataset consisted of 248 sequences including 70 specimens sampled for Na + /K + ATPase, 60 specimens for ITS1 and COI, 58 specimens for 16S genes. The list of genetic markers and GenBank accession numbers is presented in Table 1.

Sequence alignment and population genetic diversity
Sequences were aligned using MEGA ver. 6.00 with MUSCLE tool and default parameters (Tamura et al., 2013). Alignment lengths were 198, 1150, 665 and 875 bp for Na + /K + ATPase, ITS1, COI and 16S, respectively. Between-group mean distances (year 1994/2004) were computed using p-distance in MEGA ver. 6.00. To estimate the genealogical relationships among haplotypes for each gene, a maximum-parsimony haplotype network was inferred using the software TCS version 1.21 (Clement, Posada & Crandall, 2000).

Species identification
The specimens analyzed in this study had a homozygous pattern (T-T) in the last valine codon using the Na + /K + ATPase α-1 subunit with the exception of a single specimen collected in 1994 which showed a heterozygous pattern (T-G). Our phylogenetic trees (ML and BI) for COI showed that all analyzed specimens clustered with the reference sequence of A. urmiana (Maniatsi et al., 2011: HM998991) (Fig. S1). Only one specimen from 1994 placed in the diploid parthenogenetic clade also revealed a heterozygous pattern in valine codon; this sample was removed from the dataset.

ISSR Profiling
A summary of population genetic indices of A. urmiana for all observed ISSRs loci between 1994 (rainy period) and 2004 (drought period) is listed in Table 2. ISSR profiling generated  Table 3). The first and second PCoA coordinates contained 18.94% and 13.80% of the variance, respectively (overall 32.74% of total variation). PCoA demonstrated that 1994 and 2004 were distinct groups, since there was only a narrow overlap between them (Fig. 2), with 77.42% and 68.75% specimens from the rainy and dry year being distinguished, respectively (Table 4).
The COI sequences for 60 specimens contained 43 haplotypes (H1-H43). Haplotype 1 was the major haplotype that was found in ten (16.7%) specimens. H12 and H23 were found in two and four specimens from 2004, respectively. H31 and H43 were shared by two specimens belonging to the two sampling years, while all other haplotypes came from a single sampling year (Fig. 5).
The 16S sequences for 58 specimens revealed 56 haplotypes (H1-H56), which showed high variation in comparison with the other markers. The central haplotype (H1) covered  only two specimens (3.5%) including a single specimen from each year. With the exception of H39 that was shared by two specimens from 2004, other haplotypes were unique to a single specimen (Fig. 6).

Genetic variation and neutrality tests
Genetic indices and allele frequency estimated for the aforementioned four markers are presented in Table 1. of polymorphic sites, total number of mutations and number of haplotypes of the 16S marker showed lower values in 1994. Similar to COI there was no notable dissimilarity for the haplotype diversity, nucleotide diversity and average number of nucleotide differences in the 16S marker between the two periods. Haplotype frequencies of all mitochondrial and nuclear markers represented non-significant difference between rainy and drought periods. Though there was no significant difference in the amount of genetic variation of COI and 16S markers between rainy and drought periods, each marker presented different distinct of haplotypes. 16S presented higher expected heterozygosity in 1994 while ITS1 and COI showed higher values in 2004. The values of the pairwise genetic differentiation index (F ST ) were not significant. The minimum and maximum between-group distances were detected in Na + /K + ATPase (0.001) and 16S (0.011), respectively. Neutrality tests yielded negative values with different significant and non-significant levels.

DISCUSSION
The effect of ecological disturbance, especially short-term regional climate changes, on genetic diversity is not well understood (Bálint et al., 2011;Banks et al., 2013). Ruediger et al. (2012) demonstrated that inter-annual and seasonal changes in water temperature produced significant variation in the genetic structure of Daphnia populations. The impact of salinity changes on genetic variation has been considered less frequently in aquatic organisms (Stoks, Geerts & De Meester, 2014). Genetic variation and genetic diversity are important parameters to conserve biodiversity at all levels including population variabilities, individual fitness and adaptability of species to the environmental conditions (Amos et al., 2001;Hughes et al., 2008). Recently, Avolio, Beaulieu & Smith (2013) showed that the genotypic diversity of Andropogon gerardii (big bluestem grass) was significantly reduced after 10 years of an increase of experimentally-driven intra-annual precipitation variation. Brown et al. (2013) suggested habitat conflagration as a major critical process to reduce allelic richness of the mallee emu-wren Stipiturus mallee by reducing population size. Similar studies in zooplankton species such as Artemia have not been done. Environmental instability directly affects the F ST (genetic differentiation among populations) through its impact on immigration and genetic drift combined with population reduction (Banks et al., 2013). Genetic diversity has been indicated to be important to population fitness since low levels of genetic variation may decrease the ability of population to adapt to the environmental crisis (Chapman et al., 2009;Pauls et al., 2013). Although Urmia Lake is a wetland of international importance and is facing an acute ecological threat, few studies have assessed risks to its biodiversity. Asem et al. (2010) reported that in a rainy period (1994) the Artemia of this lake had a higher cysts size variation, significantly larger average egg size, and a thinner chorion than in a dry period (2004). The smaller cysts and thicker chorion produced during the dry period were attributed to decreasing food availability and to an acclimation mechanism, respectively, to increase the survivorship of the diapausing embryo under ecological crisis (Asem et al., . Sankian, Heydari & Manaffar (2011) showed that A. urmiana hatching from cysts collected in 1998 (salinity = 180 g/l) had lower mortality but higher RNA content than those from 2003 (salinity approximately 300 g/l; saturated). Our ISSR fingerprint analysis on samples collected in 1994 and 2004 showed that each group had a single unique ISSR band. AMOVA analysis showed that 21% genetic variation occurred between the two periods; the drought period had lost 5.5% of polymorphic loci in comparison with the rainy period (Tables 2 and 3). Furthermore, the 1994 and 2004 collections were divided as two distinct groups, with 77.42% of the 1994 specimens and 68.75% of the 2004 specimens separated by PCoA. These results suggest that one decade of environmental changes has caused genetic structure and biometrical variation of cyst (see Asem et al., 2010) in this population.
Theoritically a decreasing population size in response to unfavorable ecological changes is expected to lead to a reduction of genetic diversity (Bálint et al., 2011;Cobben et al., 2011;Pauls et al., 2013). But in our study, the population genetic indices generated different patterns of genetic variation. Generally, results have demonstrated that the genetic variation of ITS1 in the drought period was reduced when the salinity of the lake was increased near saturation (300 g/l). In contrast, the genetic diversity of COI and 16S was not significantly different between the two periods. Additionally Na + /K + ATPase revealed a remarkable increase in variation in the drought period. These conflicting results might be attributed to the difference in the potential of gene variability that might be confirmed by further experimental evidence. Our study showed a negative and significant Tajima's D value for both examined periods in COI and ITS1 (Table 1), which indicated an excess of rare haplotypes resulting from population expansion or from selective sweeps (Nei & Kumar, 2000;Swanson, Aquadro & Vacquier, 2001;Akey et al., 2004;Cruciani et al., 2008;Levitan & Stapper, 2009). The negative values of Tajima's D should be referred to the demographic expansion of these markers in both ecological periods with regard to developed haplotype networks of ITS1 and COI markers. Additionally, Fu and Li's D* and Fu's Fs tests showed a negative departure from the neutrality test. Given that Fu and Li's D * test and Fu's Fs test recognize an excess of rare historical mutations (Fu & Li, 1993;Fu, 1996;Zhao et al., 2008), and rare recent mutations (Fu, 1997;Ramos-Onsins & Rozas, 2002;Zhao et al., 2008), respectively. Therefore the excess of both old and novel mutations were confirmed in the gene pool of COI and ITS1 during the rainy and drought periods. Consequently, the results of neutrality tests of Fu and Li's D* and Fu's Fs could explain the patterns of haplotype networks and population expansion of ITS1 and COI.
The major influence of short-term ecological disturbance was observed in the demographic history of Na + /K + ATPase and 16S markers. Neutrality tests resulted in a non-significant value for Na + /K + ATPase in the rainy period, which supported the demographic equilibrium. While significantly negative Tajima's D, Fu and Li's D* strongly supported a demographic expansion and an excess of rare historical mutations in the drought period. These results are consistent with the pattern of haplotype distributions developed during this period. In addition, the non-significant value of Fu's Fs suggested the absence of recent mutations in Na + /K + ATPase in both periods.
Another major alteration was observed in the 16S structure. Tajima's D and Fu and Li's D* were non-significant in the rainy period which could indicate 16S marker was at demographic equilibrium without selection in the rainy period. In contrast, a negative and significant neutrality value and expanded haplotype network indicated that 16S is involved in recent expansion in the drought period (Fig. 6). The negative and significant value of Fu Fs test suggested the excess of new mutations in the gene pool of 16S in both normal and drought periods.
Overall, our results have demonstrated that ecological disturbance should be considered in hypotheses about effects of short-term environmental changes on genetic variation. The rapid genetic changes that we found has also been demonstrated in some other species of animals and plants experiencing environmental crises (Avolio, Beaulieu & Smith, 2013;Brown et al., 2013;Banks et al., 2013). This could be attributed to the hereditary potential of populations respond to immediate ecological changes (Thompson, 2009). Although previous studies have shown a decreasing genetic variation in response to ecological disturbance, we found that Artemia urmiana shows dissimilar responses to environmental changes. Consequently, changes in genetic diversity and the pathway of variation are controlled by interaction between ecological conditions and the ability of genes to vary. Rogers (2015) suggested phenotypic patterns can be affected by ecological conditions which may cause genetic variation within an anostracan population during different periods. It is evident that the ecological crisis at Urmia Lake has had a meaningful influence on Artemia urmiana genetic structure, especially reducing genetic diversity, which ultimately could risk the survival of this crustacean.