New Insights Into the Genetic Variability of Fasciola Hepatica (Trematoda) in Algeria and Relationships with Other Geographic Regions Revealed By Mitochondrial DNA.

This study aims to investigate the level of genetic variability of Fasciola hepatica flukes isolated from cattle in Algeria and to determine the phylogenetic and phylogeographic relationships with sequences isolated worldwide. Mitochondrial (Cytochrome c Oxidase subunit I gene - COI) and nuclear markers (Internal Transcribed Spacers of nuclear ribosomal DNA - ITS) for 24 F. hepatica flukes isolated from 12 cattle in North Algeria were characterised. Only two haplotypes were obtained for the COI gene, resulting in a low level of genetic variation. The analysis of variation among the COI sequences isolated from around the world did not show high levels of genetic divergence, and the phylogenetic analysis revealed a genetic similarity among F. hepatica isolates from different areas of the world. The analysis of the ITS region showed a low level of variability, which prevented obtaining informative phylogenetic and phylogeographic results. The present study also revealed that specimens of F. hepatica are genetically similar in different hosts, indicating that the genetic structure among populations of this parasite is not influenced by the host species. The low levels of genetic variation for COI and ITS regions among fluke isolates from all continents are consistent with a common origin for the flukes' worldwide distribution.

Introduction

Fasciolosis is a parasitic disease caused by two trematode species of the genus Fasciola, F. hepatica Linnaeus, 1758 and F. gigantica Cobbold, 1856, which affect the liver of herbivores and humans. The Fasciola infection has a negative impact on public health (Mas-Coma et al., 2009), and it is also responsible for significant economic losses in livestock production due to the reduction of meat and milk production, slowing growth and reducing fertility. Moreover, infected livers are condemned at meat inspection and removed from sale (Kaplan, 2001). Fasciolosis has a worldwide distribution due to the occurrence of F. hepatica in the temperate zone and F. gigantica in the tropical zone (Mehmood et al., 2017). Both species coexist in subtropical areas (Mas-Coma et al., 2005). Fasciola hepatica is the only endemic species causing fasciolosis in Algeria and represents one of the most common helminths infecting ruminants in this country. However, the occurrence of F. gigantica was recently reported for the first time in sheep from southern Algeria, close to the border with Mali (Chougar et al., 2020). The hosts had a trans-Saharan geographical origin, with introduction from Ghana, through the Sahel countries of Burkina Faso and Mali into Algeria (Chougar et al., 2020).

The prevalence of F. hepatica infection in cattle varies from one Algerian region to another. Infection levels of up to 27.6 % (Ouchene-Khelifi et al., 2018) and 52.4 % (Boucheikhchoukh et al., 2012) were in general reported for the eastern area but in this same region, Chougar et al. (2019) reported prevalence reaching 22.3 %. A lower prevalence (from 6 % to 18 %) was reported in the north-central area (Aissi et al., 2009; Chaouadi et al., 2019), while bovine fasciolosis is less frequent in western and southern Algeria, where prevalences from 2.3 % to 4.4 % (Chougar et al., 2019) and 1.7 % (Ouchene-Khelifi et al., 2018) were reported.

The morphological identification of the two species of Fasciola is based on key characteristics (Periago et al., 2006; Valero et al., 2018): leaf-shaped body, evident shoulders, and oblique body angle for F. hepatica, and indistinct shoulders, parallel margins, and rounded posterior end for F. gigantica. However, several studies have shown the characteristics of the two species overlapping (Ashrafi et al., 2006; Periago et al., 2008) and the existence of intermediated forms in localities where both species coexist, due to phenomenon of hybridisation (Alasaad et al., 2007; Amer et al., 2016). Using molecular tools is considered the best method to provide a clear differentiation between F. hepatica and F. gigantica (Mas-Coma et al., 2009; Itagaki et al., 2005; Le et al., 2008; Shoriki et al., 2016) and contribute to a better understanding of genetic diversity, origin, evolution, and expansion of this parasite (Semyenova et al., 2006; Itagaki et al., 2009; Ai et al., 2011; Thang et al., 2020). Internal transcribed spacers of nuclear ribosomal DNA (ITS) have proven to be useful in performing a proper taxonomic discrimination among Fasciola species. Mitochondrial genes may be highly polymorphic even at an intraspecific level and represent powerful markers to infer the evolutionary processes and phylogenetic relationships within populations of Fasciola species (Itagaki et al., 2005; Semyenova et al., 2006; Cwiklinski et al., 2015).

Some studies have explored the genetic variability of F. hepatica in Algeria (Farjallah et al., 2009; Chougar et al., 2019; Laatamna et al., 2021); and until now, the occurrence of two ITS haplotypes (Chougar et al., 2019), from 7 to 20 Cytochrome c Oxidase subunit I (COI) haplotypes (Chougar et al., 2019; Laatamna et al., 2021) and from 5 to 24 NADH dehydrogenase subunit I (NADH) haplotypes (Chougar et al., 2019; Laatamna et al., 2021) revealed a weak population structuring worldwide for F. hepatica. Haplotypes from Algeria are closely related to those from Europe and Africa (Chougar et al., 2019; Laatamna et al., 2021).

A recent study on Fasciola flukes from several populations located within the Tunisian-Algerian border on combined novel (pepck and pold) and common molecular markers (ITS, COI, NADH and COI trnT-rrnL) revealed a relevant gene flow between Tunisian and Algerian populations of F. hepatica (Amor et al. 2020; Chougar et al., 2020).

The present study aims to investigate the levels of genetic variability among specimens of F. hepatica in Algeria, and it uses two molecular markers with different systems of transmission to the offspring. The first marker is a mitochondrial encoding gene (COI) that does not recombine, is uniparentally inherited, and it is extensively used to depict the phylogeographic patterns of the distribution of species. The second marker is a nuclear noncoding region (ITS) that undergoes a rapid concerted evolution via unequal crossing-over and gene conversion, and it could be highly variable and effective in depicting genetic structuring among groups of species.

Materials and Methods

Sampling

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.

Table 1

Sampling plan. The table reports data on the sampling collection, the GenBank accession numbers of the sequences obtained in the present study, and the haplotypes (Hap. type) found among individuals. The unique allelic variant isolated for the ITS fragment in all the samples analysed in the present study was deposited in GenBank under the accession number MZ292402. The presence of an identical host code for different samples indicates that flukes were isolated from the same host.

Sample code	Host code	Area	Site	Host	Sampling date	GenBank COI #	Hap. type
C1_142	142	Algeria	Tipaza	Cattle	03-29-2016	MT920965	2
C2_142	142	Algeria	Tipaza	Cattle	03-29-2016	MT920966	1
C1_170	170	Algeria	Tipaza	Cattle	04-06-2016	MT920980	1
C1_858	858	Algeria	Blida	Cattle	10-03-2016	MT920969	1
C2_858	858	Algeria	Blida	Cattle	10-03-2016	MT920970	1
C1_903	903	Algeria	Boumerdes	Cattle	10-10-2016	MT920978	1
C1_995	995	Algeria	Algiers	Cattle	11-02-2016	MT920967	1
C2_995	995	Algeria	Algiers	Cattle	11-02-2016	MT920968	2
C1_1000	1000	Algeria	Algiers	Cattle	11-02-2016	MT920974	1
C2_1000	1000	Algeria	Algiers	Cattle	11-02-2016	MT920975	1
C1_1110	1110	Algeria	Boumerdes	Cattle	11-26-2016	MT920982	1
C1_1211	1211	Algeria	Algiers	Cattle	12-31-2016	MT920976	1
C2_1211	1211	Algeria	Algiers	Cattle	12-31-2016	MT920977	1
C1_1215	1215	Algeria	Algiers	Cattle	12-31-2016	MT920981	1
C1_1230	1230	Algeria	Algiers	Cattle	12-31-2016	MT920983	1
C2_1230	1230	Algeria	Algiers	Cattle	12-31-2016	MT920984	1
C3_1230	1230	Algeria	Algiers	Cattle	12-31-2016	MT920985	1
C4_1230	1230	Algeria	Algiers	Cattle	12-31-2016	MT920986	1
C5_1230	1230	Algeria	Algiers	Cattle	12-31-2016	MT920987	1
C6_1230	1230	Algeria	Algiers	Cattle	12-31-2016	MT920988	1
C1_1279	1279	Algeria	Algiers	Cattle	01-18-2017	MT920979	1
C1_1342	1342	Algeria	Algiers	Cattle	02-04-2017	MT920971	1
C2_1342	1342	Algeria	Algiers	Cattle	02-04-2017	MT920972	1
C4_1342	1342	Algeria	Algiers	Cattle	02-04-2017	MT920973	1
FHLAC1	LAC	Italy	Laconi	Cattle	12-09-2013	MT920989	1
FHLAC2	LAC	Italy	Laconi	Cattle	12-09-2013	MT920990	1
FHLAC3	LAC	Italy	Laconi	Cattle	12-09-2013	MT920991	1
FHLAC4	LAC	Italy	Laconi	Cattle	12-09-2013	MT920992	1
FHLAC5	LAC	Italy	Laconi	Cattle	12-09-2013	MT920993	1
FHLAC6	LAC	Italy	Laconi	Cattle	12-09-2013	MT920994	1
FHLAC7	LAC	Italy	Laconi	Cattle	12-09-2013	MT920995	1
FHLAC8	LAC	Italy	Laconi	Cattle	12-09-2013	MT920996	1
FHGIA1	GIA	Italy	Giara di Genoni	Goat	11-28-2013	MT920997	1
FHGIA2	GIA	Italy	Giara di Genoni	Goat	11-28-2013	MT920998	1
FHGIA3	GIA	Italy	Giara di Genoni	Goat	11-28-2013	MT920999	1
FHGIA4	GIA	Italy	Giara di Genoni	Goat	11-28-2013	MT921000	1
FHGIA5	GIA	Italy	Giara di Genoni	Goat	11-28-2013	MT921001	1
FHGIA6	GIA	Italy	Giara di Genoni	Goat	11-28-2013	MT921002	1
FHGIA7	GIA	Italy	Giara di Genoni	Goat	11-28-2013	MT921003	1
FHGIA8	GIA	Italy	Giara di Genoni	Goat	11-28-2013	MT921004	1

Fig. 1

Map of the Algerian sample collection sites. The map shows the geographical origin of the sequences from Algeria isolated in the present study and those from GenBank and BOLD.

DNA extraction, PCR and sequencing

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 MgCl2 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).

Phylogeographic and phylogenetic analyses

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.

Ethical Approval

The manuscript does not contain clinical studies or patient data. Sampling of parasites was not performed on live animals but only on tissues collected post-mortem in a slaughterhouse.

Results

COI

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).

Table 2

Indices of genetic variation. The table reports the estimates of genetic variation for the mitochondrial COI gene dataset. N: sample sizes; bp: fragment size; S: number of polymorphic sites; H: number of haplotypes; hd: haplotype diversity; π: nucleotide diversity.

	N	bp	S	H	hd	π
Samples from Algeria – present study	24	441	1	2	0.159	0.00036
Samples from Algeria – whole dataset	32	441	7	7	0.393	0.00271
Total COI dataset	187	441	42	32	0.753	0.00664

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.

Haplotype network analysis and neutrality tests

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

Median-joining network analysis. The network includes COI sequences from Algeria obtained in the present study along with those from GenBank. The small red plot on one node shows a median vector representing the hypothetical connecting sequence that was calculated using the maximum parsimony method. The number of mutations between sequences that are greater than 1 are reported on network branches. As well, the number of individuals showing the same haplotype that is greater than 1 is reported inside the spot.

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

Median-joining network analysis. The network includes all COI sequences from the present study along with those from GenBank. The small red plots on one node show a median vector representing the hypothetical connecting sequence that was calculated using the maximum parsimony method. The number of mutations between sequences that are greater than 1 are reported on network branches. As well, the number of individuals showing the same haplotype that is greater than 1 isreported inside the spot. The MP calculation post-processing option, that uses only the shortest trees sufficient to generate the graphic output, has been applied for drawing the network. This option allows to obtain a network without showing the reticulations. All Sardinian sequences in the network are form the present study.

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).

Phylogenetic analyses and PCoA

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

Principal coordinates analysis performed on the COI gene dataset. Bi-dimensional plots show the genetic differentiation among populations due to the base differences per site found in the dataset. Percentage of variation explained by the first three axes for the COI dataset: 1st = 40.96, 2nd = 13.80, 3rd=9.31.

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.

ITS

For both Algerian and Sardinian flukes, 32 identical sequences (see Table 1 for details) were obtained for the nuclear ITS fragment in the present study (GenBank accession number: MZ292402). It was not possible to obtain good and scorable sequences for 8 samples, which were not included in the analyses.

A dataset was constructed that included all 32 identical ITS sequences obtained in the present study, along with those (137) from GenBank that exactly matched the ITS fragment used in the present research (see supplementary Table S3 for the GenBank accession numbers and further details). The dataset included 169 sequences of F. hepatica (905 bp) belonging to 14 countries from every continent, except Oceania. Within the dataset, 6 polymorphic sites were retrieved, resulting in a very low level of genetic divergence (hd: 0.058, π: 0.00044) that corresponded to 5 allelic variants. All but 5 sequences (97 %) of the dataset belonged to the most frequent worldwide diffused allelic variant. The only exceptions were represented by 5 sequences from Asia (4 from China and 1 from Iran).

Because of the low level of genetic variability found among ITS variants, neither the network nor phylogenetic tree and neutrality tests analyses were informative (data not shown).

One sequence of F. jacksonii (GB# MN970006) and three of F. gigantica (GB# MW793531, JF432073, MW793533) were included in the ITS dataset for species delimitation analyses. Every species delimitation method showed that all ITS sequences of F. hepatica belonged to a unique, worldwide-distributed taxonomic entity. The only exceptions were represented by two sequences from China included in the dataset, which showed a divergent haplotype belonging to the same taxonomic unit found for F. gigantica.

PCoA was performed on 167 sequences of F. hepatica, excluding the two outlier sequences from China, which likely fell within the variability of the species F. gigantica. PCoA explained 60.08 % of the variability (PCoA1/x-axis: 36.71 %, PCoA2/y-axis: 23.37 %). The results (Fig. 5 and Table S4) showed the occurrence of three genetic groups, including 91 % (Group A), 6 % (Group B), and 3 % (Group C) of the sequences, respectively. A genetic similarity was found along the x-axis between Groups A and B. Group A was the most common and it is present on all continents; individuals of F. hepatica included in this group were isolated from different species of ruminants without a specific pattern of structuring related to the hosts. Group B included sequences isolated from several hosts that spread in the Iberian Peninsula (Spain and Andorra), North Africa (Tunisia), Central America (Mexico), and South America (Bolivia). Group C was exclusive to flukes isolated in cattle from China. No evidence of relations was found between genetic groups and sample collection dates, or host species.

Fig. 5

Principal coordinates analysis performed on the ITS fragment dataset. Bi-dimensional plots show the genetic differentiation among populations due to the base differences per site found in the dataset. Percentage of variation explained by the first three axes for the ITS dataset: 1st = 36.71, 2nd = 23.37, 3rd=11.22.

Discussion

The mitochondrial and nuclear markers used in this study identified all flukes from Algeria as F. hepatica. Although Chougar et al. (2020) recently revealed the presence of F. gigantica in Algerian sheep, present results confirmed the dominance of F. hepatica in Algeria, in accordance with previous molecular studies that found only F. hepatica in this country (Chougar et al., 2019; Farjallah et al., 2009; Laatamna et al., 2021; Farjallah et al., 2013). In accordance with previous studies focused on the genetic variation of this parasite in Algeria (Chougar et al., 2019), the analyses of the COI sequences’ variability showed a generally low level of genetic divergence among individuals in Algeria, with traces of a recent population expansion, as suggested by the Tajima’s D neutrality test results (Fu, 1997). This trend is supported by the low level of genetic variability that was found among sequences that is usually observed in areas recently colonised by these parasites (Mas-Coma et al., 2009; Hewitt, 2000; Robinson & Dalton, 2009). The possible recent introduction of F. hepatica in the Algerian sites is consistent with the fact that this species is native to Europe and expanded its geographic distribution quite recently after Europeans operated a global colonisation with livestock movement over the past five centuries (Mas-Coma et al., 2009).

In this context, the most common COI haplotype found in Algerian isolates (from GenBank and the present study) may correspond to one of the oldest mitochondrial variants present in this country and may be representative of the first lineages introduced from Europe (Walker et al., 2007). The few haplotypes found to be exclusive to single individuals may have recently derived in situ from European founders. They might have differentiated because of silent or neutral nucleotide mutations that originated synonymous codons or amino acids with similar chemical structures.

Interestingly, according to what other authors have already reported (Walker et al., 2007; Elliot et al., 2014), the two new COI haplotypes found in the present study for F. hepatica in Algeria and Sardinia co-occurred in hosts where other different mitochondrial lineages were also present. Consistently, Walker et al. (2007) found several mitochondrial composite PCR–restriction fragment length polymorphism haplotypes in the same cattle, and Elliot et al. (2014) found sheep and cattle with up to ten different mtDNA genotypes. These authors explained their findings by considering the possible occurrence of host infections with diverse fluke metacercariae coexisting in the geographical area where animals usually live or otherwise considering that livestock may have occasionally moved and grazed towards areas where F. hepatica individuals are also characterised by rare mitochondrial haplotypes. For these reasons, plants growing on the edges of wades, rivers, marshes, and irrigation canals could be sources of infection in Algeria (Massot & Senouci-Horr, 1983).

Since a low level of genetic divergence was found for the COI fragment analyzed among isolates on every continent, the present study suggests a common origin of flukes sharing the same haplotype, as it was also suggested by other authors (Le et al., 2000; Lotfy et al., 2008; Amor et al., 2011; Simsek et al., 2011; Martinez-Valladares & Rojo-Vazquez, 2014; Mucheka et al., 2015). Accordingly, a unique taxonomic entity corresponding to the monophyletic clade of F. hepatica, evidenced by phylogenetic tree, further supports the genetic affinity among F. hepatica isolates from different parts of the world. Consistently, Semyenova et al. (2006) also reported low levels of genetic variability of the COI gene in F. hepatica in several countries, with 10 haplotypes found and only 2.3 % of polymorphic sites. Although these authors considered a different fragment of the COI gene from the one the present study analysed, they demonstrated reduced levels of genetic structuring among hosts and geographic regions in Russia, Belarus, Ukraine, Bulgaria, Armenia, Azerbaijan, Georgia, Turkey, Turkmenistan, and China.

The results obtained in the present study also suggest a lack of association between the genetic structuring of the COI gene and host species. Similarly, Santos (2012), who analysed a different and more variable fragment of the COI gene than the present study examined, reported the same haplotype diffused in different hosts (cattle and sheep) from the same geographic Portuguese region. Additionally, Elliot et al. (2014) revealed that many haplotypes are shared between cattle and sheep from Australia, indicating that there is no host selection. Furthermore, present results are also consistent with Bozorgomid et al. (2019), who used the mitochondrial NADH gene to demonstrate low levels of gene flow between Fasciola species isolated from different hosts in Iran (cattle, sheep, and goats), thus suggesting that differences in host species cannot influence the genetic structure of F. hepatica.

The intergenic spacers (ITS1 and ITS2) located between the 18S, 5.8S and 28S rRNA regions generally showed a low level of genetic variability among F. hepatica isolates from almost every continent, confirming that - despite some reported differences in restricted geographical localities - this molecular marker is usually monomorphic within each trematode species (Nolan & Cribb, 2005). For this reason, the ITS molecular marker is rarely used in phylogeographic studies of trematodes, but it is useful and effective in the taxonomic attribution of these parasites (Nolan & Cribb, 2005), being a reliable genetic tool for identifying and differentiating species belonging to the genus Fasciola (Itagaki et al., 2009; Amor et al., 2011; Amer et al., 2016). For both COI and ITS fragments, some F. hepatica isolated from China diverged from the others in the present study, clustering with sequences of F. gigantica. Although the occurrence of a new cryptic species for the genus Fasciola in China cannot be ruled out, this finding more likely suggests the occurrence of hybrids between these species. More generally, introgressed forms of Fasciola are frequently reported in temperate and subtropical regions of eastern Asia based on mitochondrial and nuclear-ribosomal marker identifications (see Le et al. (2008) and references therein).

In conclusion, the present study reports two new mitochondrial COI lineages for F. hepatica identified in cattle from Algeria and Sardinia. The presence of a low number of COI haplotypes among Algerian samples may be the consequence of the recent introduction of a few founders from Europe and the possible occurrence of a high number of clonal parasites, as already reported for other geographic areas (Beesley et al., 2017).

The general low level of genetic variation retrieved for COI and ITS fragments is a frequent genetic pattern of F. hepatica (Beesley et al., 2017) and may suggest a common worldwide origin for this species. Even considering that the short length or the uniparental inheritance (for COI) of the molecular markers might have partially hindered the actual level of genetic variation, this trend can be also explained by taking into account that F. hepatica is a hermaphrodite capable of both cross- and self-fertilisation (although cross-fertilisation is most common), and the occurrence of self-fertilisation may have prompted the loss of genetic diversity in nuclear regions (Cwiklinski et al., 2015). Furthermore, it should be also taken into account that the limited number of sequences analysed for Algeria in the present study might have affected outputs thus introducing a bias due to the occurrence of genetic drift among samples. To solve this matter and better understand the level of genetic variation of F. hepatica in Algeria, additional nuclear (microsatellites) and mitochondrial (whole genomes) genetic data based on a larger sample set would be necessary to depict higher levels of polymorphism and shed new light on the phylogeographic patterns of this species.

Supplementary material