Ethics statement

Ethical approval for the cluster randomized trial was obtained from the Zanzibar Medical Research Ethics Committee (ZAMREC 0003/Sept/011), the “Ethikkomission beiber Basel” (EKBB) in Switzerland (reference no. 236/11), and the Institutional Review Board of the University of Georgia (project no. 2012-10138-0). The cluster randomized trial is registered at the International Standard Randomised Controlled Trial, Register Number (ISRCTN48837681). Further approval for the population genetic analyses was obtained from the St Mary’s Hospital Local Ethics Research Committee, R&D office, part of the Imperial College, London Research Ethics Committee (ICREC; (EC NO: 03.36. R&D No: 03/SB/033E)).

All aspects of sample collections for the cluster randomized trial were carried out in collaboration with staff of the Neglected Diseases Programme of the Zanzibar Ministry of Health (MoH) in Unguja and the Public Health Laboratory—Ivo de Carneri in Pemba. Shehia and school authorities were informed about the purpose and procedures of the study. Written informed consent was obtained from parents for all minors (children ≤16 years old) prior to recruitment. All participating adults provided their own signed consent. Participation was voluntary and individuals could withdraw or be withdrawn from the study at any time without obligation. All infected participants were offered treatment with praziquantel (40 mg/kg), donated by Merck KgA and facilitated by the WHO and Zanzibar Ministry of Health, as part of the MDA implemented by the Schistosomiasis Control Initiative (SCI) and MoH.

Shehia selection and intervention trial

The total population size of Zanzibar at 2012 was estimated at around 1.3 million, with 896,721 people residing in Unguja and 406,848 in Pemba. Unguja was composed of six districts further divided into 210 shehias, and Pemba had four districts divided into 121 shehias. Each shehia contains several villages, which varied considerably in population size (ranging from 482–26,275) and household number [23]. As described in detail elsewhere [18,19], 90 shehias (45 per island) were selected as intervention-clusters of the randomized trial, where 15 shehias per island were assigned to receive one of the three treatment arm interventions. Participants were randomly-selected, asked for consent to participate, and enrolled into the study from each school/shehia. Both populations, comparable across Unguja and Pemba, had received extensive MDA before the ZEST trial began in 2012 and thus “baseline” here does not represent PZQ-naïve populations [11] (and as detailed below).

Sample collection and microsatellite analyses

S. haematobium miracidia were collected from infected participants at baseline in January-March 2012 (before the onset of the ZEST interventions), with community-wide praziquantel treatment then implemented in April 2012. Miracidial collections were again collected from participants found infected four years into the ZEST programs in Jan-April 2016 (Year 5), following biannual treatment rounds in each of 2012–2015 (8 prior rounds in total), with CWT treatment provided in May 2016. Mid-morning urine samples were collected from each individual included in the randomized controlled trial and processed the same day, within the framework on the main ZEST programme. Samples were visually inspected for blood (macrohematuria) using a colour chart and for microhaematuria using reagent strips (Hemastix; Siemens Healthcare Diagnostics GmbH, Eschborn, Germany). Schistosoma haematobium eggs were identified and quantified by filtering 10ml of urine through a polycarbonate filter (Sterlitech, Kent, United States of America) and subsequently examined under a microscope by experienced laboratory technicians. Eggs were then concentrated from the remainder of all infected urine samples by filtration using a Pitchford funnel, rinsed and transferred into a clean Petri dish containing mineral water and exposed to light to facilitate hatching of miracidia. Individual miracidia were captured in mineral water and pipetted onto Whatman-FTA Indicating Classic cards (Whatman, Part of GE Healthcare, Florham Park, USA) for storage of DNA at room temperature (see [24,25]).

For available samples from each infected participant, the DNA from 24 randomly-selected miracidia (or all miracidia available where <24) was individually alkaline-eluted from the Whatman-FTA Indicating Classic cards (Whatman, Part of GE Healthcare, Florham Park, USA) as described in [24,25], and each sample was given a unique ID (anonymised personal identifier number). Each individual miracidium was genotyped using two previously developed S. haematobium multiplexed microsatellite PCR panels (Panel 1 and Panel 2) as described in [22] [25–27] (Table A in S1 Text). Initially the Panel 1 PCR was conducted on all DNA elutions, in a 96 well plate. Four microlitres of each reaction was visualized on a 2% gel red agarose gel and all positive reactions were cherry picked, using MICROLAB STAR Liquid Handling Workstation (Hamilton Robotics), for further fragment analysis. DNA elutions that gave successful positive panel 1 reactions were then also cherry picked, in a 96 well format, again using the MICROLAB STAR Liquid Handling Workstation (Hamilton Robotics) and analysed using the Panel 2 multiplex microsatellite PCR as described in [22]. This two-step approach was utilized to reduce the processing of samples with no or inadequate DNA, due to collection and processing errors inherent with these methodologies. Positive and negative (no DNA) controls were included with each set of 96 reactions. All positive microsatellite panel 1 reactions were diluted 1 in 10 before being denatured and injected at an optimal speed of 12 seconds into the Applied Biosystems 3130xl DNA analyser for analysis. Genotype data was exported and allele peaks were checked and edited using Geneious 6.1.4 before being placed into amplicon size “bins” and alleles exported for analysis. (Full details of the protocols followed are provided in S1 Text).

Data analyses

Microsatellite multilocus genotypes (MLGS) for each miracidium were compiled by combining data from microsatellite panels 1 and 2 using R (v.3.2.1.) [28]. Data for each miracidium were associated with data from the parasitological surveys including ID, year of collection (baseline or year 5), whether the individual was a child (6–8 or 9–12 years old) or adult, infection intensity as S. haematobium eggs per 10ml of urine, sex, age, and shehia of residence (adult) or schooling (child), island (Pemba or Unguja) and intervention arm (MDA only, MDA + snail control or MDA + behaviour change). Data on MDA coverage, shehia prevalence and infection intensity were available from the main trial [20]. Linear mixed-effects models with random intercepts (using the intensity in eggs per 10ml as the response variable, intervention arm and age as fixed effects, and shehias as random effects) were fitted to the 2012 baseline and 2016 Year 5 follow up data to explore differences in infection intensity between intervention arms, age and sex within and between islands.

Genetic diversity

Summary statistics for allelic richness (Ar), expected heterozygosity (He) and observed heterozygosity (Ho), which are measures of genetic diversity, and the inbreeding coefficient (F_ST) for each infra-population (the parasite population of each child or adult at each time point) were created using Powermarker [29] and FSTAT [30]. Allelic richness is suitable for the comparison of samples of different sizes since it rarefies to the smallest sample size in the dataset, and is thus particularly useful for schistosome population genetics where the number of miracidia collected is partly dependent on infection intensity. The alternative approach of excluding small sample size datasets would bias against lighter intensity infections, which are likely following successful PZQ treatment. The inbreeding coefficient is a measure of the degree of similarity between parasites within an individual where a lower value indicates less relatedness. Differences in Ar, He, Ho and F_ST between treatment arms and between baseline and Year 5 were investigated using linear regression models implemented in Stata (StataCorp LP, USA), controlling for age and sex of individuals, shehia/school, island and miracidia sample size. Models used robust standard errors to account for clustering of infra-populations/individuals within shehias. Two approaches to modelling were taken. Since there was minimal evidence of difference in genetic diversity indices between intervention arms at baseline, effects of the different intervention could be seen by direct comparison of the arms in the Year 5 dataset. Following this, the datasets from the two years were combined and change was modelled over time, using the baseline population as the reference population, allowing for interactions between time, island and arm to be investigated. Post-hoc testing of marginal means investigated whether there was a change for each of the arm-island populations between their own 2012 baseline and 2016 follow-up collections.

Parentage analysis, estimates of adult worm burden and fecundity

The estimation of full-relationships between the miracidia of each infra-population was carried out using COLONY version 2.0.6.1 [31]. COLONY implements a maximum likelihood algorithm when comparing different sibship configurations and also allows for genotyping error [32]. Sibship relationships were used to estimate the number of unique female parental genotypes present and the number of miracidia per sibship cluster within each individual host [33,34]. Unique parental genotypes were subsequently related to the minimum probable number of female schistosomes within a host using our recently developed statistical approach [35]. Inferred female worm burdens (where all outputs statistically accounted for the limited sample size of miracidia available), together with their individual infection intensity, were used to estimate the parasite reproductive success, measured as the mean number of eggs per 10 ml urine per female worm [35]. Linear mixed-effects models with random intercepts (using the inferred number of female worms within hosts as the response variable, island and age as fixed effects, and shehias as random effects) were fitted to the 2012 baseline and 2016 Year 5 data to explore differences in inferred numbers of female worms between islands and treatment arms. Linear mixed-effects models with random intercepts (using the intensity of infection as the response variable, inferred number of female worms within hosts, island and age as fixed effects, and shehias as random effects) were also fitted to explore differences in worm fecundity between different islands and between different shehias within the same island.

Persistent hotspots

There is no general consensus on how to define biological ‘hotspots’, although definitions which account for the failure to change over time, and/or high intensity of infection, are favoured [22,36,37]. We therefore defined biological hotspots here as either: (a) shehias where prevalence had not declined by at least 20% of its starting prevalence between 2012 baseline and 2016 Year 5 amongst the matched 9–12 year old age group (excluding three shehias with a very low starting prevalence (<3%) where stochastic events could have a large impact on definition due to the small number of infected individuals involved), or (b) shehias where the infection prevalence was  ≥ 15% in 9–12 year-old schoolchildren in at least one of the three cross-sectional parasitological surveys conducted in primary schools in 2012, 2013 and 2014 [36]. Using definition (a), nine shehias (out of 90) were defined as hotspots; using (b), 22 shehias (out of 90) were defined as hotspots. Our third categorization included shehias categorized as a ‘hotspot’ under both definitions, (a) and (b). Since all nine shehias defined as hotspots in definition (a) were included within definition (b) as well, our third categorization (c) included a total of nine shehias. Regression models were again used to determine if there were differences in mean allelic richness and inbreeding coefficients within shehias classified as hotspots or not, with robust standard errors to account for the clustering of infra-populations within shehias.

Population structure by island and persistent hotspot status

To gain further insights into the population structure of the miracidia, the microsatellite data were subjected to Principal Component Analysis (PCA) using the “adegenet” v2.1.1 package [38] in R (version 3.5.2). Each miracidium MLG was assigned to its island of origin and as either from a hotspot versus non-hotspot using the definitions as above.

Intensity of infection

Within those individuals where both epidemiological and parasitological material for population genetic analyses were obtained, infection intensity was significantly higher on Pemba compared to Unjuga across all age-groups and time points, most notably amongst the children (t = -8.28, p<0.0001, Fig 1). For urogenital schistosomiasis, infection intensities can be further classified into ‘high’, where ≥50 eggs/10ml of urine, and ‘low’ where 1–49 e/10 ml [39]. Thus the prevalence of heavy intensity infections was significantly higher on Pemba compared to Unjuga across all age-groups and time points (z = 5.63, p<0.0001). A significantly higher prevalence of heavy intensity infections was also found in adults compared to children of 6–8 years old within Pemba (z = 3.484, p<0.05). Although there was a general decrease in mean intensity of infection with increasing age, significant mean differences were found solely between children in Pemba of 6–8 years old and adults at baseline (t = -4.97, t = -4.97 p<0.0001; Fig 1), with even a non-significant upward trend amongst adults within Pemba. There was a significant decrease in infection intensity across all age groups in Pemba by the end of the study compared to baseline (t = -2.40, p = 0.01), consistent with the main ZEST findings, but this failed to reach significance in Unguja.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

A-B. Mean (± 95% confidence interval) intensity of infection (eggs/10ml urine) amongst individuals from different age-groups at baseline and Year 5, pooled across treatment arms but within islands for Pemba (a) and Unguja (b). In Pemba, the baseline data includes 53 children aged 6–8, 31 children aged 9–12 and 26 adults aged 13–75, and the Year 5 data includes 102 children aged 9–12 and 12 adults aged 20–30. Unguja’s baseline data consists of 32 children aged 7–8, 57 children aged 9–12 and 18 adults aged 20–54. Unguja’s Year 5 data includes 44 children aged 9–12 and 13 adults aged 20–46.

https://doi.org/10.1371/journal.pntd.0010419.g001

Genetic diversity

There was no difference in allelic richness between arms at year 1 baseline (MDA + behaviour change vs MDA only t = 0.67, p = 0.51; MDA + snail control vs MDA only; t = 1.58, p = 0.12) or in Year 5 (MDA + behaviour change vs MDA only t = 0.85, p = 0.40; MDA + snail control vs MDA only; t = 0.75, p = 0.46) (see Fig A in S1 Text). Over time, there was no evidence of an overall difference in allelic richness between the baseline and Year 5 collections (t = 1.44, p = 0.16), nor between the intervention arms (MDA + behaviour vs MDA only (t = 0.48, p = 0.63); MDA + snail control vs MDA only (t = 1.62, p = 0.11) or in their allelic diversity at Year 5 in comparison to the baseline MDA only population (Year 5 MDA + behaviour, t = 0.17, p = 0.87; Year 5 MDA + snail control, t = 1.20, p = 0.23) (Fig A in S1 Text). However, there was weak evidence of an interaction term between arms and their response to time (t = -1.77, p = 0.08) and post-hoc comparisons of the marginal means demonstrated that this was due to a three-way interaction between time, island and intervention arm (F_{2, 61} = 3.51, p = 0.04) with a trend towards a reduction in allelic diversity only being present in the Unguja MDA+snail control arm between baseline and Year 5 (t = -1.88, p = 0.07). Similar results were seen for He and Ho (Figs B and C in S1 Text).

Inbreeding was lower (i.e. outcrossing was higher) in Unguja (t = -2.03, p = 0.05), although there was no overall change in inbreeding coefficients over time (t = 1.19, p = 0.24), nor between the arms with additional interventions in comparison to the MDA only arm (MDA + behaviour, t = 0.66, p = 0.51; MDA + snail control, t = 0.82, p = 0.41). However, post-hoc comparisons of marginal means demonstrated a difference in the response of the intervention arms to time (F_2,72 = 5.19, p = 0.008), when each arm was compared to its own baseline population, with increased outcrossing evident in the snail arms. This was confirmed by analyses of each year independently, where there were no differences between arms or islands in Year 1, but a significantly reduced inbreeding coefficient amongst parasites from the Year 5 MDA+snail arm on Unguja (t = -0.23, p = 0.03) (Fig D in S1 Text).

Parentage analysis and estimates of female worm burden

Numbers of female worms per infrapopulation were estimated using sibship reconstruction and our recently developed statistical approach [35]. Despite the higher baseline infection intensities in Pemba island relative to Unjuga, across all arms, there were no significant difference in the mean inferred number of adult worms per infrapopulation (between the islands across both baseline and Year 5 (t = -0.894, p = 0.37) (Fig 2). In both islands, there was a significant decrease in the inferred numbers of female worms from baseline to Year 5 across all arms, indicative of successful control efforts reducing the mean adult worm populations (t = -10.21, p<0.001). However, there were no significant differences in the inferred minimum number of female worms between different treatment arms in both islands at either baseline or Year 5 (t = 0.50, p = 0.62).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

A-B. Mean (±95% confidence intervals) estimated number of Schistosoma haematobium female worms per host for (a) Pemba and (b) Unguja by treatment arm; where Pemba MDA arm n = 41 baseline, n = 22 year 5; Snail arm n = 40 baseline, n = 55 year 5; Behaviour arm n = 29 baseline, n = 37 Year 5; and where Unjuga MDA arm n = 20 baseline, n = 18 year 5; Snail arm n = 46 baseline; n = 8 year 5; Behaviour arm n = 43 baseline, n = 31 Year 5.

https://doi.org/10.1371/journal.pntd.0010419.g002

Parasite fecundity

From the estimation of adult worm burden per host, the inferred minimum mean number of eggs produced per female worm (per 10 ml urine) was significantly higher in Pemba than Unguja at baseline (32.9 (13.9–79.6) v 10.5 (2.7–42.2) in Unguja: t = -5.355, p<0.001) and Year 5 (14.6 (7.9–27.1) v 4.4 (1.5–12.8): t = -5.113, p<0.001) Fig 3). This is consistent with the observation that infection intensity was consistently higher in Pemba (Fig 1), despite the inferred number of female worms being similar between the two islands (Fig 2). There were not, however, any significant differences in fecundity within islands, in terms of either treatment arms or years, where all reduced in reponse to control irrespective of treatment arm (Fig 3), although there were significant differences between adults and the 5–8-year-old children at baseline (t = 4.769, p<0.001).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

A-B. Inferred mean (±95% confidence interval) fecundity (number of eggs per 10ml urine per Schistosoma haematobium female worm) in (a) Pemba and (b) Unguja by treatment arm and year; where Pemba MDA arm n = 41 baseline, n = 22 year 5; Snail arm n = 40 baseline, n = 55 year 5; Behaviour arm n = 29 baseline, n = 37 Year 5; and where Unjuga MDA arm n = 20 baseline, n = 18 year 5; Snail arm n = 46 baseline; n = 8 year 5; Behaviour arm n = 43 baseline, n = 31 Year 5.

https://doi.org/10.1371/journal.pntd.0010419.g003

Persistent hotspots

Fig 4A and 4B shows the mean fecundity (number of eggs per 10 ml urine per S. haematobium female worm) for each shehia at baseline and Year 5 estimated from the mixed effects model. Intensity levels at baseline were significantly higher in Pemba hotspots relative to Unjuga hotspots (t = -6.20, p<0.001), although there was no evidence of an interaction between islands and hotspots (t = -1.26, p = 0.22). We did not find significant differences in mean intensity levels (t = 0.66, p = 0.514) or fecundity (t = 0.227, p = 0.821) at baseline between individuals in shehias that went on to be designated as persistent hotspots. However, there was a non-significant trend in the inferred numbers of females at baseline between individuals in shehias that wernt on to be considered persistent hotspots (according to definition (b), which also includes all the hotspots in definition (a), t = 1.948, p = 0.057). Across both islands and both years, there was no difference in the mean number of adult worms per infrapopulation estimated from parentage analysis (t = 0.08, p = 0.94) and fecundity (t = 0.77, p = 0.44) in shehias designated as hotspots compared to responders. However, the estimated mean number of eggs per female worm was higher in some hotspots compared to the other shehias in both years (Fig 4A and 4B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4.

A-B. Mean fecundity of schistosomes (eggs per 10ml urine per female worm) within each shehia at Year 1 baseline (a) and Year 5 follow up (b) estimated using a mixed effects model with shehias as random effects. On the x-axis, shehias labelled in brown (bold) correspond to persistent hotspots defined by definition (b) where the infection prevalence was  ≥ 15% in 9–12 year-old schoolchildren in at least one of the three cross-sectional parasitological surveys conducted in primary schools in 2012, 2013 and 2014 [36]. Shehias labelled in black (bold) correspond to hotspots defined by both definition (b) and definition (a) (i.e., category (c) hotspots, see Material and Methods, Persistent hotspots) where prevalence had not declined by at least 20% of its starting prevalence between baseline and Year 5 amongst the matched 9–12-year-old age group. (Error bars not included for visual clarity and because, as a mixed effects model was used, few efficient methodologies to incorporate uncertainty into random effects parameters are available).

https://doi.org/10.1371/journal.pntd.0010419.g004

In Pemba, at baseline and Year 5, shehias Kangagani, Uwandani, Vitongoji and Pujini (located in the eastern part of the island), Wambaa (in the western part of the island) and Shumba Viamboni (at north), were found to have mean minimum worm fecundities ranging from 29.8 to 48.5 eggs/10ml urine/female worm, and were considered hotspots by both definitions (Fig 5). In Unguja, parasite mean fecundities ranging from 13.6 to 14.7 eggs/10ml urine/female worm were found in Kitope and Koani, which were also defined as hotspots by both definitions (c). These shehias are located in the centre and northwest of the island (Fig 5).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5.

A–D. Mean fecundity of schistosomes (eggs per Schistosoma haematobium female worm) and location of each shehia in 2012 Year 1 baseline (a and b) and 2016 Year 5 follow up (c and d), in Pemba (a) and (c) and Unguja (b) and (d). Shehias labelled in brown correspond to persistent hotspots by definition (b) and shehias labelled in black correspond to hotspots given the definitions (a) and (b) (see Material and Methods, Persistent hotspots for definitions). Note that, whilst the colour schemes are the same in both figures, the scales are different between different years, as fecundity ranged from 6.5–62 in Year 1 and 4–17 in Year 5). The maps were created in RStudio, using the tmap package: https://data.humdata.org/m/dataset/451bdd28-d06d-46ea-91c0-2e081f884395?force_layout=light No changes were made. The License for these data can be found at: https://data.humdata.org/m/dataset/451bdd28-d06d-46ea-91c0-2e081f884395?force_layout=light and are freely-available under the CC BY-IGO license.

https://doi.org/10.1371/journal.pntd.0010419.g005

Parasites also showed a clear genetic distinction between islands (Fig 6A), with also a potential separation by hotspot status within Pemba (Fig 6D and 6E), although with more overlap within Unguja (Fig 6B and 6C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6.

A-D. Population structure by island and persistent hotspot status using Principal Component Analysis (PCA) where each point represents an individual miracidium. Each miracidia MLG was assigned as either being from island and/or hotspots versus non-hotspots populations from baseline using the definitions as above. Scale represented as Euclidean distance, d = 0.5. Inertia ellipses indicate the distribution of the individuals from different groups (island or hotspot status). 6b –Unguja hotspot definition A. 6c –Unguja hotspot definition B. 6d –Pemba hotspot definition A. 6e –Pemba hotspot definition B.

https://doi.org/10.1371/journal.pntd.0010419.g006