# Discordant population structure inferred from male- and female-type mtDNAs from Macoma balthica, a bivalve species characterized by doubly uniparental inheritance of mitochondria

Sabrina Le Cam, Julie Brémaud, Vanessa Becquet, Valérie Huet, Emmanuel

Dubillot, Pascale Garcia, Amélia Viricel, Sophie Breton, Eric Pante

► To cite this version:

Sabrina Le Cam, Julie Brémaud, Vanessa Becquet, Valérie Huet, Emmanuel Dubillot, et al. Discordant population structure inferred from male- and female-type mtDNAs from Macoma balthica, a bivalve species characterized by doubly uniparental inheritance of mitochondria. Peer Community Journal, 2025, 5, pp.e29. 10.24072/pcjournal.529. hal-05465467

HAL Id: hal-05465467

https://hal.science/hal-05465467v1

Submitted on 19 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Research article

Published
2025-03-11

**Cite as**
Sabrina Le Cam, Julie Brémaud, Vanessa Becquet, Valérie Huet, Emmanuel Dubillot, Pascale Garcia, Amélia Viricel, Sophie Breton and Eric Pante (2025) *Discordant population structure inferred from male- and female-type mtDNAs from Macoma balthica, a bivalve species characterized by doubly uniparental inheritance of mitochondria*, Peer Community Journal, 5: e29.

**Correspondence**
sablecam15@gmail.com
pante@univ-brest.fr

**Peer-review**
Peer reviewed and
recommended by
PCI Evolutionary Biology,
https://doi.org/10.24072/pci.evolbiol.100540

<figure><img src="image_2.png" /></figure>

This article is licensed
under the Creative Commons
Attribution 4.0 License.

# Discordant population structure inferred from male- and female-type mtDNAs from *Macoma balthica*, a bivalve species characterized by doubly uniparental inheritance of mitochondria

Sabrina Le Cam<sup>1,2</sup>, Julie Brémaud<sup>1,3</sup>, Vanessa Becquet<sup>1,1</sup>, Valérie Huet<sup>1,1</sup>, Emmanuel Dubillot<sup>1</sup>, Pascale Garcia<sup>1,1</sup>, Amélia Viricel<sup>1,4</sup>, Sophie Breton<sup>1,3</sup>, and Eric Pante<sup>1,4</sup>

Volume 5 (2025), article e29

https://doi.org/10.24072/pcjournal.529

## Abstract

Doubly Uniparental Inheritance (DUI) of mitochondria is a remarkable exception to the Strictly Material Inheritance (SMI) in metazoans. In species characterized by DUI – almost exclusively gonochoric bivalve mollusks – females (F) transmit mitochondria to offspring of both sexes, while males (M) pass on their mitochondria exclusively to their sons. Under DUI, males are heteroplasmic, somatic tissues containing F-transmitted mtDNA and gametic cells containing M-transmitted mtDNAs. The aforementioned transmission routes make M- and F-transmitted mtDNA interesting as sex-specific markers which can differ in their effective population sizes, mutation rates, and selective constraints. For these reasons, looking at both markers can provide significant insights into the genetic structure of populations and investigate its determinants. In this study, we document differences in genetic diversity, divergence, inter-populational differentiation, and biogeographic structure between M- and F-type *cox1* mt genes in the Baltic tellin (*Macoma balthica*) to test whether *cox1m* and *cox1f* genes bear the marks of similar phylogeographic histories. These markers were sequenced for 302 male individuals sampled from the North Sea to the Gironde Estuary (Southern France) encompassing the intra-subspecific *M. b. rubra* hybrid zone in the Gulf of Saint-Malo. Both genes supported a scenario of cladogenesis of *M.b. rubra* clades prior to the last glacial maximum. Nucleotide diversity and net divergence were over twice higher in *cox1m* compared to *cox1f*. Genetic differentiation between northern and southern populations was nearly 3 times higher at *cox1m* compared to *cox1f* (global *F*ST = 0.345 and 0.126 respectively) and the geographic localization of the strongest genetic break significantly differed between the markers (Finistère Peninsula at *cox1f*; Cotentin Peninsula at *cox1m*, ~250 km apart). A higher mutation rate, relaxed negative selection, and differences in effective population sizes (depending on locations) at *cox1m* could explain differences in population genetic structure. As both F- and M-type mtDNAs interact with nuclear genes for oxidative phosphorylation and ATP production, geographical discordances in genetic clines in a context of secondary contact could be linked to mito-nuclear genetic incompatibilities.

<footnote n="1">¹Laboratoire Littoral Environnement et Sociétés (LIENSs), UMR 7266 CNRS-LRU, La Rochelle, France, ²Département de biologie, Université du Québec à Rimouski, Rimouski, QC, Canada, ³Département de sciences biologiques, Université de Montréal, Montréal, QC, Canada, ⁴Univ Brest, CNRS, IRD, Ifremer, UMR 6539, LEMAR, Plouzané, France</footnote>
Introduction

Some species show a remarkable exception to the maternal inheritance of mitochondria in metazoans: the doubly uniparental mode of inheritance (DUI). In this system, both males and females can transmit their mitochondria. The former transmits female-inherited F-type mitochondria to all their progeny and the latter pass on "male-inherited" (M-type) mitochondria to their male offspring, where the male mitogenomes (mt) are quartered in the germ line and gametes (reviewed in Zouros, 2013). To date, DUI has only been discovered in the class Bivalvia, with over 100 DUI species (Gusman et al., 2016) among the c.a. 11,000 (Huber, 2010; Huber et al., 2015). They are all gonochoric, except for the hermaphroditic mussel *Semimytilus algosus* (Lubośny et al., 2020). More than a simple peculiarity, DUI is suspected to play a role in sex-determination and gonad differentiation (Zouros, 2000; Breton et al., 2011; Guerra et al., 2017; Capt et al., 2018, 2019) and could well be involved in population structure through intrinsic (e.g. genetic incompatibilities; Saavedra et al., 1996) and extrinsic (e.g. selection and demography, Stewart et al., 1996) factors.

In DUI species, divergence between F-type and M-type mitogenomes is variable but generally high, ranging from 6 to over 50% (reviewed in Breton et al., 2007; Gusman et al., 2016), which questions the maintenance of mito-nuclear genetic coadaptation. Indeed, both F- and M-type mitochondria can be found in males and females but in majority, females are homoplasmic for the F-type mtDNA whereas males are heteroplasmic, accommodating two highly divergent mitogenomes (F-type in somatic tissues and M-type in sperm). The presence of the M-type mtDNA in somatic tissues is considered as a paternal leak due to elimination or segregation failure of sperm mitochondria in female or male embryos, respectively (Milani et al., 2012). Both F- and M-type mt lineages show rapid molecular evolution compared to other animals, the M-type mtDNA usually evolving faster than the F-type mtDNA (Zouros, 2013). Coevolution and coadaptation of mitochondrial and nuclear genes are required for efficient cellular energy production (i.e. oxidative phosphorylation OXPHOS) and mito-nuclear genetic incompatibilities (MNIs) can lead to a desynchronization of this machinery (Burton & Barreto, 2012; Burton et al., 2013). DUI could, therefore, bear on the maintenance of genetic structure among populations of highly dispersive bivalve species at small spatial scales and provide key insight into the establishment and maintenance of local adaptation.

Indeed, barriers to gene flow can arise and be maintained by a multitude of environmental and/or intrinsic factors (Barberousse et al., 2010), from ecological isolation to genetic incompatibilities. Hybrid zones, which correspond to regions spatially separating genetic stocks, are "natural laboratories" to study the interactions between intrinsic barriers and the environment, and the processes of adaptation and speciation.

*Macoma balthica*, a species in which DUI has recently been detected (Pante et al., 2017; Capt et al., 2020), is a noteworthy model species to study hybrid zones in marine ecosystems (Riginos & Cunningham, 2007; Strelkov et al., 2007). It has a wide distribution range spanning from the west pacific coasts, in Japan and from Alaska to Oregon (USA, Luttikhuizen et al., 2003) to the North Atlantic, where the species is found in the west from Arctic to Newfoundland (Layton et al., 2016) and in the east from the north of Russia (Hummel et al., 1997) to the Arcachon Basin (Hily, 2013 and this publication). The succession of glaciation and inter-glaciation periods has resulted in colonization of the Atlantic marked by repeated isolation, invasion, and re-colonization events (Nikula et al., 2007). These episodic colonization events have created multiple opportunities for secondary contacts between different genetic stocks and the establishment of several hybrid zones in the North Atlantic, where two subspecies of *M. balthica* co-occur. A Pacific lineage (*M. b. balthica*) is present in the Baltic Sea and the White Sea, and an Atlantic lineage (*M. b. rubra*) is found in the Norwegian Sea, the North Sea and along the British coasts, down to the southern range limit of the species (Luttikhuizen et al., 2003; Nikula et al., 2008; Väinölä, 2003). In Europe, two genetic breaks were described: an inter-subspecific break occurs between *M. b. balthica* and *M. b. rubra* in the transition zone from the North Sea to the Baltic Sea (Luttikhuizen et al., 2012;
Nikula et al., 2007, 2008), and an intra-subspecific break within *M. b. rubra* in the Gulf of Saint-Malo (Brittany, France). In the latter case, southern populations of *M. b. rubra* (distributed between the southern range limit and the aforementioned break) exhibit F-type mtDNA signatures consistent with long-term isolation in the glacial refugium of the Bay of Biscay: high genetic diversity relative to previously glaciated areas, high prevalence of private alleles, and sharp genetic differentiation from northern populations (Becquet et al., 2012). Deciphering signal of primary intergradation from secondary contact in hybrid zones is challenging; we have therefore gathered sequences from both F and M mitogenomes to estimate the timing of cladogenesis in the *Macoma* species complex, particularly in the *M. b. rubra* hybrid zone. Clarifying the phylogeographic context will allow us to focus on the effects of DUI on the evolution of *M. b. rubra* lineages. Indeed, multiple nuclear genes involved in the oxidative phosphorylation (OXPHOS) system of mitochondria were detected as significantly differentiated among southern and northern populations (Pante et al., 2012, 2019); DUI could therefore play a role in the maintenance of genetic barriers in *M. b. rubra*.

# Materials and Methods

## Sampling

Individuals were collected from a total of 14 sampling sites ranging from Arcachon (southern range limit of the species, France) to Le Croyo (Somme Bay, northern France) and from Kruiningen (the Netherlands) to Sylt (Germany) (Table S1). Individuals from sampling sites on the French coasts were treated as follows: 70 to 100 adults from 11 mm to 23 mm were randomly collected live at sexual maturity between 4^th^ and 23^rd^ of April 2018 at 10 locations ranging from the Bay of Biscay in southern France to Somme Bay (Table S1). Individuals were then held in aquaria until dissection at the LIENSs laboratory (Littoral Environnement et Sociétés) in La Rochelle, France, with water temperature maintained at 10°C. They were fed with a multispecific microalgal mixture every other day and dead or dying individuals were removed daily. For each individual, the adductor muscle was carefully severed to separate the two valves, without damaging the gonad, and a sample of the mantle was taken. Sex and gonadal maturation stage were then determined with a dissecting microscope (x100 to x400). Two types of tissue samples were collected: a gonadic sample and a somatic sample (mantle). All tissue samples were flash-frozen in liquid nitrogen and then stored at -80°C until DNA extraction. At each step of this protocol, dissecting tools and bench surfaces were thoroughly cleaned in successive baths of 10 % bleach, demineralized water, and 95 % EtOH as to avoid DNA cross-contamination. Using this dataset, we checked that male mitochondria were limited to male gonads (i.e. absent from other male tissues and from females); thus making molecular sexing permissible (Le Cam et al., 2023).

For the other sampling sites, whole individuals had been collected for a previous study (Becquet et al., 2012 and the BIOCOMBE project) and were preserved in 95% EtOH. A 3 mm³ piece of tissue was sampled from the gonad, gDNA was purified and quantified as detailed below, and PCR amplifications were attempted with both the *cox1m* and *cox1f* primers. Individuals for which both gene regions could be amplified were considered as males and included in subsequent analyses.

For phylogenetic inferences, additional samples were sequenced to grasp as much diversity as possible within the *Macoma* species complex, from the NW, NE Atlantic and the NE Pacific (Table S1). They were combined to available data from Nikula et al., (2007), Becquet et al., (2012), Layton et al., (2016) and Yurchenko et al., (2018).

## Molecular Analyses

Total DNA was extracted with the Nucleospin Tissue Kit (Macherey Nagel), following the manufacturer's instructions. In Saint-Vaast-la-Hougue (VAA) and Mont-Saint-Michel (MSM), most specimens were sexually undifferentiated at the time of sampling. Therefore, for these samples and for the samples from all the other "non-French" localities, sex was determined by checking for the amplification of both male and female *cox1* markers in gonad DNA. DNA purity and potential contaminants were checked using a NanoDrop2000.
Specific primers were designed to independently amplify the 5' portion of the *cox1* gene of the M-type and F-type mitogenomes (*cox1m* and *cox1f* respectively). To do that, both mitogenomes published in Capt et al. (2020) were used and they came from *M. b. rubra* specimens sampled in Aytré, France (clade b1b). We also used *cox1* sequences retrieved from transcriptomic data (Pante et al., 2012) from specimens sampled in Aytré, Somme Bay (France, b1a clade) and Gdansk (Poland, b2/b3 clade). Polymerase chain reaction (PCR) amplifications of *cox1f* and *cox1m* mitochondrial DNA regions were performed using the following primers (*cox1m*: cox1m14641F ATAGCTGGCCTGGTGTTTAGG, cox1m15560R TTGGACCCTTTCGAGCCAAG; *cox1f*: cox1f5434F TTAGTGACTTCACACGGTTTGC, cox1f6032R TGGGAAATAATCCCAAACCCG). Amplifications were realized in a total volume of 25 µL, with 0.1 µL of Taq Polymerase 5 U.µL⁻¹ (TaKaRa Ex Taq® Kit MgCl₂ Free Buffer) (TaKaRa Ex Taq, TaKaRa Bio, Shiga, Japan), 2.5 µL PCR Buffer 10X, 1.5 µL MgCl₂ 25 mM, 1 µL dNTP 2.5 mM each, 0.6 µL each primer 10 µM and 17.7 Milli-Q water, and 1 to 20 ng of template DNA. SensoQuest thermal cyclers (Göttingen, Germany) were used to perform the following PCR cycling profiles: (i) for the *cox1f*, 2 min of initial denaturation at 94°C, then 30 cycles consisting of 45 sec at 94°C followed by 30 sec of annealing at 57°C and 40 sec of elongation at 72°C, and a 5 min final step of elongation at 72°C, and (ii) for the *cox1m*, 2 min at 94°C, then 30 cycles of 30 sec at 94°C, 30 sec at 60°C and 55 sec at 72°C, and a final elongation for 5 min at 72°C. PCR success (i.e. amplicon concentration, specificity and absence of amplicons for PCR and extraction negative controls), was tested on 1% agarose gels. PCR products were then sent to Eurofins GATC Biotech GmbH (Konstanz, Germany) to be purified and Sanger sequencing was performed in both directions.

## Sequence Analyses

All sequences of *cox1f* sequences and *cox1m* fragments were quality-controlled and aligned in Geneious Prime 10.2.6 (Kearse et al., 2012). Female clade names follow Nikula et al., (2007). Consensus sequences were trimmed to 479 bp and 676 bp for *cox1f* and *cox1m*, respectively. *cox1m* sequences, trimmed at the same length as *cox1f* (479bp), were used for further comparative analyses. Haplotype sequences and their frequencies were computed in R v3.3.3 (R Core Team, 2017) using the *pegas* package v0.11.(Paradis, 2010). In order to model intraspecific relationships among haplotypes, for both *cox1f* and *cox1m*, median-joining networks (Bandelt et al., 1999) were constructed with PopART v1.7 software (Leigh & Bryant, 2015). For each marker, we constructed a maximum-likelihood phylogenetic tree using IQ-TREE software (Nguyen et al., 2015) with 1,000 bootstrap replicates, and nucleotide acid-substitution models were inferred using ModelFinder with Bayesian information criteria (Nguyen et al., 2015; Kalyaanamoorthy et al., 2017). For *cox1f*, data from Nikula et al., (2007), Becquet et al., (2012), Layton et al., (2016) and Yurchenko et al., (2018) were used for phylogenetic inference. For *cox1m*, we also sequenced data from the *Macoma* species complex, that is *M. petalum* (clade A) and *M. balthica balthica* (clades C and D). From these genealogical relationships, mitochondrial lineages were defined (haplogroups) and their population relative frequency estimated.

We used the net distance approach to estimate the timing of cladogenesis events within the *Macoma* tree. We used the same dataset as described for phylogenetic reconstruction, estimated the most likely model of nucleotide substitution in IQTREE. In MEGA 11 (Tamura et al., 2021), we calculated net distances using the Tamura-Nei model (Tamura & Nei, 1993) as it was the closest model identified by IQTREE. The rate variation among sites was modeled with a gamma distribution (shape parameter = 0.17). The differences in the composition bias among sequences were considered in evolutionary comparisons (Tamura & Kumar, 2002). All codon positions were included, using the pairwise deletion option for missing data. We used 115, 479 bp-long sequences for *cox1f* and 102, 676 bp-long sequences for *cox1m*. The substitution rates for *cox1f* and *cox1m* were then estimated using the net distance described above and the calibration dates proposed in Nikula et al., (2007) for the trans-Atlantic divergence between *rubra* and *balthica*. 3.5My corresponds to the opening of the Bering Strait, the earliest time of invasion of *Macoma* into the Atlantic. 2My corresponds to the earliest known fossils of *M. balthica* in Europe (Norton & Spaink,
1973). This strategy was also used by Wares & Cunningham, (2001) and Riginos & Cunningham, (2007), and Riginos et al., (2004) in the context of DUI.

Genetic diversity indices and population differentiation analyses at *cox1f* and *cox1m* were estimated only for sampling sites with a minimum sample size of 10 individuals. Genetic diversity indices were calculated using Arlequin v3.5.2.2 (Excoffier et al., 1992): haplotype number (H), haplotype diversity (H<sub>d</sub>) and nucleotide diversity (π Tajima, 1983). The differences between the distributions of each genetic diversity statistic (π, H<sub>d</sub>) were statistically evaluated using Wilcoxon rank tests for paired data. Genetic differentiation was estimated in Arlequin using the *Φ*<sub>ST</sub> fixation index and its statistical significance was tested by performing 10,000 permutations (Excoffier et al., 1992). We performed maximum-likelihood geographical cline analyses with the R package HZAR 0.2-9 (Derryberry et al., 2014) together with cline centers and cline widths estimates. Clines were fitted on haplogroup frequencies and geographical distances from the most northern sample sites (SYL). For each haplogroup frequencies, two models were fitted and compared to the null model (no cline): I (fixed scaling based on minimum and maximum observed frequencies and no exponential tails) and II (free scaling and no exponential tails). We performed model selection using AICc scores to determine the preferred model and then extracted the ML model parameters. Geographical distances along the coast were computed using the marmap 1.0.10 R package (Pante & Simon-Bouhet, 2013).

An analysis of molecular variance (AMOVA) was carried out to evaluate the genetic structure among geographic regions and overall sampling sites using 10,000 permutations with Arlequin v3.5.2.2 (Excoffier et al., 1992). We used jModeltest2 (Guindon & Gascuel, 2003; Darriba et al., 2012) to choose the Tamura-Nei model of nucleotide substitution (with gamma shape parameters of 0.1150 and 0.1450 for *cox1f* and *cox1m*, respectively) as a measure of haplotype distances in Arlequin v3.5.2.2. Based on discordances in the distribution of the genetic diversity at *cox1f* and *cox1m*, several hierarchical groupings of subpopulations were tested for each marker to determine the best geographical delineation of genetic structure (Table 2: 3 groups (CRI,KRU / CRO,VAA,MSM,BRI / MAH,BRE,AYT,FOU,ARC); 2 groups A (CRI,KRU,CRO,VAA / MSM,BRI,MAH,BRE,AYT,FOU,ARC); 2 groups B (CRI,KRU,CRO,VAA,MSM,BRI / MAH,BRE,AYT,FOU,ARC)). Different measures of populations differentiation (pairwise and AMOVA) were also carried out using conventional F-statistics.

For the *cox1f* and *cox1m* markers, mismatch distribution analyses were performed to infer
deviations from long-term population stability. We considered two genetically-homogeneous
groups (as in Marko et al., 2010): a “North” group composed of individuals from KRU, SYL, WIL
and CRI, and a “South” group with individuals from AYT, ARC, BRE and FOU. As detailed in the
Results, no genetic sub-structure was detected within these two groups using pairwise *Φ*<sub>ST</sub>
Observed distributions were compared to distributions expected under constant size and
growth/decline models. Raggedness r (Harpending, 1994) and Ramos-Onsins and Rozas’ R<sup>2</sup>
(Ramos-Onsins & Rozas, 2002) statistics were estimated, and their significance along with
confidence intervals were computed using coalescent simulations (10 000 replicates). These
analyses were carried out with DNAsp6 (Rozas et al., 2017). To test for differences in selective
pressure on the female and male mitochondrial sequence sets, Fu’s Fs (Fu, 1997) and Tajima’s D
(Tajima, 1989) test for deviations from the neutral theory model for a population of constant size
was calculated in Arlequin v3.5.2.2 (Excoffier et al., 1992) for each marker and group. Their
statistical significance was tested with 10,000 permutations. Finally, relative differences in effective
population size between both markers were estimated with the θ ratio (Watterson, 1975) among
genetically homogenous groups of samples (as described above for mismatch distributions). θ
indices, based on the number of segregating sites were estimated using the *theta.s* function of the
pegas 1.3 R package (Paradis, 2010).

Results

While this study focuses on the differences in genetic diversity and population structure at the
*cox1f* and *cox1m* markers, from the North Sea to the southern limit of the distribution range, we
first investigated the overall phylogeographic context of these lineages. Phylogenetic inferences from the *cox1f* dataset combined with Nikula et al., (2007), Becquet et al., (2012), Layton et al., (2016) and Yurchenko et al., (2018) data revealed that our English Channel / Bay of Biscay dataset is only composed of the *M. balthica rubra* lineage (Figure S1, locations of the samples are presented in Figure S3). The phylogenetic tree and the haplotype network were congruent: the "northern" haplogroup corresponded to the b2/b3 *M. b. rubra* lineage and our data revealed 2 sub-lineages within the b1 lineage, b1a and b1b. The phylogeny of *cox1m* revealed a similar pattern, with three clades m1, m2, m3 (Figure S2) and a fourth clade m4 encompassing mainly individuals from the Baltic Sea and originating at the base of the *M. b. rubra* branch. Based on net divergence, we were able to compare the timing of cladogenesis of *cox1m* and *cox1f* clades (Table S2). In particular, the divergence of *cox1m* and *cox1f* clades encountered in the *M. b. rubra* hybrid zone all seem to have occurred before the last glacial maximum (*cox1f* clades b1a and b1b between 106 and 185 ky, *cox1m* clades m1 vs m2/m3 between 30 and 50 ky).

Sexing with optical microscopy on 655 individuals from 14 sampling sites of the French Atlantic coast resulted in the detection of 248 males, 185 females and 222 undifferentiated individuals. Undifferentiated specimens might have released their gametes or resorbed them, as lipid droplets were observed in some of the undifferentiated gonads. Molecular sexing was used to detect males in undifferentiated individuals. Both *cox1m* and *cox1f* were successfully sequenced for 302 male individuals from 14 sampling sites (Table 1). Sampling size per site ranged from 3 to 33 and from 28 to 33 when only considering sampling sites in France, apart from Arcachon (ARC, N=10). Among the 302 studied male individuals, 34 and 61 haplotypes were found for *cox1f* and *cox1m*, respectively. Among these, 59% were singletons at *cox1f* and 80% at *cox1m*.

Median joining networks showed a similar geographical pattern with three major haplogroups encompassing mainly the northern/central and southern samples (Figure 1A and C). Nonetheless, the networks revealed important differences between markers. First, *cox1m* presented a higher diversity, but a similar divergence level: a maximum of 25 and 26 mutation steps separated the most divergent haplotypes for *cox1f* and *cox1m*, respectively. We estimated maximum divergence rates of 0.027 s/s and 0.028 s/s for the m and f datasets, respectively.

<figure><img src="image_3.png" /><figcaption>**Figure 1** - Median joining networks constructed from the haplotypes found among the 302 individuals distributed in 14 sampling sites for (A) the *cox1f* gene and (C) the *cox1m* gene. Haplotype associations among haplogroups are represented with a tanglegram (B) Haplotypes are sorted by haplogroup for each marker and a line represent the connection of the *cox1m* and *cox1f* haplotypes within one individual. The color of the line represents the population where it was sampled and bold lines and arrows highlighting individuals sharing distant haplogroups.</figcaption></figure>
<table><caption>**Table 1** - Genetic diversity indices at *cox1f* (479bp) and *cox1m* (479bp) loci for each sampling site: n: number of individuals, nsites: number of polymorphic sites, H: number of haplotypes, π: nucleotide diversity, Hd: Haplotype diversity and D: Tajima's D coefficient.</caption><thead><tr><th rowspan="2">Geographic zone</th><th rowspan="2">Sampling site</th><th colspan="7">cox1f</th><th colspan="3">cox1m</th></tr><tr><th>n</th><th>nsites</th><th>H</th><th>π(10<sup>2</sup>)</th><th>Hd</th><th>D</th><th>nsites</th><th>H</th><th>π(10<sup>2</sup>)</th><th>Hd</th><th>D</th></tr></thead><tbody><tr><td rowspan="6">North</td><td>Sylt (SYL)</td><td>7</td><td>7</td><td>4</td><td>0.91</td><td>-</td><td>-</td><td>5</td><td>4</td><td>0.38</td><td>-</td><td>-</td></tr><tr><td>Wilhelmshaven (WIL)</td><td>3</td><td>0</td><td>1</td><td>0.00</td><td>-</td><td>-</td><td>1</td><td>3</td><td>0.13</td><td>-</td><td>-</td></tr><tr><td>Crildumersiel (CRI)</td><td>10</td><td>6</td><td>2</td><td>0.77</td><td>0.53</td><td>2.10</td><td>14</td><td>9</td><td>0.70</td><td>0.98</td><td>-1.47</td></tr><tr><td>Kruiningen (KRU)</td><td>22</td><td>11</td><td>7</td><td>0.77</td><td>0.69</td><td>0.11</td><td>10</td><td>22</td><td>0.23</td><td>0.65</td><td>-2.00*</td></tr><tr><td>all</td><td>42</td><td>11</td><td>7</td><td>0.75±0.43</td><td>0.64±0.06</td><td>0.75</td><td>23</td><td>16</td><td>0.36±0.24</td><td>0.79±0.05</td><td>-2.26*</td></tr><tr><td>Le Croyot (CRO)</td><td>28</td><td>13</td><td>7</td><td>0.77</td><td>0.58</td><td>-0.24</td><td>14</td><td>14</td><td>0.38</td><td>0.88</td><td>-1.68*</td></tr><tr><td rowspan="2">Center</td><td>Saint-Vaast-la-hougue (VAA)</td><td>31</td><td>11</td><td>5</td><td>0.33</td><td>0.24</td><td>-1.63*</td><td>11</td><td>9</td><td>0.27</td><td>0.75</td><td>-1.69*</td></tr><tr><td>all</td><td>59</td><td>14</td><td>9</td><td>0.52±0.31</td><td>0.42±0.08</td><td>-0.79</td><td>20</td><td>20</td><td>0.32±0.21</td><td>0.81±0.04</td><td>-1.99*</td></tr><tr><td rowspan="4">Hybrid/Finistère</td><td>Mont-Saint-Michel (MSM)</td><td>31</td><td>8</td><td>4</td><td>0.22</td><td>0.24</td><td>-1.60*</td><td>5</td><td>7</td><td>0.11</td><td>0.35</td><td>-1.51*</td></tr><tr><td>Saint-Brieuc (BRI)</td><td>30</td><td>15</td><td>9</td><td>0.76</td><td>0.99</td><td>0.17</td><td>4</td><td>5</td><td>0.08</td><td>0.36</td><td>-1.58*</td></tr><tr><td>all</td><td>61</td><td>16</td><td>10</td><td>0.74±0.42</td><td>0.58±0.07</td><td>-0.28</td><td>8</td><td>10</td><td>0.09±0.09</td><td>0.35±0.08</td><td>-1.91*</td></tr><tr><td>Pont-Mahé (MAH)</td><td>30</td><td>8</td><td>6</td><td>0.18</td><td>0.36</td><td>-1.86*</td><td>5</td><td>6</td><td>0.07</td><td>0.31</td><td>-2.00*</td></tr><tr><td rowspan="3">South</td><td>Saint-Brevin (BRE)</td><td>29</td><td>9</td><td>8</td><td>0.71</td><td>0.49</td><td>-0.25</td><td>8</td><td>9</td><td>0.13</td><td>0.53</td><td>-2.14*</td></tr><tr><td>Noirmoutier (NOI)</td><td>8</td><td>3</td><td>4</td><td>0.23</td><td>-</td><td>-</td><td>1</td><td>2</td><td>0.05</td><td>-</td><td>-</td></tr><tr><td>Aytré (AYT)</td><td>33</td><td>8</td><td>6</td><td>0.43</td><td>0.55</td><td>-0.18</td><td>5</td><td>6</td><td>0.07</td><td>0.28</td><td>-1.88*</td></tr><tr><td rowspan="2">Fouras (FOU)</td><td>30</td><td>8</td><td>7</td><td>0.29</td><td>0.50</td><td>-1.12</td><td>6</td><td>6</td><td>0.08</td><td>0.31</td><td>-2.10*</td></tr><tr><td>all</td><td>63</td><td>11</td><td>10</td><td>0.36±0.23</td><td>0.52±0.07</td><td>-0.89</td><td>9</td><td>9</td><td>0.08±0.08</td><td>0.29±0.07</td><td>-2.15*</td></tr><tr><td>Arcachon (ARC)</td><td></td><td>10</td><td>5</td><td>3</td><td>0.23</td><td>0.38</td><td>-1.74*</td><td>1</td><td>2</td><td>0.04</td><td>0.20</td><td>-1.11</td></tr><tr><td></td><td>all</td><td>302</td><td>34</td><td>34</td><td>0.81±0.45</td><td>0.76±0.02</td><td>-1.06</td><td>53</td><td>61</td><td>0.29±0.20</td><td>0.68±0.03</td><td>-2.44*</td></tr></tbody></table>

Note: * indicates significant values (<0.05, 10,000 permutations), -: not available
The haplogroup relative frequencies showed different spatial patterns (Figure 2A). The *cox1m* haplogroup m3 was present only in the North Sea and east Channel samples, while the *cox1f* haplogroup b2/b3 encompassed individuals from the North Sea to the west English Channel (hybrid zone). Also, the haplogroup b1a of *cox1f* was shared among the North Sea, the English Channel and Atlantic samples, whereas the central haplogroup m2 of *cox1m* was found only in the North Sea and the Channel Sea. Interestingly, clear genetic breaks were revealed for each marker, but they were not located at the same geographical location. For *cox1m*, the Cotentin peninsula represented a sharp transition zone between the haplogroups m3/m2 and m1, which was confirmed by pairwise *Φ*<sub>ST</sub> values (Table S3). Two homogenous populations were observed north of the Cotentin peninsula and south of the Finistère (within group pairwise *Φ*<sub>ST</sub> ranged from -0.05 to 0.04) and these populations were highly differentiated (between group pairwise *Φ*<sub>ST</sub> ranged from 0.38 to 0.61). In between, the BRI and MSM sampling sites were very differentiated (between group pairwise *Φ*<sub>ST</sub> ranged from 0.50 to 0.61) from the north of the Cotentin but they were also weakly yet significantly different from some sites at the south of the Finistère (between group pairwise *Φ*<sub>ST</sub> ranged from 0.02 to 0.04). For *cox1f*, the pattern was more complex. The haplogroup b1a was composed of haplotypes found in majority at the south of the Finistère (Figure 2) and samples from either side of the Finistère were highly divergent (pairwise *Φ*<sub>ST</sub> ranging from 0.41 to 0.79, Table S4). On both sides of the Finistère peninsula, significant population substructure was detectable: the Pont-Mahé Bay (MAH) sampling site was mildly differentiated from Saint-Brevin-les-Pins (BRE) and Aytré (AYT) (pairwise *Φ*<sub>ST</sub> of 0.13 and 0.08 respectively, Table S4). Also, the relative frequencies of the 3 haplogroups were variable in the north of the Finistère peninsula, resulting in a heterogeneous distribution of the genetic diversity (Figure 2).

Saint-Brieuc (BRI) for instance was highly differentiated from all other northern sites except from MSM (Mont-Saint-Michel). Based on these results, genetically homogenous sampling sites were merged in 3 geographical groups (Table 1) and used in subsequent analysis (north, center and south). Sampling sites within the hybrid/Finistère group showed genetic differentiation only for *cox1f*.

Cline analyses performed on haplogroup frequencies consistently revealed clinal variations with the null model being rejected for each haplogroup (Figure 2B and Table S5). For the major northern clades (b1a/m2) and the southern clades (b1b/m1), a geographical discordance in the center of hybrid zone was observed, ranging from 178 (clades b1a vs. m2) to 209 (clades b1b vs. m1) km. In both cases (major northern clades and southern clades), the genetic break was located on the Cotentin Peninsula for *cox1m*, and the western tip of the Finistère Peninsula for *cox1f*. The width estimated from the maximum likelihood cline of the best-fit models were 4 to 33 times narrower at *cox1m* (16 to 80 km), compared to *cox1f* (317 and 540 km for clades b1b and b1a, respectively). Concerning the minor northern clades (b2/b3 and m3), the null model was rejected but the clinal model might not be the best fit to explain its distribution. Indeed, b2/b3 was absent in the south but minor in the northern and central sites (from 6 to 27%), except for CRI and BRI where it was dominant. BRI actually accounts for a large part of b2/b3's diversity, with 5 of the 10 haplotypes in this group, four of which are private to this sampling site. This geographic pattern contrasts with m3, which followed a clinal model: a frequency ranging between 18 and 42% in the north and absent in the south.

Haplotypes from the extreme haplogroups (S and N) were rarely shared within one individual (Figure 1B), and while 17 individuals exhibited the haplogroup b2/b3 at *cox1f* and the haplogroup m1 at *cox1m*, only 3 individuals shared the *cox1f* haplogroup b1b and the *cox1m* haplogroup m3. Also, only 1 individual presented both the *cox1f* haplogroup b1b and the *cox1m* haplogroup m2, suggesting non-random mating or lethal haplotype combinations. Noteworthily, most of these individuals were found in the hybrid zone sampling sites (18 out of 22). These results were confirmed by the significant linkage disequilibrium test carried out among haplogroups for phased data in the hybrid zone and the south (T² test df=6, p<0.001).
<figure><img src="image_4.png" /><figcaption>**Figure 2** - Geographical distribution of the *cox1f* and *cox1m* haplogroups: (A) haplogroup frequencies are represented for each marker in all the sites with a pie chart. The size of the pie chart indicates the number of individuals analyzed, (B) cline analyses comparing male- and female-type haplogroups. The transparent colored area around the fitted curve represents the 95% confidence interval. The haplogroup m4 is not presented here as it is only found in one individual from GRI.</figcaption></figure>
The level of genetic diversity also differed between the two markers (Table 1 and Figure 3). While the average nucleotide diversity (π) and its range across sampling sites were similar (a mean value of 0.0047 and a range of 0 to 0.0091 for *cox1f*, a mean of 0.0019 and a range 0.0004 to 0.007 for *cox1m*), the median values of the nucleotide diversity (0.38 and 0.12 for *cox1f* and *cox1m* respectively) showed that the nucleotide diversity is significantly greater for *cox1f* than for *cox1m* (Wilcoxon rank test V= 65 p= 0.002). We observed the lowest nucleotide diversities in southern sampling sites for *cox1m* compared to *cox1f*. The haplotype diversity was not significantly different between the two markers (Wilcoxon rank test, V = 35, p = 0.89). Finally, Tajima's D also differed among the two markers. When considering the three geographical groups, significant negative values were found only for *cox1m* indicating the possibility of population expansion after a bottleneck event or selective sweep. For *cox1f*, on the other hand, no sign of deviation from neutral expectations was detected at the population level. The proportion of singletons at the scale of sample sites (Figure 3) showed an overall higher number of rare frequency mutations in *cox1m* compared to *cox1f*, which combined with negative D values, could indicate a recent expansion in population size. Mismatch distribution analyses and drift-mutation equilibrium tests suggest different demographic histories at *cox1f* and *cox1m*, with bimodal pairwise nucleotide distributions at the former, and unimodal right-skewed distributions at the latter (Figure S4 and Table S5). Based on the non-significant raggedness indices (r), for both markers and groups, the hypothesis of sudden population expansion cannot be rejected; this result is only substantiated at *cox1m* (in both groups), by the R² statistic. This index estimated as the difference between the number of singleton mutations and the average number of nucleotide differences among haplotypes has been shown as the most powerful to detect population expansion in case of small sample size and few segregating sites (Ramos-Onsins & Rozas, 2002). While molecular signatures of the mutation/drift (dis)equilibrium are similar between southern and northern groups regardless of the marker, the mismatch distributions show more divergence among haplotypes in the north (pairwise nucleotide differences >6), compared to the south (pairwise nucleotide differences <6).

Consistently with the haplotype networks, the results of the AMOVAs suggested different geographic structuration patterns for *cox1m* and *cox1f* (Table 2). The best hierarchical model based on the F<sub>SC</sub>/F<sub>CT</sub> ratio (least variation within and highest between groups) was different for the two markers: for *cox1f* the best model was two groups (B) separated by the Finistère peninsula (KRU, CRI, CRO, VAA, MSM, BRI vs MAH, BRE, AYT, FOU, ARC) whereas for *cox1m*, the best model was also two groups of subpopulations but (A) with the Cotentin peninsula as a major geographical break (KRU, CRI, CRO, VAA) vs (MSM, BRI, MAH, BRE, AYT, FOU, ARC). The result of the AMOVAs (Table 2) showed similar level of among subpopulation genetic variation relative to the total variation for both markers (*cox1m*: 42% and Φ<sub>ST</sub>= 0.58, p<0.001; *cox1f*: 44% and Φ<sub>ST</sub>= 0.55, p<0.001). The level of differentiation among groups was higher for *cox1m* compared to *cox1f* (Φ<sub>CT</sub>= 0.57, p<0.001 and Φ<sub>CT</sub>= 0.49, p<0.001 for *cox1m* and *cox1f*, respectively). Finally, within group genetic structure was detected only for *cox1f* (Φ<sub>SC</sub>= 0.01, p=0.062 and Φ<sub>SC</sub>= 0.13, p<0.001 for *cox1m* and *cox1f*, respectively). Genetic differentiation analyses based only on haplotype frequencies (conventional F-statistics) were also carried out and gave similar results (not shown).

Finally, our results suggest that either Ne and/or μ differed geographically. The theta estimates were twice larger at *cox1m* compared to *cox1f* in the north (θf/θm =0.48) and similar in the south (θf/θm =0.98).
<figure><img src="image_5.png" /><figcaption>**Figure 3** - Genetic diversity indexes for both markers (*cox1m* and *cox1f*) along the coast. The gray area indicates the geographic range between the Cotentin and the Finistère Peninsulas.</figcaption></figure>
<table><caption>**Table 2** - Analysis of molecular variance (AMOVA) to partition the genetic variation among hierarchical geographical scales. Results are presented for both markers, with different options of hierarchical grouping of subpopulations (as detailed in the methods).</caption><thead><tr><th rowspan="2">Number of<br>groups</th><th rowspan="2">Groups</th><th colspan="4">Cox1f</th><th colspan="4">Cox1m</th></tr><tr><th>Among groups<br>relative to total<br>variance</th><th>Among<br>subpopulations<br>within groups</th><th>Among<br>subpopulations<br>relative to total<br>variance</th><th>Among groups<br>relative to total<br>variance</th><th>Among<br>subpopulations<br>within groups</th><th>Among<br>subpopulations<br>relative to total<br>variance</th><th>Among<br>subpopulations<br>within groups</th></tr></thead><tbody><tr><td rowspan="2">3</td><td>(CRI, KRU) (CRO,<br>VAA, MSM, BRI)<br>(MAH, BRE, AYT,<br>FOU, ARC)</td><td>%var</td><td>44.71</td><td>7.00</td><td>48.29</td><td>22.2</td><td>22.54</td><td>55.26</td><td></td></tr><tr><td>Fixation<br>indices</td><td></td><td>Φ<sub>CT</sub>=0.447***</td><td>Φ<sub>SC</sub>=0.517***</td><td>Φ<sub>ST</sub>=0.126**</td><td>Φ<sub>CT</sub>=0.221<sup>ns</sup></td><td>Φ<sub>SC</sub>=0.290***</td><td>Φ<sub>ST</sub>=0.447***</td><td></td></tr><tr><td rowspan="2">2(A)</td><td>(CRI, KRU, CRO,<br>VAA) (MSM, BRI,<br>MAH, BRE, AYT,<br>FOU, ARC)</td><td>%var</td><td>16.98</td><td>31.02</td><td>52.00</td><td>57.36</td><td>0.56</td><td>42.08</td><td></td></tr><tr><td>Fixation<br>indices</td><td></td><td>Φ<sub>CT</sub>=0.170<sup>ns</sup></td><td>Φ<sub>SC</sub>=0.374***</td><td>Φ<sub>ST</sub>=0.480***</td><td>Φ<sub>CT</sub>=0.574***</td><td>Φ<sub>SC</sub>=0.013<sup>ns</sup></td><td>Φ<sub>ST</sub>=0.579***</td><td></td></tr><tr><td rowspan="2">2(B)</td><td>(CRI, KRU, CRO,<br>VAA, MSM, BRI)<br>(MAH, BRE, AYT,<br>FOU, ARC)</td><td>%var</td><td>49.14</td><td>6.40</td><td>44.46</td><td>21.88</td><td>4.18</td><td>53.94</td><td></td></tr><tr><td>Fixation<br>indices</td><td></td><td>Φ<sub>CT</sub>=0.491***</td><td>Φ<sub>SC</sub>=0.126***</td><td>Φ<sub>ST</sub>=0.555**</td><td>Φ<sub>CT</sub>=0.212*</td><td>Φ<sub>SC</sub>=0.309***</td><td>Φ<sub>ST</sub>=0.460***</td><td></td></tr></tbody></table>

Note: *P>0.05, **P>0.01, ***P < 0.001, ns: non significant
# Discussion

## Phylogeography of *Macoma* lineages and F/M geographic discordance

Recovering the M- and F-type *cox1* of *M. b. balthica* from the transcriptomic study of Yurchenko et al., (2018) allowed the application of the net divergence method used in Nikula et al., (2007) on both *cox1f* and *cox1m* datasets. As detailed above, this clock hypothesizes that *rubra* and *balthica* diverged between 3.5-2 Mya (bounded by the opening of the Bering Strait and the discovery of first Atlantic fossils, respectively) and was applied by Nikula et al., (2007) on *cox3f*. At the scale of the genus, *cox1f* and *cox1m* data confirmed that the divergence between *M. balthica* and *M. petalum* occurred at the end of the Miocene, early Pliocene at the latest (*cox1f*: 4.9-8.6 My; *cox1m*: 5.4-9.5 My), as inferred with *cox3f* (5.2-9.3 My; Nikula et al., 2007). Within *M. balthica*, estimates at *cox1f*, *cox1m* and *cox3f* congruently suggest a split between clades C and D during the Calabrian subdivision of the Pleistocene (0.83-1.76, 0.98-1.72, 0.8-1.4 My, respectively), again confirming previous results (Nikula et al., 2007).

In the north (Cildumersiel, Germany) and Baltic Seas (Mecklenburg bight, Germany; Lomma, Sweden, in the Öresund Strait; Umea, Sweden, in the Gulf of Bothnia), a distinct *M. b. rubra* clade here called m4 was sampled at *cox1m* (Table S2, Figure S2). This clade diverged from m3 about 220 to 380 kya and is probably the M-type counterpart of b2/b3, which diverged from b1a about 251 to 440 kya. Males with m4 haplotypes possessed either b2/b3 or D at *cox1f*. The b2/b3 clade is common in North Sea, the Öresund and Western Baltic (Nikula et al., 2007), making its association with m4 expected. However, six m4 specimens were sampled from Umea, five of which bearing a *M. b. balthica* signature at *cox1f* (clade D), showing that *rubra* m4 introgressed far north into the Baltic Sea, as does b2/b3 (Nikula et al., 2007).

The split between the southern and northern *M. b. rubra* clades, currently separated by the Gulf of Saint-Malo where they hybridize, occurred much more recently. At *cox1f*, divergence was estimated between 106 and 185 kya, in opposition with the much later split estimated at *cox1m* (30-50 kya). These different split times may be due to the different forces acting upon the M- and F-type mtDNA, namely the selection / mutation equilibrium (known to differ between M- and F-type mtDNA, see below), drift (which may have been substantial during glacial cycles), assortative mating (a pre-zygotic barrier, see Bierne et al., 2002 for an example in *Mytilus*), sex-biased migration (see Arnaud-Haond et al., 2003 for an example in oysters), and different hybrid fitness between sexes. Nevertheless, both split times are consistent with a scenario of pre-LGM vicariance (the LGM is currently estimated to have occurred between 26.5 and 19 kya; Clark et al., 2009), in opposition to a post-glacial scenario of primary intergradation.

Zooming in on the Gulf of Saint-Malo hybrid zone revealed discordant phylogeographic patterns for the F- and M-type mitochondrial *cox1* gene. Both mtDNAs displayed significant segregation of haplotypes in space. As previously described by Becquet et al., (2012), a geographical break at the Finistère peninsula was observed at *cox1f*. Sampling done for the present study narrows the gap between northern and southern populations by about 50 km thanks to the addition of the site of Saint-Brieuc Bay, however we failed at sampling any large population between the latter and Pont-Mahé Bay, south of Finistère (Figure 2). While anecdotal data support the presence of *M. b. rubra* in the Bay of Brest (Hily, 2013 and our observations), the Finistère intertidal seem largely inhospitable to this species due to the scarcity of mud flats (Becquet et al., 2012, Le Mao et al., 2020 and this study). It therefore appears as a credible physical barrier to genetic connectivity, and it has been indeed identified in the past as a genetic breaking zone for other marine benthic-pelagic species confined to estuaries and bays (muddy-fine sediment species, Jolly et al., 2005, 2006).

On the other hand, the Cotentin peninsula located about 250 km farther East of the Finistère peninsula separated southern and northern populations *cox1m* haplotypes at a much narrower hybrid zone (an observation that might be interpreted by a different selection / migration balance compared to *cox1f*, or the heterogeneous spacing of sampling points along the coast). Sandy mud flats are numerous along western Cotentin but do not seem hospitable to *M. balthica*. Our field surveys failed at identifying populations between Mont-Saint-Michel Bay (MSM) and Saint-Vaast-la-Hougue (VAA), although anecdotal presence of the species was recorded farther along the
Cotentin coast and was reported to be present but rare on the Gulf of Saint-Malo (Le Mao et al., 2020). Along with the scarcity of seemingly appropriate habitat and patchiness of small populations in the intertidal zone of Cotentin, we observed significant differences in gonadal development and spawning phenology between populations on either side of the peninsula. While individuals with full gonads and well-developed gametes were sampled between the 4th and the 23rd of April 2018 from Fouras to Mont-Saint-Michel Bay (some individuals having even spawned in the field at MSM and Saint-Brieuc), gonads from Saint-Vaast-la-Hougue (VAA) and Somme Bay individuals were still in early stages of development by April 19, 2018 (most gametes being still undifferentiated). These observations were on par with previous reports suggesting spawning asynchrony across the Cotentin Peninsula, with spawning occurring in early April in Aytré (Saunier, 2015) and August in Somme Bay (Ruellet, 2013). In the great scallop, *Pecten maximus*, different spawning phenologies were described across the Cotentin peninsula (Lubet et al., 1991). As for the Finistère, the Cotentin peninsula was previously identified as a barrier to gene flow (Jolly et al., 2005; Quéméré et al., 2016; Handal et al., 2020). It is also the southern biogeographical boundary of the northern cold temperate Boreal region (Dinter, 2001) with distinct hydrologic and oceanographic features from the western part of the English Channel (Dauvin, 2012). There are, therefore, multiple exogenous (habitat availability; oceanographic currents, e.g. Dupont et al., 2007; Fievet et al., 2007; Hily, 2013) and endogenous (asynchrony in reproductive phenology and genetic incompatibilities, discussed below) barriers that could contribute to the genetic differentiation of populations across these two peninsulas.

The fact that the geographic discordance of genetic clines is sex-specific is quite original for a marine invertebrate. In the DUI literature, studies have reported that genetic differentiation across populations is stronger for M-type mtDNA than for F-type mtDNA (Liu et al., 1996; Riginos et al., 2004; Skibinski et al., 1999), and some authors reported discordant genetic structure between markers. The freshwater mussels *Lampilsis siliquoidea* displayed population genetic structure at *cox1f* but none at *cox1m* (Krebs et al., 2013). Conversely, Riginos et al., (2004) found weak connectivity at F-type mtDNA but none at M-type mtDNA in *Mytilus edulis* across the Atlantic Ocean. At the entrance of the Baltic, rampant introgression of *M. edulis* F-type mtDNA was observed, compared to a sharp cline at M-type mtDNA concordant with M7 lysin (a nuclear gene involved in fertilization; Stuckas et al., 2009). This makes the population structure described here for *M. b. rubra* noteworthy in the sense that the phylogeography of F- and M-type mtDNA markers do not simply differ in amplitude or resolution (due to significant differences in mutation rates, µM being about twice µF) but differ in the geographical position of the haplotypic cline separating southern and northern *M. b. rubra*.

While a genetic break located around the Finistère peninsula corresponds to previous results by Becquet et al., (2012) for *cox1f* and eight nuclear microsatellites, the Cotentin was not detected as a barrier in that study. However, spatial structure was detected at the nuclear *atp5c1* gene (encoding the gamma subunit of the FO/F1 ATP synthase protein complex) between Mont-Saint-Michel and Somme Bay (Saunier, 2015). Non-synonymous point mutations separating southern and northern populations are located in an inter-peptide interaction domain of the gamma subunit, suggesting mito-nuclear incompatibilities (Pante et al., 2019). The concordance of *atp5c1* and *cox1m* calls for further investigation of incompatibilities involving the male mitogenome specifically.

Interestingly, significant differentiation with *cox1f* was detected between Pont-Mahé Bay (MAH) and the Saint Brévin (BRE) and Aytré (AYT) populations (Table S4). These two groups, separated by the Vilaine River, were identified as significantly differentiated based on microsatellites, but not *cox1f* in the study of Becquet et al., (2012), which was based on individuals from both sexes, while our study focused on males exclusively. We may therefore benefit from higher statistical power to detect differences in differentiation across the Vilaine River. This implies that genetic structure among males at *cox1f* is higher than in females. Although we did not include females in our study, we did compare pairwise *Φ*<sub>ST</sub> values between sampling sites that were common between our study and Becquet et al.’s, (2012). *F*<sub>ST</sub> and *Φ*<sub>ST</sub> were both significantly higher in the dataset composed of males exclusively, compared to the dataset composed of males and females (Figure 4). This pattern can be attributed to stochastic genetic variation (effects of mutations and drift) and/or
differences of gene flow between sexes. The geographical discordance observed at the Finistère and Cotentin peninsulas for *cox1f* and *cox1m* are also reminiscent of sex-level differences. Indeed, in studies using female-transmitted mtDNA and male-transmitted Y chromosome to investigate population structure, sex-biased asymmetries are often cited as a determinant of discordant geographic patterns (e.g. Boissinot & Boursot, 1997; Trejo-Salazar et al., 2021). Sex-level differences in gene flow are scarcely reported but exist in bivalves, as in the protandrous pearl oyster *Pinctada mazatlanica*, which effective sex-ratio is strongly biased towards males (Arnaud-Haond et al., 2003).

<figure><img src="image_6.png" /></figure>

**Figure 4** - Comparison of pairwise FST estimated at the *cox1f* gene for study sites that were both sampled by Becquet et al. (2012) (male and female individuals) and for the present study (male individuals only). The plain line represents the x=y axis and the dash line, the linear regression (y = 1.27 x + 0.03, df = 62, R²adj=0.74, p<0.001). A Mantel test revealed a significant correlation between the two datasets (r=0.786, p=0.0084).

Other factors could explain the M- vs. F-mtDNA discordance, such as mtDNA introgression and heterogeneity in its rate, hybrid zone spatial instability and drift (Barton & Hewitt, 1985; Bierne et al., 2003; Toews & Brelsford, 2012). Patterns of asymmetric introgression of alleles are commonly observed in hybrid zones and were observed in the Bay of Biscay between *Mytilus edulis* and *M. galloprovincialis* (Rawson & Hilbish, 1995). In *M. b. rubra*, asymmetrical gene flow from north to south was detected at *cox1f* (Becquet et al., 2012, this study), and *atp5c1* (Saunier, 2015; Pante et al., 2019). At that same locus, Gagnaire et al., (2012), looking at divergence and selection among populations of the eels *Anguilla rostrata* and *A. anguilla* (two hybridizing sister species) also found evidence of unidirectional introgression from the former to the later. At *cox1m*, little introgression was detectable in *M. b. rubra*, but nevertheless sensibly higher from north to south.

Hybrid zone movement has been recognized for long, as a manifestation of variation in population density, dispersal rate, individual fitness, effects of allele frequencies on population structure, and spatially- or frequency-dependent selection (Barton & Hewitt, 1985). It has been proposed that the northeastern Atlantic *Mytilus* hybrid zones move parallel to warming sea surface temperatures (Hilbish et al., 2012). While a similar phenomenon may occur in Baltic tellins, suggesting that genotype-environment correlations could affect cline geography (Ducos et al.
unpublished), neutral forces may be sufficient to result in discordant M-type and F-type mtDNA cline centers and width. Indeed, as detailed above, both the Finistère and the Cotentin peninsulas are characterized by low population densities and oceanic currents disrupting along-coast dispersal. While mitochondrial recombination is expected to be rare (or at least, rarely transmitted to offspring; Passamonti et al., 2003), looking at whole mitogenomes as well as nuclear genes of mitochondrial function vs. other (putatively neutral) nuclear markers may help shed light on the degree to which the patterns detected at *cox1* and *atp5c1* depart from neutral expectations.

Whether the asymmetric discordant clines in *M. balthica* observed here are due to cline movement, differential gene flow, population structure, or selection (Barton & Hewitt, 1985) remains a fascinating avenue for future research. In particular, introgression across the hybrid zone in the English Channel calls for a test of asymmetric fitness of inter-lineage crosses (e.g., Turelli & Moyle, 2007). An involvement of mito-nuclear incompatibilities could explain the observed asymmetric allele frequencies if northern hybrids with southern mitochondria are less fit than southern hybrids with northern mitochondria.

## Genetic diversity differences between *cox1m* and *cox1f*

F- and M- mtDNA diversity was compared among homologous regions within heteroplasmic males, and therefore variability can be strictly attributed to heredity rather than inter-individual or inter-genic differences. M-type mtDNA was almost up to twice as polymorphic as its female counterpart, with more haplotypes, more singletons (i.e. haplotypes represented by a single male) and more segregating sites. Higher polymorphism in the M-type mtDNA, both within and between species (Skibinski et al., 1999), is commonly reported in DUI species, as in *Polititapes rhomboides* (Chacón et al., 2020), *Pyganodon grandis* (Krebs, 2004) and *Mytilus* sp. (Ort & Pogson, 2007; Śmietanka et al., 2009, 2013; Śmietanka & Burzyński, 2017). Higher frequency of rare M-type haplotypes (80% of singletons compared to 57% among F-type haplotypes in *M. b. rubra*) is another largely shared feature among DUI species (Ort & Pogson, 2007; Śmietanka et al., 2009, 2013; Śmietanka & Burzyński, 2017; Chacón et al., 2020). Major determinants of genetic diversity in animals are linked selection, mutation, and effective population size (Ellegren & Galtier, 2016). In DUI bivalves, polymorphism might result from different selection pressures on male- and female-type mitogenomes (e.g. Milani et al., 2012): female-type mitochondria ensure proper mitochondrial function at the scale of the individual whereas male-type mitochondria, which solely co-occurs in heteroplasmic individuals with the F-type mitochondria (Skibinski et al., 1994; Zouros et al., 1994) have to effectively function exclusively in the male germ line and gametes (different ‘arenas of selection’ sensu, Stewart et al., 1996). According to the male-driven evolution hypothesis, higher germ cell divisions in males (Shimmin et al., 1993) could result in an enhanced replication rate of the M-type mitochondria during spermatogenesis, providing more opportunities for mutations to accumulate (Rawson & Hilbish, 1995; Stewart et al., 1995). In addition, oxidative damage might be important in sperm (Zouros, 2013), and/or selection pressures on chaperonins and DNA repair genes might be relaxed. We can also expect sex-linked demographic differences to influence polymorphism, if the effective population size differs between males and females; we did observe larger theta at *cox1m* in northern populations, suggesting that the effective population size and/or mutation rate is larger in males. Most hypotheses on the role of selection, mutation and effective population size in shaping genetic diversity in DUI bivalves remain to be empirically tested.

However, counter examples do exist. For instance, Beauchamp et al., (2020) found more F-type haplotypes than M-type haplotypes in *Pyganodon grandis* (11 and 7, respectively) and *P. lacustris* (20 and 6, respectively). Śmietanka et al., (2013) (87 F- and 76 M-type haplotypes in *Mytilus trossulus*) and Riginos et al., (2004) (132 F- and 56 M-type haplotypes in *Mytilus edulis*) reported the same pattern, although the latter results might be due to the difficulty of amplifying the male-type mtDNA. Although often presented as such in the literature, higher polymorphism in males is therefore not a universal feature of the population genetics of DUI species.

It is worth noting that, among the studies reporting M/F differences in genetic polymorphism, all of those that explicitly tested for non-equilibrium population state found support for population expansion (Chacón et al., 2020; Ort & Pogson, 2007; Riginos et al., 2004; Śmietanka et al., 2009,
2013; Śmietanka & Burzyński, 2017, and our study at *cox1m*), except for one finding support for purifying selection (Skibinski et al., 1994). Non-equilibrium states following a glacial period can indeed affect genetic diversity, notably through drift, and confound studies focusing on the comparative molecular evolution of F and M mtDNA. It would therefore be interesting to look for biological models at equilibrium.

**Data, scripts, and supplementary information availability**

*Cox1f* and *cox1m-long* sequences were submitted to Genbank (accession N° OM855617 - OM855929 and OM856027 - OM856339 respectively). Data, scripts, and supplementary information are available online (https://doi.org/10.5281/zenodo.14849343, Le Cam et al., 2025).

**Acknowledgements**

The authors thank Mélanie Rocroy, Patrick Triplet, Thierry Ruellet, Rolet Céline (Association GEMEL, Picardie), Anthony Sturbois (RNN Baie de Saint-Brieuc, VivArmor Nature) and Xavier de Montaudouin (Station Marine d'Arcachon, UMR 5805 EPOC, Université de Bordeaux/CNRS) for their help with sample collection. We thank the participants of the Biocombe project ("The impact of BIOdiversity changes in COastal Marine Benthic Ecosystems;" contract n. EVK3-CT-2002-00072; Research n. OND1290723). Lab work was performed at the Molecular Core Facility at the University of La Rochelle. We thank the bioinformatics platform Toulouse Midi-Pyrénées (Bioinfo GenoToul). Preprint version 4 of this article has been peer-reviewed and recommended by Peer Community In Evolutionary Biology (https://doi.org/10.24072/pci.evolbiol.100540; Castilho, 2025).

**Funding**

This work was funded by the ANR (DRIVE project, grant n. ANR-18-CE02-0004-01) and by the contrat de plan Etat-Région (CPER/FEDER) ECONAT (RPC DYPOMAR).

**Conflict of interest disclosure**

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.

**References**

Arnaud-Haond S, Monteforte M, Blanc F, Bonhomme F (2003) Evidence for male-biased effective sex ratio and recent step-by-step colonization in the bivalve *Pinctada mazatlanica*. *Journal of Evolutionary Biology*, **16**, 790–796. https://doi.org/10.1046/j.1420-9101.2003.00603.x

Bandelt H, Forster P, Rohl A (1999) Median-joining networks for inferring intraspecific phylogenies. *Molecular Biology and Evolution*, **16**, 37–48. https://doi.org/10.1093/oxfordjournals.molbev.a026036

Barberousse A, Bierne N, Britton-Davidian J, Capy P, Desdevises Y, Giraud T, Jousselin E, Moulia C, Ravigné V, Smadja C, Samadi S (2010) Chapitre 4. La spéciation (dir. Ravigné, V). In: *Biologie Evolutive* (eds Thomas F, Lefèvre T, Raymond M), p. . De Boeck Supérieur.

Barton NH, Hewitt GM (1985) Analysis of hybrid zones. *Annual Review of Ecology and Systematics*, **16**, 113–148. https://doi.org/10.1146/annurev.es.16.110185.000553

Beauchamp KC, Beyett TW, Scott MW, Zanatta DT (2020) Detection of hybrid *Pyganodon grandis* and *P. lacustris* (Bivalvia: Unionidae) using F-and M-type mitochondrial DNA sequences and geometric morphometrics. *Journal of Molluscan Studies*, **86**, 233–239. https://doi.org/10.1093/mollus/eyaa013

Becquet V, Simon-Bouhet B, Pante E, Hummel H, Garcia P (2012) Glacial refugium versus range limit: conservation genetics of *Macoma balthica*, a key species in the Bay of Biscay (France).
Journal of Experimental Marine Biology and Ecology, **432**, 73-82.
https://doi.org/10.1016/j.jembe.2012.07.008

Bierne N, Borsa P, Daguin C, Jollivet D, Viard F, Bonhomme F, David P (2003) Introgression patterns in the mosaic hybrid zone between *Mytilus edulis* and *M. galloprovincialis*. *Molecular Ecology*, **12**, 447-461. https://doi.org/10.1046/j.1365-294X.2003.01730.x

Bierne N, David P, Boudry P, Bonhomme F (2002) Assortative fertilization and selection at larval stage in the mussels *Mytilus edulis* and *M. galloprovincialis*. *Evolution*, **56**, 292-298.
https://doi.org/10.1111/j.0014-3820.2002.tb01339.x

Boissinot S, Boursot P (1997) Discordant phylogeographic patterns between the Y chromosome and mitochondrial DNA in the house mouse: selection on the Y chromosome? *Genetics*, **146**, 1019-1034.
https://doi.org/10.1093/genetics/146.3.1019

Breton S, Beaupre HD, Stewart DT, Hoeh WR, Blier PU (2007) The unusual system of doubly uniparental inheritance of mtDNA: isn't one enough? *Trends in genetics*, **23**, 465-474.
https://doi.org/10.1016/j.tig.2007.05.011

Breton S, Stewart DT, Shepardson S, Trdan RJ, Bogan AE, Chapman EG, Ruminas AJ, Piontkivska H, Hoeh WR (2011) Novel protein genes in animal mtDNA: a new sex determination system in freshwater mussels (Bivalvia: Unionoida)? *Molecular Biology and Evolution*, **28**, 1645-1659.
https://doi.org/10.1093/molbev/msq345

Burton RS, Barreto FS (2012) A disproportionate role for mtDNA in Dobzhansky-Muller incompatibilities? *Molecular Ecology*, **21**, 4942-4957.
https://doi.org/10.1111/mec.12006

Burton RS, Pereira RJ, Barreto FS (2013) Cytonuclear genomic interactions and hybrid breakdown. *Annual Review of Ecology, Evolution, and Systematics*, **44**, 281-302.
https://doi.org/10.1146/annurev-ecolsys-110512-135758

Capt C, Bouvet K, Guerra D, Robicheau BM, Stewart DT, Pante E, Breton S (2020) Unorthodox features in two venerid bivalves with doubly uniparental inheritance of mitochondria. *Scientific Reports*, **10**, 1087. https://doi.org/10.1038/s41598-020-57975-y

Capt C, Renaut S, Ghiselli F, Milani L, Johnson NA, Sietman BE, Stewart DT, Breton S (2018) Deciphering the Link between Doubly Uniparental Inheritance of mtDNA and Sex Determination in Bivalves: Clues from Comparative Transcriptomics. *Genome Biology and Evolution*, **10**, 577-590.
https://doi.org/10.1093/gbe/evy019

Capt C, Renaut S, Stewart DT, Johnson NA, Breton S (2019) Putative mitochondrial sex determination in the Bivalvia: insights from a hybrid transcriptome assembly in freshwater mussels. *Frontiers in genetics*, **10**, 840. https://doi.org/10.3389/fgene.2019.00840

Castilho R (2025) Unraveling the complexities of mitochondrial inheritance in *Macoma balthica*: Insights from doubly uniparental inheritance. *Peer Community in Evolutionary Biology*, 100540.
https://doi.org/10.24072/pci.evolbiol.100540

Chacón GM, Arias-Pérez A, Freire R, Martínez L, Nóvoa S, Naveira H, Insua A (2020) Evidence of doubly uniparental inheritance of the mitochondrial DNA in *Polititapes rhomboides* (Bivalvia, Veneridae): Evolutionary and population genetic analysis of F and M mitotypes. *Journal of Zoological Systematics and Evolutionary Research*, **58**, 541-560.
https://doi.org/10.1111/jzs.12267

Clark PU, Dyke AS, Shakun JD, Carlson AE, Clark J, Wohlfarth B, Mitrovica JX, Hostetler SW, McCabe AM (2009) The last glacial maximum. *science*, **325**, 710-714.
https://doi.org/10.1126/science.1172873

Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: more models, new heuristics and parallel computing. *Nature Methods*, **9**, 772. https://doi.org/10.1038/nmeth.2109

Dauvin J-C (2012) Are the eastern and western basins of the English Channel two separate ecosystems? *Marine Pollution Bulletin*, **64**, 463-471.
https://doi.org/10.1016/j.marpolbul.2011.12.010

Derryberry EP, Derryberry GE, Maley JM, Brumfield RT (2014) HZAR: hybrid zone analysis using an R software package. *Molecular Ecology Resources*, **14**, 652-663.
https://doi.org/10.1111/1755-0998.12209
Dinter WP (2001) Biogeography of the OSPAR maritime area. *Bonn, Germany: Federal Agency for Nature Conservation*, **167**.

Dupont L, Ellien C, Viard F (2007) Limits to gene flow in the slipper limpet *Crepidula fornicata* as revealed by microsatellite data and a larval dispersal model. *Marine Ecology Progress Series*, **349**, 125–138. https://doi.org/10.3354/meps07098

Ellegren H, Galtier N (2016) Determinants of genetic diversity. *Nature Reviews Genetics*, **17**, 422–433. https://doi.org/10.1038/nrg.2016.58

Excoffier L, Smouse PE, Quattro JM (1992) Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. *Genetics*, **131**, 479–491. https://doi.org/10.1093/genetics/131.2.479

Fievet V, Touzet P, Arnaud J, Cuguen J (2007) Spatial analysis of nuclear and cytoplasmic DNA diversity in wild sea beet (*Beta vulgaris ssp. maritima*) populations: do marine currents shape the genetic structure? *Molecular Ecology*, **16**, 1847–1864. https://doi.org/10.1111/j.1365-294X.2006.03208.x

Fu Y-X (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. *Genetics*, **147**, 915–925. https://doi.org/10.1093/genetics/147.2.915

Gagnaire P-A, Normandeau E, Bernatchez L (2012) Comparative genomics reveals adaptive protein evolution and a possible cytonuclear incompatibility between European and American eels. *Molecular Biology and Evolution*, **29**, 2909–2919. https://doi.org/10.1093/molbev/mss076

Guerra D, Plazzi F, Stewart DT, Bogan AE, Hoeh WR, Breton S (2017) Evolution of sex-dependent mtDNA transmission in freshwater mussels (Bivalvia: Unionida). *Scientific Reports*, **7**, 1–13. https://doi.org/10.1038/s41598-017-01708-1

Guindon S, Gascuel OA (2003) A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. *Systematic Biology*, **52**, 696–704. https://doi.org/10.1080/10635150390235520

Gusman A, Lecomte S, Stewart DT, Passamonti M, Breton S (2016) Pursuing the quest for better understanding the taxonomic distribution of the system of doubly uniparental inheritance of mtDNA. *PeerJ*, **4**, e2760. https://doi.org/10.7717/peerj.2760

Handal W, Szostek C, Hold N, Andrello M, Thiébaut E, Harney E, Lefebvre G, Borcier E, Jolivet A, Nicolle A, Boyé A, Foucher E, Boudry P, Charrier G (2020) New insights on the population genetic structure of the great scallop (*Pecten maximus*) in the English Channel, coupling microsatellite data and demogenetic simulations. *Aquatic Conservation: Marine and Freshwater Ecosystems*, **30**, 1841–1853. https://doi.org/10.1002/aqc.3316

Harpending H (1994) Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. *Human biology*, **591–600**. https://www.jstor.org/stable/41465371

Hilbish TJ, Lima FP, Brannock PM, Fly EK, Rognstad RL, Wethey DS (2012) Change and stasis in marine hybrid zones in response to climate warming. *Journal of Biogeography*, **39**, 676–687. https://doi.org/10.1111/j.1365-2699.2011.02633.x

Hily C (2013) Redécouverte de *Macoma balthica* (Linnaeus, 1758) en rade de Brest. *An aod - les cahiers naturalistes de l'Observatoire marin*, **2**, 31–35.

Huber M (2010) Compendium of bivalves. A full-color guide to 3,300 of the world's marine bivalves. A status on Bivalvia after 250 years of research. *ConchBooks*.

Huber M, Langleit A, Kreipl K (2015) Compendium of bivalves. In: (ed M. Huber), p. 907. *ConchBooks*, Harxheim.

Hummel H, Bogaards R, Bek T, Polishchuk L, Amiard-Triquet C, Bachelet G, Desprez M, Strelkov P, Sukhotin A, Naumov A, Dahle S, Denisenko S, Gantsevich M, Sokolov K, Wolf L de (1997) Sensitivity to stress in the bivalve *Macoma balthica* from the most northern (Arctic) to the most southern (French) populations: low sensitivity in Arctic populations because of genetic adaptations? *Hydrobiologia*, **355**, 127–138. https://doi.org/10.1007/978-94-017-1907-0_13

Jolly MT, Jollivet D, Gentil F, Thiébaut É, Viard F (2005) Sharp genetic break between Atlantic and English Channel populations of the polychaete *Pectinaria koreni*, along the North coast of France. *Heredity*, **94**, 23–32. https://doi.org/10.1038/sj.hdy.6800543
Jolly MT, Viard F, Gentil F, Thiébaut E, Jollivet D (2006) Comparative phylogeography of two coastal polychaete tubeworms in the Northeast Atlantic supports shared history and vicariant events. *Molecular Ecology*, **15**, 1841–1855. https://doi.org/10.1111/j.1365-294x.2006.02910.x

Kalyaanamoorthy S, Minh BQ, Wong TK, Von Haeseler A, Jermiin LS (2017) ModelFinder: fast model selection for accurate phylogenetic estimates. *Nature methods*, **14**, 587–589.
https://doi.org/10.1038/nmeth.4285

Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thierer T, Ashton B, Meintjes P, Drummond A (2012) Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. *Bioinformatics*, **28**, 1647–1649. https://doi.org/10.1093/bioinformatics/bts199

Krebs R (2004) Combining paternally and maternally inherited mitochondrial DNA for analysis of population structure in mussels. *Molecular Ecology*, **13**, 1701–1705.
https://doi.org/10.1111/j.1365-294X.2004.02133.x

Krebs R, Borden WC, Evans NM, Doerder FP (2013) Differences in population structure estimated within maternally- and paternally-inherited forms of mitochondria in *Lampilis siliquoidea* (Bivalvia: Unionidae). *Biological Journal of the Linnean Society*, **109**, 229–240.
https://doi.org/10.1111/bij.12025

Layton KK, Martel AL, Hebert PD (2016) Geographic patterns of genetic diversity in two species complexes of Canadian marine bivalves. *Journal of Molluscan Studies*, **82**, 282–291.
https://doi.org/10.1093/mollus/eyv056

Le Cam S, Brémaud J, Malkócs T, Kreckelbergh E, Becquet V, Dubillot E, Garcia P, Breton S, Pante E (2023) LAMP-based molecular sexing in a gonochoric marine bivalve (*Macoma balthica rubra*) with divergent sex-specific mitochondrial genomes. *Ecology and Evolution*, **13**, e10320. https://doi.org/10.1002/ece3.10320

Le Cam S (2025) SabLeCam/Cox1_DUI: S. Le Cam et al. 2025 PCJ Evolutionary Biology. Zenodo. https://doi.org/10.5281/zenodo.14849343

Le Mao P, Godet L, Fournier J, Desroy N, Gentil F, Thiebaut E (2020) Mollusques. In: *Atlas de la faune marine invertébrée du golfe Normano-Breton* (ed Roscoff É de la S biologique de)

Leigh J, Bryant D (2015) PopART: Full-feature software for haplotype network construction. *Methods in Ecology and Evolution*, **6**, 1110–1116. https://doi.org/10.1111/2041-210X.12410

Liu H-P, Mitton JB, Wu S-K (1996) Paternal mitochondrial DNA differentiation far exceeds maternal mitochondrial DNA and allozyme differentiation in the freshwater mussel, *Anodonta grandis grandis*. *Evolution*, 952–957. https://doi.org/10.2307/2410870

Lubet P, Devauchelle N, Muzellec ML, Paulet YM, Faveris R, Dao JC (1991) Reproduction of *Pecten maximus* from different fisheries areas Rade de Brest, Baie de St Brieuc, Baie de Seine. In: *Fisheries, Biology and Aquaculture of Pectinids* 8th International Pectinid Workshop., pp. 22–29. Cherbourg (France). https://archimer.ifremer.fr/doc/00049/15984/

Lubośny M, Przyłucka A, Śmietanka B, Burzyński A (2020) *Semimytilus algosus*: first known hermaphroditic mussel with doubly uniparental inheritance of mitochondrial DNA. *Scientific Reports*, **10**, 1–12. https://doi.org/10.1038/s41598-020-67976-6

Luttikhuizen PC, Brenton S, Baker AJ (2003) Disjunct distribution of highly diverged mitochondrial lineage clade and population subdivision in a marine bivalve with pelagic larval dispersal. *Molecular Ecology*, **12**, 2215–2229. https://doi.org/10.1046/j.1365-294X.2003.01872.x

Luttikhuizen PC, Brenton S, Peijnenburg K, Van Der Veer H, Johannesson K (2012) Genetic architecture in a marine hybrid zone: comparing outlier detection and genomic clines analysis in the bivalve *Macoma balthica*. *Molecular Ecology*, **21**, 3048–3061.
https://doi.org/10.1111/j.1365-294X.2012.05586.x

Marko PB, Hoffman JM, Emme SA, McGovern TM, Keever CC, Nicole Cox L (2010) The ‘Expansion–Contraction’ model of Pleistocene biogeography: rocky shores suffer a sea change? *Molecular Ecology*, **19**, 146–169. https://doi.org/10.1111/j.1365-294X.2009.04417.x

Milani L, Ghiselli F, Passamonti M (2012) Sex-linked mitochondrial behavior during early embryo development in *Ruditapes philippinarum* (Bivalvia Veneridae) a species with the doubly
uniparental inheritance (DUI) of mitochondria. *Journal of Experimental Zoology Part B: Molecular and Developmental Evolution*, **318**, 182–189. https://doi.org/10.1002/jez.b.22004

Nguyen L-T, Schmidt HA, Von Haeseler A, Minh BQ (2015) IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. *Molecular biology and evolution*, **32**, 268–274. https://doi.org/10.1093/molbev/msu300

Nikula R, Strelkov P, Väinölä R (2007) Diversity and trans-arctic invasion history of mitochondrial lineages in the North Atlantic *Macoma balthica* complex (Bivalvia: Tellinidae). *Evolution*, **61**, 928–941. https://doi.org/10.1111/j.1558-5646.2007.00066.x

Nikula R, Strelkov P, Väinölä R (2008) A broad transition zone between an inner Baltic hybrid swarm and a pure North Sea subspecies of *Macoma balthica* (Mollusca, Bivalvia). *Molecular Ecology*, **17**, 1505–1522. https://doi.org/10.1111/j.1365-294x.2007.03688.x

Norton P, Spaink G (1973) The earliest occurrence of *Macoma balthica* (L.) as a fossil in the North Sea deposits. *Malacologia: International Journal of Malacology*.

Ort BS, Pogson GH (2007) Molecular population genetics of the male and female mitochondrial DNA molecules of the California sea mussel, *Mytilus californianus*. *Genetics*, **177**, 1087–1099.
https://doi.org/10.1534/genetics.107.072934

Pante E, Becquet V, Viricel A, Garcia P (2019) Investigation of the molecular signatures of selection on ATP synthase genes in the marine bivalve *Limecola balthica*. *Aquatic Living Resources*, **32**, 7. https://doi.org/10.1051/alr/2019001

Pante E, Poitrimol C, Saunier A, Becquet V, Garcia P (2017) Putative sex-linked heteroplasmy in the tellinid bivalve *Limecola balthica* (Linnaeus, 1758). *Journal of Molluscan Studies*, **83**, 226–228. https://doi.org/10.1093/mollus/eyw038

Pante E, Rohfritsch A, Becquet V, Belkhir K, Bierne N, Garcia P (2012) SNP detection from *De Novo* transcriptome sequencing in the bivalve *Macoma balthica*: marker development for evolutionary studies. *PLoS ONE*, **7**, e52302. https://doi.org/10.1371/journal.pone.0052302

Pante E, Simon-Bouhet B (2013) mmapar: A package for importing, plotting and analyzing bathymetric and topographic data in R. *PloS One*, **8**, e73051.
https://doi.org/10.1371/journal.pone.0073051

Paradis E (2010) pegas: an R package for population genetics with an integrated–modular approach. *Bioinformatics*, **26**, 419–420. https://doi.org/10.1093/bioinformatics/btp696

Passamonti M, Boore JL, Scali V (2003) Molecular evolution and recombination in gender-associated mitochondrial DNAs of the Manila clam *Tapes philippinarum*. *Genetics*, **164**, 603–611. https://doi.org/10.1093/genetics/164.2.603

Quéméré E, Baglinière J-L, Roussel J-M, Evanno G, McGinnity P, Launey S (2016) Seascape and its effect on migratory life-history strategy influences gene flow among coastal brown trout (*Salmo trutta*) populations in the English Channel. *Journal of Biogeography*. https://doi.org/10.1111/jbi.12632

R Core Team (2017) *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria.

Ramos-Onsins SE, Rozas J (2002) Statistical properties of new neutrality tests against population growth. *Molecular biology and evolution*, **19**, 2092–2100.
https://doi.org/10.1093/oxfordjournals.molbev.a004034

Rawson PD, Hilbish TJ (1995) Evolutionary relationships among the male and female mitochondrial DNA lineages in the *Mytilus edulis* species complex. *Molecular Biology and Evolution*, **12**, 893–901. https://doi.org/10.1093/oxfordjournals.molbev.a040266

Riginos C, Cunningham CW (2007) Hybridization in postglacial marine habitats. *Molecular Ecology*, **16**, 3971–3972. https://doi.org/10.1111/j.1365-294X.2007.03505.x

Riginos C, Hickerson MJ, Henzler CM, Cunningham CW (2004) Differential patterns of male and female mtDNA exchange across the Atlantic Ocean in the blue mussel, *Mytilus edulis*. *Evolution*, **58**, 2438–2451. https://doi.org/10.1111/j.0014-3820.2004.tb00873.x

Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A (2017) DnaSP 6: DNA sequence polymorphism analysis of large data sets. *Molecular biology and evolution*, **34**, 3299–3302. https://doi.org/10.1093/molbev/msx248
Ruellet T (2013) Contribution à la dynamique de population de *Macoma balthica* en baie de Somme dans le cadre du projet COMORES. Rapport du GEMEL.

Saavedra C, Stewart DT, Stanwood RR, Zouros E (1996) Species-specific segregation of gender-associated mitochondrial DNA types in an area where two mussel species (*Mytilus edulis* and *M. trossulus*) hybridize. *Genetics*, **143**, 1359–1367.
https://doi.org/10.1093/genetics/143.3.1359

Saunier A (2015) Bases génétiques de la différenciation adaptative en milieu anthropisé chez *Macoma balthica*, un bivalve marin à fort flux génique. PhD Thesis Thesis. Université de La Rochelle.

Shimmin LC, Chang BH-J, Li W-H (1993) Male-driven evolution of DNA sequences. *Nature*, **362**, 745–747. https://doi.org/10.1038/362745a0

Skibinski D, Gallagher C, Beynon C (1994) Sex-limited mitochondrial DNA transmission in the marine mussel *Mytilus edulis*. *Genetics*, **138**, 801–809.
https://doi.org/10.1093/genetics/138.3.801

Skibinski OF, Gallaguer C, Quesada H (1999) On the roles of selection, mutation and drift in the evolution of mitochondrial DNA diversity in British *Mytilus edulis* (Mytilidae; Mollusca) populations. *Biological Journal of the Linnean Society*, **68**, 195–213.
https://doi.org/10.1111/j.1095-8312.1999.tb01166.x

Śmietanka B, Burzyński A (2017) Disruption of doubly uniparental inheritance of mitochondrial DNA associated with hybridization area of European *Mytilus edulis* and *Mytilus trossulus* in Norway. *Marine Biology*, **164**, 209. https://doi.org/10.1007/s00227-017-3235-5

Śmietanka B, Burzyński A, Wenne R (2009) Molecular population genetics of male and female mitochondrial genomes in European mussels *Mytilus*. *Marine Biology*, **156**, 913–925.
https://doi.org/10.1007/s00227-009-1137-x

Śmietanka B, Zbawicka M, Sanko T, Wenne R, Burzynski A (2013) Molecular population genetics of male and female mitochondrial genomes in subarctic *Mytilus trossulus*. *Marine Biology*, **160**, 1709–1721. https://doi.org/10.1007/s00227-013-2223-7

Stewart DT, Kenchington ER, Singh RK, Zouros E (1996) Degree of selective constraint as an explanation of the different rates of evolution of gender-specific mitochondrial DNA lineages in the mussel *Mytilus*. *Genetics*, **143**, 1349–1357. https://doi.org/10.1093/genetics/143.3.1349

Stewart DT, Saavedra C, Stanwood RR, Ball AO, Zouros E (1995) Male and female mitochondrial DNA lineages in the blue mussel (*Mytilus edulis*) species group. *Molecular Biology and Evolution*, **12**, 735–747. https://doi.org/10.1093/oxfordjournals.molbev.a040252

Strelkov P, Nikula R, Väinölä R (2007) *Macoma balthica* in the White and Barents Seas: properties of a widespread marine hybrid swarm (Mollusca: Bivalvia). *Molecular Ecology*, **16**, 4110–4127.
https://doi.org/10.1111/j.1365-294x.2007.03463.x

Stuckas H, Stoof K, Quesada H, Tiedemann R (2009) Evolutionary implications of discordant clines across the Baltic *Mytilus* hybrid zone (*Mytilus edulis* and *Mytilus trossulus*). *Heredity*, **103**, 146–156. https://doi.org/10.1038/hdy.2009.37

Tajima F (1983) Evolutionary relationship of DNA sequences in finite populations. *Genetics*, **105**, 437–460. https://doi.org/10.1093/genetics/105.2.437

Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. *Genetics*, **123**, 585–595. https://doi.org/10.1093/genetics/123.3.585

Tamura K, Kumar S (2002) Evolutionary distance estimation under heterogeneous substitution pattern among lineages. *Molecular biology and evolution*, **19**, 1727–1736.
https://doi.org/10.1093/oxfordjournals.molbev.a003995

Tamura K, Nei M (1993) Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. *Molecular biology and evolution*, **10**, 512–526. https://doi.org/10.1093/oxfordjournals.molbev.a040023

Tamura K, Stecher G, Kumar S (2021) MEGA11: molecular evolutionary genetics analysis version 11. *Molecular biology and evolution*, **38**, 3022–3027. https://doi.org/10.1093/molbev/msab120

Toews DPL, Brelsford A (2012) The biogeography of mitochondrial and nuclear discordance in animals. *Molecular Ecology*, **21**, 3907–3930. https://doi.org/10.1111/j.1365-294x.2012.05664.x
Trejo-Salazar R-E, Castellanos-Morales G, Hernández-Rosales D, Gámez N, Gasca-Pineda J, Garza MRM, Medellin R, Eguiarte LE (2021) Discordance in maternal and paternal genetic markers in lesser long-nosed bat *Leptonycteris yerbabuenae*, a migratory bat: recent expansion to the North and male phylopatry. *PeerJ*, **9**, e12168. https://doi.org/10.7717/peerj.12168

Turelli M, Moyle LC (2007) Asymmetric Postmating Isolation: Darwin's Corollary to Haldane's Rule. *Genetics*, **176**, 1059-1088. https://doi.org/10.1534/genetics.106.065979

Väinölä R (2003) Repeated trans-Arctic invasions in littoral bivalves: molecular zoogeography of the *Macoma balthica* complex. *Marine Biology*, **143**, 935-946. https://doi.org/10.1007/s00227-003-1137-1

Wares JP, Cunningham CW (2001) Phylogeography and historical ecology of the North Atlantic intertidal. *Evolution*, **55**, 2455-2469. https://doi.org/10.1111/j.0014-3820.2001.tb00760.x

Watterson G (1975) On the number of segregating sites in genetical models without recombination. *Theoretical population biology*, **7**, 256-276. https://doi.org/10.1016/0040-5809(75)90020-9

Yurchenko AA, Katolikova N, Polev D, Shcherbakova I, Strelkov P (2018) Transcriptome of the bivalve *Limecola balthica* L. from Western Pacific: A new resource for studies of European populations. *Marine Genomics*, **40**, 58-63. https://doi.org/10.1016/j.margen.2018.03.007

Zouros E (2000) The exceptional mitochondrial DNA system of the mussel family Mytilidae. *Genes & Genetic Systems*, **75**, 313-318. https://doi.org/10.1266/ggs.75.313

Zouros E (2013) Biparental inheritance through uniparental transmission: the doubly uniparental inheritance (DUI) of mitochondrial DNA. *Evolutionary Biology*, **40**, 1-31. https://doi.org/10.1007/s11692-012-9195-2

Zouros E, Ball AO, Saavedra C, Freeman KR (1994) An unusual type of mitochondrial DNA inheritance in the blue mussel *Mytilus*. *Proceedings of the National Academy of Sciences*, **91**, 7463-7467.