Front. Genet. Frontiers in Genetics Front. Genet. 1664-8021 Frontiers Media S.A. 854362 10.3389/fgene.2022.854362 Genetics Original Research Contrasting Patterns of Genetic Diversity and Divergence Between Landlocked and Migratory Populations of Fish Galaxias maculatus, Evaluated Through Mitochondrial DNA Sequencing and Nuclear DNA Microsatellites Astorga et al. Genetic Divergence Among Populations of G. maculatus Astorga Marcela P. 1 Valenzuela Andrea 1 Segovia Nicolás I. 2 3 Poulin Elie 3 4 Vargas-Chacoff Luis 3 5 6 González-Wevar Claudio A. 3 5 6 * 1 Instituto de Acuicultura, Universidad Austral de Chile, Puerto Montt, Chile 2 Departamento de Biología Marina, Facultad de Ciencias del Mar, Universidad Católica Del Norte, Coquimbo, Chile 3 Instituto Milenio Biodiversidad de Ecosistemas Antárticos y Sub-Antárticos (BASE), Santiago, Chile 4 Departamento de Ciencias Ecológicas, Facultad de Ciencias, Universidad de Chile, Santiago, Chile 5 Instituto de Ciencias Marinas y Limnológicas (ICML), Facultad de Ciencias, Universidad Austral de Chile, Valdivia, Chile 6 Centro FONDAP de Investigaciones en Dinámica de Ecosistemas Marinos de Altas Latitudes (IDEAL), Universidad Austral de Chile, Valdivia, Chile

Edited by: Sandra Isabel Moreno Abril, University of Vigo, Spain

Reviewed by: Fernando Alda, University of Tennessee at Chattanooga, United States

Marta Barluenga, Spanish National Research Council (CSIC), Spain

 *Correspondence: Claudio A. González-Wevar, claudio.gonzalez@uach.cl

This article was submitted to Evolutionary and Population Genetics, a section of the journal Frontiers in Genetics

 19 05 2022 2022 13 854362 13 01 2022 08 04 2022 Copyright © 2022 Astorga, Valenzuela, Segovia, Poulin, Vargas-Chacoff and González-Wevar. 2022 Astorga, Valenzuela, Segovia, Poulin, Vargas-Chacoff and González-Wevar

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Galaxias species are interesting biogeographic models due to their distribution and different types of life cycles, with migratory and landlocked populations. To obtain a better understanding of the genetic consequences of the Quaternary glacial cycles in Galaxias maculatus, in this work we compared landlocked and migratory populations collected in areas that were differentially affected by ice advances and retreats. We included nine populations of G. maculatus, four collected from lakes (landlocked) and five from their associated estuaries/rivers (migratory) in three estuary-lake systems across southern Chile. Genetic analyses were performed using the mitochondrial control region and nine microsatellite loci. Genetic diversity measured with both markers was significantly higher in migratory than in landlocked populations across the study area. The levels of genetic differentiation showed higher differentiation among lakes than estuaries. Genetic diversity was higher in migratory populations located in areas that were less impacted by ice during Quaternary glacial processes. These results may be the consequence of recent recolonization of small freshwater bodies following the Last Glacial Maximum (LGM). Finally, the greatest differentiation was observed in populations that were exposed to continental ice advances and retreats during the LGM. Thus, in the present work we corroborate a pattern of differentiation between lakes and estuaries, using mtDNA sequences and microsatellite nuclear markers. This pattern may be due to a combination of biological factors, i.e., resident non-migratory behaviour or landlocking and natal homing-in, as well as geological factors, i.e., Expansion-Contraction Quaternary glacial biogeographic processes.

 mtDNA Chile estuary lake last glacial maximum 1210787 1110798 1151375 15150003 3120075 PFB 023 Fondo Nacional de Desarrollo Científico y Tecnológico10.13039/501100002850 Fondo de Financiamiento de Centros de Investigación en Áreas Prioritarias10.13039/501100018735 Instituto Millenium10.13039/501100007397 Comisión Nacional de Investigación Científica y Tecnológica10.13039/501100002848 1 Introduction

Galaxias maculatus is one of the freshwater fish species with the broadest distribution in cold and temperate regions of the Southern Hemisphere, including populations in southern Australia, Tasmania, New Zealand, South America, and the Falkland/Malvinas Islands (McDowall, 1970; Dyer, 2000; Habit et al., 2006). This species presents amphidromous behaviour, i.e., it reproduces in lakes and rivers and the larval stages migrate to the sea, where the individual develops before returning to complete its growth in freshwater ecosystems (McDowall, 1997). In addition, the presence of resident or landlocked populations has been observed in lakes and rivers, as well as migratory populations that maintain amphidromous behaviour (Delgado et al., 2019). This species is therefore considered an interesting model for historical biogeographic studies due to its particular distribution (Waters and Burridge, 1999; Waters et al., 2000; Zattara et al., 2005; Burridge et al., 2012; Waters et al., 2020), its response to glacial events (Zemlak et al., 2010; Zemlak et al., 2011; Carrea et al., 2012, 2013; González-Wevar et al., 2015a; González-Wevar et al., 2015b; Victoriano et al., 2020), and its migratory and non-migratory behaviour in the same habitat (Delgado et al., 2019; Delgado and Ruzzante, 2020). This variability makes G. maculatus a suitable model for evaluating genetic indicators associated with these unique characteristics, as this species has been considered to represent an intermediate evolutionary step between migratory behaviours of marine and freshwater organisms (Corush, 2019).

The glacial history of southern South America is relatively well understood. Continental ice sheet advances and retreats during the Quaternary occurred across the Pacific margin of Patagonia and generated major shifts in sea level, climate and landscape (Hulton et al., 2002; Rabassa et al., 2005, 2011; Hein et al., 2010). During the Last Glacial Maximum the Patagonian Ice Sheet expanded to an area stretching from 35°S to 56°S, and covered most of the Pacific fjords and channels of Patagonia (Clapperton, 1994; McCulloch et al., 2000). Accordingly, the Quaternary geomorphology of Patagonia varied significantly during this period and major glacial changes resulted in periodic extinctions of fauna associated with these near-shore ecosystems, allowing the colonization of vacant niches and creating opportunities for geographical isolation and speciation (Valdovinos et al., 2003; Zemlak et al., 2010; González-Wevar et al., 2011; Fraser et al., 2012; González-Wevar et al., 2012; González-Wevar et al., 2016; González-Wevar et al., 2018; Fernández-Iriarte et al., 2020). The Expansion-Contraction (E-C) model of Quaternary biogeography (Provan and Bennett, 2008) proposes that species and populations contracted their distributions to glacial refugia during glacial maxima. During interglacial periods they expanded their distributions towards previously glaciated areas following the deglaciation process (Provan and Bennett, 2008; Marko et al., 2010). Accordingly, unglaciated areas are expected to harbour higher levels of genetic diversity than ice-impacted areas, or newly founded ones. In contrast, glaciated areas should exhibit evidence of recent postglacial demographic expansions including lower levels of genetic polymorphism with small numbers of haplotypes dominating disproportionally large areas because of recolonization processes (Marko, 2004; Maggs et al., 2008). Moreover, recolonized areas should exhibit low divergence among haplotypes and lower levels of genetic differentiation than refuges (Provan and Bennett, 2008; Marko, 2004; Marko et al., 2010; González-Wevar et al., 2012; González-Wevar et al., 2013). This model provides a relatively simple paradigm to test demographic hypotheses during the Quaternary, and allows us to understand the response of populations and species to major climate change by recognizing distribution range shifts, potential refuge areas and recolonization routes (Hewitt, 2000, Hewitt, 2004; Maggs et al., 2008; Marko et al., 2010; Fraser et al., 2012; González-Wevar et al., 2013). G. maculatus includes migratory and landlocked populations distributed along the Pacific margin of South America between 30°S and the southern tip of Chilean Patagonia at 56°S. Phylogeographic and population-based analyses have recorded a small influence of Quaternary glacial cycles on the demography of the species (Zattara et al., 2005; Zemlak et al., 2010; Zemlak et al., 2011; González-Wevar. et al., 2015a; González-Wevar et al., 2015b). In general, G. maculatus is characterized by its high levels of genetic diversity and strong phylogeographic structure across its distribution (Zemlak et al., 2010; González-Wevar et al., 2015a; Delgado et al., 2019). Across non-glaciated areas between 38°S and 41°S each estuarine or riverine population represents a different genetic unit characterized by divergent haplogroups (genetic pool derived from shared mtDNA haplotypes). The diversity of G. maculatus includes four main haplogroups, three of them in the ice-free zone. The fourth haplogroup is the one that dominates Patagonia (south of 42°S). In contrast, ice-impacted populations along the Pacific margin, between 43°S and 53°S, were all included in a single haplogroup, with lower levels of polymorphism and genetic structure (González-Wevar et al., 2015a). This phylogeographic pattern is consistent with the Expansion-Contraction (E-C) Quaternary biogeographic model (Provan and Bennett, 2008).

Previous studies have shown diverse patterns of genetic differentiation in this species. Carrea et al. (2012) examined spatial patterns of genetic and phenotypic variability of G. maculatus from postglacial lakes in north-western Argentinian Patagonia and recognized three different genetic clusters using microsatellite markers. Moreover, genetic analyses of migratory and landlocked populations across a latitudinal gradient recorded higher levels of genetic structure in northern Argentinian Patagonia than in southern Patagonia (Carrea et al., 2013). Northern Patagonian populations of G. maculatus could have survived Quaternary glacial cycles in situ, while higher latitude populations were probably eradicated by extensive ice-sheet advances and retreats. Recently, Delgado et al. (2019), using neutral and putatively adaptive SNP loci, found that migratory populations in estuarine zones of the Chilean coast are highly differentiated from their landlocked counterparts. Moreover, migratory populations showed higher levels of gene flow and absence of site fidelity, while landlocked populations showed evidence of different colonization events with relatively low genetic diversity and varying levels of gene flow. Victoriano et al. (2020), using mtDNA sequences, observed that G. maculatus showed higher levels of genetic diversity and structure than the species G. platei, in the same river basin on both sides of the Andes. This result was explained by the different migratory and non-migratory behaviour of the two species.

To gain a better understanding of the genetic consequences of the Quaternary glacial in G. maculatus, and corroborate the genetic diversity and structure observed in the species, we compared landlocked and migratory populations from western Patagonia in areas that were differentially affected by ice advances and retreats. The populations were evaluated using two types of molecular markers, mitochondrial DNA sequences and nuclear DNA microsatellites. Following the E-C Quaternary biogeographic model, non-glaciated areas are expected to exhibit higher levels of genetic diversity and structure than heavily ice-impacted ones, for both migratory and landlocked populations. Similarly, migratory estuarine populations are expected to harbour higher levels of genetic polymorphism associated with older demographic histories, while landlocked populations should exhibit lower levels of genetic diversity and recent population expansions probably derived from associated estuarine/riverine ecosystems.

 2 Materials and Methods 2.1 Sampling, DNA Extraction and Amplification

Our analyses included nine populations of G. maculatus collected from lakes (n = 4) and associated estuaries (n = 5) in different biogeographical regions across the species distribution in southern Chile, in three areas that were differentially affected by Quaternary glacial processes (Figure 1). Firstly, we included two estuary-lake systems (E-LS) located in non-glaciated areas between 38°S and 42°S (Figure 1), also known as the Intermediate Area following Camus (2001). Localities in the first system (E-LS1) were Moncul Estuary (Mon-E); Lingue Estuary (Lin-E) and Colico Lake (Col-L) (Table 1). Localities in the second system (E-LS2) were: Maullín Estuary (Mau-E); Llanquihue Lake (Lla-L) and Pichilaguna Lake (Pic-L). We also included a third estuary-lake system located in an area that was heavily impacted by the Patagonian Ice Sheet in Magellan Province; localities in this third system (E-LS3) were: Tortel Estuary (Tor-R); Pascua Estuary (Pas-R); and Quetru Lake (Que-L) (Figure 1). Based on previous results (Gonzalez-Wevar et al., 2015a; Gonzalez-Wevar et al., 2015b), estuary populations are assumed to represent migratory populations and those located in lakes, landlocked ones. The geo-references of the sites are shown in Table 1 and their locations in Figure 1.

Map showing Galaxias maculatus sampling locations in western Patagonia separated into the three systems analysed.

Galaxias maculatus sampling sites in the estuaries and lakes of western Patagonia, separated into the three systems analysed, with sample size for molecular analysis.

System Locality Habitat Code Latitud Longitud N Mit N Nucl
E-LS1 Moncul Estuary MON-E 38°42′S 73°24′W 24 25
Lingue Estuary LIN-E 39°26′S 73°09′W 24 23
Colico Lake COL-L 39°03′S 72°03′W 27 24
E-LS2 Maullín Estuary MAU-E 41°35′S 73°38′W 27 31
Llanquihue Lake LLA-L 41°12′S 73°02′W 32 37
Pichilaguna Lake PIC-L 41°16′S 73°03′W 29 19
E-LS3 Tortel Estuary TOR-E 47°45′S 73°32′W 27 27
Pascua Estuary PAS-E 48°15′S 73°18′W 45 25
Quetru Lake QUE-L 48°11′S 73°13′W 31 30

About 300 individuals were sampled during the year 2013, varying in a range of 30–50 individuals per locality, using whitebait nets and storage in ethanol (95%) for DNA preservation.

DNA was extracted from muscle tissue using 1) the salting-out method described by Aljanabi and Martinez (1997), which was applied to the samples analysed in the Genomics and Molecular Ecology Lab in Valdivia; or 2) the standard DNA extraction protocol with ethanol precipitation (Taggart et al., 1992), applied to the samples analysed in the Molecular Genetic Lab in Puerto Montt. Both methods showed the same quality of DNA collection and were used in both markers. A partial fragment of the mitochondrial Control Region (D-loop) was amplified using the specific primers GAL-F and GAL-R designed from the complete mitochondrial genome of the species (González-Wevar et al., 2015b). Nine microsatellite loci were amplified using specific primers developed by Carrea et al. (2009). PCR conditions for mtDNA and microsatellite amplifications were defined following González-Wevar et al., 2015a, Gonzalez-Wevar et al., 2015b and Carrea et al. (2009) respectively. Mitochondrial DNA amplification products were purified using QIAquick Gel Extraction Kit (QIAGEN) and sequenced in both directions with an Automatic Sequencer ABI3730 × 1 at Macrogen Inc. (Seoul, South Korea). The sizes of the amplified microsatellite alleles were determined using an automatic DNA sequencer (ABI Prism 377; Applied Biosystems).

 2.2 Mitochondrial DNA Analysis

Mitochondrial sequences in G. maculatus were edited using GENEIOUS v.9.0.4 (http://www.geneious.com) and aligned with MUSCLE (Edgar, 2004). New D-loop sequences were deposited in GenBank under Accession Numbers: OM743508—OM743773. DNA saturation analysis was performed using DAMBE (Xia, 2013) to determine how saturation of transitions accumulates in relation to nucleotide divergence. Levels of genetic diversity of mtDNA marker were estimated through standard indices including: the number of haplotypes (k), the number of segregating sites (S), haplotype diversity (H), the average and number of pairwise differences (∏) for each locality, each estuary-lake system, and for the whole D-loop data set using DnaSP (Librado and Rozas 2009). We performed neutrality statistical tests (Tajima’s D and Fu’s FS) in DnaSP for each locality, each E-LS and the whole data set to estimate whether D-loop sequences in the species deviate from expectations under mutation-drift equilibrium. Considering the high levels of genetic diversity recorded in G. maculatus (Zemlak et al., 2010; González-Wevar CA. et al., 2015, Gonzalez-Wevar et al., 2015b, we determined levels of genetic differentiation between the localities analysed using mean pairwise differences (NST) following Pons and Petit (1996) in ARLEQUIN v. 3.5 (Excoffier and Lischer, 2010). The statistical significance of these analyses was estimated through 20,000 permutations. We estimated phylogeographic structure using the nearest neighbour statistic (Snn), which measures how often nearest neighbours in sequences (in sequence space) are from the same locality in geographic space (Hudson, 1990). The statistical significance of Snn was determined through 10,000 permutations. We determined the spatial genetic structure in G. maculatus through the number and composition of groups that were most differentiated based on mtDNA sequence data using SAMOVA (Spatial Analysis of Molecular Variance) (Dupanloup et al., 2002). SAMOVA is a popular analysis that uses multiple spatial scales in statistical methods to characterize spatial genetic structure based on pairwise genetic differences. We reconstructed mtDNA genealogical relationships in G. maculatus using Maximum Parsimony Networks in HAPVIEW (Salzburger et al., 2011). We also estimated the patterns of demographic history in the species by comparing the distribution of pairwise differences between haplotypes (mismatch distribution) for each locality, each E-LS and for the whole data set, up to the expected distribution under the sudden expansion growth model of Rogers and Harpending (1992) in DnaSP. Moreover, we estimated the mitochondrial age of the last demographic expansion in landlocked populations (Colico, Pichilaguna, and Quetru) under the sudden population growth model using the formula τ = 2 μt, where τ = the date of growth/decline measured in units of mutational time, t = time in years, and μ = the mutation rate per sequence per year. For this, we used a specific phylogeographic mutation rate estimated for the D-loop in galaxiid fish estimated by Burridge et al. (2008).

Finally, we performed phylogenetic reconstructions including the complete D-loop data set in G. maculatus, together with a sequence of G. platei as the outgroup. Bayesian reconstructions of haplotype relationships were done in MrBAYES v.3.1.2 (Huelsenbeck and Ronquist, 2001) using the GTR + I + G substitution model as previously selected with JModeltest v.2.1.10 (Darriba et al., 2012). Bayesian posterior probabilities (BPP) were estimated using the Metropolis coupled Markov-chain Monte-Carlo algorithm (MCMC). For this, we ran four chains for 150 × 106 generations and trees were sampled every 1,000 generations. Stationarity of this analysis was inferred when the average standard deviation of split frequencies was lower than 0.01 as suggested by Huelsenbeck and Ronquist (2001). The first 5% of the parameter values were discarded as burn-in and posterior probabilities were calculated as the fraction of trees showing a particular node. Finally, posterior probability-density was summarized as a maximum clade credibility tree using TreeAnnotator v.1.6.1 (http://beast.bio.ed.ac.uk/TreeAnnotator) and edited using the FigTree v.1.4.3 programme (http://tree.bio.ed.ac. uk/software/figtree).

 2.3 Microsatellite Analyses

Microsatellite fragment analysis was performed using Peak Scanner™ Software from Applied Biosystems and allele sizes were assigned to bins using FLEXIBIN (Amos et al., 2007). Finally, the evaluation of null alleles was reviewed and corrected using MICROCHECKER software v.2.2.3 (Van Oosterhout et al., 2004).

The statistical independence between loci was assessed using GENEPOP 4.7.5 (Rousset, 2008). Genotypic linkage disequilibrium between each pair of loci within populations was tested using Fisher’s exact test with a Markov chain. The deviation from genotypic proportion expected under the Hardy-Weinberg equilibrium was tested with exact p values by the Markov chain method using GENALEX 6.501 (Smouse and Peakall, 2012). Allele richness (Rs), standardized for the number of data, was calculated with FSTAT 2.9.4 (Goudet, 2002); the observed heterozygosity (Ho) and the expected heterozygosity (He) were estimated with GENALEX 6.501 (Smouse and Peakall, 2012). Microsatellite population genetic differentiation was assessed through FST comparisons between pairs of localities based on the total of polymorphic loci (Weir and Cockerham, 1984) using the sum of squared size differences (FST-like) with p-value significance, implemented in ARLEQUIN 3.5 (Excoffier and Lischer, 2010).

We used different approaches to estimate the spatial distribution of genetic diversity using multilocus data. We established patterns of genetic differentiation among localities in G. maculatus, through Bayesian analysis implemented in STRUCTURE v.2.3.2 (Pritchard et al., 2000). To determine which K best fitted the data, we used ΔK methods described by Evanno et al. (2005); this method specifies that the most likely K is that found when this probability is plotted against successive values of K asymptotes, and does not increase significantly with increasing K. This analysis was estimated using the Structure Harvester web page (Earl and VonHoldt, 2012).

Finally, we performed a discriminant analysis of principal components (DAPC) using ADEGENET 3.1-1 (Jombart, 2008; Jombart and Ahmed, 2011) in the R platform (R Core Team, 2014). For this analysis, the number of clusters (K) was identified using the find. clusters function, which is based on the lowest values of the Bayesian information criterion (BIC). This type of analysis defines a model in which genetic variation is divided within and between groups, and which at the same time maximises the synthesis of the variables within groups and minimises it between groups using the information of the precedence of each sample as a priori information. Discriminant analysis (DA) achieves the best separation of the individuals within predefined groups and makes a probabilistic assignment of individuals to each group.

 3 Results 3.1 Mitochondrial DNA

The complete mitochondrial DNA set in G. maculatus included 266 individuals collected from nine populations, with an average of 28 individuals per locality in a range from 19 to 45 individuals per locality (Table 2). The set consisted of 928 nucleotide positions. Considering that the analyses included a non-coding region of the mtDNA, several insertions and deletions were detected but not considered for future analyses. Control region sequences in G. maculatus were A—T rich (57.6%) compared to the G—C content (42.4%). Genetic diversity in G. maculatus was high with 197 variable positions (21.22%), of which 147 (74.61%) were parsimoniously informative. Levels of genetic diversity were much higher in migratory populations (estuaries) than in most of the landlocked populations (lakes), except for Llanquihue Lake, which exhibited high levels of polymorphism (Table 2). Haplotype diversity (H) ranged between 1.0 (Moncul, Maullín, Tortel and Pascua estuaries) and 0.729 (Quetru Lake) (Table 2). Similarly, the average numbers of nucleotide differences (Π) were higher in migratory than in landlocked populations, ranging from 30.33 (Maullín Estuary) to 1.72 (Colico Lake).

Diversity indices and neutrality tests in Galaxias maculatus populations in the different E-LS across the species distribution in western Patagonia.

D-Loop mitochondrial DNA Microsatellites nuclear DNA
Locality K H S Π D FS Rs Ho He F
MON-E 24 1.00 85 20.50 −0.76 −10.09* 10.62 ± 3.46 0.45 ± 0.07 0.78 ± 0.07 0.39 ± 0.10
LIN-E 18 0.94 77 22.77 −0.44 −0.61 10.56 ± 4.36 0.42 ± 0.04 0.76 ± 0.08 0.43 ± 0.06
COL-L 17 0.84 18 1.72 −2.29* −14.48* 5.86 ± 3.67 0.45 ± 0.11 0.57 ± 0.09 0.25 ± 0.14
E-LS1 59 0.97 103 20.75 −0.57 −21.70* 9.01 ± 4.34 0.44 ± 0.06 0.78 ± 0.06 0.42 ± 0.08
MAU-E 27 1.00 109 30.33 −0.22 −9.18* 9.68 ± 3.39 0.34 ± 0.06 0.71 ± 0.08 0.47 ± 0.11
LLA-L 28 0.99 86 11.65 −1.91* −13.31* 9.23 ± 3.79 0.59 ± 0.06 0.75 ± 0.04 0.20 ± 0.09
PIC-L 10 0.77 19 3.23 −1.15 −1.26 4.67 ± 2.84 0.57 ± 0.13 0.50 ± 0.09 −0.09 ± 0.18
E-LS2 65 0.98 162 36.63 −0.14 −12.53* 7.86 ± 3.97 0.48 ± 0.05 0.79 ± 0.05 0.34 ± 0.10
TOR-E 26 1.00 71 15.88 −1.01 −11.95* 10.37 ± 3.18 0.48 ± 0.06 0.81 ± 0.03 0.39 ± 0.07
PAS-E 43 1.00 75 13.24 −1.28 −32.40* 10.04 ± 3.48 0.42 ± 0.05 0.80 ± 0.05 0.48 ± 0.06
QUE-L 12 0.73 22 3.11 −1.60 −2.02 6.68 ± 2.40 0.49 ± 0.09 0.62 ± 0.08 0.27 ± 0.06
E-LS3 81 0.98 94 14.38 −1.18 −61.83* 9.03 ± 3.39 0.46 ± 0.05 0.79 ± 0.04 0.42 ± 0.06
Total 205 0.99 184 31.71 −0.45 −164.5* 11.70 ± 3.57 0.47 ± 0.03 0.70 ± 0.03 0.31 ± 0.04

Indices for D-loop from mitochondrial DNA where K, number of haplotypes; H, haplotype diversity; S, number of polymorphic sites; Π, average number of pairwise differences; D, Tajima’s D neutrality test and FS, Fu’s FS neutrality test and indicators with standard error of the nine microsatellite loci; Rs, allele richness; Ho, observed heterozygosity; He, expected heterozygosity; F, inbreeding coefficient.

Global genetic indices for each E-LS are marked in bold. Statistical significance for Neutrality Test (Fs) are indicated with an asterisk.

Pairwise NST comparisons detected a high degree of genetic differentiation between migratory and landlocked populations of G. maculatus along a latitudinal gradient. All pairwise population comparisons with the exception of one (between Pascua and Tortel Estuaries), were statistically significant, showing the high degree of genetic differentiation reported in the species (Table 3). Across the sampling area, levels of mtDNA genetic differentiation were higher among lakes ( X ¯ = 0.947) than among estuaries ( X ¯ = 0.513), followed by the differentiation between lake and estuary populations ( X ¯ = 0.669). Specifically, migratory populations of the southern system (Tortel and Pascua Estuaries E-LS3) showed no significant genetic differentiation (Table 3). The nearest neighbour statistic Snn in G. maculatus (Snn = 0.815) showed a high and significant phylogeographic signal in the whole D-loop data set (p < 0.001). SAMOVA supported the general pattern of mtDNA structure and identified eight groups with maximum differences accounting for 66.09% of the total genetic variance, while only 0.06% of the variance was due to within-group variation among localities (Table 4). SAMOVA recognized each locality, except for Pascua and Tortel Estuaries (E-LS3), as different populations.

Matrix of genetic differentiation values between pairs of Galaxias maculatus locations, estimated by mean general pairwise values of differentiation (NST) of Mitochondrial DNA data (above diagonal) and estimated by FST obtained from the nine microsatellite loci (below diagonal).

MON-E LIN-E COL-L MAU-E LLA-L PIC-L TOR-E PAS-E QUE-L
MON-E 0.315 0.650 0.509 0.249 0.788 0.594 0.634 0.738
LIN-E 0.014 0.431 0.391 0.412 0.730 0.417 0.467 0.618
COL-L 0.216 0.446 0.683 0.750 0.949 0.807 0.805 0.944
MAU-E 0.016 0.044 0.146 0.583 0.711 0.289 0.355 0.519
LLA-L 0.100 0.124 0.264 0.139 0.848 0.664 0.686 0.796
PIC-L 0.270 0.333 0.460 0.296 0.135 0.810 0.812 0.936
TOR-E 0.033 0.007 0.430 0.077 0.135 0.386 0.000 0.463
PAS-E 0.057 0.029 0.350 0.045 0.186 0.288 0.101 0.473
QUE-L 0.086 0.051 0.507 0.156 0.123 0.399 0.018 0.161

Statistical significant pairwise comparisons are marked in bold.

Spatial Analysis of Molecular Variance (SAMOVA) showing the percentage of genetic variation explained among groups (Moncul E, Lingue E, Colico L, Maullín E, Llanquihue L, Pichilaguna L, Quetru L, Pascua/Tortel E), among populations within groups and within populations using mtDNA. Where FSC represents differentiation within populations among groups and FCT represents differentiation among groups (***p < 0.001, **p < 0.01).

Source of variation d.f. Sum of squares Variance component Percentage of variation
Among groups 7 3067.993 13.34842 Va 66.09
Among populations within groups 1 7.237 0.01188 Vb 0.06
Within populations 257 1756.879 6.83611 Vc 33.85
Total 265 4832.109 20.19641

Fixation index FSC: 0.00173*** FCT: 0.66093**.

Statistical significant pairwise comparisons are marked in bold.

Considering the levels of genetic diversity recorded in the species, haplotype network analyses were divided into the different estuary-lake systems (E-LS1—E-LS3) analysed here (Figures 2A–C). Maximum parsimony haplotype network in E-LS1 (Figure 2A) included a total of 59 different haplotypes, and Moncul Estuary showed the highest diversity with a total of 24 different haplotypes (H = 1.00). Lingue Estuary exhibited a total of 18 different haplotypes (H = 0.94), three of which were shared by more than one individual. Considering the expanded genealogies recorded at Moncul (Π = 20.50) and Lingue (Π = 22.77) estuaries, both localities exhibited multimodal distributions of pairwise differences between haplotypes. In contrast, the Colico Lake parsimony network included a total of 17 different haplotypes. A dominant haplotype was recorded in 12 individuals (44.4%), with several closely associated singletons (Figure 2A). Considering the short genealogy with a typical star-like topology recorded at Colico Lake, the distribution of pairwise differences between haplotypes presented a unimodal distribution (not shown). Maximum parsimony network in E-LS2 (Figure 2B) included a total of 65 haplotypes; Maullín Estuary showed the highest diversity with 27 different haplotypes (H = 1.00). Llanquihue Lake showed high levels of genetic polymorphism with 28 different haplotypes (H = 0.99), three of which were found in more than one individual. Considering the expanded genealogy recorded in Maullín Estuary (Π = 30.33) and Llanquihue Lake (Π = 11.65), both localities exhibited multimodal distributions of pairwise differences between haplotypes. Interestingly, two individuals collected in Llanquihue Lake are closely related to the diversity recorded in Maullín Estuary. In contrast, the parsimony network for Pichilaguna Lake included a total of only 10 haplotypes. Two dominant haplotypes were recorded in this locality, the first was found in 11 specimens (37.93%) and the second in 9 (31.39%). Another eight low-frequency haplotypes were closely associated with the dominant ones (Figure 2B). The pairwise differences in Pichilaguna Lake showed a bimodal distribution. The maximum parsimony network in E-LS3 (Figure 2C) included a total of 81 different haplotypes. Tortel and Pascua Estuaries showed the highest diversity with 26 and 43 different haplotypes respectively. Considering the expanded genealogy recorded at Tortel (Π = 15.88) and Pascua (Π = 13.24) Estuaries, both localities exhibited multimodal distributions of pairwise differences between haplotypes. In contrast, the parsimony network in Quetru Lake included a total of 12 different haplotypes. A dominant haplotype was recorded in 16 individuals (51.61%) surrounded by several singletons and low-frequency haplotypes. Interestingly, two individuals from Quetru Lake were closely associated with the diversity recorded in the associated tributaries (Pascua and Tortel Estuaries) (Figure 2C). Considering this, the distribution of pairwise differences between haplotypes from Quetru Lake exhibited a bimodal distribution. The age of population expansions estimated for landlocked populations of G. maculatus were 13,470 years for Colico lake, 16,840 years for Pichilaguna lake and 12,120 years for Quetru lake.

Maximum parsimony haplotype networks in Galaxias maculatus based on mtDNA D-loop sequences from: (A) Estuary-Lake System 1: E-LS 1 (Moncul Estuary: Mon E, Lingue Estuary: Lin E and Colico Lake: Col L); (B) Estuary-Lake System 2: E-LS 2(Maullín Estuary: Mau E, Llanquihue Lake: Llan L and Pichilaguna Lake: Pic L); and (C) Estuary-Lake System 3: E-LS3 (Tortel Estuary: Tor E, Pascua Estuary: Pas E and Quetru Lake: Que L). Each haplotype is represented by a coloured circle showing the locality where it was collected. The size of the circles is proportional to the frequency in the sampling effort.

Phylogenetic reconstructions recognized three (HI, HIII and HIV) of the four main haplogroups previously recorded in G. maculatus in southern South America (Zemlak et al., 2010; González-Wevar et al., 2015a, Gonzalez-Wevar et al., 2015b. Haplogroup I was more abundant across the Intermediate Area localities while HIII was present in Lingue and Maullín estuaries (Supplementary Figure S1). Haplogroup IV was the dominant in estuarine and lacustrine populations located within the Patagonian Ice Sheet (Tortel Estuary, Pascua Estuary, and Quetru Lake) in areas that were heavily impacted by ice advances and retreats during glacial periods (Supplementary Figure S1). Some specimens from the localities of the Intermediate Area (Moncul Estuary, Lingue Estuary, Maullín Estuary and Llanquihue Lake) fell within the diversity of HIV, confirming the asymmetrical gene flow pattern (from north to south) previously found in the species (González-Wevar et al., 2015a).

 3.2 Microsatellites

Genetic data for the microsatellite loci analysed were obtained from a total of 241 individuals collected from nine populations, of which 131 were collected in estuaries (migratory populations) and 110 in lakes (landlocked populations). Deviation from the Hardy-Weinberg equilibrium was observed in different loci combined with different populations (Supplementary Table S1). Excess homozygosity was observed, mainly through the presence of null alleles in some loci. Allele richness was higher in migratory localities (Rs = 10.25) than in landlocked populations (Rs = 6.61), with significant differences (p < 0.01) between the groups and when estuaries and lakes were compared within each E-LS (p < 0.01). In general, estuarine migratory populations showed greater allele richness than landlocked ones, except for Llanquihue Lake (E-LS2), which showed similar diversity to those recorded in migratory estuarine localities. The highest levels of allele richness were observed in migratory populations from Moncul and Lingue Estuaries (Table 2), and the lowest in landlocked populations including Pichilaguna, Colico and Quetru Lakes. The inbreeding coefficient (F; Table 2) showed no significant differences across a latitudinal gradient. When the genetic diversity was compared among systems (E-LS), no major differences were observed using microsatellites.

Overall comparison between locations showed a value of FST = 0.149. A posteriori pairwise analysis between locations based on the FST estimator (Table 3) showed significant differences between localities; this is explained mainly by the landlocked populations which differ between all the locations, except for Quetru Lake which showed no significant differences from Tortel Estuary.

Comparing only migratory populations (estuaries) with each other, we observed a low although significant genetic differentiation with FST = 0.058 (p < 0.001). When landlocked populations (lakes) were compared with each other, a higher genetic differentiation was observed (FST = 0.228), with significant differences (p < 0.001). When the estuaries were compared with the lakes within each system, the lowest genetic differentiation was observed between the localities of E-LS 3.

Bayesian analysis of genetic structuring in STRUCTURE showed two genetic groups differentiated from one another, represented by the colours distributed among the nine locations sampled (Figure 3) in which estuaries are mainly separated from lakes. The correlation proposed by Evanno et al. (2005) shows the point where ΔK reaches its maximum value, estimating that K = 2 (Ln P(K) = −9377.4) is the most probable K for the genetic groups identified in the sample (Figure 3A). Genetic differentiation was structured with two genetic clusters: one group with the individuals belonging to Colico, Llanquihue and Pichilaguna lakes and other cluster formed by all the rivers and Quetru Lake. However, when the results were observed with a higher value of K (K = 5) (Figure 3B), it was possible to identify a pattern of differentiation, consistent with all previous statistical results, in which high similarity among the estuaries except for Maullin Estuary and high differentiation between the lakes is observed, except for Lake Llanquihue.

Genetic differentiation between locations, assessed by Bayesian analysis: (A) Diagram of the estimated population structure of Galaxias maculatus for the nine locations sampled with K = 2, the colours represent the genetic groups; (B) Diagram of population structure with K = 5. The name of the localities is observed on the lower horizontal axis: MON-E: Moncul Estuary; LIN-E: Lingue Estuary; COL-L: Colico Lake; MAU-E: Maullín Estuary; LLA-L: Llanquihue Lake; PIC-L Pichilaguna Lake; TOR-E: Tortel Estuary; PAS-E: Pascua Estuary and QUE-L Quetru Lake.

The results of the discriminant analysis of principal components (DAPC) (Figure 4) showed a consistent pattern of genetic differentiation, similar to that observed in the Fst analysis. The highest differentiation was observed in Colico Lake and the lowest among estuarine localities. The results again group the individuals from Moncul, Lingue and Pascua Estuaries together as single groups, but close to Tortel Estuary. However, a clear separation is observed between estuary-lake systems. The individuals from Colico Lake clearly separate into a differentiated group, Llanquihue and Pichilaguna Lakes are seen to be segregated jointly, with these landlocked populations highly differentiated. The exception was Maullín Estuary, which is seen to form a more differentiated group.

Ordering of the Galaxias maculatus clusters identified by the discriminant analysis of principal components (DAPC) of the different locations sampled. 1) Quetru Lake (QUE-L); 2: Pichilaguna Lake (PIC-L); 3: Lingue Estuary (LIN-E); 4: Colico Lake (COL-L); 5: Llanquihue Lake (LLA-L); 6: Maullín Estuary (MAU-E); 7: Moncul Estuary (MON-E)—Tortel Estuary (TOR-E)- Pascua Estuary (PAS-E).

 4 Discussion 4.1 Genetic Diversity 4.1.1 Mitochondrial DNA

As previously demonstrated in Galaxias maculatus (Zemlak et al., 2010; Zemlak et al., 2011; González-Wevar et al., 2015a; Gonzalez-Wevar et al., 2015b; Delgado et al., 2019), this species shows high levels of genetic diversity across its distribution in western South America. Among 266 individuals analysed with mtDNA, we recorded a total of 205 different haplotypes. On the one hand, levels of genetic diversity recorded in migratory populations of G. maculatus were higher across formerly unglaciated areas, while populations within the limits of the Patagonian Ice Sheet exhibited lower levels of polymorphism; accordingly, unglaciated areas represent more stable populations that were less affected by demographic processes associated with Quaternary glacial ice advances and retreats. On the other hand, the general diversity pattern recorded in migratory populations along a latitudinal gradient was not found in landlocked populations; levels of genetic diversity in landlocked populations of G. maculatus were lower and similar across our sampling effort, and probably a consequence of analogous demographic trajectories associated with recent bottlenecks and/or founder effects. Even though Colico and Pichilaguna lakes are located in unglaciated areas of western Patagonia, these lakes were covered by ice during the LGM (Moreno et al., 2018). Accordingly, the expansions of landlocked G. maculatus populations occurred after the Last Glacial Maximum, and hence are associated with post-glacial recolonization processes following the deglaciation process. As previously stated, the case of Llanquihue Lake is interesting and requires further study. Nevertheless, studies on G. maculatus have demonstrated a direct relationship between genetic diversity and the area of the lake analysed. Llanquihue is the second largest lake in South America; it is located at the limit of the expansion of glacial ice and probably served as a genetic reservoir for the species.

 4.1.2 Microsatellites

Microsatellite analysis in G. maculatus showed similar diversity patterns to those recorded in populations of the species from lakes and estuaries across Argentinian Patagonia (Carrea et al., 2012; Carrea et al., 2013). Our mean diversity and allele richness estimations are higher than those recorded by Carrea et al. (2013) in Argentinian Patagonia (Rs = 6.45). This study showed that migratory localities exhibited higher levels of genetic diversity than landlocked ones. It is to be expected that those populations with migratory behaviour will present greater genetic diversity than closed or landlocked populations, due to higher gene flow in the former and greater probability of genetic drift in the latter. Delgado et al. (2019) observe that estuarine migratory populations exhibited higher genetic diversity than resident populations, which is explained by the (probably low) diversity of the colonizing populations. Therefore, these differences in diversity between types of samples can be explained by bottlenecks generated by population reduction during glacial maxima, similar to those recorded by Zemlak et al. (2010) and Vera-Escalona et al., 2020.

 4.2 Genetic Structure 4.2.1 Mitochondrial DNA

As previously shown in this species (González-Wevar et al., 2015a), we also recorded contrasting patterns of genetic structure in migratory populations between the main biogeographic areas analysed. Each of the estuaries analysed north of the 42°S boundary, and across an area less than 250 km wide, represents a different genetic unit; the estuaries analysed in the southern system (E-LS3) of Patagonia exhibited lower levels of genetic differentiation. Again, major differences in terms of the genetic structure between unglaciated and glaciated areas are in basic agreement with the predictions of the Expansion-Contraction Quaternary biogeographic model in the south-eastern Pacific (Provan and Bennett, 2008; Pardo-Gandarillas et al., 2018; Fernández-Iriarte et al., 2020). González-Wevar CA. et al., 2015, González-Wevar et al., 2015b analysed populations from Chilean Patagonia and recorded low degrees of genetic structuring in estuary locations south of 42°S, belonging to the Magellan Biogeographical Province. These results coincide with observations in E-LS3 in this work, where the two estuaries showed no significant genetic differentiation. Thus, the genetic differences observed between the types of localities can be explained by the life cycles of the fish and their resident or migratory behaviour observed in the estuaries and rivers (Wallis et al., 2001; Delgado et al., 2019). Our results suggest a high probability of observing migratory behaviours in the estuary samples, and resident behaviours in lakes; the exception would be Maullín Estuary (MAU-E), where the slight genetic differentiation observed can be explained by the possible presence of a high proportion of individuals from resident populations.

Our results also show that the greatest differentiation is observed among all the landlocked populations. Quetru Lake, located in the area influenced by the LGM, showed the highest differentiation; it could have been affected by a process of extinction and recolonization from possible glacial refuges (Fraser et al., 2012).

 4.2.2 Microsatellites

Comparative analyses using microsatellites in G. maculatus detected significant differences among the populations analysed. However, pairwise analyses showed that these differences are mainly explained by the high genetic differentiation recorded in the landlocked populations, particularly those from Pichilaguna and Quetru Lakes. When landlocked and migratory populations were compared separately, we observed that the estuaries exhibit lower genetic structuring; in contrast, landlocked populations showed higher levels of genetic structure and fewer migrants. These observations were corroborated by the multivariate analyses (DAPC), where migratory populations from estuaries were grouped in a big central cloud; however, there are two clusters of higher similarity, one including Moncul, Lingue and Pascua and another one Tortel, all geographically distant from each another. These multivariate analyses showed that individuals sampled from landlocked populations could be separated from each other, and from migratory populations. Colico, Llanquihue and Pichilaguna Lakes showed the greatest similarity. Llanquihue and Pichilaguna Lakes are 3.5 km apart and are probably connected by streams. Llanquihue is the second-largest lake in South America, covering 860 km2; accordingly, it may represent a source for the G. maculatus population of Pichilaguna lake, which has an area of approximately 1 km2, although mtDNA data do not show shared genotypes. Zattara et al., 2005 found a positive correlation between genetic diversity, measured by isoenzymes, and lake area and perimeter; this agrees with the observations in our work, in which the lake with the largest area and perimeter (Llanquihue) has the highest allele richness, while the smallest lake (Pichilaguna) exhibited the lowest diversity of alleles. This is mainly due to the sizes reached by the populations in each lake; it is probable that populations of smaller effective size are more exposed to genetic drift or to the combined effects of bottlenecks, founder effects, endogamy, and restricted levels of gene flows, which would increase the inter-population genetic divergence in comparison with larger populations (Hartl and Clark, 1997). Just as was observed from the sequence data, Quetru Lake (E-LS3), located in the zone affected by the Last Glacial Maximum, shows greater differentiation and genetic structure than the lakes and estuaries located in the other two systems; this can be explained by the processes of extinction and recolonization in this area.

Finally, Bayesian Structure analysis showed a high differentiation among lakes and estuaries, separating G. maculatus populations into two main differentiated genetic groups. However, when evaluating this analysis with a greater number of differentiated groups (K = 5), it is possible to detect a differentiation pattern with high separation between lakes while the similarity between estuaries is maintained, which is consistent with our observations in the other analyses. The low differentiation detected among lakes by this Bayesian analysis may result from the high genetic diversity detected by microsatellite markers, due to their high mutational rate, which may generate allele similarity by homoplasies (Jarne and Lagoda, 1996; Estoup et al., 2002; Bhargava and Fuentes, 2010). The results obtained by the two types of markers are highly consistent; however, as observed in the structure result, a difference may be observed between the two molecular markers used. This could be due to the greater genetic diversity of microsatellites, which makes them very suitable for population analysis; however, in cases where the populations present high differentiation, it may be that the use of sequences as molecular markers can better detect the differentiation pattern than that shown by microsatellites. This would be due mainly to the possibility that the latter show high genetic similarity due to alleles that have already presented various mutations; in other words, for highly differentiated groups, the microsatellites could show homoplasies (Jarne and Lagoda, 1996; Estoup et al., 2002; Bhargava and Fuentes, 2010).

The different genetic statistical analyses show a consistent pattern of separation between the lakes analysed, and between lakes and estuaries. This pattern of genetic differentiation between lakes is similar to those recorded by Zattara et al., 2005, Delgado et al. (2019) and Rojo et al. (2020) using different types of molecular markers and analysing different watershed systems, and may be explained by the non-migratory behaviour of lacustrine populations. A similar pattern of genetic differentiation was recorded by Waters and Burridge (1999) when comparing mtDNA sequences from two Chilean lakes located in the same basin: Riñihue and Saval.

This high genetic differentiation, observed mainly in lakes and not in estuaries, may be because each lake is an “island”, and the separation between them acts as a barrier to gene flow. These results show that lakes behave as landlocked zones with little capacity for admixture, as there is probably little or no fish migration between lakes and estuaries. It has also been proposed that the level of structuring found in some basins (connected lakes, rivers, and estuaries) may be due to differential effects of Quaternary climate cycles in the different locations studied, since in some cases the population adopts a resident behaviour or life history, as observed in the variation to optional diadromous behaviour in a basin in Argentinian Patagonia (Carrea et al., 2013). This contrasts with the findings of Hollin and Schilling (1981) in locations at high latitudes which were strongly affected by glaciation in the mid-Pleistocene, the last glacial period when glaciers covered vast areas of Patagonia. This factor may have generated a refuge which now presents high genetic differentiation, as observed in the most southerly basin in Argentinian Patagonia analysed by Zemlak et al. (2008) and Carrea et al. (2013). In the present work, we observed a similar pattern of differentiation increasing with the latitude of the possible refuge zones affected by glaciation (Zemlak et al., 2008), mainly in Estuary-Lake System 3 (Tortel and Pascua Estuary with Quetru Lake), which would have been more strongly affected by the last glacial period than E-LS1 and E-LS2. However, a pattern of genetic differentiation was observed in the genetic diversity of E-LS3, explained by an effect of the last glacial period. In the case of lakes, which present high differentiation between one another and as compared to estuaries, the differentiation observed may be due to the possible presence of a refuge in each. In post-glaciation periods, these refuges would reconstruct highly differentiated populations. Thus, in the present work, using mtDNA sequences and microsatellite nuclear markers, we corroborated a pattern of differentiation between lakes and estuaries, similar to that already detected on a small scale by Zemlak et al. (2008) using mitochondrial DNA sequencing, and Delgado et al. (2019) using SNP markers. This corroborated the genetic separation of the lakes, which had already been evaluated through transplants by Delgado et al. (2020), where these differences come to be explained by evolutionary processes (natural selection or genetic drift); and by the presence of two types of behaviours in the same species, resident and migratory (Delgado et al., 2019). In addition, the effect of the Quaternary climate cycles on population structure and the loss of genetic diversity can be detected with both types of markers. Therefore, these patterns may be due to a combination of biological factors, i.e., resident non-migratory behaviour and/or landlocking and natal homing-in, as well as geological factors, i.e., E-C Quaternary biogeographic glaciation. It is therefore important to establish protection measures for this species, mainly in the lake populations due to the low diversity and high genetic differentiation observed in the lakes. Our findings need to be confirmed, to ensure the protection or conservation of these isolated populations.

 Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Center for Biotechnology Information (NCBI) BioProject database under accession numbers OM743508—OM743773.

 Ethics Statement

The animal study was reviewed and approved by the Comite de Etica-Bioetica Universidad Austral de Chile.

 Author Contributions

MPA wrote the manuscript, analyzed and interpreted data; AV carried out laboratory analysis and data analysis; NS carried out data analysis and prepared figures and tables and reviewed the manuscript; EP and LV-C contributed substantially to the drafting of the manuscript; CG-W conceived the study, wrote the manuscript, and analyzed and interpreted data.

 Funding

This study was supported by the postdoctorate 3120075 fellowship. At the same time, this research was partially funded by the Projects P05002 ICM and PFB 023 (IEB, Universidad de Chile), to EP and CG-W and projects Fondecyt 1110798 and 1151375 to LV-C. Fondap IDEAL program (15150003), Fondecyt Regular project 1210787 and ANID—Millenium Science Initiative—BASE. ANID — Programa Iniciativa Científica Milenio — ICN2021_002.

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.854362/full#supplementary-material

Bayesian maximum clade credibility tree of Galaxias maculatus haplotype relationships based on D-loop sequences. Bayesian posterior probabilities are shown above the nodes. Colored circles indicate sampling localities of the specimens analysed. Main haplogroups (HI, HIII, and HIV) recorded in G. maculatus are also indicated.

Deviation from the Hardy-Weinberg equilibrium in different microsatellites loci in different localities.

 References Aljanabi S. Martinez I. (1997). Universal and Rapid Salt-Extraction of High Quality Genomic DNA for PCR- Based Techniques. Nucleic Acids Res. 25, 4692–4693. 10.1093/nar/25.22.4692 Amos W. Hoffman J. I. Frodsham A. Zhang L. Best S. Hill A. V. (2007). Automated Binning of Microsatellite Alleles: Problems and Solutions. Mol. Ecol. Notes 7 (1), 10–14. 10.1111/j.1471-8286.2006.01560.x Bhargava A. Fuentes F. F. (2010). Mutational Dynamics of Microsatellites. Mol. Biotechnol. 44 (3), 250–266. 10.1007/s12033-009-9230-4 Burridge C. P. Craw D. Fletcher D. Waters J. M. (2008). Geological Dates and Molecular Rates: Fish DNA Sheds Light on Time Dependency. Mol. Biol. Evol. 25 (4), 624–633. 10.1093/molbev/msm271 Burridge C. P. McDowall R. M. Craw D. Wilson M. V. H. Waters J. M. (2012). Marine Dispersal as a Pre-requisite for Gondwanan Vicariance Among Elements of the Galaxiid Fish Fauna. J. Biogeogr. 39 (2), 306–321. 10.1111/j.1365-2699.2011.02600.x Camus P. A. (2001). Biogeografía marina de Chile continental. Rev. Chil. Hist. Nat. 74 (3), 587–617. 10.4067/s0716-078x2001000300008 Carrea C. Barriga J. P. Cussac V. E. Ruzzante D. E. (2012). Genetic and Phenotypic Differentiation Among Galaxias Maculatus Populations in a Patagonian Postglacial lake System. Biol. J. Linn. Soc. Lond. 107 (2), 368–382. 10.1111/j.1095-8312.2012.01939.x Carrea C. Cussac V. E. Ruzzante D. E. (2013). Genetic and Phenotypic Variation among Galaxias Maculatus populations Reflects Contrasting Landscape Effects between Northern and Southern Patagonia. Freshw. Biol. 58 (1), 36–49. 10.1111/fwb.12036 Carrea C. Paterson I. G. Cussac V. E. Ruzzante D. E. (2009). Ten Novel Microsatellite Loci Characterized for a Remarkably Widespread Fish: Galaxias Maculatus (Galaxiidae). Mol. Ecol. Resour. 9, 1503–1505. 10.1111/j.1755-0998.2009.02702.x Clapperton C. M. (1994). The Quaternary Glaciation of Chile: a Review. Rev. Chil. Hist. Nat. 67, 369–383. Corush J. B. (2019). Evolutionary Patterns of Diadromy in Fishes: More Than a Transitional State between marine and Freshwater. BMC Evol. Biol. 19 (1), 168–213. 10.1186/s12862-019-1492-2 Darriba D. Taboada G. L. Doallo R. Posada D. (2012). jModelTest 2: More Models, New Heuristics and Parallel Computing. Nat. Methods 9 (8), 772. 10.1038/nmeth.2109 Delgado M. L. Ruzzante D. E. (2020). Investigating Diadromy in Fishes and its Loss in an -Omics Era. iScience 23, 101837. 10.1016/j.isci.2020.101837 Delgado M. L. Górski K. Habit E. Ruzzante D. E. (2019). The Effects of Diadromy and its Loss on Genomic Divergence: The Case of Amphidromous Galaxias Maculatus Populations. Mol. Ecol. 28 (24), 5217–5231. 10.1111/mec.15290 Delgado M. L. Manosalva A. Urbina M. A. Habit E. Link O. Ruzzante D. E. (2020). Genomic Basis of the Loss of Diadromy in cursive: Insights from Reciprocal Transplant Experiments. Mol. Ecol. 29 (24), 4857–4870. 10.1111/mec.15686 Dupanloup I. Schneider S. Excoffier L. (2002). A Simulated Annealing Approach to Define the Genetic Structure of Populations. Mol. Ecol. 11, 2571–2581. 10.1046/j.1365-294x.2002.01650.x Dyer B. (2000). Systematic Review and Biogeography of the Freshwater Fishes of Chile. Estud. Oceanol. 19, 77–98. Earl D. A. VonHoldt B. M. (2012). STRUCTURE HARVESTER: a Website and Program for Visualizing STRUCTURE Output and Implementing the Evanno Method. Conservation Genet. Resour. 4 (2), 359–361. 10.1007/s12686-011-9548-7 Edgar R. C. (2004). MUSCLE: Multiple Sequence Alignment with High Accuracy and High Throughput. Nucleic Acids Res. 32 (5), 1792–1797. 10.1093/nar/gkh340 Estoup A. Jarne P. Cornuet J.-M. (2002). Homoplasy and Mutation Model at Microsatellite Loci and Their Consequences for Population Genetics Analysis. Mol. Ecol. 11 (9), 1591–1604. 10.1046/j.1365-294x.2002.01576.x Evanno G. Regnaut S. Goudet J. (2005). Detecting the Number of Clusters of Individuals Using the Software STRUCTURE: a Simulation Study. Mol. Ecol. 14 (8), 2611–2620. 10.1111/j.1365-294x.2005.02553.x Excoffier L. Lischer H. E. L. (2010). Arlequin Suite Ver 3.5: a New Series of Programs to Perform Population Genetics Analyses under Linux and Windows. Mol. Ecol. Resour. 10 (3), 564–567. 10.1111/j.1755-0998.2010.02847.x Fernández Iriarte P. J. González-Wevar C. A. Segovia N. I. Rosenfeld S. Hüne M. Fainburg L. (2020). Quaternary Ice Sheets and Sea Level Regression Drove Divergence in a marine Gastropod along Eastern and Western Coasts of South America. Sci. Rep. 10 (1), 844. 10.1038/s41598-020-57543-4 Fraser C. I. Nikula R. Ruzzante D. E. Waters J. M. (2012). Poleward Bound: Biological Impacts of Southern Hemisphere Glaciation. Trends Ecol. Evol. 27 (8), 462–471. 10.1016/j.tree.2012.04.011 González-Wevar C. A. Rosenfeld S. Segovia N. I. Hüne M. Gérard K. Ojeda J. (2016). Genetics, Gene Flow, and Glaciation: The Case of the South American Limpet Nacella Mytilina . PLoS ONE 11, e0161963. 10.1371/journal.pone.0161963 González-Wevar C. A. Salinas P. Hüne M. Segovia N. I. Vargas-Chacoff L. Astorga M. (2015a). Phylogeography in Galaxias Maculatus (Jenyns, 1848) along Two Biogeographical Provinces in the Chilean Coast. PloS ONE 10 (7), e0131289. 10.1371/journal.pone.0131289 González-Wevar C. A. Hüne M. Cañete J. I. Mansilla A. Nakano T. Poulin E. (2012). Towards a Model of Postglacial Biogeography in Shallow marine Species along the Patagonian Province: Lessons from the Limpet Nacella Magellanica (Gmelin, 1791). BMC Evol. Biol. 12, 139. 10.1186/1471-2148-12-139 González-wevar C. A. Nakano T. Cañete J. I. Poulin E. (2011). Concerted Genetic, Morphological and Ecological Diversification in Nacella Limpets in the Magellanic Province. Mol. Ecol. 20, 1936–1951. 10.1111/j.1365-294x.2011.05065.x González-Wevar C. Segovia T. Morley S. Chown S. Poulin E. (2013). Extinction and Recolonization of Maritime Antarctica in the Limpet Nacella Concinna (Strebel, 1908) During the Last Glacial Cycle: Toward a Model of Quaternary Biogeography in Shallow Antarctic Invertebrates. Mol. Ecol. 22, 5221–5236. 10.1111/mec.12465 González-Wevar C. A. Segovia N. I. Rosenfeld S. Ojeda J. Hüne M. Naretto J. (2018). Unexpected Absence of Island Endemics: Long-Distance Dispersal in Higher Latitude Sub-Antarctic Siphonaria (Gastropoda: Euthyneura) Species. J. Biogeogr. 45, 874–884. 10.1111/jbi.13174 González-Wevar C. Salinas P. Hüne M. Segovia N. Vargas-Chacoff L. Oda E. (2015b). Contrasting Genetic Structure and Diversity of Galaxias maculatus (Jenyns, 1848) along the Chilean Coast: Stock Identification for Fishery Management. Jhered 106, 439–447. 10.1093/jhered/esv005 Goudet J. (2002). FSTAT (Version 2.9. 3.2): A Program to Estimate and Test Gene Diversities and Fixation Indices. Available at: http://www. unil. ch/izea/software/fstat. html. Habit E. Dyer B. Vila I. (2006). Estado de conocimiento de los peces dulceacuícolas de Chile. Gayana 70 (1), 100–113. 10.4067/s0717-65382006000100016 Hartl D. L. Clark A. G. (1997). Principles of Population Genetics. Second edition. Sunderland: Sinauer Associates, 682. Hein A. S. Hulton N. R. J. Dunai T. J. Sugden D. E. Kaplan M. R. Xu S. (2010). The Chronology of the Last Glacial Maximum and Deglacial Events in central Argentine Patagonia. Quat. Sci. Rev. 29, 1212–1227. 10.1016/j.quascirev.2010.01.020 Hewitt G. M. (2004). Genetic Consequences of Climatic Oscillations in the Quaternary. Phil. Trans. R. Soc. Lond. B 359, 183–195. 10.1098/rstb.2003.1388 Hewitt G. (2000). The Genetic Legacy of the Quaternary Ice Ages. Nature 405, 907–913. 10.1038/35016000 Hollin J. T. Schilling D. H. (1981). “Late Wisconsin-Wichselian Mountain Glaciers and Small Ice Caps,” in The Last Great Ice Sheets: 179-206. Editors Denton J H Hughes T J (New York: Wiley). Available at: https://oceanrep.geomar.de/id/eprint/37507. Hudson R. (1990). Gene Genealogies and the Coalescent Process. Oxf. Surv. Evol. Biol. 7, 1–44. Huelsenbeck J. P. Ronquist F. (2001). MRBAYES: Bayesian Inference of Phylogenetic Trees. Bioinformatics 17 (8), 754–755. 10.1093/bioinformatics/17.8.754 Hulton N. Purves P. McCulloch R. Sugden D. Bentley M. (2002). The Last Glacial Maximum and Deglaciation in Southern South America. Quat. Sci. Rev. 21, 233–241. 10.1016/S0277-3791(01)00103-2 Jarne P. Lagoda P. J. L. (1996). Microsatellites, from Molecules to Populations and Back. Trends Ecology Evolution 11 (10), 424–429. 10.1016/0169-5347(96)10049-5 Jombart T. (2008). Adegenet: a R Package for the Multivariate Analysis of Genetic Markers. Bioinformatics 24 (11), 1403–1405. 10.1093/bioinformatics/btn129 Jombart T. Ahmed I. (2011). Adegenet 1.3-1: New Tools for the Analysis of Genome-wide SNP Data. Bioinformatics 27 (21), 3070–3071. 10.1093/bioinformatics/btr521 Librado P. Rozas J. (2009). DnaSP V5: a Software for Comprehensive Analysis of DNA Polymorphism Data. Bioinformatics 25 (11), 1451–1452. 10.1093/bioinformatics/btp187 Maggs C. A. Castilho R. Foltz D. Henzler C. Jolly M. T. Olsen J. (2008). Evaluating Signatures of Glacial Refugia for North Atlantic Benthic marine Taxa. Ecology 89, S108–S122. 10.1890/08-0257.1 Marko P. B. (2004). 'What's Larvae Got to Do with it?' Disparate Patterns of post-glacial Population Structure in Two Benthic marine Gastropods with Identical Dispersal Potential. Mol. Ecol. 13, 597–611. 10.1046/j.1365-294x.2004.02096.x Marko P. B. Hoffman J. M. Emme S. A. McGovern T. M. Keever C. C. Nicole Cox L. (2010). The 'Expansion-Contraction' Model of Pleistocene Biogeography: Rocky Shores Suffer a Sea Change? Mol. Ecol. 19, 146–169. 10.1111/j.1365-294x.2009.04417.x McCulloch R. D. Bentley M. J. Purves R. S. Hulton N. R. J. Clapperton D. E. (2000). Climatic Inferences from Glacial and Palaeoecological Evidence at the Last Glacial Termination, Southern South America. J. Quat. Sci. 15 (4), 409–417. 10.1002/1099-1417(200005)15:4<409::aid-jqs539>3.0.co;2 McDowall R. M. (1970). Comments on a New Taxonomy ofRetropinna(Galaxioidei: Retropinnidae). New Zealand J. Mar. Freshw. Res. 4, 312–324. 10.1080/00288330.1970.9515349 McDowall R. M. (1997). The Evolution of Diadromy in Fishes (Revisited) and its Place in Phylogenetic Analysis. Rev. Fish. Biol. Fish. 7 (4), 443–462. 10.1023/a:1018404331601 Moreno P. I. Videla J. Valero-Garcés B. Alloway B. V. Heusser L. E. (2018). A Continuous Record of Vegetation, Fire-Regime and Climatic Changes in Northwestern Patagonia Spanning the Last 25,000 Years. Quat. Sci. Rev. 198, 15–36. 10.1016/j.quascirev.2018.08.013 Pardo‐Gandarillas M. C. Ibáñez C. M. Torres F. I. Sanhueza V. Fabres A. Escobar‐Dodero J. (2018). Phylogeography and Species Distribution Modelling Reveal the Effects of the Pleistocene Ice Ages on an Intertidal Limpet from the South‐eastern Pacific. J. Biogeogr. 45 (8), 1751–1767. 10.1111/jbi.13362 Peakall R. Smouse P. E. (2012). GenAlEx 6.5: Genetic Analysis in Excel. Population Genetic Software for Teaching and Research-Aan Update. Bioinformatics 28 (19), 2537–2539. 10.1093/bioinformatics/bts460 Pons O. Petit R. J. (1996). Measwring and Testing Genetic Differentiation with Ordered versus Unordered Alleles. Genetics 144 (3), 1237–1245. 10.1093/genetics/144.3.1237 Pritchard J. K. Stephens M. Donnelly P. (2000). Inference of Population Structure Using Multilocus Genotype Data. Genetics 155 (2), 945–959. 10.1093/genetics/155.2.945 Provan J. Bennett K. (2008). Phylogeographic Insights into Cryptic Glacial Refugia. Trends Ecol. Evol. 23 (10), 564–571. 10.1016/j.tree.2008.06.010 R Core Team (2014). R: A Language and Environment for Statistical Computing. Second edition. Vienna, Austria: R Foundation for Statistical Computing. Rabassa J. Coronato A. Martínez O. (2011). Late Cenozoic glaciations in Patagonia and Tierra del Fuego: an updated review. Biol. J. Linn. Soc. 103, 316–335. 10.1111/j.1095-8312.2011.01681.x Rabassa J. Coronato A. M. Salemme M. (2005). Chronology of the Late Cenozoic Patagonian Glaciations and Their Correlation with Biostratigraphic Units of the Pampean Region (Argentina). J. South Am. Earth Sci. 20, 81–103. 10.1016/j.jsames.2005.07.004 Rogers A. R. Harpending H. (1992). Population Growth Makes Waves in the Distribution of Pairwise Genetic Differences. Mol. Biol. Evol. 9, 552–569. 10.1093/oxfordjournals.molbev.a040727 Rojo J. H. Fernández D. A. Figueroa D. E. Boy C. C. (2020). Phenotypic and Genetic Differentiation between Diadromous and Landlocked puyenGalaxias Maculatus. J. Fish. Biol. 96 (4), 956–967. 10.1111/jfb.14285 Rousset F. (2008). genepop'007: a Complete Re-implementation of the Genepop Software for Windows and Linux. Mol. Ecol. Resour. 8 (1), 103–106. 10.1111/j.1471-8286.2007.01931.x Salzburger W. Ewing G. B. Von Haeseler A. (2011). The Performance of Phylogenetic Algorithms in Estimating Haplotype Genealogies with Migration. Mol. Ecol. 20, 1952–1963. 10.1111/j.1365-294x.2011.05066.x Taggart J. B. Hynes R. A. Prodöuhl P. A. Ferguson A. (1992). A Simplified Protocol for Routine Total DNA Isolation from Salmonid Fishes. J. Fish. Biol. 40 (6), 963–965. 10.1111/j.1095-8649.1992.tb02641.x Valdovinos C. Navarrete S. A. Marquet P. A. (2003). Mollusk Species Diversity in the Southeastern Pacific: Why Are There More Species towards the Pole? Ecography 26, 139–144. 10.1034/j.1600-0587.2003.03349.x Van Oosterhout C. Hutchinson W. F. Wills D. P. M. Shipley P. (2004). Micro-Checker: Software for Identifying and Correcting Genotyping Errors in Microsatellite Data. Mol. Ecol. Notes 4 (3), 535–538. 10.1111/j.1471-8286.2004.00684.x Vera-Escalona I. Delgado M. L. Habit E. Ruzzante D. E. (2020). Historical and Contemporary Diversity of Galaxiids in South America: Biogeographic and Phylogenetic Perspectives. Diversity 12 (8), 304. 10.3390/d12080304 Victoriano P. F. Muñoz-Ramírez C. P. Canales-Aguirre C. B. Jara A. Vera-Escalona I. Burgos-Careaga T. Muñoz-Mendoza C. Habit E. M. (2020). Contrasting Evolutionary Responses in Two Co-distributed Species of Galaxias (Pisces, Galaxiidae) in a River from the Glaciated Range in Southern Chile. R. Soc. Open Sci. 7 (7), 200632. 10.1098/rsos.200632 Wallis G. P. Judge K. F. Bland J. Waters J. M. Berra T. M. (2001). Genetic Diversity in New Zealand Galaxias Vulgaris Sensu Lato (Teleostei: Osmeriformes: Galaxiidae): a Test of a Biogeographic Hypothesis. J. Biogeogr. 28 (1), 59–67. 10.1046/j.1365-2699.2001.00535.x Waters J. M. Burridge C. P. Craw D. (2020). River Capture and Freshwater Biological Evolution: A Review of Galaxiid Fish Vicariance. Diversity 12 (6), 216. 10.3390/d12060216 Waters J. M. Burridge C. P. (1999). Extreme Intraspecific Mitochondrial DNA Sequence Divergence inGalaxias maculatus(Osteichthys: Galaxiidae), One of the World's Most Widespread Freshwater Fish. Mol. Phylogenet. Evol. 11, 1–12. 10.1006/mpev.1998.0554 Waters J. M. Dijkstra L. H. Wallis G. P. (2000). Biogeography of a Southern Hemisphere Freshwater Fish: How Important Is marine Dispersal? Mol. Ecol. 9, 1815–1821. 10.1046/j.1365-294x.2000.01082.x Weir B. S. Cockerham C. C. (1984). Estimating F -Statistics for the Analysis of Population Structure. Evolution 38 (6), 1358–1370. 10.1111/j.1558-5646.1984.tb05657.x Xia X. (2013). DAMBE5: a Comprehensive Software Package for Data Analysis in Molecular Biology and Evolution. Mol. Biol. Evol. 30 (7), 1720–1728. 10.1093/molbev/mst064 Zattara E. Premoli A. (2005). Genetic Structuring in Andean Landlocked Populations of Galaxias Maculatus: Effects of Biogeographic History. J. Biogeo. 32, 5–14. 10.1111/j.1365-2699.2004.01164.x Zemlak T. S. Habit E. M. Walde S. J. Battini M. A. Adams E. D. M. Ruzzante D. E. (2008). Across the Southern Andes on Fin: Glacial Refugia, Drainage Reversals and a Secondary Contact Zone Revealed by the Phylogeographical Signal ofGalaxias Plateiin Patagonia. Mol. Ecol. 17 (23), 5049–5061. 10.1111/j.1365-294x.2008.03987.x Zemlak T. S. Habit E. M. Walde S. J. Carrea C. Ruzzante D. E. (2010). Surviving Historical Patagonian Landscapes and Climate: Molecular Insights from Galaxias Maculatus . BMC Evol. Biol. 10 (1), 67. 10.1186/1471-2148-10-67 Zemlak T. S. Walde S. J. Habit E. M. Ruzzante D. E. (2011). Climate-induced Changes to the Ancestral Population Size of Two Patagonian Galaxiids: the Influence of Glacial Cycling. Mol. Ecol. 20 (24), 5280–5294. 10.1111/j.1365-294x.2011.05352.x