The genus Contracaecum Raillet and Henry, 1912 includes species maturing in pinnipeds and waterbirds worldwide, representing the most diverse group of the family Anisakidae (1, 2). Considering only the species of the genus maturing in fish-eating birds, 19 species have been so far described as infecting cormorants and pelicans (1, 3–10). Moreover, several morphospecies of the genus have been demonstrated to include sibling species (1, 2, 8, 11–13).

Among them, the Contracaecum rudolphii (s.l.) species complex has its most habitual definitive hosts in cormorants belonging to the genus Phalacrocorax (Family Phalacrocoracidae) (4, 7, 8, 12, 14, 15). The two cryptic species of the complex parasitizing cormorants belonging to the Family Phalacrocoracidae identified in Europe are C. rudolphii sp. A and C. rudolphii sp. B (7, 15–17). They are mainly parasites in the Great Cormorant, Ph. carbo sinensis. Regarding other sibling species so far included in the C. rudolphii (s.l.) complex, C. rudolphii sp. C has been detected in Ph. auritus from Florida (3). In Australia, the two cryptic species, i.e., C. rudolphii sp. D and C. rudolphii sp. E, have been found as adult parasites in Ph. carbo and Ph. varius (8). Moreover, another species indicated as C. rudolphii sp. F has been identified as parasitizing the Great Pelican, Pelecanus occidentalis, from the Gulf of Mexico (4). Other Contracaecum species identified in cormorants are as follows: C. microcephalum detected in Ph. melanoleucos and Ph. pygmaeus from Australia and Montenegro, respectively (3, 9); C. chubutensis described in Ph. atriceps from the north Patagonian coast (5); C. australe discovered in Ph. brasilianus from Chile and Argentina (18); and C. jorgei observed in Nannopterum brasilianus from Argentina (10).

The European Shag, Phalacrocorax aristotelis, is a cormorant species commonly found along European marine coasts, mainly observed on rocky seasides and islands. Its nesting habitat in Spain includes a wide range of inaccessible sites on the coast, such as cavities, fissures, and cliff ledges. It is generally sedentary, being the individuals faithful to their first breeding area (19, 20). It comprises three subspecies: Ph. aristotelis aristotelis, which inhabits the north-east Atlantic coast, from Iceland to the Iberian Peninsula, including Norway and the Kola Peninsula in Russia; Ph. aristotelis riggenbachi, which breeds in the Atlantic coast of Morocco; and finally, Ph. aristotelis desmarestii, which is endemic of the Mediterranean and the Black Sea (21). Its diet includes a wide range of fish species, mainly demersal and/or benthic, and could be influenced by their habitat characteristics or the stage of their breeding cycle (22). Studies regarding nematodes' presence in the European Shag are scarce. In Spain, a previous analysis (23) reported the presence of C. rudolphii sp. A and C. septentrionale in Ph. a. aristotelis hosts from the north-west Atlantic coast (Galician coast). In the Mediterranean Sea, C. rudolphii sp. A has been identified in Ph. aristotelis from the east coast of Sardinia (24), while C. septentrionale has been reported in the same hosts from Iceland (25, 26).

The aim of the present study was to (1) identify genetically nematodes of the genus Contracaecum parasitizing the European Shag, Ph. aristotelis desmarestii from the north-western Mediterranean Sea (coast of Spain); and (2) investigate the population genetic structure of the parasite species in the western Mediterranean Sea areas.

The parasitological analysis of the Ph. aristotelis desmarestii individuals allowed the identification of 1,517 nematode specimens which were morphologically referred to as third-stage larvae (L3), preadults, and adults of C. rudolphii (s.l.). All cormorants were found to be parasitized with nematode specimens (Prevalence, P = 100%), being the parasites found in the stomach and along the esophagus of the hosts.

Considering the obtained sequences, the best-fit substitution models were GTR+I+G (I+G = 0.5482) and TVM+G (G = 0.3353) for mtDNA cox2 and ITS rDNA regions, respectively. The mtDNA cox2 and the ITS1 and ITS2 rDNA sequences analysis of the C. rudolphii (s.l.) specimens allowed us to identify all the nematode specimens as C. rudolphii sp. A, including 33 adults, 21 preadults, and 21 larvae. DNA sequences from mtDNA cox2 and ITS1 and ITS2 rDNA regions were deposited in GenBank under the accession numbers OP690535-40 (mtDNA cox2), OP677791-96 (ITS1 region of rDNA), and OP677785-90 (ITS2 region of rDNA). Indeed, the Bayesian inference (BI) phylogenetic tree, based on mtDNA cox2 sequences (Figure 2), indicated that all the C. rudolphii sp. A sequences here obtained from specimens parasitizing Ph. a. desmarestii, clustered together, with high probability values in a well-supported clade, with other C. rudolphii sp. A sequences previously deposited in GenBank from specimens parasitizing Ph. carbo sinensis (accession numbers EF513502 and EF558891). This clade was clearly distinct from the other two clades obtained for the other sibling species, C. rudolphii sp. B and C. rudolphii sp. F, reported in Ph. carbo and P. occidentalis, respectively. Similarly, sequences of the other Contracaecum species included in the tree, parasitizing fish-eating bird hosts of the genera Phalacrocorax and Pelicanus, were grouped into distinct and well-supported clades (Figure 2). In the same vein, the BI phylogenetic tree based on rDNA ITS1 and ITS2 sequences (Figure 3) indicated that all the sequences obtained from Contracaecum spp., previously sequenced at the mtDNA cox2 locus, were grouped in a well-supported clade with a reference sequence of C. rudolphii sp. A from the same host (Ph. aristotelis) (accession number EU678869). This clade was phylogenetically distinct in comparison with the other six clades obtained for each sibling species, i.e., C. rudolphii sp. B, C. rudolphii sp. C, C. rudolphii sp. D, C. rudolphii sp. E, and C. rudolphii sp. F, as well as for the other Contracaecum species, from hosts of the Phalacrocorax and the Pelicanus genera (Figure 3).

The hierarchical gene diversity analysis by AMOVA indicated that 3.69% of the genetic variation of C. rudolphii sp. A from the western Mediterranean localities here considered was expressed among populations while 96.31% was observed within populations. Because no significant intraspecific genetic variation was found among specimens from the Catalan-Spanish coast analyzed in the present work, they were considered as a single population. Estimates of genetic differentiation among C. rudolphii sp. A population of the three analyzed areas from the west Mediterranean, i.e., the Spanish Catalan, the Tyrrhenian, and the Sardinian coasts, inferred from the Fst value and their P-values, are shown in Table 2. Significant genetic differentiation was found among C. rudolphii sp. A metapopulations from the Sardinian coast compared to those from the Spanish region examined here (Fst ≈ 0.08, p < 0.00001). Similarly, significant Fst values were found between C. rudolphii sp. A populations from the Sardinian and Tyrrhenian coasts (Fst ≈ 0.06, p < 0.00001). In contrast, no significant genetic differentiation was observed between the populations from the Mediterranean Spanish coast (Catalan region) and the Tyrrhenian one (Fst ≈ 0.005, p > 0.05). The haplotype parsimony network (TCS) analysis, inferred from cox2 mtDNA of C. rudolphii sp. A, revealed a star-like tree as shown in Figure 4. In the Spanish coast population of the parasite species, 40 haplotypes were detected, while 80 and 28 haplotypes were identified from the Tyrrhenian and Sardinian populations, respectively (Table 3). The Hap7 was the most frequent haplotype shared by all three metapopulations. Hap4 and Hap15 were also frequently observed, shared by several individuals from all geographical areas. Whereas, the Hap3 was mainly shared by individuals from the Spanish coast. However, in these four major haplotypes, the relative frequency of individuals from the Sardinian waters was very low, comprising mainly Spanish and Tyrrhenian individuals. In addition, specimens from the Sardinian coast presented a high frequency of unique haplotypes. As shown in Table 3, the metapopulation from the Sardinian area showed the highest values of nucleotide diversity (per site) (π = 0.112) and the average number of nucleotide differences (K = 5.807). The haplotype diversity (Hd) was high in all three populations, with values ranging between 0.969 and 0.986.

Several Contracaecum species included in the C. rudolphii (s.l.) complex have been identified as parasitizing cormorants, worldwide. In the present study, we report the first genetic identification of C. rudolphii sp. A in the European Shag, Ph. aristotelis desmarestii, from the western Mediterranean Spanish coasts.

The multi-marker genetic approach using the mtDNA cox2 gene and the ITS rDNA region allowed the specific identification of all the nematode samples studied. The genetic identification of the sibling species of the C. rudolphii (s.l.) complex is important in order to understand their ecology and geographical distribution, but also the phylogenetic relationship among them. Both BI phylogenetic trees, based on the mtDNA cox2 gene locus and the ITS1 and ITS2 of the rDNA region, respectively, have shown that the C. rudolphii sp. A specimens were genetically distinct from the other species of the C. rudolphii (s.l.) complex, including sequences of C. rudolphii sp. B, C. rudolphii sp. C, C. rudolphii sp. D, C. rudolphii sp. E, and C. rudolphii sp. F (3, 4, 8, 13, 44). As shown in Figure 3, and in accordance with previous studies, a closer genetic relationship was observed between C. rudolphii sp. A with the other sibling species rather than with C. rudolphii sp. B. (7, 9). In this sense, D'Amelio et al. (3), when analyzing the phylogenetic tree inferred from the rrnS rDNA region, found the closest relationship between C. rudolphii sp. A and C. rudolphii sp. C, rather than between C. rudolphii sp. A and C. rudolphii sp. B. Although heterozygote genotypes between C. rudolphii sp. A and C rudolphii sp. B have been reported, inferred from the ITS rDNA region (38, 45), their genetic analysis should be assessed by using a multi-nuclear diagnostic loci approach, as shown between sibling species of other anisakid nematodes, i.e., the species of the Anisakis simplex (s.l.) complex (46–49). Discovering novel nuclear diagnostic markers between the two species would also improve the investigation of the possible existence of hybridization/introgression events between them.

Cormorants of the Ph. a. desmarestii subspecies are sedentary animals that habit strictly marine ecosystems and are endemic to the Mediterranean Sea (21). Data regarding the nematode presence in this host are scarce. In the Mediterranean waters, C. rudolphii sp. A has been identified parasitizing European Shag specimens on the Sardinian coasts (24). In the north-western Spanish waters (Atlantic Galician coast), the presence of C. rudolphii sp. A was reported at high prevalence values parasitizing cormorants of the subspecies Ph. aristotelis aristotelis (23). Those records, combined with the present finding of C. rudolphii sp. A in Ph. aristotelis desmarestii from the western Mediterranean Spanish coast and the absence of the sibling species C. rudolphii sp. B, confirm the adaptation of the former anisakid species to marine and brackish-water ecosystems, as suggested by Mattiucci et al. (7, 15). This hypothesis was supported by the high relative proportion of C. rudolphii sp. A found in Ph. carbo sinensis habiting in brackish environments of Italy, whereas the sibling species C. rudolphii sp. B infecting the same cormorant species was mainly found living in freshwater environments, such as rivers and lakes of Italy (7). Similar proportions of infection, with the species C. rudolphii sp. A from brackish cormorants (38) and with C. rudolphii sp. B in cormorants from freshwater environments (50) have been reported later. It has been hypothesized that the observed differential spatial distribution of C. rudolphii sp. A and C. rudolphii sp. B in distinct aquatic ecosystems could also be related to the differential feeding ecology and the behavior of distinct wintering populations of Ph. carbo sinensis in European areas (7). In this sense, the sedentary pattern and the feeding behavior of the Ph. a. desmarestii, the definitive host herein studied, would be responsible for the maintenance of the life cycle of C. rudolphii sp. A along the Catalan marine coastal ecosystem. The low mobility of the Ph. a. desmarestii could also help to maintain the genetic differentiation of the C. rudolphii sp. A population observed in the present study, justifying the relatively high number of private haplotypes (n = 26) shown by the parasite species in this area (Figure 4). On the other hand, the migration pattern of the Great Cormorant, Ph. carbo, which acts as the main definitive host of C. rudolphii sp. A in European waters might explain the population genetic differentiation here found between the metapopulation of C. rudolphii sp. A from Sardinia in comparison with those conspecific ones from the Catalan and Tyrrhenian areas (Table 2). This fish-eating bird presents its high migration pattern along the brackish and freshwater environments of Europe, inhabiting, among the others, the Catalan and the Tyrrhenian Mediterranean waters during its non-breeding season (51, 52). Although the low number of samples from the Sardinian Sea herein studied (n = 33) should be considered, the existence of a population of Ph. carbo living during the winter season in Sardinia, which could have a different breeding area in comparison with those populations of Ph. carbo wintering in the Spanish and Tyrrhenian areas, might explain the parasite genetic differentiation observed in the present study. In support of this hypothesis, a higher presence through Sardinia has been observed in some Ph. carbo populations during its migration routes (53). Moreover, Marion and Gentil (54) detected the presence of two genetically distinct cormorant populations of Ph. carbo, one along the Sardinia coast and the other one along the Italic Peninsula, as inferred from the mtDNA cox2 gene, which could also explain the present finding of a certain genetic sub-structuring of C. rudolphii sp. A, harbored by cormorants analyzed from those areas, i.e., Sardinia (38) and Latium region (7). Indeed, C. rudolphii sp. A from cormorants of the Sardinian coast also showed a high number of private haplotypes (n = 23), not shared with the other populations, despite the lower number of specimens analyzed, which is in accordance with the estimation of Fst values herein obtained. In addition, the most frequent haplotypes observed in Figure 4 seemed to be shared mainly by specimens from the Spanish and the Tyrrhenian coasts. The influence of the population structure and the migration route of the definitive hosts in shaping the genetic structure of anisakid nematodes has also been observed in Anisakis spp. (55) and Uncinaria lucasi (56). In this sense, it is suggested that the interaction among definitive, intermediate, and paratenic hosts influences the genetic structure of parasite populations, being the geographic origin of the host species a major driver (55). However, the analysis of a high number of specimens of C. rudolphii sp. A from cormorants collected from different brackish and marine ecosystems of European countries, as well as knowledge on the migration patterns of different populations of its phalacrocoracid definitive host, is necessary in order to shed light on the population genetic structure of this waterbird anisakid species.

The high prevalence of C. rudolphii sp. A (P = 100%) observed in the present study from the analyzed cormorant specimens is in accordance with that previously reported by Abollo et al. (23) in Ph. a. aristotelis from the Spanish Atlantic coast and with those recorded in Ph. c. sinensis from the Tyrrhenian Sea (7, 38). Comparison with C. rudolphii sp. A from the population of the European Shag of Sardinia is not possible because of the lack of epidemiological data (24). However, despite the high prevalence values herein observed, the life cycle of C. rudolphii sp. A in the Catalan Spanish coasts is not clear. Various studies analyzing several intermediate/paratenic hosts from this area, such as Micromesistius poutassou (57), Sardina pilchardus (58), Phycis blennoides (59), and Merluccius merluccius (60), among others, have not detected Contracaecum spp. larvae in these fish. Nevertheless, along the Spanish Mediterranean coast, some studies have identified Contracaecum spp. larvae parasitizing other fish species, including Anguilla anguilla (61), Scomber scombrus, Trachurus trachurus (62), Pagellus erythrinus, Trisopterus minutus, and Engraulis encrasicolus (59). Still, specimens were identified only at the genus level, using morphological criteria and without performing the molecular identification of the parasites. In other areas of the Mediterranean, such as Sardinian waters and the Tyrrhenian Sea, specific identification of C. rudolphii sp. A parasitizing various intermediate/paratenic hosts, including Sparus aurata (63), Trachurus sp., Mullus sp., Illex coindettii (64), A. anguilla, Dicentrarchus labrax (7, 65), Gobius niger, G. paganellus, and Solea solea (65) has been reported. Studies regarding the genetic identification of Contracaecum larvae in fish species, likely included in the prey items list of the Ph. a. desmarestii, from the western Mediterranean Spanish coast, should be accomplished for a better understanding of the C. rudolphii sp. A life cycle in the studied area.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ supplementary material.

Ethical review and approval was not required for the animal study because the cormorants used in the present study were found dead and no permit and/or ethical consideration were needed.

MC, CT, and JB performed a cormorant necropsy and the collection and morphological identification of parasites. XR-G, MP, and SM contributed to the design of the study, performed molecular and phylogenetic analyses of selected parasites, and wrote the manuscript. All authors contributed to the manuscript revision, and read and approved the submitted version.