Linking Acrosomal Morphotypes and Genetic Divergence in an Inter-Oceanic Mussel from Pacific and Atlantic Coasts: A Case of Incipient Speciation?

1. Introduction

Within an evolutionary context, both the taxonomic classification of organisms and the understanding of the diversity of biological life have historically been based on the description of morphological traits [1]. In fact, morphology is a very useful way of understanding evolutionary processes and has traditionally been used in systematics and taxonomy. In bivalve molluscs such as clams [2], oysters [3,4], and mussels [5,6,7,8], sperm morphological traits have been used in systematic studies because its ultrastructure appears to be highly conserved at species level, being used to differentiate closely related species [2]. Although the role of sperm morphology in fertilization is poorly understood it is clear that spermatozoa vary widely in size, shape and ornamentation between closely related species [9]. The interest in spermatozoa principally comes from their unusual life history and unparalleled diversity, which appear to be inextricably linked. Unique among metazoan cells, spermatozoa are cast forth from the soma into a foreign environment to spend their (pre-fertilization) lives as free-living organisms [10]. Evidence has demonstrated that sperm exhibit a remarkable diversity in form, mirroring the range of fertilization environments in nature [10].

In this work, we link acrosomal morphotypes with genetic divergence at the geographic macroscale using as a model species the inter-oceanic mussel Perumytilus purpuratus. This species is broadly distributed in the Southeast Pacific coast, from Guayaquil (3°S) in Ecuador to Cabo de Hornos (53°S) in the continental Chilean coast [11], and goes up through the Patagonia (41°S) in the south Atlantic coast of Argentina [12,13]. Such as in other external fertilizing species of Mytilidae, the sperm of P. purpuratus are released to the water, and morphologically are of the primitive or ect-aquasperm type characterized by a conical or oval nucleus [5,6,14,15]; also, have acrosomes of variable complexity, an intermediate piece composed of 4 or 5 mitochondria surrounding a pair of centrioles, and a simple flagellum [5,6,14]. Between the spermatic organelles, the acrosome appears to be the most interspecifically variable organelle, displaying diversity in shape and size [5,7,8,16,17], being the principal focus of many descriptive studies in this area.

Although sperm morphology and its ultrastructure might provide useful traits for phylogenetic analyses [2,18,19], it is important to highlight that in the marine field few studies have been conducted to explore intraspecific variation of quantitative sperm traits and its association with phylogeny. An exception was the comparative study made by [6] on P. purpuratus populations, where, based on acrosome size and sequence analyses of mitochondrial rRNA gene (16S), two genetically divergent sperm morphotypes strongly related with their geographic origin were observed: a southern morphotype (previously described by [5] for Valdivia, 41°S), and a northern novel acrosomal morphotype found at Valparaíso (32°S) and Antofagasta (23°S) [6]. The genetic divergence observed by these authors agreed with later studies at the geographical macroscale that confirmed a genetic discontinuity inside P. purpuratus distribution. On the one hand, based on mitochondrial (COI) and nuclear (18S and 28S) markers, [20] observed two divergent lineages, designed as north and south clades. Using samples from Pacific (Chile) and Atlantic (Argentina) coasts, these authors found that the samples from the Chilean coast between 20°S and 36°S corresponded to the north clade, while those from localities south of 39°S on the Pacific (Chile) up to 42°S on the Atlantic (Argentina) corresponded to the south clade. Alternatively, based on five polymorphic microsatellites markers, [21] determined two genetically divergent groups for P. purpuratus from Pacific coast of Chile, which were designated by them as stocks: a southern stock south of 40°S and a north-central stock between 18°S and 40°S. These authors suggested this sharp structure as historical in nature and maintained in P. purpuratus due to its environment-related restricted larval dispersal, to its limited potential for dispersal, or to a genetic/reproductive barrier between its historical lineages (e.g., see [6,21]).

Overall, these findings suggest that, a strong genetic break between northern and southern regions of Chile could be the reflection of the postglacial colonization after the Last Glacial Maximum (LGM) [6,20,21]; and that both, morphological and genetical divergence inside Perumytilus distributional range (i.e., between north and south populations) could be, in fact, scientific evidence for two cryptic species [6,20,21]

Our aim was to deepen the study in sperm morphology and their geographic variation in P. purpuratus and its association with phylogeny and population genetics to contrast whether there is evidence of clinal variation in sperm traits of P. purpuratus or whether specimens from the northern and southern regions could correspond to cryptic species. In the last case, we expected that both sperm morphology and genetic patterns of discontinuity were strongly correlated.

To address these questions, we analysed samples of P. purpuratus covering a coastal gradient of ca 7,180 km of coastline from Pacific to Atlantic Ocean. The sperm morphology was explored by scanning electron microscopy (SEM) and statistically significant differences among localities were explored by lineal models approaches. The genetic analyses were conducted using partial mitochondrial and nuclear genes.

2. Materials and Methods

2.1. Sample collection

Freshly adults of Perumytilus purpuratus, with a maximum shell length > 15 mm, were collected from mid-intertidal rocks along Pacific and Atlantic coasts. The whole data used in this study was obtained from a total of 21 localities covering a coastal gradient of ca. 7,180 km of coastline (Table 1, Figure 1), which were sampled between 2010 and 2015.

2.2. Sperm morphology analyses

The morphological analyses were performed using quantitative sperm traits dataset of a total of sixteen localities (Table 1, Figure 1). The individuals were opened, and the sex was identified by colour and thickness of the mantle tissue, due that in mature male individuals it gets thicker and turns into a cream colour. Small pieces of the testes were stored in labelled microtubes, fixed in 2.5% glutaraldehyde in 0.2 M phosphate buffer (pH 7.4) and preserved to 4°C until the scanning electron microscopy analyses.

2.3. Scanning electron microscopy (SEM)

For the samples obtained between 2010 and 2013, 100 µL of sperm suspension from five selected individuals were mixed, and a drop of sperm suspension was placed on a cover glass, prefixed with 2.5 % glutaraldehyde in a 0.2 M phosphate buffer (pH 7.4) for 2 h at 4°C, and then post-fixed with 1 % osmium tetroxide in a phosphate buffer for 2 h at 4°C. The material was dehydrated in a graded ethanol series, critical point dried, coated with gold, and then observed and recorded with a Leo-420 scanning electron microscope.

For the sperm sampled in 2015, three individuals selected by locality were fixed in 4% glutaraldehyde in a 0.2M phosphate buffer (pH 7.4) 400 µL volume and stored at concentration of 6x10^7 sperm/mL until analysis. Prior to analysis, the sperm samples from each locality were mixed into a tube, then were washed twice with distilled water by 5 minutes, centrifuged at 1000 rpm and filtered using a Nylon Membrane Millipore (40µm). Subsequently, a droplet of sperm suspension was applied onto an aluminum double-faced adhesive tape affixed to the SEM stub. Visualization was performed using a secondary detector (SE) in Scanning Electron Microscope under the following parameters: 2.5kV, working distance 5.7mm and magnification 6,000x in Scanning SU 3500 HITACHI-Japan.

2.4. Statistical analyses

From scan electron micrographs, length of acrosome and head of thirty spermatozoa by locality were measured using the Sigma Scan Pro v.5 and ImageJ programs following to [6]. Then, acrosome/head ratio (%) was calculated, and basic descriptive statistics such as mean, standard deviation and maximum and minimum values by locality were estimated for all sperm parameters (acrosome, head and the acrosome/head ratio (%)).

To determinate differences of mean lengths of acrosome and head among localities, one-way ANOVA with Locality as fixed factor, followed by a post hoc Tukey’s test adjusted for multiple comparisons were performed using PROC MIXED model (SAS). If ANOVA assumptions were not met for raw data, the rank transformed data were used [22].

For the relationship acrosome-head, significant differences among localities were determined using a more robust statistical approximation, instead acrosome/head ratio (%) [23]. In this way, we performed an ANCOVA with Acrosome as dependent variable, and Head as covariable (to control for size effects), using Locality as fixed factor. Then, a Tukey-Kramer post hoc test for multiple comparisons was applied using the Adjusted Acrosome length. With this test we explored the homogeneity of slopes assumption, where no interaction between Locality and Head must be observed, and therefore no evidence that the linear relationships between Acrosome and Head differed between Localities. The ANCOVA was run independently for northern and southern localities (i.e., northern, and southern lineages determined by the molecular analyses). We used PROC MIXED (SAS) given that the dependent variable (acrosome) had a normal distribution.

For data exploration, statistical analysis and graphing, Minitab ® v 19.2020 (Minitab; www.minitab.com), SAS v 9.4/SAS Enterprise Guide 8.3 (Institute Inc., NC, USA), and OriginPro v 2023b SR1 (OriginLab Corporation, Northampton, MA, USA) software were used.

2.5. Molecular analyses

To explore the morpho-genetic relationships, molecular analyses were exclusively conducted on male specimens of P. purpuratus, given the previously reported double uniparental inheritance (DUI) of the mitochondrial genome in this species [24,25].

The genetic analyses were conducted at geographic macroscale (ca. 7,180 km from Pacific to Atlantic coast) to test hypotheses of genetic discontinuity [6,21,26] in relationship with the geographically differentiated sperm (acrosomal) morphotypes of P. purpuratus (e.g., north vs south morphotypes; previously described by [6]).

2.6. DNA extraction, amplification, and sequencing

Total genomic DNA of 1-10 males by location were used for molecular analyses (Table 1). DNA was extracted from small amounts of ethanol-preserved adductor muscle tissue using the E.N.Z.A. TM Tissue DNA Kit (Omega Biotek, Inc.).

For the study of intraspecific genealogies of P. purpuratus at macro-scale level, partial mitochondrial 16S rRNA and nuclear 28S rRNA genes were amplified by polymerase chain reaction (PCR) using the primers 16S-RA/16S-RB [27] and D23F/D6R [28] respectively. The PCR program included a denaturing cycle of 95°C for 5 min, followed by 35 cycles of 95°C for 45 s, an annealing temperature of 50°C for 16S marker and 52°C for 28S marker for 45 s and 60 s respectively, 72°C for 60 s and a final extension of 10 min at 72°C. The PCR reaction had a total volume of 25 µL, containing between 10 and 20 ng of purified DNA, 2.5 µL of 5X NH4 Reaction Buffer (160 mM (NH₄)₂SO₄, 670 mM Tris–HCl, pH 8.8) (1X of Green GoTaq® Flexi Buffer), 10 mM of pre-mixed dNTPs, 2.5 mM of MgCl₂, and 1.25 U of GoTaq Flexi DNA polymerase (Promega).

DNA extractions and PCR products were visualized by UV transilluminator in 2% agarose gels. The purified amplifications were sequenced by Macrogen Inc. (Seoul, South Korea) and STABVida (Lisboa, Portugal) with the reverse primers used in PCR amplifications. Sequencing was performed in capillary automatic sequencers (ABI 3730XL or ABI PRISM 3130, Applied Biosystems®).

The sequences were edited using ProSeq v 3.5 [29] and MSA was done using default parameters of the MAFFT online server platform [30], which use the multiple alignment based on the fast Fourier Transform (FFT) [31]. All sequences were deposited in GenBank (accession numbers: OR196854 - OR196919 for 16S and OR197218 - OR197341 for 28S sequences).

2.7. Bayesian reconstructions and genetic structure

To infer some phylogenetic signal and to visualize the relationship between genetic divergence and sperm morphotypes inside P. purpuratus distribution, a Bayesian Inference (BI) using the MCMC method was performed for each genetic marker in MrBayes v 3.2.7 [32]. In both, mitochondrial 16S rRNA and nuclear 28S rRNA datasets, two closely related Mytilidae species were used as outgroups: one Brachidontes rodriguezii male individual sampled for us in Las Grutas, Neuquen province of Argentina (40°48’S/65°05’W), and one sequence downloaded of NCBI database of Ischadium recurvum. This last species was also used as outgroup in studies of Brachidontes complex including P. purpuratus populations [33,34,35]. For each gene dataset, the best-fit substitution model determinate by BIC in jModelTest v 2.1.8 [36] was HKY. Each BI analysis included two independent runs, four chains and a burn-in of 25% for a total of 10,000,000 generations sampled every 1000 generations. The resulting consensus tree was then visualised and edited in FigTree v 1.4.4 [37].

To explore the most likely number of gene pools (k) and their spatial boundaries, the population genetic structure of P. purpuratus was explored using the spatial Bayesian clustering models implemented in GENELAND v 4.9.2 R package [38]. For each gene dataset, a preliminary Bayesian exploration of the number of subpopulations (k) was performed under four genetic landscape scenarios: spatial/non-spatial combined with correlated and uncorrelated frequencies. The putative number of gene pools (k) was inferred from the modal value with the highest likelihood, and the spatial boundaries were detected as geographical areas with the lowest posterior probability of membership to clusters. For each landscape scenario, ten independents MCMC chains for k = 1-10 were assayed to estimate the most likely value of k, using as settings 1,000,000 iterations, and thinning of 100. Finally, an independent run was performed with the most likely k identified in the preliminary exploration using 5,000,000 iterations and thinning of 1000.

To test the rationality of defined 16S/28S clades and k-groups from BI and GENELAND Bayesian clustering analyses, respectively, pairwise F_ST values were computed and hierarchical analyses of molecular variance (AMOVA) of defined groups were executed in Arlequin v 3.5 [39] using non-parametric permutation procedures with 10,000 iterations.

2.8. Genetic diversity and demographic history estimation

Standard genetic diversity indices such as the number polymorphic sites (S), number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (Pi) and average number of nucleotide differences (k) were calculated using DnaSp v 6 [40] by each adjusted acrosomal morphotype and lineages previously identified (i.e., clades and k-groups from Bayesian analyses).

To reveal the relationships among haplotypes, a median-joining network [41] was constructed for each genetic dataset (16S and 28S) using the default parameters implemented in PopART software [42].

Finally, the demographic history of P. purpuratus was inferred for the genetic clades of the 16S dataset, by means of the neutrality tests of Tajima’s D, Fu’s Fs and Ramos-Onsins and Rozas’s R2 which were calculated using DnaSP v 6 [40] under the model Population Growth-Decline. Another method, pairwise mismatch distributions were carried out to infer whether demographic expansions had occurred using DnaSP v 6 [40].

3. Results

3.1. Sperm morphology analyses

Along the geographical distribution range analysed in this study, sperm of P. purpuratus showed the primitive ect-aquasperm morphology [5,6] characteristic of free-spawning bivalves [14,15].

Acrosomes and heads from a total of 480 spermatozoa were measured and statistically analysed. Mean acrosome and head length varied geographically, showing the smallest sizes at northern localities up to Lota (37°S), and thereafter increasing steeply towards the southern localities (Table 2; Figure 2a,b). In fact, there were statistically significant differences among localities (ANOVA on their rank-transformed lengths: F_{15, 464} = 112.19, p < 0.0001, for acrosomes; and F_{15, 464} = 131.37, p < 0.0001, for heads).

Along our study area, a total of three sperm morphotypes were identified (Figure 3). These morphotypes varied significantly in relationship with their maximum acrosome, head length and acrosome/head ratio (%), and therefore were named by us as short, intermediate, and long morphotypes (acrosome/head ratio = 27.8%, 49.2% and 56%, respectively; Table 2).

Short morphotype was observed in the northern localities up to Lota, where the shortest acrosomal and head lengths, typical of the smaller sperm morphotype defined by [6] (Table 2; Figure 2a,b), were significantly different from all other southern localities (adjusted Tukey test: p < 0.0001; except for mean head length from Lebu that did not differ significantly of Antofagasta (p = 0.115; Tables S1 and S2). But did not differ significantly among them for acrosome (adjusted Tukey test: p > 0.298; Table 3), and for head length (adjusted Tukey test: p > 0.292; except for Antofagasta that was significantly different to Valparaíso (p = 0.009; Table S2).

For southern localities, the intermediate and long morphotypes were identified. Punta Arenas and Chiloé showed the longest mean acrosome and head lengths (Table 2, Figure 2a,b), which corresponded to the longest morphotype found by us. This morphotype was significantly different to all the other localities (adjusted Tukey test: p < 0.0001; Tables S1 and S2) and there were not significant different between them (adjusted Tukey test: p = 1.0; Tables S1 and S2). The other southern localities showed intermediate sizes for acrosome and head lengths (Table 3, Figure 2a,b).

For the ANCOVA analyses an outlier was removed to comply normality assumption, but the results were the same as the analyses with the outlier. The terms of interaction between covariable and locality (head * locality) were not significantly different for both northern localities (ANCOVA, F_{3, 112} = 1.55, p = 0.205) and southern localities (ANCOVA, F_{11, 335} = 0.68, p = 0.755). Therefore, the slopes were homogenous within each region: northern localities: acrosome = 0.843 + 0.064 head (p = 0.075; CL95% = -0.007 - 0.135) and southern localities: acrosome = 0.273 + 0.447 head (p < 0.0001; CL95% = 0374 - 0.520). And the slopes between regions were significantly different as their 95% confidence limits not overlapped.

The adjusted acrosome length for northern and southern localities did not vary significantly with latitude (Lineal Regression, adjusted acrosome length = 1.012 + 0.003 latitude, p = 0.65, northern; adjusted acrosome length = 1.440 + 0.027 latitude, p = 0.055, southern; Figure 2c).

For the northern localities, the adjusted acrosome did not differ significantly between Antofagasta, Valparaíso and Tumbes (Tukey test: p > 0.786; Table 4, Table S3; Figure 2d), but Lota was longer and significantly different (Tukey test: p < 0.0001; Table 4, Table S3; Figure 2d). All the four northern localities showed a mean adjusted acrosome = 1.098 µm (CL = 1.090 - 1.106) which corresponded to the short morphotype.

For the southern localities, the adjusted acrosome length was segregated into two different groups which differed statistically (Tukey test: p < 0.05 as their confidence limits overlapped) (Table 4, Table S4; Figure 2d). A group conformed by seven localities: three from Isla Mocha (IM): Punta Los Piures, Faro Viejo and Caleta Derrumbre, and four from continental territory: Valdivia, Pucatrihue, Punta Pirámide and Puerto Madryn, which showed a mean adjusted acrosome = 2.441 µm (CL = 2.423 - 2.460), which corresponded to the intermediate morphotype. The other group consisted of five localities Lebu, Mehuín, Chiloé, Comodoro Rivadavia and Punta Arenas with a mean adjusted acrosome = 2.729 µm (CL = 2.683 - 2.774) and corresponded to the long morphotype.

3.2. Molecular analyses

A total of one hundred ninety sequences from male P. purpuratus specimens collected at 21 different localities (refer to Table 1) were included in the following molecular analyses. This dataset encompasses approximately 74% of the species’ documented distributional range [11,12]. For 16S and 28S markers, an average of three and six individuals by location, respectively, were sequenced.

After alignment and trimming the ends, 66 and 124 sequences of 406 and 660-bp fragments were obtained for the mitochondrial rRNA 16S and the nuclear rRNA 28S genes, respectively.

3.3. Bayesian reconstructions and genetic structure and its relationship with acrosomal morphotypes

The trees of Bayesian analyses based on 16S and 28S genes showed two well bootstrap supported clades (> 95%) sisters to the Brachidontes rodriguezii and Ischadium recurvum outgroups (Figure 4). In both 16S and 28S sequences, the phylogenetic separations were similar and were related with the acrosomal morphotypes (i.e., north vs south clade; Figure 4). On the one hand, in the 16S/28S analyses, the north clade was formed by individuals from northern localities (Iquique to Lota, i.e., all localities with short morphotype), and by three localities from 38°S (i.e., IM Caleta Derrumbe and IM Punta Los Piures with intermediate morphotype for 16S/28S trees, plus three sequences of Tirúa in the 28S tree; see Figure 4). On the other hand, the south clade showed strong support with the rest of localities showing both intermediate and long morphotype (Figure 4). Even though the results of Bayesian inference were similar between both 16S and 28S datasets, the tree topologies showed some differences. While the 16S sequences showed two independent clades (north vs south), in the 28S analysis the north clade was nested within the basal south clade (Figure 4).

After exploring the four landscape scenarios of the Bayesian clustering analysis implemented in GENELAND, the model spatial combined with correlated frequency was selected for both, 16S and 28S datasets. The correlated frequency model has been described as a more biologically grounded way to make inferences, and it has been observed that using the correlated frequency model could be more powerful at detecting subtle admixture differentiations [38]. Between 16S and 28S sequences, the results of the Bayesian clustering analyses showed different k-groups. In the 16S dataset, the two main clusters observed (north vs south clusters, K = 2) were composed by the same sequences from BI analysis (Figure 5a). While in the 28S, three main clusters were observed (K = 3): I) the north cluster with the localities Iquique to Lota (20°S - 37°S) plus IM Caleta Derrumbe (38°S; Figure 5b-1), II) a transition cluster composed by IM Punta Los Piures and Tirúa (38°S; Figure 5b-2), and III) the south cluster with the rest of the southern localities (Figure 5b-3).

Despite of these differences, both Bayesian methods suggest divergent lineages (clades) inside P. purpuratus distribution, and a limited gene flow between these geographical regions (i.e., north, transition and south). In fact, the panmixia hypothesis was rejected in the hierarchical AMOVA (p < 0.05; Table 5) suggesting genetic structure inside the P. purpuratus geographical distribution. The AMOVA results were consistent using the geographical partitions observed in both BI and GENELAND analyses. For both, mitochondrial (16S) and nuclear (28S) markers, the percentage of variation were higher among groups (> 94%), and lower among samples within groups (< 1.94%) and within samples (< 6.6%; see Table 5). Accordingly, pairwise F_ST values were significantly different between the compared groups (Figure 6).

3.4. Genetic diversity

Among the 406 and 660 bp of the alignment length of the 16S and 28S sequences, respectively, 33 (8.1%) and 15 (2.3%) were polymorphic sites; of those, 14 and six were singleton variable sites (3.5% and 0.9%, respectively) and 19 and nine were parsimony informative sites (4.7% and 1.4%, respectively; Table 6). Of the total of the sequenced individuals, the number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (Pi) and nucleotide differences (k) were greater for the 16S rRNA than 28S rRNA datasets (Table 6). For 16S marker, the indices estimated by adjusted acrosome revealed a high genetic diversity for the intermediate morphotype (Table 6a). While the indices estimated by BI analyses showed higher values for the north clade (16S) and for the transition cluster (28S) (Table 6b,c).

Of the haplotypes identified for the 16S and 28S markers (h = 20 and 9, respectively), only three were shared haplotypes. The haplotype networks showed a non-random association which was related to acrosomal morphotypes and genetic clades, and therefore with their geographic origin (i.e., north vs south zones, Figure 4). Thus, the localities with the short morphotype showed haplotypes exclusively related to the north clade (northern haplotypes: 16S h = 2, 4, 5 and 28S: h = 1-2, Figure 4b,d); of the localities with the intermediate and long morphotypes, all (with exception of IM Caleta Derrumbe, IM Punta Los Piures and Tirúa - see Figure 4b,d) showed haplotypes from south clade. In accordance with the previous Bayesian results (Figure 4), IM Caleta Derrumbe and IM Punta Los Piures (38°S) showed the northern haplotypes in the 16S network (Figure 4b). In the 28S analysis, all sequences from IM Caleta Derrumbe showed north haplotypes, however, IM Punta Los Piures and Tirúa showed both northern and southern haplotypes (Figure 4d). Finally, in addition to the shared haplotype from south clade, two private haplotypes were observed in Chiloé (Figure 4d).

3.5. Demographic history

The neutrality tests performed by the 16S genetic clades suggested a recent population expansion (Table 7). For the south clade, all neutrality tests were statistically significant (p < 0.01; Table 7), however for the north clade, only Ramos-Onsins and Rozas’s R2 test was statistically significant (p < 0.01; Table 7). When all samples were regarded as a whole, Tajima’s D, Fu’s Fs and R2 statistic values were statistically not significant (p > 0.10). The pairwise mismatch distribution used to visualize evidence of population expansion suggested by the previous neutrality tests, showed for the north clade a bimodal distribution, which is consistent with allopatric divergence followed by population growth (Figure 7a). For the south clade, unimodal mismatch distribution represented expanding populations (Figure 7b).

4. Discussion

Along a latitudinal gradient of ca. 7,180 km of coastline, from Pacific to Atlantic Ocean, we have showed evidence of a genetic and spermatic morphotypes discontinuity at 37°S, with two distinct genetic lineages (clades) in males of P. purpuratus. The north clade with a shorter morphotype and the south clade with two different (intermediate and long) morphotypes. Although the presence of two lineages, could be indicative of cryptic species for this intertidal mussel [6,20,21], our observation of a hybridization zone at 38°S is interpreted by us as a strong signal of incipient or in progress (in fraganti) speciation processes occurring on P. purpuratus [6].

4.1. Sperm morphology

In Perumytilus purpuratus there was evidence of geographical variability in sperm sizes (for mean acrosome and head lengths), with shortest sizes typical for northern localities (up to Lota: 37°04’S), and intermediate and longest sizes typical for southern localities (south of Lota) (Table 2; Figure 2a, 2b). In statistical terms, the three morphotypes observed were well differentiated at geographical macroscale, and acrosomes better defined when were adjusted (or corrected) by total size (i.e., head length; see Figure 3c). Indeed, we observed that acrosome and head lengths were positively correlated (r = 0.89), which justified the ANCOVAs analyses. Likewise, due to that acrosomal slopes between short vs intermediate/long morphotypes (i.e., = between both northern and southern region/lineage) were significantly different, morphotypes were also better defined when ANCOVAs were performed independently for each region.

Thus, our morphological findings revealed: 1) three sperm morphotypes along Chilean and Argentinean distributional range of P. purpuratus; 2) a main morphological break between northern and southern localities situated in the Pacific Ocean at 37°S; and 3) along with intermediate morphotype, a novel morphotype, longest to the previously reported [5,6,43], observed in both Pacific and Atlantic southern localities.

Several researchers have suggested that sperm traits are species-specific and, therefore, are valuable for studying taxonomic affinities among species (e.g., [5,6,7,44,45]. In this same line, a detailed morpho-anatomical, ultrastructural, and molecular work in bivalves, highlights that sperm ultrastructure characters are among the best morphological features to diagnose bivalve clades [46]. Therefore, a relationship between sperm morphotypes and genetic structure was expected.

4.2. Linking sperm morphotypes with genetic divergence

The genetic divergence, and the presence of phylogeographic structure along the distributional range of P. purpuratus, with two well-differentiated north and south clades, confirmed previous results using mitochondrial (16S, COI) and nuclear (18S, 28S) molecular markers [6,20,21,35]. Our results also agreed with the regional genetic differentiation between P. purpuratus from southeastern Pacific (Punta de Tralca, 33° 26′S) and southwestern Atlantic (Puerto Lobos, 41° 42′S) using microsatellites markers [47]. Overall, these findings are evidence for restricted gene flow between north and south lineages.

In sessile marine broadcast spawner such as P. purpuratus, the connectivity level and genetic flow between distant populations can depends on multiple factors, as for example larval dispersal capacity, contemporary oceanography, and historical climatic events (see [21]). However, the relationship between genetic divergence and the acrosomal morphotypes here reported, with short morphotype associated to the north lineage, and both intermediate and long morphotypes associated to the south lineage, suggests that reproductive barriers, and therefore, reproductive isolation, could be the current mechanism preventing gene flow between lineages.

If we consider species-specificity of spermatozoa traits to differentiate closely related species (e.g., Acrosome size as a species conservative trait) then each morphotype must be a different species. However, our morpho-genetic results were contrasting because the two genetic lineages were associated to three sperm morphotypes.

On the one hand, the presence of two genetic clades (the present work and [6,20,21,35]) suggests the evolutionary scenario of cryptic speciation [6,21,43], that means two species with prezygotic isolation i.e., barriers that prevent fertilization, and therefore, genetic exchange between previously interbreeding populations [48,49]. However, the lack of genetic differentiation observed between intermediate and long morphotypes, being both grouped within the southern genetic clade, may also suggest an evolutionary scenario of incipient species within the southern region. In reproductive terms, this scenario likely involves postzygotic isolation, characterized by barriers occurring after zygote formation, leading to the production of non-viable or infertile offspring [48,49].

4.3. Cryptic and incipient species scenarios: an evolutionary perspective

According to our results, both cryptic and incipient species scenarios are probable. In this way, we hypothesize that short morphotype to represent a cryptic species, and that individuals with intermediate and/or long morphotypes could corresponded to incipient species.

In an evolutionary context, the association between genetic divergence and adjusted acrosomal morphotypes observed by us could be explained by historical events as the LGM [6,20,21], where glacial sheets most likely disrupt larval dispersal and, consequently, the genetic connectivity between southern and northern populations. For example, glaciation-deglaciation events have shaped the geomorphological configuration of the Chilean coast. In this context, the coastline between Perú and Canal de Chacao (north of Chiloé island, ca. 41°47′ S) is continuous, smooth and without breaks or major geographical features [50]. However, from Chiloé to Cabo de Hornos (~ 56°S), area known as the Chilean archipelago, the coastal geomorphology is complex and characterized by profusion of gulfs, fjords and channels resulting from the combined effect of tectonic processes and glaciation [50](Figure 1). During the Last Glacial Maximum (LGM) dated at 23 – 25 ka ago for Patagonia [51,52], icefields covered all southern Chile, from the Chilean Lake District (40°S) to Bahía Inútil in Tierra del Fuego (53.5°S) [51], creating a great geographical barrier or breaks in habitat continuity between southern and northern populations of P. purpuratus.

In hypothetical terms, evidence indicate that north clade to correspond to isolated populations that remained unaffected by icefields during LGM period and, therefore, they kept larval connectivity and gene flow between their populations. This could explain the lower sperm morphological variability observed for the north lineage (clade), where short morphotype remained relatively constant along ca. 1,744 km of coastline (from Antofagasta to Lota; Figure 2c). In Antofagasta area (23°S), fossil records of P. purpuratus observed in molluscan assemblages of the last interglacial period (early Pleistocene ca. 120 ka - [53,54]) support this hypothesis. Consequently, our demographic analyses using mtDNA showed a bimodal mismatch distribution for northern lineage (Figure 7a). [55] suggested that a bimodal mismatch distribution is attributable to a historically differentiated allopatric population. Besides, bimodal shape may be sensitive to the age of the expansion, with the right peak representing an older expansion and the other peak a recent expansion [56,57]. In this way, northern lineage of P. purpuratus might have a more complex evolutionary history than lineage south. Instead, the mismatch distribution for southern lineage was unimodal in shape, and closely fitted the expected distribution under the sudden expansion model (Figure 7b), that may be attributable to a more recent population expansion compared with northern lineage. This finding is consistent with analyses conducted by [20] using COI molecular marker, who (assuming a mutation rate of 0.19 substitutions/Myr), estimated different population expansion timing for each P. purpuratus lineage, with a northern expansion that could have started around 15 ka (end of the Pleistocene), and a more recent expansion for southern lineage, where the largest change in population size could have occurred during the Holocene, 11.5 to 3.5 ka BP. This could explain the sperm morphotypes variability assigned to the south clade, where the two morphotypes were distributed along ca. 5,436 km of coastline, without a significative relationship with latitude (Figure 2c).

Furthermore, our comprehensive sampling effort enabled us to precisely determine the latitudinal position of the morpho-genetic break, which was identified at 37°S on the Pacific coast. This location is notably farther north than the area covered during the Last Glacial Maximum (LGM). This evidence, and the observation of the long morphotype in the break zone (Lebu), strongly suggests postglacial recolonization of the southern lineage of P. purpuratus and supports hypothesis of ice-free refugia or suitable niches within these quaternary glacial areas [58]. Therefore, we hypothesized that in the Pacific coast, the postglacial recolonization of the south lineage of P. purpuratus reached 37°S; and that the different sperm morphotypes of the south lineage could have originated in distinct glacial refugium during the LGM period.

In other marine species inhabiting southern Chile, evidence of historical influences (i.e., LGM) on cladogenesis also have been reported, as for example in red algae [59], macroalgae [60,61], fishes [62], arthropods [63] and gastropods [58,64]. For example, in the red algae Mazzaella laminarioides, [59] also localize a genetic break at 37°S that, according to these authors, could have originated due to transient habitat discontinuities driven by episodic tectonic uplifting of the shoreline around the Arauco region (37°-38°S). Besides, three divergent lineages: northern (28°55’S to 32°37’S), central (34°05’S and 37°38’S) and southern (39°40’S to 54°03’S), and evidence of postglacial recolonization from a northern glacial refugium area was observed in this species [59]. In kelps, genetic disjunction between Patagonian (49°-56°S) and northern populations (32°-44°) was observed in Durvillaea antarctica [60]; and in Macrocystis pyrifera, [61] reported a genetic break at 42°S (Chiloé Island) and shared haplotypes among some of the subantarctic islands and southern-central Chile, suggesting a recent colonization of the subantarctic region. In P. purpuratus, shared haplotypes were observed by us at intra-clade level, that is a signal of gene flow disruption between north and south lineages. However, an exception was observed for IM Caleta Derrumbe, IM Punta Los Piures and Tirúa (localities situated at 38°S, see Figure 1), where some individuals showed the typical southern haplotypes, but others showed northern haplotypes with both mitochondrial (16S) and nuclear (28S) molecular markers. As a result, these individuals were positioned inside both north and south genetic clades (see Figure 4). We interpreted this outcome as evidence of a local hybridization zone at 38°S, suggesting the presence of an incipient speciation process within the southern lineage.

4.4. Incipient species and postzygotic isolation inside southern lineage of P. purpuratus?

Our study allowed us to determine a possible hybridization zone in the 38°S of the Pacific, nearby to the morpho-genetic boundary situated at 37°S. In the 38°S two principal localities were sampled: Tirúa (38°20′S) in the continental territory and Isla Mocha (Punta Los Piures, Caleta Derrumbe and Faro Viejo), placed 35 km in front of Tirúa, in the insular territory. With exception of IM Faro Viejo, the individuals from these localities showed haplotypes typical from both north and south genetical clades (see Figure 4 and Figure 5). It should be noted that we do not have sperm morphologic information for Tirúa, but according with our findings we hypothesize that intermediate morphotype should be present at this locality. In this way, individuals at these localities showed intermediate sperm morphotype, that was the morphotype attributed to the southern clade, however, we consider that finding as the first evidence of a secondary contact zone in P. purpuratus, with individuals showing interspecific hybridization, where northern alleles probably introgressed into gene pool of southern clade.

In theory, secondary contact zone is formed when two populations or lineages that have diverged by genetic drift or selection during a period of geographic isolation come into contact [65,66,67]. As was discussed previously, this finding is consistent with the hypothetical origin of divergence between P. purpuratus lineages, and the later demographic expansion of southern lineage until 37° south latitude. In this way, when secondary contact zones are established, genetic isolation can be maintained by prezygotic and/or postzygotic mechanisms [67].

In the case of postzygotic isolation, the mechanisms for reproductive isolation have not been completed, therefore hybrid offspring will be produced by introgressive hybridization processes; and maintained by intrinsic genetic incompatibilities or through extrinsic causes of selection against hybrids [67]. According to introgressive hybridization, the movement of alleles from one species into the gene pool of another divergent species occurs by means of repeated backcrossing of an interspecific hybrid with one of its parent species [68,69]. In mussels, hybridization zones have been well-studied, especially in the Mytilus species complex that hybridize naturally [70,71,72,73,74], and the general expectation is that hybridization could lead to offspring with reduced fitness, as for example, sterility or inferior viability, probably due to the assumptions that crosses between divergent genotypes will always disrupt co-adapted genomes; however, fitness for hybrids can range to highest to lowest (see [75]). In fact, several studies in Mytilus have showed reduced hybrid fitness, such as larval inviability and early heterosis rate [76], abnormal larvae [77,78] and high levels of sterility [79]. In our work, the lowest frequency of “hybrid” specimens could be a signal of reduced fitness.

Regarding with this, the genetical analyses performed using GENELAND package, that is suitable for hybrid zone inference [80], showed a transition cluster with IM Punta Los Piures and Tirúa (Figure 5b), i.e., the localities where individuals with both northern and southern haplotypes were observed. Whether individuals from these localities to correspond to hybrids must be determined using specific markers for introgression and hybridization detection, such SNP, which should evaluate our hypothesis of hybridization zone at 38°S.

5. Conclusions

Along a latitudinal gradient of ca. 7,180 km of coastline, from Pacific to Atlantic Ocean, our work showed correlations between the interspecific variation in spermatozoa traits and the genetic structure on males of P. purpuratus. Additionally, the evidence of genetic admixture between clades (i.e., in Isla Mocha-Tirúa) and a probable hybridization zone at this latitude (38°S) strongly suggests that P. purpuratus is under incomplete reproductive isolation, where the presence of mechanisms of reproductive isolation could be being manifested between P. purpuratus lineages, and therefore be under an incipient speciation process. In conclusion, our findings contribute to a better understanding of the evolutionary history of P. purpuratus, nevertheless, more experimental studies are needed to confirm the occurrence of incipient speciation and a hybridization zone for this species, including a careful assessment of the morpho-genetic variation, by means of additional analyses such as laboratory crosses at intra/interspecific level, and the study of reproductive-related proteins, that in an evolutionary perspective are valuable because of their contribution to reproductive isolation between species and for the formation of new species [81,82,83,84].

Finally, we consider that despite uncertainties in the taxonomy of P. purpuratus, the morpho-genetic divergence between north and south P. purpuratus lineages indicate that they do not constitute the same Evolutionary Significant Unit and, therefore, should be considered as separate entities.