Bird migratory flyways influence the phylogeography of the invasive brine shrimp Artemia franciscana in its native American range

Samples, laboratory procedures, and sequences

We obtained samples from 39 A. franciscana populations across its American geographical distribution, from Canada to Chile and Argentina, including Caribbean islands and a population from Cape Verde (Table 1 and Fig. 1). Most samples were cysts obtained from the ‘cyst-bank’ of the Instituto de Acuicultura de Torre de la Sal (CSIC, Castellón, Spain), collected between 1984 and 2000. Four Canadian samples were collected in the field in 2009, two of them (MANW and CHAP) as adults, which were preserved in absolute ethanol until needed. An additional cyst sample from Mono Lake (USA), collected in the 1970s was kindly provided by the Artemia Reference Centre (ARC 270). A few cyst samples that yielded poor quality DNA extractions were subject to hatching and the resulting nauplii used for DNA extractions (i.e., MexCB and GUA samples).

DNA extractions were carried out on individual cysts (previously rinsed in distilled water), whole nauplii or partial adults using a HotSHOT protocol (Montero-Pau, Gómez & Muñoz, 2008). We used Artemia-specific primers 1/2COI_Fol-F and 1/2COI_Fol-R (Muñoz et al., 2008) to amplify and sequence a 709 bp fragment of the mitochondrial Cytochrome c Oxidase Subunit I gene (COI). We also amplified and sequenced a 535 bp fragment of the 16S ribosomal RNA gene for a subset of individuals carrying different COI haplotypes to facilitate comparison with other published sequences using primers 16Sar-5′/16Sbr-3′ (Palumbi, 1996). PCRs were performed in 20 µL total volume containing 1× reaction buffer, 2.0 mM MgCl₂, 0.2 mM dNTPs, 0.6 units Taq DNA polymerase (Bioline, London, UK) and 0.5 µM of each primer. PCR conditions were as follows: 94°C for 3 min, followed by 35 cycles of 45 s at 94°C, 60 s at 45°C (60–64°C for 16S locus), and 60 s at 72°C, followed by 5 min at 72°C. PCR products were purified for sequencing using ExoSAP-IT^® (Exonuclease I and Shrimp Alkaline Phosphatase in buffer; USB Corp., Ohio, USA), cleaned with Sephadex^®-G50 (GE Healthcare Corp.), and labelled using the BigDye Terminator Sequencing Ready Reaction v3.1 kit (Applied Biosystems). The resulting fragments were separated on an ABI 3130xl genetic analyzer. Sequences were checked, edited, and aligned using Sequencher^® v4.5 (Gene Codes Corporation, Ann Arbor, MI, USA). All sequences were deposited in GenBank [accession numbers KF662951–KF663043 and KF725843–KF725869]. Available published sequences of the same gene fragments, to which we could assign a geographic origin, were also included in our phylogenetic analyses [GenBank: DQ401259–DQ401278, GU248382–GU248387, FJ007820–FJ007834, AF202735–AF202753].

Genetic analyses

Neighbour Joining (NJ) and Maximum Likelihood (ML) phylogenetic trees were inferred for both COI and 16S gene fragments. NJ trees were constructed using evolutionary distances computed using the Maximum Composite Likelihood method and 1000 bootstrap replicate tests in MEGA5 (Tamura, Peterson & Peterson, 2011). The best-scoring ML trees for COI were estimated using RAxML-VI-HPC v. 7.2.8 (Stamatakis, 2006) on the CIPRES portal at the San Diego Supercomputer Center (http://www.phylo.org), optimising free model parameters and executing 1000 rapid bootstraps. Average genetic distances between the main COI lineages – corrected by the K2P + G substitution model – were carried out using MEGA5. Additionally, to identify lineages in our COI dataset, we used the Automatic Barcode Gap Discovery (ABGD) approach (Puillandre et al., 2012) using the webtool (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html).

Intra-population gene diversity Hs (standardized haplotype diversity) for COI was computed using the program RAREFAC (Petit, 1998) to account for population differences in sampling size. Nucleotide diversity and pairwise Φ_ST values for COI from all 39 sampled populations were estimated using ARLEQUIN v. 3.1. (Excoffier, Laval & Schneider, 2005).

Testing isolation-by-distance patterns and effect of bird flyways.

The significance of correlations between pairwise genetic and geographic distances (isolation-by-distance or IBD patterns) was tested using Mantel tests on IBDWS v.3.21 (Jensen, Bohonak & Kelley, 2005). Prior to analyses in this section, populations inferred to be introduced intentionally by humans (see Table 1) were removed. For all sampled locations, precise decimal longitude and latitude coordinates were obtained using Google Earth (http://earth.google.com). A geographic distance matrix was then computed using Geographic Distance Matrix Generator v.1.2.3 (Ersts, American Museum of Natural History, Center for Biodiversity and Conservation, http://biodiversityinformatics.amnh.org/open_source/gdmg). We used a population geographic distance matrix (Table S1) and a population genetic distance matrix for COI data (Φ_ST values using Kimura 2-Parameter as the evolutionary model, Table S2). The 99% confidence intervals for the slope and intercept were estimated using Reduced Major Axis (RMA) regression with 30,000 bootstrap randomizations using log km geographical distance.

We used canonical Redundancy Analysis (RDA) in CANOCO (ter Braak & Šmilauer, 2002) to estimate the relative contribution of geographic distances, migratory flyways and human introductions on genetic distance between populations. RDA, a multiple linear regression method widely used in community ecology, has recently been applied to infer the role of spatial versus environmental variables in structuring population genetics data (e.g., Orsini et al., 2013). As the dependent matrix we used the sample loadings of a Principal Components Analysis calculated on Φ_ST values using Kimura 2-Parameter as the evolutionary model for COI data. Environmental variables were whether the flyway overlapped with the Artemia population (0 or 1 depending on the presence of birds from the Atlantic, Central or Pacific flyways in the area) and introduction history (0 or 1). We modelled spatial variables using latitude and longitude (x and y).

We used the overlap of sampled populations with the three main migratory flyways in America based on two sources (1) Boere & Stroud (2006) for shorebirds and (2) Birdlife International (extracted from http://www.birdlife.org/datazone/userfiles/file/sowb/flyways/). These data were used as a proxy for bird movements between locations. This approach is likely to be a rough approximation to the probability of bird migration between locations, given that more precise estimates cannot be obtained due to the absence of sufficient shorebird ringing data for the whole of the Americas. Even if extensive ringing data was available, this would only estimate current bird movements, whereas climate changes since the Pleistocene are likely to have strongly affected bird movements (Alerstam, 1993).

RDA analyses were carried out using each set of migratory flyway data (i.e., Boere & Stroud (2006), and Birdlife International information) and the variance partitioning calculated according to Borcard, Legendre & Drapeau (1992) when the model was significant. All environmental variables contributed to the full model, so we constructed two additional RDAs considering only flyway or introduction history as environmental variables (data files used in RDA are deposited in Dryad DOI http://dx.doi.org/10.5061/dryad.7kb11).

Phylogenetic relationships and geographic distribution of lineages

Once PCR primers were removed and sequences trimmed to the same length, the 604 bp COI alignment contained 603 individuals newly sequenced in this study (see Table 1), which collapsed into 93 haplotypes. No ambiguities, insertions, deletions or stop codons were present in the alignment. There were 121 variable sites, 86 of them parsimony informative, and 104 synonymous and 15 non-synonymous substitutions.

Using the default parameters, ABGD did not find any partitions in our dataset, so we reduced X, (i.e., the minimum barcode gap width) as suggested by the program, to 1.0 and used Kimura 80, identifying 12 groups with a prior intraspecific distance of = 0.0028 (we use the term ‘lineages’ thereafter, see below for their geographic distribution). When the prior intraspecific distance increases to 0.0046 (not high by any standards) then the number of partitions reduces again to 1.

Both phylogenetic approaches ML and NJ were highly consistent, with two highly supported main branches displaying a geographic distribution along the continent, one mainly Atlantic (lineages 9–12) and the other split between two sub-branches along the Pacific Rim (lineages 1–7) and in Central Canada (lineage 8). Overall, there were at least 12 mtDNA lineages, most of them well supported (Fig. 2). With the exception of lineage 1, each of the lineages showed a restricted geographic distribution indicating a high level of regional endemism (Fig. 3 for the geographic distribution of the lineages). Lineage 2 was found in a single coastal site in NE Mexico, lineage 3 in the five locations in Peru, lineage 4 in two high altitude populations from Central Mexico, lineages 5 and 7 in Central Chile, lineage 6 in a single sulphate-rich location in NE Mexico, lineage 8 in four locations in Saskatchewan, Canada, lineage 9 in Puerto Rico, lineage 10 in Cape Verde, lineage 11 in several locations around the Caribbean, Baja California and SW Mexico, and lineage 12 in Argentina and Chile. In stark contrast to the rest, lineage 1 was genetically diverse and geographically widespread, found across a large geographical area across both sides of the continent, including Brazil, Chile, Cuba, Jamaica, Mexico, USA and Colombia. Out of the 27 haplotypes in lineage 1, seven haplotypes were detected in SFB and GSL, the two commercialised populations in the USA that have been sources for the invasion in the Mediterranean (J Muñoz, A Gómez, J Figuerola, F Amat, C Rico, AJ Green, 2013 unpublished) (see section Introduced populations below). In addition, lineage 1 contained six closely related haplotypes only found in Mono Lake – which harbours a population often considered as a separate species, A. monica – and haplotypes found in three Colombian populations (CSC, CT and CGZ), five Mexican populations, and the Jamaican population. Two other lineages were distributed across the continental E-W divide, creating contact zones between lineages. In lineage 12, most of the haplotypes were found in Argentinean populations, but four of these haplotypes (three of them shared across populations) were found in a Chilean Altiplano population (Salar de Llamara). Lineage 11, mainly Caribbean, has a few haplotypes in two Mexican populations from the Pacific side (Las Coloradas and Faro San José, where they coexist with lineage 1). Finally, although both lineages 1 and 11 are found in Colombia, they were not found together in any of the populations sampled. The genetic divergence between the 12 main lineages ranged from 1.8% (between lineages 1 and 4) to 6.0% (between lineages 2 and 9) (Table 2).

The 16S alignment contained 408 bp from 122 individuals, which collapsed into 59 haplotypes. There were two singleton indels, 63 variable sites and 43 parsimony informative sites. In contrast to the COI analysis, the NJ and ML reconstructions were poorly resolved, especially the basal branches, but most of the lineages recovered by the COI analyses were also recovered for the 16S data, with variable levels of support (see Fig. 4). COI lineages 3, 5, 7 and 9 were highly supported for both ML and NJ analyses in the 16S analysis, whereas lineages 1 and 2 on the one hand and lineages 4 and 6 on the other, collapsed into poorly supported branches. The 16S analyses allowed us to assign several previously sequenced populations, which we were unable to sample, to COI lineages, particularly in NW America and the Caribbean. In addition, the 16S analysis revealed the presence of two new lineages in the Caribbean, one in the Virgin Islands, related to lineage 9 (Puerto Rico) and another in Inague Island (Bahamas) related to lineage 12 (Argentina, Chile). Note that these maintain the relationship with Atlantic lineages. Regarding North-Western American populations – extensively sampled by Prof Sarane Bowen’s group – in New Mexico, Nevada and British Columbia, they hold private haplotypes which appear in a poorly supported branch with Mono Lake and Mexican haplotypes. Other populations (Clinton, Basque Lake, Baja California and Carrizo Soda lake) also appear in the composite lineage 1 and 2, underscoring the diversity of USA A. franciscana populations. Interestingly, 16S haplotypes from Jesse Lake (Nebraska) belong to the Central Canadian lineage (lineage 8) together with Little Manitou, Muskiki, Meacham and Chaplin haplotypes.

Intra- and inter-population genetic diversity in COI

The number of individuals sequenced per population ranged between 4 and 37, depending on available material (average 15.97; see Table 1 for estimates of π and H, and details of the haplotypes present in each population). The number of haplotypes per population ranged from one to seven. The highest standardized gene diversity (Hs) was found in the Mexican population MexH, whereas five populations (MexCe and Mex99 from Mexico, BRM from Brazil, PV from Peru, AMC, from Argentina, and CHLC from Chile) were fixed for a single haplotype. Most haplotypes were found within single countries, except for several haplotypes shared between some countries and the two commercial USA populations SFB and GSL (see Table 1).

Populations were strongly structured genetically (global Φ_ST = 0.92; 0.94 when putative introduced populations were removed), with Φ_ST being highly significant between all populations except for three lakes from Central Canada, plus one pair from Chile/Argentina (see Table S2).

Identification of introduced populations

We found genetic evidence for putative non-native populations originating from the commercialised SFB or GSL in nine sites from Mexico, Cuba, Jamaica, Brazil and Chile. In these populations, at least one sampled individual shared a haplotype with SFB and/or GSL populations (see Table 1). Those populations showed three patterns: (1) all individuals sampled shared haplotypes with SFB and/or GSL (BRM, GUA, and CHLC); (2) populations had haplotypes shared with the commercialised populations and a further haplotype (Af19 for MEXLC and JAYA) which differs from Af18 (a common haplotype in SFB and GSL) by a single substitution; and (3) populations sharing some haplotypes with SFB and/or GSL, but which also had unrelated additional/private haplotypes (MexT, MexSQ, MexFSJ, and CHPI). For four of the nine putative introduced populations, the occurrence of intentional introductions had been previously reported either in the same or nearby sites (see references in Table 1). Note that introduced A. franciscana populations are likely to further expand into nearby suitable habitats due to passive dispersal by birds.

Isolation-by-distance pattern and the role of American migratory flyways

Mantel tests on pairwise genetic and geographic distances for populations ranging from Chile-Argentina to Canada, excluding those inferred to be introduced, revealed a strong IBD pattern (Fig. 5) with a highly significant correlation between pairwise geographic and genetic distances, indicating that geographic distance between populations explains a large proportion of the genetic variability in the sampled area (R²-value = 0.323, p-value < 0.001).

RDA showed that both flyway and introduction history were significantly associated to population genetic distance (p < 0.02 for all the correlations with genetic distance, whether or not geographical distance was controlled for). The effect of using different flyway data was quite small. Flyway explained 18.7–21.2% of variance in population genetic distance, depending on the flyway dataset used (Table 3). Genetic distance was also affected by historical anthropogenic introductions (15.6% of variance explained) and both factors together (i.e., flyway and anthropogenic introduction) explained 30.6–31.2% of variance in genetic distance. In comparison, geographic distance only explained 8.2–8.7% of variance.

The role of bird migratory flyways

Our study provides new evidence supporting the key historical role of waterbirds as the main factor shaping the population genetic structure of continental aquatic invertebrates at an intra-continental scale. The patterning of the main phylogenetic lineages, with an Atlantic, Central and Pacific distribution – instead of a North American (Nearctic) vs. South American (Neotropical) division reflecting the recognized zoogeographic regions and the long isolation of the continents (Lomolino et al., 2010; Holt et al., 2013) – strongly suggests that historical bird migratory flyways, which occur alongside both the major coasts of this continent, determined the historical spread of A. franciscana genetic lineages. Bird movements also might have allowed the subsequent persistence of this structure by facilitating colonisation along each migratory flyway, which shaped the main East-West division in mitochondrial lineages. RDA showed that the effect of migratory flyways was highly significant, and accounted for 20% of the genetic variation between populations once geographic distance was taken into account, suggesting that the distribution of genetic lineages in A. franciscana is likely to reflect the impact of historical bird flyways on native phylogeographic patterns. In addition, the strong detected IBD pattern suggests that the chances of bird-mediated colonisation are highly distance dependent (see Viana et al., 2013), although instances of long-distance dispersal and colonisation, for example from Argentina to Chile, or to Colombia from northern Caribbean populations, are also apparent from our data. A corollary of our results is that bird movements must have shown some stability, forming parallel N-S flyways during the time frame of A. franciscana population diversification (i.e., throughout the Pleistocene) extending into new breeding areas becoming available in the north and new wintering sites in the southern extreme of the continent (Buehler, Baker & Piersma, 2006). These results are no surprise given the high transport rates of Artemia cysts by waterbirds (Green et al., 2005; Sánchez et al., 2007) and the lack of other dispersal vectors. Migratory waterbirds have existed since the Early Cretaceous (Lockley et al., 2012). Although now-extinct migratory mammals were once major vectors of plant dispersal in the Americas (Janzen, 1984), the hypersaline habitats used by Artemia are not conducive to dispersal by large mammals.

Our results for A. franciscana agree with other studies in North America pointing to an effect of bird movements on the genetic structure of passively dispersed aquatic invertebrates (Figuerola, Green & Michot, 2005), but extends these findings to the whole of the Americas, and emphasizes the major role of bird movements in facilitating colonisation into new suitable habitats – in agreement with the patterns found in A. franciscana in the invaded Mediterranean range (J Muñoz, A Gómez, J Figuerola, F Amat, C Rico, AJ Green, 2013 unpublished), and those for A. salina in its native range (Muñoz et al., 2008). Although at a continental scale our results suggest that bird movements do not promote contemporary gene flow between Artemia populations, such gene flow may still have a role at a local scale or when new areas suitable for colonization become available.

Phylogeographic patterns

The high level of endemism and population structure in native A. franciscana populations, with low ongoing gene flow and occasional long-distance migration resembles the patterns found in Artemia salina (Muñoz et al., 2008), a sexual native Mediterranean brine shrimp species, and in a range of other aquatic passively dispersed taxa including sexual reproduction in their life cycles (Gómez, Carvalho & Lunt, 2000; Edmands, 2001; De Gelas & De Meester, 2005; Mills, Lunt & Gómez, 2007; Ketmaier et al., 2008; Korn et al., 2010). Our results also expand and confirm the deep phylogenetic breaks found by Maniatsi et al. (2009) in a mtDNA and nDNA study of a more limited number of populations revealing only three lineages. As populations of these passively dispersed organisms can be founded by a small number of propagules, followed by rapid population growth and establishment of large diapausing egg banks, this favours the presence of persistent/long term founder effects, thus reducing the effect of gene flow, possibly reinforced by the development of local adaptation, according to the Monopolisation hypothesis (De Meester et al., 2002). Cysts are undoubtedly still regularly dispersed between suitable habitats by waterbirds, but they are unlikely to become established owing to the Monopolisation effects.

Given the range of genetic divergence between lineages (from 2 to 6%) the time frame of their fragmentation can be approximated roughly using a COI molecular clock for other shrimp taxa (1.4% sequence divergence per million year; Knowlton & Weigt, 1998), which translates into 1.4–4.3 million years of divergence (reaching the Pliocene), between A. franciscana lineages. Even a faster rate of 2% per million years will result on pre-Pleistocene divergence times between the main lineages. Therefore, a contribution of Pliocene/Pleistocene climatic oscillations to population fragmentations after range expansions across the continent can be inferred from our data, possibly allowing survival of lineages in separate geographical areas including Caribbean islands and areas in North and South America. Mexico has the highest lineage richness, with five out of the 12 COI lineages being native to this country – including lineages from both Pacific and Atlantic clades. These findings suggest that this region is likely to have supported separate refugia during climatically adverse periods. The occurrence of a highly divergent Central Canadian prairie lineage was unexpected, as an ice sheet covered this area during the last glacial maximum (Ehlers & Gibbard, 2004). However, the 16S data from GenBank suggest that this lineage also occurs in more southern areas in central USA (Nebraska), where it may have survived south of the ice sheets, and then undergone postglacial colonisation of Central Canada.

Following climate-driven turnover of hypersaline habitats, migratory shorebirds would be involved in expanding the lineages into newly appearing suitable habitats and the chances of successful spread would be strongly distance dependent. However, long-distance colonisation events must have also occurred. For example, assuming that the ancestor of lineages 9 and 10 inhabited the Caribbean area, the colonisation of Argentina by the ancestor of lineage 12 must have entailed successful transfer of some cysts between these distant areas (see Fig. 2). The genetic composition of lineage 12, where most of the haplotypes were found in Argentinean populations (see also Fig. 3) but four of them (three shared across populations) were found in a Chilean Altiplano population (Salar de Llamara), suggests recent colonisation from Argentina to Chile. Finally, the colonisation of Cape Verde Islands, with private haplotypes distantly related to Caribbean lineages, must have involved long-distance transport, possibly from birds accidentally landing there after storms sent them off course, although we cannot rule out the possibility of an unreported human-mediated introduction involving cysts from a Caribbean population not included in our study. Shorebirds are likely vectors for such long-distance dispersal events, and have often been implicated in the dispersal of plants between North and South America, or to oceanic islands (Cruden, 1966; Proctor, 1968).

The 16S analysis also shows higher richness of lineage 12 in Argentina with two Chilean populations (Convento and Salar de Llamara), for which evidence of nuclear DNA (nDNA) introgression among lineages exists (Maniatsi et al., 2009). Furthermore, as migratory flyways overlap on some areas, this could have occasionally resulted in transfers from the Atlantic to the Pacific coasts, as has been suggested for other passively dispersed aquatic species (Miura et al., 2011).

Natural spread of lineages from refugial areas is likely to result in contact zones between lineages, which we expect to be sharp, as we found between lineages 1 and 11 in Colombia, where despite two lineages being present in the area, there are no sites where both co-occur.

Taxonomic considerations

The COI gene is one of the most widely used tools for species delineation (Hebert et al., 2004; Costa et al., 2007; Gómez et al., 2007b). Sequence divergence of 3% have been proposed as a threshold for species delimitation in crustaceans (Costa et al., 2007, but see Lefébure et al., 2006), but other approaches are also used, such as GMYC (Pons et al., 2006) or automatic barcode delimitation, which we used here (Puillandre et al., 2012). Our analysis revealed 12 lineages in A. franciscana, some of them, like lineages 9 (Puerto Rico) plus 10 (Cape Verde) compared to all the others, or lineages 11 (Circum Caribbean) plus 12 (Argentina-Chile) compared to all the others, have genetic divergences of over 5%. Surprisingly, reproductive isolation – due to ecological specialisation and local adaptation – has only been reported between Mono Lake and San Francisco Bay populations (Bowen, Buoncristiani & Carl, 1988) due to the inability of individuals of each of these populations to survive in each others’ ecological conditions, but our data show they are very closely related. Indeed cross-fertility has been observed in the laboratory between the San Francisco Bay population and 15 other populations from the whole range of the species, including some populations included in our study that belong to very divergent mtDNA lineages such as Inague Saltern (Bahamas), Little Manitou (Canada, lineage 8), and Puerto Rico (lineage 9) (Clark & Bowen, 1976). Therefore, we concur with Bowen, Buoncristiani & Carl (1988) in regarding A. franciscana as a very diverse “superspecies”, where reproductive isolation mediated by habitat adaptation might occur, but populations in intermediate habitats could act as venues for genetic exchange between ecological isolates.

Impact of introductions on native population structure and management implications

Given the strong phylogeographic structure of A. franciscana, and the high level of private haplotypes found, we used haplotype sharing between the commercialised populations (SFB and GSL) and distant populations as a criterion of recent human mediated introduction. Using this criterion, we identified nine populations where genetic evidence pointed to putative human introductions into Mexico, Cuba, Jamaica, Brazil and Chile. For four of these populations, the occurrence of intentional introductions in the same or nearby sites during the 1970s could be documented. Therefore, our genetic data confirms that the established A. franciscana populations in these locations are, at least partially, of introduced origin, and validates our criterion. As for the impact of introductions on native populations, four Mexican populations contained additional private haplotypes not closely related to the introduced ones, which suggests the presence of a pre-existing population before the introduction and the likely introgression of both populations with persistence of native haplotypes, although nuclear loci would be necessary to confirm this. In two populations, we found a discrepancy between our genetic results and previously published ones. The first case is the MEXLC population, which we regard as introduced (lineage 1), whereas in Tizol-Correa et al. (2009) the haplotype found is shared with populations from our lineage 11. The second case is CHLC (Laguna Cejas, Chile), which we regard as introduced whereas Maniatsi et al. (2009) found that their mtDNA was most closely related to native populations from Central Chile (which belong to our lineage 7). While these authors concluded that the discrepancies between their nDNA and mtDNA data from the latter population pointed to incomplete lineage sorting, an alternative explanation is that they are due to the population being admixed with native and introduced ancestry, as introduced A. franciscana is known to hybridize even with the genetically divergent A. persimilis (Kappas et al., 2009). This failure to detect introduced mtDNA lineages in this population might arise from reduced sample sizes, or indicate real temporal changes in these populations reflecting recent introductions, as samples are likely to have been obtained in different years (unfortunately, no collection dates for these population is reported in Maniatsi et al., 2009). Despite human introductions, haplotypes presumably from pre-existent populations have survived and coexist with introduced ones, although a possible loss of genetic diversity due to introductions cannot be ruled out, and should be investigated in the future, perhaps using sediment cores (a method widely used in the crustacean model organism Daphnia) from populations where commercialised strains were introduced.

The high genetic richness found in our COI analyses, and the presence of private haplotypes belonging to lineage 1 in distant populations away from SFB and GSL, including three Colombian and five Mexican populations likely to be native, indicates that the natural distribution of lineage 1 extended further than the two commercialized populations before any human introductions took place. This was also in evidence from our 16S analyses, which included more extended geographic coverage of NW America, and revealed further sites within lineage 1 in the W USA and British Columbia harbouring private haplotypes not found in SFB or GSL. This is also consistent with the presence of private and closely related haplotypes at Mono Lake. The peculiarity and fragmentation of the habitats used by the species, the potential of salinity and varying ionic composition to act as a strong selective agent, the capacity to produce massive quantities of resting eggs that can be readily dispersed by birds, combined with the apparent limitation in modern gene flow, makes this group an ideal system for further studies testing the role of local adaptation and mass effects on reducing gene flow between populations.

Given the impact of the invasive A. franciscana across the world, and the high genetic and ecological richness of its native populations, further population translocations should be highly discouraged, and the use of native strains as a source of cysts should be encouraged even within the Americas.