Rapid Sequential Spread of Two Wolbachia Variants in Drosophila simulans

The maternally inherited intracellular bacteria Wolbachia can manipulate host reproduction in various ways that foster frequency increases within and among host populations. Manipulations involving cytoplasmic incompatibility (CI), where matings between infected males and uninfected females produce non-viable embryos, are common in arthropods and produce a reproductive advantage for infected females. CI was associated with the spread of Wolbachia variant wRi in Californian populations of Drosophila simulans, which was interpreted as a bistable wave, in which local infection frequencies tend to increase only once the infection becomes sufficiently common to offset imperfect maternal transmission and infection costs. However, maternally inherited Wolbachia are expected to evolve towards mutualism, and they are known to increase host fitness by protecting against infectious microbes or increasing fecundity. We describe the sequential spread over approximately 20 years in natural populations of D. simulans on the east coast of Australia of two Wolbachia variants (wAu and wRi), only one of which causes significant CI, with wRi displacing wAu since 2004. Wolbachia and mtDNA frequency data and analyses suggest that these dynamics, as well as the earlier spread in California, are best understood as Fisherian waves of favourable variants, in which local spread tends to occur from arbitrarily low frequencies. We discuss implications for Wolbachia-host dynamics and coevolution and for applications of Wolbachia to disease control.


Introduction
The vertically-transmitted intracellular bacterium Wolbachia may be the most widespread [1][2][3] and evolutionarily significant endosymbiont [4,5] of insects and other arthropods. Wolbachia induce many reproductive manipulations within hosts that increase their chance of spreading through females. Their ability to suppress other microbes in their hosts provides a novel method to control human vector-borne diseases such as dengue [6,7] and malaria [8,9]. Yet despite their ubiquity and potential importance in vector-borne disease control, there are few documented examples of Wolbachia infections spreading in natural host populations [10,11].
The most commonly documented host reproductive manipulation by Wolbachia is cytoplasmic incompatibility (CI) [12], which reduces hatch rates for embryos produced when sperm from an infected male fertilizes uninfected ova or ova that carry an incompatible Wolbachia strain. CI provides a reproductive advantage to infected female hosts (whose infected eggs are protected from CI), and it can drive the spread of a Wolbachia infection within and among host populations [11,[13][14][15]. Several factors can affect the spread of CI-inducing Wolbachia infections, including the strength of incompatibility (quantified by H, the relative hatch rate of embryos from incompatible fertilizations), the maternal transmission frequency (12m, where m is the frequency of uninfected ova produced by an infected female), and the fecundity of infected females relative to uninfected females (F, which can also approximate viability effects) [12].
Mathematical analyses [14][15][16] show that irrespective of the level of CI, as measured by H (0#H#1), if a CI-inducing infection satisfies F(12m),1, it will tend to decrease in frequency when very rare. The reason is that the CI-induced reproductive advantage for infected females depends on the frequency of infected males, whereas disadvantages attributable to fitness costs, F,1, or imperfect maternal transmission, m.0, are frequency independent. If the level of CI is sufficient that F(12m).H (i.e., when the CI-inducing infection is very common, more infected offspring are produced by infected mothers than viable uninfected offspring are produced by uninfected mothers), the infection will tend to spread locally once it becomes sufficiently common. Hence for Wolbachia strains that satisfy H,F(12m),1, we expect ''bistable'' dynamics [13][14][15], where there is an unstable equilibrium infection frequency, denoted p u , that satisfies 0,p u ,1; p u separates two stable equilibria, one at 0, the other a high frequency denoted p s (which always exceeds K, [17]). If maternal transmission is imperfect (m.0), then p s ,1 because uninfected individuals are continuously introduced. The unstable equilibrium frequency, p u , plays a central role in both local dynamics and spatial spread. Local dynamics depend on the initial infection frequency, denoted p 0 . If p 0 .p u , the infection frequency tends to increase towards p s ; conversely if p 0 ,p u , the frequency tends to decrease to 0. Spatial dynamics are more complex, and predictions depend on the initial frequencies over an extended area. Nevertheless spatial spread from a localized introduction cannot occur if p u is too large; roughly, bistable infections do not spread unless p u ,K (the exact condition is described in [13] and [18]). If a CI-inducing infection with p u ,K becomes established in a sufficiently large region (see Fig. 3 of [13]), it will tend to spread spatially at a rate determined by average dispersal distance of females and the intensity of CI. These bistable waves can be stopped by regional variation in population density or dispersal barriers [19].
In contrast to bistable infections, Wolbachia strains that provide a frequency-independent fitness advantage to infected hosts, such that F(12m).1, would increase locally, even from low initial frequencies and regardless of whether they cause CI. Local spread would be followed by a ''Fisherian'' spatial wave [20,21] that, unlike a bistable wave, is unlikely to be halted. For such infections, 0 is the unstable equilibrium and the unique stable equilibrium satisfies p s ,1 if maternal transmission is imperfect [12]. As in [13], we use ''Fisherian'' to describe the spatial dynamics of variants that tend to increase even when locally very rare, unlike bistable variants that tend to increase only once they exceed a critical frequency p u .0. The mechanism that reduces the unstable equilibrium to zero may involve interactions between Wolbachia, D. simulans and various pathogenic microbes, as postulated by Fenton et al. [22]. Alternatively it may involve nutritional provisioning [23,24] or other effects as yet unknown. In the absence of relevant field data, we approximate the effects of Wolbachia on D. simulans by a frequency-independent advantage that produces F(12m).1. Nothing about our analyses or conclusions, which rest on the rapid sequential spread of two Wolbachia variants, requires a more detailed model. The key feature of Fisherian variants is that they spread much more rapidly than bistable variants, because small advance propagules can catalyze invasions of new areas [25,Ch. 5], in contrast to the steadily moving wave fronts expected with bistable dynamics. Unfortunately, the temporal and spatial resolutions of our data are insufficient for fitting mechanistic models of spatial spread.
Turelli & Hoffmann [11,17] first documented a Wolbachia infection spreading within and among natural populations. A strain of Wolbachia (wRi) was initially found infecting Drosophila simulans populations in southern California. Turelli and Hoffmann monitored the northward spread of wRi over ten years (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994). Because they repeatedly found both imperfect maternal transmission in nature and reduced relative fecundity for infected females in the laboratory [15,17], they assumed F(12m),1 and inferred bistable dynamics (i.e., H,F(12m),1). Using field-based estimates of m, F and H, their mathematical analyses implied p s <0.94, in close agreement with relatively constant infection frequencies observed in several natural populations for over 20 years [26]. However, as emphasized by Jaenike [27] and elaborated by Carrington et al. [26], the CI-based prediction for p s is quite insensitive to variation in the relative fecundity of infected females, F. Based on field estimates of m near 0.045, F(12m).1 requires only that infected females have a fecundity advantage on the order of 5%. We present new mtDNA data from recent California samples of D. simulans and from a 1961 southern California collection which suggest that wRi invaded California less than 25 years before Hoffmann et al. [28] found it. This supports our new interpretation of Fisherian versus bistable spread in both California and Australia.
Other studies of Wolbachia in natural host populations have also assumed bistable dynamics [10,29], in one case [29] based on observing very rare imperfect maternal transmission (m<0.014), in the other [10] based on an indication of spatial spread analogous to that seen in D. simulans. A definitive example of bistable dynamics comes from recent field releases of Wolbachia-infected Aedes aegypti mosquitoes. The wMel Wolbachia infection was transferred in the laboratory from its native host D. melanogaster to Ae. aegypti as part of a novel strategy for blocking transmission of the human dengue virus [6,7]. The observed rate of infection frequency increase within populations was consistent with significant fitness costs, and hence bistable dynamics; and lowfrequency introductions into neighbouring populations have not led to Wolbachia establishment outside the release areas, as predicted with bistability (Eliminate Dengue Team, unpubl. data).
These mosquito data provide experimental support for bistable Wolbachia dynamics in nature. However, other infections including wMel in D. melanogaster cause minimal CI and persist despite incomplete maternal transmission. These infections may show Fisherian dynamics [12,30,31], consistent with the apparent global spread of alternative forms of wMel, none of which causes appreciable CI [32].
Here we present Wolbachia and mtDNA data that document the sequential recent spread of two Wolbachia infections in eastern Australian D. simulans populations. Previous research [33] indicated that D. simulans were polymorphic for a novel Wolbachia infection (wAu) present at a low frequency that induced no detectable CI. The persistence of this infection is most easily understood if it increases fitness consistent with F(12m).1. Although various potential fitness advantages have been associated with Wolbachia infections in the laboratory, including virus protection [34,35] and fecundity increases [36,37], they have not been demonstrated in nature for wAu or any other Wolbachia strain [38].
We show that the temporal and spatial data for wAu, and for wRi that has recently invaded Australia [39], are best explained by postulating that each variant increases host fitness in nature. We also present a new observation relevant to the history of wRi spread in Californian populations of D. simulans. These new data suggest that wRi invasion may be consistent with a Fisherian rather than a bistable wave, increasing our understanding of Wolbachia infection dynamics in nature and supporting theoretical analyses [22] suggesting that the global distribution of Wolbachia

Author Summary
Wolbachia are bacteria that live within the cells of arthropod hosts and are widespread in many groups of insects. These bacteria can rapidly spread through a population through a process of cytoplasmic incompatibility whereby females uninfected by Wolbachia show embryo death when they mate with males carrying the bacteria. Because the infected females pass on Wolbachia to their offspring, this places them at a reproductive advantage, ensuring that the infection spreads through insect populations once it reaches a high enough frequency to overcome any negative fitness effects on its host. Yet while such a rapid spread has been predicted, it has rarely been observed in nature. Here we show that a Wolbachia infection of Drosophila simulans flies has spread very rapidly in eastern Australia, replacing another Wolbachia infection that has also spread in recent years. These invasions appear to have taken place from a very low frequency, implying that both infections are likely to have had a benefit to their hosts rather than a cost. These results have implications for the spread of Wolbachia infections currently being introduced into populations of mosquitoes and other insects for disease suppression.
throughout insects and other arthropods is most easily understood if these infections tend to spread deterministically from low initial frequencies.

Wolbachia Infection Frequency Estimates
We assayed individual D. simulans samples for Wolbachia infection status and strain type using a real-time PCR/high resolution melt (RT/HTM) method designed to amplify a fragment of the wsp gene [39]. Samples which successfully amplified using this method yielded amplicons with melting temperature peak rates (T m ) that clustered into two distinct groups which showed a consistent difference of ,0.5uC. Wolbachia strain type for positively infected D. simulans individuals was assigned on this basis. Positive controls where DNA extracted from both a wAu-infected and a wRi-infected individual fly was combined and amplified in the same reaction produced intermediate and less distinct T m peaks. This pattern was never observed for amplicons derived from individual specimens, suggesting that double infections are unlikely to occur amongst Australian field populations of D. simulans at any appreciable frequency. Sequencing of the amplicons derived from 13 of these samples using a standard PCR method with primers wsp_81F_Fwd: 59-TGGTCCAA-TAAGTGATGAAGAAAC-39 and wsp_691_Rev: 59-AAAAAT-TAAACGCTACTCCA-39 [40] confirmed that 9 samples from the high-T m cluster were the wRi haplotype (cf. GI: 225591853), whilst the remaining 4 from the low-T m cluster were wAu (cf. GI: 2687519). For a further 5 high-T m and 5 low-T m cluster individuals we also obtained sequence data using standard PCR methods for each of the Wolbachia MLST genes gatB, coxA, hcpA, ftsZ and fbpA [41], and for the sucB gene [42]. These samples were drawn from isofemale lines derived from diverse geographic locations including, for the high-Tm cluster, lines established in 2011 from Brisbane, Melbourne and Gosford, New South Wales; and for the low-T m cluster, recently sourced lines from Perth and Geraldton (Western Australia), Y6 originally from Cameroon, West Africa (E.A. McGraw, pers. comm.), and Coff1 sourced from Coffs Harbour, New South Wales in the mid-1990s bearing the originally designated wAu infection [33]. A single haplotype was obtained for the high-T m samples which showed 100% sequence identity with the wRi strain of Wolbachia (cf. GI:225629872). A single distinct haplotype was also obtained for the low-T m samples which was identical to previous data for the wAu strain at the sucB locus (cf. GI:84028372).
Further validation of these results was obtained using the RT/ HTM method with two sets of strain-specific primers (wRi_wsp_Fwd: 59-TGATGTTGAAGGGCTTTATTCACAG-39; wRi_wsp_Rev: 59-GTATCTGGGTTAAATGCTGCACCTG and wAu_wsp_Fwd: 59-TGATGTTGAAGGAGTTTATTCA-TAC-39; wAu_wsp_Rev: 59-TTTGCTGGGTCAAATGTTA-CATCTT-39). 34 of the original samples from the high-T m cluster each amplified successfully using the wRi-specific primers but did not amplify with the wAu primers, whilst the reverse was true for 27 samples from the original low-T m cluster.
Results for samples collected in 2004 (n = 162), 2008 (n = 499), and 2011/12 (n = 799) are summarised in Figure 1 (full details in Table S1). In contrast to the infection frequencies found for samples collected in 1993/94 [33; Figure 1], a second Wolbachia strain (wRi) was already prevalent amongst some east Australian D. simulans populations by 2004, both in the far north (Cairns) and the south (Melbourne); it occurred at a lower frequency further down the Queensland coast (Maryborough) and in northern Tasmania, but was not yet present in the central coastal area of northern New South Wales or in southern Tasmania. Additionally, for three central populations (Maryborough, Kingscliff and Red Rock), wAu infection frequencies increased significantly compared to geographically equivalent populations sampled previously (P,0.001, G-test). Ballard [43] also reports significantly higher wAu frequencies for samples collected in December 1999 from both Coffs Harbour and Brisbane than those previously found by Hoffmann et al. [33].
The 2008 data show wRi spreading considerably since 2004, being at high frequency throughout Queensland, and now detected in both southern Tasmania and 6 of 7 populations from New South Wales. By 2011/12, wRi was at high frequency at every location we sampled in mainland eastern Australia and Tasmania (n = 794) with an overall frequency of 0.929 (0.909, 0.946); whilst the wAu infection was only detected from an island (Hayman Island) ,29 kms offshore from Airlie Beach in Queensland. We did not detect wRi infection from two sites in Western

CI Assay
Laboratory crosses indicate that the wRi-infected Australian flies have a CI phenotype indistinguishable from California wRi. Hatch rates (Table 1) from crosses between infected females from four north Queensland lines established from females collected in 2011 and five-day-old males from a wRi-infected laboratory line provided no evidence for CI and were similar to the average hatch rate when females from these lines were crossed to males from an uninfected and a wAu-infected laboratory line. Females from these lines therefore behaved like wRi lines in crossing type. In addition, when males from these lines were crossed, they were compatible with the wRi line females, but generated incompatibility in crosses with females from both the uninfected and wAu-infected lines. All differences between reciprocal crosses were significant (Mann-Whitney U-tests) at the table wide a' = 0.05 level after corrections for multiple comparisons.

Maternal Transmission
To test maternal transmission of wAu in D. simulans in Western Australia, we established 55 isofemale lines from field-caught females obtained from two locations near Perth and from Geraldton in Western Australia in 2011/12. We assayed the Wolbachia infection status of the mothers and several F1 progeny, including between 10 and 12 progeny for each of 34 lines infected with wAu. From these lines, 350 progeny were tested, and all but 8 were found to be infected. This yields a field estimate for m of 0.023 with 95% bootstrapped confidence interval (bCI) of 0.003, 0.049 based on 10,000 replicates). The bCI is broader than the 95% binomial confidence interval based on total numbers (0.01, 0.045), but is more appropriate because of heterogeneity in transmission frequencies across females [cf. 26].
We also tested maternal transmission of wRi in eastern Australia using isofemale lines established from two locations near Melbourne in 2013. Ten F1 progeny were tested for each of 23 of these lines where the mother was found to be infected. All but 7 of the F1s were also infected, yielding a field estimate of m of 0.026 (95% bCI 0.004, 0.057) which is consistent with previous field estimates of m for wRi in California [17,26].

Fecundity
The wRi infection has previously been found associated with significant fecundity increases in D. simulans lines from California [37]. We performed a three-way comparison for fecundity using Australian wRi-infected, wAu-infected and uninfected laboratory lines which had been outcrossed to a common genetic background for 5 generations. Results (Figure 2) indicate that females from the wRi-infected line showed a marginally significant fecundity advantage relative to the uninfected line (p = 0.045, t-test) but strain differences were marginally non-significant in the three-way comparison (F 2,37 = 2.7, p = 0.082).

Mitochondrial Haplotypes
To examine the association between Australian wAu and wRi infections across time, we obtained sequence data for a 1256 bp mtDNA fragment from 7 wAu-infected and 11 wRi-infected samples, and found no variation apart from a single G to A base substitution previously identified by Ballard [43] (position 457 on our fragment, position 8201 in [43], see Table 2) with the former sequence (here denoted ''A haplotype'') found only in wAupositive samples and the latter (denoted ''R haplotype'') only in wRi-positive samples. Additional mtDNA haplotype sequence data obtained for 45 Wolbachia uninfected (w-) samples from 2004, 39 wsamples from 2008 and 6 wsamples from 2011 were similarly invariant, comprising only A or R haplotypes. The mtDNA haplotype for a further 9 wsamples from 2004, and 55 wsamples from 2008, was determined by treatment of PCR amplicons with the restriction enzyme HinfI. Full details of the mtDNA haplotypes determined for wsamples are summarised in Table S1. Both A and R haplotypes were detected amongst the 2004 and 2008 eastern Australian wsamples (simultaneously for some locations), whereas sequence data for wsamples collected in 2011 (n = 5) all confirmed the presence of the R haplotype. A single 2011 wsample from Perth, Western Australia had the A haplotype.
To better understand the history of association of wRi with D. simulans, we also examined mtDNA variation in California D. simulans stocks by sequencing two regions which included those identified in Table 3 of Ballard (2004) as showing the greatest variation among siII mtDNA haplotypes, positions 1175 to 1626 and 7838 to 8287. We examined mtDNA from a stock founded in 1961 in Nueva, California, ,25 km SE of Riverside, the type locality for wRi [28]. This stock is not infected with Wolbachia, and its origin pre-dates the initial 1985 Wolbachia analyses. It yielded a novel haplotype, denoted CA61, different from those described in Ballard (2004) by two base pairs from both the canonical ''R haplotype,'' which Ballard calls DSR, and the canonical haplotype (DSW) associated with northern California D. simulans prior to the invasion of wRi. To better understand the novel CA61 haplotype, we examined mtDNA from two reference wRi stocks, Riv84 and Riv88, collected in Riverside in 1984 and 1988, respectively, and from seven stocks (denoted Y26, Y35, Y36, Y42, Y46, Y54 and Y62) collected from an orchard approximately 5 km east of Winters, California in August 2011. All stocks apart from Y35 and Y62 were wRi infected. The Y54 stock produced a novel haplotype which differs from the canonical DSR by one nucleotide ( Table 1). As expected, the mtDNA from Riv84, Riv88, and all the other stocks matched DSR.
Given the paucity of sequence variation in siII D. simulans mtDNA sequences, and in the mtDNA associated with wRi [43,44], there is little statistical power to infer phylogenetic relationships among these sequences, apart from the fact that the sequences from the Indian Ocean seem to be part of a distinct clade (see Fig. 3 of [43]). However, with only one exception, all wRi-infected stocks share the nucleotide A at position 8201, in contrast to the G found in DSW and CA61 (as well as AU23 and all other wAu-infected stocks examined by Ballard [43]. This substitution corresponds to the HinfI polymorphism described by Hale and Hoffmann [45] differentiating mtDNA of wRi-infected stocks sampled worldwide from DSW. Nucleotide A was found at position 8201 in all wRi-infected stocks examined worldwide (Ballard [43], 92 wRi-infected stocks from California [46] and 29 of 30 additional wRi-infected stocks [17]. The only exception may correspond to rare paternal transmission of mtDNA [cf. 47]. Given that uninfected flies from wRi-infected populations quickly become associated with the mtDNA from their infected maternal ancestors [46], it seems unlikely that wRi was present at an appreciable frequency near Riverside in 1961.

Mathematical Inferences
As noted by Hoffmann and Turelli [12], for a Wolbachia variant such as wAu to persist despite imperfect maternal transmission and little CI, it must increase fitness sufficiently to offset imperfect transmission, i.e., F(12m).1. Let I A (I R ) denote individuals infected with wAu (wRi), and let U denote uninfected individuals. We measure fitness of I A and I R females relative to uninfected females, and denote their relative fitnesses F A and F R . Let m A (m R ) denote the frequency of U ova produced by I A (I R ) females. The condition for wAu to increase when rare is F A (12m A ).1. When this is satisfied, wAu will reach a stable equilibrium frequency of The wAu frequencies observed in southern Queensland and northern New South Wales in 1999 and 2004, as well as the 2011/ 12 wAu frequencies in Western Australia, all suggest an equilibrium frequency for wAu of roughly 0.6. Using (1) and our  assuming roughly 20 generations per year for 6 years. With F A = 1.061 and m A = 0.023, p A would increase from 0.054 to 0.529 in 120 generations, consistent with our data. If m A were significantly smaller, say 0.01, (1) withp p A = 0.6 implies F A = 1.026. Assuming m A = 0.01 and F A = 1.026, p A would increase in 120 generations from 0.054 to only 0.229, far below our 1999 estimate. However, if m A were appreciably larger, say 0.04, (1) witĥ p p A = 0.6 implies F A = 1.11. With those values, p A would increase from 0.054 to 0.597 in 120 generations. Thus, our frequency estimates and maternal transmission data are consistent with a positive fitness effect for wAu on the order of 5-10%, but inconsistent with appreciably smaller fitness advantages of 3% or less (Figure 3). Given the uncertainty of our estimates for both the apparent stable equilibrium frequency,p p A , and the maternal transmission rate, m A , for wAu, our data documenting the frequency increase of wAu from 1993 to 1999 provide the most robust estimate for the fitness advantage it seems to induce.
As shown in Materials and Methods, once wAu is established, wRi will tend to increase when very rare only if The roughly simultaneous spread of wRi from three geographically disparate foci in north Queensland, Victoria and Tasmania   To compare these data to the predictions of our simple model (see Eqs. 3 below), we need estimates of fitness effects, maternal transmission rates and CI. As above, we assume F A = 1.061 and m A = 0.023. Assuming that H, the relative hatch rate from incompatible fertilizations (i.e., sperm from wRi-infected males fertilizing either uninfected or wAu-infected ova), and the fidelity of maternal transmission were as observed in California, e.g., H = 0.55 and m R = 0.045 [17,26], we predict that starting with p A = 0.54 and p R = 0.09, p A will fall below 0.01 within 40 generations only if F R is on the order of 1.08 (Figure 4). After 45 generations, these parameter values imply that p R should rise from 0.09 to very near its predicted stable equilibrium value of 0.94, consistent with our data.
The stable equilibrium frequency of wRi is very insensitive to F R [27]. For instance, with H = 0.55 and m R = 0.045, as F R increases from 0.95 to 1.1,p p R increases from 0.93 to only 0.94. In contrast, the unstable equilibrium for wRi (when entering an uninfected population) plummets from 0.22 to 0, converting the predicted spatial dynamics from bistable to Fisherian. Our temporal data are too coarse to approximate accurately the speed of wRi's spatial spread. However, between 2008 and 2011, wRi spread southward from Coffs Harbour and northward from Bega (Figure 1), filling in roughly 1000 km of the New South Wales coast in three years. The speed of this bidirectional spread is comparable to the northward spread of wRi observed in California, roughly 100 km/ year [11]. Such rapid spatial spread can be explained by humanmediated, long-distance dispersal and Fisherian local dynamics that allow rare long-distance migrants to greatly accelerate spatial advance [25,Ch. 5].

Discussion
Despite the ubiquity of Wolbachia in natural populations of arthropods, there are very few documented examples of Wolbachia spread [10,11]. This has limited our understanding of Wolbachia spatial dynamics, with bistable spread for CI-causing Wolbachia suggested by imperfect maternal Wolbachia transmission by wildcaught females [15,29] and demonstrated fecundity costs for infected females in the laboratory [17,48,49]. However, bistable dynamics cannot explain the persistence of Wolbachia that cause little or no reproductive manipulation [12]. Moreover, as emphasized by Fenton et al. [22], bistable dynamics are difficult to reconcile with molecular data indicating that the widespread occurrence of Wolbachia is attributable to rare horizontal transmission events involving one or very few infected founders [3,44,50,51]. These phenomena are more easily explained by Fisherian dynamics [22], in which the frequency dynamics of rare Wolbachia infections are dominated by positive fitness effects. Our temporal and spatial data, describing the sequential spread of two Wolbachia strains (wAu and wRi) through natural populations of D. simulans along the east coast of Australia over the last two decades, support the view that both strains may have spread under Fisherian dynamics.
The recent spread of wRi is particularly compelling, with apparently independent introductions in both the north and south of the country between 1994 and 2004, leading to near-fixation of wRi in all populations sampled on the east coast mainland of Australia by late 2011/early 2012. The appearance of wRi in three geographically separated locations -Queensland, Victoria and Tasmania -in 2004 cannot be explained as a single bistable wave that originated and spread from elsewhere. In 1994, before the arrival of wRi, wAu was rare in several populations.  data from 1999 and our samples from 2004 and 2008 ( Figure 1) document wAu frequency increases across the Australian east coast. Given that wAu does not induce CI [33], its spread from low frequency in multiple populations implies that it confers a fitness benefit. The rapid displacement of wAu by wRi suggests that wRi must enhance fitness even more than wAu.
Both wAu and wRi have been detected in D. simulans populations around the world [43,52,53], and these infections have co-occurred previously (e.g. Sangoqui and Rocafuerta, Ecuador in 2000; Brooksville, FL, USA in 2002 [43], and possibly also Japan [17]). The wAu infection in D. simulans may have a Neotropical origin as a result of a relatively recent horizontal transfer event from Drosophila willistoni [54]. The origin of wRi in D. simulans is less certain, but very closely related Wolbachia have been found in various Drosophila that co-occur with this cosmopolitan species [e.g., 54,55,56]. It is not known how either strain was introduced into Australian D. simulans populations, but our data suggest that both introductions occurred within the last few decades. The pooled frequency of wRi (92.7%, Table S1) in our 2011/12 eastern Australian samples (including Tasmania) is consistent with the stable equilibrium frequency found in California over the past two decades [,93%; 17,26,37]. This equilibrium frequency is determined primarily by significant CI and imperfect maternal transmission. Our data from crossing experiments, and mtDNA and Wolbachia sequencing suggest that the wRi strain now found throughout eastern Australia is essentially identical to that initially identified in California [28] and subsequently found in Africa, Asia, Europe and South America [17,43]. Genome comparisons support this (M. Turelli, unpublished data). Our mtDNA data from uninfected individuals in 2011/12 indicate that the R haplotype is being swept to fixation in eastern Australia as a consequence of the spread (with imperfect maternal transmission) of wRi [cf . 46]. A previous mitochondrial sweep of the A haplotype associated with wAu may have also occurred, but the ancestral haplotypes in eastern Australian D. simulans are not known.
Wolbachia and mtDNA are expected to be co-inherited maternally. Even though Wolbachia exhibit imperfect maternal transmission (m.0), mtDNA haplotypes associated with a given infection will be driven to fixation, by CI or other Wolbachiaassociated positive fitness effects, if all offspring of infected females inherit maternal mtDNA [46]. Eventually, all individuals in the population have infected maternal ancestors who passed on their mtDNA even if not their Wolbachia. Although rare paternal transmission of Wolbachia has been detected in laboratory D. simulans [48,49], Turelli et al. [46] and Turelli & Hoffmann [17] inferred that paternal transmission or horizontal transfer must be rare in nature. Similarly, there is evidence for rare paternal inheritance of mtDNA in D. simulans [47,57,58], but it is apparently too rare to create incongruity between Wolbachia and mtDNA lineages [43]. Drosophila simulans populations worldwide carry three distinct mtDNA haplotype groups, with very little intra-haplogroup variation [44,59,60]. Ballard [43] found both wAu and wRi only in flies carrying siII haplogroup mtDNA. He identified a single (synonymous) G to A base substitution (see Table 2) that was consistently different between wAu-infected versus wRi-infected flies respectively, and matched a restriction enzyme polymorphism identified by Hale and Hoffmann [45]. In our analyses, all wRi flies had the same base substitution as identified by these researchers, consistent with the hypothesis of a unique origin for wRi in D. simulans.
Turelli and Hoffmann [11,17] had previously suggested bistable spatial dynamics of wRi in California populations of D. simulans. However, the rate at which wRi spread spatially, on the order of 100 km per year, is more consistent with Fisherian dynamics whose deterministic spread from very low frequencies allows rare long-distance dispersal to greatly accelerate spatial advance. Assuming that the mtDNA haplotype of the 1961 strain we characterized was widespread in southern California, our new analysis suggests that wRi was not yet at appreciable frequency at this time, and that the invasion detected by Turelli and Hoffmann [11] may reflect a rapid spatial spread from populations farther south.
The speed of wRi spread throughout eastern Australia is comparable to the spread of wRi in California, on the order of 100 km/year [11]. The infection was common in two populations sampled in 2004 (Melbourne and Cairns), then swept to near fixation in all east Australian mainland populations and Tasmania by 2011/12. Although some genes and traits exhibit clinal variation over this range [61,62], several molecular markers indicate that D. simulans is effectively a single panmictic unit throughout eastern Australia, including Tasmania [63]. Thus, there are no apparent barriers to dispersal that would isolate populations from the spread of wRi. We did not, however, detect wRi in Western Australia, with the wAu infection persisting at relatively high frequency (<60%). Western Australian populations of D. simulans may be geographically isolated from eastern Australian populations, although once wRi is introduced it is likely to sweep through this area in the coming decades.
The Drosophila-Wolbachia dynamics have been used to interpret and model the spread of Wolbachia-infected Aedes aegypti mosquitoes in Australia [7]. In contrast to the apparent Fisherian spatial dynamics of wAu and wRi, these releases seem to follow bistable dynamics, with strong evidence for a non-trivial unstable equilibrium frequency, on the order of 10-20% [6,7]. After wMel was driven to a high frequency in two relatively isolated populations by 10 weeks of releases, it increased to near-fixation, and has remained at over 90% frequency for the past 2 years [7; unpublished data]. Although the infection was detected in disjunct residential areas near the release sites soon after the initial releases, it has not persisted or spread in these areas [7; unpublished data], as would be expected under Fisherian dynamics. Given the evolutionary pressure for Wolbachia to evolve towards mutualism with its hosts [37,64], a decreasing unstable point may alter these dynamics in the future [13] particularly if fitness costs produced in novel hosts are overcome.
In summary, our results suggest that the wAu and wRi Wolbachia infections of D. simulans spread at least part of the time in a Fisherian manner. The wAu infection increased in frequency despite having imperfect transmission and not causing CI. Its spread and apparent equilibrium make sense only if it increased host fitness. Over the past decade, wAu was rapidly displaced by wRi throughout eastern Australia. The rapidity of this replacement starting from three locations in 2004 indicates that wRi also increases host fitness in nature, at least under some circumstances. These data plus evidence that wRi entered southern California after 1960 lead us to reinterpret that canonical example of Wolbachia spread in California as Fisherian rather than bistable. Depending on the rapidity with which Wolbachia adapt to novel hosts, its spread after artificial introductions to target species may be greatly accelerated. . Further samples were collected from mainland eastern Australia at similar coastal and some inland localities, and from the Perth region of Western Australia during 2011, and from both Geraldton (Western Australia) and Tasmania in early 2012. All samples were preserved in 100% ethanol and stored at 220uC. Specimens subsequently used in PCR analysis were males except where a few isofemale lines had been established prior to preservation of the mother. Species identity was established by male genital arch morphology and confirmed by molecular methods [39].

Field Collections
The new California samples were from isofemale lines (Y26, Y35, Y36, Y42, Y46, Y54 and Y62) established in August 2011. The collection site was a peach orchard approximately 5 km east of Winters, CA. The stocks were maintained as isofemale lines until analysed. The 1961 Nueva, California stock was obtained from the Drosophila Species Stock Center at University of California, San Diego (stock #14021-0251.006).

Wolbachia Infection Status
DNA extractions were performed using a standard Chelex based method, and assays for Wolbachia infection status and strain type were performed with a RT/HRM method using the Roche LightCyclerH 480 system as previously described [39]. Briefly, a conserved set of primers (wsp_validation_Fwd: 59-TTGGTTA-CAAAATGGACGACATCAG-39 and wsp_validation_Rev: 59-CGAAATAACGAGCTCCAGCATAAAG-39) were used to target a variable ,340 bp region of the wsp gene. The wRi and wAu alleles are expected to differ by 22 SNPs and a single 3 bp indel over this wsp region, potentially yielding a significant difference in observed T m for the respective PCR amplicons, with the wRi allele expected to have the higher average T m [39].
Samples from D. simulans isofemale lines previously identified as wRior wAu-infected respectively, were included on each PCR plate as Wolbachia strain-type positive controls. Further positive controls entailed equal volumes of DNA extracted from an individual fly of each type being combined in a single reaction. Separate D. simulans specific primers (Dsim_RpS6_Fwd: 59-CCAGATCGCTTCCAAGGAGGCTGCT-39; Dsim_RpS6_Rev: 59-GCCTCCTCGCGCTTGGCCTTAGAT-39) were used as a host species-specific control for each sample.
The RT-PCR/HRM conditions were as follows: 10 minutes at 95uC; 45 cycles of 10 seconds at 95uC, 15 seconds at 58uC, 30 seconds at 72uC; 1 minute at 95uC; 1 minute at 40uC. High resolution melting; temperature ramped up to 95uC at 0.02u/second with continuous fluorescence acquisition; 30 seconds at 40uC.
Individual samples were scored as Wolbachia infected, uninfected or inconclusive based on machine reported wsp and RpS6 crossing point (Cp) relative and absolute value criteria and examination of individual T m profiles for evidence of non-specific products. A relatively small number of samples, which initially yielded inconclusive results, were retested.
Additional standard PCR amplification of wsp, gatB, coxA, hcpA, ftsZ, fbpA and sucB gene fragments for selected individual samples followed [41] with annealing temperatures of 54uC for gatB, coxA, hcpA, ftsZ and sucB, and 59uC for wsp and fbpA. PCR products were checked using a 2% agarose gel for presence of an unambiguous single band of expected size. Amplified DNA for selected samples which yielded a clear and positive PCR result were sent to Macrogen (Korea) for purification and sequencing. Sequence data obtained was analyzed using Geneious v6.1 (Biomatters, Auckland, NZ).

mtDNA Analyses
Assays for mitochondrial DNA haplotype for a total of 172 individual Australian samples were performed using primers (Dsim_siII_7765_Fwd: 59-ATTTAATATTCAAGCAATAGC-39; Dsim_siII_8981_Rev: 59-TTCTGGTTCTATAATTTTAGC-39) designed to target a 1256 bp region of the D. simulans mitochondrial genome previously identified as containing at least one (synonymous) G to A base substitution at position 457 (outer strand) that was a characteristic difference between the haplotypes found to be associated with either the wAu or wRi strains of Wolbachia, respectively [43].
Amplification of mtDNA was performed using a standard PCR method with the following conditions: 5 minutes at 94uC; 38 cycles of 30 seconds at 94uC, 1 minute at 55uC, 1 minute at 72uC; 5 minutes at 72uC; hold at 4uC. PCR products were checked using a 2% agarose gel for the presence of an unambiguous single band of expected size. Amplified DNA for 108 samples that yielded a clear and positive PCR result were sent to Macrogen (Korea) for purification and sequencing. Sequence data obtained was inspected using Sequencher v4.5 (Gene Codes, Ann Arbor, Mi). Subsequently PCR amplicons for a further 64 samples were checked on a 2% agarose gel after digestion with the restriction enzyme HinfI [45].
For each California line, DNA was extracted from 30 flies as described in [65]. Following the position descriptions in [43], two sets of primers were used to amplify two mtDNA regions from positions 1034 to 1718 and from 7797 to 8425. The second region contains the G to A base substitution at position 8201 that distinguished the R haplotype from the A haplotype, associated with wRi and wAu infection, respectively. (The primers used were: region 1:  DSRmt1033+:  59-CCAAAATGACTTGTAATCCA-39  and  DSRmt17192: 59-GCACCTAATATTAAAGGCACT-39, region  2: DSRmt7796+: 59-GCTACATCTCCAATTCGATTA-39 and  DSRmt84262: 59-TTTATATTCTTTTAGACAACATGG-39.) The PCR conditions were: 3 minutes at 94uC; 34 cycles of 30 seconds at 94uC, 30 seconds at 55uC, 75 seconds at 75uC; 8 minutes at 72uC; hold at 10uC. The PCR products were checked on a 2% agarose gel for an unambiguous band of the correct size. The PCR products were purified using the QIAquick PCR purification kit. The amplified fragments were Sanger sequenced by the UC Davis College of Biological Sciences UCDNA Sequencing Facility, using the same primers used for the PCR reaction. The sequencing results were analysed using ApE alignment tool (A plasmid Editor v2.0.45, by M. Wayne Davis).
To confirm the sequence coordinates we realigned the mtDNA genomes described by Ballard in [43,44,66] also including the D. yakuba mtDNA reference genome [67], and confirmed the positions of polymorphisms described in Table 3 of [43]. We then aligned our new sequences against this reference.

CI Assay
Three separate laboratory lines previously confirmed to be wRiinfected (Riv88), wAu-infected (Coffs 1) and uninfected (W88) respectively were used in reciprocal crosses with each of four separate lines (CBQ40, CBQ46, CBQ72 and CBQ80) established from field caught females collected in 2011 from north Queensland to test for CI. The level of CI was determined by mating virgin 5 dold males to virgin females (5-8 d old). Males were mated once, and females were placed after mating in a vial with a spoon containing 5 ml of agar-treacle-yeast medium and left for 24 h at 25uC. The number of unhatched eggs was counted 24-32 h later. CI data (egg hatch rates) were angular transformed prior to analysis. Mann-Whitney U-tests were used to compare CI levels between the different lines. Multiple comparisons were corrected at the table wide a' = 0.05 level using the Dunn-Sidak method [68].

Fecundity
Separate mass bred wAu-infected and uninfected fly lines were established from 20 to 30 Western Australian isofemales lines for which infection status had previously been determined. Virgin females from each of these mass bred lines were then mated with a similar number of field caught Brisbane males which had been aged in the laboratory for at least 5 d to reduce the impact of CI. An equivalent mass bred wRi-infected line was established from virgin female F1 progeny of field collected Brisbane females mated with males from the uninfected mass bred line. Each mass bred line was then outcrossed to field collected aged Brisbane males which were also aged for at least 5 d for a further 4 generations to homogenise the genetic backgrounds. Lines were then retested to confirm infection status. Flies were reared at low densities by transferring 25 eggs into vials on ,20 ml of agar-cornmeal-yeast medium. Pairs of virgin females and males were transferred to fresh medium vials for 2 d, then females were transferred to vials with spoons containing 5 ml of medium with 20% fresh live yeast paste added to the medium surface. Spoons were replaced every 24 h for 5 d and eggs counted. Twelve to 13 females were assayed for each line. Model I ANOVA (analysis of variance) and t-tests were used to compare fecundity between lines.

Mathematical Analyses
We address two issues: first, the conditions for the initial spread and maintenance of wAu in an isolated population, then the conditions for its displacement by wRi. Let I A (I R ) denote an individual infected with wAu (wRi), and let U denote an uninfected individual. Assuming discrete generations, we denote the frequency of these three types of adults in generation t by p A,t , p R,t and p Ø,t . No individuals doubly infected for wAu and wRi have been found (and each infection is associated with distinct mtDNA haplotypes, [43]), hence we assume p A,t +p R,t +p Ø,t = 1. We measure fitness of I A and I R females relative to uninfected females, and denote their relative fitnesses F A and F R . In principle, these Wolbachia may affect viability or fecundity. Given that our data indicate relatively small fitness effects (see Results), we can simplify the analysis by approximating fitness differences as female fecundity variation with no viability effects.
Field-collected I A and I R females have both been shown to produce uninfected offspring. Let m A (m R ) denote the frequency of U ova produced by I A (I R ) females. Before wRi arrives in the population, p A,t +p Ø,t = 1; so we need only keep track of p A,t . Because wAu causes no CI, Assuming F A (12m A ).1, this produces the stable equilibrium (1) in Results. (Equilibrium (1) is equivalent to the standard mutationselection equilibrium for haploids, in which the uninfected, less-fit type is maintained at frequency m/s, where m = m A is the fraction of uninfected ova produced by infected mothers and 12s = 1/F A is the relative fitness of U females.) Equilibrium (1) implies that ifp p A and m A ,,1 are known, F A <1+m A /(12p p A ).
To understand the dynamics of wRi entering a population with wAu, note that both Hoffmann et al. [33], and James and Ballard [53] demonstrated that wAu infection is not able to rescue CI caused by wRi. Indeed, wRi induces the same level of CI in matings with both uninfected (U) and wAu-infected (I A ) females. The joint dynamics of p A,t , p R,t and p Ø,t follow a simplified version of the recursions for double and single infections derived by Hoffmann and Turelli [ Note that H H in (4d) is just the sum of the right hand sides of (4a-c). Each of these three expressions has a simple interpretation as the product of two terms, the first (in square brackets) proportional to the fraction of ova of each type, the second proportional to the fraction of fertilized ova that hatch. (In (4b) this second term is one, because wRi-infected ova are compatible with all sperm.) Given that the three frequencies sum to 1, there are only two independent variables. Now consider the conditions for wRi to increase when it is extremely rare and wAu is at equilibrium with U, as described by Eq. (1). When p R,t <0, H H<p Ø,t +F A p A,t . When wAu and U are at equilibrium, Eq. (3) implies that p Ø,t +F A p A,t = F A (12m A ). Hence, from (4b), we see that the condition for wRi to increase when very rare is condition (2) provided in Results. This is just the condition for spread of a Wolbachia variant that is completely compatible with the existing type [64]. Because essentially no CI occurs when wRi is very rare, the incompatibility of wRi with wAu does not enter condition (2).
When (2) is not met, p R,t increases only once it exceeds a frequency threshold so that the fitness advantage I R individuals receive from CI overcomes their frequency-independent disadvantage relative to variant A. This condition is a simple generalization of the unstable equilibrium that arises for a single Wolbachia infection that satisfies H,F(12m),1.