Skip to main content

Advertisement

  • Research article
  • Open Access

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

BMC Evolutionary Biology201919:118

https://doi.org/10.1186/s12862-019-1442-z

  • Received: 1 February 2019
  • Accepted: 27 May 2019
  • Published:

Abstract

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.

Keywords

  • Southeast Pacific
  • Chile
  • Phylogeography
  • Population genetics
  • Speciation
  • Genital characters
  • Morphological variation

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 [812].

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 [1724]. 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, 2427].

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 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 [3842] and phylogeographic discontinuities coincident with biogeographic breaks [23]. 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 [4547]. 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.

Results

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).
Table 1

Summary of samples of Excirolana hirsuticauda used for COI analyses, and genetic diversity and neutrality tests for COI sequences of Excirolana hirsuticauda

    

Diversity indices

Neutrality tests

Demo. expansion

Geo. expansion

  

Latitude

N

S

H

Rh

π

k

D

F

P (SSD)

P (Hr)

P (SSD)

P (Hr)

Taltal

TAL

25°70′S

38

6

7

4.303

0.0010

0.609

−1.548

−4.280

0.532

0.505

0.319

0.533

Caldera

CAD

26°59′S

20

3

4

3.000

0.0005

0.300

−1.723

−2.749

0.455

0.595

0.450

0.611

Playa Blanca

PBL

28°11′S

36

12

10

5.000

0.0011

0.667

−2.419

−8.953

0.333

0.514

0.168

0.517

Coquimbo

COQ

29°54′S

20

9

5

4.000

0.0017

0.989

−2.098

−1.142

0.198

0.450

0.569

0.655

Los Vilos

LVI

31°51′S

38

21

11

6.841

0.0078

4.700

−0.200

0.172

0.449

0.614

0.573

0.856

Maitencillo

MAI

32°38′S

39

27

16

9.848

0.0087

5.225

−0.622

−2.640

0.009

0.001

0.060

0.081

Pichilemu

PMU

34°23′S

20

16

12

11.000

0.0071

4.242

−0.221

−3.411

0.291

0.418

0.192

0.658

Pangua

PAN

34°29′S

36

25

14

8.174

0.0068

4.073

−1.118

−2.725

0.163

0.166

0.289

0.564

Purema

PUR

36°26′S

20

18

12

11.000

0.0036

2.142

−2.174

−7.309

0.459

0.593

0.498

0.674

Tranaquepe

TRA

38°10′S

36

26

18

10.460

0.0041

2.462

−2.101

−11.785

0.881

0.981

0.884

0.981

Queule

QUE

39°22′S

20

18

11

10.000

0.0040

2.384

−1.994

−5.080

0.572

0.862

0.610

0.861

Calfuco

CAF

39°46′S

25

14

8

5.600

0.0039

2.320

−1.299

−1.159

0.001

0.999

0.713

0.782

Puñihuil

PUÑ

41°55′S

20

17

12

11.000

0.0028

1.700

−2.414

−8.927

0.968

0.894

0.976

0.890

Cucao

CUC

42°40′S

36

25

22

12.175

0.0039

2.333

−2.113

−20.55

0.360

0.554

0.292

0.575

Total

  

404

114

128

13.757

0.0041

2.439

      

Genetic diversity measures for 14 local populations used for COI sequences analysis: Samples size (N), number of segregating sites (S), number of haplotypes (H), haplotype richness after rarefaction (Rh), nucleotide diversity (π), average number of nucleotide differences between pairs of sequences (k). Demographic inference analyses: neutrality tests [Tajima’s D (D) and Fu & Li’s F (F)] (significant values are in bold), probability values associated with Sum of Squares Deviation of the mismatch frequency distribution of the number of pairwise nucleotide differences [P (SSD)] and probabilities of Harpending’s Raggedness index [P (Hr)]. Both SDD and Hr were calculated based on expected mismatch distribution according to the demographic (Demo. expansion) and geographic expansion (Geo. expansion) models

Fig. 1
Fig. 1

Genetic structure based on COI haplotypes of 404 individuals of 14 populations of E. hirsuticauda. (a) Maximum likelihood phylogram of the relationships between some Excirolana species based on COI sequences. Sequences for E. braziliensis and E. chiltoni were obtained from GenBank (accession numbers at tip labels). (b) Maximum likelihood phylogram of COI haplotypes of E. hirsuticauda. Clade in red was unique to the sites to the north of 30°S (marked with red square on map) and showed reciprocal monophyly with sites to the south of 30°S. Grey branches and labels indicate a derived linage, albeit without reciprocal monophyly. In (a) and (b) numbers along the branches represent bootstrap support values. (c) Haplotype network of COI sequences. Circles represent each haplotype with size proportional to each colour representing a locality as designated in the contiguous map. (d) Geographic distribution of haplotypes. Each pie represents a haplotype, and the subdivisions represent relative frequencies of each haplotype per site. Colours represent shared haplotypes according to haplotype network in grey inset box. White portions denote haplotypes not shared with other localities

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).
Table 2

AMOVA analysis for COI sequences of 14 local populations of Excirolana hirsuticauda

 

DF

SS

Var. Comp.

% of Var.

Stats.

Among groups

2

2177.93

7.93

78.06

ΦCT = 0.78

Among pops. Within groups

11

304.02

0.94

9.21

ΦSC = 0.42

Within pops

390

503.89

1.29

12.73

ΦST = 0.87

The three groups considered were the North, Center, and South groups. All analyses yielded P values lower than 0.001. DF degrees of freedom, SS sum of squares, Var. Comp variance component, % of Var percentage of variance explained, Stats statistics

Table 3

Estimates of onset of splitting time (t) across the two detected genetic discontinuities (30° and 35°S) and migration rates (m) from North to South (N > S) and S to N (S > N) for COI sequences of Excirolana hirsuticauda based on the isolation-with-migration model implemented in IMa2

Model

Value

mN > S

mS > N

t (years)

North-Center

HiPt

0.0005

0.0005

1,365,000

Mean

0.0428

0.01815

1,338,333

95%HPD

0.1205

0.0515

4,165,000

LRT

0.000a

0.000a

 

Center-South

HiPt

0.205

0.005

72,083

Mean

0.299

0.0532

73,925

95%HPD

0.675

0.175

98,750

LRT

2.178a

0.000a

 

For each model and parameter, the high point (HP), mean, 95% highest posterior density (95% HPD) of the marginal posterior probabilities are shown, as well as the Likelihood Ratio Test (LRT) performed on migration rates. anon-significant values

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).
Table 4

Demographic inference analyses for COI sequences of 14 local populations of Excirolana hirsuticauda

 

Neutrality Tests

Demographic Expansion

 

D

F

P (SSD)

P (Hr)

τ

North

−1.403

−12.913

0.120

0.099

4.316

Center

−1.148

− 17.989

0.180

0.217

7.961

South

−2.243

−26.169

0.625

0.848

6.711

Neutrality tests [Tajima’s D (D), Fu & Li’s F (F)] (significant values are in bold), probability values associated with Sum of Squares Deviation of the mismatch frequency distribution of the number of pairwise nucleotide differences [P (SSD)] and probabilities of Harpending’s Raggedness index [P (Hr)], and Tau (τ) statistic, under the demographic expansion model

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 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).
Fig. 2
Fig. 2

Bayesian assignment of individuals of E. hirsuticauda to three detected genetic groups based on microsatellite data. (a) Assignment of each individual to genetic clusters is represented with 3 colours. (b) Best-k, optimal number of clusters, according to the estimation based on Evanno’s (Evanno et al. 2005) method

Table 5

AMOVA analysis for microsatellite loci of eight local populations of Excirolana hirsuticauda

 

DF

SS

Var. Comp.

% of Var.

Stats.

Among groups

2

954.87

1.79

79.39

ΦCT = 0.14

Among pops. Within groups

5

61.89

0.16

7.03

ΦSC = 0.08

Within populations

225

133.84

0.31

13.58

ΦST = 0.21

The three groups considered were the North, Center, and South groups. All analyses yielded P values lower than 0.001. DF degrees of freedom, SS sum of squares, Var. Comp. variance component, % of Var percentage of variance explained, Stats. statistics

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).
Fig. 3
Fig. 3

Body size distribution and morphometric differentiation analyses of Excirolana hirsuticauda. (a) Latitudinal variation in total body length, showing the mean (circles) and 95th and 5th percentiles (broken lines). Colours in the graph represent the three groups of populations detected with genetic analyses (North, Center and South). (b) Plot of principal component analysis for 90 males. Symbols indicate from which group individuals were collected. (c) Results of Hotelling’s discriminant test (T2) corroborating the ordination of PCA analysis. (d) ANCOVA results showing the relationship between body size and length of projection of appendix masculina for each group of populations. Insert pictures correspond to detail of the morphology of the appendix masculina of males from North and South groups. Scale bars represent 0.5 mm

Table 6

Summary of nested ANOVA on total body length of 3362 individuals of Excirolana hirsuticauda

Factor

DF

SS

MS

F

P

%Var

Group

3

223.09

74.36

2182.4

< 0.01

61

Group:Site

4

28.78

7.20

211.2

< 0.01

12

DF degrees of freedom, SS sum of squares, MS mean of squares, %Var % of variance explained. Group represents the three groups (regions) as fixed variable, and Group:Site, corresponds to sites nested within regions, with Site as random variable

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).
Table 7

PERMANOVA analysis testing morphological differences between the three detected groups of Excirolana hirsuticauda in the PCA analysis (see Additional file 13)

Model

Source

DF

SS

MS

Pseudo-F

P-value

Total Model

Groups

2

35.696

17.848

39.872

0.0001*

Residual

87

38.943

0.045

  

Total

89

74.639

   

Pairwise Comparisons

Groups

   

t

P-value

N, C

   

68.829

0.0001*

N, S

   

7.373

0.0001*

C, S

   

22.552

0.0023*

Total PERMANOVA model including the three groups (North, Center, and South groups) and pairwise comparisons between groups. Table include degrees of freedom (DF), sum of squares (SS), mean of squares (MS), Pseudo-F ratio and p-values. P-values were calculated with 9999 permutations of the residual model [48]. * significant values

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).
Table 8

Results from ANCOVA on Log10 transformed variables evaluated on four selected morphometric characters of Excirolana hirsuticauda

Relationship

Intercept

 

F

DF

MS

P

TBL / Length of the lateral projection of AM

62.28

1

0.99

< 0.01

TBL/ AM length

3.098

1

0.02

> 0.05

TBL/ Interocular distance

2.75

1

0.004

> 0.05

TBL/ Length of peduncle A2

10.37

1

0.007

< 0.01

Table includes F-ratio (F), degrees of freedom (DF); mean of squares (MS); total body length (TBL); appendix masculina (AM)

Discussion

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 [5052]. 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 [2]), creating an important habitat discontinuity. Habitat discontinuity has been invoked to explain phylogeographic breaks in other marine organisms such as algae [5355], 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 [212], 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 [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 at ~ 35°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 [6264]. 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 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.

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 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 [76] and Fu & Li’s F [77] were calculated using DNASP 5 [78]. Haplotype richness (Rh), after rarefaction to the minimum sample size of N = 20, was estimated using CONTRIB 1.4 [79].

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 × 106 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 (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].

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 FIS (FSTAT 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 × 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 [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 (pseudo-r2 = 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 (T2) 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].

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.

Abbreviations

AMOVA: 

Analysis of Molecular Variance

ANCOVA: 

Analysis of Covariance

ANOVA: 

Analysis of Variance

COI: 

Cytochrome Oxidase I

ENSO: 

El Niño South Oscillation

HCS: 

Humboldt Current System

HKY: 

Hasegawa-Kishino-Yano model

HWE: 

Hardy-Weinberg Equilibrium

LTR: 

Likelihood Ratio Test

MCMC: 

Markov Chain Monte Carlo

Myr: 

Million Years

PCA: 

Principal Component Analysis

PERMANOVA: 

Permutation Analysis of Variance

PIC: 

Polymorphic Information Content

SAMOVA: 

Spatial Analysis of Molecular Variance

tMCRA: 

Time of the Most Recent Common Ancestor

Declarations

Acknowledgements

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.

Funding

This research was primarily funded by FONDECYT grants 1051076, 1090670, and 1140862.

Authors’ contributions

P.A.H. & M.T. conceived the study, M.T. and R.R. obtained the samples, R.R., and A.I.V. performed laboratory work and collected the data. M.M.R., R.R., P.A.H., N.I.S and A.I.V. analysed the data. P.A.H. led the writing. M.T., M.M.R., contributed substantially to the writing process. All authors read and approved the manuscript.

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.

Competing interests

The author(s) declare(s) that they have no competing interests.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Departamento de Biología Marina, Facultad de Ciencias del Mar, Universidad Católica del Norte, Larrondo, 1281 Coquimbo, Chile
(2)
Núcleo Milenio de Ecología y Manejo Sustentable de Islas Oceánicas (ESMOI), Universidad Católica del Norte, Coquimbo, Chile
(3)
Centro de Estudios Avanzados en Zonas Áridas (CEAZA), Coquimbo, Chile

References

  1. Camus PA. Biogeografía marina de Chile continental. Rev Chil Hist Nat. 2001;74:587–617.View ArticleGoogle Scholar
  2. 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.View ArticleGoogle Scholar
  3. 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.Google Scholar
  4. 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.View ArticleGoogle Scholar
  5. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  6. Hormazábal S, Shaffer G, Leth O. Coastal transition zone off Chile. J Geophys Res. 2004;109:C01021.View ArticleGoogle Scholar
  7. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  8. Broitman B, Navarrete S, Smith F, Gaines S. Geographic variation of southeastern Pacific intertidal communities. Mar Ecol Prog Ser. 2001;224:21–34.View ArticleGoogle Scholar
  9. 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.View ArticleGoogle Scholar
  10. 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.PubMedView ArticleGoogle Scholar
  11. 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.PubMedView ArticleGoogle Scholar
  12. 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.View ArticleGoogle Scholar
  13. 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.PubMedView ArticleGoogle Scholar
  14. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  15. 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.View ArticleGoogle Scholar
  16. Kuo C-H, Avise JC. Phylogeographic breaks in low-dispersal species: the emergence of concordance across gene trees. Genetica. 2005;124:179–86.PubMedView ArticleGoogle Scholar
  17. Dawson MN. Phylogeography in coastal marine animals: a solution from California? J Biogeogr. 2001;28:723–36.View ArticleGoogle Scholar
  18. Wares JP. Community genetics in the northwestern Atlantic intertidal. Mol Ecol. 2002;11:1131–44.PubMedView ArticleGoogle Scholar
  19. Dawson MN. Incipient speciation of Catostylus mosaicus (Scyphozoa, Rhizostomeae, Catostylidae), comparative phylogeography and biogeography in Southeast Australia. J Biogeogr. 2005;32:515–33.View ArticleGoogle Scholar
  20. Pelc RA, Warner RR, Gaines SD. Geographical patterns of genetic structure in marine species with contrasting life histories. J Biogeogr. 2009;36:1881–90.View ArticleGoogle Scholar
  21. 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.PubMedView ArticleGoogle Scholar
  22. Altman S, Robinson JD, Pringle JM, Byers JE, Wares JP. Edges and overlaps in Northwest Atlantic phylogeography. Diversity. 2013;5:263–75.View ArticleGoogle Scholar
  23. 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.View ArticleGoogle Scholar
  24. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  25. Hellberg ME, Burton RS, Neigel JE, Palumbi SR. Genetic assessment of connectivity among marine populations. Bull Mar Sci. 2002;70:273–90.Google Scholar
  26. Irwin DE. Phylogeographic breaks without geographic barriers to gene flow. Evolution. 2002;56:2383–94.PubMedView ArticleGoogle Scholar
  27. 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.View ArticleGoogle Scholar
  28. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  29. 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.View ArticleGoogle Scholar
  30. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  31. 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.View ArticleGoogle Scholar
  32. Nosil P, Vines TH, Perspective FDJ. Reproductive isolation caused by natural selection against immigrants from divergent habitats. Evolution. 2005;59:705–19.PubMedGoogle Scholar
  33. Wu CI. The genic view of the process of speciation. J Evol Biol. 2001;14:851–65.View ArticleGoogle Scholar
  34. Norris RD, Hull PM. The temporal dimension of marine speciation. Evol Ecol. 2012;26:393–415.View ArticleGoogle Scholar
  35. Villamor A, Constantini F, Abbiati M. Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates. PLoS One. 2014;9:e101135.PubMedPubMed CentralView ArticleGoogle Scholar
  36. Nosil P, Harmon LJ, Seehausen O. Ecological explanations for (incomplete) speciation. Trends Ecol Evol. 2009;24:145–56.PubMedView ArticleGoogle Scholar
  37. Martin P, Zuccarello GC. Molecular phylogeny and timing of radiation in Lessonia (Phaeophyceae, Laminariales). Phycol Res. 2012;60:276–87.View ArticleGoogle Scholar
  38. 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.View ArticleGoogle Scholar
  39. 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.View ArticleGoogle Scholar
  40. 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.PubMedView ArticleGoogle Scholar
  41. 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.View ArticleGoogle Scholar
  42. 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.View ArticleGoogle Scholar
  43. 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.View ArticleGoogle Scholar
  44. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  45. 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.PubMedView ArticleGoogle Scholar
  46. 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.View ArticleGoogle Scholar
  47. 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.View ArticleGoogle Scholar
  48. Anderson MJ. Permutation test for univariate or multivariate analysis of variance and regression. Can J Fish Aquat Sci. 2001;58:626–39.View ArticleGoogle Scholar
  49. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  50. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  51. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  52. 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.View ArticleGoogle Scholar
  53. 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.PubMedView ArticleGoogle Scholar
  54. 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.PubMedView ArticleGoogle Scholar
  55. 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.View ArticleGoogle Scholar
  56. Lancellotti DA, Vásquez JA. Biogeographical patterns of benthic invertebrates in the southeastern Pacific littoral. J Biogeogr. 1999;26:1001–6.View ArticleGoogle Scholar
  57. 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.View ArticleGoogle Scholar
  58. Brattström H, Johanssen A. Ecological and regional zoogeography of the marine benthic fauna of Chile. Sarsia. 1983;68:289–339.View ArticleGoogle Scholar
  59. 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.View ArticleGoogle Scholar
  60. 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.View ArticleGoogle Scholar
  61. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  62. Provan J, Bennett KD. Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008;23:564–71.PubMedView ArticleGoogle Scholar
  63. Hellberg ME, Balch DP, Roy K. Climate-driven range expansion and morphological evolution in a marine gastropod. Science. 2001;292:1707–10.PubMedView ArticleGoogle Scholar
  64. Cox LN, Zaslavskaya NI, Marko PB. Phylogeography and trans-Pacific divergence of the rocky shore gastropod Nucella lima. J Biogeogr. 2014;41:615–27.View ArticleGoogle Scholar
  65. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  66. 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.View ArticleGoogle Scholar
  67. GDF W. Functional morphology and evolution of isopod genitalia, Columbia University press. Crustacean Sexual Biology, New York/Oxford. 1991;228245.Google Scholar
  68. 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.Google Scholar
  69. Prada C, Hellberg ME. Long prereproductive selection and divergence by depth in a Caribbean candelabrum coral. Proc Natl Acad Sci. 2013;110:3961–6.PubMedView ArticleGoogle Scholar
  70. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  71. Simmons LW. Sexual selection and genital evolution. Austral Entomol. 2014;53:1–17.View ArticleGoogle Scholar
  72. 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.PubMedView ArticleGoogle Scholar
  73. 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.PubMedGoogle Scholar
  74. 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].
  75. Swofford DL. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). 2002. Sunderland, MA, USA: Sinauer Associates.Google Scholar
  76. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1998;123:585–95.Google Scholar
  77. Fu YX, Li WH. Statistical tests of neutrality of mutations. Genetics. 1993;133:693–709.PubMedPubMed CentralGoogle Scholar
  78. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.PubMedView ArticleGoogle Scholar
  79. Petit RJ, El Mousadik A, Pons O. Identifying populations for conservation on the basis of genetic markers. Conserv Biol. 1998;12:844–55.View ArticleGoogle Scholar
  80. 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.View ArticleGoogle Scholar
  81. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.PubMedPubMed CentralView ArticleGoogle Scholar
  82. Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.PubMedPubMed CentralView ArticleGoogle Scholar
  83. Darriba D, Taboada G, Doallo R, Posada D. jModelTest 2: more models, new 612 heuristics and parallel computing. Nat Methods. 2012;9:772.PubMedPubMed CentralView ArticleGoogle Scholar
  84. Nielsen R, Wakeley J. Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics. 2001;158:885–96.PubMedPubMed CentralGoogle Scholar
  85. 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.PubMedView ArticleGoogle Scholar
  86. 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.View ArticleGoogle Scholar
  87. Malay MC, Paulay G. Peripatric speciation drives diversification and distributional pattern of reef hermit crabs (Decapoda: Diogenidae: Calcinus). Evolution. 2010;64:634–62.PubMedView ArticleGoogle Scholar
  88. 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.View ArticleGoogle Scholar
  89. 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.View ArticleGoogle Scholar
  90. Brookfield J. A simple new method for estimating null allele frequency from heterozygote deficiency. Mol Ecol. 1996;5:453–5.PubMedPubMed CentralView ArticleGoogle Scholar
  91. Rousset F. Genepop 007: a complete reimplementation of the Genepop software for windows and Linux. Mol Ecol Res. 2008;8:103–6.View ArticleGoogle Scholar
  92. Kalinowski ST. HP-RARE 1.0: a computer program for performing rarefaction on measures of allelic richness. Mol Ecol Notes. 2005;5:187–9.View ArticleGoogle Scholar
  93. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  94. Goudet J. FSTAT: a program to estimate and test gene diversities and fixation indices. 2002; Version 2.9.3.2. Available from URL: https://www2.unil.ch/popgen/softwares/fstat.htm.
  95. Prichard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.Google Scholar
  96. 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.PubMedView ArticleGoogle Scholar
  97. 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.PubMedPubMed CentralView ArticleGoogle Scholar
  98. 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.View ArticleGoogle Scholar
  99. 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.PubMedView ArticleGoogle Scholar
  100. 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.Google Scholar
  101. 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.View ArticleGoogle Scholar
  102. Reist J. An empirical evaluation of several univariate methods that adjust for size variation in morphometric data. Can J Zool. 1984;63:1429–39.View ArticleGoogle Scholar
  103. Hammer Ø, Harper DAT, Ryan PD. PAST: paleontological statistics software package for education and data analysis. Palaeontol Electronica. 2001;4:4–9.Google Scholar
  104. 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/ Google Scholar
  105. Anderson MJ, Gorley RN. Clarke KR. PERMANOVA + for PRIMER: guide to software and statistical methods. UK: PRIMER-E, Plymouth; 2008.Google Scholar

Copyright

© The Author(s). 2019

Advertisement