Phylogenetic analysis of brine shrimp (Artemia) in China using DNA barcoding.

DNA barcoding is a powerful approach for characterizing species of organisms, especially those with almost identical morphological features, thereby helping to to establish phylogenetic relationships and reveal evolutionary histories. In this study, we chose a 648-bp segment of the mitochondrial gene, cytochrome c oxidase subunit 1 (COI), as a standard barcode region to establish phylogenetic relationships among brine shrimp (Artemia) species from major habitats around the world and further focused on the biodiversity of Artemia species in China, especially in the Tibetan Plateau. Samples from five major salt lakes of the Tibetan Plateau located at altitudes over 4,000 m showed clear differences from other Artemia populations in China. We also observed two consistent amino acid changes, 153A/V and 183L/F, in the COI gene between the high and low altitude species in China. Moreover, indels in the COI sequence were identified in cyst and adult samples unique to the Co Qen population from the Tibetan Plateau, demonstrating the need for additional investigations of the mitochondrial genome among Tibetan Artemia populations.

Introduction

Brine shrimps (Artemia) live in hypersaline environments, such as salt lakes and lagoons with a worldwide distribution. They are often used as model organisms to evaluate management policies in aquatic resources 1., 2. and also in biodiversity studies in inland hypersaline lakes, since they show enormous diversity at the genus level in terms of their ability to survive under different ionic compositions, climatic conditions, and altitudes 3., 4., 5.. In addition, these morphological variations provide excellent materials for studying adaptive genetic polymorphisms at the molecular level. During the past two decades, the phylogenetic relationships among Artemia species have been established by combined studies based on crossbreeding, morphological differentiation, cytogenetics, nuclear markers (including allozymes and other nuclear DNA sequences), and mitochondrial DNA sequence markers 6., 7., 8., 9., 10..

DNA barcoding is a technique that uses a short genetic marker in an organism’s mitochondrial DNA to identify if it belongs to a particular species. It has been used as a powerful approach for characterizing species of organisms, especially those with almost identical morphological features (. Here we report a phylogenetic study on Artemia samples obtained from different habitats around the world based on DNA barcoding procedures. A 648-bp segment of the mitochondrial gene, cytochrome c oxidase subunit 1 (COI), was chosen as a standard barcode region (. Our initial study focused on the plateau regions in China, particularly the isolated salt lakes in the Tibetan Plateau. We also investigated population diversity among the isolated salt lakes as well as the evolutionary relationship among Artemia populations in the Himalayas.

Results and Discussion

Phylogenetic analysis based on COI divergence

We compiled sequences from the mitochondrial COI gene of Artemia samples and yielded a high-quality alignment of a 474-bp segment. Most of our samples gave rise to consistent high-quality sequences without any background sequences (if exist, they are often observed as minor peaks bellow a series of major peaks in the sequence chromatographs) except for the samples ARC1526 and 1612, where some minor but above-background sequence heterogeneities were constantly observed. There are two possibilities for the background signal; one is the heterogeneity of mitochondrial sequences and the other is a result of co-amplification of a potential nuclear pseudogene. We prefer the first case for the following reasons. First, when we cloned the PCR product as a library and sequenced a few dozens of clones from it, we found that a small fraction of the clones contain indels (insertions or deletions), suggesting that the indels exist among the samples and are not sequencing errors. Second, we repeated DNA extraction from two batches of cyst samples and the result confirmed the initial observation. In addition, the same indels were also found in some adult samples. Third, if these indels are from a potential pseudogene in a nuclear genome, their existence would suggest either multiple integrations with sequence loss or multiple deletions of the same inserted fragment in this particular population. Such assumption seems less possible than a simple mechanism of sequence deletions and is an extremely rare event in the mitochondrion of this unusual population due to functional selection, since we have not found any other indels in all sampled populations. Therefore, we chose the most frequent alleles from these sequences for the construction of phylogenetic trees (see the last section for a more detailed discussion).

The COI sequence fragment among the Artemia groups collected around the world displayed considerable divergence with intraspecific Kimura-2-parameter (K2P) distances at an average of 14% (ranging from 0 to 26%). Based on the K2P distance and bootstrap values, both neighbor-joining (NJ) and maximum parsimony (MP) trees were built with effective sequence alignments of the 474-bp fragment within the barcode region (Figure 1), giving rise to a stable and comparable topology of five clades. The first clade mainly contains genotypes from populations collected near the Bohai Bay area in East China, which showed a high sequence similarity to A. franciscana. Because the samples were collected from different lakes with multiple sampling over time, a previous large-scale invasion of A. franciscana in Bohai Bay was confirmed (. The second clade mainly contains genotypes from populations in the Tibetan Plateau in Southwest China with high sequence similarity to A. urmiana, including samples from five major salt lakes located at altitudes over 4,000 m, representing the high altitude populations of Artemia in China. Since A. urmiana and parthenogenetic Artemia are genetically similar, and samples ARC1526 and 1612 have only females, it suggests that parthenogenetic Artemia may be present in the sample. Two consistent unique sequence variations resulting in amino acid alterations, 153A/V and 183L/F, were observed when we compared the high and low altitude species from China. These sequence alterations are candidates for further functional studies such as to determine if the adaptation to high altitude had resulted in the fixation of such mutations. The third clade corresponds to A. sinica, containing populations from Inner Mongolia in North China. The fourth and fifth clades correspond to A. persimilis and A. salina, respectively; however, no sample populations in China are found in the two clades in this study.

Fig. 1

Phylogenetic trees built by NJ method (A) and MP method (B). The trees were constructed based on the COI sequence fragment showing phylogenetic relationships among Artemia species using Daphnia pulex as outgroup.

Positive selection based on Ka/Ks analysis

The rate of substitution in synonymous and nonsynonymous sites (Ka/Ks ratio) provides a powerful tool for understanding DNA sequence evolution (. We used the GY method 14., 15. to estimate Ka/Ks, hoping to reveal sequence signatures of natural selection in the COI gene (Table 1). Compared with populations in the first clade, significantly higher Ka/Ks ratios were found among the Tibetan populations; their COI genes have relatively higher Ks values, indicating that they possess more synonymous mutations. The two variations (153A/V and 183L/F) shared by the high altitude group that alter amino acid sequences between the high and low altitude populations were also detected, suggesting that they were positively selected for better fitness in the Tibetan Plateau. Representative sequence of sample ARC1612 has the highest Ka/Ks ratio, suggesting a relatively stronger selection posed on this population, and synonymous mutations provide clues that the populations are diverging, most likely due to the changing environment during the last three million years rather than genetic drift (. In addition, we also found that populations from Inner Mongolia (ARC1188 and 1317) have significantly higher Ka values than those in the first clade, implying relatively stronger selective pressure on this species.

Table 1

Ka/Ks analysis results among different Artemia populations

Sequence	Ka	Ks	Ka/Ks
Daphnia pulex	8.00E-02	3.4873	0.0229
1612	1.22E-02	0.9703	0.0126
1526	9.17E-03	1.0071	0.0091
1524	6.14E-03	1.1643	0.0053
1348	6.14E-03	1.1643	0.0053
1610	6.12E-03	1.2177	0.005
1609	6.15E-03	1.3616	0.0045
A. urmiana	6.16E-03	1.386	0.0044
A. salina	1.28E-02	3.0613	0.0042
1611	6.20E-03	1.703	0.0036
1188	6.08E-03	2.7791	0.0022
1317	5.48E-03	3.5289	0.0016
1321	2.93E-03	2.0111	0.0015
1256	2.94E-03	2.1369	0.0014
A. persimilis	2.95E-03	2.2582	0.0013
A. sinica	3.03E-03	2.5426	0.0012
1588	7.12E-06	0.0071	0.001
1590	7.12E-06	0.0071	0.001
1592	7.12E-06	0.0071	0.001
1125	2.86E-05	0.0286	0.001
1233	2.86E-05	0.0286	0.001
1258	2.86E-05	0.0286	0.001
1261	2.86E-05	0.0286	0.001
1349	2.86E-05	0.0286	0.001
1520	2.89E-05	0.0289	0.001
1613	2.86E-05	0.0286	0.001
1580	1.43E-05	0.0143	0.001
1619	1.43E-05	0.0143	0.001

Divergence of the Tibetan populations

The Tibetan populations showed a clear difference from populations elsewhere in China. We obtained high-quality sequences from individual adults of the six Tibetan populations and calculated their K2P distances (Table 2). For the tree construction, we used the consensus sequences when sequence heterogeneities were encountered among a minor set of samples. Among the six populations, individuals from ARC1524 (Jingyu Lake) have an average distance of 0.51%, the lowest in comparison to the other populations, with the distance ranging from 0 to 1.31%. Individuals from ARC1526 and ARC1612, both from Co Qen, have an average distance of 4.17% and 7.07%, respectively, with the distance ranging from 0.38% to 11.07% for ARC1526 and 2.55% to 12.01% for ARC1612. Individuals from ARC1348 (Lagkor Co) and ARC1609 (Nima) showed similar degrees of divergence to each other. The higher diversity in Co Qen may reflect a mixed species status in the population as a result of either a significantly longer genetic history or higher environmental pressures. Based on phylogenetic and divergence analyses of the selected samples from different regions of the world, it is possible that the high altitude group of Artemia are descendents of a local ancestral species in the Himalayas, which diverged genetically as the Tibetan Plateau arose stepwise over approximately the last three million years (.

Table 2

K2P distance showing the COI sequence divergence of six Tibetan populations

ARC No.	Sample quantity	Min	Max	Mean	S.D.
1348	20	0	7.76%	2.30%	0.34%
1524	18	0	1.31%	0.51%	0.17%
1526	9	0.38%	11.07%	4.17%	0.52%
1609	20	0	4.92%	2.49%	0.45%
1610	20	0	9.56%	3.42%	0.50%
1612	8	2.55%	12.01%	7.07%	0.79%

Sequence indels in the COI gene unique to the Co Qen population

Several sequence indels in the COI gene were identified in the cyst and adult samples from the Co Qen population (Figure 2). Such indels certainly abolish the function of the gene as null mutations, even a triplet in frame. Since performing sequencing directly on PCR product (using it as sequencing template) is incapable of providing information on allelic frequencies or ratios, we cloned the PCR product and sequenced several clones from each PCR-product-derived library. Initial analysis of DNA from adult samples revealed 2 single-nucleotide deletions, 2 three-nucleotide insertions, 27 single-nucleotide insertions, and 25 normal alleles (polyC with n=7) among 56 clones. These sequence indels are not likely all to be sequencing errors since they were found only in this particular population and a track of seven guanines are usually not a difficult task for thermo-sequencing. Furthermore, we found that almost equal proportions of the alleles with or without single-nucleotide indels exist within the population. We subsequently examined the cyst samples using the same PCR product cloning procedure. Surprisingly, the percentage of indel-bearing alleles vs. normal alleles was 10%, which is quite different than what was seen in adults. It is suggested that the indels, most of which are single-nucleotide insertions, should already exist within mitochondrial populations in an individual Artemia cyst or adult. It is not known if this is a result of natural selection or genetic drift, depending on how they alter the overall functionality of the COI gene. Although one may argue that there is still a possibility that these indels are a result of potential co-amplification of nuclear pseudogenes due to very recent mitochondrial gene transfers, our experimental data are not in favorite of this far-fetched assumption. Nevertheless, to understand why this phenomenon only happens to the Co Qen population and how the indels are generated in COI gene, a full-scale sequence analysis on the Tibetan populations as well as functional studies on the deleterious protein are of necessity.

Fig. 2

Alignment of partial COI sequences including the indels among adult samples from the Co Qen population (ARC1526 and 1612). Initial analysis of DNA from adult samples revealed 2 single-nucleotide deletions, 2 three-nucleotide insertions, 27 single-nucleotide insertions, and 25 normal alleles (polyC with n=7) among 56 clones.

Materials and Methods

Sample collection

All Artemia strains used in the present study are kept as cysts at the Laboratory of Aquaculture & Artemia Reference Center (ARC), Ghent University, Belgium, with ARC code numbers (Table 3). Among them, six populations were obtained from five salt lakes in the Tibetan Plateau, namely Lagkor Co (ARC1348), Jingyu Lake (ARC1524), Co Qen (ARC1526 and 1612), Nima (ARC1609), and Yangnapeng Co (ARC1610). Cysts of these six populations were hatched and the resulting larvae were allowed to develop into adulthood. The hatching and survival rates are shown in Table 4. Ten female and ten male adults were chosen randomly from each of the four populations (ARC1348, 1524, 1609, and 1610). Because the two populations from Co Qen (ARC1526 and 1612) produced only females, twenty female adults were used for our experiments. All adult samples were stored in ethanol, allowing distant transportation before DNA extraction.

Table 3

List of Artemia samples in this study with collection location and time

ARC No.	Location	Harvesting year
1125	Tanggu, Tianjin, China	1989
1188	Yimeng, Inner Mongolia, China	1991
1233	Gangzhou, Hebei, China	1991
1256	La Pampa, Argentina	1993
1258	San Francisco Bay, California, USA	1993
1261	Luannan, Hebei, China	1993
1317	Bameng, Inner Mongolia, China	1995
1321	Argentina	1996
1348	Lagkor Co, Tibet, China	1996
1349	Vinh Chau, Vietnam	1997
1520	Great Salt Lake, Utah, USA	2001
1524	Jingyu Lake, Xinjiang, China	2001
1526	Co Qen, Tibet, China	2001
1580	Wudao, Shangdong, China	2002
1588	Beidaba, Shangdong, China	2002
1590	Huangnigou, Shangdong, China	2002
1592	Leguantai, Shangdong, China	2002
1609	Nima, Tibet, China	1999
1610	Yangnapengco, Tibet, China	2002
1611	Qixiang Lake, Tibet, China	2002
1612	Co Qen, Tibet, China	2001
1613	Luannan, Hebei, China	1992
1619	Chengkou, Shangdong, China	1991

Table 4

Hatching and survival rates of samples from Tibetan Plateau

ARC No.	Hatching rate	Survival rate
1348	29.72%	74.29%
1524	23.14%	67.68%
1526	19.00%	62.50%
1609	32.84%	74.59%
1610	28.72%	67.58%
1612	21.40%	64.89%

DNA extraction

DNA from cyst and adult samples were extracted respectively. On the one hand, 0.1 g cysts was ground into fine powder in the presence of liquid nitrogen and then 1 mL of 2X CTAB (2% CTAB, 0.1 M Tris-HCl at pH 8.0, 20 mM EDTA, and 1.4 M NaCl) was added immediately. On the other hand, adult samples were homogenized in 200 μL CTAB buffer and incubated at 65°C for 1 h. The viscous aqueous supernatant was carefully collected after a brief centrifugation at 9,500 g for 30 min at ambient temperature. DNA was extracted from the supernatant with chloroform/isoamylalcohol (24:1), precipitated with isopropanol, washed with ethanol, and re-suspended in double-distilled water containing 20 μg/mL RNase.

PCR amplification

We used a pair of primers based on the Artemia mitochondrial genome sequence (GenBank Accession No. NC_001620) (COI-upper: 5’-CTGATTGGCGATGAACA-3’, COI-lower: 5’-ATGCTTCCTTCTTACCTCTT-3’), which generated a 648-bp DNA product. PCR amplifications were carried out in a final reaction volume of 25 μL containing 2.5 μL of 10X PCR buffer, 1.8 μL of 25 mM MgCl2, 1.2 μL of 2.5 mM dNTP, 10 pM of each primer, and 0.5 unit of Taq DNA polymerase. The thermocycling profile was set as: 1 cycle of 5 min at 94°C, 35 cycles of 30 sat 94°C, 30 s at 50°C, and 1 min at 72°C, followed by a final extension of 10 min at 72°C. PCR products were analyzed based on agarose gel electrophoresis.

Sequencing of PCR products

DNA sequencing was carried out with 2 μL (50-100 ng) of PCR products after a desalting step with the Minipore’s 96-well cleanup plate. Sequencing reactions were set up in a 10-μL mixture with dyeterminator chemistry (ABI) and 10 pmol of the sequencing primer; in this case, we used the two PCR primers in separate sequencing reactions to obtain DNA sequences from both the 5’ and 3’ directions (or the strands). The sequencing amplification protocol is set up as follows: 1 cycle of 15 s at 95°C, 35 cycles of 15 s at 95°C, followed by a final annealing of 15 s at 50°C, and a 1.5-min extension at 60°C. Sequences were acquired with an ABI 3730XL sequencer. Raw sequence data were analyzed with the Phred-Phrap-Consed package 17., 18..

Phylogenetic and Ka/Ks analyses

The newly acquired sequences of A. franciscana, A. urmiana, A. sinica, A. persimilis, and A. salina (Artemia_franciscana_EF615571 to Artemia_persimilis_EF615594), together with their corresponding regions using Daphnia pulex as outgroup (available in GenBank database with gi number of gi|71482580|, gi|5835848|, gi|5835051|, gi|71482578|, gi|71482572|, and gi|90019508|, respectively), were aligned using the ClustalW package (, followed by manual inspections and quality assessments. A 474-bp consensus was used for the phylogenetic analysis. We used the Ka/Ks calculator ( for Ka/Ks analysis based on the K2P model and the GY method 14., 15.. NJ and MP analyses implemented in MEGA3.1 were employed to examine relationships among all the populations. In order to reduce the model effect, the NJ tree was built based on the K2P model as well as the Kimura-3-parameter (K3P), Jukes-Cantor (JC), and Tamura-Nei (TN) models for COI DNA sequences. We used the heuristic searches (close-neighbor interchange) with 500 random addition replications to construct the MP tree. The bootstrap test with 1,000 replicates was applied to assess the reliability of phylogenetic trees.

Authors’ contributions

WW, QL and HG carried out the study and prepared the manuscript. PB, GVS, PS and NX helped with sample collection and co-wrote the manuscript. QS, SH and JY supervised the research and revised the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors have declared that no competing interests exist.