Trematode genetic patterns at host individual and population scales provide insights about infection mechanisms

Abstract

Multiple parasites can infect a single host, creating a dynamic environment where each parasite must compete over host resources. Such interactions can cause greater harm to the host than single infections and can also have negative consequences for the parasites themselves. In their first intermediate hosts, trematodes multiply asexually and can eventually reach up to 20% of the host's biomass. In most species, it is unclear whether this biomass results from a single infection or co-infection by 2 or more infective stages (miracidia), the latter being more likely a priori in areas where prevalence of infection is high. Using as model system the trematode Bucephalus minimus and its first intermediate host cockles, we examined the genetic diversity of the cytochrome c oxidase subunit I region in B. minimus from 3 distinct geographical areas and performed a phylogeographic study of B. minimus populations along the Northeast Atlantic coast. Within localities, the high genetic variability found across trematodes infecting different individual cockles, compared to the absence of variability within the same host, suggests that infections could be generally originating from a single miracidium. On a large spatial scale, we uncovered significant population structure of B. minimus, specifically between the north and south of Bay of Biscay. Although other explanations are possible, we suggest this pattern may be driven by the population structure of the final host.

Introduction

With about 45 000 species described in a wide range of ecosystems, trematodes are one of the most common and widespread group of parasites (Carlson et al., 2020). They can be found at almost across all trophic levels of dynamic food chains (Bartoli and Gibson, 2007). Trematodes are an important component of ecosystem biodiversity with significant impacts at the host individuals (SchulteOehlmann et al., 1997; Curtis et al., 2000; Thieltges, 2004), host populations (Fredensborg et al., 2005) and ecosystem communities (Poulin, 1999; Mouritsen and Poulin, 2002; Goedknegt et al., 2016). Besides, by contributing to the nutrient cycle, acting as indicators of environmental changes or as proxy of environmental diversity (due to multi-host life cycles), trematode presence may indicate a healthy and resilient ecosystem (Johnson et al., 2010; Hatcher and Dunn, 2011).

Trematodes have a complex life cycle that alternates between free-living and parasitic stages. The miracidium, the trematode larva hatched from an egg, infects the first intermediate host (usually a mollusc) and transforms into sporocysts or rediae (parasitic stage). At this stage, sporocysts or rediae, through asexual multiplication, produce cercariae (free-living stage) that emerge from the first host to infect the second intermediate host (a vertebrate or invertebrate) where they settle as metacercariae. In the trematode's final host (a vertebrate), after ingestion of the second host, metacercariae develop into adult flukes, reproduce sexually, and complete the life cycle (Cribb et al., 2003; Bartoli and Gibson, 2007). In the first intermediate host, sporocyst stages are overtly destructive, replacing host tissue and reaching up to 20% of the host's biomass (Dubois et al., 2009; Preston et al., 2013), with direct consequences for host reproduction (Carballal et al., 2001), growth (Bowers, 1969) and energy demand (Jokela et al., 1993), leading to eventual host death (Thieltges, 2006). In most species, it is currently unknown whether this sporocyst biomass results from a single miracidium, that excludes other miracidia by predation or intraspecific competition, or from co-infection. If co-infection is the rule, trematode invasion may result in a burden that the host might not be able to bear (Fredensborg and Poulin, 2005; Mideo, 2009). On the other hand, co-infection can, occasionally, benefit the host by lessening the overall burden of infection, by reducing parasite infection success, or by strengthening the host's immune response, leading to higher resistance to infection (Dumont et al., 2007; Balmer et al., 2009).

Determining the genetic diversity of trematode sporocysts within and among individual first intermediate hosts is therefore important to understand the biology, behaviour and evolutionary patterns of these parasites. Nonetheless, current knowledge regarding biology of these parasites, and particularly about the sporocyst life stage, is still very scarce, despite trematodes' wide distribution and importance to the ecosystem. Few studies have focused on the conspecific diversity of trematode sporocysts within the first intermediate host (Rauch et al., 2005; Keeney et al., 2007; Lagrue et al., 2007), showing that the likelihood of infection by conspecifics increased with the prevalence of the parasite in the community (Keeney et al., 2008; Louhi et al., 2013). Moreover, at the population level, the host has a significant impact on genetic diversity and population structure of parasites, with substantial gene flow occurring in parasite species with efficient dispersion mechanisms (Agola et al., 2009; Feis et al., 2015). However, most studies of trematode genetic diversity have focused on intermediate and final hosts, with particular emphasis on trematodes with harmful impacts on human (Theron et al., 2004; Bell et al., 2006; Balmer et al., 2009) or socio-economically important species, namely fish (Vilas et al., 2003; Criscione and Blouin, 2006).

Bucephalus minimus is a marine trematode parasite that occurs in several aquatic systems along the Northeast Atlantic coast and Mediterranean Sea (Magalhães et al., 2015). In the Atlantic area, this parasite infects the European edible cockle, Cerastoderma edule, which serves as the first intermediate host when a miracidium penetrates its tissue. In cockles, the prevalence of this parasite varies greatly among coastal systems and season; depending on the time since infection, the parasite's dry mass in infected cockles can range from 1 to 20% of the total living tissue within the cockle shell (Magalhães et al., 2015; de Montaudouin et al., 2021). Bucephalus minimus initially infects the cockle's gonad and digestive gland but promptly spreads to other parts of the host, eventually invading the entire body (Desclaux et al., 2002; de Montaudouin et al., 2009). Infection by B. minimus results in castration (Carballal et al., 2001) and energy consumption (Dubois et al., 2009), leading to starvation and autolysis of the host's digestive tract. This parasite is considered as one of the most harmful trematode parasites infecting C. edule (Magalhães et al., 2015; de Montaudouin et al., 2021). Inside cockles, B. minimus produces sporocysts and cercariae through asexual multiplication, which emerge and infect the goby Pomatoschistus spp. (second intermediate host), where they encyst and develop into metacercariae. The final host, Dicentrarchus labrax, the European seabass, is infected after consumption of parasitized gobies. Metacercariae develop into adult flukes and produce eggs, through sexual reproduction, to complete the cycle (Pina et al., 2009; Magalhães et al., 2015).

Due to the limited information regarding this trematode's sporocyst stages, the primary goal of the present study was to assess the genetic variability of the cytochrome c oxidase subunit I (COI) region of B. minimus sporocyst DNA within and among the first intermediate host, C. edule. We tested the hypothesis that higher genetic variability at the host individual scale (resulting in 2 or more haplotypes among sporocysts in a single cockle) are more common in localities with high prevalence of infection, where joint infections should be more frequent by chance alone. As our second goal, a phylogeographic study of B. minimus was also carried out combining information available in the literature and from samples taken on cockle beds that were examined in this study for the first time (i.e. Iberian Peninsula and Great Britain).

Materials and methods

Bucephalus minimus samples

The genetic variability of B. minimus sporocyst haplotypes within the same host was studied by collecting specimens present in the first intermediate host, the edible cockle, at 3 different beds with different prevalence along the European Atlantic coast: Ria de Aveiro, Aveiro, Portugal (lowest prevalence [Magalhães et al., 2018]); de la Ramallosa Lagoon, Baiona, Spain (moderate prevalence [Intecmar, 2021]); and Île aux Oiseaux, Arcachon, France (with high levels of prevalence [Magalhães et al., 2015]) (Fig. 1). Adult cockles (between 20 and 30 mm shell length) were haphazardly collected at low tide and dissected in the laboratory to morphologically identify B. minimus infection. The flesh was then transferred and observed under a stereomicroscope by carefully compressing it between 2 sterilized glass slides. Four sporocyst replicates were extracted per cockle (in a total of 5 cockles per sampling site) using forceps, and preserved separately in 100% ethanol at −20°C. All material was sterilized between samples. To enhance the possibility of different clones, sporocysts were taken from different infected tissues of the same cockle (i.e. 1 from the foot, 1 from the gills and 2 from the digestive gland).

Figure 1. Geographical location of the Cerastoderma edule cockle beds sampled for the study of Bucephalus minimus genetic variability at host level (Aveiro, Baiona and Arcachon [in italics]) and for phylogeographic analysis. *Cockle beds sampled for the first time in this study.

DNA isolation, amplification and sequencing

Genomic DNA extraction from B. minimus specimens was performed using E.Z.N.A Mollusc DNA kit (Omega Bio-Tek, Norcross, GA, USA) in accordance with the manufacturer's instructions. Nanodrop was used to assess DNA concentration, and, if needed, aliquots were created to dilute DNA to approximately 30 ng μL⁻¹.

The mitochondrial COI fragment was amplified using the MplatCOX1-dF (5′-TTW CIT TRG ATC ATA AG-3′) and MplatCOX1-dR (5′-TGA AAY AAY AII GGA TCI CCA CC-3′) primers (Moszczynska et al., 2009), resulting in sequences of 587 bp. The polymerase chain reaction (PCR) was carried out in a final volume of 20 μL composed of 1X reaction buffer, 2.5 mm MgCl₂, 100 μm deoxynucleotide triphosphates (dNTPs), 0.5 μm of forward and reverse primers, 0.65 units of ThermoFisher AmpliTaq Gold DNA polymerase (ThermoFisher, Waltham, MA, USA) and 60 ng of DNA. The PCR programme employed had an initial denaturation step for 10 min at 95°C, followed by 35 cycles of 30 s at 94°C, 30 s at 50°C and 1 min at 72°C and a final extension of 10 min at 72°C. Following this initial PCR, a second PCR was carried out using the same conditions previously described but using 2 μL of PCR product instead of DNA. The amplified PCR products were analysed by electrophoresis through a 1% agarose gel dyed with SYBR safe DNA Gel Stain (Invitrogen, Carlsbad, CA, USA) and visualized under UV light.

PCR products were enzymatically purified with ExoSAP mix (10 μL of PCR product, 0.6 units of EXO I [DNA nuclease] and 0.3 units of shrimp alkaline phosphatase [SAP] for a final volume of 12 μL) under the following conditions: 60 min at 37°C and 15 min at 85°C. Purified PCR products were sequenced using the ABI Prism BigDye^TM Terminator v3.1 Cycle Sequencing Kit protocol on an ABI Prism 3730 xl automatic sequencer (Applied Biosystems, Foster City, CA, USA). All sequences obtained in this study were deposited in GenBank (accession numbers: OQ625925–OQ625936; Table 1). Variable sites were manually checked using the SEQSCAPE 2.5 program (Applied Biosystems) and aligned using ClustalW algorithm implemented on BioEdit v.7.2.5 (Hall, 1999) with LaB haplotype (GenBank accession number: KF880429.1) as reference. Identification of the different haplotypes in the B. minimus specimens analysed was carried out using the software DNAsp v5.10 (Librado and Rozas, 2009). Finally, to check for the presence of premature stop codons, COI haplotypes identified were translated to amino acid sequences using the flatworm mtDNA code in the online software EMBOSS Transeq (Rice et al., 2000; Goujon et al., 2010).

Table 1. Accession number for B. minimus COI gene DNA sequences downloaded (with reference) and deposited (in bold) in GenBank

Data analysis

Bucephalus minimus genetic variability at host level

To determine the genetic variability of B. minimus specimens within the same cockle and among different cockles, haplotypes found in the different tissues of each analysed cockle were compared using BioEdit. A haplotype network was built using the obtained data, identifying haplotypes for each cockle with a different colour. The haplotype network was constructed by calculating the distance (based on number of base pair differences) between DNA sequences and determining the number of mutations between haplotypes using ‘pegas’ and ‘ape’ packages of R Statistical Software v.4.2.2 (Paradis, 2010; Paradis and Schliep, 2019).

Phylogeographic analysis

For phylogeographic analyses, together with the specimens collected in the present study, DNA extractions of 13 specimens from 7 cockle beds sampled as part of the COCKLES Interreg project (http://cockles-project.eu/) were sequenced as previously described. Samples were available for Aveiro, Portugal (1 sample), Noia, Spain (1 sample), Arcachon, France (1 sample), Bay of Somme, France (4 samples), The Dee, Wales (2 samples), Burry Inlet, Wales (2 samples) and Wadden Sea, the Netherlands (2 samples). Additionally, the analysis included another 54 COI sequences available in the GenBank database retrieved in January 2023 (see Feis et al., 2015 and Table 1). In total, sequences represented specimens from 11 different cockle beds located in 7 countries, covering a large part of the natural distributional range of cockles in the Atlantic area (Fig. 1).

Phylogenetic relationships were studied using different and complementary approaches. First, a haplotype network was computed for the full dataset as previously described, which included information regarding the cockle bed in which the haplotypes were found and the frequency of occurrence. Identification of unique haplotypes present in the dataset was carried out with DnaSP v.5.10 (Librado and Rozas, 2009). Using MEGA X software (Kumar et al., 2018), the Hasegawa–Kishino–Yano nucleotide substitution rate (HKY) with a γ value of 0.655 and invariable sites of 0.728 was identified as the most probable nucleotide substitution model for our data. Phylogenetic trees were constructed using exclusively the different haplotypes identified and the nucleotide substitution model described above. Maximum likelihood (ML) and neighbour-joining rooted and unrooted trees were constructed using the R Statistical Software v.4.2.2. The rooted tree was created using the COI sequences of Rhipidocotyle sp. (a trematode from the same family as B. minimus, Bucephalidae, GenBank accession number: KM538111.1) and Himasthla quissetensis (a trematode that infects cockles as second intermediate host but from a different family, Himasthlidae; GenBank accession number: MN272732.1) as outgroups. The sequences were trimmed to 540 bp to remove missing data. The unrooted trees were constructed using the full 587 bp sequences of B. minimus. The robustness of the branches for the phylogenetic trees was estimated with 1000 bootstrap replicates and a likelihood ratio test was performed based on the minimum Akaike information criterion values for ML. All phylogeographic analyses were performed with the ‘ape’, ‘pegas’, ‘ggtree’ and ‘phangorn’ packages of R Statistical Software v.4.2.2 (Paradis, 2010; Schliep, 2011; Yu et al., 2017; Paradis and Schliep, 2019).

Genetic diversity parameters, calculated as haplotype diversity (h) and nucleotide diversity (π), were estimated within each cockle bed studied using Arlequin v.3.5.1.3 (Excoffier et al., 2005). The HKY model is not available in this software. For this reason, the Tamura–Nei (TN) model with a γ value of 0.639 (similar to the HKY and identified as the third best option for MEGA X) was used for π estimations. Genetic structure and population differentiation were assessed with global and pairwise coefficients of population differentiation applying the TN with γ value of 0.639 substitution rate (ϕ_ST values) with Arlequin. Analysis of molecular variance (AMOVA) applying different models of a priori clustering (based on host population genetics [Souche et al., 2015; Vera et al., 2022] and observed data – see results) was carried out to study the distribution of genetic variation within (ϕ_SC) and among (ϕ_CT) bed groups using Arlequin. The significance for all the ϕ statistics was evaluated with 10 000 permutations.

Results

Bucephalus minimus genetic variability at host level

During this study, a total of 210 cockles were analysed, of which 17 were found to be infected with B. minimus (5 each in Aveiro and Baiona and 7 in Arcachon). The prevalence of B. minimus in the cockle beds sampled varied from 3.3% in Aveiro (Portugal) to 23.3% in Arcachon (France). In Baiona (Spain), B. minimus was present in 16.7% of the sampled cockles. Five infected cockles per bed were used to extract 4 sporocysts per cockle, yielding a total of 60 samples. Fifty-six samples were successfully sequenced, while 4 samples, from a single cockle from Arcachon (France), were not successfully sequenced due to DNA extraction problems. From the 56 sequenced sporocysts, belonging to 14 infected cockles, 12 different haplotypes were identified, with 5, 3 and 4 haplotypes found in Aveiro, Baiona and Arcachon, respectively (Table 2). Six of the identified haplotypes were characterized for the first time (named as BmA–BmF; see Table 1).

Table 2. Coordinates of each cockle bed, the number of analysed cockles (N _cockles), number of parasites sequenced (N_{B. minimus}), prevalence of B. minimus, number of haplotypes (k), number of polymorphic sites (PS) and haplotype composition (between parentheses the number of individuals bearing the same haplotype when different from one)

All sporocysts of B. minimus from the same cockle had identical haplotype, however B. minimus haplotypes identified in different cockles from the same bed were different, except in Baiona where the haplotype found in 3 different cockles was identical (haplotype LaE; see Table 2 and Fig. 2). Moreover, haplotypes were not shared among the 3 beds (Fig. 2).

Figure 2. Haplotype network of Bucephalus minimus samples from Aveiro (Portugal), Baiona (Spain) and Arcachon (France) based on genetic distance (number of base pair differences) of cytochrome c oxidase subunit 1 (COI) gene sequences. Different haplotypes with respective names are represented by circles, with circle size proportional to observed frequency. Inferred mutation steps are shown by black dots. Colours depict samples taken from the same cockle.

Phylogeographic analysis

From the 69 DNA sequenced samples (56 for the B. minimus haplotype genetic variability at host level study, and 13 from different European sites selected from the COCKLES project), 12 resulted in novel haplotypes (Table 1). No premature STOP codons were identified in these sequences (data not shown). Thus, when GenBank resources were included, a total of 162 COI gene sequences of B. minimus specimens from 11 cockle beds were analysed. From these available sequences, 66 represented unique haplotype sequences. Shared haplotypes (i.e. those found in more than 1 cockle bed) accounted for 17% of the total. The LaAQ haplotype was the most prevalent and abundant haplotype, occurring 42 times across 6 different beds. On the other hand, 83% of the haplotypes were exclusively found in a single bed, with several reported only once (i.e. singletons). Arcachon presented the highest number of different B. minimus haplotypes detected (19, see Table 3).

Table 3. Haplotype frequency for B. minimus COI gene per cockle bed including haplotype diversity (h ± standard deviation) and nucleotide diversity (π ± standard deviation)

Due to the high variability found in the COI region, the analysis of B. minimus haplotypes across the various cockle beds produced a complex network made up of several closely connected haplotypes and associated mutational steps, with no more than 5 mutations separating any 2 successive haplotypes identified. A common haplotype (LaAQ), observed in several beds situated north of Arcachon (44°N), was located in the centre of the network, from which numerous other haplotypes diverged in a star-like pattern. These haplotypes were exclusively found in a single bed or shared between relatively close beds (Fig. 3). Nonetheless, haplotype clusters (i.e. haplogroups) were identified in specific geographic areas, and beds from the South (Merja Zerga, Aveiro, Baiona and Noia) and North (Bay of Somme, English Channel, Celtic Sea, Burry Inlet, The Dee and Wadden Sea) did not share any haplotype, with the exception of the LaE haplotype identified in Baiona, Arcachon and the English Channel. Phylogenetic relationships observed in the network were also confirmed with the phylogenetic trees (Fig. 4).

Figure 3. Haplotype network of Bucephalus minimus based on genetic distance (number of base pair differences) of cytochrome c oxidase subunit 1 (COI) gene sequences across 11 beds of Cerastoderma edule where Bucephalus minimus specimens were collected. Different haplotypes with respective names are represented by circles, with circle size proportional to observed frequency. Inferred mutation steps are shown by black dots. Colours depict sample location.

Figure 4. Maximum likelihood phylogenetic tree based on cytochrome c oxidase subunit I (COI) gene sequences of Bucephalus minimus haplotypes from this study and retrieved from GenBank database. Rhipodocotyle sp. and Himasthla quissetensis haplotypes were used as outgroups. Numbers at nodes represent the percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) using neighbour-joining and maximum likelihood methods (NJ/ML). The scale bar indicates the distances in substitutions per nucleotide.

The Arcachon bed, located in the centre of cockle's distributional range between the northern and southern geographic areas, exhibited the highest number of detected haplotypes (19), sharing haplotypes with locations from both regions (Table 3). Excluding Noia, where only 1 individual was analysed, haplotype diversity ranged from 0.3846 in Merja Zerga to 1.0000 in beds where all individuals analysed had a distinct haplotype (Aveiro, Burry Inlet and The Dee). Nucleotide diversity ranged from 0.0007 in Merja Zerga to 0.0076 in Baiona (Table 3). The high diversity and the haplotype distribution among locations were also reflected in the ϕ_ST values. Global ϕ_ST for the whole region was 0.2922 (P value < 0.001). Many pairwise ϕ_ST values resulted in significant differences, although many comparisons between close locations were non-significant, mainly among those involving northern beds (Supplementary Table S1). Moreover, all pairwise ϕ_ST values involving Merja Zerga were high and highly significant (P value < 0.001), suggesting the singularity of this bed (global ϕ_ST in the whole region excluding Merja Zerga = 0.1706, P value < 0.001). These results suggest the presence of one northern group (composed by Bay of Somme, English Channel, Celtic Sea, Burry Inlet, The Dee and Wadden Sea) more homogeneous genetically (ϕ_ST = 0.0295, P value = 0.073) than the southern one (composed by Merja Zerga, Aveiro, Baiona and Noia; ϕ_ST = 0.6694, P value < 0.001), with Arcachon representing a potential contact region between both geographic areas. The ϕ_ST value in the northern group increased up to 0.0336 (P value = 0.018) when Arcachon was included, while this value decreased in the southern group when this location was included although it remained quite high (ϕ_ST = 0.4184, P value < 0.001). Hence, these results suggest a closer relationship of Arcachon with the northern group. AMOVA analysis assigned 34.58% of the genetic differentiation to differences between northern and southern groups (ϕ_CT = 0.3458, P value = 0.008), this percentage being three times higher than those assigned to differences among beds within groups (ϕ_SC = 0.2109, P value < 0.001, percentage of genetic differentiation = 13.80%). The AMOVA model including Arcachon in the northern group yielded similar values (ϕ_CT = 0.3100, P value = 0.006, percentage of genetic differentiation = 31.00%; ϕ_SC = 0.1567, P value < 0.001, percentage of genetic differentiation = 10.81%). This model assigned a higher percentage of genetic differentiation among groups and a lower percentage to differences among beds within groups than the model including Arcachon in the southern group (ϕ_CT = 0.1327, P value = 0.048, percentage of genetic differentiation = 13.27%; ϕ_SC = 0.2299, P value < 0.001, percentage of genetic differentiation = 19.13%), suggesting a more coherent grouping of the beds in the former model.

Discussion

Bucephalus minimus genetic variability at host level

Co-infection by multiple parasites, from the same or different species, within the same host is a well-recognized phenomenon in the parasitological literature (Poulin, 2001; Read and Taylor, 2001). This pattern has been extensively studied for several parasite species, namely with an impact on human health (Theron et al., 2004; Bell et al., 2006). For example, in the case of malaria, more than 5 strains have been found to be infecting the same host (Bell et al., 2006). The same trend was observed for the trematode parasite Schistosoma mansoni within their second and final host (Theron et al., 2004). Similar to what is observed for metacercarial or adult stages of trematode parasites, it would be anticipated that different clones would infect the same first intermediate host when thousands of eggs per infected definitive host are shed into the water column, i.e. thousands of miracidia hatching within metres of each other. This was observed for some trematode species (Rauch et al., 2005; Keeney et al., 2007; Lagrue et al., 2007). In the present study, only 1 COI haplotype was found inside each infected cockle (regardless of the samples' origin), in contrast with what has been previously recorded. Nevertheless, it should be noted that in the present study only the COI region (maternally inherited) was sequenced, while for previous studies, microsatellite markers were used to identify individual variability. In fact, microsatellite markers are more accurate for population structure analysis and individual identification since they are highly variable polymorphic regions (Abdul-Muneer, 2014) and inherited from both parents. Despite its limitations, the high genetic variability of the COI region found in this species (12 haplotypes out of 14 analysed cockles), as well as in each of the studied beds individually (no repeated haplotypes in 2 out of 3 studied cockle beds – see Table 2), suggests that possibly only 1 individual (i.e. 1 miracidium) of B. minimus infects the host and/or prevails inside it. Unfortunately, no microsatellite markers are currently developed for this trematode species or other closely related species (which could have been used by cross-validation). Therefore, further studies using nuclear markers that are either highly polymorphic (such as microsatellites) or in a high number (single nucleotide polymorphisms, SNPs) – which would reduce the probability of random sharing of multilocus genotypes among individuals – will be necessary to confirm our results.

In the Ria de Aveiro, given the low prevalence of B. minimus, the presence of only 1 parasite haplotype per host was not surprising. There is a well-known upwelling mechanism offshore of this coastal lagoon (Queiroz et al., 2012), resulting in low water temperature and consequently lower prevalence and abundance of trematode parasites compared to other coastal systems where cockles are distributed (Correia et al., 2020). Trematodes are highly sensitive to temperature, both in their free-living and parasitic stages (Thieltges and Rick, 2006; Selbach and Poulin, 2020). For example, the production and hatching rate of trematode eggs are positively correlated with temperature, peaking under ideal thermal conditions (Morley, 2012; Morley and Lewis, 2017). The same happens with cercarial multiplication within and emergence from the first intermediate host (Poulin, 2006; de Montaudouin et al., 2016). Hence, Ria de Aveiro may have fewer free-living stages (miracidia) in the water and might take longer for the cycle to complete. Adding to it the high density of the host found in the area, a parasite intensity dilution effect (as occurs in other regions, e.g. Magalhães et al., 2017) might also contribute to a single conspecific parasite infection in each individual host, as we observed. However, the same pattern (i.e. 1 haplotype per host) was found for cockles from de la Ramallosa lagoon (Baiona) and Île aux Oiseaux (Arcachon), where the prevalence of B. minimus can exceed 20%. Double infection by trematode sporocysts in cockles is rare (Magalhães et al., 2015, 2020), most likely due to rare exposure of the host to a second miracidium. However, occurrence of co-clone infection has been shown to arise when the prevalence of hosts with trematode sporocyst infection rises (Keeney et al., 2008; Louhi et al., 2013). Therefore, the rationale above that relates temperature and low prevalence as causes of single haplotype infection lacks support in Baiona (Spain) and Arcachon (France) and suggests again that B. minimus infection may originate from a single miracidium.

Alternatively, the genetic diversity at host level may be determined by the infection mechanisms of these parasites, such as intraspecific competition. The presence of a single B. minimus haplotype per host could also be explained by the production of substances that could change host chemical attractiveness (Baiocchi et al., 2017) or be toxic against competitors (Burman, 1982; Selva et al., 2009). This is true for nematode parasites, but there is no information on trematodes. However, synthesis of harmful chemicals seems highly unlikely as it would impact the trematode's own clones. Nonetheless, to test these predictions, specific experiments would need to be conducted.

Another explanation may be that cockles with multiple infection (i.e. those with more than 1 haplotype) are rare because they incur higher mortality rates than single infections. This stage of the trematode life cycle is highly deleterious for the host and can lead to mass mortality events during outbreaks (Thieltges, 2006; de Montaudouin et al., 2021). Periods of high prevalence may also result in a rise in co-infections, which may heighten the host's susceptibility and mortality. Theoretical and empirical studies provide conflicting findings, and it is therefore unclear whether susceptibility and subsequent host death are directly connected to co-infection by conspecific parasites (Read and Taylor, 2001; Davies et al., 2002; Alizon and van Baalen, 2008). Regardless, a study conducted using different strains of the trematode S. mansoni in their host showed an increase in overall pathogenicity (Davies et al., 2002).

As an alternative and excluding the above scenarios, it is possible that after parasite settlement, co-infection becomes unlikely due to alterations in host characteristics. For instance, the nematode Acanthocheilonema viteae induces an immune response by the host to further infection by free-living larvae, with no effect on the nematode's persistence within the host (Rajakumar et al., 2006). The extensive use of host tissues, particularly the gonad and digestive gland (Dubois et al., 2009), could be another reason for the absence of co-infection. This exhaustive use of resources might prevent the settlement of other parasite individuals. Interestingly, co-infection with another trematode species at the same sporocyst stage is exceedingly rare and less frequent than expected (in terms of probability), without predominance by either species (Magalhães et al., 2020). The exhaustive use of tissue preventing the establishment of further miracidia seems to be the most likely explanation for our findings. Thus, the present results cannot fully demonstrate whether co-infection by conspecifics is present or not due to the limited sample size and the genetic marker used. They do, however, suggest that co-infections might be rare and not the rule, and that, for a co-infection to occur, different miracidia must infect simultaneously or within a short period of time. Further studies using codominant nuclear DNA molecular markers, such as microsatellites or SNPs (unavailable to date for B. minimus), or attempts at experimental infection of previously infected cockles, would help to determine whether co-infections are possible and confirm our results.

Bucephalus minimus phylogeography

The COI sequences amplified in this study, as well as those available from the literature (Feis et al., 2015), were all grouped together, with short branches separating each sequence, indicating that the sequences all belong to the same species throughout the entire cockle distributional range, with no evidence of any subspecies.

The haplotype network constructed resulted in a complex pattern with several different haplotypes, generally with few mutational steps between them (at most 5 steps between 2 adjacent haplotypes) in accordance with the previous study by Feis et al. (2015). These results are typical of a population with relatively stable age structure and a fixed ratio of natural growth (i.e. stable demography) (Loewe and Hill, 2010). The absence of shared haplotypes among the various beds (most haplotypes were present in a single bed) supports this conclusion. These results were expected and consistent with previous reports on this species (Feis et al., 2015). Trematodes have a complex life cycle that involves a vertebrate, usually a fish or a bird, as definitive host (Cribb et al., 2003). In the case of B. minimus, its final host is the European seabass, which has limited migratory capacity compared to birds (de Pontual et al., 2023). This limited mobility results in a higher degree of isolation between populations of B. minimus.

In our study, 2 haplogroups were identified, one located south of the Bay of Biscay (including Merja Zerga, Aveiro, Baiona and Noia) and the other to the north (including Bay of Somme, English Channel, Celtic Sea, Burry Inlet, The Dee and Wadden Sea). These groups did not share haplotypes except one (LaE). However, this division also coincides with the cockle beds where a lower number of individuals were assessed (see Table 3). Besides, it is noteworthy that there is a core haplotype present in 6 different beds from which most of the other haplotypes diverge. This haplotype has only been identified in the northernmost beds. However, it may represent an ancestral haplotype that has spread across the different regions, and then undergone evolutionary mutations that are specific to each bed in which it occurs. For example, the haplotype network created for the B. minimus genetic variability at host level study (Fig. 2) revealed higher similarity of Aveiro (south group) and Arcachon (central/north group) haplotypes than Aveiro and Baiona (south group) haplotypes due to a parsimonious position (position 50 of our alignment) in the COI sequences where a cytosine replaces a thymine in all Baiona haplotypes. Consequently, the presence or absence of haplotypes in a particular bed should be interpreted cautiously, especially for southern beds.

Despite the limitations previously described, the population genetics of B. minimus exhibited similar geographic clusters as those of its first intermediate host, C. edule [structured as northwards and southwards of French Brittany (Vera et al., 2022)]. Cockle dispersal mainly occurs during their larval stages (Martel and Chia, 1991), when B. minimus cannot infect them (Magalhães et al., 2015). It may also occur through human-mediated movements, although no information is available regarding this possibility. Therefore, it is more logical to expect a correlation with the population genetic structure of seabass, the most mobile host in the trematode life cycle (Zemmer et al., 2020). Two distinct genetic population units, one in the Atlantic area and another in the Mediterranean, have been previously identified for seabass (Souche et al., 2015), with a slight genetic differentiation in the Atlantic area observed south of the Strait of Gibraltar (Morocco) attributed to a hybrid zone between 2 evolutionary lineages (Lemaire et al., 2005; Vandeputte et al., 2019). A genomic survey conducted with more than 2700 molecular markers (i.e. SNPs) throughout the Atlantic area revealed the presence of 3 different groups weakly differentiated and geographically distributed. Thus, all Atlantic wild fish belonged to a single group, except specimens from the northern North Sea (i.e. Norway) and the Strait of Gibraltar (AQUATRACE, 2017), matching the results of Souche et al. (2015). This structure (see Fig. S1) could explain the high B. minimus population differentiation found when Merja Zerga (located on the Moroccan coast) was compared with the remaining beds, although more beds south of the Strait of Gibraltar need to be studied, along with greater sample sizes from southern beds, especially in the Iberian Peninsula, to corroborate this explanation. Following seabass genetic structure, the other cockle beds analysed may be included within the same genetic group. Despite this homogeneity, Quéré et al. (2010) identified genetic differentiation between seabass from the Bay of Biscay and North Sea using a molecular marker (microsatellite) under selection associated to the somatolactin gene, although this pattern was not confirmed either by other candidate genes or adaptive variation screenings (Souche et al., 2015; AQUATRACE, 2017). Moreover, an electronic tagging study in seabass has identified the Bay of Biscay as a potential hybridization zone for several subpopulations of seabass from various Atlantic coast locations (de Pontual et al., 2023), which could account for the region's high haplotype diversity and differentiation in the population genetics of B. minimus. In any event, the effect of selective processes on seabass population structure cannot be ruled out and could explain the current genetic structure in B. minimus, with the Bay of Biscay (here represented by Arcachon bed) being a possible hybrid region between northern and southern groups. Thus, the genetic structure among B. minimus population uncovered in our study may reflect local adaptation of the parasite to the most common host genotypes occurring in each of the regions sampled (e.g. Sasal et al., 2000). One possible approach to further investigate this possibility would be to conduct co-phylogeographic analyses between the genetic structure of the parasite and that of their cockle and seabass hosts (e.g. Nieberding et al., 2004; Nieberding and Olivieri, 2007).

In conclusion, with a single COI haplotype observed per host, B. minimus genetic variability at the host level was very limited. This suggests that a single miracidium may be infecting C. edule, possibly indicating the existence of strong mechanisms operating in the background to reduce the likelihood of multiple infections. However, the processes involved are currently poorly understood, with a long way to go to fully comprehend trematode host–parasite interactions. Particularly, when discussing the early stage of the parasite's life cycle (sporocysts), laboratory studies are required to confirm the likelihood of co-infections and the defence mechanisms involved. Similar to the population structure of its first host, C. edule, 2 B. minimus groups (north and south of Bay of Biscay) were found. It is probable that this pattern derives from D. labrax, this parasite's final host, when specific genomic areas under selection are investigated. However, this structure may also be the result of sampling limitations mostly in the southern beds. The full parasite population connectivity will only be revealed by additional research, specifically by increasing the number of samples of understudied beds, such as in the Mediterranean.