When homoplasy mimics hybridization: a case study of Cape hakes (Merluccius capensis and M. paradoxus)

Sampling

Sampling took place in the years 2012–2013, from the Cunene River mouth in northern Namibia, to Port Elizabeth in South Africa, covering the entire distribution range of both species (Fig. 1). Samples were obtained from Namibian and South African governmental research surveys and commercial fisheries by trawling. A piece of muscle was collected immediately after capture and stored in 95% ethanol. DNA was extracted using a standard chlorophorm:isopropanol method (Winnepenninckx, Backeljau & Dewachter, 1993).

Preliminary data analyses

DNA of 1,200 individuals (300 per year and per species) was Polymerase Chain Reaction (PCR) amplified for a fragment of the mitochondrial DNA control region (CR—Quinteiro, Vidal & Rey-Mendez, 2000), and for ten nuclear microsatellite loci. Although more loci should ideally be used (Cornuet et al., 1999; Vähä & Primmer, 2006; Väli et al., 2010) empirical studies suggest that 100% of correct assignments can be achieved with ten loci, as long as they have an average genetic diversity of H ∼ 0.6 and 30–50 individuals per population are analysed (Cornuet et al., 1999). Moreover, the usage of ten loci allows for direct comparison with previous published work on hybridization in the marine environment (Table S1), and in particular Miralles, Machado-Schiaffino & Garcia-Vazquez (2014). Seven newly developed for M. paradoxus (MP51, MP318, MP374, MP8450, MP8478, MP8494, MP8448—Hoareau et al., 2015) and three developed for M. merluccius also used in the previous hybridization study (Mmerhk-3b, Mmerhk-20, Mmerhk-29—Moran et al., 1999).

All PCR methods followed the original protocols. The mtDNA PCR products were sequenced (CAF, Stellenbosch, South Africa) and final alignment was conducted using CLUSTAL X (Thompson et al., 1997) in GENEIOUS 7.1.4 (Biomatters, Auckland, New Zealand). Microsatellite fragments were genotyped in an ABI377 (CAF, Stellenbosch, South Africa) using LIZ as an internal size marker, and scored based on size in GENEIOUS 7.1.4. Accuracy of scores was ensured by using a reference individual as positive control. The microsatellite dataset was evaluated for quality of amplification by estimating the occurrence of null alleles, large allele drop out and stuttering in MICROCHECKER (Van Oosterhout, Weetman & Hutchinson, 2006). Deviations to the expectation of outcrossing and linkage disequilibrium were tested in FSTAT (Goudet, 1995). Statistical significance was tested with 10,000 permutations. One locus was found to be in linkage with locus MP374, and the latter was therefore removed from further analyses. The possibility of temporal differentiation was tested by performing pairwise genetic divergence analyses between sampling years for both species for the CR dataset in ARLEQUIN (Excoffier, Laval & Schneider, 2005), and for the nine microsatellite loci in FreeNA (Chapuis & Estoup, 2007), with statistical significance assessed after 10,000 iterations. As no evidence of genetic differentiation was observed, sampling years were pooled together and treated as one single population. Genetic diversity levels were measured as haplotype (h) and nucleotide (π) diversity for mtDNA in ARLEQUIN, and expected (H_E) and observed heterozygosity (H_O), number of alleles (Na), allelic richness (AR), percentage of null alleles (NAL) and inbreeding index (F_IS) for the microsatellite dataset in ARLEQUIN, FSTAT and FreeNA. Deviations of neutral expectations were tested within geographical regions (Namibia vs. South Africa) and between species (M. capensis vs. M. paradoxus) using the F_ST-based approach implemented in LOSITAN (Antao et al., 2008), under the Stepwise Mutation Model (SMM) and run for 50,000 simulations.

The ability of the microsatellite dataset to distinguish between M. paradoxus and M. capensis was investigated by: (a) constructing a distribution of allelic frequencies, to assess the presence and frequency of fixed alleles, (b) estimating pairwise F_ST values in FreeNA (Chapuis & Estoup, 2007), with statistical significance assessed after 10,000 iterations, and (c) performing a Principal Component Analysis (PCA) based on allelic frequencies in adegenet 1.3.1 (Jombart & Ahmed, 2011; this method is less influenced by deviations to Hardy-Weinberg equilibrium and the presence of null alleles). Furthermore, as different age classes of both species are commonly caught together and external morphological identification is difficult, preliminary analyses of individual assignment were conducted using both mtDNA and nDNA datasets to validate morphological identification. Medium-joining networks were reconstructed for the CR sequences of both species combined in NETWORK (Bandelt, Forster & Rohl, 1999). The coalescent-based assignment software STRUCTURE (Pritchard, Stephens & Donnelly, 2000) was used to assign individuals to species, based on microsatellite data. Five preliminary runs for K = 2 were performed under a strict assignment model, with no-admixture and independent allelic frequencies, for 250,000 MCMC iterations, after an initial burnin of 50,000 MCMC steps. The five runs were combined in CLUMPP (Jakobsson & Rosenberg, 2007), and final results visualized in DISTRUCT (Rosenberg, 2004).

Assignment analyses and identification of hybrid states

Assignment tests and identification of hybrids were performed in a two-step manner. First, a simulation study was conducted in Hybridlab (Nielsen, Bach & Kotlicki, 2006) to assess if the microsatellite dataset could accurately differentiate among multiple hybridization events, i.e. that hybrids were viable and capable of backcrossing with their parental species. Four hybrid states were simulated: first generation (F₁), second generation (F₂), backcrosses of F1 with M. paradoxus (BcP) and backcrosses of F1 with M. capensis (BcC). A total of 80 simulated hybrids were generated based on the allelic frequencies of 160 randomly chosen individuals identified as “pure” M. paradoxus and M. capensis in the previous STRUCTURE analysis, using a cut-off value of 0.9 for the ancestry coefficient (q). This simulated dataset was then used in Bayesian clustering analyses implemented in STRUCTURE and NEWHYBRIDS (Anderson & Thompson, 2002). Both approaches rely on the assessment of an ancestry coefficient (q) for each individual based on the allelic frequencies of the identified clusters. We ran STRUCTURE analyses for four different admixture models mimicking a range of evolutionary scenarios: (i) no interbreeding and two independently evolving species (no admixture and independent allelic frequencies); (ii) no interbreeding, but species sharing a recent common ancestor (no admixture and correlated allelic frequencies); (iii) interbreeding but two distantly related species (admixture and independent allelic frequencies) and (iv) interbreeding between recently evolved species/populations (admixture and correlated allelic frequencies) (Pritchard, Stephens & Donnelly, 2000). Each model was assessed with five preliminary runs, for K = 2, with an initial burnin of 50000 MCMC steps, followed by 250000 MCMC iterations. The five runs were combined in CLUMPP, and results visualized in DISTRUCT. There is no general consensus regarding the most likely q to distinguish between “pure” and admixed individuals, with the choice of q greatly depending on the purpose of the study (Vähä & Primmer, 2006). A q = 0.2 (20% of admixture) increases accuracy: the ability to correctly assign individuals to their classes (“pure” vs hybrids), while q = 0.1 (10% of admixture) will decrease the probability of incorrectly assigning “pure” as hybrids (Vähä & Primmer, 2006). As one of the main aims of this study is to detect admixed individuals and assess if the microsatellite dataset can accurately distinguish between different hybrid states, we employed both q levels to understand how this metric influences the number of putative hybrids identified (Bohling, Adams & Waits, 2013; Vähä & Primmer, 2006; Väli et al., 2010). Only individuals with non-overlapping credibility intervals were retained as putative hybrids. Finally, identification of hybrids states was conducted in NEWHYBRIDS. Five runs were performed with different initial seed values, using a Jeffrey prior, for 250,000 MCMC iterations, after an initial burnin of 50,000 MCMC steps. Runs were combined in CLUMP and individuals identified based on q = 0.5 (Anderson & Thompson, 2002). The accuracy of both approaches was recorded as a percentage of correct assignments to each state. Once the simulation study was performed, the same analyses were conducted for the true microsatellite dataset.

Preliminary analyses revealed significant deviations to the expectations of outcrossing, mainly in M. capensis, the presence of null alleles for both species and deviations from neutrality in three loci (MP318, MP8450 and Mmerhk-3b). Although these assignment tests are generally robust and can be performed even in non-ideal conditions (Carlsson, 2008; Cornuet et al., 1999; Hauser et al., 2006), we created two additional sub-datasets: (a) removing loci with a null allele frequency above 15% and (b) without the three loci under selection, in order to understand how these features influence the identification of putative hybrids in the Cape hakes. Both datasets were run in STRUCTURE (only with model iii—admixture with independent allelic frequencies for computational purposes) and NEWHYBRIDS, with the same run conditions described above.

Verification of hybrids

A total of 38 “pure” individuals and 19 putative hybrids identified with more than one assignment method were further PCR amplified and sequenced for Mmer-hk3b, following the original protocol (Moran et al., 1999). At the microsatellite level, the sequenced hybrids were either heterozygous, exhibiting allele sizes from both species (e.g., 327/341—M. capensis/M. paradoxus), or homozygous for the alleles of the other species (M. capensis: 341/341; M. paradoxus: 327/327). Furthermore, individuals with mito-nuclear discordance were re-sequenced for mtDNA, to confirm their status. The obtained sequences were aligned using CLUSTAL X, in GENEIOUS 1.7.4, and used for hybrid identification. Individuals with admixed origin are expected to be heterozygous for sites with species-specific mutations, and alleles of F1 individuals should thus appear in both species specific clades after a gene tree is constructed for this data. By comparing direct sequences of the flanking regions, it is also possible to identify potential homoplasy due to insertions or deletions either in the core sequence of the microsatellite, or in the flanking region (Estoup, Jarne & Cornuet, 2002; Grimaldi & CrouauRoy, 1997; Van Oppen et al., 2000). In order to understand if admixed individuals resulted from interbreeding, incomplete lineage sorting or homoplasy, we used the available sequence of Mmerhk-3b for M. merluccius in Genbank (Accession number: AF136627.1) in the phylogenetic analyses. Reconstruction of phylogenetic relationships was performed with the Maximum Likelihood algorithm implemented in PhyML (Guindon et al., 2009), using the most suitable nucleotide substitution model as indicated by jModelTest (Posada, 2008). Nodal support was obtained employing the X²-arLT statistics (Anisimova & Gascuel, 2006), and trees visualized in Figtree (Rambaut, 2009).

Preliminary analyses

Of 1,181 individuals, 1,137 were successfully PCR amplified and sequenced/genotyped for both mtDNA and the microsatellite datasets. As extensive introgression (via mito-nuclear discordance) was described in Miralles, Machado-Schiaffino & Garcia-Vazquez (2014), only individuals amplified for both datasets were retained for further analyses.

Assessment of amplification quality of microsatellite loci revealed that although there was no evidence of large allele drop out or stuttering, both species’ datasets contained null alleles and did not conform to Hardy-Weinberg expectations, due to heterozygotes deficit (Table 2). Overall, M. paradoxus had a lower frequency of null alleles (∼3%) at the cross-specific loci Mmerhk-20 and Mmerhk-29, and at locus MP8450. In contrast, the null allele frequencies in M. capensis ranged from 1–22.9% across six loci (Table 2). Locus MP8448 and Locus MP374 were discovered to be in linkage, and thus the first was removed from the final dataset. All subsequent analyses were performed with nine microsatellite loci. Neutrality tests performed within and between species revealed outlier markers, which appeared to be under positive selection. Loci MP318, MP8450 and Mmerhk-3b were flagged as outlier loci between M. paradoxus and M. capensis (F_ST > 0.1, p < 0.05). Loci MP8894 and Mmerhk-3b were revealed to be non-neutral in M. capensis (F_ST > 0.1, p < 0.05). No significant deviations to neutrality were observed for the M. paradoxus dataset.

A fragment of 406 bp was sequenced for the CR, yielding ten haplotypes in M. paradoxus, and 160 haplotypes in M. capensis. Genetic diversity levels varied with species and dataset used (Tables 1 and 2). Overall, M. paradoxus had lower genetic diversity for the mtDNA and higher for the microsatellite dataset: h = 0.541, π = 0.002, H_O = 0.689; while M. capensis showed an inverse pattern with higher genetic diversity observed for the mtDNA dataset and lower for the microsatellite loci h = 0.893, π = 0.008, H_O = 0.619 (Tables 1 and 2). Number of alleles ranged from 5–45 per locus in M. paradoxus and 5–33 per locus in M. capensis (Table 2). No significant genetic differentiation was observed between sampling years for either dataset or species (M. paradoxus: φ_ST = 0.002, F_ST = 0.000, p > 0.05; M. capensis: φ_ST = 0.001, F_ST = 0.003, p > 0.05), and the presence of null alleles did not significantly influence differentiation measures in M. capensis (uncorrected F_ST = 0.0034, corrected F_ST = 0.0032). Thus, both sampling years were pooled.

The microsatellite dataset could accurately distinguish between the two species with different loci exhibiting non-overlapping allele size ranges for M. paradoxus and M. capensis (Fig. S1). In particular, locus Mmerhk-3b was the most divergent, with two sets of alleles fixed in each taxon, and only a few individuals of either species exhibiting shared alleles (Fig. S1). The PCA plot clearly separated the two species along the first two axes, with limited overlap between them (Fig. 2). Pairwise F_ST values were estimated at F_ST = 0.126 (uncorrected, p < 0.05) and F_ST = 0.114 (corrected, p < 0.05) between M. capensis and M. paradoxus.

Preliminary assignment tests using both datasets revealed the presence of 47 individuals that were misidentified during sampling (4.13% of total samples): 42 for M. capensis and five for M. paradoxus, which were reassigned to represent the correct species (Fig. S2). Reconstruction of the phylogenetic relationships between individuals clearly separated M. capensis and M. paradoxus with only two individuals identified as having mito-nuclear discordance (after correcting for sampling mislabels; Fig. S2).

Assignment analyses and identification of hybrids

Simulation analyses showed that identification of simulated hybrids significantly varied with the model and approach used (Fig. 3). Assignment analyses performed in STRUCTURE were able to detect F1 hybrids in the majority of cases, but accurate identification of further hybrid states varied significantly with the chosen admixture model (Fig. 3). Overall, assignment models iii and iv identified a higher number of hybrids than models i and ii, while q = 0.1 had a more accurate detection of multiple hybridization events (BcC and BcP) than q = 0.2 (Fig. 3). Successful detection of simulated F2 hybrids ranged from 21.3% (q = 0.2, model i) to 93.7% (q = 0.1, model iv); BcC hybrids ranged from 6.3% (q = 0.2, model i) to 33.8% (q = 0.1, model iv) and BcP hybrids ranged from 6.3% (q = 0.2, model i) to 33.8% (q = 0.1, model iii). However, while individuals considered “pure” had non-overlapping credibility intervals, implying that the nine microsatellites can accurately assign individuals to either species, the majority of identified hybrids exhibited overlapping intervals regardless of the method employed. The NewHybrids approach performed better and allowed to correctly identify the hybrid status in 96% of the simulated cases for F1, 87.5% for F2, 85% for BcC and 87.5% for BcP (Fig. 3).

Similarly, analyses of the complete dataset revealed that the number of putative hybrids identified varied significantly with the admixture model used and the ancestry coefficient threshold (Table 3 and Fig. 4). As expected, q = 0.1 identified a higher number of admixed individuals when compared to q = 0.2: 46 vs. 16 (Table 3). Assignment models i and ii (no admixture) identified up to 14 putative hybrids (1.14%), while models iii and iv (admixture) retrieved up to 42 hybrids (3.78%) (Table 3 and Fig. 4). The NEWHYBRIDS approach identified 21 individuals as F2 (1.93%), one M. capensis as pure M. paradoxus and the remaining individuals as pure species (Table 3). Only three individuals were identified as putative hybrids by all models and approaches (0.27%), 19 (1.67%) by two or more models and of these only six had non-overlapping Bayesian credibility intervals (Table 3). The majority of admixed individuals were M. capensis (26), and all hybrids were found throughout the distribution range of the species, from northern Namibia to the East Coast in South Africa (Table 3).

Assignment tests using the two additional datasets identified different numbers of putative hybrids (Table S3). While dataset 2 (8 loci, null allele frequencies <15%) identified approximately the same individuals as before (Table S3), removing the two loci with fixed differences between species (dataset 3, 6 loci) resulted in a significant increase of the number of admixed individuals identified (Table S3), resulting in a decrease in accuracy, with a higher number of “pure” individuals classified as hybrids (Table S3). However, out of the 108 individuals identified as hybrids, only four exhibit non-overlapping Bayesian credibility intervals. Therefore, the majority of admixed individuals identified reflected the lack of statistical power of the microsatellites.

In the real dataset, only two individuals appeared to exhibit a mito-nuclear discordance: one M. capensis sample appeared to have the nuclear genotype of M. paradoxus and mtDNA haplotype of M. capensis, and another M. capensis sampled was identified as M. paradoxus based on mtDNA, but had the genotype of M. capensis (Table 3).

Verification of hybrids

Amplification and sequencing of the flanking region of the microsatellite Mmerhk-3b resulted in a fragment of 190 bp for 38 pure individuals (20 M. paradoxus and 18 M. capensis), with five fixed nucleotide differences and one 10 bp indel observed between the two species (Table 4). The sequence from M. merluccius exhibited the same 10 bp deletion than found in M. capensis, and had a mix of the fixed positions between M. capensis and M. paradoxus (Table 4). A total of 19 putative hybrids, identified based on multiple methods and models, were sequenced for the same region (including the two samples with mito-nuclear discordance). No heterozygotic sites were observed for the five fixed positions that differentiate the two hake species, and reconstruction of phylogenetic relationships, using K80 as the model for nucleotide substitution, retrieved two monophyletic clades, one for M. paradoxus and one for M. capensis (Fig. 5). In particular, individuals that had allele 227 (M. capensis) and 241 (M. paradoxus) showed no heterozygote sites in the flanking region. The two M. capensis individuals identified as mito-nuclear discordances grouped either with M. paradoxus (13_C9SW_15) or with M. capensis (13_C104W_5) for both re-sequenced mtDNA and Mmerhk-3b markers, suggesting that the original discordance was likely the result of a mislabelling error (13_C9SW_15) and microsatellite artefact (13_C104W_5).

Interbreeding, ancestral polymorphism or homoplasy?

The microsatellite dataset used in this study revealed the presence of multiple individuals with admixed ancestry, 19 of which were consistently identified as F₂ hybrids by different models and approaches. The majority of hybridization studies would thus accept these 19 individuals as hybrids, with an interbreeding rate of 1.58% between M. paradoxus and M. capensis. However, there is a distinct caveat in using microsatellites when making such assumptions in species that have evolutionary divergent histories, as microsatellite loci do not evolve in the same manner as other genomic sequences (Calabrese & Sainudiin, 2005; Ellegren, 2004). For example, a study between two sister-species of coral trouts (Plectropomus) revealed a complex pattern of hybridization, from ancestral introgression, and incomplete lineage sorting at one microsatellite locus, to complete reproductive isolation (Van Herwerden et al., 2006). In this case, the two species had been isolated for only 500,000 years. In another case, cichlid species isolated for 700,000 years exhibited evidence of homoplasic mutations in one locus (Van Oppen et al., 2000). These cases illustrate the need to assess the accuracy of microsatellites for distinguishing between interbreeding, retention of ancestral polymorphisms through incomplete lineage sorting and homoplasy.

The sequencing of the most divergent microsatellite loci (Mmerhk-3b) revealed five fixed nucleotide substitutions (all transitions) and one 10 bp insertion in M. paradoxus. Reconstruction of phylogenetic relationships clearly distinguished between the two hake species, and all putative hybrids were assigned to either clade, with no heterozygotes observed. These findings support a hypothesis of no hybridization between M. paradoxus and M. capensis, corroborating earlier allozyme studies (Grant & Leslie, 2001). We then used the available sequence of Mmerhk-3b of M. merluccius, the sister-species of M. capensis, in order to understand if retention of ancestral polymorphism could explain the presence of same-size alleles in the two species. Merluccius merluccius grouped with M. capensis, shared nucleotide changes with both species, but did not have the 10 bp insertion observed in M. paradoxus. Rare genomic changes such as indels are considered good indicators of phylogenetic relationships, as their presence is highly conserved across groups of animals (Matthee et al., 2007). The absence of the 10 bp insertion in both M. merluccius and M. capensis suggests that these species share a common ancestor (as observed by their close phylogenetic relationship). Therefore, it is not likely that the presence of same-size alleles in M. paradoxus and M. capensis results from the retention of ancestral polymorphisms in this case since they share no homologous alleles at this locus (Campo et al., 2007; Grant & Leslie, 2001; Quinteiro, Vidal & Rey-Mendez, 2000). It is more likely that complete lineage sorting has occurred, at least for this marker. Furthermore, as the changes in the flanking region were observed in both putative hybrids and pure individuals, it is not likely that the alleles were obtained via heterospecific crossing. These findings point to size homoplasy as the most probable mechanism behind the identification of putative hybrids between these species.

Homoplasy in microsatellites is generally expressed as changes in the repeat unit (point mutations, compound microsatellites) or mutations in the flanking region (Angers & Bernatchez, 1997; Estoup, Jarne & Cornuet, 2002; Grimaldi & CrouauRoy, 1997; Van Oppen et al., 2000), as observed in Mmerhk-3b. Nevertheless, the use of only one marker under selection (Hoareau et al., 2015), to make inferences about genome-wide homoplasy has in itself several drawbacks. First, recombination and unequal gene conversion can, by chance, mimic the effects of homoplasy and thus influence the conclusions (Balaresque et al., 2014; Ziegler et al., 2009). However, the mutation rate within the repeat unit is generally higher than the recombination rate between the core sequence and the flanking region (10⁻⁶ for sequences 100 bp long—Hilliker, Clark & Chovnick, 1991), making it unlikely for recombination or gene conversion to be responsible for the observed differences (Estoup, Jarne & Cornuet, 2002). Second, the ability to detect multiple generation hybrids using nuclear sequences is somewhat limited, as the probability of identifying them decreases with an increase in repeated events (F1–100%; F2–33%; backcrosses–16.67%). For that reason, we specifically targeted individuals that were either heterozygotes or homozygotes for alleles of the other species. If these admixed individuals resulted from heterospecific crossings then we would expect the flanking region to either show heterozygotes for the fixed positions (as they had one allele from each species), or the diagnostic SNPs from the other species (for homozygotes). That was not the case, which further supports the hypothesis of homoplasy instead of hybridization between the Cape hakes. Finally, Mmerhk-3b exhibited significant deviations to neutrality, and thus selection against hybrids might have influenced our findings (by removing heterozygotes from the gene pool). As mentioned above, we selected individuals exhibiting the alleles common to the other species, and so, would still expect to retrieve heterozygotes for the fixed positions. Therefore, it is likely that the observed changes in Mmerhk-3b indeed reflect homoplasy and not heterospecific crossings between the two Cape hake species.

Size homoplasy has long been known to affect microsatellite loci, as it is deeply linked to the most common mutational model for these markers (SSM: Single Stepwise Model) (Angers & Bernatchez, 1997; Estoup, Jarne & Cornuet, 2002; Grimaldi & CrouauRoy, 1997; Van Oppen et al., 2000). In addition, evolutionary factors such as mutation rate, effective population size and time since divergence can also influence homoplasic loci (reviewed in Estoup, Jarne & Cornuet, 2002). Empirical studies suggest that for an average mutation rate of 5 × 10⁻⁴ all modelled alleles are homoplasic after 6,000 generations, regardless of effective population sizes and mutational models (Estoup, Jarne & Cornuet, 2002). With an average generation time of 3.5 years (Botha, 1985) 6,000 generations would correspond to 21,000 years for the Cape Hakes. As the species have been isolated for 3–4 Ma, it is thus likely that similar allelic forms result from widespread homoplasy, and not heterospecific crossings (Estoup, Jarne & Cornuet, 2002).

Nevertheless, a large number of studies in population genetics in the marine environment continues to rely on these markers to describe hybridization and often the analytical methods to make inferences are not fully explored (Table S1). This is particularly concerning when the species are commercially exploited since accurate evolutionary histories are required for adequate management. In conclusion, our study provides empirical data showing that homoplasy has likely mimicked the effects of heterospecific crossings, and resulted in false positive signals of hybridization between Cape hakes. This has important implications for the management of these two co-distributed species, as hybridization plays no obvious part in the genetic structuring of these species.