Patterns of admixture and introgression in a mosaic Mytilus galloprovincialis and Mytilus edulis hybrid zone in SW England

The study of hybrid zones offers important insights into speciation. Earlier studies on hybrid populations of the marine mussel species Mytilus edulis and Mytilus galloprovincialis in SW England provided evidence of admixture but were constrained by the limited number of molecular markers available. We use 57 ancestry‐informative SNPs, most of which have been mapped genetically, to provide evidence of distinctive differences between admixed populations in SW England and asymmetrical introgression from M. edulis to M. galloprovincialis. We combine the genetic study with analysis of phenotypic traits of potential ecological and adaptive significance. We demonstrate that hybrid individuals have brown mantle edges unlike the white or purple in the parental species, suggesting allelic or non‐allelic genomic interactions. We report differences in gonad development stage between the species consistent with a prezygotic barrier between the species. By incorporating results from publications dating back to 1980, we confirm the long‐term stability of the hybrid zone despite higher viability of M. galloprovincialis. This stability coincides with a dramatic change in temperature of UK coastal waters and suggests that these hybrid populations might be resisting the effects of global warming. However, a single SNP locus associated with the Notch transmembrane signalling protein shows a markedly different pattern of variation to the others and might be associated with adaptation of M. galloprovincialis to colder northern temperatures.

. Hybrid zones are considered excellent natural scenarios for investigating the evolutionary process and understanding how genetic and environmental factors contribute to reproductive isolation (Gompert et al., 2017;Taylor & Larson, 2019).Knowledge of hybridization and speciation is less advanced in marine compared with freshwater or terrestrial ecosystems but is increasing with recent developments in genomic analysis and the integration of genomic results with other research fields (Faria et al., 2021;Gardner, 1997).
Hybrid zones may be of particular value in assessing the effects of climate change which may modify opportunities for hybridisation and affect the relative fitness of parental species and hybrids (Chunco, 2014;Taylor et al., 2015;Wielstra, 2019).New adaptive alleles introduced by hybridisation and introgression might also help species deal with climate change through evolutionary rescue resisting both extinction and hybrid zone movements (Baskett & Gomulkiewicz, 2011;Becker et al., 2013;Brauer et al., 2023).
Hybrid zones in mussels are characterized as mosaic with parental and admixed populations mixed up in a patchy distribution (e.g.Bierne, Borsa, et al., 2003;Fraïsse et al., 2014;Skibinski et al., 1983).Mosaic hybrid zones whose position might be influenced by climate were observed in a variety of species including insects, amphibians, plants, and birds (Taylor et al., 2015), raising the possibility that mussel hybrid zones might be similarly affected.In the present study, we use 57 ancestry-informative SNP markers to analyse interbreeding and admixture between M. edulis and M. galloprovincialis from the SW England hybrid zone.We demonstrate significant differences in genetic structure between hybrid populations and provide support for two distinct hybrid zone regions in SW England.M. galloprovincialis has higher viability than M. edulis consistent with earlier studies for other fitnessrelated traits.By contrast, the results suggest asymmetrical gene flow from M. edulis into M. galloprovincialis.Despite these factors and warming sea temperatures, the zone appears not to have moved over the 35 years that it has been studied genetically.One SNP locus, situated in the gene Notch coding for a transmembrane signalling protein, showed a pattern of geographic variation markedly different from the others and might reflect introgression from M. edulis into M. galloprovincialis and local adaptation.

| Sampling of mussels, tissue collection and phenotypic traits
Marine mussels (Mytilus spp.) were collected in March 2016 from the intertidal rocky shore in south-west (UK) England for six populations of Mytilus edulis and Mytilus galloprovincialis including localities where these species live in sympatry and hybridize (Skibinski, Ahmad, et al., 1978;Skibinski et al., 1983) (Figure 1).Two of these, Bude and Padstow, are considered regional reference populations for M. galloprovincialis.The other four are mixed populations of both species.Samples were also taken from two allopatric reference populations, Swansea in South Wales (M.edulis) and Vigo in Spain (M.galloprovincialis).The word 'reference' is not intended to imply genealogical relationships between populations.For some analyses, additional mussels from spawning experiments and also from published studies were included to increase power (see File S1).
Mussels were transported to the marine station at the University of Vigo, Spain (ECIMAT, CIM-UVigo) for processing.Mussels were sorted according to shell length, into small (20-29 mm) and large (30-40 mm) classes.This size range was used because in the UK mussels <20 mm are usually immature and cannot be assayed easily for gonad development stage, and 40 mm represents an upper size limit naturally.Size classes were used, rather than size itself as a continuous variable, to allow comparison with early studies where such classes were used (see Section 4).For each population, 60 individual mussels (50% for each size class) were taken randomly.Mantle edge tissue was dissected and preserved in 95% ethanol for DNA extraction and genetic analyses.Gonad tissue was collected for histological analysis.
Gonad development stage was assigned to seven categories (0-6) using a histological test based on haematoxylin-eosin stain (Seed, 1969), with the classification scheme of Kenchington et al. (2020).Stages 0-3 do not exhibit mature gametes and were pooled together for analysis (abbreviated 'Not Ripe') as were stages 4-6 which exhibit mature gametes (abbreviated 'Ripe').Sex was de-
The alleles at each locus were recoded as G or E from the original DNA base assignation according to whether they were at higher frequency in M. galloprovincialis or M. edulis, respectively.The sequences amplified by KASP™ which contained the 57 SNPs were used for blastX against the nr NCBI protein database.This was done to assess whether the KASP™ sequences were located in coding or non-coding regions, and whether the SNP substitutions were synonymous or non-synonymous (File S2-wsA and wsB).

| Data analysis
ENTROPY v.2.0 (Gompert et al., 2014;Shastry et al., 2021), a hierarchical Bayesian modelling software, was run to estimate ancestry coefficients q and Q 12 for each individual over loci using genotype likelihoods assuming k = 2, which corresponds with the two Mytilus species.1,000,000 MCMC steps following a burn-in period of 100,000 steps were used, with all other parameters as the default values.The coefficient q estimates the admixture proportion with reference to one of the two species and Q 12 estimates the interspecific heterozygosity (also called interspecific ancestry).In a non-modelling context, q is analogous to hybrid index, the frequency of the diagnostic alleles for one of the species, and Q 12 is analogous to the frequency of heterozygotes for the diagnostic alleles.'Triangle' plots of interspecific heterozygosity (Q 12 ) against admixture proportion (q) are informative in assessing genetic structure (Barton, 2000;Fitzpatrick, 2012;Gompert & Buerkle, 2016).
A genetic linkage map of 36 SNP markers was constructed for Mytilus by Simon et al. (2020).They comprise six groups of linked markers and a further nine unlinked markers.Thirty-two of those SNP markers were used in the present study (File S3-wsA).As admixture proceeds, linkage disequilibrium decays at a rate which is a function of the recombination frequency (Barton, 2000;Barton & Hewitt, 1985).Using the map distances, STRUCTURE vs.2.3.4 program (Pritchard et al., 2000) was used to estimate the number of generations since a recent admixture event in hybrid populations following a linkage model (Falush et al., 2003).Swansea and Bude were used as reference populations for M. edulis and M. galloprovincialis, respectively.The POPFLAG option was set as ON.Run parameters were set to k = 2, assuming a linkage and correlated allele frequencies model, with 100,000 burn-in and 1,000,000 iteration steps.The results of this analysis are supported by simulation studies of interbreeding between M. edulis and M. galloprovincialis using a program (HybSim) written in Fortran 95 (File S4).
NewHybrids vs.1.1 (Anderson & Thompson, 2002) was used to calculate the posterior probabilities that individuals in the population samples belong to different genotypic frequency categories for two generations of interbreeding (G2).They are P0 (M.edulis), P1 (M.galloprovincialis), F1, F2, and backcrosses BC0 and BC1 and are distinguished by having different expected genotypic proportions for EE, EG, GE, and GG (File S3-wsB).They are distinct from genealogical classes which can be named p0, p1, f1, f2, bc1, and bc2.Thus, a cross between an M. galloprovincialis parent (p1) and a first-generation hybrid (f1) would produce a backcross individual (bc1).An individual belonging to a particular genealogical class (e.g.bc1) might not generate a high posterior probability to the corresponding genotypic frequency category (i.e.BC1), but to another.This may happen when the individual belongs to a genealogical class whose corresponding genotypic frequency category is not amongst those presented in the analysis.Individuals were assigned to the genotypic category for which they had the highest posterior probability.The application to three generations of interbreeding (G3) was also investigated.In all the NewHybrids runs, Jeffreys-like priors for both mixing proportions and allele frequencies, and 100,000 burn in and 1,000,000 iterations were used.Unless stated otherwise results are presented by default for the G2 analysis.Simulation was also carried out to assess results for three generations of interbreeding using the program HybSim (File S4).Bayescan v.2.1 (Foll & Gaggiotti, 2008) was used to check for loci with outlying F ST values.Prior odds for the neutral model of 10, the default, were used.Because the SNP loci were chosen a priori because they discriminated between M. edulis and M. galloprovincialis, the analysis was applied separately only to individuals classified as P0 and P1 in the G2 NewHybrids analysis from the six populations in the hybrid zone (Swansea and Vigo excluded).A test for loci showing atypical patterns of variation was also carried out using the Mahalanobis distance.The calculation was carried out in SPSS vs.28, and the distance was converted to a p-value which in turn was converted to a q-value (Storey, 2002) employing the Bioconductor vs.2.28.0 q-value package.
Contingency analyses were used to examine associations between NewHybrids genotypic categories, populations and the measured phenotypic characters (sex, mantle edge, and gonad colour and gonad developmental stage) using SPSS vs.28 with Fisher's Exact test with the Monte Carlo method based on 100,000 sampled tables.Significance of individual cells was determined by converting the adjusted standardized residual to a two-tail p-value using the cumulative distribution function (CDF) (function norm.s.dist in Excel).
As a measure of association, Cramer's V test was calculated for each contingency table .F IS and p-values were calculated in Genepop vs.4.7.5 (Rousset, 2008).Two four-point correlation coefficients (R1 and R2) were computed from the observed numbers of the nine possible twolocus genotypes formed from diagnostic E and G alleles.R1 measures the association of the four two-locus homozygote genotypes.

| Allele and genotype proportions
Raw genotypic data for each SNP marker and for Glu-5, KASP™ sequences, blastX results, information for the phenotypic traits studied, derived statistics including allele and genotype frequencies, and deviations from Hardy-Weinberg equilibrium (F IS ) are provided in supplementary material (File S2-wsA to wsC and wsE to wsM).
The EE, EG, and GG genotype proportions, G allele proportion, and F IS were computed for each SNP over all mussels within each population.These summary values were averaged over all SNP loci (Figure 2).The M. edulis reference population (Swansea, SW) has much lower G allele frequency than the M. galloprovincialis reference populations Bude (BD), Padstow (PD), and Vigo (VG).The four populations Croyde (CR), Porthcurno (PC), Carlyon (CB), and Whitsand (WB) referred to as hybrid populations are intermediate.Deficit of heterozygotes (positive F IS ), consistent with the presence of two types between which there is a restriction on interbreeding, is higher for the hybrid populations and highly significant (File S2-wsF to wsM).The DNA marker Glu-5 shows a pattern across populations which correlates with the mean G allele proportion across the SNP loci (R 2 = 0.997, File S2-wsQ).
There are some noticeable population-dependent differences between large and small mussels.For example, the G allele is greater in frequency in large mussels in the four hybrid populations.At Croyde, F IS is higher in large than small mussels, whereas F IS is higher

| Interspecific heterozygosity and admixture proportion
Values of admixture proportion (q) and interspecific heterozygosity (Q 12 ) (File S2-wsE) were used to construct triangle plots for individual mussels, coloured according to the G2 NewHybrids genotypic frequency categories (Figure 3 have a wider range of q-values than at Vigo.These might be the result of a third generation of interbreeding.This is supported by analysis consistent with the occurrence of individuals belonging to G3 NewHybrids categories and supported by simulation results (Files S1 and S3-wsP).Another possibility is that allopatric reference populations outside the hybrid zone (e.g.Vigo) are less introgressed with M. edulis genes than those within the hybrid zone (e.g.Bude).This might result in the identification of a higher proportion of admixed individuals than would have been observed if local reference populations had been used (Bierne, Borsa, et al., 2003;El Ayari et al., 2019).This is tested by repeating the NewHybrids analysis excluding the Vigo sample (File S1).
There is a small reduction in the number of admixed individuals, but the conclusion that there are many admixed individuals still stands.
An alternative perspective tested is that admixture would be observed to be greater in UK M. galloprovincialis if only Mediterranean M. galloprovincialis were included in the dataset as reference populations outside the UK (Files S1 and S3-wsM).A greater number of admixed individuals are detected, but a large proportion are still identified as M. galloprovincialis (P1).In an additional test, STRUCTURE was used to assess similarity of UK M. galloprovincialis samples to Atlantic and Mediterranean M. galloprovincialis using seven of the 57 SNP loci that discriminate between the latter two.The results show a high similarity of UK to Atlantic rather than to Mediterranean M. galloprovincialis (File S3-wsN).
NewHybrids was developed primarily for unlinked markers (Anderson & Thompson, 2002).Analysis using 15 unlinked markers from amongst the 57 studied (File S3-wsA) and derived from Simon et al. (2020) give results consistent with those obtained for the 57 loci (File S1-wsG).
In the triangle plots, there are more BC1 individuals on the lefthand side than BC0 individuals on the right-hand side (Figure 3).
The greater frequency of BC1 suggests greater backcrossing to M.
F I G U R E 3 Plots of interspecific heterozygosity (Q 12 ) against admixture proportion (q) for samples from eight populations.The plotted points for individual mussels are coloured according to the NewHybrids G2 analysis genotypic frequency categories.
galloprovincialis than to M. edulis.As BC1 has a significantly higher frequency in the large than small-size class, a contributing fact might be higher viability of BC1.

| Extent of admixture in SW England
Estimates of the number of generations since admixture from STRUCTURE range between 0.2 generations at Croyde and 4.3 generations at Whitsand with Porthcurno (2.0) and Carlyon (2.3) intermediate.These low values could suggest ongoing interbreeding between parental M. edulis and M. galloprovincialis to repeatedly generate individuals from the early stages of admixture such as F1, G2, and G3.This is supported by simulation results demonstrating that when there are high numbers of parental types, few individuals that are admixed beyond G3 are produced (File S3-wsO).An equilibrium is established in which highly admixed individuals are at very low frequency or absent.Ongoing early-stage admixture is supported by a pattern of decline in correlation between loci measured by R1 and R2 as the map distance between the SNP loci increases (File S3-wsK).At lower map distance, the values are high suggesting that there has not been sufficient time for breakdown of linkage disequilibrium (R1).There are distinctive differences between the hybrid populations indicative of differences in the amount of admixture.For example, the values of R1 and R2 are high at Croyde consistent with the more limited amount of interbreeding evident in the triangle plots (Figure 3).

| Population and size variation in genotypic categories
A summary of the population frequencies of the G2 genotypic categories indicates marked differences between populations (Figure 4).
The underlying data and individual-based histograms (File S3-wsC and wsD(A)) confirm this variation.The frequencies of G2 categories were compared between the four hybrid populations using three Population x Category contingency tests.The tests were carried out separately for all mussels, and for small-and large-size classes, and in a hierarchical manner allowing comparison between different categories or groups of categories.In eight of the nine tests, the overall contingency tables are significant indicating differences between hybrid populations (File S3-wsH).In accord with Figures 3 and 4, P0 and P1 differ in frequency between the hybrid populations, and there is a significant excess of F1 at Croyde and deficit of BC1 at Croyde and Carlyon.
There are significant differences in the genotypic category frequency distributions between large-and small-size classes for the individual hybrid populations (File S3-wsI).In accord with Figure 4, P0 and P1 show a significant excess in small and large mussels, respectively, at Croyde, Carlyon, and Whitsand, and BC1 has a significantly higher frequency in large mussels at Carlyon and in all mussels pooled.Figures 3 and 4 both suggest a higher frequency of BC1 than BC0 individuals.Over all mussels, the ratio of BC1 to BC0 is 54 to 16.This is significant in a binomial test against a 1:1 ratio (p = .000).

| SNP loci showing atypical variation
The BayeScan analysis of all loci for the six populations within the hybrid zone reveals no outlying loci for P0 individuals.For P1 individuals, there are two outlier loci, SNP18 and SNP32 with high Eighteen of the mussels are classified as P0 by NewHybrids confirming that this is an M. edulis region.However, one P1 and one BC0 mussel were found.This contrasts with previous reports that found only M. edulis in this region (see Gilg & Hilbish, 2003a).Mantle colour, gonad colour, and gonad development stage show significant differences between populations (File S1 and S5).There is no variation in sex ratio between populations.In addition, a binomial test against a 1:1 expectation is not significant over all populations pooled (p = .922).There is a small but significant association between genotypic category and sex for the reference populations (Figure 6) with females more frequent in P0 and males in P1 (File S5-wsB).

| Phenotypic characters
There is a strong and significant association of genotypic category with mantle colour as expected (Figure 6a), and also a strong association between gonad colour and sex (panel b).Gonad development stage shows less strong but still significant associations with genotypic category and with sex (panels a and b).For mantle colour (panel c), a significant excess of purple is observed for BC1 as well as P1.Of note is the high and often significant excess of brown mantle colour in the F1 and admixed genotypic frequency categories compared with P0 and P1.For gonad colour (panel d), there is a significant excess of orange and white in females and males, respectively, in the hybrid populations, and also in the reference populations.For the latter, association between genotypic category and gonad development is towards a significant excess F I G U R E 4 NewHybrids genotypic category assignations (G2) summed over mussels within populations for all mussels and separately for small and large mussels. of 'ripe' and 'not ripe' gonads in P0 and P1, respectively (panel e), and a significant excess of 'ripe' and 'not ripe' in male and female, respectively (panel f).
A higher frequency of 'not ripe' mantles is observed in the M. galloprovincialis at Padstow and Bude compared with Swansea where 'ripe' is dominant (File S5-wsA).The hybrid populations do not show well-defined spawning periods as 'ripe' is significantly more frequent at Porthcurno and 'not ripe' more frequent at Whitsand.There is also a sex difference with a significant excess of 'ripe' and 'not ripe' in male and female, respectively (Figure 6f and File S5-wsB).When values for 'ripe' (as a %) are plotted as a line graph from west to east, separately for the north and south coast populations in File S5-wsE, there is an interaction such that 'ripe' increases from west to east on the north coast but decreases from west to east on the south coast.In regression analysis, the north-south difference is significant (p = .034)as is the interaction (p = .017).It is surprising that the frequency of F2 is so high (9/239 or 3.8%) given the low frequency of F1 individuals in hybrid populations in SW England (3.8%) and that f2 individuals must be produced from f1 × f1 crosses.A spawning time for f1 intermediate between the parental species would favour a high frequency of such crosses.

| Confirming earlier hybridization findings and UK Mytilus galloprovincialis status
The results of the present study based on 57 SNP loci are consistent with earlier allozyme work (e.g.Skibinski et al., 1983) suggesting hybridisation between Mytilus edulis and Mytilus galloprovincialis in SW England.However, they extend knowledge of this hybrid zone.
The clear differentiation of M. galloprovincialis into genetically distinct types occurring in the Atlantic and Mediterranean was not well understood at the time of the allozyme studies.Failure to recognize this, particularly in wide geographic surveys, creates the possibility of incorrect species identification at some localities (Kotwicki et al., 2021;Simon et al., 2020;Wenne et al., 2020).We show that the UK M. galloprovincialis studied here are closely similar to the Atlantic rather than Mediterranean type.While not excluding the possibility of some introgression in SW England of genes from Mediterranean M. galloprovincialis, the interbreeding and admixture in SW England should focus on that between M. edulis and Atlantic M. galloprovincialis.

| Genetic variation and admixture patterns in SW England
In SW England, the parental species and hybrids are distributed in patches characteristic of a mosaic rather than clinal hybrid zone.The significant differences in the frequency of the G2 genotypic categories between the hybrid populations establish that they are not simply replicate samples from a homogeneous extended hybrid zone.
Studies in SW England with Glu-5 provided evidence of a localized hybrid region from Start Point to St Ives on the south coast (Gilg & Hilbish, 2000, 2003a, 2003b;Hilbish et al., 2002).Parental M. edulis occurs to the east of Start Point and M. galloprovincialis to the northeast of St Ives.Dispersal from the hybrid region into the parental populations was suggested, but there was little evidence of dispersal into the region.Self-recruitment, which may be important in mussel populations (Bierne, Bonhomme, et al., 2003;Gilg et al., 2007;Gosling et al., 2008), is likely to be important in this region.A comparison of allele frequencies at the SNP18 locus in admixed individuals from northern and southern populations supports the hypothesis of local recruitment.The idea that dispersal is occurring northeast of St Ives out of the southern hybrid region (Gilg & Hilbish, 2003a, 2003b) is not supported.The difference in allele frequencies at SNP18 between the north and south coasts suggests two distinct hybrid regions with different recruitment patterns.

The pattern of rising and falling of the G allele frequency from
Start Point to Croyde is similar to that for the hybrid zone between M. edulis and M. galloprovincialis on the Atlantic coast of France from the Bay of Biscay to Normandy (Bierne, Borsa, et al., 2003).This is characteristic of mosaic hybrid zones in Mytilus.Assortative mating, the adaptation of parental species to patchy habitats, and lower fitness of hybrids may be important in the development and maintenance of mosaic hybrid zones (Daguin et al., 2001;M'Gonigle & FitzJohn, 2010).Differences in habitat preference might contribute to reproductive barriers in Mytilus species (Bierne, Bonhomme, et al., 2003;Bierne, David, Langlade, & Bonhomme, 2002;Gosling & Mcgrath, 1990) and could lead to assortative fertilization (Bierne, David, Boudry, & Bonhomme, 2002) as well as long-term maintenance of both species in separate  (Edwards & Skibinski, 1987).This would be consistent with gene flow from Mytilus edulis to M. galloprovincialis and the occurrence of M. edulis in SW England and the absence of M. galloprovincialis in South Wales.Asymmetric introgression of E alleles into M. galloprovincialis has also been observed in French hybrid populations (Fraïsse et al., 2014;Simon et al., 2021).This might reflect an adaptive advantage of M. edulis alleles in other species (Fraïsse et al., 2016;Simon et al., 2021).
The two-locus correlation approach using R1 and R2 statistics contributed to the first genetic evidence of widespread admixture between M. edulis and M. galloprovincialis which varied between localities (Skibinski et al., 1983;Skibinski & Beardmore, 1979).
However, the small number of allozyme loci did not allow individuals to be assigned to different genealogical classes.By contrast the availability of a much larger number of SNP loci, together with more recent dedicated software such as ENTROPY, NewHybrids, and STRUCTURE, allows resolution of a variety of admixed genotypes.There might be confidence that the NewHybrids G2 analysis NewHybrids G2 analysis are reassigned here as genotypic categories unique to the G3 analysis such as P1_BC1.Whether these would belong to the genealogical class p1_bc1 is less certain though our simulations suggest they might be.
The STRUCTURE analysis of generations since admixture and analysis of R1 and R2, both employing mapped SNP loci, suggests that most of the highly admixed individuals might be from the first few generations of interbreeding.The simulation results also support this idea given that M. edulis and M. galloprovincialis are maintained at high frequency within the hybrid zone.This contrasts with dock populations where complete mixing of gene pools can occur and estimates of the number of generations since admixture are higher ranging from 4 to 14 (Simon et al., 2020).Measurement of generation since admixture and the genealogical classes of admixed individuals are relevant to assessment of the strength of the barrier to long-term introgression.At Croyde, F1 individuals predominate, and admixture is low, contrasting with the results for the other hybrid populations.This suggests a stronger local barrier to introgression at Croyde.This has not prevented the high frequency of the SNP18 E allele at Croyde, although such individuals may originally have been produced elsewhere.Local variation in introgression has been observed in other Mytilus species, for example, higher introgression rates in European compared to American M. edulis × M. trossulus hybrid zones (Fraïsse et al., 2016;Simon et al., 2021) and higher introgression of M. trossulus into M. galloprovincialis in America (Saarman & Pogson, 2015).In the Baltic, hybridisation between M. edulis and M. trossulus is so extensive that a hybrid swarm results (Stuckas et al., 2017).

| Factors affecting genetic and phenotypic trait variation
Explanation of the variation in frequency of genotypic categories and extent of admixture between populations, the mosaic structure, and asymmetric admixture can be sought in a variety of fitness-related factors.There is evidence that in SW England M. galloprovincialis has higher selective viability (Gardner et al., 1993;Gardner & Skibinski, 1988, 1991a;Skibinski, 1983;Skibinski & Roderick, 1991;Wilhelm & Hilbish, 1998), fecundity (Gardner & Skibinski, 1990), growth rate (Gardner et al., 1993;Seed, 1971), and feeding and metabolic rate (Hilbish et al., 1994).A low frequency of first-generation hybrids (f1) in a population such as Croyde suggests the involvement of prezygotic and or postzygotic reproductive barriers (Bierne et al., 2006;Brannock & Hilbish, 2010;Diz & Skibinski, 2007;Romero et al., 2019;Wilhelm & Hilbish, 1998).Since admixed individuals are frequent overall in the hybrid populations, low viability of f1 seems unlikely.This is consistent with the studies of age-dependent variation in the zone which suggest that hybrids are not inferior in viability to the parents (Gardner & Skibinski, 1988, 1991a).A low rate of generation of f1 would be favoured by differences in spawning time between M. edulis and M. galloprovincialis.
Through analysis of gonad index or histological tests, overlapping reproductive cycle differences have been reported between sites and seasons with M. edulis spawning earlier in the year (Bierne, Bonhomme, & David, 2003;Bignell et al., 2008;Gardner, 1992;Gilg & Hilbish, 2003c;Secor et al., 2001;Seed, 1971).The observation here of a greater frequency of 'ripe' mantles at Swansea in March compared with Padstow and Bude is consistent with this.
The observed sex difference in spawning time is consistent with a previous report that M. galloprovincialis males spawn earlier than females, although synchronous spawning occurred in M. edulis (Secor et al., 2001).
The ability to allocate individual mussels to NewHybrids genotypic categories allows assessment of phenotypic and adaptive differences between these categories.Mantle edge colour is an important distinguishing feature between the two species (Lewis & Seed, 1969).Darker mantles in M. galloprovincialis may be an evolutionary response to more intense sunlight in the Mediterranean (Seed, 1971).The excess of brown mantle in the F1 suggests an allelic or non-allelic interaction associated with high heterozygosity.
The brown colour may have some as yet unidentified adaptive significance for hybrid individuals.The results confirm previous observations that female and male gonads tend to be pinky-orange and creamy-white, respectively (e.g.Seed, 1969), though more recent work suggested that gonad colour is not a reliable trait to determine sex (Mikhailov et al., 1995;Petes et al., 2008).

| Implications of the results for selective breeding in aquaculture
Because consumers might have a preference for mussels with particular gonad colour, for example, orange, which is due to carotenoids, the results may have important applications in aquaculture.Doubly uniparental inheritance of mitochondrial DNA (DUI) in Mytilus (Fisher & Skibinski, 1990;Skibinski et al., 1994aSkibinski et al., , 1994b;;Zouros et al., 1994aZouros et al., , 1994b) ) is associated with biasing of sex ratios in Mytilus (Kenchington et al., 2002;Saavedra et al., 1997).This introduces the possibility that mussels with particular gonad colour might be generated in hatcheries with sex-biased crosses.A heritability of 0.9 for this character in M. galloprovincialis (Pino-Querido et al., 2015) suggests possible modification by conventional selective breeding or strain selection.The sex biasing in individual crosses seems to average out in population samples to give a 1:1 sex ratio (Fisher & Skibinski, 1990;Seed, 1969).Thus, the cause of the association between genotypic category and sex for the reference populations reported here is unclear.
Techniques for interbreeding mussels in the laboratory are improving (Kenchington et al., 2020).Thus, if a high-density linkage map can be obtained for Mytilus in future, better information on genealogical classes and the ancestral relationship between individual mussels should be accessible as in human genetics (Henn et al., 2012).Such a linkage map would allow great progress towards determining the genetic basis and adaptive significance of variation of mantle colour, gonad colour, and other phenotypic attributes through genome-wide association studies (Luo et al., 2023;Peñaloza et al., 2022;Uffelmann et al., 2021).

| Evidence of introgression from M. edulis into M. galloprovincialis at SNP18
The SNP loci studied here as well as Glu-5 show a similar pattern of geographic variation.However, SNP18 shows a strikingly different pattern.Outlier loci, which may be caused by selection (Lotterhos & Whitlock, 2014;Morin et al., 2004), were observed in previous studies on Mytilus in France and Spain (Bierne, David, Langlade, & Bonhomme, 2002;Fraïsse et al., 2014;Gosset & Bierne, 2013).In these studies, as with SNP18, an M. edulis allele had an unexpectedly high frequency in M. galloprovincialis.The interpretation of outlier loci and their causes is complex (Fraïsse et al., 2016).However, the possibility that introgression and selection of a specific region of the M. edulis Notch gene are adapting M. galloprovincialis to colder northern temperatures could be pursued.Future work could derive three-dimensional structures for the Notch protein to try to infer function as has been used for a mtDNA protein (Skibinski et al., 2017).The function of the Notch protein in ectothermic organisms might be influenced by temperature (Shimizu et al., 2014).Thus, different alleles could have adaptive significance by promoting optimal development or growth at different seawater temperatures.Evidence in support of the hypothesis that the allele present in M. galloprovincialis could have evolved after its evolutionary split from M. edulis is supported by the observation that mussel samples from the conspecific M. trossulus are fixed for the same allele as in M. edulis (Simon et al., 2021), as well as by BlastX analysis (nrNCBI) indicating that Notch sequences from the more distantly related M. californianus and M. coruscus also have the same allele as in M. edulis and M. trossulus.

| Temporal stability of SW England hybrid zone and climate change implications
Mytilus edulis and M. galloprovincialis mussels have been recognized as such in SW England over many decades of research, dating back to the beginning of the nineteenth century (see Hepper, 1957).In mussels collected in SW England from 1980-1981to 1996-1998, the frequency of M. galloprovincialis was greater in larger and older than smaller and younger mussels (Gardner et al., 1993;Gardner & Skibinski, 1988;Hilbish et al., 2002;Skibinski, 1983;Wilhelm & Hilbish, 1998).This phenomenon was also reported in another mosaic hybrid zone region in NE England (Skibinski & Roderick, 1991).
The relationship between shell length and M. galloprovincialis allele frequencies has been updated to 2016 incorporating the allozyme Est-D, Glu-5, and the SNP loci for samples from Whitsand, Croyde, and Porthcurno (Figure 7).The positive relationship between shell length and the G allele frequencies is generally similar in these different studies over a period of 35 years showing little evidence of substantial shifts or oscillations during the period.Three hypotheses were suggested to explain the apparent temporal stability, that is, differential growth, differential mortality or viability, and a historical change involving replacement of M. edulis by M. galloprovincialis within the hybrid zone (Gardner & Skibinski, 1988;Skibinski, 1983).Higher mortality of M. edulis and viability of M. galloprovincialis is likely the preferred hypothesis (Gardner et al., 1993;Gardner & Skibinski, 1988;Wilhelm & Hilbish, 1998).Higher strength of attachment to the substrate or higher resistance to thermal stress favouring M. galloprovincialis were proposed as selective factors (Gardner & Skibinski, 1988;Hilbish et al., 1994Hilbish et al., , 2002;;Popovic & Riginos, 2020;Willis & Skibinski, 1992).The stability of the hybrid zone might be due to the compensating immigration of M. edulis (Edwards & Skibinski, 1987;Gardner & Skibinski, 1988;Skibinski & Roderick, 1991;Wilhelm & Hilbish, 1998).Sea surface temperature (SST) in the UK showed a small gradual increase of 0.2-0.3°Cfrom 1870 to around 1980, but from 1980 to 2020 showed a more dramatic increase of 0.7-0.8°C(Kendon et al., 2021).Thus, these hybrid populations do not seem to have responded to global warming by a historical change shifting towards higher G allele frequencies.There may, however, be highly localized effects.For example, although early work reported typical M. edulis at Padstow (e.g.Seed, 1971;Skibinski, Ahmad, & Beardmore, 1978), such individuals were not identified during collection for the present study.
A distinction can be drawn between the movement of a hybrid zone through asymmetric introgression and the displacement of one species by another through invasion.Mytilus, particularly M. galloprovincialis, is renowned as an invasive species and has displaced native mussels at many localities worldwide (Branch & Steffani, 2004;Pickett & David, 2018).For example, M. galloprovincialis is steadily displacing native M. trossulus in California because of ecological factors (Lockwood & Somero, 2011;Saarman & Pogson, 2015).An increase of M. galloprovincialis on the Irish coast was also reported during a time interval of 25 years (Gosling et al., 2008).High levels of genetic variation in Mytilus may contribute to its invasive success (Han & Dong, 2020).The northern expansion of M. galloprovincialis might be limited by cold winter temperatures affecting reproductive development and the preference of M. galloprovincialis for warmer temperatures (Fly et al., 2015;Fly & Hilbish, 2013).In studies of three hybrid zones of M. edulis and M. galloprovincialis in France over a period of two decades, one of the zones changed its position into warmer water (Hilbish et al., 2012).
Many factors might influence the stability and movement of clinal and mosaic hybrid zones including the fitness and adaptation of parents and hybrids, variation in dispersal, changes in population density, assortative mating, and asymmetric introgression (Barton, 1979;Buggs, 2007;M'Gonigle & FitzJohn, 2010;Wielstra, 2019).We have cited studies that suggest that in SW England M. galloprovincialis is superior to M. edulis in several fitness-related variables.However, there are also potentially balancing factors such as asymmetric admixture of M. edulis genes into M. galloprovincialis and immigration of juvenile M. edulis into the zone.Despite such factors that might predispose the zone to move, there is no clear evidence that this is happening in response to global warming either through invasion or introgression.

| CON CLUDING REMARK S
The growing literature on hybrid zones suggests an increasing appreciation of their general importance for the study of evolution (Gompert et al., 2017;Johannesson et al., 2020).This includes the realization that hybrid zones, in both animals and plants, are perhaps more frequent than was once thought (Taylor & Larson, 2019;Zbinden et al., 2023).We pooled results for two reviews reporting hybrid zones potentially influenced by climate change (Chunco, 2014;Taylor et al., 2015).A total of 40 examples were for animals and five for plants.There were only two marine examples (mussels and sole fish).While recognizing that there may be many ascertainment biases in drawing conclusions about natural hybrid zones, the results do suggest that there might be a relatively low level of activity of marine research in this area.Demonstrating a clear link between climate and hybrid zone movement is difficult because long-term studies are needed involving the monitoring of genetic, phenotypic, environmental, and ecological variables.
There is growing understanding of the importance of hybridisation in generating novel genetic adaptations, and that the porous nature of many hybrid zones allows introgression of genes which might undergo positive selection in recipient species (Harrison & Larson, 2014).Advanced in genomics and genome sequencing has greatly advanced such knowledge (Moran et al., 2021).For example, several percent of the genome of some modern humans derives through introgression from close relatives such as the Neanderthals and appears to have potential adaptive significance (Dannemann et al., 2017;Dannemann & Kelso, 2017).Such approaches are developing in Mytilus augmented by large-scale studies of the correlation of genetic and environmental variables (Kijewski et al., 2019;Wenne et al., 2020).Given that global warming does increase, invasive species such as Mytilus may become increasingly important in monitor- termined with this test and recorded as Female or Male.Mantle edge colour, gonad colour, and the shell length and height of each mussel were also recorded.Mantle colour was assigned as Purple, Brown, or White/Straw (abbreviated 'White'), and gonad colour as deep orange (abbreviated 'Orange'), intermediate, or Pale Orange/White (abbreviated 'White').
For genotyping, PCR products together with the NZYDNA ladder V (NZYtech, Lisbon, Portugal) were separated in 2% agarose gels supplemented with 2.5 μL of SYBR Safe dye at 90 volts for 40 min.Individuals were genotyped for two alleles: 126 (G) and 180 (E) base pairs long.The G or E alleles are at high frequency (>95%) in the reference populations M. galloprovincialis and M. edulis, respectively, used in this study.SNP genotyping was carried out by LGC genomics (Hoddesdon, UK) following KASP™ (Kompetitive Allele Specific PCR; Semagn et al., 2014) for 57 SNP markers which show high F ST values between allopatric reference populations from M. edulis and M. galloprovincialis
R2 measures the relative frequency of backcross (EE EG, EG EE, GG EG, and EG GG) compared with F1 (EG EG) and other genotypes.The methods of computation of R1 and R2 are shown in File S3-wsK.R1 and R2 were calculated for individual populations for pairs of SNPs on the genetic linkage map of Simon et al. (2020) for a range of map distances.

F
I G U R E 2 Genotype proportions (EE, EG, and GG), G allele proportion, and F IS averaged over 60 individuals and 57 SNPs.Values are given for all mussels and separately for small (20-29 mm) and large (30-40 mm) mussels.The corresponding numerical values are given in File S2-wsP.The underlying raw data are given in File S2-wsF to wsM which also give calculated p-values for F IS .Population abbreviations are defined in the text. in smaller mussels at Whitsand.These contrasts between large and small are apparent in the pattern of F IS values for the individual SNP loci (File S2-wsN).
, and all samples pooled in File S3-wsJ).The plotted points are in different regions of the triangle as expected.In reference or mixed populations of M. edulis and M. galloprovincialis without interbreeding, individuals are clustered at low Q 12 to the bottom left (e.g.Vigo, Bude, and Padstow) and right (e.g.Swansea) corners of triangle.As interbreeding proceeds, F1 individuals appear at the apex (e.g.Croyde on the north coast).Porthcurno, Carlyon, and Whitsand on the south coast show a greater level of interbreeding and admixture than Croyde.Though the number of F1 individuals identified is small in the four hybrid populations (7/239 or 2.9%) (File S3-wsE), there is a significant excess at Croyde relative to the others (File S3-wsH).Croyde has no other admixed individuals identified by NewHybrids (Figure 3, File S3-wsF(A)).Individuals in the Swansea and Vigo reference populations have q-values within a narrow range.By contrast, the M. galloprovincialis populations Bude and Padstow have a few individuals assigned to the F1 and BC1 categories suggesting some interbreeding and admixture.Some of the individuals assigned to P1

F
ST (File S3-wsL(A)).Mahalanobis distance revealed SNP18 alone as being a significant outlier (File S3-wsL(B)).Plots of G allele proportions for the P1 individuals for the six populations within the hybrid zone region demonstrate the contrast between SNP18 and the other SNP loci (File S3-wsL(C)).The E allele proportion is high in P0 individuals in both the north and south coasts of SW England (Figure 5).However, in P1 individuals, the G allele though having the expected high frequency in the southern populations Porthcurno, Carlyon, and Whitsand occurs at a low frequency in the northern populations Padstow, Bude, and Croyde.The high frequency of the SNP18 E allele in P1 individuals in northern populations (Figure 5) allows two tests for dispersal of admixed individuals out of the hybrid zone region that stretches from Start point to St Ives discussed by Gilg and Hilbish (2003a) (Figure 1).These tests are set out in detail in File S1.The first test shows that seven hybrid individuals at Padstow, Bude, and Croyde have SNP18 E allele frequencies characteristic of M. galloprovincialis from the north but significantly different from M. galloprovincialis from the southern hybrid zone region.The second test shows that individuals classified as the progeny of P1 and BC1 (P1_BC1 in the G3 NewHybrids analysis) have a high frequency of the SNP18 E allele in line with its overall high frequency in M. galloprovincialis in the north (File S1).Recently, we have genotyped 20 mussels collected in 2018 at Dartmouth east of Start Point for the same SNPs (File S2-wsO).
The SNP18 locus is a good candidate to investigate possible selection acting on it.BlastX with the sequence containing the SNP18 polymorphism gave best hits for the Notch signalling transmembrane protein from Mytilus.The polymorphism occurs in a codon where a cytosine in M. edulis is substituted by an adenine in M. galloprovincialis associated with a non-synonymous change from proline to glutamine.The Gene Ontology term associated with Notch (GO: 0007219) links the 'Notch signaling pathway' to 'a series of molecular signals initiated by the binding of an extracellular ligand to the receptor Notch on the surface of a target cell, and ending with regulation of a downstream cellular process, e.g.transcription'.
The classification of individual mussels for sex, mantle edge, and gonad colour, and gonad development stage are recorded in File S2-wsE.Photos of mussels differing in gonad and mantle colour are shown in File S2-wsD.

F
I G U R E 5 Map showing values for the proportion of E allele and G allele in P0 and P1 individuals in the eight populations for SNP18 Notch gene.patches important for maintenance of mosaic patterns of variation.A theoretical model suggested that mosaic zones can be stable despite long-distance dispersal and an application of the model to the hybrid zone in France confirmed high scores of mosaicity (M'Gonigle & FitzJohn, 2010).The excess of BC1 over BC0 in the triangle plots suggesting asymmetric admixture has a parallel in earlier studies with mitochondrial DNA.Haplotypes at high frequency in M. galloprovincialis from Bude and Padstow were absent from Swansea though haplotypes at high frequency in Swansea were present in Bude and Padstow picks out individuals in the specific genealogical classes, for example, that F1 individuals are f1.Many individuals assigned as P1 in the F I G U R E 6 Results of tests of association for pooled reference and pooled hybrid populations.Cramer's V for association between NewHybrids genotypic categories and phenotypic characters (a) and between Sex and the other phenotypic characters (b).(c-f) Distributions for characters as shown in the figure.Significance levels are obtained from adjusted standardized residuals for individual cells in the contingency tables.*p < =.05; ** < =.01; *** < =.001.Further details are given in File S5-wsB, wsC, and wsD.
ing work.Over four decades ago when our genetic work on Mytilus in SW England began we could not have foreseen these developments and their importance.AUTH O R CO NTR I B UTI O N S A.P.D. Conceptualization, Methodology, Investigation, Formal analysis, Writing -original draft, review and editing, Funding acquisition.

F
Historic evolution of G allele frequencies obtained from highly ancestry-informative markers in different studies where mussel samples analysed were collected from 1980 to 2016 in the following hybrid populations in SW England: Croyde (CR), Porthcurno (PC), and Whitsand Bay (WB).The molecular markers used (in parenthesis) and studies from which the G allele frequencies were derived as follows: Skibinski (1983) (Est-D allozyme marker; sample collection in 1980-1981), Gardner & Skibinski (1988) (Est-D; sample collection in 1986-1987), Wilhelm & Hilbish (1998) (Est-D; sample collection in 1990-1991), Hilbish et al. (2002) (Glu-5 nuclear marker; sample collection in 1996-1998), and results from the present study deriving from the genotyping of Glu-5 and 57 SNPs (average G allele frequency) in mussel samples collected in 2016.