Population Genetic Structure and Hybridization of Schistosoma haematobium in Nigeria

Background: Schistosomiasis is a major poverty-related disease caused by dioecious parasitic flatworms of the genus Schistosoma with a health impact on both humans and animals. Hybrids of human urogenital schistosome and bovine intestinal schistosome have been reported in humans in several of Nigeria’s neighboring West African countries. No empirical studies have been carried out on the genomic diversity of Schistosoma haematobium in Nigeria. Here, we present novel data on the presence and prevalence of hybrids and the population genetic structure of S. haematobium. Methods: 165 Schistosoma-positive urine samples were obtained from 12 sampling sites in Nigeria. Schistosoma haematobium eggs from each sample were hatched and each individual miracidium was picked and preserved in Whatman® FTA cards for genomic analysis. Approximately 1364 parasites were molecularly characterized by rapid diagnostic multiplex polymerase chain reaction (RD-PCR) for mitochondrial DNA gene (Cox1 mtDNA) and a subset of 1136 miracidia were genotyped using a panel of 18 microsatellite markers. Results: No significant difference was observed in the population genetic diversity (p > 0.05), though a significant difference was observed in the allelic richness of the sites except sites 7, 8, and 9 (p < 0.05). Moreover, we observed two clusters of populations: west (populations 1–4) and east (populations 7–12). Of the 1364 miracidia genotyped, 1212 (89%) showed an S. bovis Cox1 profile and 152 (11%) showed an S. haematobium cox1 profile. All parasites showed an S. bovis Cox1 profile except for some at sites 3 and 4. Schistosoma miracidia full genotyping showed 59.3% of the S. bovis ITS2 allele. Conclusions: This study provides novel insight into hybridization and population genetic structure of S. haematobium in Nigeria. Our findings suggest that S. haematobium x S. bovis hybrids are common in Nigeria. More genomic studies on both human- and animal-infecting parasites are needed to ascertain the role of animals in schistosome transmission.


Introduction
Schistosomiasis is one of the major neglected tropical diseases with public and veterinary health concerns and is endemic in tropical and subtropical regions. With a global burden of about 1.4 million disability-adjusted life years (DALYs), the disease is ranked second after malaria based on morbidity [1]. The Schistosoma genus shows a wide definitive host spectrum that ranges from humans to domestic and wild animals. Humans could be infected with one or more of the six human-infecting Schistosome species, and this may lead to combined disease symptoms and co-morbidities. Four human-infecting parasites (Schistosoma haematobium, S. mansoni, S. intercalatum, and S. guineensis) are common in Africa, while two (S. japonicum and S. mekongi) are prevalent in Asia [2]. Except for

Cox1 Phylogenetic Trees
On the 59 Cox1 sequences we have sequenced 21 different haplotypes: 2 S. haematobium and 19 S. bovis haplotypes. An S. bovis Cox1 phylogenetic tree performed with the haplotypes recovered from Nigeria and several other haplotypes from Cameroon, Benin, Senegal, Ivory Coast, Kenya, and Tanzania did not reveal any spatial structuration (see Supplementary Figure S2).

Microsatellite Analysis
No significant deviation from HW equilibrium or linkage disequilibrium was observed across loci. The genetic variability indices (He, A, Ar, and Fis) are shown in Table 3 from the 14 microsatellite loci. Mean heterozygosity across the population ranged from 0.527-0.598 but we did not observe any significant difference (p > 0.05). For allelic richness (Ar), the mean values ranged from 4.718-6.929. Significant differences were observed between site 9 and sites 2 to 4 (p < 0.05).

Population Genetic Structure
The pairwise genetic differentiation estimates (F ST ) between most of the sampling sites are statistically significant after Bonferroni correction, except for between sites 1 and 2, 8 and 9, 8 and 10, 9 and 10, 7 and 11, and 7 and 12 ( Table 4). The PCA ( Figure 1) showed a structuration among populations with the population from the west (1-4) separated from the population in the east (sites 7-12). Random sampling of two miracidia by patient does not change the latter result (Supplementary Figure S3) This genetic structure was confirmed using Structure software that revealed the highest probability for two clusters (K = 2) ( Figure 2). Table 4. Pairwise genetic differentiation estimate (F ST -above the diagonal) and the Euclidian geographic distances (Km-below the diagonal) between the sampling sites. Most F ST values are statistically significant (marked with an asterisk (*)) with the level of significance adjusted with Bonferroni correction (p < 0.0011). No link was observed between the geographical and genetic distances (Mantel test; p > 0.05).

Discussion
We report for the first time the population genetic structure and hybridization of S. haematobium in Nigeria. Based on the Cox1 profile, our study revealed a country-wide minimum proportion of 89% prevalence of S. bovis x S. haematobium hybrids and almost equal repartition among the study sites. Most studied sites revealed a hybrid prevalence of 100% except for sites 3 and 4. S. bovis x S. haematobium hybrid prevalence obtained from other West African countries: Cote d'Ivoire 57.5% [9] and Senegal 9-72% [7,8,12,[28][29][30] revealed that Nigeria has the highest prevalence of hybrids. An important variation in hybrid frequency, ranging from 2% to 26%, between different villages has been evidenced in Senegal [29]. The authors have positively associated this variation with the prevalence of S. mansoni. They hypothesized that a first schistosome infection would favor ongoing infections and subsequently hybridizations. Because hybrid prevalence is 100% in the majority of the sites we have sampled, we cannot test for an eventual link with proximal factors such as prevalence and socio-demographic factors we have measured [31].
We have obtained full genotypes (both Cox1 and ITS2) for a sub-sample of 59 parasites. Interestingly, we found a high percentage (39%) of Sb x SbSb genotypes. This genotype has not been found in Senegal [8,12,28,29] and is in a very low percentage in Cote d'Ivoire [32] and in Corsica [18]. The high percentage of Sb x SbSb genotype we found is associated with a preponderance of both S. bovis Cox1 haplotype and S. bovis ITS2 alleles compared to S. haematobium. Concerning the ITS2, in Cote d'Ivoire and in Senegal, the frequency of the S. haematobium allele is 87% and more than 88%, respectively [8,12,28,29]. We found only 40% of S. haematobium ITS2 alleles in Nigeria. When a population is at equilibrium, the ITS is expected to harbor a single allele from one of the parents resulting in a concerted evolution [33]. This supposes that in Nigeria, contrary to other countries the population of hybrid schistosomes is not stabilized.
Concerning Cox 1, in Cote d'Ivoire and Senegal, the frequency of the S. haematobium haplotype is 46% and more than 77%, respectively [8,12,28,29]. We found only 11% of the S. haematobium haplotype in Nigeria, and these haplotypes were restricted to two sites. As previously proposed by Boon et al. [29], two main factors could explain a variation in Cox1 haplotype frequency: genetic drift and/or selection. Because the mitochondrial genes are only maternally inherited, they are more prone to genetic drift compared to bi-parentally (i.e., nuclear) inherited markers. Mate choice or mate competition could select for a given mitochondria species. Recently, this has been shown in random mating between S. haematobium and S. bovis excluding the selection of mitochondria through sexual selection [34]. Boon et al. [29] proposed that the environment could select for different mitochondrial haplotypes. For instance, these authors hypothesized that the snail strain host could select for hybrid parasites in a given area. This interesting snail driver selection hypothesis needs to be tested.
Concerning Nigeria, the high frequency of S. bovis genes could be explained by active zoonotic transmission and ongoing gene flow between animal (i.e., S. bovis) and human parasites. Recent genomic studies have shown that the S. haematobium x S. bovis hybrid is certainly the result of an ancient introgression event [35,36]. The age of the hybridization does not exclude ongoing zoonotic transmission. This zoonotic transmission has been evidenced in Benin with cows and rodents [20,21] and only with rodents in Senegal [19]. S. haematobium x S. bovis hybrids have not been evidenced in cows in Cameroon [6]. Consid-ering the high prevalence of S. bovis genes in parasite-infecting humans in Nigeria, looking for the presence of hybrid schistosomes in animals (rodents or cows) seems necessary.
To determine the genetic structure among the populations, we measured the pairwise genetic estimates (F ST values), for all pairs. Generally, values < 0.05, 0.05-0.15, 0.15-0.25 and >0.25 indicate low, moderate, high and very high genetic differentiations respectively [6]. Our study revealed F ST values of 0.0104-0.1688 which is an indication of low to very high genetic differentiation among the populations. Few population genetic studies involved S. haematobium compared with S. mansoni [36], and the studies involved local scales, between 8 and 45 km distances between sites for Gower et al. [37] and Boon et al. [28], respectively. Our study proposed a wider range from local (10's of kilometers) to regional (10s to 100s of kilometers) scale. When populations are separated by a few kilometers, paired F ST values are in agreement with previous studies and range from 0.01 to 0.04 [37]. At the regional scale, the Fst values for S. haematobium are similar to S. mansoni [38,39].
Regardless of the method used (PCA or Bayesian analyses using Structure software), we showed a clear clustering into two groups of populations: one from the west (populations 1-4) and one from the east (populations 7-12). S. haematobium populations are usually not well structured compared to S. mansoni [36]. These two parasite species possess similar transmission dynamics that could influence the parasite's genetic structure. For instance, for both species, the transmission is focused on water bodies, the intermediate and definitive hosts have similar mobility, and the number of intermediate host species is restricted. It is well established that S. haematobium has less genetic diversity than S. bovis or S. mansoni [6,36]. This low genetic diversity reduces the power of determination of structuring units. Our study shows that no structuring units are detectable under around 250 km distances between populations. Fst values are also lower among populations 1-4 or among populations 5-10 than between the two clusters of populations. In comparison, it has been shown that clear population structures between S. mansoni populations are separated by a 127 km distance in Ethiopia [39].
Various factors including hybridization can favor genetic structure [27]. Introgression through hybridization can influence the genetic structure by adding new alleles in a given area and in turn favor population clustering. We have obtained 100% of hybrids in the majority of sites we have sampled. Nowadays, this does not exclude the influence of hybridization in genetic structuring. Indeed, the molecular barcoding method we have used only infer the presence/absence of hybrids and not the genomic introgression level. Molecular markers such as SNP are needed to infer the role of hybridization in genetic structuring.

Study Area and Study Population
This study was carried out in twelve sites in Nigeria, West Africa ( Figure 3). This study was integrated into a survey carried out on prevalence and risk factors associated with urinary schistosomiasis among primary school-age pupils in Nigeria (Onyekwere et al. submitted).

Urine Sample Collection and Miracidia Sampling
A labeled, clean, and sterile plastic container with an "identification code" for anonymity of 20 mL was given to each patient whose parents or legal guardians gave oral consent. Each participant whose urine sample was positive for the parasite was treated with a single oral dose of 40 mg/kg body weight of praziquantel (600 mg, Biltricide, Bayer, Leverkusen, Germany) through their Primary Health Center (PHC). Individual miracidium was harvested using a P10 Gilson micropipette in 3 µL of water under a 20× or 40× magnification binocular microscope. About 20-25 miracidia were individually captured for each participant with each miracidium being checked in the pipette tip before placing on Whatman FTA ® cards (GE Healthcare Life Sciences; Amersham, UK). Each FTA ® card filled with miracidia was stored at room temperature while on the field and transferred to "Laboratoire Interactions Hotes-Pathogenes-Environnements" (IHPE), France, for genetic analysis. Table 5 shows the number of miracidia collected from participants and genotyped with Cox1 and microsatellites for each of the sampling site.

DNA Extraction
Genomic DNA from Schistosoma randomly selected miracidia were individually extracted from FTA ® cards using the Chelex method [40]. Harris-Micro-Punch (VWR; London, UK) was used to perforate a 2 mm disc at the center where the sample was placed. The disc was washed in 50 µL ultra-pure water for 10 min, the water discarded, and the disc was incubated in 80 µL of 5% Chelex ® solution (Bio-Rad; Hercules, CA, USA) at 65 • C for 30 min with agitation. This was incubated again at 99 • C for 8 min without agitation. The solution was centrifuged at 14,000 rpm for 2 min and 60 µL of the supernatant was transferred into a 96-wells micro-plate and stored at −20 • C for genomic analysis.

Estimation of Hybrid Prevalence by Mitochondrial DNA Identification
Hybrid schistosomes are generally characterized by the combination of the maternal DNA (mt-DNA) from Cox1 and the nuclear DNA (rDNA) from ITS2 [7,9]. The results will be used to assign each parasite a genetic signature based on the haplotype-alleles combinations: Sb x SbSb, Sb x ShSh, Sb x SbSh, Sh x ShSh, Sh x SbSb or Sh x SbSh. We obtained this full genotyping only for a subsample of miracidia (see below). However, a basic estimation of the hybrid frequency could be assessed only by the Cox1 gene characterization considering that human infected by an animal parasite (i.e S. bovis) gene is a hybrid parasite [29]. Hence, the frequency of hybrids is a synonym for the frequency of miracidia with an S. bovis Cox1 profile. This method can lead to the underestimation of the frequency of hybrids because Sh x SbSb or Sh x SbSh genotypes are considered as pure S. haematobium instead of hybrids. For this purpose, each miracidium was molecularly characterized by rapid diagnostic multiplex PCR (RD-PCR) on Cox1 [41]. The number of miracidia tested per site is presented in Table 1. We used species-specific primers to amplify the region to discriminate each Schistosoma species fragment: S. haematobium 120 bp, S. mansoni 215 bp, and S. bovis 260 bp [41]. The primers we used were a single universal reverse primer; (Shmb.R, 5-CAA GTA TCA TGA AAY ART ATR TCT AA-3 ) and three species-specific forward primers; (Sh. , and 4 µL of DNA extract, making a total volume of 10 µL for the PCR mix. Thermal cycling was performed in (plate thermal cycler) a PerkinElmer 9600 Thermal Cycler (PerkinElmer, Waltham, MA, USA) and the PCR conditions used were: pre-denaturing at 95 • C for 3 min; 45 cycles of 10 s at 95 • C (denaturing), 30 s at 52 • C (annealing), and 10 s at 72 • C (extending). This was followed by a final extending period of 2 min at 72 • C. The PCR product was stored in the refrigerator at 4 • C until use. The PCR products (Cox1) were visualized on 2% agarose gel stained with 8 µL Midori dye. Nine microliters of the PCR product was loaded into each well using a multi-channel micro-pipette (including wells for positive controls; S. haematobium, S. mansoni, S. bovis, and water for a negative control) and 4 µL for size standard 100 bp (base-pair) ladder. The PCR products in the gel were analyzed by electrophoresis at 135 V for 30-35 min and transferred to the UV trans-illuminator where gel images were taken.

Mitochondrial DNA (Cox1) and Nuclear Internal Transcribed Spacer II (ITS2) Sequencing
Full genotyping on a sub-sample was assessed by SANGER sequencing of the two genes (Cox1 and ITS2). Cox1 and ITS2 genes of six to seven miracidia harboring S. bovis Cox1 RD-PCR profile were sequenced for all sites. Seven more miracidia were sequenced for sites 3 and 4, the only sites harboring the S. haematobium Cox1 RD-PCR profile (see results). S. haematobium Cox1 PCR mix was performed in 96 wells with a single forward primer COI1_F: 5 -GGGGGTTTTATTGGTTTAGGTT-3 and a single reverse primer COI1_R: 5 -CCAATTATAAAAGGCCATCACC-3 , while S. bovis COI1 PCR mix was performed in 96 wells with a single forward primer COI1_F: 5 -GAGGTGGTTTTATTGGTCTTGG-3 and a single reverse primer COI1_R: 5 -GGCCACCATCATACCAACAT. Schistosome ITS2 PCR mix was performed with a single forward primer ITS4_F: 5 -TAACAAGGTTTCCGTAGGTG AA-3 and a single reverse primer ITS5_R: 5 -TGCTTAAGTTCAGCGGGT-3 (Kane and Rollinson, 1994). The PCR mix was made up of 17.35 µL of ultra-pure water, 6 µL of buffer (colorless GoTag Flexi buffer 5X; Promega; Madison, Wisconsin, USA) 1.8 µL of 25 mM MgCl2 (Promega), 0.6 µL of 10 mM dNTPs mix (Promega), 1 µL of each 10µM primer, 0.25 µL GoTaq G2 Hot Start polymerase (Promega), and 2 µL of DNA extract, making a total volume of 30 µL for each PCR mix. The PCR thermal cycling conditions used was the same for all markers and was performed in (plate thermal cycler) a PerkinElmer 9600 Thermal Cycler (PerkinElmer, Waltham, MA, USA): pre-denaturing at 95 • C for 3 min, 45 cycles of 30 s at 95 • C (denaturing), 40 s at 56 • C (annealing), and 80 s at 72 • C (extending). This was followed by a final extending period of 2 min at 72 • C. The PCR product was stored in the refrigerator at 4 • C until used. Then, 4.5 µL of the product was mixed with 1.5 µL of a green loading dye to make 6 µL which was loaded into each well of a 1% agarose gel with 8 µL Midori dye using a multi-channel micro-pipette and 5 µL for size standard 100 bp ladder. This was analyzed by electrophoresis at 135 V for 30 min and transferred to the UV trans-illuminator where gel images were taken. The expected band size was 1000-1100 bp. Fifty-nine (59) samples were selected based on the quality of the amplicons. These successfully amplified PCR products were purified and sequenced on an Applied Biosystem Genetic Analyzer at Genoscreen, Lille, France.

Sequence Analysis
The Cox1 and ITS2 sequences were assembled separately and edited with a 4.5 sequencer version: (Gene Codes Corporation; Ann Arbor, MI, USA). The sequences were aligned using BioEdit Version 7.0.9 and ClustalW software. The aligned sequences were compared with the sequences in the GenBank Nucleotide Database for species designation: (https://www.ncbi.nlm.nih.gov/nucleotide/ accessed on 30 March 2022). The nuclear ITS2 region between S. haematobium and S. bovis differs at five polymorphic sites, hence the sequence chromatograms were checked at these SNPs to identify any possible heterozygosity (Supplementary Figure S2). We constructed a Cox1 gene phylogenetic tree only using S. bovis sequences because the Cox1 S. haematobium gene is known to be poorly variable and a phylogenetic study revealed only two clusters in all the areas of repartition of the parasites [42]. The phylogenetic tree was constructed using MEGA version 6.0.6 (Pennsylvania State University, Philadelphia, PA, USA) using an HKY + G nucleotide substitution model identified as the best model describing data. The support for tree nodes was calculated with 1000 bootstrap iterations. The phylogenetic analysis includes S. bovis sequences from various African countries retrieved from GenBank databases with a minimum length of 778 bp (see Supplementary Table S1). The tree was rooted in the S. haematobium haplotypes of the present study. All sequences were uploaded onto the GenBank database (OL840258-OL840278).

Microsatellite Genotyping
Microsatellite genotyping was performed on parasites from all sites except 5 and 6. A total of 1136 samples (Table 1) were individually genotyped with 18 microsatellite markers divided into two panels of 9 loci [43]. The multiplex PCR mix for each panel in two tubes was performed using the Qiagen ® multiplex PCR Kit (Qiagen, Hilden, Germany) according to the manufacturer's standard amplification protocol. The forward primers were fluorescently labeled using 6-FAM, VIC, NED, and PET dyes (Applied Biosystems, Foster City, California, USA). The PCR mix consists of 5 µL Qiagen MM 2X, 1 µL of 10X microsatellite primer mix, and 4 µL DNA extract making a final volume of 10 µL. The thermal cycling was performed in a plate thermocycler, PerkinElmer 9600 Thermal Cycler (PerkinElmer, Waltham, MA, USA): pre-denaturing at 95 • C for 15 min, 40 cycles of 30 s at 94 • C (denaturing), 90 s at 56 • C (annealing), and 60 s at 72 • C (extending). This was followed by a final extending period of 30 s at 60 • C [43]. The microsatellite PCR products were sent to Genoscreen, Lille, France for genotyping. Each microsatellite locus was visibly peak called with GS500Liz size standard (Applied Biosystem) and GeneMarker software. Eighty percent of our samples were successfully amplified by 14 loci and were used for result analysis while 4 markers (C131, Sh4, Sh8, and Sh15), which amplified less than 20% of the samples, were excluded.

Population Genetic Structure
Linkage disequilibria and departures from Hardy-Weinberg expectations were tested using exact tests (1200 permutations) adjusted for multiple tests using Bonferroni's correction as implemented in the FSTAT software 2.9.3.2 [44]. We analyzed the genetic variability of schistosomes from each study site by computing the expected heterozygosity (He), number of alleles (A), allelic richness (Ar), and the inbreeding coefficient (Fis) in each microsatellite's locus with FSTAT v.2.9.3.2 [44]. Heterozygosity (He) and allelic richness (Ar) between the populations were compared using the pairwise Friedman rank test followed by Nemenyi post hoc test.
Genetic structure was first assessed by calculating pairwise F ST values between sites according to [45] using FSTAT version 2.9.3.2. A possible link between geographic (in Km) and genetic distances (Fst) was assessed using the Mantel test. Second, we used the principal components analysis (PCA) implemented in Genetix [46]. Because we sampled several miracidia per patient, these miracidia are related, which could influence the genetic structure. In order to assess a possible bias of our sampling strategy, we performed PCA by randomly sampling two miracidia per patient. Third, we used the Bayesian clustering approach implemented in the Structure software to determine the uppermost level of genetic structure [47]. We tested the number of clusters from K = 1 to K = 12, by computing three runs for each cluster which is made up of 10 6 iterations after a "burn-in" period of 250,000 iterations with other parameters set by default and an admixture model. The mean logarithm probability for each cluster (K) was taken for the three runs with the corrsieve package in R. The ∆K-values were then computed in R to determine the probable cluster number from the total clusters (∆K) tested according to Evanno et al. [48], from which we identified K = 2 as the most probable genetic clusters. Lastly, an additional 10 runs were computed for K = 2 using 10 6 iterations and setting the same parameters as earlier described. The mean probability for a miracidium to belong to each cluster over the 10 runs was taken as Q-values, and we used Clumpp version 1.1.2 according to Francis R.M. [49], and Distruct version 1.1 according to Rosenberg N.A [50].

Conclusions
This study revealed that S. haematobium-bovis hybrids are predominant in Schistosoma eggs isolated in the urine samples of primary school-aged pupils in Nigeria. Our findings provide evidence that S. haematobium x S. bovis hybrids are common in Nigeria. Based on the high prevalence of S. haematobium x S. bovis hybrids, we advocate research priority on domestic and wild animals to investigate the role of zoonotic transmission.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/pathogens11040425/s1, Table S1. Accession numbers for Supple-mentary Figure S1; Table S2. Microsatellite database; Figure S1. Maximum likelihood phylogenetic tree built with 21 (2 S. haematobium NigeriaHap1&2, and 19 S. bovis NigeriaHap3-21) haplotypes from the present study and haplotypes from Cameroon, Benin, Senegal, Cote d'lvoire, Kenya and Tanzania from Genbank database. See supplementary Table S1 for AN database; Figure S2. The sequence chromatograms show the pure and mixed signal in the nuclear ITS2 marker. The double sequence chromatogram (heterozygous) showing bi-parental inheritance of the nuclear DNA; Figure S3. Population genetic structure graph assessed by principal component analysis using 2 miracidia by patient. Each sampling site is represented by a dot. The first and second axis of the PCA represent 46.2% and 16.9% respectively of the total variation in allele frequency.  Institutional Review Board Statement: This study was part of an epidemiological surveillance not submitted to an ethics committee.
Informed Consent Statement: Urine sample were collected only from children whose parents or legal guardians gave oral consent. Each participant whose urine sample was positive for the parasite was treated with a single oral dose of 40 mg/kg body weight of praziquantel.
Data Availability Statement: Datasets generated for this report can be found on NCBI database for sequences and for microsatellite database (Supplementary Table S2).