Local introduction and heterogeneous spatial spread of dengue-suppressing Wolbachia through an urban population of Aedes aegypti

Dengue-suppressing Wolbachia strains are promising tools for arbovirus control, particularly as they have the potential to self-spread following local introductions. To test this, we followed the frequency of the transinfected Wolbachia strain wMel through Ae. aegypti in Cairns, Australia, following releases at 3 nonisolated locations within the city in early 2013. Spatial spread was analysed graphically using interpolation and by fitting a statistical model describing the position and width of the wave. For the larger 2 of the 3 releases (covering 0.97 km2 and 0.52 km2), we observed slow but steady spatial spread, at about 100–200 m per year, roughly consistent with theoretical predictions. In contrast, the smallest release (0.11 km2) produced erratic temporal and spatial dynamics, with little evidence of spread after 2 years. This is consistent with the prediction concerning fitness-decreasing Wolbachia transinfections that a minimum release area is needed to achieve stable local establishment and spread in continuous habitats. Our graphical and likelihood analyses produced broadly consistent estimates of wave speed and wave width. Spread at all sites was spatially heterogeneous, suggesting that environmental heterogeneity will affect large-scale Wolbachia transformations of urban mosquito populations. The persistence and spread of Wolbachia in release areas meeting minimum area requirements indicates the promise of successful large-scale population transformation.


Introduction
Dengue fever is the most common arboviral disease affecting humans [1]. Over 2,500,000,000 people live in dengue-afflicted regions, and dengue incidence is increasing at an alarming rate in tropical and subtropical countries [2]. A number of other arboviruses also represent emerging disease risks, including chikungunya and Zika, the latter being associated with a recent explosive epidemic in South America [3,4]. The main approach to controlling these diseases has been suppression of the principal mosquito vector, Ae. aegypti, either through source reduction or insecticide-based control programs. Given the increasing incidence of Ae. aegypti-associated human disease, it is clear that current control measures are insufficient. In response to this problem, a number of new control approaches are currently being developed and tested [5,6,7,8,9,10].
In contrast to control efforts that require repeated population suppression, the Eliminate Dengue Program (http://www.eliminatedengue.com/program) aims to modify populations using long-lasting local introductions of a dengue-inhibiting Wolbachia into naturally uninfected populations of Ae. aegypti. The strain, wMel, was transferred from Drosophila melanogaster into laboratory-raised Ae. aegypti, who inherit the infection maternally [11,12]. Following introgression of the infection into a native genetic background, Wolbachia-infected mosquitoes are released into the field to mate with wild uninfected mosquitoes, and wMel frequency increases through cytoplasmic incompatibility (CI) [13,14]. CI describes the fact that uninfected females mated with Wolbachia-infected males produce inviable embryos. In Ae. aegypti, this is believed to occur in 100% of these incompatible crosses [11]. In contrast, infected females can mate with either infected or uninfected males and produce almost 100% infected progeny. CI greatly reduces the relative fitness of uninfected females when infected males are common and drives rapid establishment of Wolbachia in isolated mosquito populations [14], given that there is no mating bias against wMel-infected Ae. aegypti [15].
Although wMel-infected females receive a frequency-dependent relative fitness advantage from CI, they also suffer from frequency-independent fitness costs, including decreases in fecundity and larval competitive ability [16,17,18,19]. Thus, CI does not produce a net fitness advantage while wMel is rare, resulting in dynamics analogous to those produced by an Allee effect in ecology [20,21] and by natural selection on a locus (or alternative karyotypes) in which heterozygotes are less fit than either homozygotes (i.e., underdominance, [22,23,24]). The interaction of the frequency-dependent advantage associated with CI and the frequencyindependent cost(s) produces "bistable dynamics" with a threshold frequency of infection (denotedp) below which the infection will be locally eliminated and above which frequencies systematically increase [25,26,27].
Curtis [23] first proposed transforming pest populations by introducing translocations that are expected to show bistable dynamics (cf. [28]). The bistable model for Wolbachia spread was introduced by Turelli and Hoffmann [29] to explain the rapid spread of wRi, a CI-causing Wolbachia variant, through California populations of Drosophila simulans. Although this interpretation of wRi dynamics has now been challenged by more recent data on the spread of natural Wolbachia infections [30], 3 lines of evidence nevertheless support bistability of wMeltransinfected Ae. aegypti [31]: (1) frequency dynamics from the original field releases [14], (2) direct experimental evidence for lower fecundity and viability [19,32], and (3) new data showing that persistent influx over 2 years of wMel-infected Ae. aegypti into a relatively isolated population has not led to establishment of wMel there [31].
In order for the invasion to spread spatially under bistability, new uncolonised areas must receive infected immigrants at a rate high enough to be pushed past the threshold frequency, p. Under the dynamics produced by CI-inducing Wolbachia, spatial spread is expected in a habitat with relatively homogeneous population densities ifp is below a critical value near 0.5 [20,29]. For wMel in Ae. aegypti near Cairns,p is thought to be moderate (p % 0:2 À 0: 35) because of its relatively low fitness costs and near-perfect maternal transmission [11,14,31].
Previously, wMel-infected Ae. aegypti released in 2 relatively isolated communities in Northern Queensland, Australia (Gordonvale and Yorkeys Knob), colonised each area rapidly [14], and the infection has persisted at high frequency (>90%) at both sites [18]. Moreover, wMel continues to show strong blockage of dengue transmission in laboratory-challenged mosquitoes derived from field collections [33]. Here we present data from 3 subsequent releases of wMel-infected Ae. aegypti in Cairns, Northern Queensland, a city with about 150,000 residents that is located between the communities of Gordonvale and Yorkeys Knob. These releases followed protocols similar to those of [14], but the release zones were centred within suburban landscapes, providing a continuous habitat for Ae. aegypti. This study investigates the capability of the wMel infection to spread spatially through urban Ae. aegypti populations and the stability of the infection in invaded regions over time.
Spread from localized releases to surrounding uninfected areas depends on mosquito dispersal and relative population densities. Spatial spread can be slowed or stopped if densities are higher in surrounding uninfected areas [20]. Dispersal of Ae. aegypti varies with local environmental conditions. Poor habitats generally induce larger dispersal distances as gravid females must travel further to find the relatively rare oviposition sites [34,35,36]. Despite its global success as an invasive species in tropical habitats, presumably through dispersal of eggs and larvae [37], adult Ae. aegypti are generally considered weak dispersers. Females usually remain within 50-150 m of their eclosion site [34,38,39,40,41,42]. They appear to disperse poorly across highways [31,42,43] and through vegetated parkland [44]. Occasional longrange dispersal, on the order of 0.5-1 km, has been observed [45,46,47,48]. However, given the bistable dynamics of wMel in Ae. aegypti, rare long-range dispersal will not accelerate Wolbachia spread because the infection will not increase locally from low initial frequencies [20,31].
We document local wMel establishment and heterogeneous spatial spread from the 2 relatively large release areas. Our new data demonstrate that local Wolbachia introductions can succeed, persist for at least 2 years, and produce slow spatial spread. Using graphical summaries, we approximate the rate of spatial spread and the width of the spreading wave. We also show that our field data are broadly consistent with simple mathematical models that depend critically on bistable frequency dynamics for wMel transinfected into Ae. aegypti. These models involve only 2 parameters, one describing the position of the unstable threshold point,p, and the other, σ, describing average Ae. aegypti dispersal distance. Both parameters can be estimated independently of spread data [31]. We also present likelihood-based data analyses that fit simple curves to estimate the shape and speed of Wolbachia spread. The shape of the advancing wave is summarized by wave width, defined as the inverse of the maximum slope in infection frequencies, averaged over the wave front [49]. As discussed below, wave width provides an estimate of dispersal distance averaged over time. Wave speed is defined as the average rate of movement of an intermediate infection frequency (e.g., 0.5.) The theory of bistable waves leads to a simple prediction for wave speed in terms of wave width andp, the threshold infection frequency above which local increases in infection frequencies (p) are expected [20,24,31,50]. The observed speed of wMel spread in Cairns is broadly compatible with this prediction, and the estimated wave width is also consistent with independent estimates of dispersal. Moreover, the lack of clear establishment or spread from our third, significantly smaller, release area (only 0.11 km 2 ) is consistent with the prediction for bistable dynamics that releases must be conducted over sufficiently large areas to initiate spatial spread.
Our likelihood analyses also quantify significant heterogeneity in rates of spatial spread that is apparent from our graphical representations. We attempt to link this heterogeneity to easily measured habitat variables. Heterogeneity in host population density is expected to strongly influence Wolbachia invasions subject to bistable dynamics, especially affecting wave speed and potentially restricting the extent of spread [20]. Even if Ae. aegypti disperse equally in all directions, heterogeneities in population density produce asymmetries in net migration. This asymmetry accelerates spread from high-density patches to low-density patches and decelerates-or halts-spread out of low-density patches [20]. Habitat variables such as shade, yard condition, and abundance of oviposition sites have been correlated with Ae. aegypti abundance [41,51]; the frequency of Wolbachia infection within the release zone in Gordonvale, Queensland, was higher in neighbourhoods with more brick and screened houses, which are associated with lower Ae. aegypti abundance [32]. This motivates our attempts to understand patterns of local spread by inferring local densities from easily measured habitat variables. However, the variables we assessed did not predict observed heterogeneities in spread beyond the release zones. aegypti adults infected with the wMel strain of Wolbachia between January 10th and April 24th, 2013. The release zones, located in the suburbs of Edge Hill/Whitfield (EHW), Parramatta Park (PP), and Westcourt (WC), were within 2 km of each other and encompassed 0.97 km 2 , 0.52 km 2 , and 0.11 km 2 , respectively. Mosquitoes were released evenly throughout each release zone at weekly intervals. Total BG-Sentinel trap collections for EHW, PP, and WC are summarised in S1 Table. Our collections continued for about 2 years and are summarized in 4 time intervals. Weekly trap yields at EHW and PP decreased progressively from the onset of each dry season but rose again sharply at the beginning of each wet season (Fig 2). Mosquitoes were caught in consistently higher numbers at PP than at EHW (two-tailed Student t test: P < 0.001), and onsite traps (traps within the release zone) collected mosquitoes at a faster rate than offsite traps (traps outside the release zone; two-tailed Student t test: P < 0.001). When accounting for seasonal changes, yields of uninfected mosquitoes caught in offsite traps at both sites tended to decrease over time (Fig 2 panel A). At PP, there was a corresponding increase in infected mosquito numbers, while at EHW, infected mosquito yields were relatively consistent throughout.

Mosquito collections and abundance
Among onsite traps, yields of infected mosquitoes were consistent with seasonal expectations and were stable over time (Fig 2 panel B). The higher local infection frequencies onsite might lead to an assumption that uninfected mosquito numbers would decline more rapidly than those offsite, but this was not observed, with uninfected mosquito yields at both sites increasing sharply in W2. This proliferation was particularly surprising considering the 2-month-long periods at EHW in the previous season, D2, during which no uninfected mosquitoes were caught onsite (Fig 2 panel B). Wolbachia frequencies: Onsite EHW and PP were both invaded quickly, and by the time releases had finished, Wolbachia infection frequencies within each release zone had reached p = 0.85 (S1 Fig). Following the final releases, p remained relatively stable and near fixation within each release zone. However, in W2, onsite p at EHW dropped from 0.96 to 0.84, the lowest recorded since monitoring began. Considering Fig 2 panel B, it appears that this was due to neither imperfect maternal transmission of Wolbachia [26] nor increased mortality among infected mosquitoes, as their Spread of dengue-suppressing Wolbachia in Aedes aegypti numbers increased to levels similar to those observed in W1. Rather, a sudden influx of uninfected mosquitoes seems most plausible. Averaged across all 4 seasons, 0.88 of mosquitoes within the EHW release zone were infected, while 0.90 were infected at PP.
The WC release zone was invaded as quickly as EHW and PP. However, beginning in September 2013, onsite p dropped sharply to p < 0.7, after which frequencies fluctuated. While onsite p never dropped below any plausible value forp, at no point did the invasion at WC exhibit either the near-fixation values of p or the temporal stability observed at both EHW and PP (see figures below and compare panel C of S2 Fig with panels A and B).
Spread of wMel at EHW and PP but apparent collapse at WC: Graphical analysis The changes of p with time at EHW, PP, and WC between 7 May 2013 and 30 April 2015 are displayed in Figs 3, 4 and 5, respectively, along with trap locations and yields. The plots, based on spatial averaging (ordinary Kriging as described in the Methods section, performed using ArcMap 10.2.2 [52]), show considerable seasonal heterogeneity in the spatial structure of the invasions at EHW, PP, and WC.
At EHW (Fig 3) after D1, the infection was confined largely to the north and northeast, but by the end of W1, the invasion had spread to the east, northeast, and southwest. This pattern persisted through D2, with a small retraction in the north and expansion in the east, though for this season, Kriging was affected by a small sample size (N = 31). Kriging on W2 trap data At PP (Fig 4), spread through D1 was confined mostly to the southeast, from the edge of the release zone up to Mulgrave Road. In the following season (W1), infected mosquitoes were found south across Mulgrave Road and north of the release zone. The infection persisted south of Mulgrave Road but only at below-threshold (p 0.3) frequencies. Over D2, the invasion expanded in range, with high frequencies observed in the north and the southeast and moderate frequencies in the northwest.
At both EHW and PP, the area covered by the infection tended to increase over time (Figs 3 and 4; summarized in panels A and B of S2 Fig), except for PP in W1, in which the area within the p ! 0.8 contour decreased by 3% from D1, and for EHW in D2, in which the area within the p ! 0.8 and p ! 0.5 contours decreased by 7% and 1%, respectively. Nevertheless, from D1 to W2, the area enclosed by the p ! 0.8 contours grew by 85% at EHW and 77% at PP.
At WC (Fig 5; S2 Fig panel C), establishment or spread of the infection was not observed. Following D1, the Wolbachia invasion failed to expand to the south or west of the release zone. Within the release zone, a gradual retreat of the p ! 0.8 contour was observed, with several onsite traps in W2 registering p 0.3. The area covered by the infection at WC reached a peak at W1, but by W2 the area enclosed by the p ! 0.8 and p ! 0.5 contours had decreased by 52% and 44%, respectively, from this maximum (S2 Fig panel C). Wolbachia also failed to spread from WC into PP (or vice versa), but these areas are separated by parkland, which is likely to act as a barrier to movement and prevents ongoing monitoring there.

Spread of wMel at EHW and PP: Heuristic analysis of graphical data summaries
The time interval from D1 to W2 is around 1.5 years, which can be approximated as 15 generations, assuming about 10 generations per year (explained below), or simply viewed as 548 days. When containers are initially colonised and food is available for larvae, developmental time is likely to be rapid at 7-10 days. However, larval populations can rapidly exceed the carrying capacity of the container and its food source (typically leaves), and development is then slowed to 20-50 days [53]; these variable conditions produce a range of adult body sizes that is typically found in field samples from Cairns [54]. If we assume an intermediate value of 20 days in the field, along with time for adult maturation to mating and blood feeding (2-3 days post eclosion), blood-meal digestion and egg formation and oviposition (4 days), and egg embryonation (3 days) [55], this adds another 10 days of adult and egg developmental time. In Cairns, a cooler winter period will lengthen developmental periods, while dry periods delay hatching. Overall, 10 generations per year is likely to be a reasonable estimate.
S2 Fig provides approximations for the areas covered by wMel in different seasons after the releases. For EHW, the area covered in which wMel has at least frequency 0.5 is about 1.3 km 2 in D1, and this rises to about 2.2 km 2 in W2. We can calculate wave speed per generation (assuming 10 generations a year) or per day using alternative geometric approximations described in the Methods section: approximation 4 assumes a circular release area, approximation Eq (5) assumes a rectangle (for which we approximate parameter y = 2 [i.e., a release area twice as long as wide]), or approximation 6 assumes a rectangle in which spread does not occur (or is not monitored) in one direction. (Note that very little spread occurred to the south at EHW.) The resulting estimates of wave speed per day are, respectively, c d = 0.35 m per day, 0.31 m per day, and 0.45 m per day. If we assume 10 generations per year (and so 15 generations separating D1 from W2), the corresponding wave speeds per generation are: c = 12.9, 11.2, and 16.6 m/gen.
Assuming dispersal parameter σ % 100 m/(gen) 1/2 and unstable equilibriump % 0:3 (see Turelli and Barton [31]), the cubic diffusion approximation for wave speed (see Eq 2), c ¼ sð 1 = 2 ÀpÞ, predicts roughly 20 m/gen. As discussed in the context of our likelihood analyses below, the discrepancy between the estimated speeds and this analytical prediction can be resolved by assuming longer generations, a higher unstable point, and/or long-tailed dispersal [31].
For PP, the area covered in which wMel has at least frequency 0.5 is about 0.65 km 2 in D1, and this rises to about 1.17 km 2 in W2. Using our geometric Models 4, 5 and 6 (with y = 2, as for EHW), the resulting estimates of wave speed per day are, respectively, c d = 0.28 m per day, 0.25 m per day, and 0.37 m per day. If we assume 10 generations per year (and so 15 generations separating D1 from W2), the corresponding wave speeds per generation are: c = 10.4, 9.0, and 13.4 m/gen. The speed estimates for PP are systematically smaller than for EHW. As discussed in the Methods section, both wave speed and wave width (describing the distance over which infection frequencies change appreciably) are proportional to average dispersal distances. Thus, slower wave speed is expected if the higher adult densities observed at PP versus EHW translate into a more desirable habitat and consequently smaller average dispersal distances (lower σ). Consistent with this, we find a sharper wave at PP as quantified by smaller average distances between the 0.3 and 0.8 contours at PP than EHW; these distances average 326 m at EHW and only 252 m at PP.

Spread of the infection at EHW and PP: Likelihood analysis
Our likelihood analyses are independent of the graphical summaries produced by Kriging. They rely on an approximate description of the expected shape of local spread and/or collapse (see Eq 7 in the Methods section). We present several successive analyses that summarize the rate and pattern of spatial spread of Wolbachia at EHW and PP. Our summaries focus on 2 statistics: wave width and wave speed. We start by analysing the data averaged over space and time, then present more detailed analyses that document heterogeneous spread. We begin by analysing the data assuming that observed frequencies deviate from deterministic expectations only because of binomial sampling variation. We then use a more complex probability model that accounts for additional sources of heterogeneity. Finally, we explicitly test for directional heterogeneity in rates of spread, as documented visually in Figs 3 and 4. The details of the likelihood analyses are relegated to S1 Text.
Analysis of pooled data. S2 Table provides the maximum likelihood estimates of r 0 and w of the Gaussian/logistic Model 7 for each time interval with the corresponding likelihood Log(L). Note that r 0 ! 0 is measured relative to the centre of the release areas, so that a value near 340 m (220 m) corresponds to the edge of the EHW (PP) releases. When r 0 << w, the Gaussian/logistic model approximates a Gaussian, centred at 0, as illustrated by the figure in our Methods section. For EHW, the increase in r 0 (wave position) becomes roughly linear with time after the first 2 intervals (i.e., after day 120), and the estimated values of w (wave width) settle down to approximate constancy. In this initial phase of spread, r 0~w and the fitted model approximates a logistic. Dropping the first 2 time intervals, the regression of r 0 on time has slope c d = 0.474 m per day. This estimate is broadly consistent with our heuristic graphical wavespeed estimates of 0.31-0.45 m per day. The mean wave width (again, dropping the first 2 intervals) is 439 m.
For PP, the increase in r 0 becomes roughly linear with time after the first 3 intervals (i.e., after day 150), with r 0~w , and the fitted model approximates a logistic. Dropping the first 3 Spread of dengue-suppressing Wolbachia in Aedes aegypti time intervals, the regression of r 0 on time has slope c d = 0.289 m per day at PP; the mean wave width is 366 m. Again, the likelihood estimate of wave speed is comparable to our heuristic graphical estimates of 0.25-0.37 m per day. Moreover, our statistical results agree qualitatively with our graphical analyses in showing slower spatial spread at PP than EHW, with the slower speed accompanied by a sharper cline in frequencies (i.e., smaller w).
Before comparing these results to our theoretical predictions, we consider 2 statistical refinements, whose methods are described in S1 Text. The likelihood analysis summarized in Fig 6 assumes that binomial sampling is the only source of variation in infection frequencies observed at fixed distances from the release areas. It also fits frequency data averaged over time and space. In S1 Text section 1.1, we first generalize our statistical model to allow for nonbinomial variation, implicitly accounting for factors such as environmental heterogeneity that contribute to differences in infection frequencies among sampling locations equidistant from the release areas. The additional stochasticity is summarized by a parameter F, where F > 0 accounts for greater-than-binomial variation. The likelihood analyses for EHW and PP incorporating this factor appear in S3 and S4 Tables, respectively. Our second refinement fits the data directly without spatial or temporal averaging. This involves estimating a constant, denoted R, which describes how far the wave has moved from the centre of each release area when approximately linear spread was observed (see S1 Text section 1.2 for details). Table 1 shows that our likelihood estimates of wave speed, c d (measured as meters per day), and wave width, w, are relatively insensitive to these more refined analyses.
Analyses of heterogeneous spatial spread at EHW. Fig 3 shows the heterogeneous spread of wMel at EHW, with more rapid spread to the north and east than to the south and west. Here we summarize a likelihood analysis that quantifies the heterogeneity. The details of the analysis are provided in S1 Text section 1.3. For simplicity, we divided the samples into 4 equal-angle triangular sectors, centred on the middle of the release area. We fit Model 7 to the pooled data from each sector separately and sought the orientation of the sectors that produced the best fit to the data (allowing for additional nonbinomial variance as discussed above). This allows us to estimate 4 separate wave speeds and apply a likelihood test for heterogeneity. The results of the analysis, demonstrating statistically significant spatial heterogeneity in the rates of spread, are summarized in S5    As demonstrated by Turelli and Barton [31], even with fast local dynamics and long-tailed dispersal, we can accurately approximate average local dispersal as σ = w/4 m/(gen) 1/2 (Eq 2). From this we infer σ % 115 m/(gen) 1/2 at EHW; in contrast, we obtain σ % 95 m/(gen) 1/2 at PP. Given that the support intervals for the estimates of w at EHW and PP do not overlap (Table 1), we expect these results reflect differences in local dispersal. Given that PP has consistently higher population densities, this difference may reflect less dispersal in a habitat where mosquito densities are higher. However, this needs further testing against alternative hypotheses, such as more dispersal barriers surrounding the PP versus the EHW release areas. It is notable that both the EHW and PP estimates of dispersal are consistent with values obtained from release-recapture experiments (reviewed in [31]). If we assume that the wave speed follows the cubic diffusion approximation c ¼ sð 1 = 2 ÀpÞ, per generation and that generations are T days long, we can in principle reconcile observed wave speeds with expected wave speeds at each release site by choosingp and T appropriately, namely

Comparison of observed spread at EHW and PP to theoretical predictions
where σ is the local dispersal estimate and c d is the observed wave speed per day. For instance, if we assume that at both EHW and PP,p ¼ 0:3, the observed and expected wave speeds can  be reconciled if we assume that T = 46 days for EHW, whereas T = 63.3 days for PP. Given that population densities are higher for PP, increased crowding may indeed produce longer generation times [53]. These times are systematically larger than our conjecture of 10 generations per year, which we supported by an informal data review above. These inferences assume that the cubic-diffusion prediction for wave speed (c ¼ sð 1 = 2 ÀpÞ per generation) is accurate for these field populations. However, as shown by Turelli and Barton [31], long-tailed dispersal with fast local frequency dynamics (as expected with complete cytoplasmic incompatibility, corresponding to s h = 1 in the models of [31]), can slow the expected wave speed by 20%-40% below the cubic-diffusion prediction. If the expected wave speed is reduced by 30%, the observed wave speeds match the modified expectations with generation times reduced to 32.2 and 44.3 days at EHW and PP, respectively. These times are closer to our conjecture of 10 generations per year. In general, there seems to be reasonable quantitative agreement between the slow observed wave speeds and the predictions of simple models using parameter values that are consistent with the poorly known field biology of Ae. aegypti and the deleterious fitness effects of wMel in Ae. aegypti. Despite many caveats, including uncertainty about parameter values and the imprecise meaning of the one-dimensional unstable pointp for populations with overlapping generations and complex ecology [27], the observed spread rates at EHW and PP are clearly consistent with approximation Eq (1) using plausible estimates of dispersal distance, the unstable point, and generation time.

Comparison of apparent collapse at WC to theoretical predictions
In contrast to EHW and PP, the releases at WC did not lead to clear establishment and certainly did not produce spatial spread (see Figs 5 and 7). Turelli and Barton [31] provide conditions on minimum release areas (and maximum dispersal distances) consistent with spatial spread, allowing for long-tailed dispersal and rapid local dynamics. We expect that p % 0:25 À 0:3 and σ % 100 m/gen 1/2 . If these parameter estimates are accurate, the release area at WC is likely to be just below the minimum needed to produce successful local establishment and spread (see Table 2 of [31]). Moreover, the fact that the apparent collapse at WC is extremely slow is consistent with the slow dynamics expected near that critical size threshold for wave-establishing releases [31]. Overall, the bistable dynamics of wMel in Ae. aegypti will impose some minimum release size, and only WC is near a plausible minimum. To rigorously test the minimum-release-area predictions of Barton and Turelli [20] and Turelli and Barton [31], several more replicate releases in small areas would be needed.

Discussion
Our data demonstrate that wMel can be stably established locally within urban areas surrounded by uninvaded but suitable habitat. Hence, stable population replacement is not limited to small isolated habitats such as those where the initial releases and establishment of wMel in Ae. aegypti took place (cf. [14]). Moreover, the temporal increase in infection frequency within the EHW and PP release zones was comparable to that seen in the isolated areas. In contrast, the smallest release area, WC, did not show stable invasion. This suggests that there is little impediment to the local establishment of Wolbachia in urban areas, provided the releases are conducted over sufficiently large areas (e.g., on the order of 0.5 km 2 when dispersal distances are comparable to those in Cairns [31]). These findings highlight the feasibility of patchy releases across large cities, suggesting that area-wide replacement can be produced gradually, with patchy releases complemented by natural local spread. At EHW and PP, the area in which Wolbachia persists at high frequency roughly doubled after 2 years (Figs 3  and 4).
The failure of wMel to establish and spread at WC seems attributable to the small area of the release zone, as the habitat conditions in and around WC are similar to EHW and PP. This is consistent with mathematical predictions concerning the minimum release zone radius, R crit [20,31]. Based on the wider advancing wave front seen at EHW versus PP, we infer greater average dispersal distance at EHW (which is likely to provide fewer feeding and breeding opportunities than PP). Mosquito dispersal differences probably explain the faster spread observed at EHW versus PP. In contrast, the slow temporal and spatial dynamics of local infection frequency at WC suggests that 0.11 km 2 , the area of the WC release zone, may be very close to the minimum size needed to initiate spread, at least for the levels of dispersal typical of Cairns. When contrasted against the successful spread at PP, we conclude that the critical release area under Cairns conditions is somewhere between 0.11 km 2 and 0.52 km 2 . In tropical regions that support denser Ae. aegypti populations, we expect lower dispersal distances. This would allow successful local establishment using smaller release areas, but spatial spread would also be expected to be even slower than the 100-200 m per year observed at EHW and PP.
The heterogeneity in both the speed and patterns of the spatial dynamics at EHW and PP suggests that local environmental factors greatly influence the spread of Wolbachia transinfections (such as wMel in Ae. aegypti) that produce significant fitness costs. Spread at each site exhibited strong spatial structure throughout the study, and the structure persisted across the monitoring period. Areas that were easily invaded during the first dry season after the releases (D1, see Figs 3 and 4) generally stayed invaded in successive seasons, and the autocorrelation among mosquito numbers and infection frequencies increased as the study progressed (S7 Table). The invasion spread well beyond the initial release zones at EHW and PP, and our likelihood analyses (Fig 6) suggest that slow but steady spread would continue in the absence of further releases until significant barriers to dispersal are encountered.
Barriers to spread can include both barriers to Ae. aegypti dispersal and variation in Ae. aegypti population density [20]. At PP, the invasion spread south from the release zone immediately but never established to a high frequency south of Mulgrave Road. Nevertheless, infected mosquitoes were caught at low frequencies south of Mulgrave Road from season W1 onwards. These observations are consistent with the demonstration in Trinidad that roads represent partial barriers to Ae. aegypti dispersal [43]. At the very least, such barriers slow wave propagation [20]. It remains unclear whether Mulgrave Road provides a sufficient barrier to stop the wave of Wolbachia, as is the case of the Bruce Highway at Gordonvale. There, Wolbachia have failed to invade an area adjacent to the 2011 release zone for several years, despite persistent migration across the highway [31]. Other evidence from mark-release experiments and genetic studies have pointed to potential barriers (roads, rivers, forests) to movement of Ae. aegypti at a local scale [43].
In W2 at EHW, there was an apparent drop in p in the southern half of the release zone. This was unexpected given that in previous seasons traps in this region had recorded Wolbachia frequencies close to fixation. It appears that the drop was due to a sudden increase in uninfected mosquito numbers onsite, which may represent the hatching of dormant uninfected eggs or an early influx of uninfected mosquitoes from an external source at the start of W2. One possibility is that Wolbachia infected larvae experienced a fitness cost under highstress conditions prevailing at that time; such costs have been recently documented under stressful conditions that produce a range of adult sizes [19,56] similar to those seen under field conditions [54], even though earlier studies suggested only modest fitness costs associated with wMel [11,16]. The openness of the EHW study area may also make it more vulnerable to reinvasion, as immigration into the release zone was possible from 360˚of the surrounding area. In contrast, one of the long edges of the PP release zone was bounded by parkland that blocked immigration and may help explain the slow but uninterrupted spread observed there (S2 Fig). Seasonal variations in invasion dynamics are expected when mosquito abundance varies throughout the year. Models of Wolbachia population dynamics show that when mosquito host abundance fluctuates seasonally, p is expected to decline in the low seasons because of slow recruitment and comparatively high mortality among infected imagoes [57]. This was not observed in PP but was in EHW, where during D2 the area covered by the infection shrank and offsite p dropped considerably, only to recover the following wet season. PP may have been shielded from these effects by its apparent abundance of good mosquito habitat, reflected in its high trap yields throughout the study. Very few mosquitoes were caught in EHW during D2, and the sluggish recruitment there was a likely cause of the retraction of the invasion during that season.
Initial onsite infection frequencies at EHW correlated positively with window screens and negatively with habitat quality (S8 Table). This corroborates the findings of Hoffmann et al. [32] that following mass release of infected mosquitoes within an area, Wolbachia frequencies are highest in areas of poor mosquito habitat. However, no relationship was found between trap yields and simple measures of habitat quality (S2 Text). While BG-Sentinel traps can pick up on seasonal changes in mosquito abundance [58,59], they may not be able to give precise estimates of local Ae. aegypti densities at the scale of deployment used in this study.
Modelling offsite spread as a function of easily observed habitat variables was inconclusive (S2 Text). No variables were consistently predictive across seasons or sites, and in some cases variables that were expected to encourage spread (i.e., areas of low Ae. aegypti density: those with window screens, low-set dwellings, poor habitat quality) were found to deter it. However, predictions based on variation in population density and uniform dispersal (e.g., [20]) may be confounded by active searching for favourable oviposition sites [60], if density variation is driven by local habitat quality.
The lack of any discernible predictor variables, the strongly heterogeneous spread, and the drop in infection frequencies at the centre of the EHW release zone during W2 suggest that stochastic processes may have played a role in the invasions of EHW and PP. This is surprising considering that the scale of the Cairns invasions was much larger than those thought to be susceptible to stochastic effects associated with very small numbers of infected individuals [61]. Alternatively, a series of highly localised processes may have influenced the heterogeneity of the spread. If this is the case, our BG-Sentinel traps may be too dispersed to pick up on local variability that could inform future releases [62]. Spatial structure in Ae. aegypti populations has been observed at the house scale in Cairns [63]; in ensuing releases of Wolbachia-infected Ae. aegypti, a more clustered placement of traps within and around the advancing wavefront may provide a clearer picture of the processes at work. Despite the heterogeneity, our simple 2-parameter model seems to plausibly account for the slow rates of spread observed at our larger release sites. Bistability for the wMel transinfection, versus the apparent tendency for successful natural Wolbachia infections to spread even when very rare, accounts for the fact that spread in Cairns is orders of magnitude slower than observed Wolbachia spread in natural Drosophila populations [30,31].
In summary, we have found rapid local establishment of wMel Wolbachia in the Ae. aegypti populations of urban release areas, with an adjacent suitable habitat available for mosquito dispersal. In the 2 release areas that exceeded the predicted minimum size threshold for local establishment, the infection remained at a moderately high frequency for 2 years. Moreover, wMel spread slowly outward at a rate consistent with theoretical predictions, based on realistic estimates of local dispersal and the position of the unstable equilibrium frequency. While this rate of Wolbachia spread is extremely slow, these findings indicate that large urban areas can be transformed gradually with patchy local Wolbachia releases [31]. Local information about barriers to dispersal can inform the minimum number of releases required, but it remains a challenge to understand the heterogeneity of spatial spread in terms of easily obtained data concerning habitat quality.

Study area
Ae. aegypti adults infected with the wMel strain of Wolbachia were released by Eliminate Dengue staff at 3 zones in Cairns, Queensland, from January to April 2013. Overall, 131,420 mosquitoes were released at EHW, 286,379 at PP, and 35,196 at WC. More mosquitoes were released at PP because of its denser local Ae. aegypti population, and the 3 sites began with comparable proportions of wild and introduced mosquitoes.
The largest of the zones, EHW, was also the most open, situated amid residential suburban development with no potential dispersal barriers in its vicinity. PP and WC were in contrast both semiclosed, with each having 1 side of the release zone bounded by parkland so that from the centrum only 268˚of the release zone at PP and 254˚at WC was connected directly to urbanised Ae. aegypti habitat. Additionally, both PP and WC were near major roads, specifically Mulgrave Road to the southeast and Captain Cook Highway to the northeast (Fig 1), which could act as barriers to mosquito dispersal. These roads each consisted of 6-10 lanes totalling >50 m throughout. They were flanked by commercial buildings, for the most part, interspersed with apartment complexes.
Due in part to an abundance of modern single-storey housing, EHW was known to support a lower density of Ae. aegypti than the other sites. In contrast, PP was adjacent to Cairns' central business district, and its household size of 2.00 per dwelling-smaller than that of either EHW (2.29) or WC (2.11) (http://profile.id.com.au/cairns/population)-reflects a larger number of multistorey apartment complexes, fewer bungalows and a higher density of unscreened older houses. At each of the 3 locations, the area enclosed by the release zone was thought to support a higher density of Ae. aegypti than the area surrounding the release zone, where there tended to be a higher density of modern houses. Cairns experiences a tropical monsoon climate, with a wet season running from November to April.
The successful establishment and spatial spread of Wolbachia following releases requires that the release area exceed a theoretical minimum, described by a critical radius R crit [20,31]. While both EHW and PP clearly surpass the minimum, the small area of WC makes establishment there uncertain. Numerical analyses show that R crit depends on both the shape of the dispersal function and the average dispersal distance, with small releases more likely to be successful with lower dispersal that is highly leptokurtic (i.e., showing both more long-distance and more short-distance dispersal that expected under a Gaussian function [31]). Dynamics of establishment and spread likewise depend on the size of the release zone relative to R crit , with faster dynamics predicted when the release zone is very large or very small relative to the minimum size. From this, we expect EHW and PP to display relatively rapid spread. In contrast, given that the release area at WC is close to the minimum, the failure or success of establishment is expected to take on the order of 2 years (about 20 generations) to ascertain [31].

BG-Sentinel trap collections
BG-Sentinel mosquito traps (Biogents AG, Regensburg, Germany) were set up within and around each release zone at fixed positions in the yards of consenting householders, covering a distance of 25-530 m (EHW), 35-670 m (PP), and 30-520 m (WC) from the release zones in every available direction. The exact number of trapping sites varied over time as traps were moved from households whose residents had either moved or decided to terminate participation in the project. New traps were also added occasionally. By April 2015, data had been collected from 182 traps at EHW (44 onsite traps within the release zone, 138 offsite traps outside the release zone), 142 traps at PP (42 onsite, 100 offsite), and 74 traps at WC (20 onsite, 54 offsite). BG-Sentinel traps catch mosquito adults by means of a visual lure (black entry cup) and suction fan. In Cairns, they capture Ae. aegypti with high specificity and in large numbers [59].
Traps were checked weekly from 7 May 2013 to 30 April 2015. Traps that failed because of malfunction, invasion by predators (ants, spiders) or physical disturbance were scored as null observations for that time point. Adult mosquitoes from each trap were stored in ethanol at -20˚C. Traps at WC ceased being checked after 1 April 2015, as new releases began in the area.
Ae. aegypti abundances are known to vary considerably throughout the year in Cairns [64]. This led us to partition the trap data into seasonal units, reflecting the 6-month wet and dry seasons of Northern Queensland. The first dry season D1 (May 2013-October 2013) began immediately after the releases, followed by the first wet season W1 (November 2013-April 2014), the second dry season D2 (May 2014-Oct 2014), and the second wet season W2 (Nov 2014-Apr 2015). This allowed comparisons across time as well as space.
For our graphical analyses, we aggregated data for each trap to give the following: (1) total mosquito abundance per season, (2) total number infected with Wolbachia per season, (3) an average number of mosquitoes observed per week (total and infected), and (4) a seasonal infection frequency, p. As in Hoffmann et al. [18], we checked species identity and Wolbachia infection status by PCR. For each mosquito, PCR was performed using 3 primer sets, Aedes universal primers (mRpS6_F/mRpS6_R), Ae. aegypti-specific primers (aRpS6_F/aRpS6_R), and Wolbachia-specific primers (w1_F/w1_R).

Analyses of spatial spread: Mathematical background
Mathematical background. Before describing the graphical and statistical analyses, we review key theoretical results summarized in Barton and Turelli [20] and Turelli and Barton [31] concerning spatial spread with bistable dynamics. Letp denote the threshold frequency above which the wMel infection frequency is expected to locally increase. A realistic model of wMel frequency dynamics must minimally include overlapping generations with age structure [27] and density dependence [53,65]. Hence,p is best viewed as a convenient summary statistic indicating the relative fitness costs associated with wMel in comparison to the frequency-dependent advantage produced by CI [27]. Turelli and Barton [31] summarize empirical evidence suggesting that near Cairns,p for wMel in Ae. aegypti is likely to be in the range 0.25-0.35. Dispersal behaviour is summarized by the parameter σ, the standard deviation of dispersal distances per generation along any given axis. Assuming Gaussian dispersal, the mean Euclidean distance between the birthplaces of mothers and daughters is s ffiffiffiffiffiffiffiffi p=2 p % 1:25s. As reviewed by Barton and Turelli [20], a partial differential equation model for spatialtemporal dynamics, with an idealized cubic approximation for local infection-frequency dynamics, provides analytical predictions for the rate of spatial spread and the shape of the spreading wave [66]. Assuming complete CI (i.e., all embryos produced from incompatible crosses die), the predicted wave speed is per generation. The spreading wave assumes a characteristic asymptotic shape. If we define wave width, w, as the inverse of the maximum slope of infection frequency [49], the cubic-diffusion approximation produces Prediction 3 provides an estimate of average dispersal distance from spatial infection frequency data once steady spatial spread is initiated. Using the asymptotic wave formula that generates Eqs 2 and 3, we expect infection frequencies to change from about 0.18 to 0.82 over 3σ. (Similarly, infection frequencies are expected to change from 0.3 to 0.8 over about 2.23σ.) Using estimates of c and σ with Eq 1, we can approximatep from joint estimates of speed and width--assuming these analytical predictions are robust.
Turelli and Barton [31] examined the robustness of predictions 2 and 3 to both long-tailed dispersal and the rapid local changes in infection frequency expected with complete CI. Relation 2 between wave width and σ is quite robust, with maximum departures on the order of 10% for 0:2 p 0:35. In contrast, over the same range of parameters and models, longtailed dispersal (corresponding to higher frequencies of both long-distance and short-distance dispersal than expected under a Gaussian) with complete CI can reduce wave speed by 10%-40%. Hence, estimates of σ from observed wave width and prediction 3 are likely to be quite robust, whereas wave speed may be systematically overestimated by prediction 2.
Barton and Turelli [20] presented conditions for wave initiation and wave stopping. To avoid being swamped by immigration of uninfected individuals, releases must cover a sufficiently large area to initiate an expanding wave. Assuming thatp 0:35, releases within circles of a radius greater than 3σ should suffice to initiate spatial spread. As shown by Turelli and Barton [31], even smaller releases should initiate spread with long-tailed dispersal, especially if Wolbachia-induced fitness reductions mainly involve fecundity. However, bistable waves can be relatively easily stopped by environmental heterogeneities and barriers to dispersal such as roads. This phenomenon is illustrated by the wMel frequency data from Pyramid Estate, a small suburb separated by a highway from the 2011 release site in Gordonvale, about 20 km south of central Cairns. As reported by Hoffmann et al. [18], the wMel infection frequency has remained stable at over 95% in Gordonvale since 2011. Occasional wMelinfected Ae. aegypti migrate across the highway to Pyramid Estate. Yet, the infection frequency has never increased appreciably in Pyramid Estate and the long-term average is only 0.106. As noted by Turelli and Barton [31], this provides an approximate lower bound ofp ! 0:21.
Bistable waves can be slowed or stopped by increases in population density. In general, increases in population density slow wave speed. However, asp approaches 0.5, even very small increases suffice to stop wave spread [20]. Hence, natural heterogeneity in Ae. aegypti population density is likely to produce heterogeneous rates of spatial spread.

Analyses of spatial spread: Graphical analyses of infection-frequency data
Graphical analyses of infection-frequency data. Ordinary Kriging [67] was performed to interpolate data on a map and visualize the patterns of spread based on the seasonal p of traps from which at least 4 mosquitoes had been collected in the season. Interpolative maps were created predicting Wolbachia frequencies throughout each site, from which the direction and extent of the invasion could be inferred. Kriging was performed in ArcMap 10.2.2 [52] using an exponential semivariogram model and a 24-point nearest-neighbour search function for EHW and PP and a 16-point nearest-neighbour search function for WC (because of fewer traps at WC).
Kriging maps of infection frequencies averaged over 6-month windows allowed for the rough approximation of c and σ. To approximate c, we first calculated the areas enclosed by the p ! 0.5 and p ! 0.8 map contours for each season, using ArcMap 10.2.2. These contours defined areas within which Wolbachia frequency was greater than the highest possible estimates ofp consistent with spatial spread (roughlyp ! 0:5, [20]) and areas within which Wolbachia was near fixation. With each plot representing 6 months of spread, the increase in area can be used to estimate wave speed using formulas provided below. To estimate σ, we calculated for each season the average distance between the p ! 0.3 and p ! 0.8 map contours. This was achieved by plotting 36 lines at 10˚intervals from the centre of each of the EHW and PP release zones. From where each of these lines intersected the p ! 0.8 contour, we calculated the shortest distance to the p ! 0.3 contour.
To translate the estimates of expanding Wolbachia-infected areas into heuristic approximations of wave speed, we used simple geometric models. Consistent with our likelihood analyses below, we assume that 6 months after the releases began, the infection approached its asymptotic rate of spatial spread, c. The relationship between changes in area and speed of the wave front, c, depends on the shape of the release area and whether spread can occur in all directions. The closed curve of fixed length that encloses the largest area is a circle, hence we can approximate the maximum rate of uniform spread in all directions from a release by assuming a circular release area with symmetric spread in all directions. If the area in which the infection has frequency ! 0.5 increases from A 0 to A t over t generations, the wave speed per generation is To examine the effect of release-area shape, we can instead assume that the release occurs in a rectangle, with the long sides y times longer than the short (so y = 1 with a square release area). As the infection spreads, the area covered has quarter-circle corners, asymptotically approaching a circle through time. With a rectangular release, the relationship between area change and wave speed is Finally, suppose that the infection can spread in only 3 directions, with no expansion possible in the direction of one of the longer sides (so that the expanding wave asymptotically becomes a semicircle). Then wave speed can be approximated as We apply these approximations to the EHW and PP data below and compare them to model-based likelihood estimates. Note that c and t can be measured in generations or in days.
Likelihood analyses. To estimate wave speed and wave width, we reduce the two-dimensional data to one dimension by calculating distances from the edge of the release area and averaging over time and space.
Defining from the nearest edge. The effective radial distance was taken to be r = r Ã − r min , so that the sample point within the release area farthest from the release area edge has r = 0.
Statistical model for estimating wave speed and width: Instead of fitting a mechanistic model of temporal and spatial dynamics, we fit a statistical model that approximates both initial conditions after our releases and the spreading wave. The model assumes the spreading wave shape approximates the one-dimensional asymptotic solution of the cubic diffusion approximation (see Eq. 13 of [20]). For each time, we approximate the spatial distribution of infection frequencies as a function of distance from the release area, r, using where a ¼ pwr 0 2 Exp 2r 0 w and n ¼ wr In p(r), r 0 indicates the position of the wave when r 0 > w. For small r 0 (r 0 << w), Model 7 is approximated by a Gaussian distribution scaled to total mass α, i.e., where α and v are given by Eq (7b). For large r 0 , Model 7 approaches a logistic cline with width w centred at r 0 , i.e., where w is cline width, as defined in Eq 2. Fig 8 shows p(r) from Model 7 for various r 0 with w = 1. Pooling data. As a first approximation, data were pooled into 9 time intervals with boundaries at 90, 110, 120, 150, 250, 300, 400, 550, 700, and 900 days and distances pooled at 100-m intervals, counting relative to the edge of the release area (r Ã = 0). Maximum likelihood estimates of the parameters r 0 and w of the Gaussian/logistic Model (7) were obtained separately for each time interval, initially assuming that binomial sampling is the only source of variation in frequencies at each distance from the release boundary. Other sources of variation are considered below.
To estimate wave speed, we calculated the slope of the regression of r 0 against time once rough linearity was achieved. Similarly, the wave width was estimated as the average value of w once spread had become roughly linear (see Results).
In S1 Text, we describe 3 additional likelihood analyses. The first provides a simple model to account for variance in infection frequencies above that expected from binomial sampling alone. The second describes an analysis of wave speed and width without pooling the data over time or space. The third is an explicit analysis of the spatially heterogeneous spread at EHW, which supports the heterogeneity apparent from  show the pooled data from EHW and PP; panels B and D show the maximum likelihood fits of those data to the model described by Eq (7). The x-axis in each panel is distance (in meters) from the edge of the release area, the y-axis is infection frequency. Each of the nine colored lines shows the spatially spreading infection, the first centred on t = day 100 (releases began on day 99, blue), the last centered on day 800, dark red), nearly two years after the releases were  Table. Likelihood analyses of wave position and width for EHW and PP using a binomial model. Likelihood analyses of infection frequencies using data pooled over time and by distance from release sites. These analyses account for only binomial sampling variance. For each time interval, the pooled data are used to estimate the parameters r 0 , describing the position of the wave, and w, the wave width, in Eq (7). (XLSX) S3 Table. Likelihood analyses of wave position and width for EHW allowing for non-binomial sources of variation. EHW estimates of the parameters in Model (7) obtained with model (S1) to account for non-binomial sources of variation in infection frequencies via the parameter F. For each time interval, we provide Log(L) and the MLEs for r 0 and w with interval-specific F (left), or with a common F = 0.226 (right). (XLSX) S4 Table. Likelihood analyses of wave position and width for PP allowing for non-binomial sources of variation. PP estimates of the parameters in Model (7) obtained with model (S1) to account for non-binomial sources of variation in infection frequencies via the parameter F. For each time interval, we provide Log(L) and the MLEs for r 0 and w with interval-specific F (left), or with a common F = 0.167 (right). (XLSX) S5 Table. Heterogeneity of spatial spread at EHW. Estimates of wave speed per day (c d ) and wave width (w) in four directions. The four rows correspond roughly to south, east, north and west in Fig 3. (XLSX) S6 Table. Likelihood analyses of wave position and width for WC allowing for non-binomial sources of variation. WC estimates of the parameters in Model (7) obtained with model (S1) to account for non-binomial sources of variation in infection frequencies via the parameter F. For each time interval, we provide Log(L) and the MLEs for r 0 and w with interval-specific F (left), or with a common F = 0.1 (right).