Mitochondrial sequence diversity reveals the hybrid origin of invasive gibel carp (Carassius gibelio) populations in Hungary

Sample collection and DNA extraction

One hundred and thirty-two gibel carp samples were collected from six locations in Hungary including Lake Balaton (Siófok; n = 29; N46°54′24 E18°02′41), two reservoirs of Kis-Balaton Water Protection System (KBWPS) I stage (n = 17; N46°36′02 E17°09′01) and II (n = 18; N46°39′47 E17°07′23), Hőgyész (n = 30; N46°28′34 E18°26′02), Siófok-Töreki fish pond system (n = 19; N46°52′32 E18°00′14) and Őszödi-berek wetland (n = 19; N46°49′02 E17°48′12) (Fig. 1). Fish collection for laboratory examinations was authorized by the Government Office of Pest county (Permit no.: XIV-I-001/2302-4/2012). Fishes were collected with electrofishing and from the eel trap of Balaton Fish Management Nonprofit Ltd. at Siófok. Since gibel carp is an invasive fish species, all collected individuals were euthanized with clove oil before tissue samples were collected for DNA sequencing. Tissue samplings have been authorized by the Minister of Agriculture (Permit no.: HHgF/122-1/2018). Tissue samples were collected from the caudal fin and then stored in 96% ethanol at −20 C° before use. DNA was isolated from the tissue samples by E.Z.N.A Tissue DNA Kit (Omega Bio-tek, Norcross, GA, USA) following the producer’s protocol. DNA concentration was measured by spectrophotometer (IMPLEN, NanoPhotometer™ Uv/Vis) and the quality was tested by running 250 ng of each DNA sample on 1.5% agarose gel.

PCR amplification and sequences

The control region of the mitochondrial genome (D-loop) was amplified with primers from the Cyprinidae family, Carp-pro2-F (5′-TCACCCCTGGCTCCCAAAGC-3′) and Carp-phe2-R (5′-CTAGGACTCATCTTAGCATCTTCAGTG-3′) (Wang et al., 2010). PCR reaction final volume was 25 µl, contained 1  × PCR buffer with (NH₄)₂SO₄ (Fermentas; Thermo Fisher Scientific, Waltham, MA, USA), 2000 µM dNTP mix, 250 nM for each primer, 1.5 mM MgCl₂, 100 ng template and 1 U Taq polymerase (Fermentas). To collect more information about the populations, the resulting haplotypes were analysed by using the same protocol with primers from the Cytochrome c oxidase I gene (CoI): CO1_FF2d_F (5′-TTCTCCACCAACCACAARGAYATYGG-3′), CO1_FR1d_R (5′-CACCTCAGGGTGTCCGAARAAYCARAA-3′) (Ivanova et al., 2007) and with primers from the Cytochrome b (Cytb): Cytb_H_2_R (5′-GTTTGTTTTCTAACCCGATCAATG-3′) Cytbas_F (5′-GAAGGCGGTCATCATAACTAG-3′) (Xiao, Zhang & Liu, 2001). PCR temperature profiles were the following for all three markers: Denaturation at 95 °C for 2 min., then 30 s. at 94 °C, 20 s. at 52 °C and 1 min. at 72 °C for 35 cycles. Final elongation was at 72 °C for 10 min. PCR product quality was assessed on a 1.5% agarose gel, then purified by NucleoSpin Gel and PCR Clean-up kit (Macherey-Nagel, Germany). The purified product was sequenced from both ends using a Big Dye Terminator v. 3.1 Cycle Sequencing kit (Applied Biosystem), based on the method of sanger sequencing with an ABI 3130 Genetic Analyzer machine.

Bioinformatic and statistical analysis

Chromatograms were converted to FASTA file format with BioEdit Sequence Alignment Editor (Hall, 1999). The 699 bp long sequences of the D-loop region were aligned and analysed by the MegaX software (Kumar et al., 2018) with the ClustalW algorithm. Haplotypes diversity was calculated based on the polymorphic sites and haplotypes by using the DnaSP version 5 (Librado & Rozas, 2009) software. Pairwise FST values were calculated by using the MegaX software (Kumar et al., 2018). Different haplotypes were checked and compared with outgroup sequences from the National Centre for Biotechnology Information (NCBI GenBank) standard database (https://www.ncbi.nlm.nih.gov/) by nucleotide BLAST (Basic Local Alignment Search Tools). CoI results were compared to the sequences of the BOLD (Barcode of Life Data System) system (http://v3.boldsystems.org/). The network of the D-loop haplotypes was built using a median joining algorithm by using PopART software (Rohl, Bandelt & Forster, 1999; Clement et al., 2002). We used the closest related species (C. auratus, C. a. buergeri, Carassius carassius, C. cuvieri, Cyprinus carpio carpio) as reference sequences along with other C. gibelio haplotypes found in the literature. In case of Cytb based on the genetic distance and the phylogenetic tree was reconstructed by MegaX software with Neighbour Joining fitting with Kimura-2 parameter and 1,000 bootstrap replication. The references were the haplotypes described in C. gibelio neotypes by Kalous et al. (2012), as well as European sequences described by Takada et al. in the Carassius auratus-complex, where a European clade was identified based on the Cytb sequence (Takada et al., 2010).

D-loop

Based on the comparisons of 699 bp long sequences of the D-loop region, 22 haplotypes were identified in 132 individuals (Table 1). Three of these (HapDl_1, HapDl_7, HapDl_21) were identified in more than 15 individuals while the others were less frequent. These haplotypes were present in maximum 6 individuals. Within the haplotypes, the number of polymorphic sites was 43 (Fig. 2). For each population, the highest haplotype diversity was in Siófok 0,83 (HD) ± 0,04 (SD). The second was the KWBPS I with 0,80 (HD) ± 0,08 (SD). The lowest value was 0,19 (HD) ± 0,11 (SD) for the population of Őszödi-berek. In the population of Hőgyész HD was 0,72 ± 0,05 (SD), in the population of Siófok-Törek 0,66 ± 0,07 and in the population of KWBPS II. HD was 0,63 ± 0,09 (SD). Thirty-five sites were parsimony informative, based on the DnaSP analyses. Samples from all locations represented a mixture of different haplotypes except those for the sequences from Őszödi-berek, which separated into two individual groups (HapDl_21, HapDl_22) and the HapDl_2 group with seven individuals from Lake Balaton. Twelve other haplotypes were found (HapDl_3, HapDl_5, HapDl_6, HapDl_10, HapDl_11, HapDl_12, HapDl_14, HapDl_15, HapDl_17, HapDl_18, HapDl_19, HapDl_20), of which each contained only one or two samples from one location. The identified haplotypes were divided into three major groups separated on the network figure (Fig. 3). Most of the haplotypes, included the two largest groups (HapDl_1, HapDl_7), clustered together in the first group (Dl_Group 1.). The second group (Dl_Group 2.) contained five haplotypes (HapDl_3, HapDl_8, HapDl_11, HapDl_14, HapDl_20) and the third group contained only one haplotype (HapDl_4) closer to the Carassius auratus buergeri (AB377291.1) sequence, than the other Carassius gibelio haplotypes did. The rest of the outgroup sequences (AC.: LC019787.1, JN117597.1, JX122531.1, JQ390593.1) were separated according to the genetic distance.

Six groups (HapDl_3, HapDl_4, HapDl_8, HapDl_9, HapDl_11, HapDl_15) with the total of 12 individuals were identified differently as gibel carp (Table 1). Two of them (HapDl_11, HapDl_15) contained only samples from Hőgyész and were identified as C. auratus. Samples from haplotype 8, 9, and 3 were identified as C. auratus too, along with samples from Lake Balaton and KBWPS Stage I. Haplotype 4 has shown the highest similarity with C. a. buergeri instead of gibel carps. Samples in this group were from Lake Balaton and the ponds of Siófok-Töreki. KBWPS Stage II. was the sole population with only C. gibelio haplotypes (EF633617.1, MF083605.1, MF036180.1, MF036179.1). Lake Balaton samples have the highest genetic diversity of nine haplotypes. This was followed by the samples from KBWPS Stage I. with eight haplotypes (Table 1). Őszödi-berek had only two haplotypes (HapDl_ 21, HapDl_22). The two D-loop sequence variants, detected in the case of 19 individuals characterized only this population.

Pairwise FST value (Table 2) among the populations has shown a moderate difference between each pair but most of the values were lower than 0.1, with the exception of one population (Őszödi-berek) which resulted in higher values. The lowest difference (−0.015) was identified between the populations of Hőgyész and KBWPS Stage I.

Cytochrome b

Cytb sequences were analysed to confirm the results of 22 D-loop haplotypes. Only six Cytb haplotypes were identified. The complete length of the sequence region was 1057 bp. Within the haplotypes the numbers of polymorphic sites were 27 (Fig. 4). NCBI database was used to check the evolutionary origin of the sequences. All six haplotypes were detected as Carassius gibelio (Table 3). None of the first 100 listed results are recognized as C. auratus, C. a. buergeri or other CAC species. Based on the previous work of Kalous et al. (2012) Cytb haplotypes were compared to the neotype sequences they described (AC: HM000009, HM0000020, GU170378, FJ822041, FJ478019, Ab368700, HM000008, HM008678, JN402305, HM008684, HM008685.1, DQ868924, DQ868925, DQ868926, HM008690) and the sequences described by Takada et al. (2010) (DQ399926.1., DQ399929.1) and were presented together on the phylogenetic tree (Fig. 5). Two of the six Hungarian haplotypes (HapCb_1, HapCb_4) were integrated into the first neotype group with other European sequences and two haplotypes to neotype II. (HapCb_3, HapCb_5) which contain only reference sequences from Mongolia. The third group on the top of the three contained the sequences described by Takada et al. (2010), and the Hungarian HapCb_6 and HapCb_2.

Cytochrome c oxidase 1

CoI gene sequences of the 22 D-loop haplotypes were also determined. The length of the analysed fragment was 562 bp long. Four haplotypes (Fig. 6) were identified. Within the haplotypes the number of polymorphic sites was three. Sequences were analysed by NCBI BLAST and BOLD system. According to the BLAST standard nucleotide database two groups (HapCoI_2 and HapCoI_4) out of four showed the highest similarity with Tatia intermedia, one with Cyprinus carpio haematopterus (HapCoI_1) and one with C. auratus (HapCoI_3) (Table 4).

The group containing most of the samples (HapCoI_2) had 100% of identity coverage of the Tatia intermedia query sequence. Furthermore, based on the BOLD system HapCoI_2 and HapCoI_4 were identified as C. gibelio. HapCoI_1, as well as HapCoI_3 have shown the highest similarity with C. a. auratus. In the case of HapCoI_3, not even the first 20 hits contained C. gibelio sequences. The phylogenetic tree formed three clearly visible branches (Fig. 7). The first branch included the HapCoI_1 and HapCoI_3. The second branch contained the HapCoI_2 and the HapCoI_4 was in the third group.