Introduction
Contracaecum sp. nematodes are parasites of aquatic and terrestrial animals with a global distribution (Hartwich, Reference Hartwich1964; Shamsi, Reference Shamsi2019). This genus is considered the most diverse within the family Anisakidae (Hartwich, Reference Hartwich1964), with over one hundred described species (Shamsi , Reference Shamsi2019), and has a significant pathogenic impact on hosts (Shamsi, Reference Shamsi2019). Molecular studies are providing a continuous update of Contracaecum nematode species composition. Mattiucci et al. (Reference Mattiucci, Paoletti, Solorzano and Nascetti2010) reported two new species C. gibsoni and C. overstreeti from Greece, Garbin et al. (Reference Garbin, Mattiucci, Paoletti, González-Acuña and Nascetti2011) described a new taxa C. australe from Chile and D'Amelio et al. (Reference D'Amelio, Cavallero, Dronen, Barros and Paggi2012) reported two new taxa C. fagerholmi and C. rudolphii F from Mexico.
Contracaecum rudolphii Hartwich, Reference Hartwich1964 is worldwide distributed and can be found as an adult in diverse seabirds, as cormorants (Phalacrocorax spp.) and pelicans (Pelecanus spp.) (Hartwich, Reference Hartwich1964; Nadler et al., Reference Nadler, D'Amelio, Fagerholm, Berland and Paggi2000, Reference Nadler, D'Amelio, Dailey, Paggi, Siu and Sakanari2005; Szostakowska et al., Reference Szostakowska, Myjak and Kur2002; Barson and Marshall, Reference Barson and Marshall2004; Li et al., Reference Li, D'Amelio, Paggi, He, Gasser, Lun, Abollo, Turchetto and Zhu2005; Mattiucci et al., Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008, Reference Mattiucci, Sbaraglia, Palomba, Filippi, Paoletti, Cipriani and Nascetti2020; D'Amelio et al., Reference D'Amelio, Cavallero, Dronen, Barros and Paggi2012; Cole and Viney, Reference Cole and Viney2018; Shamsi , Reference Shamsi2019). Epidemiological and molecular studies can help to explain the population genetics and the dispersal pattern of the different Contracaecum and host species, and provide useful data for the conservation of host seabird populations (Cole and Viney, Reference Cole and Viney2018). Hosts can play an important role in increasing the population genetic variability of parasites, by mixing parasites from multiple localities. Indeed, host migration may lead to important gene flow among parasite populations, and so to a poor genetic structuring (Cole and Viney, Reference Cole and Viney2018).
The variety of hosts and the wide range for spicule length of C. rudolphii suggested, since its first description (Hartwich, Reference Hartwich1964), that this taxon included several sibling species. The application of allozyme analysis of C. rudolphii collected from Phalacrocorax aristotelis and P. carbo sinensis in Europe revealed the co-existence of two sibling species, C. rudolphii A and C. rudolphii B (Bullini et al., Reference Bullini, Nascetti, Paggi, Orecchia, Mattiucci and Berland1986; D'Amelio et al., Reference D'Amelio, Nascetti, Mattiucci, Cianchi, Orecchia, Paggi, Berland and Bullini1990; Mattiucci et al., Reference Mattiucci, Cianchi, Nascetti, Paggi, Sardella, Timmi, Webb, Bastida, Rodriquez and Bullini2003, Reference Mattiucci, Sbaraglia, Palomba, Filippi, Paoletti, Cipriani and Nascetti2020; Li et al., Reference Li, D'Amelio, Paggi, He, Gasser, Lun, Abollo, Turchetto and Zhu2005). Subsequently, based on molecular markers, several new species were characterized from different regions and hosts: D'Amelio et al. (Reference D'Amelio, Barros, Ingrosso, Fauquier, Russo and Paggi2007) described C. rudolphii C from P. auritus in west-central Florida; Shamsi et al. (Reference Shamsi, Norman, Gasser and Beveridge2009) described C. rudolphii D from P. carbo and P. varius and C. rudolphii E from P. varius in Australia; D'Amelio et al. (Reference D'Amelio, Cavallero, Dronen, Barros and Paggi2012) described C. rudolphii F from Pelecanus occidentalis in Mexico.
In Europe, several studies confirmed P. carbo as the main host of C. rudolphii (Szostakowska et al., Reference Szostakowska, Myjak and Kur2002; Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2012). Moreover, investigations on the parasites of brackish and freshwater fish showed that C. rudolphii B exists in both ecosystem, while C. rudolphii A occurs mainly in marine and brackish environments (Szostakowska et al., Reference Szostakowska, Myjak and Kur2002; Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2007, Reference Szostakowska and Fagerholm2012; Culurgioni et al., Reference Culurgioni, Sabatini, De Murtas, Mattiucci and Figus2014; Mattiucci et al., Reference Mattiucci, Sbaraglia, Palomba, Filippi, Paoletti, Cipriani and Nascetti2020). However, Szostakowska and Fagerholm (Reference Szostakowska and Fagerholm2012) reported mixed infections in P. carbo, and a rare ‘hybrid’ of C. rudolphii A and B in the Baltic Sea.
Two sub-species of P. carbo are recognized in Europe, differing mainly in their geographical distribution. The sub-species P. carbo carbo (L.) lives in Europe, North Africa and on the Atlantic coast of North America, whereas P. carbo sinensis ranges from western Europe to China, India and North Africa. In Italy, several studies have shown that the breeding colonies of the great cormorant have disappeared in the past centuries, but not in Sardinia (Carpegna et al., Reference Carpegna, Grieco, Grussu, Veronesi and Volponi1997), where between 2006 and 2012, an increase of about 82% in the reproductive number of cormorants (P. carbo sinensis) has been reported, corresponding to a total of 48 colonies.
The aim of the present study is to contribute to extending the knowledge on the Contracaecum taxa that infect the great cormorant in the Mediterranean Sea (Sardinia), using both morphological and molecular data (four mitochondrial and one nuclear genes). This information will allow to provide a phylogenetic analysis of the studied taxa and their congeners, to elucidate the pattern of dispersion of Contracaecum spp. populations, useful to evaluate the possible impact on the host.
Materials and methods
Sampling
The stomach contents of 60 specimens of P. carbo sinensis from six localities of Sardinia were examined for Contracaecum sp. The hosts came from the controlled abatement campaign carried out from 5 to 27 February 2009 in the Oristano coastal brackish water ponds, in order to protect fish production (Decree Regione Autonoma della Sardegna, No 2225/DecA/3, 30/01/2009). The first sorting of the stomach contents was carried out by Merops s.r.l. and the Oristano Department of the Istituto Zooprofilattico Sperimentale della Sardegna. The localities (n = 6) and the number of hosts examined are shown in Fig. 1. Collected specimens of Contracaecum sp. were sorted by locality, host, stage and sex, then fixed in ethanol 70% prior to morphological and molecular analyses.
Fig. 1. Sampling localities of Phalacrocorax carbo sinensis in the Oristano coastal brackish water ponds. In brackets number of hosts examined.
Morphological identification
Nematodes were identified based on the main characters with taxonomic significance, i.e. shape of interlabia (both sexes), number and pattern of pre- and post-cloacal papillae, length and shape of the extremity of spicula (males) (Hartwich, Reference Hartwich, Anderson, Chabaud and Willmott1974; Barus et al., Reference Barus, Sergeeva, Sonin and Ryzhikov1978; Abollo et al., Reference Abollo, Gestal and Pascual2001). Specimens were deposited in the collection of the Department of Zoology, King Saud University.
DNA extraction, amplification and sequencing
From each locality, five hosts were randomly selected, and from each host, six specimens of C. rudolphii (two fourth-stage larvae, two adult males and two females) were stored in 70% ethanol for molecular studies (N total = 180). Genomic DNA was extracted using SDS/proteinase K protocol (Gasser et al., Reference Gasser, Chilton, Hoste and Beveridge1993).
Polymerase chain reactions (PCR) were performed in a 50 μL reaction mix containing 100 ng DNA, ddH2O, 1U Taq polymerase (Bioline, USA), 5X buffer and 5 pmol from each primer (further information on the primer sets are listed in the Supplementary Table S1). One nuclear locus, the ribosomal internal transcribed spacer (ITS) including ITS1 and 5.8S ribosomal RNA gene and ITS2, and four mitochondrial loci, the cytochrome c oxidase 1 (cox1), cytochrome c oxidase 2 (cox2), the NADH dehydrogenase subunit I (nad1) and the small subunit of rRNA (rrnS), were PCR-amplified. Amplifications comprised a denaturation at 95°C (4 min), 35 cycles at 95°C for 60 s, 55/°C (ITS), 50°C (cox1, cox2, nad1 and rrnS) for 90 s and 72°C for 60 s with a final extension at 72°C for 8 min. Negative controls were added to check for eventual contamination. PCR products were sequenced at Macrogen sequencing facility (Macrogen Inc, Seoul, Korea).
Data analysis
The obtained sequences were checked and aligned using Unipro UGENE 1.3 (Okonechnikov et al., Reference Okonechnikov, Golosova and Fursov2012). Representative sequences of Contracaecum sp. (Supplementary Table S2) were included in the datasets from Genbank using BLAST algorithm (Altschul et al., Reference Altschul, Gish, Miller, Myers and Lipman1990). Mixed Sequencer Reader (MSR program) was used to examine mixed chromatograms (Chang et al., Reference Chang, Tsai, Tang, Chen, Lian, Hu, Tsai, Chao, Lai, Wang and Lee2012), useful to identify and compare heterozygous ITS sequences. For each base position, MSR calculates the log ratio of the chromatographic intensities (LRi). The LRi value varies from 0, for heterogeneous chromatographic traces with equivalent intensities, to infinite, for position containing only one main band.
DnaSP v 5.10 (Librado and Rozas, Reference Librado and Rozas2009) was used to estimate the number of haplotypes (h), nucleotide (π) and haplotype (Hd) diversities. Network 5 (Fluxus Technology Ltd., United Kingdom) was used in order to generate Median-Joining network (Bandelt et al., Reference Bandelt, Forster and Rohl1999) depicting evolutionary relationships between haplotypes obtained from each single marker and from the concatenated mitochondrial set. Analysis of molecular variance (AMOVA) was estimated between populations with sufficient sample size, using Arlequin 3.5 (Excoffier and Lischer, Reference Excoffier and Lischer2010). Contracaecum sp. populations were clustered according to locality into six groups (Fig. 1).
Fu's Fs (Fu, Reference Fu1997), Tajima's D (Tajima, Reference Tajima1989) and Mismatch distributions were performed for all mitochondrial markers to evaluate changes in population size using DnaSP v 5.10 (Librado and Rozas, Reference Librado and Rozas2009).
Models of evolution and the best partitioning scheme were inferred, for each marker and for the concatenated sequences, using PartitionFinder 2.1 (Lanfear et al., Reference Lanfear, Frandsen, Wright, Senfeld and Calcott2016). The appropriate models were selected under the best value for the Bayesian Information Criterion (BIC) (Lanfear et al., Reference Lanfear, Frandsen, Wright, Senfeld and Calcott2016). Sequences of Toxocara cati (KY003086, AM411622) and Hysterothylacium auctum (AF115571) were used as outgroup (Supplementary Table S2). Bayesian (BI) and maximum likelihood (ML) phylogenetic trees were constructed using MrBayes v3.0 (Huelsenbeck and Ronquist, Reference Huelsenbeck and Ronquist2001) and RAXML (Stamatakis, Reference Stamatakis2006), respectively. BI analyses were performed with two independent runs of 1 × 108 generations. Trace plots were used to test the number of burn-in fraction (20%). The chains were sampled every 1000 generations. For ML analyses, bootstrap was assessed by 2000 pseudoreplicates.
Results
Morphological identification
The morphological analysis of 15375 nematodes, based on the size and shape of interlabia (both sexes), and the number and pattern of the pre- and post-cloacal papillae, length and shape of spicula (males), allowed to identify adult males (n = 5845) and females (n = 8312) as C. rudolphii sensu lato, and the fourth stage larvae as Contracaecum sp. (n = 1218). The main morphological characters were the following. Adults, cuticle brownish-yellowish and striated transversely, lips and interlabia well developed, interlabia with bilobed distal end. Males, smaller than females, 20 mm (14–25 mm) long, tail conical, curved at the tip, with 27–40 pairs of pre-cloacal papillae and seven pairs of post-cloacal papillae, spicules sub-equal with longitudinal alae and rounded tips. Females, 30 mm (18–48 mm) long, vulva located on the second quarter of body length, tail conical, rounded at the tip. Fourth stage larvae, 7.5 mm (5–12 mm) long, cuticle whitish, with striations at the anterior extremity, oral opening with three lips, interlabia absent, tail conical without papillae.
Table 1 shows the levels of infection of fourth stage larvae and adult males and females according to sex, weight and locality of the host. All the examined specimens of P. carbo sinensis were infected with adults, and 88% with fourth-stage larvae. Total intensity of infection ranged from 15 (Santa Giusta) to 1075 (Cabras). Hosts from Cabras and Mistras showed the highest total mean intensity of infection (334.9 and 332.2, respectively), those from Is Benas and Santa Giusta the lowest (143.6 and 179.7, respectively). In all localities, the intensity of infection of fourth stage larvae was the lowest, and that of females the highest, showing an unbalanced sex ratio (0.7:1.0).
Table 1. Levels of infection of Contracaecum rudolphii in Phalacrocorax carbo sinensis from the Oristano ponds according to sex, weight and locality.
P%, prevalence; Mean I, mean intensity of infection.
Molecular characterization
Sequences of the ribosomal ITS [ITS1 (434 bp) and ITS2 (286 bp)] revealed the presence of C. rudolphii A (n = 157, accession numbers: MT341236, MT341180), C. rudolphii B (n = 22, accession numbers: MT341234, MT341178) (Table 2, Supplementary Fig. S1A). Pairwise comparison revealed seven variable positions within ITS1 and 4 substitutions and an indel (9 bp) within ITS2. One sample from Is Benas showed a mixed chromatogram (accession numbers: MT341235, MT341179). The analysis of the mixed chromatograms with MSR (Supplementary Fig. S1A–C) exhibited a shift of signal intensity at 635 bp, corresponding to an ‘indel’-type heterozygosity. Eleven LRi line drops were observed, corresponding to the heterogeneous fluorescence traces within the chromatogram (Supplementary Fig. S1B, C).
Table 2. Number of identified individuals of Contracaecum rudolphii A, and C. rudolphii B, and heterozygote of these (C. rudolphii A/B), in cormorants from the different localities of the Oristano ponds
All the fourth-stage larvae were identified as C. rudolphii A. Fifteen hosts from all the six localities showed a mixed infection. Of the five P. carbo from Mistras, only one was infected with both Contracaecum species, whereas hosts from Is Benas showed the highest rate of coinfection (4/5).
The length of the PCR products of the mitochondrial markers cox1, cox2, nad1 and rrnS was, approximately, 450, 560, 400 and 530 bp, respectively. For all mtDNA region, no size variation was detected on agarose gel, but among nad1 amplicons, two different patterns were observed, 200 and 400 bp (Fig. 2). No pseudo-genes were identified after the examination of mitochondrial genes translation. The nucleotide variation of the cox1, cox2, nad1 and rrnS was mainly related to changes at the third codon position, while few changes were detected at the first or second position. The G + C contents were 36.8% (cox1), 36.2% (cox2), 30.0% (nad1) and 30.5% (rrnS), respectively (Table 3).
Fig. 2. Gel showing the size variation in PCR products of the mitochondrial nad1, the two different patterns, 200 bp and 400 bp. lanes 1, 3–6, 10 and 11 representative samples of C. rudolphii A and lanes 2, 7–9 of C. rudolphii B.
Table 3. Genetic variability and G + C content of four mitochondrial markers, cox1 (450 bp), cox2 (560 bp), nad1 (400 bp) and rrnS(530 bp), within C. rudolphii A (n = 158) and C. rudolphii B (n = 22), from the Oristano ponds
H, haplotype number; Hd, haplotype diversities; π, nucleotide diversity.
Phylogenetic relationships
For all mitochondrial markers, both C. rudolphii A and C. rudolphii B were characterized by a very high haplotype diversity whereas their respective nucleotide diversity was very low (Table 3). Most haplotypes were unique and represented by a single individual (Table 4).
Table 4. Observed mitochondrial haplotypes of two Contracaecum spp. used in this study, with their numbers, codes and locality
Alignment of the partial cox1 sequences (440 bp) of C. rudolphii A and C. rudolphii B generated 64 and 14 haplotypes, respectively (Table 4, Fig. 3a) (accession numbers: MT338400-77). Contracaecum rudolphii A and C. rudolphii B differed by 33 mutations. For C. rudolphii A, only five haplotypes (HA1, HA2, HA4, HA21 and HA36) were shared among localities. Contracaecum rudolphii A network consisted of two star-like haplogroups, the first included 53 haplotypes with three major haplotypes (two from Is Benas, one from Cabras) and the second one was more aggregated, with one haplotype from S'Ena Arrubia grouping 11 haplotypes. Both network patterns were not dependent on the samples' origin. The heterozygote sample showed the C. rudolphii B pattern with the following haplotype combination: H6cox1-H3cox2-H1nad1-H7rrnS.
Fig. 3. Haplotype network obtained from (a) cox1 (450 bp), (b) cox2 (560 bp), (c) nad1 (400 bp) and (d) rrnS (530 bp). (Ha) C. rudolphii A haplotype, (Hb) C. rudolphii B haplotype. The area of each circle is proportional to the haplotype frequency and mutation numbers are indicated in each branch. The heterozygote sample showed the C. rudolphii B pattern with the following haplotype combination: H6cox1-H3cox2-H1nad1-H7rrnS.
Alignments of the mtDNA cox-2 sequences (540 bp) of C. rudolphii A and C. rudolphii B included 61 and 40 variable sites, respectively (accession numbers: MT338999-038). Network analysis revealed a star-like haplogroup for each studied species. Contracaecum rudolphii A haplogroup consisted of 33 haplotypes aggregated around three major haplotypes (two from Is Benas, one from Cabras) and differing from C. rudolphii B haplogroup by 33 mutations (Table 4, Fig. 3b).
Mitochondrial nad1 sequencing generated variable sequence sizes, 400 bp for all C. rudolphii A haplotype (n = 56) and for one C. rudolphii B haplotype, and 200 bp for the four remaining C. rudolphii B haplotypes (accession numbers: MT338938-99). Multiple alignment of nad1 dataset revealed a homologous region to the forward primer starting from the position 166. The network analysis yielded one star-like haplogroup for each studied species differing by 32 mutations. Contracaecum rudolphii A haplotypes were aggregated around three haplotypes, whereas C. rudolphii B haplotypes were derived from a major haplotype (Table 4, Fig. 3c).
Analysis of rrnS sequences (530 bp) yielded two star-like haplogroups differing by 12 mutations, including C. rudolphii A (n = 39) and C. rudolphii B (n = 14) haplotypes (accession numbers: MT341181-233). Contracaecum rudolphii A samples were divided into two main groups differing in one site from major haplotypes HA1 and HA22 (Table 4, Fig. 3d). The second haplogroup included C. rudolphii B samples, differing from the major haplotype HB1 by one to three mutations.
The obtained matrix of the four concatenated mitochondrial genes consists of 1940 bp after alignment (450 bp for cox1, 560 bp for cox2, 400 bp nad1 and 530 bp for rrnS). In the total matrix, 292 positions were variable of which 231 were parsimony informative. Multiple alignment of the concatenated set generated 75 and 17 mitochondrial haplotypes for C. rudolphii A and C. rudolphii B, respectively (Fig. 4a, b). Network analysis generated two star-like haplogroups, supporting the results of single markers analysis. One major haplotype was shared by three populations, Cabras, Mistras and S'Ena Arrubia, in the case of C. rudolphii A and one was observed in two populations, Cabras and Marceddi, for C. rudolphii B (Fig. 4a, b).
Fig. 4. Haplotype network obtained from concatenated mitochondrial loci [cox1 (450 bp), cox2 (560 bp), nad1 (400 bp) and rrnS(530 bp)]. (a) C. rudolphii A haplotype network, (b) C. rudolphii B haplotype network. The area of each circle is proportional to the haplotype frequency and mutation numbers are indicated in each branch.
AMOVA analysis revealed that the highest amounts of total genetic variation originated from differentiation within populations. Within populations proportions varied from 95% (P < 0.0001) for rrnS to 100% (P < 0.0001) for cox1, cox2 and nad1 (Table 5).
Table 5. Analysis of molecular variance (AMOVA)
Percentage of variation explained by different hierarchical levels for cox1, cox2, nad1 and rrnS loci. Degrees of freedom (d.f.).
Neutrality tests Tajima's D and Fu's Fs revealed for C. rudolphii A significant negative values (Supplementary Table S3) indicating an excess of rare polymorphisms in the studied population suggesting a recent expansion and/or positive selection (Tajima, Reference Tajima1989). Overall Tajima's D and Fu's Fs statistics in C. rudolphii B datasets were not statistically significant and consistent with a population at drift-mutation equilibrium (Pichler, Reference Pichler2002). Mismatch distributions for all markers in C. rudolphii A were unimodal (Supplementary Fig. S2), patterns in concordance with the above neutrality test values supporting the assumption of a sudden expansion model. Mismatch distribution plot C. rudolphii B showed a multimodal shape, revealing demographic equilibrium or a stable population.
The best evolution models and the optimal partition schemes selected for each single marker and for the concatenated mitochondrial set are shown in Table 6. BI analysis and ML yielded trees with similar topology for the nuclear ITS, the four mitochondrial loci (cox1, cox2, nad1 and rrnS) and the concatenated dataset (Fig. 5). All the obtained trees showed two highly supported clades including C. rudolphii A and C. rudolphii B samples, which only differed from the resolution within both clades. Contracaecum ogmorhini (parasite of the otariid Arctocephalus pusillus) appeared within the highly-supported clade of C. rudolphii (avian hosts), as sister group to C. rudolphii B in all mitochondrial trees, or sister group to [C. rudolphii D + (C. rudolphii E + C. rudolphii F)] in ITS tree. Contracaecum eudyptulae and C. septentrionale are sister group to the clade containing C. rudolphii complex (A, B, C, D, E and F) + C. ogmorhini. Only in cox1 tree, C. septentrionale was more divergent with a basal position. Cox2 topology revealed a highly supported clade including C. ogmorhini, to the avian C. margolisi (host: California sea lion). In ITS and cox2 trees, C. osculatum complex, Phocascaris cystophorae C. radiatum and C. miroungae (all from phocid hosts) appeared as sister group to C. microcephalum, C. micropapillatum and C. multipapillatum species that mature in birds. Trees of both BI and ML analyses are shown in Fig. 5, with node support including BI posterior probability and ML bootstrap values.
Fig. 5. Inferred phylogenetic relationships among Contracaecum spp. based on sequences of (a) ITS, (b) cox1, (c) cox2, (d) nad1, (e) rrnS and (f) concatenated mitochondrial loci. The numbers along branches indicate bootstrap probability (BP) values resulting from ML analysis and posterior probability (PP) from the MrBayes analysis.
Table 6. Evolutionary models by gene and codon position found in PartitionFinder.
Discussion
Combination of morphological and molecular data allowed the identification of Contracaecum specimens from P. carbo sinensis of the Oristano ponds as C. rudolphii A and C. rudolphii B. It is confirmed that sequence analysis of both nuclear and mitochondrial markers effectively discriminates both Contracaecum spp. The present study confirms that the two Contracaecum spp. can coexist in a single host (bird), as previously suggested by Szostakowska and Fagerholm (Reference Szostakowska and Fagerholm2007, Reference Szostakowska and Fagerholm2012). Contracaecum rudolphii A ranges mainly in marine and brackish waters of Europe, whereas C. rudolphii B in the freshwater habitats of the central Europe (Abollo et al., Reference Abollo, Gestal and Pascual2001; Szostakowska et al., Reference Szostakowska, Myjak and Kur2002; Li et al., Reference Li, D'Amelio, Paggi, He, Gasser, Lun, Abollo, Turchetto and Zhu2005; Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2007, Reference Szostakowska and Fagerholm2012; Farjallah et al., Reference Farjallah, Merella, Ingrosso, Rotta, Slimane, Garippa, Said and Busi2008a; Mattiucci et al., Reference Mattiucci, Sbaraglia, Palomba, Filippi, Paoletti, Cipriani and Nascetti2020). The identification of all the fourth stage larvae as C. rudolphii A and none as C. rudolphii B suggests that only specimens of C. rudolphii A had recently infected the cormorants in the same sampling locality, while the infections of C. rudolphii B probably occurred sometime before, and not in the sites of sampling. Indeed, P. carbo sinensis is a migratory species, widespread in the inland and coastal wetlands of central Europe and the Mediterranean Sea. Thus, it is probable that C. rudolphii B could have infected the host in the central European region, where this species is common in brackish and freshwater habitats (Szostakowska et al., Reference Szostakowska, Myjak and Kur2002; Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2007, Reference Szostakowska and Fagerholm2012; Mattiucci et al., Reference Mattiucci, Sbaraglia, Palomba, Filippi, Paoletti, Cipriani and Nascetti2020). In the same Island of Sardinia (Archipelago di La Maddalena), Farjallah et al. (Reference Farjallah, Merella, Ingrosso, Rotta, Slimane, Garippa, Said and Busi2008a) reported only C. rudolphii A from the European shag P. aristotelis. Unlike its congener, P. aristotelis is a typical marine species, feeding mainly on marine fish, and with a more restricted geographical distribution (western Palearctic).
Several studies confirmed the important role of Phalacrocorax spp., and several other bird species, for the dispersion of Contracaecum spp., and in modelling their genetic variability within the distribution areas (Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2007, Reference Szostakowska and Fagerholm2012; Farjallah et al., Reference Farjallah, Merella, Ingrosso, Rotta, Slimane, Garippa, Said and Busi2008a; Mattiucci et al., Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008, Reference Mattiucci, Sbaraglia, Palomba, Filippi, Paoletti, Cipriani and Nascetti2020). The pattern of distribution of C. rudolphii A and C. rudolphii B has long been debated, and although several studies have demonstrated that both species may occur in fresh, brackish and marine water systems (Nadler, Reference Nadler2000; Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2007, Reference Szostakowska and Fagerholm2012; Farjallah et al., Reference Farjallah, Merella, Ingrosso, Rotta, Slimane, Garippa, Said and Busi2008a; Mattiucci et al., Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008), present results show that C. rudolphii A appears to be more seaward and C. rudolphii B more freshwaterward, as recently suggested by Mattiucci et al. (Reference Mattiucci, Sbaraglia, Palomba, Filippi, Paoletti, Cipriani and Nascetti2020). In fact, the outnumbering of C. rudolphii A compared to C. rudolphii B (157 specimens vs 22), as well as the identification of all fourth stage larvae such as C. rudolphii A suggest that C. rudolphii A was acquired by the host in the studied brackish water ponds; on the other hand, the absence of larval stages of C. rudolphii B suggests that this species has infected the host in environments other than those of the studied lagoons. Considering the migratory habits of P. carbo sinensis, which in Sardinia is mainly present as wintering, it is likely that C. rudolphii B could have infected the host in the breeding grounds of central Europe, where the third-stage larvae of this species have been identified in fish (Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2007; Mattiucci et al., Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008; Molnár et al., Reference Molnár, Székely, Baska, Müller, Zuo, Kania, Nowak and Buchmann2019).
Among all C. rudolphii examined, multiple alignment of ITS sequences combined to chromatogram analysis distinguished one heterozygote C. rudolphii A/B. Previously, only Szostakowska and Fagerholm (Reference Szostakowska and Fagerholm2012) reported a similar heterozygote in the Baltic Sea. Few other studies reported ‘hybrids’ within anisakid species, i.e. A. simplex/A. pegreffii in the Atlantic, Mediterranean Sea and in Western Pacific Ocean (Abollo et al., Reference Abollo, Paggi, Pascual and D'Amelio2003; Marques et al., Reference Marques, Cabral, Busi and D'Amelio2006; Farjallah et al., Reference Farjallah, Busi, Mahjoub, Slimane, Paggi, Said and D'Amelio2008b; Cavallero et al., Reference Cavallero, Ligas, Bruschi and D'Amelio2012, Reference Cavallero, Costa, Caracappa, Gambetta and D'Amelio2014); however, such observations are still very rare. Even when parental species are reproductively isolated, the coexistence of different anisakid congeneric species, with high intensity in single hosts, increases the probability of observing such a rare heterozygote. It should be emphasized that the heterozygote was found in Is Benas, where the hosts showed the highest coinfection rate (4/5). The present results confirm that in the brackish-marine waters of the southwestern Atlantic Ocean and of the southern Mediterranean Europe C. rudolphii A is the most abundant species in Phalacrocorax spp., whereas in the same regions C. rudolphii B is dominant in freshwater habitats (Abollo et al., Reference Abollo, Gestal and Pascual2001; Li et al., Reference Li, D'Amelio, Paggi, He, Gasser, Lun, Abollo, Turchetto and Zhu2005; Farjallah et al., Reference Farjallah, Merella, Ingrosso, Rotta, Slimane, Garippa, Said and Busi2008a, Reference Farjallah, Busi, Mahjoub, Slimane, Paggi, Said and D'Amelio2008b; Mattiucci et al., Reference Mattiucci, Sbaraglia, Palomba, Filippi, Paoletti, Cipriani and Nascetti2020). On the other hand, in the cormorants of northeastern Europe, C. rudoplhii B is always more abundant than C. rudolphii A, although the second species is significantly more abundant in the brackish water habitats compared to freshwater ones, where it is practically absent (Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2012). The coinfection rate of C. rudolphii A and C. rudolphii B ranged from 20 to 80% (mean 50%), a similar value (42%) was reported by Szostakowska and Fagerholm (Reference Szostakowska and Fagerholm2012) in brackish water habitats of northern Europe, whereas the same authors found a much lower coinfection in the freshwater habitats of the same region (17%).
For mitochondrial markers, the AT content was higher than 65%, which agrees with the previous studies of Contracaecum species (Nadler et al., Reference Nadler, D'Amelio, Fagerholm, Berland and Paggi2000; Dzido et al., Reference Dzido, Kijewska and Rokicki2012; Lin et al., Reference Lin, Liu, Zhang, D'Amelio, Zhou, Yuan, Zou, Song and Zhu2012; Delgado and Garcia, Reference Delgado and García2015). Polymorphism analyses of C. rudolphii A and C. rudolphii B based on several markers revealed high Hd and low Pi, which are evidence of an important genetic diversity (Ferreri et al., Reference Ferreri, Qu and Han2011). Similar values have also been reported for other parasite species (Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2007, Reference Szostakowska and Fagerholm2012; Dzido et al., Reference Dzido, Kijewska and Rokicki2012; Lin et al., Reference Lin, Liu, Zhang, D'Amelio, Zhou, Yuan, Zou, Song and Zhu2012; Szudarek et al., Reference Szudarek, Kanarek and Dabert2017). Such pattern of diversity, coupled to the ‘star-like’ haplotype networks with a high ratio of singletons, indicates that the studied C. rudolphii A populations had undergone a bottleneck or founder effect event subsequent to a rapid growth and expansion of the population (Ferreri et al., Reference Ferreri, Qu and Han2011; De Jong et al., Reference De Jong, Wahlberg, van Eijk, Brakefield and Zwaan2011; Levin and Parker, Reference Levin and Parker2013; Dabert et al., Reference Dabert, Coulson, Gwiazdowicz, Moe, Hanssen, Biersma, Pilskog and Dabert2015). Neutrality test values and Mismatch distributions observed in C. rudolphii A samples provide strong evidence for past population expansion. These parameters pointed to different demographic history in the case of C. rudolphii B. However, it is important to point out that in the case of C. rudolphii B, small numbers of samples originated from different localities and hosts were studied. This may conduct to misevaluate the polymorphism level for these samples (Rand et al., Reference Rand, Dorfsman and Kann1994). Thus, it will decrease the ability to detect evolutionary forces acting at the population level (Schmid et al., Reference Schmid, Nigro, Aquadro and Tautz1999).
Population genetic analysis of Contracaecum samples revealed high genetic diversity, and no genetic structure among different localities. The lack of genetic structure is common within parasitic nematodes and several hypotheses have been suggested to clarify such pattern. One hypothesis is the high intensity of infection of these nematodes in the definitive host (mammals and birds) leading to a high effective population size (Ne), which slows down the genetic drift (Cole and Viney, Reference Cole and Viney2018). A second hypothesis is the parasite dispersal model, actually the long distance gene flow mediated by host dispersal dynamic was suggested as the cause of lack of differentiation between parasite populations (Slatkin and Hudson, Reference Slatkin and Hudson1991; Li et al., Reference Li, D'Amelio, Paggi, He, Gasser, Lun, Abollo, Turchetto and Zhu2005; Levin and Parker, Reference Levin and Parker2013; Dabert et al., Reference Dabert, Coulson, Gwiazdowicz, Moe, Hanssen, Biersma, Pilskog and Dabert2015; Cole and Viney, Reference Cole and Viney2018).
Phylogenetic reconstructions yielded two well-supported monophyletic clades corresponding to C. rudolphii A and C. rudolphii B, and thus confirm the status of the studied taxa. This topology was concordant with the haplotype network results. In all trees, C. eudyptulae and C. septentrionale are sister group to the clade containing [C. rudolphii complex (A, B, C, D, E and F) + C. ogmorhini] and appeared as the most basal species within this clade. These findings agree with results based on protein electrophoretic data (Orrechia et al., Reference Orecchia, Mattiucci, D'Amelio, Paggi, Plotz, Cianchi, Nascetti, Arduino and Bullini1994), 28S (Nadler et al., Reference Nadler, D'Amelio, Fagerholm, Berland and Paggi2000) and mitochondrial markers (Lin et al., Reference Lin, Liu, Zhang, D'Amelio, Zhou, Yuan, Zou, Song and Zhu2012). Within the obtained topologies, strong BI inferences and bootstrap values supported the position of C. ogmorhini (parasite of an otariids) as sister group of the clades of avian parasites. Moreover, the C. osculatum complex, P. cystophorae C. radiatum and C. miroungae (all from phocid hosts) appeared sharing a common ancestor with avian species namely C. microcephalum, C. micropapillatum and C. multipapillatum.
More sequences were available for ITS and cox2 trees construction, and the obtained topology showed a major clade including several sister groups: (C. rudolphii complex + C. variegatum + C. ogmorhini + C. eudyptulae + C. septentrionale) and (C. osculatum complex + P. cystophorae C. radiatum + C. miroungae + C. microcephalum + C. micropapillatum + C. multipapillatum). This topology suggests the paraphyly of the C. rudolphii complex and the avian Contracaecum species, which is consistent with the findings inferred from the nuclear large subunit ribosomal DNA analysis (Nadler et al., Reference Nadler, D'Amelio, Fagerholm, Berland and Paggi2000).
In conclusion, the present results support the paraphyly of the avian Contracaecum species, with in addition C. ogmorhini. This is consistent with previous molecular studies (Nadler et al., Reference Nadler, D'Amelio, Fagerholm, Berland and Paggi2000; Szostakowska and Fagerholm, Reference Szostakowska and Fagerholm2007, Reference Szostakowska and Fagerholm2012; Farjallah et al., Reference Farjallah, Merella, Ingrosso, Rotta, Slimane, Garippa, Said and Busi2008a; Dzido et al., Reference Dzido, Kijewska and Rokicki2012; Lin et al., Reference Lin, Liu, Zhang, D'Amelio, Zhou, Yuan, Zou, Song and Zhu2012; Delgado and Garcia, Reference Delgado and García2015; Szudarek et al., Reference Szudarek, Kanarek and Dabert2017) and allozyme data (Mattiucci et al., Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008). These results were previously proposed by Fagerholm and Gibson (Reference Fagerholm and Gibson1987) on the basis of morphological characters, especially the similarities of caudal papillae patterns, which suggested C. ogmorhini as closely related to the avian Contracaecum, in contrast to Berland (Reference Berland1964), which stated avian Contracaecum as a monophyletic clade.
As previously suggested by Nadler (Reference Nadler2000), these phylogenetic results, based on nuclear and mitochondrial markers, support the urgent need for a nomenclatural revision of these ascaridoid nematodes. This task requires sampling covering the maximum of species, from various hosts and localities, and the analysis of several nuclear and mitochondrial markers that can provide a significant number of characters necessary to clarify the taxonomic status of each group and the phylogenetic relationships between them.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182020001341
Acknowledgement
Thanks to Paolo Briguglio and Merops s.r.l. for making available the stomach contents of the cormorants.
Financial support
The project was financially supported by the Vice Deanship of Research Chairs, Deanship of Scientific Research of the King Saud University.
Conflict of interest
None.
Ethical standards
None.
