- Research article
- Open Access
Genetic and morphological divergence at a biogeographic break in the beach-dwelling brooder Excirolana hirsuticauda Menzies (Crustacea, Peracarida)
BMC Evolutionary Biologyvolume 19, Article number: 118 (2019)
There is a biogeographic break located at 30°S in the southeast Pacific, in a coastal area of strong environmental discontinuities. Several marine benthic taxa with restricted dispersal have a coincident phylogeographic break at 30°S, indicating that genetic structure is moulded by life history traits that limit gene flow and thereby promote divergence and speciation. In order to evaluate intraspecific divergence at this biogeographic break, we investigated the genetic and morphological variation of the directly developing beach isopod Excirolana hirsuticauda along 1900 km of the southeast Pacific coast, across 30°S.
The COI sequences and microsatellite data both identified a strong discontinuity between populations of E. hirsuticauda to the north and south of 30°S, and a second weaker phylogeographic break at approximately 35°S. The three genetic groups were evidenced by different past demographic and genetic diversity signatures, and were also clearly distinguished with microsatellite data clustering. The COI sequences established that the genetic divergence of E. hirsuticauda at 30°S started earlier than divergence at 35°. Additionally, the three groups have different past demographic signatures, with probable demographic expansion occurring earlier in the southern group (south of 35°S), associated with Pleistocene interglacial periods. Interestingly, body length, multivariate morphometric analyses, and the morphology of a fertilization-related morphological character in males, the appendix masculina, reinforced the three genetic groups detected with genetic data.
The degree of divergence of COI sequences, microsatellite data, and morphology was concordant and showed two geographic areas in which divergence was promoted at differing historical periods. Variation in the appendix masculina of males has probably promoted reproductive isolation. This variation together with gene flow restrictions promoted by life history traits, small body size, oceanographic discontinuities and sandy-beach habitat continuity, likely influenced species divergence at 30°S in the southeast Pacific coast. The degree of genetic and morphological differentiation of populations to the north and south of 30°S suggests that E. hirsuticauda harbours intraspecific divergence consistent with reproductive isolation and an advanced stage of speciation. The speciation process within E. hirsuticauda has been shaped by both restrictions to gene flow and a prezygotic reproductive barrier.
Biogeographic breaks are areas of shifts in species composition and represent the limits of biogeographic regions. These discontinuities are often located in areas with prominent topographical features or steep environmental discontinuities. Along the southeast Pacific coast, biogeography of marine communities is characterized by two major biogeographic breaks. The southern break is located at 42°S and corresponds to the northern limit of the Magellan Province that is dominated by sub-Antarctic species and is characterized by a topography of fjords, islands and channels, with strong influence of rivers and glaciers . Also occurring at 42°S is the bifurcation of the West Wind Drift Current into the Humboldt Current System (HCS) to the north, and the Cape Horn Current towards the south. Additionally, the northern limit of the ice sheet during the last glacial period was located at 42°S.
To the north of 42°S the coast is mainly linear and under the influence of the HCS [1, 2]. The second biogeographic break is along the linear portion of the coast, at 30°S. At this latitude, there are no physical barriers or divergence of currents that may explain the presence of a biogeographic break. The 30°S biogeographic break divides to the north, the Peruvian Province, and to the south, the Intermediate Area (between 30°-42°S), and represents the limits of the geographic range of distribution of temperate and warm-affinity taxa, respectively [1, 2].
The HCS is an eastern boundary system with coastal wind-driven upwelling . Wind, upwelling, and offshore transport of upwelled waters are enhanced at 30°S, which corresponds to the average location of the subtropical anticyclone [2, 4, 5]. These features create extensive changes in kinetic energy of surface waters, temperature, salinity, oxygen concentration, and in thermal seasonal variations [6, 7]. Associated with the discontinuities in environmental characteristics at 30°S, there are important changes in the structure of marine communities and in larval recruitment patterns, including a significantly lower larval recruitment and adult abundance [8,9,10,11,12].
A common pattern along biogeographic breaks is a concordant genetic discontinuity, or phylogeographic break, of intraspecific lineages [13, 14]. Phylogeographic breaks result from ancient and sometimes ongoing barriers to gene flow that promote divergence . Concordance in the location of intraspecific genetic discontinuities suggest that there are common and persistent factors that promote intraspecific genetic structure , particularly for species with low dispersal capability [17,18,19,20,21,22,23,24]. Even though many factors influence intraspecific patterns of genetic diversity, including historical factors, developmental mode, planktonic larval duration, habitat availability and environmental characteristics, the intrinsic dispersal potential of a species often has the strongest influence on intraspecific genetic structure. The fact that dispersal potential is so relevant emphasizes the role of barriers to gene flow at biogeographic breaks as promoters of intraspecific divergence and speciation [20, 24,25,26,27].
In concordance with patterns detected for other marine biogeographic breaks, the 30°S break in the southeast Pacific also harbours phylogeographic breaks in species with low dispersal potential, including algae , intertidal barnacles [29, 30], gastropods [24, 31], and the amphipod Orchestoidea tuberculata . All of these species have limited dispersal potential, either because they have short-lived spores such as algae, or relatively short-lived dispersive larval stages, or completely lack dispersive stages, as is the case of O. tuberculata. Concordant intraspecific phylogeographic breaks in species with low dispersal potential suggest that the 30°S biogeographic break is a barrier to gene flow that originated in the past .
Strong environmental discontinuities, as those encountered at some biogeographic breaks, promote speciation because they maintain isolation of populations to each side of the biogeographic break . In the context of marine speciation, phylogeographic breaks can be considered as a stage of the speciation process, which can be understood as a continuum between a panmictic population and the complete reproductive isolation between two or more populations [33, 34]. In a scenario where gene flow restrictions are persistent in time, genetic drift will lead to increasing intraspecific divergence over time on either side of the barrier [24, 35]. Along the speciation divergence gradient, advanced stages are characterized by strong reproductive isolation, genetic discontinuities and lineage sorting with reciprocal monophyly . Intraspecific divergence and speciation are more likely to occur with more extreme environmental discontinuities, with greater temporal persistence of physical discontinuities, and/or with lower dispersal potential . Estimates of intraspecific divergence of marine benthic taxa at the 30°S biogeographic break indicate that it may have started as early as the Pleistocene, over 200,000 years ago [24, 37]. Taxa that have larvae that spend less than two weeks in the water column also display divergence and the onset has been dated during or soon after the last glacial period . The 30°S phylogeographic break is deeper in taxa with very low dispersal potential suggesting that for those taxa the break’s onset was thousands of years ago. The ancient onset of the break and the current environmental discontinuities that occur at 30°S, make this phylo- and bio-geographic break a natural scenario for the evaluation of intraspecific divergence processes.
Excirolana hirsuticauda Menzies (Crustacea, Peracarida) is a beach-dwelling isopod of > 12 mm of body length with low dispersal potential and a wide geographic distribution (ca. 22°S-43°S), crossing the 30°S biogeographic break. Shore-dwelling peracarids generally display strong spatial genetic structure [38,39,40,41,42] and phylogeographic discontinuities coincident with biogeographic breaks . The sympatric congener of E. hirsuticauda, E. braziliensis, harbours three reciprocally monophyletic lineages, two of which are parapatric at ~ 30°S [43, 44].
In this study, we focus on the 30°S latitude in the southeast Pacific as a natural laboratory to explore intraspecific lineage divergence in E. hirsuticauda at a biogeographic break. We analysed genetic variation of mitochondrial DNA sequences and microsatellite loci, and morphological and morphometric variation along most of the distributional range of this species (25°S - 42°S). Previous studies in marine peracarids that incorporate morphological analyses in addition to the genetic counterpart suggest that morphology usually correlates with genetic divergence [45,46,47]. Using genetic and morphological divergence, herein we test the hypothesis that E. hirsuticauda, a species with limited dispersal, has a phylogeographic break coincident with the 30°S biogeographic break, and that the degree of genetic and morphologic divergence represents an advanced stage of speciation at 30°S.
The analysed Cytochrome Oxidase I (COI) sequence dataset was composed of 404 sequences of 600 base pairs of length from 14 localities, which corresponded to 128 distinct haplotypes (Table 1, Additional file 1). The inferred best model of molecular evolution for COI was the GTR + G + I model. Phylogenetic reconstructions supported the monophyly of E. hirsuticauda (Fig. 1a) and revealed two reciprocally monophyletic intraspecific lineages; one to the north and one to the south of 30°S (Fig. 1b). A third group of haplotypes, albeit not reciprocally monophyletic, was detected embedded within the clade located to the south of 30°S (Fig. 1b). The haplotype network revealed the same two main clades, separated by at least 24 mutational steps (Fig. 1c). To the south of 30°S there was the third group of divergent haplotypes occurring mainly in the central study area, between 30°S and 35°S (Fig. 1c, d). Given that most analyses were performed considering these groups, we will refer to them as: North group (north of 30°S), Center group (between 30 and 35°S) and the South group (to the south of 35°S).
Most population pairwise ΦST values were high and significant (Additional file 2). Analysis of Molecular Variance (AMOVA) determined that 78.06% of the genetic variance was explained by a priori defined groups (North, Center and South) (Table 2). Migration rates from IMa2 analyses revealed that the North and Center groups were strongly isolated and started diverging ~ 1.4 million years ago. The Center and South groups started diverging more recently, ~ 72,000 years ago (Table 3, Additional file 3). No migration was detected between North and Center groups and very low and asymmetrical migration between Center and South groups, with a slight southern direction, but not statistically different from the null model of no migration (Table 3, Additional file 3).
Genetic diversity differed between groups. The North group had low haplotype richness, low nucleotide diversity, and a low number of shared haplotypes. The Center group had high haplotype richness, high nucleotide diversity, and numerous shared haplotypes, while the South group had high haplotype richness, low nucleotide diversity, and low number of shared haplotypes (Table 1, Fig. 1d). Likewise, indicators of departures of DNA variability from neutrality processes expectations of the three groups differed. Tajima’s D and Fu & Li’s F (Table 1), were mostly significant in the North and South groups and not-significant in the Center group. Tajima’s D and Fu & Li’s F, calculated per group including all sites, showed that the North and South groups had highly significant values while the Center clade was not-significant (Table 4). Mismatch frequency distributions of the number of differences between pairs of sequences had higher significance in the North and South groups (Table 1). Temporal reconstructions of effective population depicted as Bayesian Skyline plots, inferred using HKY + I and TN93+ G models respectively, were similar in shape and indicative of a past population expansion for the North and South groups (Additional file 4). Altogether data show that the North and South groups underwent a past demographic expansion event. Based on a τ = 4.316 (Table 4), the estimated time since the last population expansion of the North group was 240,000 years ago, while the South group had an earlier expansion, dating back to 370,000 years ago based on τ = 6.711 (Table 4). These estimations were in agreement with the median of the estimated date for the Most Recent Common Ancestor (tMRCA) calculations carried out in BEAST. Estimations located the tMRCA of the North group at 240,000 years ago and of the South group at 340,000 years ago approximately (Additional file 4).
Eight loci were genotyped for 317 individuals of eight localities (Additional file 5), however, three loci (Ehir1, Ehir6 and Ehir49) were excluded from further analyses after corroborating that they were in significant linkage disequilibrium and deviated from Hardy-Weinberg Equilibrium (HWE) across all populations (data not shown).
The Brookfield method in MICRO-CHECKER found no scoring errors or linkage disequilibrium for the five remaining loci (Additional file 6). HWE analyses using the 5 remaining loci at the population level indicated that most localities were out of HWE. Inspection by locus and population shows that loci had less than 50% of localities out of HWE (Additional file 7). Loci were highly polymorphic, with a total Polymorphic Information Content (PIC) of 0.810, and with similar observed and expected heterozygosities, ranging from 0.438 to 0.855 and 0.467 to 0.923, respectively (Additional file 8). Overall, locus-by-locus analyses showed non-significant FIS values, while population-by-population analyses showed significant FIS values, indicating population processes may be determining the heterozygosities (Additional file 9).
Individual-based Bayesian clustering analysis without a priori information of origin location, performed in STRUCTURE, grouped genotypes in three population clusters (Fig. 2) that corresponded to the three groups detected with COI sequences (North of 30°S, Center between 30 and 35°S, and South of 35°S) (Fig. 2). Population pairwise FST values were mostly significant (Additional file 10). AMOVA for three groups (the same ones tested for COI data) revealed that most of the molecular variance of microsatellite loci was explained by groups (79.39%) (Table 5).
Morphology and morphometrics
A total of 3362 isopods were measured from seven sites (Additional file 11). Total body length ranged from 2.5 mm to 13.2 mm (Fig. 3a, Additional file 12). Individuals of the North group were smaller than individuals belonging to the Center and South groups (P < 0.01). Total body length was significantly different both between the three groups, with groups explaining 61% of variation in total body length, and between sites within groups, although sites explained only 12% of the variation (Nested ANOVA, F: 211.2; df: 4 74.36; P < 0.01) (Table 6).
Principal component analysis (PCA) of adjusted morphometric data of males revealed a sharp difference between the North, Center and South groups (Fig. 3b) that was significant according to Hotelling’s discriminant test (P < 0.01) (Fig. 3c). The first principal component (PC) explained 55.1% of the variation and combined with PC 2 they explained 77.3% of the variation. In PC 1, the character with the highest contribution was the lateral projection of the appendix masculina (− 0.742) (Additional file 13). The analysis was redone after eliminating the discrete variables, but the results were the same (data not shown). Permutation Analysis of Variance (PERMANOVA) performed with the three detected groups revealed significant differences between them (df = 2; Pseudo-F = 39.872; P(perm) = 0.0001). Pairwise comparisons showed the same pattern between North-Center groups (P(perm) = 0. 0001), North-South (P(perm) = 0.0001) and with a lower distance, between Center-South groups (P(perm) = 0.0023) (Table 7).
Analysis of Covariance (ANCOVA) on Log10-transformed variables evaluated on selected morphometric characters demonstrated that, accounting for body length, the lateral projection of the appendix masculina was the only character that was different between groups (Fig. 3d); other slopes showed no differences between groups. Individuals belonging to the northern group showed relatively larger projections than individuals from the southern groups, despite their smaller body length (ANCOVA on Log10-transformed data, F: 62.28; df: 1; P < 0.01) (Table 8).
Divergence of Excirolana hirsuticauda
Phylogeographic and morphological analyses of the beach-dwelling isopod E. hirsuticauda along the southeast Pacific coast showed concordant genetic and morphological divergence at 30°S and at 35°S. The strongest divergence, both with genetic and morphological markers, was detected across the 30°S biogeographic break, where the degree of divergence is congruent with a very advanced stage of the speciation process, that started over a million years ago. The significant differentiation of a reproductive character suggests that there may be reproductive isolation and complete speciation to the north and south of 30°S.
Intraspecific lineage divergence at 30°S has been likely maintained by the oceanographic discontinuities at 30°S that promote ecological shifts [5, 8, 11] and divergent selection with local adaptation . Another variable that may enhance divergence is the detected habitat discontinuity across the study area [50,51,52]. Habitat connectivity values, determined from analysis of the length and distance between sandy beaches distributed between 25°S and 42°S, showed lower habitat connectivity at ~ 30–31°S, with a value up to eight orders of magnitude lower than in the rest and with overall less and lower values (Additional file 14). Indeed, there are very few and only small sandy beaches between 30 and 31°S (see also ), creating an important habitat discontinuity. Habitat discontinuity has been invoked to explain phylogeographic breaks in other marine organisms such as algae [53,54,55], which have very low dispersal potential, and also for the isopod Jaera albifrons (38). Further studies that directly evaluate habitat connectivity along the HCS will allow determining the contribution of habitat continuity in the phylogeographic structure of marine taxa at 30°S. The major constrains to gene flow detected herein are likely a combination of low habitat connectivity, oceanographic discontinuities [2,3,4,5,6,7,8,9,10,11,12], the very limited dispersal potential and small body size of E. hirsuticauda.
A second discontinuity in population divergence of E. hirsuticauda was detected between 34°23′S and 36°26′S (~ 35°S). Intraspecific divergence at 35°S was weaker than at 30°S, and started more recently (less than 100,000 years ago). At the 30°S biogeographic break, genetic data detected reciprocally monophyletic groups to the North and South, while at 35°S the divergence did not lead to reciprocal monophyly. In spite of the weaker overall signal, the 35°S break was evident with both genetic and morphological data.
Although the presence of a genetic discontinuity at 35°S in a marine species of the coast of Chile is novel, it is in keeping with species boundaries known to occur at this location. The closest described genetic discontinuity is between 32°S and 34°S for the seaweed Mazzaella laminarioides . Based on the distribution of almost 1000 species, Lancellotti & Vásquez [56, 57] described a discontinuity in species distribution where the northern limit of a transition zone was at 35°S. This transition zone extends south to 41°S and is characterized by gradual species replacement and changing environmental conditions. In particular, short-ranged, warm-temperate species have their southern limit of distribution at ~ 35°S . Analysing peracarid diversity, Rivadeneira et al.  found that species diversity decreases to the north of 35°S. The direct effects of El Niño Southern Oscillation (ENSO), as well as the warm Peruvian Countercurrent that reaches from the south up to 35°S , have been invoked to explain species diversity changes at 35°S [57, 59]. Recently, Lara et al.  reported a biogeographic break at 35°S for marine rocky intertidal species of the study area, particularly species with planktotrophic larvae. They found the break to be consistent with the spatial structure of sea surface temperature, chlorophyll-a, and river outflow. All of the above conforms relevant evidence about the existence of a genetic discontinuity and a biogeographic break at 35°S on the southeast Pacific coast for some marine species.
Direction of population expansion and genetic diversity of genetic groups of E. hirsuticauda
A past population expansion event could be suggested for the North and South, but not the Center group. The demographic expansions of the North and South groups date earlier than the divergence of the Center group from the South group. The expansions date during two Pleistocene interglacial periods, 240,000 and 240,000–370,000 years ago, for the North and South groups, respectively.
The genetic diversity and likely demographic differences between North and South groups are likely outcomes of historic signatures of isolation in glacial refugia during low sea-level glaciation periods of the Pleistocene. The regional persistence of species in refugia has been commonly described in marine species (e.g. [21, 61]). Subsequently, during interglacial periods, there may have been a northern coastal recolonization explaining the more recent population expansion detected for the North clade. Lower diversity in the North clade suggests that the direction of geographic and demographic expansions of E. hirsuticauda were from south to north, similar as suggested for other species expanding from southern glacial refugia [62,63,64]. Gene flow between Center and South groups is low but asymmetric with gene flow occurring only in the south-north direction. In addition to the effects of the direction of colonization, sufficient time has passed to allow drift within populations to enhance differentiation.
Species from sandy beaches such as cirolanid isopods may locally disappear after major changes in beach morphology and height, but some species, including E. hirsuticauda can also rapidly recolonize sandy beaches . Environmental changes such as the ones caused by more intense effects of ENSO at lower latitudes [57, 58] may have also shaped the lower genetic diversity in the north. The effects of ENSO may have led to repeated bottlenecks during El Niño phases , further reducing genetic divergence in the northern area.
Speciation within Excirolana hirsuticauda
The degree of genetic and morphological divergence of E. hirsuticauda to the north and south of 30°S, with reciprocal monophyly and lack of gene flow, suggests an advanced stage of the speciation process occurring at a biogeographic break. The morphological trait that explained most of the divergence at 30°S was the appendix masculina of males. Isopods have internal insemination and sperm transfer is often mediated by the appendix masculina, a structure of the pleopod that shows great variation between species . Detected groups of E. hirsuticauda populations showed significant divergence in the relative length (and not the absolute length) of the projection of the appendix masculina in relation to the rest of the appendix, which should lead (if it has not already) to a prezygotic reproductive isolation, similar as reported for E. braziliensis  and other marine invertebrates (e.g., [69, 70]). Divergence in genital characters is commonly observed in arthropods, likely enforcing reproductive isolation and speciation [71, 72]. Herein we demonstrate that the shape of the appendix masculina of E. hirsuticauda clearly distinguishes and allows the correct assignment of individuals from the north and south of the 30°S biogeographic break. Additionally, the appendix masculina also allows identification of linages separated at 35°S albeit with less divergence than at 30°S. There is consistency in the degree of genetic and morphological divergence in each of the two discontinuities detected, at 30°S and 35°S.
The evolutionary history of E. hirsuticauda led to concordant signals of divergence in genetic diversity and morphology, specifically of a trait related to reproductive isolation, at two geographic areas with differing time since onset of the divergence. This study highlights that in E. hirsuticauda, divergence in the appendix masculina of males, morphometry and genetic diversity co-vary at the same latitudes. Additionally, it exemplifies the result of divergence across biogeographic breaks that can shape species diversity by restricting gene flow. Restricted gene flow resulting from the combination of small body size, mode of development without a pelagic larval stage, oceanographic and habitat discontinuities have likely shaped the divergence processes in E. hirsuticauda along the Humboldt Current System. The evolution of a morphological trait that could be a reproductive barrier suggests that even in the absence of ecological or habitat constrains to gene flow, reproductive isolation would persist. Thus, both restricted gene flow and the prezygotic reproductive barrier were likely factors that promoted speciation within E. hirsuticauda.
Mitochondrial COI gene
Individuals of E. hirsuticauda were collected from each of 14 localities between 25°S and 42°S (Table 1). Isopods were washed out of bulk sand samples using 2 mm mesh bags, and were preserved in 95% ethanol and stored at − 20 °C. Genomic DNA was isolated using the QIAamp DNA kit (Qiagen Inc., Valencia, CA, USA).
We obtained partial mtDNA COI gene sequences using the procedures described by Varela & Haye  using primers HCO and LCO . Inspection and alignment of chromatographs were carried on with the software GENEIOUS R8 .
Phylogenetic analyses were carried out using Maximum Likelihood with MODELTEST in PAUP* 4.0b10 . Support for nodes was estimated with 1000 non-parametric bootstrap replicates.
A median-joining haplotype network was constructed using NETWORK 4.6 (Fluxus Technology, 2010). Number of segregating sites, number of haplotypes, nucleotide diversities, Tajima’s D  and Fu & Li’s F  were calculated using DNASP 5 . Haplotype richness (Rh), after rarefaction to the minimum sample size of N = 20, was estimated using CONTRIB 1.4 .
To detect historical population expansions, mismatch frequency distributions of the number of nucleotide differences between pairs of sequences and Raggedness index were performed in Arlequin 3.5  and Bayesian Skyline plots, to infer population dynamics through time, in BEAST v 2.4.8 [81, 82]. In BEAST, we used an uncorrelated log normal relaxed clock model for each of the detected groups with five independent runs using 100 × 106 generations sampled every 1000 generations. For each group, we performed separately a JModelTest v2  in order to use the best fit substitution model for the demographic reconstruction. Conversed chains were checked with Tracer v 1.6 (http://beast.community/tracer). Tracer v 1.6 was used to depict the Bayesian skyline plots. For this, we used the median values and corresponding 95% HDP (highest density probabilities) confidence intervals for effective population size (Ne) dynamics trough time and the time of the most recent common ancestor (tMRCA) of the demographic expansion for each group.
Pairwise ΦST values were estimated using ARLEQUIN 3.5. ARLEQUIN 3.5 was used to perform AMOVA using a priori groups detected by previous analyses using the pairwise difference as distance method with 10,000 permutations.
To estimate the splitting time (t) of the detected genetic groups and the degree of isolation of each group (migration patterns), we used the isolation-with-migration model implemented in the software IMa2 [84, 85]. For these purposes, we carried out several preliminary runs in the M mode of the software to determinate the best set of priors that ensure mixing and convergence of the Markov Chains (MCMC). Uniform priors were used to estimate splitting time (t = 100), whereas an exponential prior (mean = 1) for gene flow (m) was adopted. We executed 1 × 108 MCMC iterations sampling every 100 generations, with a burn-in period of the first 25%. According to the authors’ recommendation, we assumed that the mutation rate is under a Hasegawa-Kishino-Yano (HKY) model. After achieving convergence in the M-mode, we used the simulated genealogies using the L-Mode (Load Tree mode) to calculate the log maximum-likelihood and credibility intervals (95% under HPD) estimates for migration parameters using a Likelihood Ratio Test (LRT). Splitting times were rescaled into years (t/μ) and effective rate of migration (Θx/ mx)/2) (rate at which genes come into population, per generation) using a mutational rate of 2% per Myr [86, 87].
For microsatellite DNA analyses, eight primer pairs  were used to assess the degree of population differentiation of E. hirsuticauda from eight localities. Amplicons were run on an automated sequencer and allele sizes were scored using GENEMAKER 1.80 (SoftGenetics). MICRO-CHECKER 2.2.3  was used to detect scoring errors and null alleles using the Brookfield 1 method .
GENEPOP 4.0  was used to test for linkage disequilibrium and deviations from HWE expectations within each population and per locus, followed by sequential Bonferroni correction. Allele richness estimates were obtained with a rarefaction analysis using HP-RARE v1.0 . The following analyses were also performed: estimation of heterozygosities and proportion of polymorphic loci (GENALEX ); estimation of inbreeding coefficients FIS (FSTAT 184.108.40.206 ); estimation of number of groups with SAMOVA 1.0. Finally, the Bayesian method of Prichard et al.,  implemented in STRUCTURE 2.3.4 was used to investigate the spatial genetic structure. The most likely value for K was assessed using Evanno’s method  by comparing the likelihood of the data for different values of K. Calculations were conducted from K 1–9 with 10 replications using an ancestry model that incorporates admixture, each one with a burn-in period of 25%, followed by 5 × 105 iterations.
Morphology and morphometrics
At least 200 individuals were collected at each of seven beaches between 26.99°S and 42.66°S. Since a recent study had confirmed that body sizes of E. hirsuticauda increased with latitude , herein we determined the body length for both males and females. For each specimen, total body length from the tip of the rostrum to the posterior end of the pleotelson was measured using a dissecting microscope with a graded ocular (Additional file 15).
In a subset of the collected males (10 per population, but 20 males at the populations just north and south of 30°S), 19 measures of body morphology (Additional file 15) were taken. Sample sizes of 10–25 adults have been determined as adequate to characterize the morphometry of populations of E. braziliensis Weinberg & Starczak. Especially since the appendages involved in sperm transfer might vary among populations (e.g. [72, 98]), we dissected the second pleopod and measured the length of the endopod, appendix masculina, and the lateral projection (i.e. “deep notch” in: Ribetti & Roccatagliata ). Dissected appendages were placed in a drop of glycerine on a microscope slide and covered with a cover slide before taking digital images through a microscope. Images were analysed with the program IMAGE-PRO PLUS 6.0 (Media Cybernetics, Inc.).
Nested ANOVA was used to evaluate changes in body length, using sites nested within three regions. The regions were defined with a multivariate classification tree [99, 100] that was based on measurements, using an Euclidian distance matrix with latitude as predictor. The resulting classification showed high accuracy (pseudo-r2 = 0.78) , defining the presence of three groups coincident with those detected with genetic divergence: North of 30°S, Center (one site at 31°51′S), and South (34°-42°30′S) (Additional file 16).
Principal component analysis (PCA) of morphological variables was used to evaluate geographic differences in morphology, accounting for allometric effects on body size using Thorpe’s method [101, 102]. Hotelling’s discriminant test (T2) was used to determine differences in shape among groups. Analyses were carried out using the software PAST . We examined the allometric relationship between selected morphological variables (i.e. those contributing significantly to the PCA) and body length (Log-transformed) using an ANCOVA, with individuals pooled according to the three detected groups. Analyses were carried out using the software R .
PERMANOVA [48, 105] analyses were performed in order to identify significant differences between the detected groups in the PCA. The analyses were based on a Euclidian distance matrix type III (partial) and the sum of squares were calculated with 9999 permutations in the original data set.
Availability of data and materials
Sequences were deposited in GenBank with accession numbers FJ532097-FJ532224. All data generated or analysed during this study are included in this published article [and its Additional files].
Analysis of Molecular Variance
Analysis of Covariance
Analysis of Variance
Cytochrome Oxidase I
El Niño South Oscillation
Humboldt Current System
Likelihood Ratio Test
Markov Chain Monte Carlo
Principal Component Analysis
Permutation Analysis of Variance
Polymorphic Information Content
Spatial Analysis of Molecular Variance
Time of the Most Recent Common Ancestor
Camus PA. Biogeografía marina de Chile continental. Rev Chil Hist Nat. 2001;74:587–617.
Thiel M, Macaya E, Acuña E, et al. The Humboldt current system of northern-Central Chile. Oceanographic processes, ecological interactions and socio-economic feedback. Oceanogr Mar Biol. 2007;45:195–344.
Strub PT, Mesias JM, Montecino V, Rutllant J, Salinas S. Coastal Ocean circulation off western South America. The sea, vol. 11; 1998. p. 273–313.
Marín VH, Delgado LE, Luna-Jorquera G. S-chlorophyll squirts at 30°S off the Chilean coast (eastern South Pacific): feature-tracking analysis. J Geophys Res. 2003;108(C12):3378.
Aravena G, Broitman B, Stenseth MC. Twelve years of change in coastal upwelling along the central-northern coast of Chile: spatially heterogeneous responses to climate variability. PLoS One. 2014;9:e90276.
Hormazábal S, Shaffer G, Leth O. Coastal transition zone off Chile. J Geophys Res. 2004;109:C01021.
Tapia FJ, Largier JL, Castillo M, Wieters EA, Navarrete SA. Latitudinal discontinuity in thermal conditions along the nearshore of Central-Northern Chile. PLoS One. 2014;9:e110841.
Broitman B, Navarrete S, Smith F, Gaines S. Geographic variation of southeastern Pacific intertidal communities. Mar Ecol Prog Ser. 2001;224:21–34.
Rivadeneira MM, Fernandez M, Navarrete SA. Latitudinal trends of species diversity in rocky intertidal herbivore assemblages: spatial scale and the relationship between local and regional species richness. Mar Ecol Prog Ser. 2002;245:123–31.
Navarrete SA, Broitman BR, Menge BA. Interhemispheric comparison of recruitment to rocky intertidal communities: pattern persistence and scales of variation. Ecology. 2008;89:1308–22.
Navarrete SA, Wieters EA, Broitman BR, Castilla JC. Scales of benthic-pelagic coupling and the intensity of species interactions: from recruitment limitation to top-down control. Proc Natl Acad Sci. 2005;102:18046–51.
Valdivia N, Aguilera MA, Navarrete SA, Broitman BR. Disentangling the effects of propagule supply and environmental filtering on the spatial structure of a rocky shore metacommunity. Mar Ecol Prog Ser. 2015;538:67–79.
Bowen BW, Gaither MR, DiBattista JD, Iacchei M, Andrews KR, Grant WS, Toonen RJ, Briggs JC. Comparative phylogeography of the ocean planet. Proc Natl Acad Sci. 2016;113:7962–9.
Marie AD, Lejeusne C, Karapatsiou E, Cuesta JA, Drake P, Macpherson E, Bernatchez L, Rico C. Implications for management and conservation of the population genetic structure of the wedge clam Donax trunculus across two biogeographic boundaries. Sci Rep. 2016;6:39152.
Avise JC, Arnold J, Ball RM, et al. Intraspecific phylogeography: the mitochondrial DNA bridge between population genetics and systematics. Ann Rev Ecol Syst. 1987;18:489–522.
Kuo C-H, Avise JC. Phylogeographic breaks in low-dispersal species: the emergence of concordance across gene trees. Genetica. 2005;124:179–86.
Dawson MN. Phylogeography in coastal marine animals: a solution from California? J Biogeogr. 2001;28:723–36.
Wares JP. Community genetics in the northwestern Atlantic intertidal. Mol Ecol. 2002;11:1131–44.
Dawson MN. Incipient speciation of Catostylus mosaicus (Scyphozoa, Rhizostomeae, Catostylidae), comparative phylogeography and biogeography in Southeast Australia. J Biogeogr. 2005;32:515–33.
Pelc RA, Warner RR, Gaines SD. Geographical patterns of genetic structure in marine species with contrasting life histories. J Biogeogr. 2009;36:1881–90.
Teske PR, Papadopoulos I, Mmonwa KL, Matumba TG, McQuaid CD, Barker NP, Behengaray LB. Climate-driven genetic divergence of limpets with different life histories across a southeast African marine biogeographic disjunction: a different process, same outcome. Mol Ecol. 2011;20:5025–41.
Altman S, Robinson JD, Pringle JM, Byers JE, Wares JP. Edges and overlaps in Northwest Atlantic phylogeography. Diversity. 2013;5:263–75.
Eberl R, Mateos M, Grosberg RK, Santamaria CA, Hurtado LA. Phylogeography of the supralittoral isopod Ligia occidentalis around the point conception marine biogeographical boundary. J Biogeogr. 2013;40:2361–72.
Haye PA, Segovia NI, Muñoz-Herrera NC, Gálvez FE, Martínez A, Meynard A, Pardo-Gandarillas MC, Poulin E, Faugeron S. Phylogeographic structure in benthic marine invertebrates of the Southeast Pacific coast of Chile with differing dispersal potential. PLoS One. 2014;9:e88613.
Hellberg ME, Burton RS, Neigel JE, Palumbi SR. Genetic assessment of connectivity among marine populations. Bull Mar Sci. 2002;70:273–90.
Irwin DE. Phylogeographic breaks without geographic barriers to gene flow. Evolution. 2002;56:2383–94.
Sherman C, Hunt A, Ayre D. Is life history a barrier to dispersal? Contrasting patterns of genetic differentiation along an oceanographically complex coast. Biol J Linn Soc. 2008;95:106–16.
Montecinos A, Broitman BR, Faugeron S, Haye PA, Tellier F, Guillermin M-L. Species replacement along a linear coast habitat: phylogeography and speciation in the red alga Mazzaella laminarioides along the south east pacific. BMC Evol Biol. 2012;12:97.
Zakas C, Binford J, Navarrete SA, Wares JP. Restricted gene flow in Chilean barnacles reflects an oceanographic and biogeographic transition zone. Mar Ecol Prog Ser. 2009;394:165–77.
Ewers-Saucedo C, Pringle JM, Sepúlveda HH, Byers JE, Navarrete SA, Wares JP. The oceanic concordance of phylogeography and biogeography: a case study in Notochthamalus. Ecol Evol. 2016;6:4403–20.
Sánchez R, Sepúlveda R, Brante A, Cárdenas L. Spatial pattern of genetic and morphological diversity in the direct developer Acanthina monodon (Gastropoda: Mollusca). Mar Ecol Prog Ser. 2011;434:121–31.
Nosil P, Vines TH, Perspective FDJ. Reproductive isolation caused by natural selection against immigrants from divergent habitats. Evolution. 2005;59:705–19.
Wu CI. The genic view of the process of speciation. J Evol Biol. 2001;14:851–65.
Norris RD, Hull PM. The temporal dimension of marine speciation. Evol Ecol. 2012;26:393–415.
Villamor A, Constantini F, Abbiati M. Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates. PLoS One. 2014;9:e101135.
Nosil P, Harmon LJ, Seehausen O. Ecological explanations for (incomplete) speciation. Trends Ecol Evol. 2009;24:145–56.
Martin P, Zuccarello GC. Molecular phylogeny and timing of radiation in Lessonia (Phaeophyceae, Laminariales). Phycol Res. 2012;60:276–87.
Carvalho G, Piertney S. Interspecific comparison of genetic population structure in members of the Jaera albifrons species complex. J Mar Biol Assoc UK. 1997;77:77–93.
Wilson A, Boates J, Snyder M. Genetic isolation of populations of the gammaridean amphipod, Corophium volutator, in the bay of Fundy, Canada. Mol Ecol. 1997;6:917–23.
De Matthaeis E, Davolos D, Cobolli M, Ketmaier V. Isolation by distance in equilibrium and nonequilibrium populations of four talitrid species in the Mediterranean Sea. Evolution. 2000;54:1606–13.
Baratti M, Goti E, Messana G. High level of genetic differentiation in the marine isopod Sphaeroma terebrans (Crustacea: Isopoda: Sphaeromatidae) as inferred by mitochondrial DNA analysis. J Exp Mar Biol Ecol. 2005;315:225–34.
Ketmaier V, De Matthaeis E, Fanini L, Rossano C, Scapini F. Variation of genetic and behavioural traits in the sandhopper Talitrus saltator (Crustacea Amphipoda) along a dynamic sand beach. Ethol Ecol Evol. 2010;22:17–35.
Varela AI, Haye PA. The marine brooder Excirolana braziliensis (Crustacea: Isopoda) is also a complex of cryptic species on the coast of Chile. Rev Chil Hist Nat. 2012;85:495–502.
Hurtado LA, Mateos M, Mattos G, Liu S, Haye PA, Paiva PC. Multiple transisthmian divergences, extensive cryptic diversity, occasional long-distance dispersal, and biogeographic patterns in a marine coastal isopod with an amphi-American distribution. Ecol Evol. 2016;6:7794–808.
Lessios HA, Weinberg JR, Starczak VR. Temporal variation in populations of the marine isopod Excirolana: how stable are gene frequencies and morphology? Evolution. 1994;48:549–63.
Liu J, Li Q, Kong L, Zheng X. Cryptic diversity in the pen shell Atrina pectinata (Bivalvia: Pinnidae): high evidence and hybridization revealed by molecular and morphological data. Mol Ecol. 2001;20:4332–45.
Yin J, Pan D, He C, Wang A, Yan J, Sun H. Morphological and molecular data confirm species assignment and dispersal of the genus Ligia (Crustacea: Isopoda: Ligiidae) along the northeastern coastal China and East Asia. Zool J Linnean Soc. 2013;169:362–76.
Anderson MJ. Permutation test for univariate or multivariate analysis of variance and regression. Can J Fish Aquat Sci. 2001;58:626–39.
Segovia NI, Gallardo-Escárate C, Poulin E, Haye PA. Lineage divergence, local adaptation across a biogeographic break, and artificial transport, shape the genetic structure in the ascidian Pyura chilensis. Sci Rep. 2017;7:44559.
Teske PR, Papadopoulus I, Newman BK, Dworschak PC, MacQuaid CD, Barker N. Oceanic dispersal barriers, adaptation and larval retention: an interdisciplinary assessment of potential factors maintaining a phylogeographic break between sister lineages of an African prawn. BMC Evol Biol. 2008;8:341.
Fraser C, Thiel M, Spencer H, Waters J. Contemporary habitat discontinuity and historic glacial ice drive genetic divergence in Chilean kelp. BMC Evol Biol. 2010;10:203.
Wang J, Tsang LM, Dong YW. Causations of phylogeographic barrier of some rocky shore species along the Chinese coastline. BMC Evol Biol. 2015;14:114.
Alberto F, Raimondi PT, Reed DC, et al. Habitat continuity and geographic distance predict population genetic differentiation in a giant kelp. Ecology. 2010;91:49–56.
Buonomo R, Assis J, Fernandes F, Engelen AH, Airoldi L, Serrao EA. Habitat continuity and stepping-stone oceanographic distances explain population genetic connectivity of the brown alga Cystoseira amentacea. Mol Ecol. 2017;26:766–80.
Durrant HMS, Barret NS, Edgar GJ, Coleman MA, Burriedge CP. Seascape habitat patchiness and hydridu¡ynamics explain genetic structuring of kelp populations. Mar Ecol Prog Ser. 2018;587:81–92.
Lancellotti DA, Vásquez JA. Biogeographical patterns of benthic invertebrates in the southeastern Pacific littoral. J Biogeogr. 1999;26:1001–6.
Lancellotti DA, Vásquez JA. Zoogeografía de macroinvertebrados bentónicos de la costa de Chile: Contribución para la conservación marina. Rev Chil Hist Nat. 2000;73:99–129.
Brattström H, Johanssen A. Ecological and regional zoogeography of the marine benthic fauna of Chile. Sarsia. 1983;68:289–339.
Rivadeneira MM, Thiel M, González ER, Haye PA. An inverse latitudinal gradient of diversity of peracarid crustaceans along the Pacific coast of South America: out of the deep south. Glob Ecol Biogeogr. 2011;20:437–48.
Lara C, Saldías GS, Cazelles B, Rivadeneira MM, Haye PA, Broitman BR. Coastal biophysical processes and the biogeography of rocky intertidal species along the south-East Pacific. J Biogeogr. 2019;46:420–31.
Haye PA, Muñoz-Herrera NC. Isolation with differentiation followed by expansion with admixture in the tunicate Pyura chilensis. BMC Evol Biol. 2013;13:252.
Provan J, Bennett KD. Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008;23:564–71.
Hellberg ME, Balch DP, Roy K. Climate-driven range expansion and morphological evolution in a marine gastropod. Science. 2001;292:1707–10.
Cox LN, Zaslavskaya NI, Marko PB. Phylogeography and trans-Pacific divergence of the rocky shore gastropod Nucella lima. J Biogeogr. 2014;41:615–27.
Jaramillo E, Dugan JE, Hubbard DM, Melnick D, Manzano M, Duarte C, Sánchez R. Ecological implications of extreme events: footprints of the 2010 earthquake along the Chilean coast. PLoS One. 2012;7:e35348.
Faugeron S, Martinez E, Correa J, Billot C. Long-term copper mine waste disposal in northern Chile associated with gene flow disruption of the intertidal kelp Lessonia nigrescens. Mar Ecol Prog Ser. 2005;288:129–40.
GDF W. Functional morphology and evolution of isopod genitalia, Columbia University press. Crustacean Sexual Biology, New York/Oxford. 1991;228245.
Weinberg JR, Starczak VR. Morphological differences and low dispersal between local populations of the tropical beach isopod, Excirolana braziliensis. Bull Mar Sci. 1988;42:296–309.
Prada C, Hellberg ME. Long prereproductive selection and divergence by depth in a Caribbean candelabrum coral. Proc Natl Acad Sci. 2013;110:3961–6.
Klibansky LKJ, McCartney MA. Conspecific sperm precedence is a reproductive barrier between free-spawning marine mussels in the Northwest Atlantic Mytilus hybrid zone. PLoS One. 2014;9:e108433.
Simmons LW. Sexual selection and genital evolution. Austral Entomol. 2014;53:1–17.
Mathews LM, Adams L, Anderson E, Basile M, Gottardi E, Buckholt MA. Genetic and morphological evidence for substantial hidden biodiversity in a freshwater crayfish species complex. Mol Phylogenet Evol. 2008;48:126–35.
Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3:294–9.
Drummond AJ, Ashton B, Buxton S, Cheung M, Cooper A, Heled J, Kearse M, Markowitz S, Moir R, Stones-Havas S, Sturrock S, Thierer T, Wilson A. GENEIOUS. 2010; version 5.4. [http://www.geneious.com].
Swofford DL. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). 2002. Sunderland, MA, USA: Sinauer Associates.
Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1998;123:585–95.
Fu YX, Li WH. Statistical tests of neutrality of mutations. Genetics. 1993;133:693–709.
Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.
Petit RJ, El Mousadik A, Pons O. Identifying populations for conservation on the basis of genetic markers. Conserv Biol. 1998;12:844–55.
Excoffier L, Lischer HEL. ARLEQUIN version 3.5: a new series of programs to perform population genetics analyses under Linux and windows. Mol Ecol Res. 2010;10:564–7.
Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.
Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.
Darriba D, Taboada G, Doallo R, Posada D. jModelTest 2: more models, new 612 heuristics and parallel computing. Nat Methods. 2012;9:772.
Nielsen R, Wakeley J. Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics. 2001;158:885–96.
Hey J, Nielsen R. Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proc Natl Acad Sci. 2007;104:2785–90.
Poulakakis N, Sfenthourakis S. Molecular phylogeny and phylogeography of the Greek populations of the genus Orthometopon (Isopoda, Oniscidea) based on mitochondrial DNA sequences. Zool J Linnean Soc. 2008;152:707–15.
Malay MC, Paulay G. Peripatric speciation drives diversification and distributional pattern of reef hermit crabs (Decapoda: Diogenidae: Calcinus). Evolution. 2010;64:634–62.
Haye PA, Marchant S. Isolation and characterization of eleven PCR primers for microsatellite loci for the Chilean marine isopod Excirolana hirsuticauda. Mol Ecol Res. 2008;8:1088–90.
van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004;4:535–8.
Brookfield J. A simple new method for estimating null allele frequency from heterozygote deficiency. Mol Ecol. 1996;5:453–5.
Rousset F. Genepop 007: a complete reimplementation of the Genepop software for windows and Linux. Mol Ecol Res. 2008;8:103–6.
Kalinowski ST. HP-RARE 1.0: a computer program for performing rarefaction on measures of allelic richness. Mol Ecol Notes. 2005;5:187–9.
Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in excel. Population genetic software for teaching and research-an update. Bioinformatics. 2012;28:2537–9.
Goudet J. FSTAT: a program to estimate and test gene diversities and fixation indices. 2002; Version 220.127.116.11. Available from URL: https://www2.unil.ch/popgen/softwares/fstat.htm.
Prichard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.
Jaramillo E, Dugan JE, Hubbard DM, Contreras H, Duarte C, Acuña E, Schoeman DS. Macroscale patterns in body size of intertidal crustaceans provide insights on climate change effects. PLoS One. 2017;12:e0177116.
Ribetti V, Roccatagliata D. Redescription of Excirolana armata (Dana, 1853) and synonymy of Cirolana argentina Giambiagi, 1930 (Crustacea: Isopoda: Cirolanidae). Proc Biol Soc Wash. 2006;119(1):14.
Cutler DR, Edwards TC, Beard KH, Cutler A, Hess KT, Gibson J, Lawler JJ. Random Forest for classification in ecology. Ecology. 2007;88:2783–92.
Fenberg PB, Menge B, Raimondi P, Rivadeneira MM. Biogeographic structure of the northeastern Pacific rocky intertidal: the role of upwelling and dispersal to drive patterns. Ecography. 2015;37:1–13.
Thorpe R. Quantitative handling of characters useful in Snake systematics with particular reference to intraspecific variation in the ringed snake Natrix natrix (L.). Biol J Linn Soc. 1975;7:27–43.
Reist J. An empirical evaluation of several univariate methods that adjust for size variation in morphometric data. Can J Zool. 1984;63:1429–39.
Hammer Ø, Harper DAT, Ryan PD. PAST: paleontological statistics software package for education and data analysis. Palaeontol Electronica. 2001;4:4–9.
R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2018. URL http://www.R-project.org/
Anderson MJ, Gorley RN. Clarke KR. PERMANOVA + for PRIMER: guide to software and statistical methods. UK: PRIMER-E, Plymouth; 2008.
The authors would like to thank Iván Hinojosa and Leonardo Miranda for help with collection of individuals of Excirolana hirsuticauda along the coast of Chile. The authors would also like to thank Raúl Vera, Natalia Muñoz-Herrera and Sergio Marchant for assistance in laboratory procedures.
This research was primarily funded by FONDECYT grants 1051076, 1090670, and 1140862.
Ethics approval and consent to participate
Procedures of collecting and post-collection treatment of samples was approved by ethics committee of the Universidad Católica del Norte.
The author(s) declare(s) that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Raw COI haplotype frequencies of the 128 haplotypes of Excirolana hirsuticauda. (XLSX 44 kb)
Pairwise population ΦST values for COI mtDNA of Excirolana hirsuticauda among 14 sites. (DOCX 79 kb)
Marginal posterior probability distributions of migration rates and time since splitting between groups of Excirolana hirsuticauda using IMa2 with COI sequences. (DOCX 714 kb)
Demographic reconstructions in COI sequences of Excirolana hirsuticauda using BEAST. (DOCX 183 kb)
Microsatellite data of Excirolana hirsuticauda in GeneAlex format. (CSV 22 kb)
Frequency of null alleles per locus and population for microsatellite loci of Excirolana hirsuticauda. (DOCX 80 kb)
Locus-by-locus probability of deviation of Hardy-Weinberg Equilibrium after Bonferroni correction. (DOCX 72 kb)
Summary statistics for microsatellite loci of Excirolana hirsuticauda. (DOCX 117 kb)
FIS values for microsatellite loci of Excirolana hirsuticauda. (DOCX 69 kb)
Microsatellite pairwise population FST of Excirolana hirsuticauda. (DOCX 61 kb)
Data of samples used for body length and morphometric analysis. (DOCX 62 kb)
Body length range distribution of Excirolana hirsuticauda. (DOCX 137 kb)
Eigenvectors of morphometric analysis. (DOCX 97 kb)
Connectivity of sandy beaches across the study area based on beach length and beach separation. (DOCX 190 kb)
List of the 19 morphometric measurements analysed in males of Excirolana hirsuticauda. (DOCX 89 kb)
Multivariate classification tree based on morphological traits of Excirolana hirsuticauda. (DOCX 106 kb)