Phylogenetic analysis of Fasciola spp. isolated from slaughtered cattle in KwaZulu-Natal and Mpumalanga provinces of South Africa based on the cytochrome c oxidase subunit I mitochondrial marker

Fasciola spp. are the causative agents of fascioliasis in humans and livestock. Before the development of control and management measures, the geographical distribution of the species and patterns of infection must be considered. Because of difficulties in the phenotypic differentiation and morphometric classification of Fasciola spp., DNA molecular markers have become more useful for fluke differentiation and description of phylogenetic patterns. This study aimed to differentiate and describe the phylogenetic background of Fasciola spp. isolated from cattle slaughtered at three abattoirs in the Mpumalanga and KwaZulu-Natal provinces of South Africa. The cytochrome c oxidase I (COI) – FHCO1 (forward: 5′-TTGGTTTTTTGGGCATCCT-3′) and FHCO1 (reverse: 5′ -AGGCCACCACCAAATAAAAGA3′) – marker was sequenced from 55 Fasciola flukes that were collected from abattoirs in catchment areas of the KwaZulu-Natal and Mpumalanga provinces. Fasciola hepatica was demonstrated to have 100% prevalence in KwaZulu-Natal and Mpumalanga (highveld), respectively, and 76% prevalence in the lowveld (Belfast area) of Mpumalanga. Two animals from the Belfast metapopulation were co-infected with both Fasciola gigantica and F. hepatica. DNA sequence analysis of all the isolates demonstrated a sequence conservation of 0.472, nucleotide diversity of 0.082 and Tajima’s D of -1.100; however, it was not statistically significant (p > 0.05). Twenty-two haplotypes were identified, with 18 novel haplotypes being unique to the isolates from South Africa. Within the study samples, 12 haplotypes were isolated to a few individuals, with a haplotype diversity of 0.8957 indicating high genetic diversity. Principal coordinate analysis supported the clustering and distribution of the haplotypes, with 11.38% of the variation being attributed to coordinate 2 and 55.52% to coordinate 1. The distribution of Fasciola spp. has been demonstrated to be related to the distribution of the freshwater intermediate host snails, Lymnaea spp., as well as the relative altitude of the localities in South Africa. Information provided by this study serves as preliminary evidence for further studies on the mapping of the distribution of F. gigantica and F. hepatica in South Africa, which is key in designing control programmes for fascioliasis in humans and livestock.


Introduction
Fascioliasis is a food-borne parasitic disease with a near global distribution and dire socioeconomic effects because of human and livestock mortalities, retarded growth, reduced meat and milk production, condemnation of infected livers, animal emaciation and cost of anthelmintics for treatment (Mas-Coma et al. 2004;Mas-Coma, Valero & Bargues 2009b;Spithill & Dalton 1998). The disease occurs as a result of infection by digenean trematode species of the genus Fasciola. Several species have been described under the genus, with the taxonomic standing still under scrutiny and debate (Agatsuma et al. 2000;Ai et al. 2011;Hashimoto et al. 1997;). However, the main causative agents of fascioliasis in livestock and humans are Fasciola hepatica and Fasciola gigantica, with hybrids being reported in areas where the distribution of the two species overlap (Ai et al. 2011;Ashrafi 2015;Walker et al. 2008;WHO 1995). Studies have reported the presence of F. hepatica in Africa, America, Asia, Europe and tropical islands of the Pacific (Mas-Coma, Bargues & Valero 2005;Mas-Coma et al. 2009b;Mucheka et al. 2015;Nguyen et al. 2009;Walker et al. 2008). Fasciola gigantica, on the contrary, has a limited geographical distribution and has been reported in Africa and Asia and some parts of southern Europe (Mas-Coma & Bargues 1997). The two Fasciola spp. are the causative agents of fascioliasis in humans and livestock. Before the development of control and management measures, the geographical distribution of the species and patterns of infection must be considered. Because of difficulties in the phenotypic differentiation and morphometric classification of Fasciola spp., DNA molecular markers have become more useful for fluke differentiation and description of phylogenetic patterns. This study aimed to differentiate and describe the phylogenetic background of Fasciola spp. isolated from cattle slaughtered at three abattoirs in the Mpumalanga and KwaZulu-Natal provinces of South Africa. The cytochrome c oxidase I (COI) -FHCO1 (forward: 5′-TTGG TTTTTTGGGCATCCT-3′) and FHCO1 (reverse: 5′ -AGGCCACCACCAAATAAAAGA3′) -marker was sequenced from 55 Fasciola flukes that were collected from abattoirs in catchment areas of the KwaZulu-Natal and Mpumalanga provinces. Fasciola hepatica was demonstrated to have trematode species are known to mainly infect humans and ruminant animals, resulting in losses in livestock production, fertility and resistance to draught in ruminants (Bernardo et al. 2011;Keiser & Utzinger 2007). Economic losses in livestock have been reported to exceed $3 billion annually, with the zoonotic aspect of the disease also making it a public health concern (Mas- Coma et al. 2009a;Spithill & Dalton 1998).
Fasciola hepatica and F. gigantica exhibit similar life cycles (Graczyk & Fried 1999;Mas-Coma & Bargues 1997). The diheteroxenous life cycle of Fasciola (Andrews 1999) is dependent on Lymnaeid freshwater snails as the intermediate hosts for both F. hepatica and F. gigantica (Bargues et al. 2001). Lymnaeids have been reported to act as intermediate hosts to numerous other trematode species depending on the geographical region, ecology and snail-parasite specificity (Bargues et al. 2001;De Kock, Wolmarans & Bornman 2003;Mas-Coma et al. 2005). The life cycle of Fasciola includes an intermediate host snail and a definitive host with complex reproductive biology, which results in high gene flow and genetic variability within Fasciola spp. populations (Cwiklinski et al. 2015). Studies have shown that multiple Lymnaeid snail species may be susceptible to Fasciola (De Kock, Joubert & Pretorius 1989;Mas-Coma et al. 2009) and this may also contribute to potential population sub-structuring as well as development of distinct genetic clusters between different geographical localities (Beesley et al. 2017;Vilas, Vázquez-Prieto & Paniagua 2012).
Over the last few decades, reports of fascioliasis have increased, resulting in the increasing need for taxonomic clarity as well as differentiation of the two fasciolid species. Fasciola gigantica and F. hepatica overlap in their geographical distribution (Mucheka et al. 2015;Robles-Pérez et al. 2015), and this has resulted in numerous discussions and notions on the taxonomic patterns and identity of Fasciola spp. as well as observed intermediate forms (Agatsuma et al. 2000;Ai et al. 2011;Nguyen et al. 2009;Periago et al. 2007). Morphological techniques, primarily used in diagnosis, take their basis on the differentiation of morphological characteristics of eggs and the adult fluke (Cuomo, Noel & White 2009;Mas-Coma et al. 2005;Thanh 2012). However, these techniques have limited sensitivity and specificity as they require skilled taxonomists who are scarcely available in many parts of the world such as South Africa. There have also been several reports of hybrids of F. gigantica and F. hepatica, which may not be morphologically differentiated from the main species Schweizer et al. 2007). With the increase in reports of fascioliasis and the rapid emergence of resistance to anthelmintics, further genetic analysis studies are therefore warranted. Resistance to triclabendazole to the recommended treatment dose has been reported from several parts of Eurasia (Brennan et al. 2007;Brockwell et al. 2014;Fairweather 2009;Kelley et al. 2016;Peng et al. 2009;Vilas et al. 2012), with a report of potential resistance of F. hepatica to treatment in a human case in the Netherlands (Winkelhagen et al. 2012).
Treatment and control of fascioliasis largely depends on early diagnosis of infection, and with advances in DNA technology, application of molecular markers for species identification and genetic characterisation has become necessary in studying Fasciola isolates (Cwiklinski et al. 2015;Elliott et al. 2014;Mucheka et al. 2015;).
Studies have shown the utility of nuclear ribosomal markers such as ITS 1 and ITS 2 in species identification; however, mitochondrial markers are highly variable and can resolve taxonomic patterns of more closely related species and/or populations (Mucheka et al. 2015;Patwardhan, Ray & Roy 2014). The cytochrome c oxidase subunit I (COI) and nicotin adenin dinucleotide dehydrogenase (NADH dehydrogenase) subunit I (NadI) have been used as molecular markers for elucidating the phylogenetic relationships of Fasciola spp. and describing to an extent the genetic diversity of Fasciola populations (Ai et al. 2011;Hashimoto et al. 1997;Mucheka et al. 2015;Peng et al. 2009;Semyenova et al. 2006). The COI marker is widely regarded as an efficient DNA barcode and has been commonly used in species identification and differentiation of Fasciola spp. in recent studies in South Africa and Zimbabwe (Mucheka et al. 2015).
Against this background, the present study was aimed at identifying Fasciola spp. isolates collected from cattle slaughtered at abattoirs located in KwaZulu-Natal (KZN) and Mpumalanga (MP) provinces using the (mtDNA) COI region and also examining interspecies genetic diversity among the isolates using COI haplotypes.

Study areas and collection of samples
A total of 55 flukes were collected from cattle at three abattoirs located in the KZN and MP provinces (Figure 1), respectively. The abattoirs served as catchment areas for cattle slaughtered in their respective areas and the exact location of origin of each animal was not determined. Twenty-one flukes were collected from cattle originating from the Swartberg and Lions River in the Underberg area of KZN (one fluke was collected from each animal), 17 were collected from the Barberton abattoir, Belfast, MP (two to three flukes collected per animal) and 17 from the Mpumalanga highveld, MP (one fluke from each animal) ( Table 1).
The Underberg region of KZN has a warm and cold season, with an average temperature of 14.1 °C. The average annual rainfal is 985 mm, with the altitude being 1540 m above sea level (Climate-Data 2017a). The climate and weather patterns in the Mpumalanga province depend on the topography of the specific region analysed. In Belfast, which is in the lowveld, the climate is warm and temperate and the average temperature is 13.2 °C and the average annual precipitation is 835 mm. In the highveld, the average annual temperature is 15.5 °C, with rainfall averaging 683 mm. The altitude of the province ranges from 1500 m to more than 2000 m above sea level. A decrease in temperature is noted from east to west, with the increase in altitude (Climate-Data 2017b).

DNA extraction, amplification and sequencing
The genomic DNA was isolated from tissue from the posterior end of the fluke using the phenol-chloroform method, which uses sodium dodecylsulfate (SDS) and proteinase K to enzymatically digest proteins and non-nucleic acid cellular components. A mixture of phenol:chloroform:isoamyl alcohol (25:24:1) was used to separate lipids and cellular debris into the organic phase, resulting in isolated DNA being in the aqueous phase (McKiernan & Danielson 2017). All DNA was stored at -20 °C until further analysis. Polymerase chain reaction (PCR) amplification of the partial COI region was done using the primers FHCO1 (forward: 5′-TTGGTTTTTTGGGCATCCT-3′) and FHCO1 (reverse: 5′ -AGGCCACCACCAAATAAAAGA3′) (Mucheka et al. 2015).
All PCR reactions were performed in 25 μL volumes, with each reaction containing 12.5 μL of DreamTaq PCR Master Mix (2X), 0.5 μL per primer of the forward and reverse primers, 10.5 μL of nuclease free water and 1 μL of template DNA. PCR was performed using a thermocycler (Bio-Rad, Hercules, CA, United States) under the following cycling conditions: initial denaturation at 95 °C for 3 minutes (min), followed by 48 cycles of denaturation at 95 °C for 30 seconds (s), annealing at 59 °C for 30 s and 72 °C for 1 min and a final elongation step of 10 min at 72 °C, with a 10 °C hold.
Polymerase chain reaction products were run on a 1% (w/v) agarose gel, which was stained with ethidium bromide. On each gel, a molecular weight marker of 100-1000 base pairs (bp) (O'GeneRuler, Fermentas, South Africa) was used to determine the size of visible bands. The gel was viewed using the ChemiDoc™ MP Imaging System (Bio-Rad) to confirm amplification.
The unpurified PCR products were subsequently sent to the Central Analytical Facilities (CAF), Stellenbosch University  (South Africa) for sequencing in the forward direction using next-generation Sanger sequencing technology. Using this technology, PCR-amplified DNA was denatured to single strands and chain-terminating dideoxynucleotides (ddNTPs) were added to the reaction. During DNA replication, the chain-terminating ddNTPs were selectively hybridized to the template strand by DNA polymerase. The chain-terminating signals were then used to determine the identity of the bases (Jamuar, D'Gama & Walsh 2016).

Sequence analysis
BioEdit version 7.2.5 (Hall 2013) was used to edit all sequences by calling up the base peaks in the chromatogram and adjusting N, S and K substitutions to appropriate the nucleotide bases. Sequence identity was then confirmed (  (Hammer, Harper & Ryan 2001) and principal coordinate analysis (PCoA) were used for multivariate ordination analysis. The PCoA plot was generated using the Gower similarity index and c = 2 transformation exponent. Phylogenetic relationships among the populations were inferred using MrBayes version 3.2.6 (Ronquist & Huelsenbeck 2003). Four Markov chains were run for 200 000 generations to ensure a standard deviation of split frequencies less than 0.01. The print frequency was set to 1000, with a sample frequency of 10 and a burn-in period of 20 000 for Bayesian inference (BI). The haplotype network was constructed using Network 5 (Bandelt, Forster & Röhl 1999) and the median joining rooting method under default (10) weight and with epsilon (ε) set to 0.

Ethical considerations
This study was approved by the Animal Research Ethics Committee of the University of KwaZulu-Natal (AREC/044/016PD).

Fluke identification
The trimmed COI sequences were BLAST searched against the NCBI GenBank database to obtain the closest matches, with percentage identities (  Figure 2). Two animals were co-infected with both F. hepatica and F. gigantica (Table 1).

Sequence variation and diversity
The aligned KZN sequences showed a sequence conservation of 0.57, indicating moderate sequence variation although conservation was high ( Table 2). The nucleotide diversity was estimated to be 0.086, indicating a low probability of difference in random sequences. Tajima's test was performed to examine the deviation of the sequences from mutation drift equilibrium. The test yielded a Tajima's D statistic of -0.75339, which is indicative of negative selection; however, the results were not statistically significant ( p > 0.10). The Belfast (MP) sequences showed a sequence conservation of 0.641, indicating high sequence conservation ( Table 2). The nucleotide diversity was estimated to be 0.069, also indicating low probability of random sequence difference. Tajima's test showed a D value of -0.84641; however, it was not statistically significant ( p > 0.10). Highveld (MP) sequences showed a sequence conservation of 0.457, indicating a higher sequence variation than conservation. The nucleotide diversity was estimated to be 0.109. A Tajima's D of -0.608 is also indicative of negative selection although not significant ( p > 0.10) ( Table 2).
The total aligned sequences showed a sequence conservation of 0.472, indicating moderate sequence variation with higher variable sites than conserved sites. The nucleotide diversity was estimated to be 0.082, showing low probability of random sequence difference. Tajima's D was estimated to be -1.100 and not statistically significant ( p > 0.05) ( Table 2). The lack of significance in Tajima's D could be attributed to the small sample size or could indicate the neutrality of the population.

Genetic distances between haplotypes
The genetic p-distances between the F. hepatica isolates ranged from 0.01 to 0.28, while a range of 0.02-0.04 was observed for the F. gigantica isolates. Interestingly, the Highveld haplotypes (H_15, H_16, H_17, H_20) showed a relatively higher genetic distance to the rest of the F. hepatica haplotypes, with a range from 0.09 to 0.28 (Table 3). Genetic distance between haplotypes identified within the two study species ranged from 0.07 to 0.35 (Table 3), indicating a close relation between isolates of the same species and a relatively distant interspecies genetic relation.

Phylogenetic analysis
The phylogenetic analysis based on the COI gene showed two main clades as shown in Figure 3. The first clade constitutes F. hepatica based on the identity of sequences (      Principal coordinate analysis supports the genetic distance and distribution of haplotypes, with the clustering of individuals being similar to the respective F. gigantica and F. hepatica haplogroups, as well as the outgroup. The PCoA plot also shows 11.38% of the variation being attributed to the left axis and 55.52% to the bottom axis ( Figure 5).

Discussion
Phenotypic identification of Fasciola isolates is routinely used although with difficulties (Ai et al. 2011). Molecular techniques in the differentiation of Fasciola isolates have been found to be effective in disease control and management (Ai et al. 2011;Beesley et al. 2017;Mucheka et al. 2015). In our study, we were able to differentiate F. hepatica and F. gigantica isolates using phylogenetic analysis of the partial COI region, which showed distinct clades of the two species. Fasciola hepatica isolates were restricted to KZN and the highveld area of the MP province and 76.5% to the lowveld area in Belfast. The F. gigantica isolates were restricted to the lowveld (Belfast) region of MP, with a 23.5% prevalence over F. hepatica in the same locality. This study has shown similar results, with respect to the distribution of F. hepatica and F. gigantica in the KZN and MP provinces of South Africa, to those recorded by Mucheka et al. (2015). The consistency in the results of this study with those observed in cited studies suggests the validity of the distribution patterns of Fasciola spp. in the study areas of South Africa.
A major driving factor of the wide distribution of F. gigantica and F. hepatica is the dependence of the life cycle on the fresh    (Ai et al. 2011). No intermediates were identified in the current study; however, the F. gigantica isolates from Belfast showed a close relationship with intermediates from China.
A link between environmental conditions and genetic diversity has been established based on Darwin's theory of adaptation to environmental stresses as well as geographical location (Walker et al. 2011). The results of this study demonstrate low interspecific genetic diversity (π range between metapopulations = 0.06901-0.10911 and 0.07784 overall) between the flukes. This was comparatively higher than that observed by Ai et al. (2011) in a study of diversity of flukes; however, their study consisted of isolates with a wider geographical distribution. Elliott et al. (2014) (Appleton 2003). Fasciola gigantica has been reported to be not as widely distributed as F. hepatica in South Africa, however (De Kock et al. 1989;Mucheka et al. 2015) unlike countries in Africa, such as Zimbabwe where it is predominant (Mucheka et al. 2015). Fasciola gigantica isolates from this study did not share any haplotypes with the F. gigantica GenBank isolates; however, they showed a relationship with only one mutational step between the haplotypes, which suggests a distant relation between the isolates. The isolates forming part of the F. gigantica haplogroup could share a distant maternal lineage with geographical isolation as well as microevolution contributing to diversity. Fasciola gigantica is reported to have spread to Africa earlier than F. hepatica (Mas-Coma et al. 2009) which could influence the rate of divergence of the populations isolated in South Africa because of the reported short time of haplotype change in fasciolids (Walker et al. 2011). The high haplotype diversity is also indicative of a recent population expansion, which can be attributed to the ability of fasciolids to adapt to a wide range of environmental conditions (Mas-Coma 2004;Mas-Coma et al. 2005). Climate change and human activities have been cited as major contributing factors of the transmission and spread of fascioliasis (Afshan et al. 2014;Mas-Coma et al. 2009). Life cycles of both the liver flukes and their intermediate host snails are highly dependent on geo-climatic parameters, with emphasis on temperature, rainfall and altitude (Mas-Coma 2004). Variations in climatic conditions can result in the variations in genetic diversity of the liver fluke between populations, displaying spatial and temporal separation (Walker et al. 2011). Studies have shown that husbandry and management of farms influence the movement of the definitive hosts, which, in turn, affects the population genetic diversity of the parasites (Beesley et al. 2017;Mas-Coma et al. 2005;Semyenova et al. 2006). Currently, the authors could not find conclusive studies on patterns of movement of cattle and other livestock in South Africa.

Conclusion
Although an increase in the importance of facioliasis as a zoonotic disease has been observed, the paucity of information in Southern Africa on the phylogenetic relationships as well as genetic diversity of Fasciola spp. is evident. This study has demonstrated the presence of F. hepatica in both the KZN and MP provinces of South Africa, confirming a recent observation by Mucheka et al. (2015) with F. gigantica only observed in the lowveld of the MP province. Fasciola hepatica was shown to be prevalent in the two provinces studied and shows similarity to studies by Mucheka et al. (2015). A link between the distribution of Fasciola spp. and that of the Lymnaeid intermediate host snails has been established in the literature (Mas-Coma et al. 2009;Mucheka et al. 2015), giving an inclination to attributing the fasciolid distribution to that of the snail intermediate hosts. The mitochondrial COI marker showed high genetic diversity between the study populations; however, further studies using more variable and neutral markers, such as microsatellites, are required to elucidate population genetic structuring. Future research should be focused on elucidating the phylogeography of Fasciola spp. in other provinces of South Africa, as well as the parasiteintermediate host snail interactions. This study has added to the current knowledge regarding the distribution of Fasciola spp. in the KZN and MP provinces of South Africa. The study was however limited by sample size and the geographical distribution of sampling sites, which could both be increased to further elucidate the genetic structuring and distribution of Fasciola on a wider scope.