Yogyakarta City in south-central Java, Indonesia, has a population of 422,732 in an area of 32 km ², with 14 administrative districts comprising 45 kelurahans (urban villages). Seven kelurahans (total population 64,599; area 4.9 km ²) on the north-western perimeter of the city were selected as the site for scaled deployment of Wolbachia-carrying mosquitoes ( Figure 1), on the basis of logistical and operational feasibility. Three kelurahans (population 33,535; 3.1 km ²) on the south-eastern perimeter were selected a priori as an untreated control area, on the basis of comparable socio-demographic characteristics ( Table 1) and historical dengue incidence.

A Public Acceptance Model ¹⁷ was applied in engaging with the local community prior to Wolbachia deployments, and throughout the release and monitoring periods. Key elements of this approach included: meetings with key stakeholders and community leaders; meetings and ongoing regular communication with existing community reference groups at the village, city and provincial level; a communications campaign through social media, traditional media, mobile billboards, and community events; a household-based survey to evaluate awareness and acceptance prior to releases; and a ‘stakeholder enquiry system’ to receive and respond to any issues arising from stakeholders or community members.

A household survey of baseline community acceptance of WMP’s Wolbachia method was undertaken in November 2015 prior to intensive community engagement, followed by a second survey in June 2016 to evaluate acceptance prior to commencing releases. Two-stage cluster random sampling was used to select respondent households from across Yogyakarta City: n=587 in the 2015 survey (of which n=96 were from the seven urban villages in the quasi-experimental intervention area) and n=862 in the 2016 survey (n=180 from the intervention area). The sample size in the 2016 pre-release survey was sufficient to estimate community acceptance in the quasi-experimental intervention area with a margin of error of 10%, under conservative assumptions of a true community acceptance rate of 50% and a design effect of 2 due to the clustered sampling. Administrative wards were the primary sampling unit, and 12 households were randomly selected per ward using a list of households provided by the ward head as the sampling frame. The survey respondent was the head of household or another adult family member; there were no other inclusion or exclusion criteria. Respondents were asked about their level support for Wolbachia releases as a method for controlling dengue in their community. The survey participation rate was > 95% in each survey. During the 2015 baseline survey 67% of respondents from the quasi-experimental intervention area were supportive of the release of Wolbachia-infected Ae. aegypti as a new technology for controlling dengue, 14% were not supportive, and 20% were neutral or did not know. In the 2016 pre-release survey, acceptance in the intervention area had increased to 79%; 9% were not supportive and 11% were neutral or did not know.

The stakeholder enquiry system was managed by a dedicated staff member and enquiries could be logged by phone call, SMS text, WhatsApp, email or direct communication with WMP field staff. In total we had 216 enquiries during the release phase, relating to: details of release activities (location, timing, methods) (35%), requests for face-to-face meetings with community leaders (25%), expressions of support for the Wolbachia technology (18%), reports of dengue cases in the release area (7%), reports from householders on the condition of mosquito release containers (MRCs) (6%), and refusal to host an MRC (6%).

An existing colony of local Ae. aegypti containing the wMel Wolbachia strain, created in 2013 for pilot releases in Bantul District of Yogyakarta Province ¹⁹ , was used as the founder colony for the releases described here. It was backcrossed for three generations with wild-type males collected from the study intervention area. Four cages with 150 wMel-carrying female Ae. aegypti and 150 wild-type males were maintained for two gonotrophic cycles, and blood-fed by human volunteers. One hundred larvae were randomly selected from each cage and screened for Wolbachia, and only cages with 100% Wolbachia prevalence were maintained for the second backcrossing. An open colony was maintained by adding 10% of wild-type Aegypti after the second backcrossing. Insecticide resistance testing was conducted as described in 22, except that test specimens were adult females derived from egg, rather than larval, collections. This demonstrated equivalent insecticide resistance profiles between the colony and wild material ( Figure 2).

Mosquitoes were blood-fed by human volunteers once a week for two gonotrophic cycles, as per previous protocols ¹⁷ . All volunteer blood-feeders were afebrile and free of clinical signs or symptoms of any arbovirus infection both at the time of blood-feeding and for three days thereafter. Four days after blood feeding, mosquito eggs were collected by placing oviposition strips in each cage for 3 days. Eggs were dried slowly and stored in sealed containers until releases.

For quality assurance (QA) of the mosquito colony, a sample of ~1% of eggs from the parent generation of the release material was screened for Wolbachia by qualitative PCR Taqman assay on a Roche LightCycler 480. The pre-specified minimum acceptable Wolbachia prevalence was 97%. In each gonotrophic cycle, nine blood-fed female mosquitoes per human volunteer were tested for infection with DENV-1–4, chikungunya and Zika virus by qRT-PCR as previously described ²³ .

For QA of the release material, eggs from 10 oviposition strips were hatched each release week, with a minimum threshold hatch rate of 80%. Each release week, 10% of all release containers were selected randomly for QA, to estimate numbers of mosquitoes released. Eggs from these QA containers were manually counted under a microscope prior to field deployment. Release success was evaluated in the same QA containers at the time of servicing two weeks later (see below). A container was classified as ‘failed’ if: 1) it was lost; 2) it was dry; 3) pupae skin n<25; 4) dead adult mosquitoes >10; 5) wing carcass >20; or 6) remaining larvae and pupae >30.

Wolbachia-carrying mosquitoes were released as eggs using mosquito release containers (MRCs). These were 2-litre plastic buckets each containing one oviposition strip with 100–150 eggs, Tetra Pleco Wafers fish food (Tetra GmbH, Germany), and 1 litre of water. MRCs were covered and placed outside houses, protected from direct sun and rain. Holes drilled near the top of the bucket walls allowed adult mosquitoes to escape. Releases occurred between August 2016 and March 2017, with 13–15 rounds of releases in each kelurahan. Releases stopped in each kelurahan when the prevalence of Wolbachia in field-caught mosquitoes was >60% for three consecutive weeks releases ( Table 2). MRCs were reset every two weeks. An MRC was placed in 1–2 randomly selected locations within each 50×50 m grid square across the intervention area. Permission was obtained from property owners to place MRCs on private property.

Prevalence of Wolbachia in the local Ae. aegypti population was monitored by weekly collection of adult mosquitoes via a network of 89 BG Sentinel traps (Biogents, Germany). The median (range) trap density was 15.7 (13.0–20.0) BG/km ² in the intervention area and 3.6 (2.0–4.0) BG/km ² in the non-release area. Mosquitoes were demobilised at -20°C for ≥1 hour, then identified by morphological features. The number of mosquitoes caught in each BG trap was recorded by species, sex, and in total. Ae. aegypti were stored at -20°C in 80% ethanol until Wolbachia screening.

Field-caught Ae. aegypti were screened for wMel Wolbachia by qualitative PCR Taqman assay on a Roche LightCycler 480. The qPCR conditions consisted of a denaturation step at 95°C for 5 minutes followed by 45 cycles of PCR (denaturation at 95 °C for 10 seconds, annealing at 60 °C for 15 seconds, and extension at 72 °C for 1 second with the single acquisition) followed by a cooling down step at 40°C for 10 seconds. Specific primers targeting the gene encoding Ae aegypti Rps17 and wMel WD0513 were used as previously described ²⁴ , but with replacement of the Cy5-BHQ3 fluorophore-quencher pair in the wMel probe with the fluorophore-quencher LC640-IowaBlack (Integrated DNA technologies) ²⁵ . Testing was weekly when Wolbachia prevalence was <80%; two-weekly when ≥80% and four-weekly when ≥90%.

Data on hospitalized dengue haemorrhagic fever (DHF) cases (ICD-10 code A91; International Classification of Diseases 10th revision) were obtained from the Yogyakarta District Health Office (DHO), for January 2006–September 2019. Data on hospitalised dengue fever (DF) cases (ICD-10 code A90) were also obtained for January 2017–September 2019, as DF notification only began in 2017. Data were monthly case counts by kelurahan of residence. Population data were obtained from the Indonesian Bureau of Statistics. Anonymised results of dengue rapid diagnostic tests (RDTs) (SD Bioline Dengue Duo, Abbott, USA) performed in 18 primary health care clinics across Yogyakarta City from March 2016–September 2019 were obtained directly from the clinics. RDT results were extracted for participants resident in the quasi-experimental study intervention or control area, based on the kelurahan of residence where available, otherwise the kelurahan of clinic location.

The crude dengue incidence rate ratio in the pre-intervention and post-intervention periods was calculated as the aggregate number of DHF cases divided by the aggregate person-months, in the intervention area versus the control area. Crude dengue incidence was similarly calculated using an endpoint that included both DF and DHF case notifications, but only in the post-intervention period due to the lack of DF reporting prior to January 2017. Spearman’s rank correlation test was used to evaluate the correlation in monthly notified dengue incidence between the intervention and control areas, in the pre-intervention period. The Wolbachia intervention effect was estimated using controlled interrupted time series analysis ^{17, 18, 20} . Negative binomial regression was used to model monthly DHF case counts in the aggregate intervention and control areas, with an offset for population size. The primary analysis, defined a priori, included data from January 2006 until March 2019, two years post-release. A secondary analysis included an additional six months of data to September 2019. Annual population estimates were used for 2006–2013 and 2015–2017; 2013 population estimates were used for 2014 and 2017 estimates were used for 2018 and 2019, due to unavailability of data for those years. Seasonal variability in dengue incidence was controlled using 6-monthly flexible cubic splines. A binary ‘group’ variable indicated the study arm (intervention or control). A binary ‘treatment’ variable distinguished the pre-intervention period (up to the end of Wolbachia deployments in the last release area) and the post-intervention period, and the intervention effect was estimated from the interaction between the ‘group’ and ‘treatment’ variables. This allows explicitly for a level change in the outcome (dengue case incidence) in both intervention and control arms in the post-intervention period, for example in a scenario where other secular effects coincident with the Wolbachia deployments may have influenced dengue incidence independently of Wolbachia. The unavailability of DF data prior to January 2017 precluded an ITS analysis with a combined DF/DHF endpoint. Robust standard errors were used for all models. Fisher’s exact test was used to compare the proportion of patients with an NS1-positive RDT result between intervention and control areas, in the pre-intervention period and post- Wolbachia deploymnets. Analyses were performed using Stata® statistical software package version 14.2 (StataCorp, College Station, TX).

Power was estimated using 1000 simulated datasets drawn from a negative binomial distribution fitted to an 11-year time series (2006–2016) of monthly DHF case notifications from the intervention and control areas prior to Wolbachia deployment. Post-intervention time periods of 1, 2, 3 or 4 years were simulated, with the pre-intervention period fixed at 7 years (maximum total simulated time series of 11 years). Dengue case numbers were not modified for the untreated control area, and for the Wolbachia intervention area were either kept at baseline values (for the simulation at the null; i.e. relative risk (RR)=1) or reduced proportionately (for simulations of intervention effects of RR=0.7, 0.6, 0.5, 0.4, 0.3). For each of these five ‘true’ effect sizes, applied to each of the 1000 simulated time series, the ‘observed’ effect size was calculated from a negative binomial regression model of monthly case counts in the intervention and control areas, as described above. Power was estimated from the proportion of 1000 simulated scenarios in which a significant intervention effect (p<0.05) was observed ( Figure 3). With seven years of pre-intervention data and one year post-intervention follow-up, there was 80% power to detect a ≥30% reduction in dengue incidence (incidence rate ratio (IRR)=0.7). There was only a marginal gain in power from >1 year of follow-up, if the true reduction in dengue incidence is greater than 30%. Power calculations were performed using R version 3.5.0 (R Foundation for Statistical Computing, Austria) with the “MASS” package ²⁶ used for estimation of the regression models.

An independent risk assessment conducted by the Indonesia Ministry of Research, Technology and Higher Education in 2016 concluded that over the next 30 years, there is a negligible risk of harm as a result of releasing Wolbachia-infected Ae. aegypti ²⁷ .

Approval to release Wolbachia mosquitoes was obtained from the provincial and city governments of Yogyakarta prior to releases. Ethical approval was obtained from the Universitas Gadjah Mada Faculty of Medicine, Public Health and Nursing Ethics Committee (approval number KE/FK/105/EC/2016) to release Wolbachia mosquitoes, for blood-feeding the mosquito colony on human volunteers, and to access non-identifiable aggregate data on monthly notified DHF and DF case numbers from the Yogyakarta DHO. Approval to conduct the household surveys of community acceptance was also obtained from the UGM Ethics Committee (approval number KE/FK/1023/EX/2015). Written agreement to share dengue surveillance data and non-identifiable individual dengue RDT results from puskesmas clinics was obtained from Yogyakarta DHO.

Verbal and written consent was obtained from heads of households for participating in the baseline surveys, and for hosting a BG trap or MRC. BG hosts were compensated 50,000 IDR for the cost of powering traps.

Longitudinal entomological monitoring demonstrated rapid establishment of Wolbachia in the intervention areas, a continuous increase in Wolbachia prevalence in trapped Ae. aegypti throughout the first year post-release, and persistence at a very high prevalence ever since ( Figure 4). The median Wolbachia prevalence was 73% (range 67–92%) one week after releases stopped, and 100% (96–100%) two years post-deployment. In the control areas, single Wolbachia-infected Ae. aegypti mosquitoes were detected on 11 occasions, but there has been no evidence of Wolbachia establishment.

In the decade prior to Wolbachia releases, dengue outbreaks occurred annually in both the intervention and control areas ( Figure 4). A median of 125 DHF cases were notified each year in the intervention area (range 65–308), corresponding to a median annual incidence of 169 per 100,000 population (range 98–477 per 100,000). In the control area, a median of 44 DHF cases were notified each year (range 19–159) corresponding to a median annual incidence of 130 per 100,000 (range 59–477 per 100,000). Per-capita dengue incidence in the intervention area was on average 15% higher than the control area during the pre-intervention period (crude IRR 1.15, p = 0.002). Monthly dengue incidence in the intervention and control areas was highly correlated over time (Spearman’s rho = 0.75, p<0.001), supporting the validity of the control area as a counterfactual for evaluating the epidemiological impact of Wolbachia releases ( Figure 5).

During two years following the completion of Wolbachia deployments (April 2017–March 2019), 34 DHF cases were notified from the intervention area and 53 from the control area, corresponding to 67% lower crude incidence in the Wolbachia-treated area (26 vs 79 cases per 100,000 person-years; IRR = 0.33, p<0.001).

In an interrupted time series (ITS) analysis of monthly DHF case notifications January 2006 – March 2019, which adjusts for baseline differences between the intervention and control areas and seasonal and inter-annual time effects, this translates to a 73% reduction in notified DHF incidence (95% CI 49%,86%; p<0.001) associated with the Wolbachia intervention.

In an exploratory analysis including six months of additional post-intervention data to September 2019, the DHF incidence in intervention vs control areas was 30 vs 115 per 100,000 person-years (crude IRR = 0.26, p<0.001) which translates in the regression model to a 76% reduction in dengue incidence (95% CI 60%,86%; p<0.001) associated with the Wolbachia intervention

During the first 2.5 years post-intervention (April 2017–September 2019), 97 cases (DF + DHF) were notified from the intervention area and 137 from the control area (60 vs 163 per 100,000 person-years; IRR = 0.37, p<0.001) corresponding to 63% lower crude incidence in the Wolbachia-treated area, using a combined endpoint of DF and DHF case notifications. The lack of pre-release data on hospitalised DF cases precluded an ITS analysis with this combined endpoint.

Prior to Wolbachia deployments (March 2016 – March 2017), there was no difference between intervention and control areas in the proportion of tested patients with an NS1-positive RDT result (127/447 (28%) in intervention area vs 106/408 (26%) in control area; Fisher’s exact p=0.44; Figure 5; raw data available on Figshare ³² ). By comparison, in the post-intervention period (April 2017–September 2019) NS1 positivity was significantly lower in the intervention area than the control area (6/568 (1%) vs 61/429 (14%); Fisher’s exact p<0.001; Figure 6).

Figshare:

Figure 2_IRtesting.xlsx https://doi.org/10.6084/m9.figshare.12026943.v3 ³³

Figure 3_PowerPlot.xlsx https://doi.org/10.6084/m9.figshare.12026937.v3 ³⁴

Figure 4_MosquitoCollections.xls https://doi.org/10.6084/m9.figshare.12026940.v3 ³⁵

Figure 5_6_DengueCaseData.xls https://doi.org/10.6084/m9.figshare.12199688.v3 ³⁶

Wolbachia community acceptance survey. https://doi.org/10.6084/m9.figshare.12199754.v1 ³⁷

The dengue case notification data and NS1 rapid diagnostic test results were provided to us by the Yogyakarta District Health Office and individual Yogyakarta puskesmas clinics for the purpose of the current study, and the conditions of the release of the raw data do not permit further sharing to a third party. Data aggregated by month and kelurahan is available from the Figshare link above. Applications for access to raw data for research purposes or data reanalysis can be made directly to the Head of the Yogyakarta City Health Office at kesehatan@jogjakota.go.id.

Data deposited with Figshare are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).

STROBE checklist for 'Reduced dengue incidence following deployments of Wolbachia-infected Aedes aegypti in Yogyakarta, Indonesia: a quasi-experimental trial using controlled interrupted time series analysis' https://doi.org/10.6084/m9.figshare.12026955.v2 ³⁸