Introduction
The Contracaecum osculatum species, members of the family Anisakidae, are mainly parasitic in marine crustaceans and vertebrates. Adult C. osculatum inhabit the alimentary tract of pinnipeds, whereas the larvae (L3) are found primarily in the fish body cavity. Allozyme data analysis resulted in the identification of five C. osculatum sibling species: C. osculatum A, C. osculatum B and C. osculatum C ( = C. osculatum s.s.) (Nascetti et al., Reference Nascetti, Cianchi, Mattiucci, D'Amelio, Orecchia, Paggi, Brattey, Berland, Smith and Bullini1993) from the Arctic-Boreal regions, and C. osculatum D and C. osculatum E from the Antarctic (Orecchia et al., Reference Orecchia, Mattiucci, D'Amelio, Paggi, Plotz, Cianchi, Nascetti, Arduino and Bullini1994). Contracaecum osculatum A and C. osculatum B were reported from the eastern and western parts of the North Atlantic as well as from the North Pacific (Mattiucci & Nascetti, Reference Mattiucci and Nascetti2008). In contrast, C. osculatum C was found mainly in the north-eastern Atlantic (Mattiucci & Nascetti, Reference Mattiucci and Nascetti2008). In addition, C. osculatum C is the only species of the sibling complex present in the Baltic Sea. Due to the morphological similarity of these sibling species, the only appropriate way to differentiate between them, in both the larval and adult stages, is via molecular identification. Allozyme analysis was the first molecular technique successful in identification of the Contracaecum species. The introduction of DNA-based methods has initiated a search for new markers. Anisakids are often identified by means of non-coding parts of the ribosomal DNA (rDNA-ITS), but its relatively high conservatism renders appropriate identification of the closely related species difficult or even impossible. For example, internal transcribed spacer (ITS) sequence analysis cannot be applied to distinguish between C. osculatum D and C. osculatum E because of their identical composition in the ITS1 and ITS2 fragments, as reported by Zhu et al. (Reference Zhu, D'Amelio, Paggi and Gasser2000). Similar difficulties were encountered when identifying other cryptic species, e.g. Pseudoterranova decipiens. Only exceptionally has the analysis of an ITS fragment succeeded in discriminating between species in a sibling species complex (Anisakis simplex complex; D'Amelio et al., Reference D'Amelio, Mathiopoulos, Santos, Pugachev, Webb, Picanco and Paggi2000; Pontes et al., Reference Pontes, D'Amelio, Costa and Paggi2005).
Correct identification of genetically similar species requires the analysis of more variable DNA fragments. Mitochondrial DNA (mtDNA), in which substitutions accumulate much faster than in the nuclear ribosomal DNA, is one of the most promising candidates for a source of species markers. This was borne out by analyses of C. osculatum mitochondrial genes for which Hu et al. (Reference Hu, D'Amelio, Zhu, Paggi and Gasser2001) and Mattiucci et al. (Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008) observed a relatively high variation in the sequences of cytochrome c oxidase subunits I and II (CO1 and CO2) as well as in the rRNA subunits (ssrRNA and lsrRNA). Nevertheless, C. osculatum D and C. osculatum E could be distinguished only by the variation in the oxidase subunit II (CO2) sequence. Studies on mtDNA of other anisakids maturing in marine mammals confirmed the differences in mitochondrial DNA to be relatively higher than those in ribosomal DNA. Cross et al. (Reference Cross, Collins, Campbell, Watts, Chubb, Cunningham, Hatfield and MacKenzie2007) reported relatively high variation within A. simplex s.s. in a large CO1 subunit. Valentini et al. (Reference Valentini, Mattiucci, Bondanelli, Webb, Mignucci-Giannone, Colom-Llavina and Nascetti2006) investigated the genetic divergence in CO2 within the A. simplex complex. A mitochondrial DNA fragment was also successfully used to study differentiation within A. simplex s.s. (Kijewska et al., Reference Kijewska, Dzido and Rokicki2009).
The present study was aimed at comparing selected mtDNA protein-coding genes of C. osculatum A, B and C as potentially useful markers for species identification. The sequences analysed, with their different variation levels, can be useful in studying the complex phylogeny of Contracaecum and in population analyses, including environment-induced parasite–host interactions during the parasite's life cycle.
Materials and methods
Parasites and DNA extraction
The third-stage larvae (L3) of C. osculatum complex were collected from Gadus ogac (Davis Strait) and from Reinhardtius hippoglossoides (Barents Sea). Five specimens of C. osculatum A (vouchers: 27.7, 8.4, 67.3, 27.2, 27.15), C. osculatum B (61.9, 24.3, 27.13, 116.1, 25.10) and C. osculatum C (27.3) were collected from the body cavity of R. hippoglossoides. Two more specimens of C. osculatum C (74a, 105a) were collected from the body cavity of G. ogac. The nematodes were washed in physiological saline and fixed in 70% ethanol. Total genomic DNA was extracted using a modified hexadecyltrimethylammonium bromide (CTAB) protocol (De Jong et al., Reference De Jong, Van Der Wurff, Stam and Olsen1998). Each specimen was placed in 800 μm extraction buffer containing 2 μl of β-mercaptoethanol and 8 μl of proteinase K (20 mg/ml). After incubation (4 h at 56°C) double chloroform extraction (with 700 μl of chloroform each) and alcohol precipitations were carried out to obtain the DNA pellet which was air dried and resuspended in 30 μl of TE buffer (10 mm Tris–HCl, pH 7.5 and 1 mm EDTA).
Species identification and mtDNA amplification
Identification of representatives of the C. osculatum complex was performed based on ITS sequences. The nuclear rDNA regions containing the ITSs (ITS1 and ITS2) were amplified using the primers Nc5 (forward: 5′-GTA GGT GAA CCT GCG GAA GGA TCA TT-3′) and Nc2 (reverse: 5′-TTA GTT TCT TTT CCT CCG CT-3′) (Zhu et al., Reference Zhu, Gasser, Podolska and Chilton1998). Two additional internal primers were used for sequencing: Nc13 (reverse: 5′-ATC GAT GAA GAA CGC AGC-3′; Zhu et al., Reference Zhu, Gasser, Podolska and Chilton1998) and Anc13 (forward: 5′-CTT AGT GCT CAA TGT GTC TG-3′). Polymerase chain reactions (PCR) were performed in a total volume of 50 μl containing 1 × PCR buffer, 1.5 mm MgCl2, 1.25 U of OptiTaq DNA polymerase (EURx, Poland), deoxynucleoside triphosphates (dNTPs; at 300 μm each), and primers at 0.3 μm each. The PCR amplification conditions were as follows: 30 cycles with 30 s denaturation at 94°C (2 min in the first cycle), 30 s annealing at 57°C and 45 s extension at 72°C (7 min in the last cycle). Sequencing was performed by Macrogen Inc. (Seoul, South Korea). The sequences obtained for ITS1–5.8S–ITS2 were compared with those deposited previously in GenBank by Zhu et al. (Reference Zhu, D'Amelio, Paggi and Gasser2000) for species identification. The mtDNA fragments were amplified using AT3Fcon (forward: 5′-TAA GCG CAA AAG ATA TAA AAC AAC TGA C-3′) and RNA2R (reverse: 5′-TTC CCC TAC CTC TAC TTT ACT ACA ACT TAC TC-3′) (Kijewska et al., Reference Kijewska, Dzido and Rokicki2009) primers. For amplification Expand Long Range, dNTPack was used (Roche, Mannheim, Germany): 1 × PCR buffer 2, 3.75 U/50 μl of polymerase, 1.25 mm MgCl2, 3% dimethyl sulphoxide (DMSO), 500 μm dNTPs and 0.3 μm primers. PCR was performed under the following conditions: 10 cycles of 10 s denaturation at 92°C, 15 s annealing at 47°C and 12 min at 68°C, and 20 cycles of 10 s at 92°C, 15 s at 47°C and 12 min 30 s at 68°C. Additionally there was 2 min denaturation in the first cycle and 7 min annealing in the last cycle. PCR products were used for sequencing of C. osculatum complex genes using the primers given in table 1.
Table 1 List of genes analysed, with primers (name, sequence and orientation) used for sequencing, and the GenBank accession numbers of the Arctic Contracaecum osculatum species sequences obtained.
aO, Primer orientation (F, forward; R, reverse); bKijewska et al. (Reference Kijewska, Dzido and Rokicki2009).
Sequence analysis
Sequences were aligned in the same reading frame using ClustalW software (Thompson et al., Reference Thompson, Higgins and Gibson1994) and cut to an unambiguous length (table 2). Stop codons were removed from the analysis. Genetic nucleotide divergence among the Contracaecum taxa was estimated by the Kimura-2-parameter (K2P) method and by p-distance values. The pairwise distance calculation, nucleotide frequencies, ratios of transitions/transversions (Ts/Tv) and variable nucleotides at codon positions were conducted in MEGA5 (Tamura et al., Reference Tamura, Peterson, Peterson, Stecher, Nei and Kumar2011). Calculations of dN/dS were performed using KaKs_Calculator 2.0 (Zhang et al., Reference Zhang, Li, Zhao, Wang, Wong and Yu2006) under the GY model which takes account of sequence evolutionary features, such as transition/transversion rate ratio and nucleotide frequencies. Tajima's D tests were calculated in the DnaSP program (Librado & Rozas, Reference Librado and Rozas2009).
Table 2 The overall intraspecific differences in nucleotide (nt) and amino acid (aa) sequences of mtDNA protein-coding genes of the Arctic Contracaecum osculatum species.
ant, Nucleotide; aa, amino acid; bTs, transition; Tv, transversion.
Results and discussion
The mtDNA fragments (~ 11,550 bp) of three Arctic C. osculatum species were amplified and used as templates for sequencing. No size variations of amplicons were observed. Sequencing was performed using ten primers, 1–3 for each gene (table 1). The sequences of eight mtDNA protein-coding genes were deposited in GenBank (table 1). The length of the fragments varied from 231 to 1104 bp (without stop codons) (table 2). The set of fragments consisted of complete sequences of CYTB, ND2, ND3, ND4L and ND6, and partial sequences of CO1 (42%), CO3 (73%, N-terminus) and ND5 (42%, C-terminus).
All the genes studied proved AT-rich. The AT content in the first, second and third position averaged 67.8, 70.9 and 72.8%, respectively (table 3). The overall AT content across all the genes examined was 70.5%, and was very similar to that (69.5%) reported by Kim et al. (Reference Kim, Eom and Park2006) for protein-coding genes from A. simplex. The highest AT contents were recorded in ND2, ND6 and ND4L in all the positions, and the lowest contents were observed in CYTB, CO3 and CO1. The sequence data were skewed towards thymine (50.3% in all the genes studied). The highest thymine contents were observed in ND2, ND4L and ND6 (table 3). The high AT content did not reflect the variability in nucleotides or amino acids in the genes. The transition (A ↔ G or C ↔ T) to transversion (A ↔ C, A ↔ T, C ↔ G, G ↔ T) ratios ranged from 4.9 in ND3 to 11.8 in ND6 (table 2) and were at a level typical of mitochondrial DNA reported in other studies (e.g. CO1; Valentini et al., Reference Valentini, Mattiucci, Bondanelli, Webb, Mignucci-Giannone, Colom-Llavina and Nascetti2006).
Table 3 Nucleotide composition of mtDNA protein-coding genes in the Arctic Contracaecum osculatum species (C. osculatum A, C. osculatum B and C. osculatum C) arranged by codon position.
The pairwise comparison of the species revealed K2P to range within 0.06–0.12 (table 4). The highest pairwise divergence in all the genes (4830 bp) occurred between C. osculatum C and the two other species, the lowest divergence being recorded between C. osculatum A and C. osculatum B (the B versus C and A versus C differences were almost equal). The CO3 and ND3 genes proved most different, but the difference between two pairs of species never exceeded 0.02 (table 4). Zhu et al. (Reference Zhu, D'Amelio, Paggi and Gasser2000) reported a similar pattern in the ribosomal DNA they analysed in C. osculatum: the A versus B pairwise difference was lower than that between C and A, and between C and B. The overall intraspecific differentiation observed in this study for C. osculatum A and C. osculatum B was similar (K2P = 0.02), whereas the sequence polymorphism calculated for C. osculatum C (K2P = 0.03) was significantly higher (P < 0.001). The hypothesis that the sequence variation found in this study would be due to the geographic origin of the individuals examined (C. osculatum C collected from the north-east and north-west Atlantic) should be treated with caution. We observed an inconsistency between the data yielded by our set of CO1 sequences and a set of C. osculatum collected by Hu et al. (Reference Hu, D'Amelio, Zhu, Paggi and Gasser2001). In the case of C. osculatum C, the sequence variation calculated for the Hu et al. (Reference Hu, D'Amelio, Zhu, Paggi and Gasser2001) sample (from the Baltic Sea) was slightly lower than that in this study (0.03 versus 0.02). In contrast, a comparison of C. osculatum B samples revealed intraspecific variation (0.02 versus 0.03) which was lower than that produced by the set of sequences analysed by Hu et al. (Reference Hu, D'Amelio, Zhu, Paggi and Gasser2001) (from the Canadian Atlantic).
Table 4 Nucleotide differences (Kimura-2-parameter below the diagonal and p-distance above the diagonal) per sequence between sequences. Calculations were based on pairwise comparisons of the Arctic Contracaecum osculatum species (CosA, C. osculatum A; CosB, C. osculatum B; CosC, C. osculatum C).
The K2P and p-distance were at their highest in the CYTB gene (K2P = 0.11, table 4). The intraspecific variation (fig. 1) observed in this gene is three and six times lower than the interspecific differences between C. osculatum C (0.03) and C. osculatum A and B (0.02). The comparison of the taxa revealed 101 variable nucleotide positions in the CYTB sequence (35, 40 and 38 fixed nucleotide changes typical only of C. osculatum A, B and C, respectively). This nucleotide variation represented the domination of transitions (n = 76) over transversions (n = 21) and multiple substitution events (n = 6). There were also six fixed amino acid changes which could be species-specific (fig. 2). CO3 was the other variable gene (K2P = 0.10 and 0.12), its intraspecific variation being at a level comparable to that observed in CYTB. There were 11, 17 and 13 fixed nucleotide changes typical only of C. osculatum A, B and C, respectively. These nucleotide variations represented transitions (n = 34) and transversions (n = 7). There were five variable species-specific amino acids (highlighted in fig. 2). The overall differentiation in ND5 (7%, table 2) was slightly lower than in CO3 and CYTB (7.5 and 7.6%, respectively). The observed pairwise interspecific K2P values varied from 0.09 to 0.11 in ND5. This gene showed 15, 20 and 24 nucleotide changes specific only for C. osculatum A, B and C, respectively. These nucleotide variations represented transitions (n = 47), transversions (n = 11) and a single multiple substitution event. The CO1 nucleotide distances obtained in this study were somewhat lower than the distance observed for CO3 (table 2) and similar to that found for CO2 (Mattiucci et al., Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008). In CO2, all the nucleotide substitutions were silent (Mattiucci et al., Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008). In this study, CO1 showed one amino acid substitution which seemed to be C. osculatum C specific (isoleucine in position 131 instead of valine) (fig. 2). The relatively low amino acid differentiations in the cytochrome c oxidation genes (CO1 and CO2) confirmed that most of the nucleotide changes in these genes are silent. A different situation was observed in ND6, the overall nucleotide differentiation (6%) was translated into amino acid differentiation at the level of 1.2%. The gene showed four fixed amino acid changes. A similar pattern was observed in ND3 which showed overall 2.1% amino acid variation with four fixed amino acid changes (fig. 2). Relatively low interspecific differentiation was observed in ND2 and ND4L (table 4). There were no fixed amino acid changes in ND4L and only 17 species-specific nucleotide positions. Only one amino acid change (position 72 in a single individual of C. osculatum B) was recorded in ND2.
Fig. 1 Interspecific nucleotide divergence (Kimura-2-parameter) in mtDNA protein-coding genes of three Arctic Contracaecum osculatum species (diagonal lines): C. osculatum A (grey), C. osculatum B (white) and C. osculatum C (black). The plot is based on 13 individuals examined.
Fig. 2 Characteristic amino acid substitutions in mtDNA protein-coding genes of the Arctic Contracaecum osculatum species: C. osculatum A (grey), C. osculatum B (black) and C. osculatum C (bold black). The number indicates the position in each protein; the amino acid change is shown in brackets.
The overall variation in the C. osculatum genes examined was four times that observed in the ITS and rDNA. The average ω (dN/dS) was estimated at 0.008, whereas species-specific estimates ranged from 0.010 in C. osculatum B to 0.021 in C. osculatum A. At the interspecific level, the estimated values of ω were lower (0.005–0.007) than at the intraspecific level. The ω levels suggest a strong purifying selection. The variation at synonymous and non-synonymous sites within (polymorphism) and between species (divergence) showed a substantial difference between the intraspecific and interspecific levels, the dN/dS ratios being lower at the intraspecific level. The homogeneity of the observed values of ω in each species demonstrated their dN and dS distributions to be similar – evidence of a similar selection pressure.
As described above, the highest differences in DNA sequences between the Arctic C. osculatum species were observed in CYTB, CO3 and ND5 (table 4). The interspecific dS and dN levels for these three genes was higher than the intraspecific ones. Of all the genes tested in this study, the three mentioned above appear to be the best candidates for species-specific markers in the Arctic C. osculatum complex, and probably in the Antarctic species (C. osculatum D and C. osculatum E) as well. The higher interspecific differentiation found in the genes mentioned above is very promising from the standpoint of identifying the two Antarctic species, and seems to be even better than the interspecific differentiation in CO1, lsrRNA and ssrRNA (Hu et al., Reference Hu, D'Amelio, Zhu, Paggi and Gasser2001), and perhaps CO2 (Mattiucci et al., Reference Mattiucci, Paoletti, Webb, Sardella, Timi, Berland and Nascetti2008). So far, CO2 has been – except for allozymes – the only marker with which to distinguish successfully between the Antarctic C. osculatum species.
The data collected for C. osculatum C confirmed the taxon to be most different among the entire Contracaecum complex. Due to the small sample size, interpretation of results obtained in this study was very difficult. However, considering the already published data on the phylogeny of C. osculatum C (Mattiucci & Nascetti, Reference Mattiucci and Nascetti2008), the higher K2P observed in C. osculatum C could, in our opinion, be a result of the effective total population consisting of a number of different geographic subpopulations. This conjecture is supported by a number of findings. The species diverged during the Pleistocene refuge (the Baltic Sea), whereas the other species of the complex started diverging when the complex was distributed over both polar regions (Mattiucci & Nascetti, Reference Mattiucci and Nascetti2008). Contracaecum osculatum C is characterized by a narrow specificity and relatively limited distribution, compared to other Arctic species, which contributes to the geographical separation of the population, particularly in regions such as the Baltic Sea, with a limited distribution of definitive hosts. Tajima's D test (Tajima, Reference Tajima1989) of CO1 sequences deposited in GenBank yielded positive values for C. osculatum C only. This, according to Anderson et al. (Reference Anderson, Blouin and Beech1998), could confirm the effect of an admixture of a subpopulation.
The results obtained showed the within-species variation to be lower than the interspecific variation. In our opinion, increasing the number of individuals examined will not affect the observed intraspecific polymorphism and interspecific variability in a way that prevents species identification. On the other hand, however, increasing the number of individuals examined will allow a more comprehensive interpretation of evolutionary processes.
Acknowledgements
We are grateful to anonymous referees for their valuable suggestions and comments. The research was supported by the University of Gdańsk grant no. L115-5-0102-9 and the Polish Ministry of Education and Science grant no. 1182/IPY/2007/01. We would like to thank Dr Lone Nukaaraq Møller for parasitic nematode material from Gadus ogac.
