Genetic and morphological divergence at a biogeographic break in the beach-dwelling brooder Excirolana hirsuticauda Menzies (Crustacea, Peracarida)

Background 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. Results 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. Conclusions 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. Electronic supplementary material The online version of this article (10.1186/s12862-019-1442-z) contains supplementary material, which is available to authorized users.


Background
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 [1]. 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 [3]. 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 [15]. Concordance in the location of intraspecific genetic discontinuities suggest that there are common and persistent factors that promote intraspecific genetic structure [16], 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 [28], intertidal barnacles [29,30], gastropods [24,31], and the amphipod Orchestoidea tuberculata [24]. 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 [20].
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 [32]. 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 [36]. 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 [36]. 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 [24]. 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 phyloand bio-geographic break a natural scenario for the evaluation of intraspecific divergence processes.
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.

Mitochondrial DNA
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) (    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 notsignificant (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).

Microsatellite loci
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 F IS values, while population-by-population analyses showed significant F IS 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 F ST 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 Log 10 -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 Log 10 -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 [49]. Another variable that may enhance divergence is the detected habitat discontinuity across the study area [50][51][52]. Habitat  (see also [2]), 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  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 [28]. 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 at3 5°S [58]. Analysing peracarid diversity, Rivadeneira et al. [59] 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 [58], have been invoked to explain species diversity changes at 35°S [57,59]. Recently, Lara et al. [60] 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 [65]. 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 [66], 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 [67]. 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 [68] 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.

Conclusion
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 covary 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.

Methods
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 [43] using primers HCO and LCO [73]. Inspection and alignment of chromatographs were carried on with the software GENEIOUS R8 [74].
Phylogenetic analyses were carried out using Maximum Likelihood with MODELTEST in PAUP* 4.0b10 [75]. Support for nodes was estimated with 1000 nonparametric bootstrap replicates.
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 [80] 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 × 10 6 generations sampled every 1000 generations. For each group, we performed separately a JModelTest v2 [83] 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 (N e ) 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 isolationwith-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 × 10 8 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 / m x )/2) (rate at which genes come into population, per generation) using a mutational rate of 2% per Myr [86,87].

Microsatellite loci
For microsatellite DNA analyses, eight primer pairs [88] 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 [89] was used to detect scoring errors and null alleles using the Brookfield 1 method [90]. GENEPOP 4.0 [91] 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 [92]. The following analyses were also performed: estimation of heterozygosities and proportion of polymorphic loci (GENALEX [93]); estimation of inbreeding coefficients F IS (F STAT 2.9.3.2 [94]); estimation of number of groups with SAMOVA 1.0. Finally, the Bayesian method of Prichard et al., [95] 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 [96] 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 × 10 5 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 [97], 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 [98]). 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 (pseudor 2 = 0.78) [99], 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 (T 2 ) was used to determine differences in shape among groups. Analyses were carried out using the software PAST [103]. 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 [104].