# Brine chemistry matters: isolation by environment and by distance explain population genetic structure of Artemia franciscana in saline lakes

Dagmar Frisch, Christophe D Lejeusne, Masaki Hayashi, Mark T Bidwell, Javier Sánchez-fontenla, Andy J Green

<figure><img src="image_1.png" /></figure>

## To cite this version:

Dagmar Frisch, Christophe D Lejeusne, Masaki Hayashi, Mark T Bidwell, Javier Sánchez-fontenla, et al.. Brine chemistry matters: isolation by environment and by distance explain population genetic structure of Artemia franciscana in saline lakes. Freshwater Biology, 2021, 66 (8), pp.1546-1559. 10.1111/fwb.13737. hal-03249374

HAL Id: hal-03249374

https://amu.hal.science/hal-03249374

Submitted on 4 Jun 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Brine chemistry matters: isolation by environment and by distance explain population genetic structure of *Artemia franciscana* in saline lakes

2 population genetic structure of *Artemia franciscana* in saline lakes

3

4 Dagmar Frisch¹,*

5 Christophe Lejeusne²,³

6 Masaki Hayashi⁴

7 Mark T. Bidwell⁵

8 Javier Sánchez-Fontenla³

9 Andy J. Green³

10

11 ¹University of Birmingham, School of Biosciences, Birmingham UK

12 ²Aix Marseille Univ, CNRS, IRD, Avignon Université, IMBE, UMR 7263, Station Marine

13 d'Endoume, Rue de la Batterie des Lions, 13007 Marseille, France

14 ³Department of Wetland Ecology, Estación Biológica de Doñana (EBD-CSIC), C/Américo

15 Vespucio 26, 41092 Sevilla, Spain

16 ⁴Department of Geoscience, University of Calgary, 2500 University Drive NW, Calgary,

17 Alberta, T2N 1N4, Canada

18 ⁵Canadian Wildlife Service, Environment and Climate Change Canada, Saskatoon, SK,

19 Canada

20

21 * corresponding author, email: dfrisch@ou.edu

22

23 Running title: Population genetics of native *Artemia*

24
25 Keywords: brine shrimps, native region, hydrogeochemistry, colonization route, dispersal
26 **Abstract**

27 1. The American brine shrimp *Artemia franciscana* is important in aquaculture and has
28 become invasive in other continents, aided by dispersal via waterbirds. However, little
29 is known about processes underlying its genetic diversity and population structure in its
30 natural habitat in North America. These processes, including dispersal and local
31 adaptation, are pivotal drivers of species distribution and community structure, and
32 therefore central to aquatic biodiversity.

33 2. We studied 15 populations in natural saline lakes of Saskatchewan, Canada to
34 determine the influence of variation in geological history, water chemistry, lake size and
35 location. We aimed to determine the relative importance of isolation by distance (IBD)
36 and isolation by environment (IBE) using the cytochrome c oxidase subunit 1 gene
37 (CO1) as a mitochondrial marker and five nuclear microsatellite markers.

38 3. Geographic patterns for CO1 and microsatellites differed, with lakes clustering in
39 different groups based on genetic distances according to the marker used. CO1 better
40 indicated historical colonization processes, suggesting potential routes of initial
41 colonization when lakes were formed after deglaciation 11,000-15,000 years ago.

42 4. Differentiation between lakes based on nuclear markers was strongly related to
43 variation in hydrochemistry, suggested by distance-based redundancy analysis
44 (dbRDA), but there was no indication of IBD. The ratio between alkalinity and the sum of
45 Ca and Mg concentrations was particularly important, although a lake with a high Cl
46 concentration caused by Potash mining also had a unique *Artemia* population.

47 5. Geochemistry is important in the adaptive radiation of anostracan crustaceans. Our
48 study suggests it also underlies intraspecific genetic variation between populations,
49 promoting IBE, and making dispersal ineffective when cysts are moved by birds
50 between lakes with different hydrochemistry.
51 **Introduction**
52
53 The role of dispersal and environmental filtering in structuring communities and
54 populations is central to research that seeks a better understanding of species
distributions and biogeographical patterns (D'Amen, Mod, Gotelli, & Guisan, 2018;
55 Ovaskainen et al., 2017). Dispersal connects populations across spatial distances,
56 leading to gene flow if followed by successful establishment in the new habitat
57 (Bohonak, 1999). Such successful establishment may be hindered by unsuitable local
58 environmental conditions or by a resident population with a numerical or fitness
59 advantage resulting from priority effects, or local adaptation (De Meester et al 2016).
60
61 Several patterns of genetic differentiation are predicted to result from three main
processes or their interaction: isolation by distance (IBD (Wright, 1943)), isolation by
62 environment (IBE (Wang & Bradburd, 2014) and isolation by adaptation (IBA (Nosil,
63 Egan, & Funk, 2008)). The processes underlying these spatial genetic patterns can be
64 identified by testing relationships between genetic differentiation and geographic and
65 ecological distances (Orsini, Vanoverbeke, Swillen, Mergeay, & De Meester, 2013).
66
67 Here, we examine drivers of biogeographic patterns of the brine shrimp *Artemia
franciscana* Kellog, 1906 (Branchiopoda: Anostraca) in its native range of North
68 America. This species is highly invasive in other continents, to which it was originally
69 exported from the Great Salt Lake and San Francisco Bay as resistant eggs (cysts) used
70 widely as food in aquaculture (Amat et al., 2005; Horváth et al., 2018). It has become a
71 model species in toxicology (Ruebhart et al., 2008; Libralato et al., 2016). In addition, it
72 is also probably the best known example of an aquatic invertebrate that disperses via
73 migratory waterbirds, which consume *Artemia* and disperse their cysts via gut passage
74 owing to their ability to resist digestion (Muñoz, Amat, Green, Figuerola, & Gómez, 2013;
75 Reynolds, Miranda, & Cumming, 2015; Sánchez, Hortas, Figuerola, & Green, 2012).
We focus our study on saline lakes of the Saskatchewan prairies, within the Great Plains bioregion for anostracans (Rogers, 2014a). This area was covered by an ice sheet during the last glacial period until the onset of deglaciation of southern Saskatchewan about 17,000 years BP, which gradually retreated in a north-east direction (Christiansen, 1979).

The export of cysts from San Francisco Bay and the Great Salt Lake in the USA has largely determined the genetic structure of *A. franciscana* in Mexico, as well as across continents other than North America (Eimanifar, Van Stappen, Marden, & Wink, 2014; Horváth et al., 2018). However, populations in central Canada, where the lakes freeze in winter and the hydrochemistry is markedly different, are unaffected by this export (Muñoz et al., 2013). Population genetic studies of this and other *Artemia* species have previously been concentrated in solar salt works or aquaculture ponds (Muñoz et al., 2008; Muñoz et al., 2013, 2014), and studies in natural lakes are rare. Naihong et al. (2000) found greater allozyme heterozygosity for *A. sinica* in larger lakes. Panmixia has been reported for *A. franciscana* in Great Salt Lake (Eimanifar, Marden, Braun, & Wink, 2015) and for *A. urmiana* in Lake Urmia (Eimanifar & Wink, 2013).

Within North America, the distribution of different species of anostracan crustaceans is strongly related to regional differences in geochemistry, and hence hydrochemistry (Rogers, 2014a). This suggests that variation in geo- and hydrochemistry may also be a strong selective force for local adaptation amongst widely distributed species such as *A. franciscana*, and may potentially result in population genetic patterns that suggest IBE.

This is supported by differences observed under laboratory conditions between *A. franciscana* populations in their tolerance of different ionic concentrations (Bowen, Buoncristiani, & Carl, 1988). Saskatchewan has an estimated 500 saline lakes greater than one square kilometer in area, mainly dominated by sodium sulfate and magnesium
101 sulfate salts (Hammer, 1993). Beginning in 1918 in Muskiki Lake, sodium sulfate has
102 been industrially extracted from 20 lakes, which supplied about 50% of the North
103 American demand (Last & Ginn, 2005). However, the hydrochemistry varies
104 considerably between different saline lakes in this region (Bowman & Sachs, 2008).

105 Many of these lakes are important for the conservation of migratory waterbirds,
106 including two of our study sites: Chaplin Lake and Old Wives / Frederick Lakes
107 (Important Bird Areas, https://www.ibacanada.ca/).

108 We would expect that *Artemia* populations colonized the natural saline lakes of
109 Saskatchewan following glacial retreat 11,000-15,000 years ago (Christiansen, 1979),
110 with migratory waterbirds facilitating long-distance dispersal and hence colonization of
111 newly formed lakes (Muñoz et al., 2013). Contemporary gene flow via bird movements
112 may be expected to increase the similarity between lakes, particularly those that are
113 close together, which may be visited repeatedly by the same individual birds even on a
114 daily basis (Demers, Colwell, Takekawa, & Ackerman, 2008; Green, Figuerola, & Sánchez,
115 2002) and thereby increase the effect of IBD. Larger lakes are also likely to hold greater
116 abundance and diversity of waterbirds (Bidwell, Green, & Clark, 2014; Sebastián-
117 González & Green, 2014), increasing the arrival and diversifying the origin of cysts
118 brought from other lakes, and hence promoting a correlation between lake size and
119 genetic diversity. On the other hand, differences in hydrochemistry and resulting local
120 adaptation may mean that only a subset of lakes provide suitable habitat for a given
121 genotype, promoting IBE and reducing the chances that *Artemia* brought by birds as
cysts may survive.

123 For this study, we apply the mitochondrial CO1 gene and microsatellite markers to study
124 historic and current drivers of *Artemia franciscana* biogeography in its native
125 distribution. In particular, we test the hypotheses that (1) genetic diversity is positively
126 correlated with lake area, and that (2) spatial genetic structure of populations in the
127 Saskatchewan saline lakes follows a pattern of isolation by environment (IBE) rather
128 than a pattern of isolation by distance that would result from dispersal limitation (IBD).

129 **Methods**

130 **Study area**

131 From June to August 2011, 50 saline lakes were visited within the Saskatchewan plains
132 and searched for *Artemia*, although many of these lakes were too low in salinity to
133 support *Artemia franciscana*, and some of them only held other anostracans associated
134 with lower salinities (*Branchinecta* spp.). These lakes were partly selected on the basis
135 of existing literature on saline lakes, as some of them were included in broader studies
136 of the chemistry and invertebrate communities of saline lakes in Saskatchewan (Bentley,
137 Hayashi, Zimmerman, Holmden, & Kelley, 2016; Bowman & Sachs, 2008; Hammer, 1993;
138 Last & Ginn, 2005; Wissel, Cooper, Leavitt, & Pham, 2011). A minority of our study lakes
139 are well known for their *Artemia* populations, which have been used in previous
140 laboratory studies (e.g. Chaplin and Little Manitou; Bowen et al., 1988). Additional lakes
141 were sampled because they were encountered during ground surveys of Piping Plover
142 (*Charadrius melodus*) and other waterbirds for a separate study. Finally, other lakes
143 were visited because their saline nature was evident from inspection on Google Earth
144 (https://earth.google.com), with distinctive white soils around the water edge
145 indicating crystallization of saturated salts. Many of our study lakes are previously
146 undocumented, despite previous inventories of the worldwide distribution of *Artemia*
147 (Muñoz & Pacios, 2010). Our sampling was not exhaustive and additional unknown
148 populations of *Artemia* are likely to occur in Saskatchewan. Saline lakes of the Canadian
149 prairies are unaffected by commercial exploitation and export of *A. franciscana* cysts,
150 and the harsh climate makes aquaculture unviable, although adult *Artemia* have been
151 harvested in Chaplin Lake and some other lakes (Hammer, 1986).

152

153 **Sampling**

154 Shorelines were visited during daylight by public access roads and the water column
155 was inspected for brine shrimps, especially on the side where winds concentrated them
156 along the shoreline by wave action. Where present, shrimps were collected with a hand
157 net then placed in alcohol for later confirmation of species identification. At the same
158 time, the shoreline was inspected for concentrations of cysts that form masses at the
159 water's edge or at the high water mark. These cysts were collected in 200 ml plastic
160 bottles, then later washed in distilled water before desiccation for better storage. A
161 water sample was also collected close to the shoreline in 200 ml plastic bottles which
162 were stored in a refrigerator for 1-6 weeks before they were analyzed for major ion
163 chemistry. Lake surface area was estimated from Google Earth (2012). Waterbirds were
164 present in all lakes at the time of visit (details available on request). Lakes used in this
165 study were separated from each other by up to 413 km.

166

167 **Water chemistry analysis**

168 Alkalinity (Alk) was determined by sulfuric-acid titration, and concentrations of other
169 major ions were measured by ion-exchange chromatography after diluting brine
170 samples to a suitable concentration range.

171 Based on these measured parameters, we selected key descriptors of the variation
172 between lakes based on process-based understanding of hydrology and geochemistry of
173 saline lakes in the region (e.g., Bentley et al.(2016); Last & Ginn (2005)). These
174 descriptors were SO₄ concentration, and the ratios Cl/SO₄, Alk/SO₄, and Alk/(Ca +Mg).
175 Lake salinity is primarily controlled by the degree of evaporative enrichment of sulfate,
176 which is sourced by the weathering of sulfide minerals (e.g. pyrite) contained in glacial
177 till (Nachshon, Ireson, van der Kamp, & Wheater, 2013; Van Stempvoort, Hendry,
178 Schoenau, & Krouse, 1994). Evaporative enrichment is greater in closed-basin lakes with
179 higher water inputs by surface runoff and groundwater discharge (Wood & Sanford,
180 1990). Sulfate concentration is a useful indicator of enrichment to a point where the
181 saturation of mirabilite (Na₂SO₄·10H₂O) is reached. Further enrichment increases
182 relative concentration of Mg compared to other cations. If the original composition of
183 water has relatively high alkalinity, the geochemical evolution of water takes a path that
184 is on the alkalinity side of the chemical divide (Hardie & Eugster, 1970). In contrast,
185 originally alkalinity-poor water will take a path on the other side of the chemical divide
186 characterized by high alkaline earth concentrations. Therefore, alkalinity/(Ca + Mg)
187 ratio and alkalinity/SO₄ ratio may be useful for describing the geochemical evolution of
188 saline lakes. Some saline lakes in Saskatchewan (including one in our study, PAT)
189 receive brines from potash mines that contain high concentrations of chloride.
190 Therefore Cl/SO₄ ratio is useful for distinguishing the effects of artificial inputs of
191 chloride.

192

193 **mtDNA genotyping**

194 DNA was extracted from cysts sampled from 12 lakes (Table S1) in the study area using the
195 Hotshot protocol (Montero-Pau, Gomez, & Muñoz, 2008). Two additional populations from
196 two British Columbia lakes sampled in 2011 that belong to a lineage from the other side of
197 the Continental Divide that includes populations of the Great Salt Lake, Utah, and San
198 Francisco Bay (BCO and PCB, Table S1, Fig. S1) were included as outgroup in the haplotype
199 network construction. We used *Artemia*-specific primers 1/2CO1_Fol-F and 1/2CO1_Fol-R
200 (Muñoz et al., 2008) to amplify a 709 bp fragment of the mitochondrial cytochrome c oxidase subunit 1 gene (CO1). Amplifications were performed in 20 µL total volume 201 containing 1X reaction buffer, 2.0mM MgCl2, 0.2mM dNTPs, 0.6 units Taq DNA 202 polymerase (Bioline Corp.) and 0.5µM of each primer. PCR conditions were as follows: 203 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 204 locus), and 60 s at 72°C, followed by 5 min at 72 °C. Quality of sequences was checked 205 manually and consensus sequences were obtained and aligned in BioEdit v. 7.1.9 (Hall, 206 1999). The total number of newly obtained sequences varied between 16 and 24 per lake for 207 12 lakes. All corresponding haplotypes were deposited in NCBI GenBank (Accession 208 numbers: MW799827-MW799878). 209

210

211 **Nuclear microsatellite (nSSR) genotyping**

212 DNA extractions were carried out on individual cysts (previously rinsed in distilled 213 water) from each population making use of the HotSHOT protocol (Montero-Pau et al., 214 2008). The cysts were collected from *Artemia franciscana* populations of 14 lakes in 215 Saskatchewan, Canada (Table 1, Table S1). We amplified eight nuclear microsatellite loci 216 (in the following nSSR or simple sequence repeats) using the primer pairs described in 217 Muñoz, Green, Figuerola, Amat, & Rico (2009) (Af_B109, Af_A108, Af_B9, Af_B139, 218 Af_B117, Af_A136, Af_A104). The locus Af_B11 (Muñoz et al., 2009) was initially 219 included but was dropped because its amplification failed in > 200 individuals. The 220 results from the genotype accumulation curve computed with the R package *poppr* 2.8.3 221 (Kamvar, Tabima, & Grünwald, 2014) suggest that a full resolution of the multi-locus 222 genotypes (MGLs) was already achieved by including only five loci (Fig. S2). We thus 223 decided to use a higher number of genotyped individuals with the trade-off of removing 224 locus Af_B11, yielding data from 592 individuals in total. Following tests for Hardy-
225 Weinberg equilibrium (HWE) (see below) we excluded two additional loci (Af_B117, Af_A136) which deviated from HWE in all or almost all studied populations (Table S3).
226
227 Therefore, throughout the paper, we report results for the five nSSR loci Af_B109, Af_A108, Af_B9, Af_B139, Af_A104). All included loci were polymorphic and their allelic richness was between 21 and 48 alleles (Table S4).
228
229
230

231 **Statistical analysis**

232 All R packages used for data analyses described in this section were run in R version 3.6.2 (R Core Team, 2019).
233
234

235 *Environmental variables.* In order to quantify and visualize how the environment available to *Artemia franciscana* differed between individual lakes, Principle Component Analysis (PCA) was performed with R package *vegan* 2.5.6 (Oksanen et al., 2019) and plotted with the autoplot() function of *ggplot2* (Wickham, 2016). The PCA included the following variables: lake surface area, SO₄, and the ratios Cl/SO₄, Alk/SO₄, and Alk/(Ca +Mg). 236
237
238
239
240

241 *mtDNA.* To analyze mtDNA, we supplemented the haplotype dataset newly sequenced in this study, with haplotypes that were previously reported from four of our study lakes (MAN, CHA, MEA and MUS), based on earlier sampling in 2009 (Muñoz et al., 2013). The entire dataset used for the current study therefore included 65 haplotypes in total: 52 haplotypes newly detected in the present study, three haplotypes recorded both in the present study and in Muñoz et al. (2013) and 10 additional haplotypes only reported in Muñoz et al. (2013). Details on sample size and haplotypes can be found in Table S1 in supplementary material. 242
243
244
245
246
247
248
A haplotype network was generated using TCS 1.21 (Clement, Posada, & Crandall, 2000) at the 95% connection limit. This network includes 12 lake populations from Saskatchewan and two lake populations from British Columbia. Genetic diversity within lakes was characterized by the number of private haplotypes (haplotypes found only in one location), standard diversity indices of haplotype diversity *h* and nucleotide diversity *π* (Nei 1987) calculated in Arlequin v. 3.5 (Excoffier & Lischer, 2010) and DnaSP v.5 (Librado & Rozas, 2009). Pairwise *Φ*<sub>ST</sub> with Tamura & Nei (TrN) distances and their statistical significance (10,000 permutations) were performed in Arlequin v. 3.5.

An analysis of molecular variance (AMOVA) was conducted to determine the hierarchical structure of CO1 variation (Excoffier, Smouse, & Quattro, 1992). This analysis was performed with the aforementioned substitution model with 10,000 random permutations in Arlequin and based on the groups of populations that were not significantly genetically differentiated according to pairwise *Φ*<sub>ST</sub> (Group 1: BEN, BES, WHI, FRE, CHA; Group 2: AKE; Group 3: EIN, LYD, MAN, MUS, PAT, MEA). We tested three hierarchical levels: between groups, between populations within groups, within populations.

*Microsatellite (nSSR) markers.* To estimate null allele frequencies, we used the software FreeNA (Chapuis & Estoup 2007) that follows the Expectation Maximization (EM) algorithm of Dempster, Laird, and Rubin (1977) (detailed in Chapuis & Estoup 2007).

Hardy-Weinberg equilibrium was tested with the function hw.test(), computing an exact test based on 5000 Monte Carlo permutations of alleles. This function is implemented in the R package *pegas* 0.14 (Paradis 2010). Results were corrected for multiple comparisons using the method of Holm (1979), a less conservative modification of the Bonferroni test. Total allelic richness per populations was computed by rarefaction implemented in the R package *PopGenReport* 3.0.4 (Adamack & Gruber 2014). The population with the smallest number of individuals genotyped for nSSRs (SHO, 8
274 individuals) was excluded from the rarefaction procedures and the second smallest

275 population size (GUL, n = 21 individuals) applied as the smallest size for rarefaction.

276 Multilocus genotype diversity H (Shannon-Wiener diversity (Shannon, 1948)) was

277 computed with the R package *poppr* 2.8.3 (Kamvar et al., 2014). Locus- and population-

278 specific allele number and number of private alleles per population was calculated with

279 *poppr* 2.8.3 (Kamvar et al., 2014).

280 For pairwise population genetic distances, we computed Weir's (1996) *F*<sub>ST</sub> using FreeNA

281 (Chapuis & Estoup 2007). This software provides a correction of the bias induced by null

282 alleles following the method described in Chapuis & Estoup (2007). We calculated

283 pairwise population *D*<sub>EST</sub> with GenAlEx 6.502 (Peakall & Smouse, 2012) with associated

284 p-values based on 999 permutations.

285 To further explore the population structure suggested by pairwise population distance,

286 we applied discriminant analysis of principal components (DAPC) implemented in the

287 package adegenet 2.1.3 (Jombart 2008), a model-free analysis that makes no prior

288 assumptions about population genetic processes (such as HWE). We used the function

289 *xvalDAPC()* to find the optimum number of PCs to be retained prior to their use in the

290 final discriminant analysis, which we performed with the first 40 PCs and the first five

291 discriminant functions. We specified the 14 study populations as *a priori* groups to

292 analyse population membership probabilities for each individual.

293 *Distance-based redundancy analysis (dbRDA)*: To test for the importance of isolation by

294 environment (IBE) vs. isolation by distance (IBD), we performed dbRDA with functions

295 implemented in the R package *vegan* 2.5.6 (Oksanen et al., 2019).

296 We used genetic distance estimators (pairwise *D*<sub>EST</sub> for nSSR markers and pairwise *Φ*<sub>ST</sub>

297 for mtDNA) as dependent variables. To test the effect of individual environmental

298 variables (the same ones as used in the above PCA) on population differentiation, we
299 built the null model (genetic differentiation ~ 1) and the full model (genetic differentiation ~ Area + SO₄ + Cl/SO₄ + Alk/SO₄ + Alk/Ca+Mg) using the function capscale(). Both the null model and the full model were used as input to the forward selection model choice (on adjusted R² and p-values) with the function ordiR2step() with 999 permutations. Model choice was repeated 100 times to obtain robust results.
300
301
302
303
304 To test the effect of spatial configuration, we computed distance-based Moran's Eigenvector Maps (dbMEMs, previously known as PCNMs) using the R package adespatial 0.3.7 (Dray et al., 2019) and chose the dbMEMs with associated Moran's I values larger than expected under the null hypothesis of no spatial autocorrelation (here: the first three dbMEMs: MEM1, MEM2, MEM3). For dbRDA, the null model was tested against the full spatial model (genetic differentiation ~ MEM1 + MEM2 + MEM3) using the forward selection model choice as described above.
305
306
307
308
309
310
311 *Network analysis.* To visualise the connectivity (as a proxy for gene flow) between populations based on nSSR loci, we constructed a network graph using EDENnetworks version 2.18 (Kivelä, Arnaud-Haond, & Saramäki, 2015) with manual thresholding, based on FST corrected for null alleles obtained by FreeNA as described above. Manual thresholding was performed at different values (0.02, 0.03 and 0.05) for visualisation of possible population clusters under different threshold scenarios. Analysis of molecular variance (AMOVA) for nSSR markers was used to test the significance of the population structure suggested by the lowest threshold (0.02). Three hierarchical levels were tested: (1) between 6 individual populations and two multi-population groups: AKE (group1), WHI (group2), MUS (group3), MAN (group4), PAT (group5), FRE (group6), BEN + BES (group7) and EIN + LYD + CHA + GUL + FRS + SHO (group8), (2) between populations within groups, and (3) within populations. AMOVA was calculated using the
323 R package *ade4* 1.7-13 implemented in *poppr* 2.8.3 (Kamvar et al., 2014). Statistical
324 significance of variance components was tested by permutation tests (999
325 permutations).

326

327

328 **Results**

329 *Study sites – spatial distribution, lake size and hydrochemical properties*

330 The 15 Saskatchewan saline lakes included in the genetic analysis spanned over a range
331 of 2.63° in latitude, 3.93° in longitude, and were separated by distances of up to 410 km.

332 Their spatial distribution was represented as Moran Eigenvectors for further analysis
333 (Fig. S3). Their surface area ranged from 0.14 to 29.6 km² (Table 1). We obtained

334 hydrochemical data for 13 of the lakes, including those for which nSSR data are available
335 (with the exception of FRE and SHO, Table S2). The total dissolved solids of the 13 lake
336 water samples ranged from 13 to 207 g L⁻¹, which were roughly proportional to SO₄
337 concentrations, except for PAT (see below).

338 A Principal Component Analysis (PCA) including SO₄ meq, the ratios Cl/SO₄, Alk/SO₄,
339 Alk/(Ca+Mg) and lake area explained 68.9% of the environmental variation between
340 lakes between the first two components (Fig. 1). Two lakes were clearly separated from
341 the other sites: PAT was the site of a Potash mine and was associated with high values of
342 Cl/SO₄, whereas AKE had particularly high Alk/(Ca+Mg). GUL was also hydrochemically
343 separated from other lakes by a positive association with Alk/SO₄. There was a tendency
344 for greater lake area to be associated with higher SO₄ (Fig. 1): MUS and WHI were
345 among the largest lakes, and also had high SO₄.

346
347 *Haplotype and allelic diversity*

348 We identified a total of 52 new haplotypes in the sampled lakes. Of these, seven were
349 unique to the British Columbia locations. Three additional haplotypes were also
350 reported in a previous study, that also recorded 10 additional haplotypes not found in
351 our study (Muñoz et al., 2013). To construct the haplotype network and for all other
352 mtDNA analyses below, we included haplotype data from these authors (details in Table
353 S6).

354 Three star-like CO1 haplogroups formed by a total of 65 haplotypes exist in the
355 Saskatchewan area, each with several additional haplotypes of low frequency (Fig. 2).

356 These three main groups were centered on two haplotypes also recorded by Muñoz et
357 al. (2013) (Af_35, AF_41, Group A and C in Fig. 2) and a newly detected haplotype
358 (Hap_05, Group B). Overall, separation from the British Columbia haplotypes involved
359 about 4-8 mutations, while haplotypes within Saskatchewan were not strongly
360 differentiated (the two main haplotypes differed only by one mutation).

361 Diversity patterns differed between CO1 and nSSR markers (Table 1): CO1 haplotype
362 diversity ranged from 0.000 (PAT) to 0.779 (AKE). For nSSR, rarefied allelic richness
363 was between 26.3 (PAT) and 54.9 (FRS). Diversity of multilocus genotypes H ranged
364 from 2.493 (GUL) to 3.493 (EIN). No linear relationships were detected between allelic
365 or haplotype diversity and lake area (Table S7).

366

367 *Population genetic structure in geographic and environmental space*

368 There was significant structure between groups of populations according to the
369 mitochondrial CO1 gene (Table 2) with two main groups (BEN, BES, FRE, WHI, CHA, and
370 LYD, EIN, MAN, MUS, MEA, PAT), each of which was dominated by one of the two main
371 haplotypes. Differentiation between members of these two groups was strong, with
372 pairwise ΦST values between 0.548 and 0.891. AKE (with the greatest haplotype diversity) was significantly differentiated from all other populations. An AMOVA of mtDNA based on the above mentioned three groups confirmed that the variation was highest between groups (68%), and much smaller between populations within groups (1.37%). Within population variation was 31% (Table 3a).
373
374
375
376
377 While the pattern of differentiation was different from that of mitochondrial markers, we also detected significant genetic differentiation between populations according to pairwise DEST based on nSSR loci (Table 2), ranged from -0.007 (between BEN and BES) and 0.418 (between MUS and AKE).
378
379
380
381 Network visualisation based on FST corrected for null alleles (Table S5) revealed a topology with two network components at a threshold of 0.02 (Fig. 3a). One included EIN, LYD, GUL, CHA, FRS and SHO (some of which were well separated in space), while the second included BES and BEN (which were closer together). The other five populations remained disconnected at this threshold. A gradual increase of the threshold led to higher complexity of the networks (Fig. 4b and c) with increasing connectivity between the two clusters identified in Fig. 4a. However, AKE, WHI, PAT, MAN and FRE remained separated even at the highest threshold. Using the two network components identified at a threshold of 0.05, we found that genetic diversity (rarefied allelic richness of nSSR loci) differed between connected and disconnected nodes (1-way ANOVA: F(1,11) = 17.3, p = 0.002) (Fig. 4d). This analysis suggested that populations within a network component were genetically more diverse (median = 50.44, SD = 4.41) than separated populations (median = 40.54, SD = 6.46).
381
382
383
384
385
386
387
388
389
390
391
392
393
394 The population groups identified from the network analysis at a threshold of 0.02 (Fig. 3a) were tested by AMOVA (Table 3b). The results were significant for all strata, but by far the highest amount of variation was found within individual populations (83.50%).
395
396
397 Variation between groups of populations was 14.09%, suggesting a different pattern
398 than that seen in mtDNA, as also suggested by the differences in group composition (e.g.
399 the nearby BEN and BES group together with the more distant CHA and FRE by CO1, yet
400 they form a separate group according to nSSR). Genetic variation of nSSR was lowest
401 between populations within groups (6.09%).

402 To analyse environmental conditions and spatial configuration of sample sites that
403 might underlie the observed genetic structure, we performed dbRDA using pairwise
404 genetic distance (Φ<sub>ST</sub> for the mitochondrial CO1 data, D<sub>EST</sub> for nSSR loci, Table 2) as
405 independent variables. No significant relationship was detected for mtDNA for either
406 spatial or environmental variables. In contrast, for nSSR markers the forward selection
407 model choice found strong and significant relationships with the ratio Alkalinity/(Ca
408 +Mg) but not with spatial configuration. Alkalinity/(Ca +Mg) was significant in 87 of
409 100 repeats of the model choice (p ≤ 0.05), explaining 30% of the observed genetic
410 variation with a mean p-value of 0.04 (Table 4).

411 A discriminant analysis of principal components (DAPC, Fig. 4) of nSSR markers
412 revealed a population genetic structure with six comparatively well-defined populations
413 in which the majority of individuals were assigned to their respective population of
414 origin (AKE, WHI, PAT, MUS, FRE and MAN), echoing the unconnected network nodes in
415 Fig. 3. A considerable amount of admixture is suggested for other populations where
416 EIN, LYD, GUL, CHA and FRS had higher membership probabilities for one of the
417 westernmost populations while MUS, BEN and BES showed a high degree of admixture
418 among themselves.

419

420 Discussion
As an important species for aquaculture, *Artemia franciscana* occurs widely outside its native distribution in North America. It is highly invasive, rapidly outcompeting native *Artemia* species. Invasive populations generally originate from commercially harvested cysts from the Great Salt Lake and San Francisco Bay, or secondary introductions from China (Amat et al., 2005; Horváth et al., 2018). In contrast, this species has rarely been studied in other areas of its native distribution, where it inhabits natural saline lakes with a wide range of hydrochemistries in particular regarding salt composition (Bowman & Sachs, 2008). In order to identify the roles of IBE and IBD in driving spatial genetic structure of *A. franciscana* in the saline lakes of Central Canada, we applied mitochondrial and nuclear markers to 14 lake populations in the Saskatchewan prairies.

We found significant genetic structure between the Saskatchewan populations in both marker systems, and a strong relationship between the divergence pattern observed by nSSR markers and the hydrochemical environment of the lakes.

The mtDNA patterns observed are likely to reflect historical patterns of lake colonization, especially when lakes were formed after deglaciation, and then became saline enough to provide suitable habitat for *A. franciscana*, a species associated with higher salinities than all other anostracans in North America (Rogers, 2014a, 2014b). In the unique case of Patience (PAT), a lake contaminated with chloride by potash mining, colonization may have occurred in modern times after mining began. The presence of only one haplotype in PAT, combined with its unique pattern of nuclear markers, suggests strong and rapid local adaptation to these unusual chemical conditions absent from our other study sites, which were lakes in a more natural condition.

*What mitochondrial markers reveal*
For lakes other than PAT, our results suggest a close relationship of colonization by *Artemia* and deglaciation after the last glacial period, which first opened up habitat alongside the southwestern border of the ice sheet (14,000y BP). Under such a scenario, the area may first have been colonized in the western area by individuals originating from ice-free saline lake populations further to the west or south, such as Alberta or the USA. From the four western-most Saskatchewan populations (AKE, EIN, WHI, LYD), which hold all three haplotype groups identified in this study, dispersal may have proceeded towards the newly opened south-eastern habitats, subsequently colonized by one of the main haplotype groups (group C). Later, as the glacier retreated northeastwards and aquatic habitat opened towards the east in Saskatchewan, these northern lakes were largely colonized by the second main haplotype group (group A).

Initial colonization by the two main haplotypes was followed by haplotype diversification, consistent with the observed star-like network pattern. After the area was ice-free in the next 1000 years (Christiansen, 1979), additional dispersal from the southern and western areas may have allowed colonization of these postglacial lakes.

While this proposed sequence for colonization is consistent with the haplotype pattern observed in this study, it would require further study by including more *A. franciscana* populations from additional lakes in Saskatchewan and to the west in Alberta. Future sampling should also include populations further south in the Great Plains bioregion of the USA (Rogers, 2014b) to test the possibility that dispersal towards the northern areas may have followed the retreating ice sheet, and originated from regions south of Saskatchewan.

*What nuclear markers reveal*
469 In contrast, nSSR markers showed significant patterns of Isolation by Environment (IBE), explaining ~30% of the variation in the genetic structure, but unrelated to geographic configuration of the sampled lakes. The observed patterns are indicative of local adaptation to the environment according to the criteria detailed in Orsini et al. (2013). Since *Artemia* produces vast egg banks in established populations, ecological and genetic priority effects occurring soon after initial colonisation and subsequent monopolisation effects (De Meester, Gómez, Okamura, & Schwenk, 2002, De Meester, Vanoverbeke, Kilsdonk, & Urban, 2016) may have reinforced such a pattern, including the lack of a spatial signal that could be related to such priority effects. Our evidence for IBE suggests that bird-mediated dispersal of cysts between lakes differing markedly in hydrochemistry may not be effective (i.e. will not lead to colonization and establishment in the new location). Laboratory tests show that nauplii larvae fail to reach adulthood when placed in brine of unsuitable ionic composition (Bowen et al., 1988), and that cyst hatching rates are sensitive to metal concentrations (Brix et al., 2006).

483 In several cases, strong differentiation is highly consistent with local adaptation to extreme hydrochemical conditions, notably in PAT and AKE, but also in WHI, which had the highest SO₄ concentrations in our study. In contrast to the mitochondrial data, nSSR markers revealed that PAT was a highly divergent population, strongly related to an elevated Cl/SO₄ ratio resulting from potash mining. Similarly, in AKE, this was related to a particularly high Alkalinity/(Ca+Mg) ratio. Adaptation of different *A. franciscana* populations to specific ion concentrations was demonstrated in the laboratory by Bowens et al. (1988), leading to ecological isolation of individual populations and an inability to colonize other lakes with divergent chemistry. This is also consistent with the apparent influence of geochemistry on the evolution and distribution of different anostracan taxa in North America (Rogers, 2014a). Since the nSSR markers are
494 essentially neutral, to the extent that lakes have similar genetic diversity due to similar hydrochemistry, this cannot simply reflect convergent selection on genes for osmoregulation or other key physiological processes. Instead, the genetic structure is likely to reflect contemporary gene flow, which is inhibited by divergent chemistry but promoted by a close spatial relationship between lakes.

499 Network analysis revealed a high connectivity between many populations, reflecting their genetic similarity for nuclear markers. Together with PAT and AKE, a strongly isolated population was FRE, which is surprising since this population is adjacent to FRS and only isolated from it by a dam which also serves as a road. This creates two of several separate sub-basins in Frederick Lake used for the extraction of sodium sulfate via evaporation. The first basin (FRE) discharges water to FRS, which is a closed basin with no surface water outlet. As a result, the water in FRE is expected to be less saline compared to the concentrated brine in FRS, although we have no chemistry data for FRE.

507 Elsewhere, BES and BEN were also separated by a road, but were not differentiated genetically, and had rather similar hydrochemistry, as they were parts of a larger lake connected by culverts under the road.

510 Other lakes that were relatively isolated according to network analysis were WHI (which was characterised by the highest SO₄ concentrations) and MAN. It is possible that connectivity (i.e. genetic similarity) is influenced by the extent of interchange of migratory birds between lakes. The most important site for waterbirds is CHA, which is protected as a Ramsar site in recognition of the major concentrations of shorebirds it holds. CHA was a central node in the network, connected with more populations in the network than any other lake. On the other hand, adequate data on bird movements between our study lakes do not exist, but our bird counts (AJ Green, unpublished data) would suggest that WHI and MAN should be reasonably well-connected via waterbirds
519 with other lakes. WHI held thousands of Wilson's phalaropes *Phalaropus tricolor* which
520 were feeding on *Artemia* and are likely to be excellent vectors for their dispersal.

521 Hatching of nauplii from faeces collected at Frederick Lake from the American avocet
522 *Recurvirostra americana* during the course of this study confirmed the internal transport
523 of viable *A. franciscana* cysts in our study area (AJG, unpublished data). We did not find
524 evidence that lake area significantly affected genetic diversity. However, network
525 analysis showed that there is greater genetic diversity in highly connected populations
526 than in the less connected populations, and future research could test whether this
527 reflects movement patterns of birds between lakes. Outside of migratory periods,
528 waterbirds will generally provide higher connectivity for sites that are close together,
529 and hence coincide within the home range of more individual birds (Green et al., 2002).

530
531 *Conclusions*
532 Saline lakes of Saskatchewan provide a fascinating insight into the population genetics
533 of *Artemia* under natural conditions. The importance of IBE in this keystone crustacean
534 was only apparent as a consequence of our unusual application of hydrochemical
535 analysis. Studies of the population or invasion genetics of aquatic invertebrates in inland
536 waters typically limit their analyses of water quality to conductivity or salinity, and few
537 if any have gone into as much detail about ionic composition as in our study. More
538 frequent collaborations between hydrogeologists and molecular ecologists are
539 recommended. Although many of the sites occupied by other *Artemia* species in
540 different continents are coastal lagoons and salt ponds that are likely to have little
541 diversity in salt composition, many others are inland salt lakes (Muñoz & Pacios, 2010)
542 that may be hotspots of genetic diversity associated with local adaptation to unique
543 chemical conditions. On the other hand, the *Artemia* populations in such inland lakes
544 (e.g. *A. sinica* and *A. tibetiana* populations, Lin et al., 2017) may be resistant to invasion
545 by *A. franciscana* populations originating from aquaculture, owing to contrasting
546 hydrochemistry.

547

548 Acknowledgements

549 Robert G. Clark hosted AJG and provided vital assistance during sampling. Cysts from BC
550 were kindly provided by Curtis A. Suttle, Amy M. Chan and Allyson L. Brady. Funding
551 was provided by Spanish National Plan projects CGL2010-16028 and CGL2016-76067-P
552 (AEI/FEDER, EU), and the PICS program funded by CSIC (PIC2015FR4) in Spain and
553 CNRS (PICS07360) in France. DF received funding from the European Union's FP6-
554 MOBILITY programme under the Marie Skłodowska Curie grant agreement No. 40285.

555 AJG was also supported by mobility grants from the Spanish Ministry of Science and
556 Innovation (PR2008-0293) and the Spanish Ministry of Education (PR2010-0602).

557 Additional support was provided by Canadian Wildlife Service.

558

559 *Data availability statement*

560 Sequence data that support the findings of this study are available in NCBI GenBank,
561 accession numbers MW799827-MW799878. Hydrochemical data is made available at
562 digital.csic.es.

563
564 **References**

565 Adamack, A.T. & Gruber, B. (2014). PopGenReport: simplifying basic population genetic analyses in R. *Methods in Ecology and Evolution*, 5(4), 384-387.

567 Amat, F., Hontoria, F., Ruiz, O., Green, A. J., Hortas, F., Figuerola, J., & Hortas, F. (2005). The American brine shrimp as an exotic invasive species in the western Mediterranean. *Biological Invasions*, 7(1), 37-47.

570 Bentley, L. R., Hayashi, M., Zimmerman, E. P., Holmden, C., & Kelley, L. I. (2016). Geologically controlled bi-directional exchange of groundwater with a hypersaline lake in the Canadian prairies. *Hydrogeology Journal*, 24(4), 877-892. doi: 10.1007/s10040-016-1368-0

574 Bidwell, M. T., Green, A. J., & Clark, R. G. (2014). Random placement models predict species-area relationships in duck communities despite species aggregation. *Oikos*, 123(12), 1499-1508. doi: 10.1111/oik.00821

577 Bohonak, A. J. (1999). Dispersal, gene flow, and population structure. *Quarterly Review of Biology*, 74(1), 21-45.

579 Bowen, S. T., Buoncristiani, M. R., & Carl, J. R. (1988). *Artemia* habitats: Ion concentrations tolerated by one superspecies. *Hydrobiologia*, 158(1), 201-214. doi: 10.1007/BF00026278

582 Bowman, J. S., & Sachs, J. P. (2008). Chemical and physical properties of some saline lakes in Alberta and Saskatchewan. *Saline Systems*, 4(1), 3. doi: 10.1186/1746-1448-4-3

585 Brix, K. V., Gerdes, R. M., Adams, W. J., & Grosell, M. (2006). Effects of Copper, Cadmium, and Zinc on the Hatching Success of Brine Shrimp (*Artemia franciscana*). *Archives of Environmental Contamination and Toxicology*, 51(4), 580-583.
588 https://doi.org/10.1007/s00244-005-0244-z

589 Chapuis, M.P., & Estoup, A. (2007). Microsatellite null alleles and estimation of population differentiation. *Mol. Biol. Evol.* 24(3): 621-631.

591 Christiansen, E. A. (1979). Wisconsinan deglaciation of southern Saskatchewan and adjacent areas. *Can J Earth Sci*, 16(4), 913-938. doi: 10.1139/e79-079

593 Clement, M., Posada, D., & Crandall, K. A. (2000). TCS: A computer program to estimate gene genealogies. *Molecular Ecology* 9(10), 1657-1659. doi: 10.1046/j.1365-294X.2000.01020.x

596 D'Amen, M., Mod, H. K., Gotelli, N. J., & Guisan, A. (2018). Disentangling biotic interactions, environmental filters, and dispersal limitation as drivers of species co-occurrence. *Ecography*, 41(8), 1233-1244. doi: 10.1111/ecog.03148

599 De Meester, L., Gómez, A., Okamura, B., & Schwenk, K. (2002). The Monopolization Hypothesis and the dispersal-gene flow paradox in aquatic organisms. *Acta Oecologica*, 23, 121-135.

602 De Meester, L., Vanoverbeke, J., Kilsdonk, L. J., & Urban, M. C. (2016). Evolving Perspectives on Monopolization and Priority Effects. *Trends in Ecology and Evolution*, 31(2),136-146. doi: 10.1016/j.tree.2015.12.009

605 Demers, S. A., Colwell, M. A., Takekawa, J. Y., & Ackerman, J. T. (2008). Breeding stage
606 influences space use of female American avocets in San Francisco Bay, California.

607 *Waterbirds*, 31(3), 365-371. doi: 10.1675/1524-4695-31.3.365

608 Dempster, A. P., Laird, N. M. & Rubin, D. B. (1977) Maximum likelihood from incomplete

609 data via the EM algorithm. *Journal of the Royal Statistical Society B* 39:1-38.

610 Dray, S., Bauman, D., Blanchet, G., Borcard, D., Clappe, S., Guenard, G., ... Wagner, H. H.

611 (2019). *adespatial: Multivariate Multiscale Spatial Analysis*. R package version 0.3-7.

613 Eimanifar, A., Marden, B., Braun, M. S., & Wink, M. (2015). Analysis of the genetic

614 variability of *Artemia franciscana* Kellogg, 1906 from the Great Salt Lake (USA)

615 based on mtDNA sequences, ISSR genomic fingerprinting and biometry. *Marine*

616 *Biodiversity*, 45(2), 311-319. doi: 10.1007/s12526-014-0256-x

617 Eimanifar, A., Van Stappen, G., Marden, B., & Wink, M. (2014). *Artemia* biodiversity in

618 Asia with the focus on the phylogeography of the introduced American species

619 *Artemia franciscana* Kellogg, 1906. *Molecular Phylogenetics and Evolution*, 79, 392-

620 403. doi: 10.1016/j.ympev.2014.06.027

621 Eimanifar, A., & Wink, M. (2013). Fine-scale population genetic structure in *Artemia*

622 *urmiana* (Günther, 1890) based on mtDNA sequences and ISSR genomic

623 fingerprinting. *Organisms Diversity and Evolution*, 13(4), 531-543. doi:

624 10.1007/s13127-013-0135-5

625 Excoffier, L, Smouse, P. E., & Quattro, J. M. (1992). Analysis of Molecular Variance

626 Inferred from Metric Distances among DNA Haplotypes - Application to Human

627 Mitochondrial-DNA Restriction Data. *Genetics*, 131(2), 479-491.

628 Excoffier, Laurent, & Lischer, H. E. L. (2010). Arlequin suite ver 3.5: A new series of

629 programs to perform population genetics analyses under Linux and Windows.

630 *Molecular Ecology Resources*, 10(3), 564-567. doi: 10.1111/j.1755-

631 0998.2010.02847.x

632 Green, A. J., Figuerola, J., & Sánchez, M. I. (2002). Implications of waterbird ecology for

633 the dispersal of aquatic organisms. *Acta Oecologica*, 23, 177-189.

634 Hall, T. A. (1999). *BioEdit: a user-friendly biological sequence alignment editor*. *Nucleic*

635 *Acids Symposium Series*.

636 Hammer, U. Theodore. (1986). Saline Lake Resources of the Canadian Prairies. *Canadian*

637 *Water Resources Journal*, 11(1), 43-57. doi: 10.4296/cwrj1101043

638 Hammer, U.T. (1993). Zooplankton distribution and abundance in saline lakes of Alberta

639 and Saskatchewan, Canada. *International Journal for Salt Lake Research*, 2(111-

640 132).

641 Hardie, L., & Eugster, H. (1970). Evolution of closed-basin brines. *Mineral. Soc. Amer.*

642 *Spec. Pap. 3*, 273-290.

643 Holm, S. (1979). A simple sequentially rejective multiple test procedure. *Scandinavian*

644 *Journal of Statistics*, 6, 65-70.

645 Horváth, Z., Lejeusne, C., Amat, F., Sánchez-Fontenla, J., Vad, C. F., & Green, A. J. (2018).

646 Eastern spread of the invasive *Artemia franciscana* in the Mediterranean Basin, with

647 the first record from the Balkan Peninsula. *Hydrobiologia*. doi: 10.1007/s10750-
648 018-3683-z

649 Jombart, T. (2008). adegenet: a R package for the multivariate analysis of genetic
650 markers. *Bioinformatics*, 24, 1403-1405.

651 Kamvar, Z. N., Tabima, J. F., & Grünwald, N. J. (2014). Poppr: an R package for genetic
652 analysis of populations with clonal, partially clonal, and/or sexual reproduction.
653 *PeerJ*, 2, e281. doi: 10.7717/peerj.281

654 Kivelä, M., Arnaud-Haond, S., & Saramäki, J. (2015). EDENetworks: A user-friendly
655 software to build and analyse networks in biogeography, ecology and population
656 genetics. *Molecular Ecology Resources*, 15(1), 117-122. doi: 10.1111/1755-
657 0998.12290

658 Last, W. M., & Ginn, F. M. (2005). Saline systems of the Great Plains of western Canada:
659 an overview of the limnogeology and paleolimnology. *Saline Systems*, 1,10. doi:
660 10.1186/1746-1448-1-10

661 Librado, P., & Rozas, J. (2009). DnaSP v5: A software for comprehensive analysis of DNA
662 polymorphism data. *Bioinformatics*, 25(11), 1451-1452. doi:
663 10.1093/bioinformatics/btp187

664 Libralato, G., E. Prato, L. Migliore, A. M. Cicero, and L. Manfra. 2016. A review of toxicity
665 testing protocols and endpoints with Artemia spp. *Ecological Indicators* **69**:35-49.

666 Lin, Q., Xu, L., Hou, J., Liu, Z., Jeppesen, E., & Han, B.-P. (2017). Responses of trophic
667 structure and zooplankton community to salinity and temperature in Tibetan lakes:
668 Implication for the effect of climate warming. *Water Research*, 124, 618-629.
669 https://doi.org/https://doi.org/10.1016/j.watres.2017.07.078

670 Montero-Pau, J., Gomez, A., & Muñoz, J. (2008). Application of an inexpensive and high-
671 throughput genomic DNA extraction method for the molecular ecology of
672 zooplanktonic diapausing eggs. *Limnology and Oceanography—Methods*, 6, 218-222.

673 Muñoz, J., Amat, F., Green, A. J., Figuerola, J., & Gómez, A. (2013). Bird migratory flyways
674 influence the phylogeography of the invasive brine shrimp *Artemia franciscana* in
675 its native American range. *PeerJ*, 1, e200. doi: 10.7717/peerj.200

676 Muñoz, J., Gómez, A., Figuerola, J., Amat, F., Rico, C., & Green, A. J. (2014). Colonization
677 and dispersal patterns of the invasive American brine shrimp *Artemia franciscana*
678 (Branchiopoda: Anostraca) in the Mediterranean region. *Hydrobiologia*, 726(1), 25-
679 41. doi: 10.1007/s10750-013-1748-6

680 Muñoz, J., Gomez, A., Green, A. J., Figuerola, J., Amat, F., & Rico, C. (2008). Phylogeography
681 and local endemism of the native Mediterranean brine shrimp *Artemia salina*
682 (Branchiopoda : Anostraca). *Molecular Ecology*, 17(13), 3160-3177.

683 Muñoz, J., Green, A. J., Figuerola, J., Amat, F., & Rico, C. (2009). Characterization of
684 polymorphic microsatellite markers in the brine shrimp *Artemia* (Branchiopoda,
685 Anostraca). *Molecular Ecology Resources*, 9(2), 547-550.

686 Muñoz, J., & Pacios, F. (2010). Global biodiversity and geographical distribution of
diapausing aquatic invertebrates: The case of the cosmopolitan brine shrimp,
*Artemia* (Branchiopoda, Anostraca). *Crustaceana*, 83(4), 465-20.A. doi:
689 10.1163/001121610X489449
690 Nachshon, U., Ireson, A., van der Kamp, G., & Wheater, H. (2013). Sulfate salt dynamics in the glaciated plains of North America. *Journal of Hydrology*, 499, 188-199. doi: 10.1016/j.jhydrol.2013.07.001

693 Naihong, X., Audenaert, E., Vanoverbeke, J., Brendonck, L., Sorgeloos, P., & De Meester, L. (2000). Low among-population genetic differentiation in Chinese bisexual *Artemia* populations. *Heredity*, 84(2), 238-243. doi: 10.1046/j.1365-2540.2000.00664.x

696 Nosil, P., Egan, S. P., & Funk, D. J. (2008). Heterogeneous genomic differentiation between walking-stick ecotypes: 'Isolation by adaptation' and multiple roles for divergent selection. *Evolution*, 62(2), 316-336. doi: 10.1111/j.1558-5646.2007.00299.x

700 Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., ... Wagner, H. (2019). vegan: Community Ecology Package. R package version 2.5-6.

702 Orsini, L., Vanoverbeke, J., Swillen, I., Mergeay, J., & De Meester, L. (2013). Drivers of population genetic differentiation in the wild: Isolation by dispersal limitation, isolation by adaptation and isolation by colonization. *Molecular Ecology*, 22(24), 5983-5999. doi: 10.1111/mec.12561

706 Ovaskainen, O., Tikhonov, G., Dunson, D., Grotan, V., Engen, S., Saether, B. E., & Abrego, N. (2017). How are species interactions structured in species-rich communities? A new method for analysing time-series data. *Proc Biol Sci*, 284(1855), 20170768. doi: 10.1098/rspb.2017.0768

710 Paradis, E. (2010). pegas: an R package for population genetics with an integrated-modular approach. *Bioinformatics*, 26, 419-420.

712 Peakall, R., & Smouse, P. E. (2012). GenALEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research-an update. *Bioinformatics*, 28(19), 2537-2539. doi: 10.1093/bioinformatics/bts460

715 R Core Team (2019). (2019). R: A language and environment for statistical computing. *Accessed* 1st April 2019.

717 Reynolds, C., Miranda, N. A. F., & Cumming, G. S. (2015). The role of waterbirds in the dispersal of aquatic alien and invasive species. *Diversity and Distributions*, 21(7), 744-754. doi: 10.1111/ddi.12334

720 Rogers, D. C. (2014a). Anostracan (Crustacea: Branchiopoda) zoogeography II. Relating distribution to geochemical substrate properties in the USA. *Zootaxa*, 3838(3), 1. doi: 10.11646/zootaxa.3856.1.1

723 Rogers, D. C. (2014b). Anostracan (Crustacea: Branchiopoda) zoogeography I. North American bioregions. *Zootaxa*, 3838(3), 251. doi: 10.11646/zootaxa.

725 Ruebhart, D. R., I. E. Cock, and G. R. Shaw. 2008. Brine shrimp bioassay: Importance of correct taxonomic identification of Artemia (Anostraca) species. *Environmental Toxicology* **23**:555-560.

728 Sánchez, M. I., Hortas, F., Figuerola, J., & Green, A. J. (2012). Comparing the potential for dispersal via waterbirds of a native and an invasive brine shrimp. *Freshwater Biology*, 57(9), 1896-1903. doi: 10.1111/j.1365-2427.2012.02852.x

731 Sebastián-González, E., & Green, A. J. (2014). Habitat use by waterbirds in relation to
732 pond size, water depth, and isolation: Lessons from a restoration in southern Spain.

733 *Restoration Ecology, 22*(3), 311-318. doi: 10.1111/rec.12078

734 Shannon, C. E. (1948). *A Mathematical Theory of Communication*. *Bell System Technical*

735 *Journal* 27(3), 379-423. doi: 10.1002/j.1538-7305.1948.tb01338.x

736 Van Stempvoort, D. R., Hendry, M. J., Schoenau, J. J., & Krouse, H. R. (1994). Sources and

737 dynamics of sulfur in weathered till, Western Glaciated Plains of North America.

738 *Chemical Geology*, **111**(1-4), 35-56. doi: 10.1016/0009-2541(94)90081-7

739 Wang, I. J., & Bradburd, G. S. (2014). Isolation by environment. *Molecular Ecology*,

740 23(23), 5649-5662. doi: 10.1111/mec.12938

741 Weir, B. S. (1996). Genetic Data Analysis II. Sinauer Associates, Sunderland, Mass.

742 Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis (Use R!). In *Springer*.

743 doi: 10.1007/978-0-387-98141-3

744 Wissel, B., Cooper, R. N., Leavitt, P. R., & Pham, S. V. (2011). Hierarchical regulation of

745 pelagic invertebrates in lakes of the northern Great Plains: A novel model for

746 interdecadal effects of future climate change on lakes. *Global Change Biology*, 17(1),

747 172-185. doi: 10.1111/j.1365-2486.2010.02291.x

748 Wood, W. W., & Sanford, W. E. (1990). Ground-water control of evaporite deposition.

749 *Economic Geology*, **85**(6), 1226-1235. doi: 10.2113/gsecongeo.85.6.1226

750 Wright, S. (1943). Isolation by distance. *Genetics*, *28*(28), 114-138.

751

752 **Figure captions**

753

754 Fig. 1 Principal Component Analysis of environmental variables in Saskatchewan saline

755 lakes. The percentage of total variation explained by the first two components is given

756 on the respective axes.

757

758 Fig. 2. a) Haplotype network (mtDNA) for Canadian *Artemia franciscana* (includes CO1

759 haplotypes from Saskatchewan and British Columbia). Grey squares highlight the four

760 main haplotype groups present in Saskatchewan. The top nodes of the network (grey

761 shades) represent members of two populations (BCO, PBC) in British Columbia ~1000

762 km west of the Saskatchewan study area. Haplotype groups A–C are based on genetic

763 differentiation between the populations that hold these haplotypes as identified by

764 pairwise ΦST and subsequent AMOVA (see Tables 2,3).

765 b) Geographic distribution of haplotype groups A–C. Numbered arrows indicate the

766 possible sequence of colonization after deglaciation (details in the text). Lines indicate

767 the glacial border at 14,000y BP (blue), 12,500y BP (green, dotted) and 11,500y BP (red,

768 hatched) (drawn after Christiansen, 1979). Scale is approximate, for exact distances between

769 lakes see Table S8. Note: to avoid confusion between patterns of mitochondrial and nuclear

770 markers, we applied different colour schemes in this figure and Fig. 4.

771

772 Fig. 3 a -d. EDENetwork for *Artemia franciscana* in Saskatchewan constructed from nSSR

773 loci. Each population is positioned according to its geographical coordinates. Edges

774 (connections between populations) are shown for three different thresholds at (a) 0.02,
<figure><img src="image_3.png" /><figcaption>b)</figcaption></figure>

<figure><img src="image_4.png" /></figure>
821

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

822
a)

<figure><img src="image_10.png" /></figure>

b)

<figure><img src="image_11.png" /></figure>

<footnote>823</footnote>

<footnote>824</footnote>

<footnote n="825">825 *Authors' contribution statement*</footnote>

<footnote>826</footnote>

<footnote n="827">827 AJG conceived the study and collected samples. MTB collected samples. CL and JSF</footnote>
828 generated the molecular data. MH generated and analysed the hydrochemical data. DF

829 analysed the molecular data with CL and wrote the paper with AJG and input from all

830 other authors.

831

832 *Conflict of interest*

833 The authors declare no conflict of interest.
Table 1. Genetic diversity of 14 *Artemia franciscana* populations from Saskatchewan. For nSSR markers (5 loci), rarefied total allelic richness (AR) is indicated as well as the number of private alleles (A<sub>p</sub>) and the Shannon-Wiener index for multilocus genotypes (H). SHO is excluded here, as only 8 individuals were genotyped. For mtDNA, haplotype diversity (*h*), nucleotide diversity (π) and the number of private haplotypes (h<sub>p</sub>) are indicated. Lake size is given as surface area in km<sup>2</sup>.

841

<table><thead><tr><th rowspan="2">Pop</th><th colspan="3">nSSR</th><th colspan="3">mtDNA</th><th rowspan="2">Area (km<sup>2</sup>)</th></tr><tr><th>AR</th><th>A<sub>p</sub></th><th>H</th><th>h</th><th>π</th><th>h<sub>p</sub></th></tr></thead><tbody><tr><td>AKE</td><td>40.9</td><td>7</td><td>2.876</td><td>0.779</td><td>0.00290</td><td>4</td><td>3.282</td></tr><tr><td>MAN</td><td>40.5</td><td>2</td><td>3.242</td><td>0.270</td><td>0.00047</td><td>4</td><td>15.599</td></tr><tr><td>PAT</td><td>26.3</td><td>4</td><td>3.445</td><td>0.000</td><td>0.00000</td><td>0</td><td>5.725</td></tr><tr><td>WHI</td><td>42.0</td><td>5</td><td>3.240</td><td>0.569</td><td>0.00129</td><td>5</td><td>29.618</td></tr><tr><td>BEN</td><td>50.6</td><td>4</td><td>3.149</td><td>0.380</td><td>0.00134</td><td>3</td><td>0.140</td></tr><tr><td>BES</td><td>53.8</td><td>1</td><td>3.442</td><td>0.411</td><td>0.00075</td><td>4</td><td>0.707</td></tr><tr><td>CHA</td><td>49.6</td><td>8</td><td>3.405</td><td>0.641</td><td>0.00246</td><td>8</td><td>23.206</td></tr><tr><td>EIN</td><td>50.3</td><td>4</td><td>3.493</td><td>0.260</td><td>0.00060</td><td>2</td><td>1.964</td></tr><tr><td>FRE</td><td>38.9</td><td>3</td><td>2.873</td><td>0.320</td><td>0.00122</td><td>2</td><td>0.890</td></tr><tr><td>FRS</td><td>54.9</td><td>9</td><td>3.441</td><td>NA</td><td>NA</td><td>NA</td><td>0.590</td></tr><tr><td>LYD</td><td>51.9</td><td>5</td><td>3.324</td><td>0.537</td><td>0.00223</td><td>5</td><td>2.061</td></tr><tr><td>GUL</td><td>50.1</td><td>4</td><td>2.493</td><td>NA</td><td>NA</td><td>NA</td><td>0.576</td></tr><tr><td>MUS</td><td>40.3</td><td>0</td><td>3.243</td><td>0.356</td><td>0.00123</td><td>7</td><td>18.284</td></tr><tr><td>MEA</td><td>NA</td><td>NA</td><td>NA</td><td>0.253</td><td>NA</td><td>4</td><td>0.900</td></tr></tbody></table>
<caption>Table 2: Pairwise genetic distances between *Artemia franciscana* populations in Saskatchewan. The lower diagonal shows Pairwise Φ<sub>ST</sub> (Tamura & Nei) based on the mtDNA marker (CO1). The upper diagonal shows D<sub>EST</sub> computed for five nSSR markers. Numbers in light grey have a non-significant associated p-value (> 0.01).</caption>

845

<table><thead><tr><th></th><th>AKE</th><th>BEN</th><th>BES</th><th>FRE</th><th>FRS</th><th>WHI</th><th>CHA</th><th>LYD</th><th>EIN</th><th>MEA</th><th>MAN</th><th>MUS</th><th>PAT</th><th>SHO</th><th>GUL</th></tr></thead><tbody><tr><td>AKE</td><td>0.341</td><td>0.321</td><td>0.326</td><td>0.374</td><td>0.316</td><td>0.316</td><td>0.332</td><td>0.296</td><td>-</td><td>0.323</td><td>0.418</td><td>0.331</td><td>0.322</td><td>0.313</td></tr><tr><td>BEN</td><td>0.476</td><td>-0.007</td><td>0.318</td><td>0.150</td><td>0.334</td><td>0.077</td><td>0.093</td><td>0.069</td><td>-</td><td>0.253</td><td>0.048</td><td>0.204</td><td>0.038</td><td>0.159</td></tr><tr><td>BES</td><td>0.574</td><td>0.016</td><td>0.294</td><td>0.155</td><td>0.296</td><td>0.075</td><td>0.086</td><td>0.058</td><td>-</td><td>0.235</td><td>0.070</td><td>0.181</td><td>0.044</td><td>0.139</td></tr><tr><td>FRE</td><td>0.493</td><td>0.000</td><td>0.038</td><td>0.188</td><td>0.140</td><td>0.159</td><td>0.233</td><td>0.213</td><td>-</td><td>0.108</td><td>0.391</td><td>0.225</td><td>0.260</td><td>0.166</td></tr><tr><td>FRS</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.305</td><td>0.015</td><td>0.072</td><td>0.086</td><td>-</td><td>0.255</td><td>0.139</td><td>0.326</td><td>0.082</td><td>0.105</td></tr><tr><td>WHI</td><td>0.521</td><td>0.013</td><td>-0.006</td><td>0.025</td><td>-</td><td>0.219</td><td>0.182</td><td>0.159</td><td>-</td><td>0.100</td><td>0.374</td><td>0.215</td><td>0.293</td><td>0.157</td></tr><tr><td>CHA</td><td>0.386</td><td>0.061</td><td>0.105</td><td>0.009</td><td>-</td><td>0.090</td><td>0.029</td><td>0.024</td><td>-</td><td>0.178</td><td>0.089</td><td>0.220</td><td>0.022</td><td>0.038</td></tr><tr><td>LYD</td><td>0.490</td><td>0.636</td><td>0.700</td><td>0.643</td><td>-</td><td>0.658</td><td>0.548</td><td>0.004</td><td>-</td><td>0.213</td><td>0.097</td><td>0.258</td><td>0.058</td><td>0.075</td></tr><tr><td>EIN</td><td>0.570</td><td>0.736</td><td>0.824</td><td>0.750</td><td>-</td><td>0.777</td><td>0.593</td><td>0.085</td><td>-</td><td>0.168</td><td>0.082</td><td>0.192</td><td>0.034</td><td>0.046</td></tr><tr><td>MEA</td><td>0.622</td><td>0.764</td><td>0.839</td><td>0.777</td><td>-</td><td>0.801</td><td>0.637</td><td>0.123</td><td>0.001</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>MAN</td><td>0.623</td><td>0.769</td><td>0.848</td><td>0.783</td><td>-</td><td>0.808</td><td>0.637</td><td>0.123</td><td>0.001</td><td>0.000</td><td>0.310</td><td>0.149</td><td>0.251</td><td>0.121</td></tr><tr><td>MUS</td><td>0.568</td><td>0.698</td><td>0.761</td><td>0.706</td><td>-</td><td>0.728</td><td>0.602</td><td>0.048</td><td>0.015</td><td>0.029</td><td>0.028</td><td>0.295</td><td>0.076</td><td>0.152</td></tr><tr><td>PAT</td><td>0.623</td><td>0.794</td><td>0.891</td><td>0.813</td><td>-</td><td>0.841</td><td>0.632</td><td>0.115</td><td>-0.007</td><td>-0.016</td><td>-0.015</td><td>0.016</td><td>0.229</td><td>0.218</td></tr><tr><td>SHO</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.081</td></tr></tbody></table>
846 Table 3 a, b. Population genetic structure of *Artemia franciscana* in Saskatchewan
847 identified by AMOVA for (a) CO1 haplotypes (mtDNA) with p-values obtained from
848 10,100 randomisations as implemented in Arlequin v3.5.2.2, and (b) five nSSR loci, with
849 associated p-values obtained by Monte Carlo permutation test with 999 repeats.
850 For a), population genetic structure was tested for the following groups of populations
851 with low genetic differentiation according to Table 2 (lower diagonal). Group1: AKE;
852 group2: BEN, BES, FRE, WHI, CHA; group3: EIN, LYD, MUS, MAN, PAT, MEA. Groups of
853 populations for b) were chosen according to connected and unconnected populations
854 resulting from network construction (Fig. 3a): AKE (group1), WHI (group2), MUS
855 (group3), MAN (group4), PAT (group5), FRE (group6), group7: BEN, BES, and group8:
856 EIN, LYD, CHA, GUL, FRS, SHO.

857 ,

<table><thead><tr><th>Source of variation</th><th>d.f.</th><th>Sum of squares</th><th>Variance components</th><th>% variation</th><th>p-value</th></tr></thead><tbody><tr><td colspan="6"><strong>a) CO1 haplotypes</strong></td></tr><tr><td>Between groups</td><td>2</td><td>135.42</td><td>0.839</td><td>68.06</td><td>&lt;10<sup>-5</sup></td></tr><tr><td>Between populations within groups</td><td>9</td><td>7.18</td><td>0.017</td><td>1.37</td><td>&lt;10<sup>-5</sup></td></tr><tr><td>Within populations</td><td>281</td><td>105.96</td><td>0.377</td><td>30.57</td><td>&lt;10<sup>-5</sup></td></tr><tr><td>Total</td><td>292</td><td>248.56</td><td>1.233</td><td></td><td></td></tr><tr><td colspan="6"><strong>b) nSSRs</strong></td></tr><tr><td>Between groups</td><td>7</td><td>176.90</td><td>0.321</td><td>14.09</td><td>0.001</td></tr><tr><td>Between populations within groups</td><td>6</td><td>24.70</td><td>0.055</td><td>2.41</td><td>0.001</td></tr><tr><td>Within populations</td><td>578</td><td>1100.03</td><td>1.903</td><td>83.50</td><td>0.001</td></tr><tr><td>Total</td><td>591</td><td>1301.62</td><td>2.279</td><td></td><td></td></tr></tbody></table>

858
Table 4. Summary of results of dbRDA for nSSR markers obtained by forward model selection (100 repeats) to test the effect of individual environmental variables (lake area, SO₄, Cl/SO₄, Alkalinity/SO₄, Alkalinity/Ca+Mg) on genetic distance (D<sub>EST</sub>, for details see Methods). Lake area, SO₄, Cl/SO₄, Alkalinity/SO₄ were not selected by forward model selection and therefore were not included in the table. "signif (n)" = the number of times that results were significant (of 100 repeats).

<table><thead><tr><th>statistics</th><th>Alk/Ca+Mg</th></tr></thead><tbody><tr><td>R2.adj</td><td>0.30</td></tr><tr><td>AIC</td><td>-16.3</td></tr><tr><td>F</td><td>4.49</td></tr><tr><td>p.mean</td><td>0.04</td></tr><tr><td>p.min</td><td>0.024</td></tr><tr><td>p.max</td><td>0.05</td></tr><tr><td>p.SD</td><td>0.007</td></tr><tr><td>signif (n)</td><td>87</td></tr></tbody></table>