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., Reference Carlson, Dallas, Alexander, Phelan and Phillips2020). They can be found at almost across all trophic levels of dynamic food chains (Bartoli and Gibson, Reference Bartoli and Gibson2007). Trematodes are an important component of ecosystem biodiversity with significant impacts at the host individuals (SchulteOehlmann et al., Reference SchulteOehlmann, Oehlmann, Fioroni and Bauer1997; Curtis et al., Reference Curtis, Kinley and Tanner2000; Thieltges, Reference Thieltges2004), host populations (Fredensborg et al., Reference Fredensborg, Mouritsen and Poulin2005) and ecosystem communities (Poulin, Reference Poulin1999; Mouritsen and Poulin, Reference Mouritsen and Poulin2002; Goedknegt et al., Reference Goedknegt, Feis, Wegner, Luttikhuizen, Buschbaum, Camphuysen, van der Meer and Thieltges2016). 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., Reference Johnson, Dobson, Lafferty, Marcogliese, Memmott, Orlofske, Poulin and Thieltges2010; Hatcher and Dunn, Reference Hatcher and Dunn2011).
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., Reference Cribb, Bray, Olson and Littlewood2003; Bartoli and Gibson, Reference Bartoli and Gibson2007). 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., Reference Dubois, Savoye, Sauriau, Billy, Martinez and de Montaudouin2009; Preston et al., Reference Preston, Orlofske, Lambden and Johnson2013), with direct consequences for host reproduction (Carballal et al., Reference Carballal, Iglesias, Santamarina, Ferro-Soto and Villalba2001), growth (Bowers, Reference Bowers1969) and energy demand (Jokela et al., Reference Jokela, Uotila and Taskinen1993), leading to eventual host death (Thieltges, Reference Thieltges2006). 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, Reference Fredensborg and Poulin2005; Mideo, Reference Mideo2009). 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., Reference Dumont, Mone, Mouahid, Idris, Shaban and Boissier2007; Balmer et al., Reference Balmer, Stearns, Schotzau and Brun2009).
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., Reference Rauch, Kalbe and Reusch2005; Keeney et al., Reference Keeney, Waters and Poulin2007; Lagrue et al., Reference Lagrue, McEwan, Poulin and Keeney2007), showing that the likelihood of infection by conspecifics increased with the prevalence of the parasite in the community (Keeney et al., Reference Keeney, Boessenkool, King, Leung and Poulin2008; Louhi et al., Reference Louhi, Karvonen, Rellstab, Louhi and Jokela2013). 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., Reference Agola, Steinauer, Mburu, Mungai, Mwangi, Magoma, Loker and Mkoji2009; Feis et al., Reference Feis, Thieltges, Olsen, de Montaudouin, Jensen, Bazairi, Culloty and Luttikhuizen2015). 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., Reference Theron, Sire, Rognon, Prugnolle and Durand2004; Bell et al., Reference Bell, De Roode, Sim and Read2006; Balmer et al., Reference Balmer, Stearns, Schotzau and Brun2009) or socio-economically important species, namely fish (Vilas et al., Reference Vilas, Paniagua and Sanmartin2003; Criscione and Blouin, Reference Criscione and Blouin2006).
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., Reference Magalhães, Freitas and de Montaudouin2015). 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., Reference Magalhães, Freitas and de Montaudouin2015; de Montaudouin et al., Reference de Montaudouin, Arzul, Cao, Carballal, Chollet, Correia, Cuesta, Culloty, Daffe, Darriba, Díaz, Engelsma, Freitas, Garcia, Goedknegt, Gonzalez, Grade, Groves, Iglesias, Jensen, Joaquim, Lynch, Magalhães, Mahony, Maia, Malham, Matias, Nowaczyk, Ruano, Thieltges and Villalba2021). 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., Reference Desclaux, de Montaudouin and Bachelet2002; de Montaudouin et al., Reference de Montaudouin, Thieltges, Gam, Krakau, Pina, Bazairi, Dabouineau, Russell-Pinto and Jensen2009). Infection by B. minimus results in castration (Carballal et al., Reference Carballal, Iglesias, Santamarina, Ferro-Soto and Villalba2001) and energy consumption (Dubois et al., Reference Dubois, Savoye, Sauriau, Billy, Martinez and de Montaudouin2009), 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., Reference Magalhães, Freitas and de Montaudouin2015; de Montaudouin et al., Reference de Montaudouin, Arzul, Cao, Carballal, Chollet, Correia, Cuesta, Culloty, Daffe, Darriba, Díaz, Engelsma, Freitas, Garcia, Goedknegt, Gonzalez, Grade, Groves, Iglesias, Jensen, Joaquim, Lynch, Magalhães, Mahony, Maia, Malham, Matias, Nowaczyk, Ruano, Thieltges and Villalba2021). 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., Reference Pina, Barandela, Santos, Russell-Pinto and Rodrigues2009; Magalhães et al., Reference Magalhães, Freitas and de Montaudouin2015).
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., Reference Magalhães, Correia, de Montaudouin and Freitas2018]); 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., Reference Magalhães, Freitas and de Montaudouin2015]) (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).
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−1.
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., Reference Moszczynska, Locke, McLaughlin, Marcogliese and Crease2009), 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 MgCl2, 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 BigDyeTM 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, Reference Hall1999) 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, Reference Librado and Rozas2009). 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., Reference Rice, Longden and Bleasby2000; Goujon et al., Reference Goujon, McWilliam, Li, Valentin, Squizzato, Paern and Lopez2010).
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, Reference Paradis2010; Paradis and Schliep, Reference Paradis and Schliep2019).
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., Reference Feis, Thieltges, Olsen, de Montaudouin, Jensen, Bazairi, Culloty and Luttikhuizen2015 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, Reference Librado and Rozas2009). Using MEGA X software (Kumar et al., Reference Kumar, Stecher, Li, Knyaz and Tamura2018), 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, Reference Paradis2010; Schliep, Reference Schliep2011; Yu et al., Reference Yu, Smith, Zhu, Guan and Lam2017; Paradis and Schliep, Reference Paradis and Schliep2019).
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., Reference Excoffier, Laval and Schneider2005). 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., Reference Souche, Hellemans, Babbucci, MacAoidh, Guinard, Bargelloni, Chistiakov, Patarnello, Bonhomme, Martinsohn and Volckaert2015; Vera et al., Reference Vera, Maroso, Wilmes, Hermida, Blanco, Fernández, Groves, Malham, Bouza, Cockle's, Robins and Martínez2022] 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).
Coordinates of each cockle bed, the number of analysed cockles (N cockles), number of parasites sequenced (NB. 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).
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).
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).
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.
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, Reference Poulin2001; Read and Taylor, Reference Read and Taylor2001). This pattern has been extensively studied for several parasite species, namely with an impact on human health (Theron et al., Reference Theron, Sire, Rognon, Prugnolle and Durand2004; Bell et al., Reference Bell, De Roode, Sim and Read2006). For example, in the case of malaria, more than 5 strains have been found to be infecting the same host (Bell et al., Reference Bell, De Roode, Sim and Read2006). The same trend was observed for the trematode parasite Schistosoma mansoni within their second and final host (Theron et al., Reference Theron, Sire, Rognon, Prugnolle and Durand2004). 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., Reference Rauch, Kalbe and Reusch2005; Keeney et al., Reference Keeney, Waters and Poulin2007; Lagrue et al., Reference Lagrue, McEwan, Poulin and Keeney2007). 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, Reference Abdul-Muneer2014) 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., Reference Queiroz, Humphries, Noble, Santos and Sims2012), 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., Reference Correia, Magalhães, Freitas, Bazairi, Gam and de Montaudouin2020). Trematodes are highly sensitive to temperature, both in their free-living and parasitic stages (Thieltges and Rick, Reference Thieltges and Rick2006; Selbach and Poulin, Reference Selbach and Poulin2020). For example, the production and hatching rate of trematode eggs are positively correlated with temperature, peaking under ideal thermal conditions (Morley, Reference Morley2012; Morley and Lewis, Reference Morley and Lewis2017). The same happens with cercarial multiplication within and emergence from the first intermediate host (Poulin, Reference Poulin2006; de Montaudouin et al., Reference de Montaudouin, Blanchet, Desclaux-Marchand, Lavesque and Bachelet2016). 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., Reference Magalhães, Freitas, Dairain and de Montaudouin2017) 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., Reference Magalhães, Freitas and de Montaudouin2015, Reference Magalhães, Daffe, Freitas and de Montaudouin2020), 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., Reference Keeney, Boessenkool, King, Leung and Poulin2008; Louhi et al., Reference Louhi, Karvonen, Rellstab, Louhi and Jokela2013). 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., Reference Baiocchi, Lee, Choe and Dillman2017) or be toxic against competitors (Burman, Reference Burman1982; Selva et al., Reference Selva, Viana, Regev-Yochay, Trzcinski, Corpa, Lasa, Novick and Penades2009). 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, Reference Thieltges2006; de Montaudouin et al., Reference de Montaudouin, Arzul, Cao, Carballal, Chollet, Correia, Cuesta, Culloty, Daffe, Darriba, Díaz, Engelsma, Freitas, Garcia, Goedknegt, Gonzalez, Grade, Groves, Iglesias, Jensen, Joaquim, Lynch, Magalhães, Mahony, Maia, Malham, Matias, Nowaczyk, Ruano, Thieltges and Villalba2021). 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, Reference Read and Taylor2001; Davies et al., Reference Davies, Fairbrother and Webster2002; Alizon and van Baalen, Reference Alizon and van Baalen2008). Regardless, a study conducted using different strains of the trematode S. mansoni in their host showed an increase in overall pathogenicity (Davies et al., Reference Davies, Fairbrother and Webster2002).
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., Reference Rajakumar, Bleiss, Hartmann, Schierack, Marko and Lucius2006). The extensive use of host tissues, particularly the gonad and digestive gland (Dubois et al., Reference Dubois, Savoye, Sauriau, Billy, Martinez and de Montaudouin2009), 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., Reference Magalhães, Daffe, Freitas and de Montaudouin2020). 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., Reference Feis, Thieltges, Olsen, de Montaudouin, Jensen, Bazairi, Culloty and Luttikhuizen2015), 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. (Reference Feis, Thieltges, Olsen, de Montaudouin, Jensen, Bazairi, Culloty and Luttikhuizen2015). 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, Reference Loewe and Hill2010). 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., Reference Feis, Thieltges, Olsen, de Montaudouin, Jensen, Bazairi, Culloty and Luttikhuizen2015). Trematodes have a complex life cycle that involves a vertebrate, usually a fish or a bird, as definitive host (Cribb et al., Reference Cribb, Bray, Olson and Littlewood2003). 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., Reference de Pontual, Heerah, Goossens, Garren, Martin, Le Ru, Le Roy and Woillez2023). 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., Reference Vera, Maroso, Wilmes, Hermida, Blanco, Fernández, Groves, Malham, Bouza, Cockle's, Robins and Martínez2022)]. Cockle dispersal mainly occurs during their larval stages (Martel and Chia, Reference Martel and Chia1991), when B. minimus cannot infect them (Magalhães et al., Reference Magalhães, Freitas and de Montaudouin2015). 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., Reference Zemmer, Detwiler, Sokol, Da Silva Neto, Wyderko, Potts, Gajewski, Sarment, Benfield and Belden2020). Two distinct genetic population units, one in the Atlantic area and another in the Mediterranean, have been previously identified for seabass (Souche et al., Reference Souche, Hellemans, Babbucci, MacAoidh, Guinard, Bargelloni, Chistiakov, Patarnello, Bonhomme, Martinsohn and Volckaert2015), 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., Reference Lemaire, Versini and Bonhomme2005; Vandeputte et al., Reference Vandeputte, Gagnaire and Allal2019). 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. (Reference Souche, Hellemans, Babbucci, MacAoidh, Guinard, Bargelloni, Chistiakov, Patarnello, Bonhomme, Martinsohn and Volckaert2015). 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. (Reference Quéré, Guinand, Kuhl, Reinhardt, Bonhomme and Desmarais2010) 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., Reference Souche, Hellemans, Babbucci, MacAoidh, Guinard, Bargelloni, Chistiakov, Patarnello, Bonhomme, Martinsohn and Volckaert2015; 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., Reference de Pontual, Heerah, Goossens, Garren, Martin, Le Ru, Le Roy and Woillez2023), 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., Reference Sasal, Durand, Faliex and Morand2000). 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., Reference Nieberding, Morand, Libois and Michaux2004; Nieberding and Olivieri, Reference Nieberding and Olivieri2007).
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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182023000987
Data availability statement
The data that support the findings of this manuscript are provided in the text and are also available from the corresponding authors (S.C. and M.V.) upon reasonable request.
Acknowledgements
The authors are grateful to Txetxu Santiago and Leslie Stout for their help in cockle sampling and Cristina Alonso for her technical support. We are also grateful to Professor J. Russell Stothard and the anonymous reviewers for constructive comments on an earlier version. Simão Correia benefited from PhD grant (2020.04688.BD) funded by National Funds through the Portuguese Science Foundation (Fundação para a Ciência e a Tecnologia, FCT). Luísa Magalhães is funded by national funds (OE), through FCT (grant 2021.01858.CEECIND). Thanks are also due to CESAM (UIDP/50017/2020 + UIDB/50017/2020 + LA/P/0094/2020). Sampling in Arcachon (France) was performed thanks to Planula 4 vessel (CNRS-INSU, Flotte Océanographique Française) and its staff.
Author contributions
S.C., S.F.-B., L.M. and M.V. conceived and designed the study. S.C. conducted data gathering. X.d.M. and G.D. contributed with data. S.C. and M.V. performed data analyses. S.C. wrote the original draft. S.F.-B., L.M., X.d.M., R.P. and M.V. reviewed and edited the draft. S.F.-B., L.M., X.d.M., G.D. and M.V. provided resources. S.F.-B., L.M., R.P. and M.V. supervised the work.
Financial support
The COCKLES project (EAPA_458/2016 COCKLES Co-Operation for Restoring CocKle SheLfisheries and its Ecosystem Services in the Atlantic Area) contributed with samples for the phylogeographic study. The project was funded by the Interreg Atlantic Area Programme through the European Regional Development Fund.
Competing interests
None.
Ethical standards
Not applicable.
Open access
