New insights into the genetic variability of Fasciola hepatica (Trematoda) in Algeria and relationships with other geographic regions revealed by mitochondrial DNA

In the present study, 24 individuals of F. hepatica were collected from 12 cattle slaughtered at the Mitidja area in the North-center of Algeria (1 to 6 flukes per each host were isolated) between March 2016 and February 2017 (Table 1 and Fig. 1). Furthermore, 16 individuals of F. hepatica, from one cattle (8 flukes) and one goat (8 flukes) in the Mediterranean island of Sardinia (Italy), in September and November of 2013 (see Table 1 for details), were collected with the aim to enlarge the dataset of isolates used for comparison with Algerian isolates.

Fig. 1

Genomic DNA of the specimens was extracted using the kit Macherey-Nagel NucleoSpin Tissue (MACHEREY-NAGEL GmbH & Co. KG), following the protocol used by Cossu et al. (2015). Sample quality and DNA concentration were determined via spectrophotometry using a NanoDrop™ Lite (NanoDrop Technologies, Thermo Fisher Scientific Inc., Wilmington, DE). The DNA mean concentration obtained for the samples was 75 ng/μL.

PCR amplification of a partial fragment of the COI gene (441 bp) for samples from Algeria and Sardinia were performed using the primers, Ita 8 (forward; 5’-ACGTTGGATCATAAGCGTGT-3’) and Ita 9 (reverse: 5’-CCTCATCCAACATAACCTCT-3’) (Itagaki et al., 2005). Furthermore, PCRs were also performed for a fragment (900 bp) of the nuclear region including ITS-1, 5.8S rDNA, and ITS-2 (ITS) using the primers BD1 (forward: 5′-GTCGTAACAAGGTTTCCGTA-3′) and BD2 (reverse: 5′-TATGCTTAAATTCAGCGGGT-3′) (Luton et al., 1992). All PCRs were carried out in a total volume of 25 μl containing 10 ng of total genomic DNA on average which was combined with 0.6 μM of each primer and one pellet of PuReTaq Ready-To-Go PCR beads (GE Healthcare; 9900 West Innovation Drive, Wauwatosa, WI, USA). Each pellet of PuReTaq Ready-To-Go PCR beads contained reaction buffer, 2.5 units of PuReTaq DNA polymerase, bovine serum albumin (BSA), deoxynucleotide triphosphates (dNTPs) and stabilizers. For each bead reconstituted to a 25 μl final volume, the concentration of each dNTP was 200 μM and of MgCl² was 1.5 mM. The PCR conditions were 4 min at 94 °C as an initial step, followed by 35 cycles of 30 sec at 94 °C, 30 sec at the annealing temperature (56 °C for COI and 57 °C for ITS), and 30 sec at 72 °C, with a final post-treatment of 5 min at 72 °C. Both positive and negative controls were used to test the effectiveness of the PCR protocols, and the absence of possible contamination. The PCR products were visualized on 2 % agarose gels (TAE 1×) and purified by ExoSAP-IT (USB Corporation). Sequencing was performed for both strands using the PCR primers by an external sequencing core service (Macrogen Europe).

The sequences obtained for specimens of F. hepatica from Algeria in the present study were merged with those available for this species in GenBank and Barcode of Life Data system (BOLD) from all over the world (see Fig. S1 in Supplementary Materials for Gen-Bank accession numbers), with the scope to perform a broader phylogeographic analysis of F. hepatica in Algeria and the rest of the world.

The sequences of F. hepatica isolated in Sardinia (Italy) were also included in the analysis, to involve data also from this poorly investigated western Mediterranean island.

Forty contiguous sequences for COI and 32 for the ITS, were aligned and inspected for errors using the package Clustal Omega (Sievers & Higgins, 2014) available at https://www.ebi.ac.uk/Tools/msa/clustalo/) and the data were deposited in the GenBank (see Table 1 for GenBank accession numbers).

The genetic variation within the datasets was assessed estimating the number of polymorphic sites (S), number of haplotypes (H), haplotype diversity (hd), and nucleotide diversity (π) using the software package DnaSP 6.12.03 (Librado & Rozas, 2009).

Median-joining networks (Bandelt et al., 1999) were constructed using the software package Network 10.0.0.0 (www.fluxus-engineering.com) to infer the genetic relationships among haplotypes and allelic variants, thus detecting the possible occurrence of evolutionary forces acting on populations. The transitions and transversions were equally weighed. Due to the lack of knowledge regarding the possible occurrence of retromutation events, the same weight (10) was assigned to all of the observed polymorphisms.

The Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) neutrality tests were performed using the software package DnaSP 6.12.03 (Librado & Rozas, 2009) to infer departures from equilibrium models of the Algerian population. Combining different neutrality tests can help to distinguish among the different evolutionary processes responsible for departures from equilibrium; Fu’s Fs can better detect demographic expansions, whereas Tajima’s D can better detect bottlenecks and populations contractions (Soriano et al., 2008).

The best probabilistic model of sequence evolution was determined using jModeltest 2.1.1 (Posada, 2008), with a maximum likelihood optimised search. The Akaike Information Criterion (AIC) found “TPM3uf+I+G” as the best-fitting model, while the Bayesian Information Criterion (BIC) found the “HKY+G” model. The parameters of the more sophisticated model between the two which were detected were used for input files (i.e. TPM3uf+I+G). Phylogenetic relationships among different taxa (if any) were investigated using a species tree based on Bayesian Inference (BI) by means of the software MrBayes 3.2.7 specifying setting as model parameters: NST = 3, rates = invgamma, ngammacat = 4. Two independent runs, each consisting of four metropolis-coupled MCMC chains (one cold and three heated chains), were run simultaneously for 5,000,000 generations, sampling trees every 1,000 generations. The first 25 % of the 10,000 sampled trees was discarded as burnin. Runs were executed by means of the CIPRES Phylogenetic Portal (Miller et al., 2010). In order to verify the convergence of chains, it was checked that the average standard deviation of split frequencies (ASDSF), approached 0 (Ronquist et al., 2012), and the Potential Scale Reduction Factor (PSRF) was around 1 (Gelman & Rubin, 1992) following Scarpa et al. (2019a).

Phylogenetic trees were visualized and edited using FigTree 1.4.1 (available at http://tree.bio.ed.ac.uk/software/figtree/).

To verify the taxonomic assessment of every sequence in the dataset, four different methods of species delimitation, which are listed below, were used.

The ST-GMYC (Single Threshold-Generalized Mixed Yule Coalescent) method (Pons et al., 2006), which follows the phylogenetic species concept to delimit species, was applied by means of the SPLITS (SPecies LImits by Threshold Statistics) package (Ezard et al., 2009) implemented in the R statistical environment (available at http://r-forge.r-project.org/projects/splits/) on the ultrametric species tree which was obtained by the software Beast 1.10.4 (Drummond & Rambaut, 2007) following Scarpa et al. (2018). The nucleotide divergence threshold (NDT) method was implemented by means of a script (Scarpa et al., 2019b) written in the R statistical environment. For the K/θ method (Birky et al., 2010), used with the corrected formula for sexual organism showed in Birky (2013), clades were selected on the topology of a mid-point rooted Neighbour-joining (NJ) tree (Saitou & Nei, 1987) obtained using the R package APE (Analysis of Phylogenetics and Evolution) (Paradis et al., 2004). Clades showing values of K/θ ≥ 4 should be considered as well-defined entities with a 95 % probability of having an independent evolutionary history. The ASAP (Assemble Species by Automatic Partitioning) method (Puillandre et al., 2020), which is fully exploratory (it does not require any kind of a priori knowledge), was performed using the p-distance model (as substitution model to calculate the distances matrix), selecting default options. Within the list of the best partitions, the species hypothesis, valuating their gap-width score, p-value and threshold distance following Puillandre et al. (2020) were chosen.

On the datasets obtained the principal coordinate analysis (PCoA) was performed using GenAlEX 6.5 (Peakall & Smouse, 2012). This analysis allows to distinguish genetic clusters running on a pairwise genetic distance matrix corrected with K2P (Kimura, 1980) model. The rate of variation among sites was modelled with a gamma distribution and all ambiguous positions were removed for each sequence pair.

Twenty-four sequences of the central portion of the COI gene (441 bp) were obtained for the samples from Algeria in the present study (Table 1). Among them, only one polymorphic site was found that defined two haplotypes (type 1 and type 2, see Table 1 for details) that were shared by 92 % and 8 % of the samples, respectively (see Table 2 for details on the genetic divergence estimates). The two haplotypes diverged from one another for one neutral point-mutation (transition A → G ), which does not affect the protein structure since it produces a change between two non-polar aliphatic amino acids (isoleucine → valine). The mutation occurred at position 799 of the COI gene nucleotide sequence (reference sequence used for the F. hepatica COI gene: NC_002546).

A COI dataset, which included the sequences from Algeria obtained in the present study and those of F. hepatica from Algeria recorded in GenBank (see Fig. 1 for details on the geographic origin of the sequences), was constructed to infer a set of sequences that could likely represents the Algerian population. The dataset showed low levels of genetic variation and it included 32 sequences (24 from the present study and 8 from GenBank), with 7 polymorphic sites that defined 7 haplotypes (see Table 2 for details on genetic divergence estimates).

A further COI dataset, including the sequences from Algeria (24) and Sardinia (16) obtained in the present study and those (147) corresponding to the same portion of the COI gene (441 base pairs of the central fragment of the gene) from F. hepatica strains isolated worldwide and deposited in GenBank and BOLD (see Table S1 for the accession numbers), was also constructed to include the data obtained in the present study on a wider geographic context. This dataset included 187 COI sequences of F. hepatica belonging to 15 countries from every continent except South America (see Supplementary Table S1 for details). Among the sequences, 42 polymorphic sites were retrieved, resulting in a good level of genetic divergence that corresponds to 32 haplotypes (see Table 2 for further details on the genetic divergence estimates). In particular, the sequences from Algeria and Sardinia obtained in the present study belong to the most frequent worldwide diffused haplotypes, except for two never-reported haplotypes found in Algeria (in two flukes from two different cattle, C1_142 from Tipaza and C2_995 from Algiers in Table 1) and one fluke in Sardinia (FHLAC 5 from Laconi in Table 1). Interestingly, these two uncommon lineages were isolated from flukes infecting cattle that were also infected with flukes characterised by the most common COI haplotype.

The network analysis based on sequences from Algeria (Fig. 2) showed a well-defined, star-like shape with a major and highly diffused haplotype that was found in 78 % of sequences and 6 derived haplotypes diverging for 1 to 4 point-mutations. Almost all the derived haplotypes were exclusive to single individuals, except for one lineage that was found in two individuals from Algiers (sample C2_995 in Table 1) and Tipaza (sample C1_142 in Table 1). Overall, 22 of the sequences from Algeria obtained in the present study belonged to the most frequent haplotype of the network.

Fig. 2

The neutrality tests performed on the same dataset of sequences from Algeria showed non-significant negative values of D (-1.604 with 0.10≥ P≥ 0.05) and Fu’s Fs (-3.567 with P≥ 0.10).

The network analysis performed on the COI dataset, including sequences from all over the world (Fig. 3), showed evidence of the occurrence of three highly diffused haplotypes, which are surrounded by many derived lineages that diverged for a single point mutation from the central ancestor and are generally exclusive to single individuals. Two of the three most frequent haplotypes of the network were diffused across nearly every country included in the analysis, while the third most frequent haplotype was exclusive to sequences from Spain, aside from one sequence from Austria. Furthermore, 10.7 % of the haplotypes included in the dataset were exclusive to single individuals, and one haplotype found in three flukes from China was highly divergent (more than 20 point-mutations) from the others.

Fig. 3

The neutrality tests performed on this dataset of sequences showed a significant departure from the equilibrium, with a significant negative value for Tajima’s D test (D=-2.065 with P ≤ 0.05) that is consistent with population expansions and a non-significant negative value for Fu’s Fs (-13.864 with P≥ 0.10).

For phylogenetic inferences, one sequence of F. jacksonii and three sequences of F. gigantica were included within the COI dataset of sequences from around the world as outgroups (see Fig. S1 in Supplementary Materials for GenBank accession numbers).

The phylogenetic tree (Fig. S1 in Supplementary Materials) showed a unique, well-supported monophyletic cluster that included all F. hepatica sequences, except for five sequences from one province in central China (Gansu), which were deposited in GenBank in 2016 and isolated from goats. In particular, two of these sequences set within the large F. hepatica clade in a poorly supported internal sub-cluster, while the three remaining sequences set in an internal well-supported sub-cluster within the F. gigantica clade.

Consistently, every species delimitation method that was used showed that all the COI sequences of F. hepatica included in the dataset belong to a unique, worldwide-distributed taxonomic entity. In accordance with results of the phylogenetic tree, the only exception was represented by three Chinese sequences from Gansu that belonged to the same taxonomic unit found for F. gigantica.

Principal coordinate analysis (PCoA) was performed on all the COI sequences of F. hepatica (n=184), except for the three sequences from China that belong to the F. gigantica variability. PCoA explained 54.76 % of the variability (PCoA1/x-axis: 40.96 %, PCoA2/y-axis: 13.80 %). Results (see Fig. 4 and Table S2) showed the occurrence of three genetic groups, including 59.24 % (Group A), 25.54 % (Group B), and 15.22 % (Group C) of the sequences, respectively, along with one sequence isolated in Algeria from one cattle that set as an outlier outside the three groups.

Fig. 4

Group A mainly spread in western Europe and Africa, with only a few sequences isolated on other continents; individuals of F. hepatica included in this group were isolated from cattle, sheep, and goats without a specific structuring pattern related to the hosts. Group B was almost exclusive to sequences isolated from sheep in Spain, with a few sequences generally isolated in cattle from China, Italy, Austria, and Algeria. Group C was less frequent and scattered across all continents, particularly in flukes isolated from Chinese goats, European cattle and sheep, and African cattle. No association was found between genetic groups and sample collection dates or host species.