Genomic data reveal a North-South split and introgression history of blood fluke populations across Africa

The human parasitic fluke, Schistosoma haematobium hybridizes with the livestock parasite S. bovis in the laboratory, but the frequency of hybridization in nature is unclear. Here, we analyze 34.6 million single nucleotide variants in 162 samples from 18 African countries, revealing a sharp genetic discontinuity between northern and southern S. haematobium. We find no evidence for recent hybridization. Instead the data reveal admixture events that occurred 257-879 generations ago in northern S. haematobium populations. Fifteen introgressed S. bovis genes are approaching fixation in northern S. haematobium with four genes potentially driving adaptation. Further, we identify 19 regions that are resistant to introgression; these are enriched on the sex chromosomes. These results (i) suggest strong barriers to gene flow between these species, (ii) indicate that hybridization may be less common than currently envisaged, but (iii) reveal profound genomic consequences of rare interspecific hybridization between schistosomes of medical and veterinary importance.

Methods
Sample collection: description, ethics, and identification: We used samples or data from three sources.  i) The first dataset was generated from samples provided by the Schistosomiasis Collection At the Natural History Museum1 which is housed at the Natural History Museum (London).  SCAN samples consisted of individual miracidia and cercariae preserved on Whatman FTA cards 2.  We analyzed 114 S. haematobium and S. bovis samples from 123 individual hosts (snails or humans) and 12 Africa countries. ii) In addition to the SCAN samples, we collected nine adult Schistosome worms, presumed to be S. bovis, from the intestines of routinely slaughtered cattle from meat vendors at three abattoirs located in Auchi, Benin City, and Enugu in Nigeria. In the laboratory, the mesenteric vessels of each purchased intestines were visually inspected for schistosome parasites. Adult schistosomes were recovered using forceps and washed in saline solution. Adult pairs were separated into males and females before being stored in 96% ethanol for subsequent DNA isolation analyses.  iii) Finally, for the third source of data we used whole genome sequence data from NCBI3-8.

Samples provided by the SCAN repository were originally collected in accordance with protocols approved by local, state, and national authorities, including the Ministry of Health. The Imperial College Research Ethics Committee (ICREC) at Imperial College London, in conjunction with ongoing Schistosomiasis Control Initiative (SCI) activities, provided additional ethical guidance for samples collected through the CONTRAST program. Ethical clearance and study protocols for Nigerian samples were approved by the National Health Research Ethics Committee of Nigeria (NHREC) (protocol number: NHREC/01/01/2007– 30/10/2020 and approval number: NHREC/01/01/2007– 29/03/2021) and the Institutional Review Board (IRB) of University of Texas Health, San Antonio Texas, United States of America (protocol number: HSC20180612H). Informed consent was obtained from all participants, with processes tailored to ensure understanding and voluntary participation. All data were anonymized to protect participant privacy, and schistosomiasis-positive individuals were treated with a single dose of praziquantel (40 mg/kg). For livestock parasite collection, approval was secured from local veterinarians. No animals were euthanized for research purposes; Schistosoma samples were collected during routine activities at abattoirs. Further details on collection methods, ethical approvals, and data availability for public samples can be found in their respective publications documented in Supplemental Data File 1.

Provisional species identifications were assigned to cercariae and miracidia based on sampled host.  For example, miracidia hatched from eggs collected from human urine samples were assumed to be S. haematobium while miracidia hatched from eggs in cattle feces were assumed to be S. bovis. Cercariae collected from snails were identified by Sanger sequencing the mitochondrial cox1 region and the ribosomal internal transcribed spacer (ITS) rDNA.  The mitochondria was genotyped at the cox1 gene using a multiplex PCR that contains a standard forward primer (Asmit1: forward [TTT TTT GGT CAT CCT GAG GTG TAT]) and species specific reverse primers (SbR: reverse [CAC AGG ATC AGA CAA ACG AGT ACC], ShR: reverse [TGA TAA TCA ATG ACC CTG CAA TAA]).  Amplicons were visualized on a 2% gel.  Larger amplicons (543 bp) indicated S. haematobium and smaller amplicons (306 bp) were diagnostic for S. bovis9.  The ribosomal ITS sequence was amplified using the ETTS1 (TGC TTA AGT TCA GCG GGT) and ETTS2 (TAA CAA GGT TTC CGT AGG TGA A) primers 10.  Amplicons were Sanger sequenced and the resulting fragments were assigned to species based on comparison to reference samples11. Downstream genetic analysis with whole genome SNVs was used to confirm and reassign species identifications where necessary.

Library prep and sequencing: DNA from single parasites stored on FTA cards was subjected to whole-genome amplification (WGA).  Single miracidia2 isolated by punching the FTA card into a sterile tube.  We used the GenomiPhi V2 DNA Amplification kit (Sigma-Aldrich: Cat. No. GE25-6600-31) and the recommended protocols to amplify the schistosome DNA. DNA was extracted from single male adult S. bovis worms using the DNeasy® Blood and Tissue kit (Qiagen: Cat. No.69504) before subsequent WGA. We quantified amount of schistosome DNA in each WGA sample by real time quantitative PCR (qPCR) reactions using the single copy gene α-tubulin 1 gene markers primers (S. haematobium: forward [GGT GGT ACT GGT TCT GGT TT], reverse [AAA GCA CAA TCC GAA TGT TCT AA]; S. bovis: forward [ATG GCC TCG TTA TCA ACC AT], reverse [TGG CCT CGT TAT CAA CCA TA]) 2.  The PCR reaction was denatured at 95 °C for 10 minutes, followed by 40 cycles of 95 °C for 15 seconds and 60 °C for 1 minute. A standard curve was generated using six dilutions of α-tubulin PCR product spanning 1.29 × 10¹ to 1.29 × 10⁷ copies/µL. DNA sequencing libraries were generated from 500 ng of DNA per sample using the KAPA Hyperplus kit (Roche: Cat. No. 07962401001) protocol with the following modifications: i) enzymatic fragmentation at 37°C for 10 minutes, ii) adapter ligation at 20°C for an hour, and iii) 4 cycles of library PCR amplification. After qPCR quantification of each library with KAPA Library Quantification kits (Roche: Cat. No. 07960140001), samples with similar concentrations were combined into pools for sequencing at 4nM, while samples with disparate concentrations were equalized in 10 mM Tris-HCl pH 8.5 before pooling. Libraries were sequenced with 150 bp paired-end reads on two Illumina NovaSeq flowcell. All resulting reads were deposited in the NCBI Short Read Archive under BioProject PRJNA636746 and are documented in Supplemental Data File 1.

Read filtering and Mapping: Raw reads were quality trimmed with trimmomatic v0.39 12 using the following parameters: LEADING:10, TRAILING:10, SLIDINGWINDOW:4:15, MINLEN:36, ILLUMINACLIP:2:30:10:1:true.  This command removed low quality bases at the beginning and ends of the reads, removed portions of the read where quality dropped below a minimum threshold, trimmed adapter sequences and discarded reads <36 nts. We then mapped the trimmed reads to the Egyptian-strain S. haematobium reference genome, GCF_000699445.35, with BBMap v38.1813. On average the S. haematobium and S. bovis (GCA_944470425.1) genome assemblies are ~97% similar across their genomes8 which should minimally affect reference biases when mapping short reads.  However, to avoid reference biases we used the ‘vslow’ and ‘minid=0.8’ options with BBMap and discarded ambiguously mapping reads (‘ambig=toss’).

Genotyping, phasing, and filtering: Mapped reads were sorted with SAMtools v1.1314 and checked for duplicates with GATK v4.2.0.0’s15 mark_duplicates. Then single nucleotide variants (SNVs) were genotyped with HaplotypeCaller and GenotypeGVCFs. To make the dataset more manageable, we genotyped each chromosome separately using the -L option.  Next, we removed all indels and hard filtered SNVs based on QualByDepth (QD < 2.0), RMSMappingQuality (MQ < 30.0), FisherStrand (FS > 60.0), StrandOddsRation (SOR > 3.0), MappingQualityRankSumTest (MQRankSum < -12.5), and ReadPosRankSumTest (ReadPosRankSum < -8.0) with GATK’s VariantFiltration. We removed multi-allelic sites, and sites with genotype quality (GQ) <20 or read depth (DP) <8 with VCFtools v0.1.1616.   After these filters were applied we removed genomic sites that were genotyped in ≤50% of individuals and then any individuals that were genotyped at ≤50% of sites.     

References

1          Emery, A. M., Allan, F. E., Rabone, M. E. & Rollinson, D. Schistosomiasis collection at NHM (SCAN). Parasites & vectors 5, 185, doi:10.1186/1756-3305-5-185 (2012).

2          Le Clec'h, W. et al. Whole genome amplification and exome sequencing of archived schistosome miracidia. Parasitology 145, 1739-1747, doi:10.1017/s0031182018000811 (2018).

3          Rey, O. et al. Diverging patterns of introgression from Schistosoma bovis across S. haematobium African lineages. PLoS pathogens 17, e1009313, doi:10.1371/journal.ppat.1009313 (2021).

4          Platt, R. N. et al. Ancient Hybridization and Adaptive Introgression of an Invadolysin Gene in Schistosome Parasites. Molecular biology and evolution 36, 2127-2142, doi:10.1093/molbev/msz154 (2019).

5          Stroehlein, A. J. et al. Chromosome-level genome of Schistosoma haematobium underpins genome-wide explorations of molecular variation. PLoS pathogens 18, e1010288, doi:10.1371/journal.ppat.1010288 (2022).

6          Comparative genomics of the major parasitic worms. Nature genetics 51, 163-174, doi:10.1038/s41588-018-0262-1 (2019).

7          Young, N. D. et al. Whole-genome sequence of Schistosoma haematobium. Nature genetics 44, 221-225, doi:10.1038/ng.1065 (2012).

8          Oey, H. et al. Whole-genome sequence of the bovine blood fluke Schistosoma bovis supports interspecific hybridization with S. haematobium. PLoS pathogens 15, e1007513, doi:10.1371/journal.ppat.1007513 (2019).

9          Webster, B. L., Rollinson, D., Stothard, J. R. & Huyse, T. Rapid diagnostic multiplex PCR (RD-PCR) to discriminate Schistosoma haematobium and S. bovis. Journal of helminthology 84, 107-114, doi:10.1017/s0022149x09990447 (2010).

10        Kane, R. A. & Rollinson, D. Repetitive sequences in the ribosomal DNA internal transcribed spacer of Schistosoma haematobium, Schistosoma intercalatum and Schistosoma mattheei. Molecular and biochemical parasitology 63, 153-156, doi:10.1016/0166-6851(94)90018-3 (1994).

11        Pennance, T. et al. Interactions between Schistosoma haematobium group species and their Bulinus spp. intermediate hosts along the Niger River Valley. Parasites & vectors 13, 268, doi:10.1186/s13071-020-04136-9 (2020).

12        Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114-2120, doi:10.1093/bioinformatics/btu170 (2014).

13        Bushnell, B. BBMap: a fast, accurate, splice-aware aligner. (Lawrence Berkeley National Lab.(LBNL), Berkeley, CA (United States), 2014).

14        Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078-2079, doi:10.1093/bioinformatics/btp352 (2009).

15        McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20, 1297-1303, doi:10.1101/gr.107524.110 (2010).

16        Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156-2158, doi:10.1093/bioinformatics/btr330 (2011).