Population genetic structure of Schistosoma haematobium and Schistosoma haematobium × Schistosoma bovis hybrids among school-aged children in Côte d’Ivoire

While population genetics of Schistosoma haematobium have been investigated in West Africa, only scant data are available from Côte d’Ivoire. The purpose of this study was to analyze both genetic variability and genetic structure among S. haematobium populations and to quantify the frequency of S. haematobium × S. bovis hybrids in school-aged children in different parts of Côte d’Ivoire. Urine samples were subjected to a filtration method and examined microscopically for Schistosoma eggs in four sites in the western and southern parts of Côte d’Ivoire. A total of 2692 miracidia were collected individually and stored on Whatman® FTA cards. Of these, 2561 miracidia were successfully genotyped for species and hybrid identification using rapid diagnostic multiplex mitochondrial cox1 PCR and PCR Restriction Fragment Length Polymorphism (PCR-RFLP) analysis of the nuclear ITS2 region. From 2164 miracidia, 1966 (90.9%) were successfully genotyped using at least 10 nuclear microsatellite loci to investigate genetic diversity and population structure. Significant differences were found between sites in all genetic diversity indices and genotypic differentiation was observed between the site in the West and the three sites in the East. Analysis at the infrapopulation level revealed clustering of parasite genotypes within individual children, particularly in Duekoué (West) and Sikensi (East). Of the six possible cox1-ITS2 genetic profiles obtained from miracidia, S. bovis cox1 × S. haematobium ITS2 (42.0%) was the most commonly observed in the populations. We identified only 15 miracidia (0.7%) with an S. bovis cox1 × S. bovis ITS2 genotype. Our study provides new insights into the population genetics of S. haematobium and S. haematobium × S. bovis hybrids in humans in Côte d’Ivoire and we advocate for researching hybrid schistosomes in animals such as rodents and cattle in Côte d’Ivoire.


Introduction
Schistosomiasis is a chronic neglected tropical disease caused by trematodes belonging to the genus Schistosoma [14,40]. The disease affects both humans and animals and is of considerable public health and veterinary concern, particularly in tropical and subtropical zones. According to the World Health Organization (WHO), schistosomiasis is transmitted in over 78 countries and territories throughout a wide belt of the tropics and subtropics [68]. More than 250 million people are infected, mostly in Africa [31], and the global burden of schistosomiasis was estimated at 1.4 million disability-adjusted life years (DALYs) in 2017 [26]. Six species of schistosomes can infect humans: S. guineensis, S. haematobium, S. intercalatum and S. mansoni in Africa and the Arabian Peninsula [60], and S. mekongi and S. japonicum in Asia; S. mansoni mainly occurs in Africa, but is also found in Brazil and some Caribbean islands [51]. Schistosoma haematobium causes the urogenital form of the disease. Classified as a group I carcinogen, urogenital schistosomiasis can lead to squamous-cell carcinoma of the bladder [35]. In most countries, while human schistosomiasis is well documented, little is known about the prevalence and transmission dynamics of animal schistosomiasis.
In Côte d'Ivoire, both S. haematobium and S. mansoni are endemic, causing urogenital and intestinal schistosomiasis, respectively [13]. The former is predominant in the central and southern parts [15,56], while the latter is mainly found in the western parts of Côte d'Ivoire [7,47,59]. We recently reported high prevalence of S. haematobium in school-aged children in South Côte d'Ivoire [5]. Limited data are available on animal-infecting schistosomes, such as S. bovis, a parasite of livestock and rodents. A previous study on post-mortem examinations of cattle in North Côte d'Ivoire reported a prevalence of 35% [1]. These findings were confirmed in a recent cross-sectional survey conducted in different parts along a transect from North to South Côte d'Ivoire, with the highest prevalence of S. bovis found in cattle in the northern parts of the country [39].
Habitat change and migration of hosts, which is also linked to climate change, can influence the epidemiology and distribution of schistosomiasis and enhance the occurrence of interspecies hybridization between human and animal schistosomes [38]. Schistosoma haematobium Â S. bovis hybrids are well documented in West Africa, particularly in Benin, Mali, Niger and Senegal [34,42], with occurrences also reported in Malawi [66] and Corsica, France [9,10]. More generally, genomic studies have highlighted the fact that several natural strains that were initially identified as "pure" S. haematobium are in fact introgressed to some extent with S. bovis [44,48]. Recently, hybrids from miracidia (although not thought to be viable) have been identified in Côte d'Ivoire between S. haematobium and S. mansoni excreted by humans [17,41], and between S. haematobium and S. bovis being transmitted by freshwater snails of the genus Bulinus [58] and also from humans [6].
Population genetic structure and genetic diversity varies by Schistosoma species [49]. Among the African species, S. mansoni is the most widely studied. Several investigations have shown high genetic diversity and strong genetic structure in several countries, including Ethiopia, Kenya, Senegal and Uganda [2-4, 57, 61, 63]. In contrast, studies on S. haematobium have shown both less genetic diversity and less population structure compared to S. mansoni [30,46]. This difference in genetic structuration between S. mansoni and S. haematobium is also apparent at the continental scale [30,62,63]. Population genetics data for S. bovis are even more limited. One study has investigated the population genetic patterns of S. bovis in Cameroon revealing an intermediate pattern with high genetic diversity (i.e. like in S. mansoni), but no genetic structuration (i.e. like in S. haematobium) at the country scale [18]. At the continental scale, a comparative genomic study revealed a significantly higher genomic diversity between S. bovis lineages than between S. haematobium lineages [48]. More recently, a study analyzed the genetic patterns of parasites collected from humans and animals in North Senegal and molecularly characterized both microsatellite markers and hybrid genotyping [11]. The authors demonstrated (i) clear genetic separation between parasites recovered from animals compared to those from humans; (ii) no genetic structuration between hybrids and pure parasites from human hosts; and (iii) significant genetic differentiation between different villages of the Senegal River basin. These authors concluded that, in Senegal, animals do not represent a real reservoir for human schistosomiasis and also, there is no cross-over of transmission between humans and animals. However, this conclusion contradicts recent observations from Benin, where S. haematobium Â S. bovis hybrids were found in cows [54] and rodents [55]. The purpose of this study was to investigate the population genetic structure and genetic diversity of S. haematobium and S. haematobium Â S. bovis hybrids that were obtained from human urine samples, in four sampling sites of Côte d'Ivoire. We molecularly characterized individual miracidia, using cox1 and ITS2 markers to identify S. haematobium and S. haematobium Â S. bovis hybrids within the different populations, alongside microsatellite markers [65] to analyze population genetic diversity and structure.

Ethics statement
Ethical clearance was obtained from the Ministère de la Sant e et de l'Hygiène Publique en Côte d'Ivoire (reference no. 003-18/MSHP/CNER-kp). School authorities, teachers, parents/guardians and children were informed about the objectives, procedures and potential risks and benefits of the study. Written informed consent was obtained from children's parents or legal guardians, while children provided oral assent. After sampling, a praziquantel treatment (40 mg/kg) was offered to children found with a Schistosoma infection.

DNA extraction for miracidia genotyping
Genomic DNA from individual miracidia was extracted using Chelex Ò beads (Bio-Rad; Hercules, CA, USA) from 2.0 mm discs containing sample that had been punched out of the FTA card using a Harris-Micro-Punch (VWR; London, UK [37]. Cox1 and ITS analysis for species and hybrid identification A cox1 multiplex PCR was performed to identify the species-specific mitotype of each sample, as described in our previous work [6]. For the analysis of the nuclear internal transcribed spacer 2 (ITS2) region, we used a restriction fragment length polymorphism (RFLP) approach ( [19] and Supplementary Fig. S1), following amplification using the forward primer Sc_ITS_F: 5 0 -GGC TGC AGC GTT AAC CAT TA-3 0 and reverse primer Sc_ITS_R: 5 0 -ACA CAC ACC ATC GGT AC AAA-3 0 , which targets 505 bp of the ITS2 [19]. We performed PCRs in a total reaction volume of 25 lL, comprising 2 lL of gDNA, 5 lL of Green GoTaq flexi buffer 5Â, 1.5 lL of 25 mM MgCl 2 , 0.5 lL of 10 mM dNTP mix, 1 lL of each 10 lM primer and 1 U of GoTaq Hot Start Polymerase (Promega; Madison, WI, USA). The reaction conditions included an activation step of 95°C for 3 min, followed by 45 cycles of 95°C for 40 s, 58°C for 40 s and 72°C for 40 s, and a final extension at 72°C for 6 min. In the subsequent step, 5 lL of PCR products were then digested at 37°C using 0.5 lL of restriction enzyme MboI (Thermo Fisher Scientific; Waltham, MA, USA), 2.5 lL of CutSmart buffer, and 17 lL of molecular water for 15 min, followed by inactivation at 60°C for 20 min for a total reaction volume of 25 lL. The cut sites of the MboI enzyme for species specific SNPs were based on the following sequences positions: ;GATC and CTAG". The PCR products were not cleaned or visualized before enzyme digestion. We used only samples that were positive with the cox1 multiplex PCR method. This enzymatic digestion cuts the 505 bp ITS2 fragment into four bands for S. haematobium (281 bp, 82, bp, 98 bp and 44 bp) and three bands for S. bovis (379 bp, 82 bp and 44 bp). All digested ITS2 PCR products were visualized using 2% agarose electrophoresis gels stained with GelRed TM (Biotium Inc.; Darmstadt, Germany). The nuclear ITS2 and mitochondrial cox1 identities for each individual miracidia were combined to produce a mito-nuclear genetic profile of S. haematobium (Shcox1 + ShITS2: Sh Â ShSh), S. bovis (Shcox1 + ShITS2: Sb Â SbSb) or S. haematobium-bovis hybrid genotypes, the latter of which is recognized by mito-nuclear discordance (Shcox1 + SbITS2: Sh Â SbSb), (Shcox1 + Sb/ShITS2: Sb Â SbSh), (Sbcox1 + ShITS2: Sb Â ShSh) or (Sbcox1 + Sh/SbITS2: Sb Â ShSb)).

Microsatellite analysis
Individual miracidia were further genotyped using a set of 18 microsatellite markers divided into two panels developed by Webster and colleagues [65]. The forward primers were fluorescently labelled with 6-FAM, VIC, NED and PET dyes (Applied Biosystems; Foster City, CA, USA) to enable identification within the multiplex PCR, as previously described [65]. Microsatellite PCRs were performed, using a Microsatellite PCR Kit (Qiagen; Hilden, Germany), in a final volume of 10 lL, including 4 lL of the DNA template, 5 lL of the 2Â microsatellite PCR Buffer Kit (Qiagen; Hilden, Germany) and 1 lL of 10Â microsatellite primer mix in two PCR panels of nine per multiplex. Thermal cycling was performed with an initial hot-start activation of 15 min at 95°C, followed by 45 cycles of 94°C for 30 s, 56°C for 90 s and 72°C for 60 s, with a final extension at 60°C for 30 min. The PCR products were sent to Genoscreen (Lille, France) for genotyping. All microsatellite loci were visually peak called, using GS500Liz size standard (Applied Biosystems), and GeneMarker software. Only 16 loci were subsequently used for analysis as allelic dropout of two loci (C131 and Sh8) occurred in at least 20% of the samples (Supplementary Tab. S1). Errors due to large allele dropout or stutter bands and evidence of the presence of null alleles at each locus were checked using Micro-checker, version 2.2.3 [43].

Genetic diversity
The number and percentages of the different cox1/ITS2 profiles obtained from each miracidium were calculated for each site and any differences using a v 2 test was observed. We also calculated the number of S. bovis and S. haematobium haplotypes for cox1 or alleles for ITS2 among sites. The mitochondrial cox1 haplotype is unique as it is inherited from the maternal line, while the ITS2 was scored as one for heterozygous (Sb/Sh) and two for homozygous (SbSb or ShSh). Differences in the relative frequencies of S. bovis vs. S. haematobium haplotypes or alleles between sites were tested using a binomial test [32].
For the microsatellite data, tests for deviation from Hardy-Weinberg equilibrium per locus and site were carried out using Genepop, version 4.0 [53]. The genotypic disequilibrium test for pairs of loci overall and the adjusted p-value for 5% nominal level was performed using FSTAT, version 2.9.4 [28]. Genetic diversity was compared between sites and between cox1/ITS2 genetic profiles. Expected heterozygosity (He), number of alleles (A), allelic richness (Ar), and the inbreeding coefficient (F IS ) of each microsatellite locus were computed per sampling site, using FSTAT version 2.9.4 [28]. He and Ar were compared between the populations using the pairwise Friedman rank test, followed by the Nemeyi multiple comparison test implemented in the R studio PCMP plus package [36].

Genotypic differentiation and population structure
Genotypic differentiation between sampling sites was assessed using pairwise F ST values [63], calculated in FSTAT, version 2.9.4 [28,67] with a threshold of significance adjusted for multiple tests using Bonferroni's standard correction [50]. Pairwise F ST values were also calculated between samples that exhibited the different cox1/ITS2 profiles. Principal component analysis (PCA) was performed to compare either the four sites or the six cox1/ITS2 profiles, using GENETIX, version 4.05 [8].
The uppermost level of genetic structure for all individuals was determined by a Bayesian clustering approach using Markov Chain Monte Carlo (MCMC) permutation analysis, implemented in STRUCTURE, version 2.3 [45], using data of children from whom we collected at least 15 miracidia. The length of the burn-in period was 250,000 with the number of MCMC replicates after burn-in at 1,000,000 and K from 1 to 8 using an admixture model. The log likelihood for each K was averaged over three runs with the CorrSieve package in R (Vienna, Austria) and the delta K-values were then computed to determine the most likely number among the K values tested [21]. For the most likely number of genetic clusters, an additional 10 runs were computed with the same initial parameters as those described. The probability of each miracidium belonging to each cluster was averaged over the 10 runs and graphically represented using CLUMPP, version 1.1.2 [23] and DISTRUCT, version 1.1 [52].
Analysis of molecular variance (AMOVA) was performed to evaluate the partitioning of the overall genetic variance to the hierarchical levels at infrapopulation "within miracidia", "among miracidia within individual children", "between miracidia from individual children within each site" and "between sites" using ARLEQUIN, version 3.5 [22]. The level "within miracidia" reflects the differences between microsatellite loci and has no population genetic meaning. Hence, we report the residual variance of the remaining three levels and normalize their relative contributions to 100%.
Population genetic structure between parasite infrapopulation was assessed independently for each sampling site using STRUCTURE, version 2.3 with the same procedure as described above. If no best K was found, K was set as the number of children at each sampling site (16 for Adzop e, 18 for Agboville, 18 for Duekou e and 22 for Sikensi). We first calculated the mean probabilities assignations of the miracidia from a given child to a given cluster. Then repartitions of the miracidia among the different clusters were tested in each participant using v 2 under the null hypothesis of random repartition (1/K).

Genetic diversity
Of the 2164 miracidia from which cox1/ITS2 profiles were obtained, 1966 (90.8%) also produced reliable allele data for at least 10 microsatellite loci for subsequent analyses. We found a significant genotypic disequilibrium (adjusted p-value at 5% level < 0.001) for 120 pairwise locus combinations across the whole dataset. Genetic diversity indices (He, A, Ar and F IS ) and the probability of deviation from Hardy-Weinberg equilibrium (P HWE ) for each microsatellite locus by sampling site are shown in Table 2. Most loci showed a deviation from Hardy-Weinberg equilibrium. All loci were highly polymorphic and the number of alleles and the allelic richness ranged from 7 to 17 and 7.0 to 15.8, respectively per locus in the whole dataset. Inter-population comparisons showed that heterozygosity (He) was relatively stable (Friedman test, v 2 = 1.33, degree of freedom (df = 3, p = 0.72) and that allelic richness was more variable ranging between 8.8 and 10.5 for Sikensi and Adzop e, respectively. A significant difference was observed for this last parameter (Friedman test, v 2 = 15.33, df = 3, p = 0.002), allelic richness from Duekou e being statistically higher compared to the remaining three sites (Nemeyi test, p < 0.05). No differences in genetic diversity indices were observed according to the different cox1/ITS2 genetic profiles.

Population genotypic differentiation
The pairwise F ST values between sites ranged from 4.1% to 6.9% with Duekou e being consistently most differentiated from the other sites (Table 3). The F ST values between the samples that presented these cox1/ITS2 profiles ranged from 0.3% to 2.7% and thus showed that the S. bovis cox1-ITS2 profiles (Sb Â SbSb) were not different from the other cox1/ITS2 profiles ( Table 4). The miracidia with the S. bovis cox1 Â S. bovis ITS2 (Sb Â SbSb) genetic profile came from three sites (7 from Adzop e, 4 from Agboville, and 4 from Sikensi) (Supplementary Tab. S2).
The PCA revealed only weak structuration among the four sites with only Duekou e being partially separated from the other three sites (Fig. 2). Comparing the different genetic profile, the PCA showed no genetic structuring according to the six different genetic profiles (Fig. 3).
Population genetic structure was assessed both at the population level (between sampling sites) and at the infrapopulation level (between children within each sampling sitesee below). K = 4 displayed maximal delta K at the population level suggesting that all miracidia were grouped into four genetic clusters. The STRUCTURE analysis showed a relatively strong genetic structure between sampling sites (Fig. 4). At the population level, an average of 51% of the genotypes were assigned to the respective dominant cluster in Adzop e, 63% in Agboville, 69% in Sikensi, and 83% in Duekou e.

Sites
Sb Â SbSb Sh Â ShSh Sb Â ShSb Sb Â ShSh Sh Â SbSb Sh Â ShSb Total hybrids All total Cox1 haplotypes ITS2 alleles belong to each cluster by child, using only those children who hosted at least 15 miracidia. The analysis of molecular variance (AMOVA) showed that most of the miracidial genetic variance was found within sites (57.8% of the variance). There is more variance among miracidia by site rather than by host (10.4%) ( Table 5).

Discussion
Our study showed a high proportion of S. haematobium Â S. bovis hybrids among school-aged children in different parts of Côte d'Ivoire. Similar results were also reported from a large study in Senegal, where most children (88%) excreted hybrid miracidia [64]. Notably, virtually all children (97.8%) were producing miracidia with a hybrid genotype (i.e. a discrepancy between the mitochondrial and the nuclear signature and/or heterozygous nuclear profiles). Moreover, we found no genetic differentiation of miracidia presenting S. haematobium, S. bovis or hybrid genetic profiles, which indicates that all genotyped parasites most likely belong to a single genetic entity. This finding is in line with a recent genomic study indicating that most, if not all S. haematobium lineages across Africa, are likely to be introgressed to some extent with S. bovis [48].
We assigned all miracidia to six distinct cox1-ITS2 profiles, of which four were unambiguous hybrid profiles. Hybrids with a mitochondrial cox1 haplotype from S. bovis and homozygous ITS2 from S. haematobium (Sb Â ShSh) were the most Table 4. Pairwise estimates of F ST (below the diagonal) and significance (above diagonal) between cox1-ITS2 profiles of the four sampling sites based on 16 microsatellite loci.    commonly identified type at 41.9%. Such hybrids have sometimes been documented in Senegal and ascribed to bidirectional introgressive hybridization [33,63]. Beside this general pattern, significant variation was observed between sampling sites with the highest frequency of hybrids (reaching 79%) identified in Adzop e. Interestingly, the genetic contribution of the S. bovis genetic profile to the genetic make-up of schistosome populations is much more evident at the mitochondrial level than at the nuclear level. The frequency of miracidia with a heterozygous ITS2, i.e. carrying both S. bovis and S. haematobium ITS alleles, was relatively low, which could be explained by the concerted evolution of ITS2 gene region [25]. The S. bovis cox1 Â S. bovis ITS2 genetic profile (SbSb) has also been reported in miracidia recovered from humans in  Corsica [10]. We only detect 0.7% of this profile in our study. Importantly, this genetic profile using only two gene genotyping methods does not mean that these parasites are "pure" S. bovis parasites at the whole genome level. The genetic signature can only represent a certain introgression level between both species, and can explain some discrepancy in parasite phenotype such us the route of excretion expected to be in the feces for the S. bovis parasite.
Data on the genetic structure of S. haematobium populations are sparse [29,30]. Previous findings based on a phylogeographic approach, using cox1 data, at the African continental scale suggested that S. haematobium is less structured than S. mansoni [63] and S. bovis [48]. In contrast, the significant differentiation between sampling sites and even within individual hosts reported here, suggests a genetic structure of miracidia at small geographical scales. The large geographical distance separating Duekou e in the western part of Côte d'Ivoire to the other three sites is likely to limit mixing between these populations and lead to genetic structuring between the schistosome populations. More surprising is the genetic structure observed at the individual level. In the majority of cases, each child hosted an infrapopulation different from the other children. Our finding does not corroborate the "genetic mixing bowl hypothesis" [16]. The authors reported that the long-living definitive hosts cumulate genotypes coming from several short-living infected intermediate host snails, which can lead to homogeneity of infrapopulations. The pattern found in our study could be explained by recent treatment of children by praziquantel one year before sample collection. In this scenario, children have only recently started to harbour parasites again and there was insufficient time for "mixing" after praziquantel treatment. This would homogenize parasite populations among children rather than promoting structuration. Children would therefore be infected with only one or a few pairs of parasites and they could also be highly exposed and have several infections. This clustering may be amplified by the genotypic specific immune responses. In a mouse model, it was found that a trickle infection (repeated light infections) can protect mice against a Schistosoma challenge infection [20]. Subsequently, the same authors have shown that this level of protection is genotype-dependent and correlated with genetic dissimilarity between the immunizing and the challenging infection [20]. The smaller the genetic distance between immunizing and challenging clones, the lower the infectivity rate of the challenging clones. This mechanism may, in turn, increase genetic diversity in each host and differentiate the parasite infrapopulation between hosts. No matter what the reasons are for this pattern, it emphasises the importance of investigating the parasites from as many different patients (or hosts in general), when investigating population genetics between different sampling sites. Increasing the number of hosts sampled, rather than the number of miracidia per host, has already been proposed after a stochastic re-sampling approach [24].
No genetic differentiation was observed according to the genetic profile. This result suggests that there are no genetic mating restrictions between the six different genetic profiles, and hence, all recombinations are possible [11,27]. Similar findings have recently been reported from Senegal [11]. So, our results together with recent genomic studies [44,48], call into question S. haematobium and whether forms of hybrids detected in humans thus far are truly different genomic entities.
The current study has some limitations, related to the reliability of the RFLP analysis, the absence of samples from animals, or the hybrid genotyping method. First, in our study we used an RFLP method based on a short 505 bp ITS2 fragment. Even though the results we obtained are robust, the RFLP method is known to be sensitive to different factors such as enzyme used, and the time or temperature of digestion. A new method for single-nucleotide polymorphism detection adapted to large scale screening needs to be developed. Second, miracidia from livestock or rodents were not collected. This is an important point to infer the possible zoonotic capacity of these hybrid pathogens. The situations are contrasted depending on the countries without current evidence of zoonotic transmission in Senegal [12] compared to rodent and cow reservoir hosts in Benin [54]. Finding a high frequency of hybrids in Côte d'Ivoire does not confirm zoonotic transmission and we advocate looking for the presence of hybrid schistosomes in rodents and livestock. Third, mito-nuclear genotyping methods are well adapted for large population screening; however, as previously mentioned, the whole genome cannot be resumed to these two gene assignations, and genomic studies using more genetic markers are desirable.

Conclusions
Our study investigated the genetic diversity and population structuring of schistosomes at the population level. Each child harboured a genetic cluster of schistosomes at the infrapopulation level, which could lead to putative resistant parasite selection due to drug pressure by mass drug administration (MDA). The high frequency of hybrids between S. haematobium Â S. bovis observed can lead to difficulties in accurately diagnosing schistosome hybrids by conventional techniques using light microscopy. This could also have adverse consequences on schistosomiasis control towards elimination and negatively impact disease transmission. More studies are need on population genetics of schistosomes at the human and animal interface to evaluate the parasite's gene flow.

Contributions of authors
EKA, HM, OB and JB conceived and designed the study. EKA, AV and JFA performed the molecular analyses. EKA, and JB performed statistical analysis. EKA wrote the first draft of the manuscript. JTC, WY, AOT, EKN, JZ, JU, OB and JB revised the manuscript. MSW helped with figure resolution. All authors read and approved the final manuscript before submission and resubmission.