Highly Connected Populations and Temporal Stability in Allelic Frequencies of a Harvested Crab from the Southern Pacific Coast.

For marine invertebrates with a benthic adult form and a planktonic larva phase, the connectivity among populations is mainly based on larval dispersal. While an extended larval phase will promote gene flow, other factors such as an intensive fishery and geographical barriers could lead to changes in genetic variability. In this study, the population genetic structure of the commercial crab Metacarcinus edwardsii was analyzed along 700 km of the Chilean coast. The analysis, based on eight microsatellite loci genotyped from megalopae and adult crabs, considered temporal and spatial patterns of genetic variation. The results showed no evidence of spatial patterns in genetic structure, suggesting high connectivity among the sampling sites. The temporal analysis showed no evidence of changes in allele frequencies and no evidence of a recent bottleneck. The lack of spatial structure and allele variation over time could be explained by the interaction of factors such as i) low reproductive variance due to the capability of females to store sperm in the seminal receptacle, which can be used for successive broods, ii) high larval dispersal and iii) high individual reproductive output. Using our data as priors, a genetic modelling approach coincided, predicting this temporal and spatial stability. The same analysis showed that a reduction in population size leads to the loss of genetic variability in populations, as well as of the genetic cohesiveness between populations, pointing out the importance management for species under exploitation, such as M. edwardsii.

Introduction

Most benthic marine invertebrates are sessile or have little movement in their adult stage; dispersion happens during the planktonic larval period [1,2]. Typically, in this process a level of larval self-recruitment occurs in populations of origin, as well as an exchange of individuals among geographically separate populations [3]. The behavior of the larvae favors horizontal migration along the coast and vertical migration, promoting larval retention in nearshore waters [4,5,6].

Population structure cannot be predicted only by the mode of larval development and the extent of the planktonic period [7]. Species without a planktonic larval phase usually show higher levels of population genetic differentiation than species with planktonic larval development. However, species with planktonic larvae may have small-scale spatial structure, for example the starfish Patiriella regularis, where upwelling acts as a barrier to gene flow [8]. On the other hand, there are species with direct development that have genetic homogeneity over significant spatial scales, such as the gastropods Littorina sitkana and Nucella lapillus, where larval transport is aided by rafting [7, 9]. In this context, larval dispersal involves interactions between larval behavior and ocean transport [10]. Therefore, the genetic structure of marine populations is not only influenced by the mode of development, but also by oceanographic conditions such as tides, coastal upwelling [11,12], winds, physical barriers and life history traits [13].

Species with long lived pelagic larvae do not always show homogenous genetic patterns in time and space. For example, temporal variations in allele frequencies have been described [14,15] and related to a phenomenon called chaotic patchiness [16,17]. Several hypotheses try to explain the factors that can cause genetic heterogeneity, for example strong genetic drift caused by high mortality during early development, the action of natural selection during the larval stage and catastrophic events in the marine environment could change the pattern of larval recruitment in different years [18]. The possibility of low larvae mixture [19] and high reproductive variance [20] may also occur.

For species facing commercial fishing pressure, population size is an important factor influencing genetic variability and population connectivity [21]. Genetic variability is related to population size in relation to the effect of the genetic drift [22]. It would be expected that populations with reduced size will have lower genetic variability and changes in allele frequencies generation by generation. Changes in allele frequencies and allele diversity have been well described in marine species under exploitation [23,24], however other species have shown temporal stability in genetic diversity despite fishing pressure [25,26]. Therefore, it is important to determine the significance of effective population size in genetic stability and population cohesiveness.

Another important point to consider for the genetic structure of populations is the presence of biogeographic barriers and the phylogeographic breaks. Well-known phylogeographic breaks in the Pacific are located at Point Conception (34.5° N) in California [27], 30–32° S in northern [28] and 42° S in southern Chile [29]. While studies performed in Chile have identified species limited by the break at 30–32° S [30,31], other species did not align with these biogeographic breaks [32,31].

While these studies focused on phylogeographic analysis with a historical component, for harvested species the actual patterns of connectivity must be understood. In this context, decapod crustaceans are a good model to study connectivity due to their high fertility and extended larval development. In this study we performed a population genetic analysis on the crab Metacarcinus edwardsii, a crustacean inhabiting from the coast of Guayaquil, Ecuador to the Beagle Canal in Chile [33]. In Chile, M. edwardsii represents around 70% of national crab landings in artisanal fisheries [34]; exploitation is concentrated mainly around Chiloe Island in southern Chile [35]. The pelagic larval development of this species lasts three months (at 14±0.5°C) with five zoea larval stages and one megalopa stage [36] that settle mainly between November and December [37]. Adults inhabit subtidal environments, feeding on carrion and bivalves [38].

While previous studies have described aspects of the species reproductive biology [39,40], larva identification [41] and the effects of the fisheries on the reproductive potential of males [35], this information must be related to population genetics studies in order to be integrated in measures of protection for exploited populations. Our study covered 700 km of coastline, considering different sites and individuals sampled over 4 different years, in order to describe the spatial and temporal genetic structure of M. edwardsii.

Materials and Methods

Sampling, DNA extraction and microsatellite amplification

To study the spatial genetic structure, 291 adult crabs were sampled from six localities: Concepción (CO), Los Molinos (LM), Ancud (AN), Calbuco (CA), Dalcahue (DA) and Quellón (QU) (Fig 1). Samples were obtained by Scuba diving during 2011 to 2014. To study the temporal genetic variability a total of 157 megalopa larvae were collected from Los Molinos (LM) during four years (2009 to 2012) using a passive larva collector, described by Pardo et al. 2010 [42]. Each megalopa was identified using the description of Pardo et al. (2009b) [41] and stored in 95% ethanol until analysis. Sample size and year of sampling are shown in Table 1. All collections and analyses were conducted in Chile and complied with its existing laws (Resolución Exenta No. 3088 Subsecretaría de Pesca).

Fig 1

Metacarcinus edwardsii sampling sites.

Table 1

Sample size by site and year.

Group	Locality	Coordinates	Year	N	Code
Adult	Los Molinos (LM)	39°51′16.7″ S; 73°23′40.3″O	2011	24	LM 2011
			2013	24	LM 2013
	Dalcahue (DA)	42°22′46.3″ S; 73°35′42.5″O	2011	24	DA 2011
			2013	24	DA 2013
	Ancud (AN)	41°50′59.8″ S; 73°51′32.5″O	2012	37	AN 2012
			2013	24	AN 2013
	Calbuco (CA)	41°45′47.1″ S; 73°05′20.1″O	2012	16	CA 2012
			2013	16	CA 2013
	Quellón (QU)	43°08′18.4″ S; 73°36′43.4″O	2012	47	QU 2012
			2013	29	QU 2013
	Concepción (CO)	36°21'6.75"S; 72°50'52.68"O	2014	26	CO 2014
Megalopae	Los Molinos (LM)	39°51′16.7″ S; 73°23′40.3″O	2009	35	Megalopae 2009
			2010	51	Megalopae 2010
			2011	37	Megalopae 2011
			2012	34	Megalopae 2012

Genomic DNA was extracted from adult crabs using the salt method described by Aljanabi and Martinez (1997) [43]. DNA was extracted from megalopae using the QIAamp DNA Mini Kit (QIAGEN) and purified with sodium acetate (3M, pH 5.2). Eight polymorphic microsatellite loci were amplified following the procedure of Rojas-Hernández et al. (2014) [44]. Fragment analysis was performed using an Applied Biosystems 3100 sequencer at the Pontificia Universidad Católica de Chile. The allelic data matrix was built using the Peak Scanner software (Applied Biosystems).

Adult data analysis

The number of alleles per locus, expected (HE) and observed (HO) heterozygosity were estimated with the GENETIX software [45]. With this same software, 5000 permutations on monolocus genotypes were used to test linkage disequilibrium and 5000 allele permutations were used to test for departures from Hardy Weinberg Equilibrium (HWE). FSTAT software [46] was used to estimate the allelic richness, while MICROCHECKER software [47] was used to identify possible null alleles in the microsatellite loci data.

To ensure that the analysis was performed with unrelated individuals, the rxy index [48] was estimated with IDENTIX software [49]. IDENTIX used multilocus data to estimate the relatedness between pairs of individuals within each locality, compared with a distribution of simulated genotypes constructed after an allele permutation among individuals. 1000 permutations were generated to build this simulated distribution.

When highly variable genetic markers are used for the study of population structure, an alternative to the FST index is recommended [50]. One such index is G”ST; which is not influenced by the heterozygosity of samples. Thus, population structure was evaluated using four different methods: i) FST with GENETIX software, ii) G”ST with GenAlEx software [51], iii) a Bayesian approximation with STRUCTURE software [52] and iv) an iterative reassignment of individuals with FLOCK software [53]. The statistical significance of the global FST was evaluated with GENETIX software using 5000 permutations. The GenAlEx software allowed the estimation of the global G”ST and its significance was evaluated with 9999 permutations. The STRUCTURE software infers the most probable number of populations; the analysis was performed using the admixture model and the correlated allele frequencies model. The procedure was run 5 times for each K estimation (from K = 1 to K = 11) with a burn-in of 200000 MCMC iterations. LnP(D) values obtained for each K value were compared using the Structure Harvester software [54]. FLOCK software randomly divides the collection of genotypes into K genetic groups; the genotypes are reassigned to the group with the highest probability of belonging at each iteration following the multilocus method of maximum likelihood described by Paetkau et al. (1995) [55]. For our analysis, 50 runs and 30 iterations were used for values of K = 2 to 6 with a LLOD value ranging from 0 to 0.6.

To assess association between genetic differentiation and geographic distance, a Mantel test was performed using the GENETIX software. The statistical significance was tested with 5000 permutations. Since temporal comparisons within each locality showed no significant difference, data from each locality were pooled for posterior spatial analysis. In order to determine the migration patterns of M. edwardsii among localities, historical migration rates were estimated using the MIGRATE software [56]. MIGRATE uses a coalescent approach to estimate mutation-scaled migration rates (M) for each population over the last 4Ne generations. The analysis used maximum likelihood, the Brownian motion mutation model and the matrix migration model containing 10 short chains of 40000 steps and three long chains of 400000 steps, after a burn-in step of 40000 and a static heating scheme of 6 chains with increasing temperatures (1, 1.25, 1.67, 2.5, 5 and 1e6; the swapping interval was 1). Mutation rates of 1 x 10−2 and 1 x 10−3 were used to estimate the migration rate.

Effective population size

The effective population size was estimated using two different methods. First, the Bayesian method was implemented in ONeSAMP software [57] with the default settings. The second was the LD method implemented in NeEstimator [58]. In addition, a Kruskal-Wallis test was performed in order to compare Nb and Ne values obtained from both software.

Searching for evidence of a recent bottleneck

BOTTLENECK software [59] was used to detect evidence of a recent genetic bottleneck. Briefly, BOTTLENECK estimates the likelihood of recent reductions in effective population size by comparing the expected heterozygosity under HWE with the expected heterozygosity under mutation-drift equilibrium. The two phase mutation model with a 70% stepwise mutation was used.

Megalopa data analysis

Four cohorts of megalopae were sampled in Los Molinos from 2009 to 2012 (Table 1). The number of alleles per locus, linkage disequilibrium, HE, HO and departures from HWE were estimated for each year with the GENETIX software. To determine differences among the cohorts analyzed, FST and G''ST were estimated in the same way as for adults. Two analyses were performed: the first comparing the four megalopa cohorts and the second including each megalopa cohort and the adults collected in the same locality. STRUCTURE and FLOCK software were also used to detect differences among megalopa cohorts with same settings as in the adult analysis.

Taking into account that larval retention near the coast and the reproductive variance could generate non-random larval settlement in the same area, the likelihood of detecting related individuals in a sample was estimated using the rxy value and the statistical significance was obtained with 1000 permutations in IDENTIX software.

For iteroparous species, estimates of effective population size need to consider samples composed of individuals of different ages. Therefore, the estimated population size of each cohort is interpreted as an approximation of the actual number of breeders per season (Nb), instead of an effective population size per generation (Ne) [60]. Considering this difference, Nb was estimated for each megalopa cohort and an approximate Ne was estimated using pooled data from these four cohorts. These calculations were done with ONeSAMP and NeEstimator software.

Simulation of the effect of population size and the migration rate

For economically important species, exploring the effects of change in both population size and migration rate is relevant in order to understand population connectivity and changes in genetic variability. Considering that no spatial structure was observed in M. edwardsii and there was no evidence for a recent bottleneck, we explored the factors that would affect this pattern in a scenario of increased fishing. Currently, around five thousand tons of this species are harvested every year from southern Chile, affecting the population size [35]. Migration rates could be influenced by changes in both oceanographic circulation and upwelling processes [61]. Thus, this simulation tested for an effect of a reduction in population sizes and migration rates among sampling sites. The simulated datasets were built using EASYPOP 2.0.1 software [62] considering the following conditions: diploid organism, two sexes, random mating system, 6 sampling sites, equal number females and males in each generation, island migration model, eight loci with free recombination, a 0.01 microsatellite mutation rate, a mixed mutation model with Kam = 0.3, allelic states = 24 (similar to the mean allele number observed in this study) and maximum variability of the initial population. Three variables were tested simultaneously in the simulation: i) population size: 100000, 50000, 10000, 5000, 1000 and 500 individuals, ii) migration rate: 10% (as observed in this study), 1%, 0.1%, 0.01% and 0.001%, and iii) three different time horizons: after 10, 50 and 100 generations.

Results

Analysis of adults

All loci were polymorphic; the mean allele richness varied from 6.29 (Cedw15) to 25.44 (Cedw4) with a mean of 24 alleles per locus (S1 Table). Departures from HWE were detected in three out of eight microsatellite loci (α = 0.01), but deviations were not associated with any specific locus, sampling period or location. Null alleles do not appear to be responsible for the observed departures from HWE. Ten out of the 168 comparisons showed evidence for linkage disequilibrium (α = 0.01), however this disequilibrium was not observed in the same pair of alleles at all sites. The mean rxy values estimated for the groups of individuals collected in the same year showed no statistical evidence (Permutation, P > 0.05) for highly related individuals. This evidence suggests that samples are composed of unrelated individuals (S1 Table).

For the temporal genetic structure, FST and G”ST estimated for samples obtained in the same locality but in different years showed no significant differences: LM (FST = 0.001, P = 0.366 and G”ST = 0.026, P = 0.325), DA (FST = 0.005, P = 0.078 and G”ST = 0.067, P = 0.067), AN (FST = -0.003, P = 0.933 and G”ST = -0.035, P = 0.924), CA (FST = -0.008, P = 0.949 and G”ST = -0.082, P = 0.942) and QU (FST = 0.005, P = 0.068 and G”ST = 0.002, P = 0.523). Considering these results, the next analyses were performed with pooled data from each locality.

For the spatial genetic structure, the global FST and G”ST showed no significant differences among sampling sites (FST = 0.001, P = 0.147 y G”ST = 0.006, P = 0.293). The Structure Harvester showed K = 1 had the highest likelihood L(K) mean value and the lowest standard deviation (Fig 2A). FLOCK software showed the same pattern; the individual assignments in two or more reference groups showed no statistical significance for all runs performed with different LLOD values. It is important to note that FLOCK does not include K = 1 in the analysis. Duchesne and Turgeon (2012) [52] indicate that an indecisive result must be interpreted as an absence of population structure or a low power of the analysis due to the low number of loci or alleles. Considering the high number of alleles found and reasonable sample sizes used in this study, the analysis suggests no spatial or temporal structure in M. edwardsii. While, with regards to the effect of distance, the Mantel test showed no evidence of an association between genetic distance and geographic distance (Mantel test, Z = 33.86, P = 0.1496).

Fig 2

Inference of the number of the genetic clusters estimated through the L(K).

(a.) adults collected from six localities. (b.) megalopae collected at Los Molinos. Each K represents the mean and standard deviation of five independent runs.

Effective population size and recent bottleneck

The Ne estimated were similar between sites, ranging from 16 (LM 2013) and 475 (DA 2011), except for QU 2013 that presented an extreme value of 1574 (Table 2). A Kruskal-Wallis test comparing Nb with Ne values from two software programs showed no significant differences (p = 0.338). Evidence for a recent bottleneck was found in one (DA 2011, P = 0.006) out of 11 possible estimations (CO 2014, P = 0.320, LM 2011, P = 0.371, LM 2013, P = 0.156, DA 2013, P = 0.273, AN 2012, P = 0.098, AN 2013, P = 0.473, QU 2012, P = 0.098, QU 2013, P = 0.156, CA 2012, P = 0.098, CA 2013, P = 0.125).

Table 2

Nb and Ne estimated with ONeSAMP and NeEstimator for adults and megalopae collected from the different sampling sites and years.

			OneSamp		NeEstimator
		n	Ne	Confidence interval	Ne	Confidence interval
Adult	DA 2011	24	475.17	264.83–1185.00	196.7	112.2–708.5
	DA 2013	19	205.45	109.81–512.36	49.4	35.4–79.2
	LM 2011	18	132.37	76.73–280.68	78.2	56.4–124.2
	LM 2013	23	16.68	12.88–22.55	44.8	33.5–65.9
	AN 2012	35	38.23	30.30–49.25	323	190.5–987.3
	AN 2013	22	20.43	15.64–27.94	155.2	94.7–400.8
	QU 2012	32	196.77	133.64–334.86	359.9	205.5–1312.8
	QU 2013	25	1574.59	643.27–5932.69	177.9	110.7–427.4
	CA 2012	16	232.31	122.54–672.29	25.1	20.4–32.1
	CA 2013	14	190.96	92.85–531.93	49.2	33.7–87.2
	CO 2014	24	148.35	89.59–289.64	40.1	34–48.4
Megalopa larvae	Megalopae 2009	25	61.28	42.39–92.39	∞	226.8 - ∞
	Megalopae 2010	38	153.30	95.09–256.62	571.2	194.7 - ∞
	Megalopae 2011	23	315.32	155.67–909.01	1345	122.1 - ∞
	Megalopae 2012	29	10.85	9.05–13.04	118.2	65.5–488.3
	All (≈Ne)	111	198.58	167.62–239.01	∞	∞

Migration between pairs of localities

The mutation-scaled historical migration rates (M) ranged from 1.1 to 3.2 (Table 3) and the estimates of Θ ranged from 0.2 to 7.4 (LM = 6.1255, DA = 5.4971, AN = 6.0260, CA = 4.2633, QU = 4.4672, CO = 6.5019). There was no tendency related to a reduction in migration rate for geographically distant pairs of localities. Considering the mutation rate of 1 x 10−2, the migration rate between localities ranged from 0.01 to 0.031, with a total immigration rate of 0.1 for each location. While, when the mutation rate was 1 x 10−3, the migration rate between localities ranged from 0.001 to 0.003, with a total immigration rate of 0.01 for each location. Considering this result, the mean for self-recruitment (estimated as 1- total immigration rate) ranged from 0.90 (with a mutation rate of 1 x 10−2) to 0.99 (with a mutation rate of 1 x 10−3).

Table 3

Number of migrants per generation (with 95% CI) of Metacarcinus edwardsii by sampling site.

From/to	LM	DA	AN	CA	QU	CO
LM		2.964 (3.175–3.590)	1.350 (1.490–1.638)	2.541 (2.736–3.083)	1.512 (1.659–1.815)	1.528 (1.671–1.823)
DA	3.246 (3.040–3.461)		2.549 (2.740–2.941)	1.091 (1.180–1.316)	2.457 (2.645–3.022)	1.434 (1.573–1.721)
AN	1.066 (0.882–1.315)	3.036 (3.249–3.472)		1.867 (2.034–2.381)	1.421 (1.564–1.716)	2.016 (2.180–2.353)
CA	2.920 (2.723–3.125)	1.276 (1.415–1.564)	1.970 (2.138–2.315)		2.279 (2.460–2.650)	1.373 (1.508–1.738)
QU	1.339 (1.208–1.479)	2.753 (2.956–3.169)	1.242 (1.377–1.630)	2.657 (2.8567–3.215)		1.666 (1.815–1.973)
CO	1.8186 (1.665–1.981)	1.200 (1.335–1.635)	2.412 (2.598–2.794)	1.704 (1.864–2.243)	1.580 (1.822–2.065)

Analysis of megalopa larvae

Seven out of eight microsatellites were used for the megalopa analysis. The locus Cedw4 showed low quality and was excluded for all analyses performed for the megalopa. The seven loci used for megalopae showed high polymorphism with values ranging from 4.92 (Cedw15) to 12.57 (Cedcrab3) (S2 Table). While 1 out of 82 comparisons between pairs of loci showed evidence for linkage disequilibrium (α = 0.01), this finding was not consistent across sampling years.

In the temporal analysis performed for the four years of samples obtained at Los Molinos, the global analysis showed no evidence for differences in allele frequencies (FST = 0.0017, P = 0.1536 and G”ST = 0.016, P = 0.180). Even when the analysis was performed with the adults collected at the same location, no temporal structure was detected (FST = 0.0016, P = 0.1218 and G”ST = 0.011, P = 0.241). This result is consistent with that obtained with STRUCTURE (Fig 2B) and FLOCK software; both suggested K = 1 as the most probable number of genetic groups. Furthermore, there was no evidence of related individuals for any of the four cohorts analyzed (S2 Table), suggesting that zoea larvae belonging to the same litter are separated early in water column.

The effective number of breeders (Nb) was calculated for each cohort (megalopae 2009 to 2012) and the effective population size (Ne) was estimated pooling all megalopa data. Nb and Ne had a similar range of variation and even some values of Nb were higher than values of Ne estimated for the adults sampled in the same locality (Los Molinos) (Table 2).

Modelling of changes in population size and reciprocal migration

The simulation started with 6 sites containing 100,000 individuals and an average of 24 alleles per locus. The estimated Ne in generation t = 1 for each of these six sampling sites had an average Ne = 476, similar to the Ne values obtained with our data.

This simulation allows for the observation of different scenarios at different times (see Figs 3 and 4). First, sites containing 100,000 individuals for 100 generations showed no significant differences in the FST values (Fig 3) and did not reveal a decreasing number of alleles even when the rate of migration was reduced from 10% (observed value for a mutation rate of 1 x 10−2) to a rate 1000 times lower. Second, sites with small population sizes showed differences in allele frequencies and lost part of their allelic diversity. After the 10th generation, the sites containing 1,000 individuals and 10% migrants had significantly different FST values. The same was observed for populations composed of 5,000 individuals and a 0.1% migration rate. Finally, in the worst simulated case (500 individuals and migration rate of 0.001%), a significant value of FST (FST = 0.0128, P < 0.05) was observed. This difference was primarily due to changes in allele frequencies, since only one allele per locus was lost during this time. The decrease in the number of alleles appears to be gradual; at the 50th generation the simulation predicted a loss of 4 alleles per locus (16% of the total alleles) and 7 alleles (29%) at the 100th generation (Fig 4).

Fig 3

Changes in F

Fig 4

Reduction in the number of alleles per locus estimated for the simulation of six populations as a function of the population size and migration rate, data after a) 100, b) 50 and c) 10 generations.

Discussion

The genetic analysis performed in Metacarcinus edwardsii along 700 km of coastline showed a population genetic structure in accordance with the duration of planktonic larvae. Genetic analyses suggest there is no genetic structure at temporal and spatial levels, with a historical mean of 10% of the progeny migrating every year among sampling sites.

Temporal and spatial genetic patterns in M. edwardsii

The traditional FST and the unbiased G”ST showed no evidence for spatial and temporal structure in M. edwardsii; these results were consistent with the analyses performed in both STRUCTURE and FLOCK software. Furthermore, the analysis showed no evidence of isolation by distance. The spatial genetic homogeneity described here for M. edwardsii is consistent with the long planktonic larval period. Other species with long planktonic development have shown a lack of spatial structure, for example invertebrates in the coast of South Africa [63] and the crustaceans Chionoecetes opilio [64] and Cancer pagurus [26]. This lack of structure is related to how far larvae can disperse in the early phases of development. Shanks et al. (2003) [65] reported the dispersal distances of different crustaceans, in particular, Carcinus maenas, a brachyuran with a larval period of 80 days, can disperse from 63 to 173km depending on the zone studied. Similarly, Kinlan and Gaines (2003) [66] reported that invertebrates with feeding larvae show an average dispersal distance of 100 km.

Due to the long planktonic period, high mortality during this stage of development could be expected [67], which along with the variation of oceanographic conditions might have caused changes in allele frequencies between cohorts. Genetic homogeneity among consecutive cohorts observed in M. edwardsii is an unexpected result considering that Los Molinos shows seasonal coastal upwelling events [32], and even had a tsunami which led to changes in coastal topography in February, 2010 [68]. While this homogeneity is an unusual pattern among marine invertebrates, the crab Carcinus maenas [69] and the bivalve Panopea generosa [70] have also shown temporal genetic stability in a similar context of complexity.

Factors explaining the pattern observed in M. edwardsii

The biology of M. edwardsii could explain the observed pattern of genetic structure.

1. Reproductive variance. Genetic homogeneity over time has been associated with species that have low reproductive variance among individuals that compose each population [69]. As was pointed out by Broquet et al. (2013) [71], high reproductive variance is the main factor explaining ephemeral genetic differentiation at small scale (chaotic genetic patchiness) observed in benthic invertebrates. Thus, this low reproductive variance can be explained by characteristics exhibited in the brachyuran crab reproductive system, with capacity to store spermatophores for several months [40]. Therefore, if females do not find mates in one season, sperm transferred during the previous season and stored in seminal receptacles could be used for fertilization. In this condition, all females are able to produce viable broods every year.

2. High larval dispersal. Species with high potential for dispersion given the duration of their pelagic larvae, usually have high levels of connectivity [72]. According to Haye et al. (2014) [31], the time spent by larvae in the water column, as a proxy of the dispersion potential, is the factor that best explains the genetic structure present in marine invertebrates in the transition area (30–33° S) in Chile. M. edwardsii, with its long larval period, is a good example of a species showing a lack of genetic structure at a wide spatial scale.

3. Complete larval mixing during the planktonic period. It has been widely assumed that during the planktonic period, larvae from different locations and/or parents are completely free to mix (e.g. [73]). However, there is evidence showing that discrete patches of larvae could be maintained in the water column on a scale of days to weeks [74,75], allowing the settlement of related individuals in the same place [76,19,77]. The factors that promote or reduce larval mix are related to oceanographic conditions and the behavior of the females and/or larvae, and not to the potential these species have to disperse [19, 20]. In this context, incomplete larval mix has been pointed to as a secondary factor affecting stochastic changes in allele frequencies [71]. The similar allele frequencies observed each year in M. edwardsii and the lack of evidence of relatedness within each megalopa cohort suggested that larvae mix completely during the planktonic period before benthic settlement. Another study performed in the brachyuran Carcinus maenas also showed genetic stability among cohorts [69].

4. Reproductive output. Species with greater reproductive output show more broad-scale connectivity among populations than species with low reproductive output [78]. Brachyuran species have the highest reproductive output among marine invertebrates. Hines (1991) [79] estimated 18,200 to 2,208,000 eggs per brood in nine Cancer crab species on the east coast of the USA. Interestingly, females of the brachyuran Cancer pagurus produced between 605,600 to 2,310,000 eggs per brood [79], and had no population structure in the Norwegian Sea [26]. Our measurements, performed on M. edwardsii, showed values of 500,000 to 1,000,000 eggs per brood.

Effective population size

Traditionally, estimates of Ne were based on samples of middle age; while samples of one cohort provide information about the Nb (model based on semelparous species) [80]. In this case, Ne per generation could be larger than Nb by a factor similar to the generation length (Ne≈G×Nb) [81]. For iteroparous species, Waples et al. (2014) [82] suggested that samples from each cohort provide information to estimate the Nb and also could be used to the estimate Ne; however, this requires taking into account a correction that considers some life history traits. In our study, the Nb and Ne had similar values, confirming the analysis of Waples et al. (2014) [82]. Similar results also have been described for the frog Bufo calamita [83], the sturgeon Acipenser fulvescens [84] and the salamanders Ambystoma opacum and A. talpoideum [85].

Genetic basis for fisheries management and resource conservation

M. edwardsii is the most important artisanal crab fishery in Chile; which makes it important to understand the current state of the species in order to understand the effect of the fishery on its future. This study was conducted in the most intense crab fishing area in the country [35] (Concepción, Los Molinos and Calbuco are low fishery intensity localities, while Ancud, Dalcahue and Quellón are high fishery intensity localities), so the conclusions made here are valid for future management plans. Our results indicate that the larvae produced in unexploited areas are naturally connected with exploited areas, aiding to avoid the reduction of stocks and irreversible loss of genetic diversity. Considering that crab fishery increases every year due to the pressure imposed by the international market and reallocation of fishery effort after the depletion of other resources (i.e. finfishes, Chilean abalone), crab landings in the southern Chile have increased in the last decade and stock seems be resilient to increment of effort. High larval subsidies of individuals from unexploited areas could be helping maintain stock. Considering our data and the simulation, the crab fishery should remain an artisanal-only activity [86] and control should be reinforced in effective management measures in order to avoid drastic changes in genetic fluxes.

Summary of genetic values by sampling site.

(DOCX)

Click here for additional data file.

Summary of genetic values for each megalopa cohort.

(DOCX)

Click here for additional data file.