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

Parasitology cambridge.org/par Research Article Cite this article: Correia S, Fernández-Boo S, Magalhães L, de Montaudouin X, Daffe G, Poulin R, Vera M (2023). Trematode genetic patterns at host individual and population scales provide insights about infection mechanisms. Parasitology 150, 1207–1220. https://doi.org/10.1017/S0031182023000987 Received: 22 May 2023 Revised: 2 October 2023 Accepted: 7 October 2023 First published online: 20 October 2023 Keywords: Bucephalus minimus;Cerastoderma edule; clonal diversity; COI; host–parasite interactions; parasite; population genetics Corresponding authors: Simão Correia; Email: [email protected]; [email protected]; Manuel Vera; Email: manuel.ver[email protected] © The Author(s), 2023. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited. Trematode genetic patterns at host individual and population scales provide insights about infection mechanisms Simão Correia1,2,3,4 , Sergio Fernández-Boo2, Luísa Magalhães1, Xavier de Montaudouin5, Guillemine Daffe6, Robert Poulin4 and Manuel Vera3 1 Department of Biology, CESAM, University of Aveiro, 3810-193 Aveiro, Portugal; 2 Aquatic and Animal Health Group, CIIMAR, University of Porto, 4450-208 Matosinhos, Portugal; 3 Department of Zoology, Genetics and Physical Anthropology, Campus Terra, University of Santiago de Compostela, 27002 Lugo, Spain; 4 Department of Zoology, University of Otago, 9054 Dunedin, New Zealand; 5 CNRS, UMR EPOC, Station Marine, University of Bordeaux, F-33120 Arcachon, France and 6 Université de Bordeaux, CNRS, Observatoire Aquitain des Sciences de l’Univers, F-33615 Pessac, France 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 multihost 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 https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press 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). 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., 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 mMMgCl 2 , 100 μMdeoxynucleotide triphosphates (dNTPs), 0.5 μMof 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 1208 Simão Correia et al. https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press 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). 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 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. Parasitology 1209 https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press 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;Yuet 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 Table 1. Accession number for B. minimus COI gene DNA sequences downloaded (with reference) and deposited (in bold) in GenBank Haplotype name GenBank accession number Reference LaA KF880428.1 Feis et al.(2015) LaB KF880429.1 Feis et al.(2015) LaC KF880430.1 Feis et al.(2015) LaD KF880431.1 Feis et al.(2015) LaE KF880432.1 Feis et al.(2015) LaF KF880433.1 Feis et al.(2015) LaG KF880434.1 Feis et al.(2015) LaH KF880435.1 Feis et al.(2015) LaI KF880436.1 Feis et al.(2015) LaJ KF880437.1 Feis et al.(2015) LaK KF880438.1 Feis et al.(2015) LaL KF880439.1 Feis et al.(2015) LaM KF880440.1 Feis et al.(2015) LaN KF880441.1 Feis et al.(2015) LaO KF880442.1 Feis et al.(2015) LaP KF880443.1 Feis et al.(2015) LaQ KF880444.1 Feis et al.(2015) LaR KF880445.1 Feis et al.(2015) LaS KF880446.1 Feis et al.(2015) LaT KF880447.1 Feis et al.(2015) LaU KF880448.1 Feis et al.(2015) LaV KF880449.1 Feis et al.(2015) LaW KF880450.1 Feis et al.(2015) LaX KF880451.1 Feis et al.(2015) LaY KF880452.1 Feis et al.(2015) LaZ KF880453.1 Feis et al.(2015) LaAA KF880454.1 Feis et al.(2015) LaAB KF880455.1 Feis et al.(2015) LaAC KF880456.1 Feis et al.(2015) LaAD KF880457.1 Feis et al.(2015) LaAE KF880458.1 Feis et al.(2015) LaAF KF880459.1 Feis et al.(2015) LaAG KF880460.1 Feis et al.(2015) LaAH KF880461.1 Feis et al.(2015) LaAI KF880462.1 Feis et al.(2015) LaAJ KF880463.1 Feis et al.(2015) LaAK KF880464.1 Feis et al.(2015) LaAL KF880465.1 Feis et al.(2015) LaAM KF880466.1 Feis et al.(2015) LaAN KF880467.1 Feis et al.(2015) LaAO KF880468.1 Feis et al.(2015) LaAP KF880469.1 Feis et al.(2015) LaAQ KF880470.1 Feis et al.(2015) LaAR KF880471.1 Feis et al.(2015) LaAS KF880472.1 Feis et al.(2015) (Continued) Table 1. (Continued.) Haplotype name GenBank accession number Reference LaAT KF880473.1 Feis et al.(2015) LaAV KF880474.1 Feis et al.(2015) LaAW KF880475.1 Feis et al.(2015) LaAX KF880476.1 Feis et al.(2015) LaAY KF880477.1 Feis et al.(2015) LaAZ KF880478.1 Feis et al.(2015) LaBA KF880479.1 Feis et al.(2015) LaBB KF880480.1 Feis et al.(2015) LaBC KF880481.1 Feis et al.(2015) BmA OQ625925 This study BmB OQ625926 This study BmC OQ625927 This study BmD OQ625928 This study BmE OQ625929 This study BmF OQ625930 This study BmG OQ625931 This study BmH OQ625932 This study BmI OQ625933 This study BmJ OQ625934 This study BmK OQ625935 This study BmL OQ625936 This study 1210 Simão Correia et al. https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press (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). 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). 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). 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). 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 (Pvalue < 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 (Pvalue < 0.001), suggesting the singularity of this bed (global ϕ ST in the whole region excluding Merja Zerga = 0.1706, Pvalue < 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, Pvalue = 0.073) than the southern one (composed by Merja Zerga, Aveiro, Baiona and Noia; ϕ ST = 0.6694, Pvalue < 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 (Pvalue = 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, Pvalue = 0.008), this percentage being three times higher than those assigned to differences among beds within groups (ϕ SC = 0.2109, Pvalue < 0.001, percentage of genetic differentiation = 13.80%). The AMOVA model including Arcachon in the northern group yielded similar values (ϕ CT = 0.3100, Pvalue = 0.006, percentage of genetic differentiation = 31.00%; ϕ SC = 0.1567, Pvalue < 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, Pvalue = 0.048, percentage of genetic differentiation = 13.27%; ϕ SC = 0.2299, Pvalue < 0.001, percentage of genetic differentiation = 19.13%), suggesting a more coherent grouping of the beds in the former model. 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) Cockle bed Coordinates N Cockles B. minimus Prevalence kPSHaplotype composition Aveiro 40.710123, −8.704596 5 20 3.3% 5 8 BmA, BmB, BmC, BmD, BmE Baiona 42.117020, −8.820283 5 20 16.7% 3 3 BmF, LaD, LaE (3) Arcachon 44.690111, −1.182944 4 16 23.3% 4 5 BmG, LaAE, LaAQ, LaAW General 14 56 8.1% 12 14 Parasitology 1211 https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press 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 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. 1212 Simão Correia et al. https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press Table 3. Haplotype frequency for B. minimus COI gene per cockle bed including haplotype diversity (h± standard deviation) and nucleotide diversity (π± standard deviation) Haplotype Merja Zerga Aveiro Baiona Noia Arcachon Bay of Somme English Channel Celtic Sea Burry Inlet The Dee Wadden Sea Sum LaA 000000000011 LaB20000200000022 LaC 000010000001 LaD 001010000002 LaE 003020100006 LaF 000000010001 LaG 000000100001 LaH000000000011 LaI 500000000005 LaJ 000000081009 LaK 000000010001 LaL 100000000001 LaM000000010001 LaN000000010001 LaO000010000001 LaP 000000000011 LaQ000000000011 LaR 000000110114 LaS 000000010001 LaT 000000010001 LaU000000000011 LaV 000000020002 LaW000010000001 LaX 000000100001 LaY 000000010001 LaZ 000000010001 LaAA 0 0 0 0 1 0 0 0 0 0 0 1 LaAB 0 0 0 0 1 0 0 0 0 0 0 1 LaAC000010000001 LaAD 0 0 0 0 0 0 0 0 0 0 1 1 LaAE 0 0 0 0 2 0 0 0 0 0 0 2 (Continued) Parasitology 1213 https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press Table 3. (Continued.) Haplotype Merja Zerga Aveiro Baiona Noia Arcachon Bay of Somme English Channel Celtic Sea Burry Inlet The Dee Wadden Sea Sum LaAF 0 0 0 0 0 0 4 1 0 0 0 5 LaAG 0 0 0 0 0 0 0 1 0 0 0 1 LaAH 0 0 0 0 0 0 0 1 0 0 0 1 LaAI 0 0 0 0 0 0 1 0 0 0 0 1 LaAJ 0 0 0 0 0 0 0 0 0 0 1 1 LaAK 0 0 0 0 0 0 0 0 0 0 1 1 LaAL 0 0 0 0 0 0 1 0 0 0 1 2 LaAM 0 0 0 0 0 0 0 1 0 0 0 1 LaAN 0 0 0 0 1 0 2 0 0 0 0 3 LaAO000000000011 LaAP 0 0 0 0 0 0 0 0 0 0 1 1 LaAQ 0 0 0 0 7 2 12 4 1 0 16 42 LaAR 0 0 0 0 0 0 1 0 0 0 0 1 LaAS 0 0 0 0 1 0 0 0 0 0 0 1 LaAT000000100012 LaAV000000000011 LaAW000020000002 LaAX 0 0 0 0 0 0 1 0 0 0 0 1 LaAY000010000001 LaAZ 0 0 0 0 1 0 0 0 0 0 0 1 LaBA 0 0 0 0 0 0 1 0 0 0 0 1 LaBB 0 0 0 0 0 0 2 0 0 0 0 2 LaBC 0 0 0 0 1 0 0 0 0 0 0 1 BmA010000000001 BmB010000000001 BmC010100000002 BmD010000000001 BmE010000000001 BmF001000000001 BmG000010000001 BmH000001000001 BmI000001000001 1214 Simão Correia et al. https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press (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 BmJ000000000101 BmK010000000001 BmL000010000001 Sum 26 6 5 1 29 4 30 27 2 2 30 162 h0.3846 ± 0.1017 1.0000 ± 0.0962 0.7000 ± 0.2184 1.0000 ± 0.0000 0.9384 ± 0.0340 0.8333 ± 0.2224 0.8299 ± 0.0632 0.9003 ± 0.0461 1.0000 ± 0.5000 1.0000 ± 0.5000 0.7241 ± 0.0922 Π0.000688 ± 0.000732 0.007572 ± 0.005013 0.002082 ± 0.001834 0.0000 ± 0.0000 0.006253 ± 0.003626 0.003468 ± 0.002899 0.003077 ± 0.002031 0.005543 ± 0003282 0.001719 ± 0.002426 0.005170 ± 0.005961 0.003172 ± 0.002079 Parasitology 1215 https://doi.org/10.1017/S0031182023000987 Published online by Cambridge University Press