Gene-drive-mediated extinction is thwarted by population structure and evolution of sib mating

Abstract Background and objectives Genetic engineering combined with CRISPR technology has developed to the point that gene drives can, in theory, be engineered to cause extinction in countless species. Success of extinction programs now rests on the possibility of resistance evolution, which is largely unknown. Depending on the gene-drive technology, resistance may take many forms, from mutations in the nuclease target sequence (e.g. for CRISPR) to specific types of non-random population structures that limit the drive (that may block potentially any gene-drive technology). Methodology We develop mathematical models of various deviations from random mating to consider escapes from extinction-causing gene drives. A main emphasis here is sib mating in the face of recessive-lethal and Y-chromosome drives. Results Sib mating easily evolves in response to both kinds of gene drives and maintains mean fitness above 0, with equilibrium fitness depending on the level of inbreeding depression. Environmental determination of sib mating (as might stem from population density crashes) can also maintain mean fitness above 0. A version of Maynard Smith’s haystack model shows that pre-existing population structure can enable drive-free subpopulations to be maintained against gene drives. Conclusions and implications Translation of mean fitness into population size depends on ecological details, so understanding mean fitness evolution and dynamics is merely the first step in predicting extinction. Nonetheless, these results point to possible escapes from gene-drive-mediated extinctions that lie beyond the control of genome engineering. Lay summary Recent gene drive technologies promise to suppress and even eradicate pests and disease vectors. Simple models of gene-drive evolution in structured populations show that extinction-causing gene drives can be thwarted both through the evolution of sib mating as well as from purely demographic processes that cluster drive-free individuals.


INTRODUCTION
Engineered gene drives offer an exciting new technology for the possible control of pests and vector-borne diseases, and which might even be used to rescue wildlife species from the edge of extinction. The selective advantage of some drives is so powerful that they can be used to cause species extinction, but many potential applications propose using them more benignly, to deliver a harmless genetic cargo throughout a species.
Gene drives are classic selfish genes that give themselves an 'unfair' advantage in meiosis, gamete competition, or offspring competition, often to the detriment of the population or family. It seems paradoxical that natural selection can favor genes with such extreme deleterious effects that they could cause extinction, but the theory indicating such possibilities is over half a century old, developed in response to a few natural systems of selfish genes [1][2][3][4]. The discovery of gene drives in natural populations was followed quickly by the realization that they could be used to transform species for various socially desirable ends [5] or even extinction [4], but engineering to implement gene drives remained the challenge for decades. Now, the insight of Burt [6] combined with CRISPR technology has led to a revolution in interest [7][8][9][10][11]; laboratory experiments have now shown the feasibility of various implementations [12][13][14]. The enthusiasm for gene-drive extinctions is reflected in such profound dreams as possibly eradicating Anopheles mosquitoes that carry malaria [14] and eliminating all invasive mammals in New Zealand [15].
But will it work? The pace at which engineering methods have enabled gene-drive construction has vastly exceeded our experience with implementations, so that we stand poised to introduce gene drives on a massive scale without appreciating how they might fail or deviate from expectations. Given the demonstrated success of engineered gene drives in experimental populations, the most obvious basis of possible failure now becomes the evolution of resistance, the focus of this article. At a minimum, resistance would limit coverage of the population by a gene drive; at worst, resistance would fully reverse a drive's effect. Furthermore, the evolution of resistance to one implementation may thwart subsequent implementations, so early failures may have longterm ramifications for later interventions. There is thus an imperative to understand resistance evolution before implementing gene drives on a wide scale.
Resistance may take many forms. Some forms may be specific to the mechanistic underpinnings of the gene drive implementation; others may operate largely independent of the drive mechanism. For a homing endonuclease gene, an obvious form of resistance is mutation in the target sequence recognized by the nuclease [6]. Resistance that, at a molecular level, blocks genedrive expression or interferes with its operation will be difficult to predict or study except empirically, in the context of specific applications. Other types of resistance, especially those that transcend mechanistic details of the drive, may be more amenable to a priori analysis. Furthermore, CRISPR is not the only method successfully used to engineer gene drives [16], further warranting the consideration of types of resistance evolution that are mechanism-independent.
Here we use mathematical models to specifically consider population structure as a foundation for resistance or extinction failure-how specific types of non-random mating will operate and evolve to thwart extinction-causing gene drives; extinctioncausing drives encounter the strongest selection for resistance and possibly face the greatest potential for failure from the demographic consequences of population reduction. Gene drives, indeed perhaps all selfish genes, benefit from mass action processes that facilitate their interactions with alternative alleles. Mass action (classic random mating) pits alleles against each other in proportion to their abundances. Population structure, which may take various forms, tends to cluster alleles of common origin to compete against themselves, not only constraining the scope of their selfish benefits, but also increasing the degree to which they compete against themselves instead of against the wild-type alleles.
Although this article relies heavily on mathematical models for insight, we recognize that it is impossible to capture the full empirical complexity of any gene-drive implementation in a mathematical model. Even such basic properties as the interplay between genetic evolution and demography are unknown and may vary through time and space in a real population. We therefore offer a variety of highly simplified models intended to capture different features of the process. Specifically, we develop two classes of models. One is a metapopulation model with simplified dynamics regulating local extinctions (as if from a gene drive), recolonization, and interactions between drive and non-drive populations. It is merely a modified version of Maynard Smith's original haystack model [17]. The second class of model is borrowed from formal population genetics and allows sib mating to evolve in response to a gene-drive introduction. Precedents for investigating inbreeding as a possible mitigation against gene drives are provided by Hamilton [4], Bull [18] and Drury et al. [19]. Our study goes beyond this prior work to offer a broad analysis of the consequences of population structure on extinction failure, whether due to evolution of resistance or just pure demographic consequences. One emphasis of specific interest is whether and how easily sib mating can evolve to block gene drives. Furthermore, we compare the effects of sib mating and its evolution for both Y-chromosome gene drive, which does not kill individuals, and for recessive-lethal gene drive, which does kill. The different models represent alternative structurebased processes that might ensue from an extinction-causing gene drive, and they are thus complementary.

EXTREME POPULATION STRUCTURE WORKS AGAINST GLOBAL EXTINCTION: THE HAYSTACK MODEL REVISITED
Many intended uses of gene drives rely on the drive suppressing or even extinguishing populations. The populations most suited to this end are unstructured with random mating. As is well known from the decades-old theoretical literature on group selection of cooperative (altruistic) traits, a structured population is protected against selfish elements, the degree of protection depending on quantitative details of fitness effects, migration rates, and group extinction rates [20][21][22]. When the selfish element is a lethal gene drive, it accelerates extinction of the subpopulations in which it resides, but the then-empty patches are recolonized disproportionately by the 'altruists' (wild-type) lacking the selfish element.
To illustrate this process in a highly simplified but intuitively tractable form, we modify the original haystack model of Maynard Smith [17], with parallels to Hamilton and May's dispersal model [23]. We imagine many small populations, each inhabiting patches (islands) in a large habitat of many patch sites; i.e. a metapopulation (Fig. 1). Migrants from one population can colonize other sites. There are two types of populations: those consisting purely of wild-type individuals (B, for beneficial) and those with at least some gene-drive individuals (S, for selfish). In the spirit of the haystack model, gene drive spread within a local population is considered so rapid that any patch with even a few gene-drive individuals is immediately converted to type S. This separation of time scales is an essential feature of metapopulation models and allows for a consolidation of state variables. Although migrants from B cannot convert S populations and thus can be ignored, the reverse migration (rate SB ) is highly effective because of the gene-drive effect. Furthermore, high values of SB represent a near absence of population structure. The model does not specifically include or even require genetics, the key within-patch process being the rapid takeover of B patches by S individuals when they invade B.
Differential equations describing the dynamics are presented in Equation (1), with variables and parameters defined in Table 1.

Parameter constraints
The very nature of extinction-causing gene drives dictates that our interest is confined to parameter values where B alone can persist but S alone cannot (resulting in the conditions B > B ; S > S ), a higher colonization rate of empty patches by B than by S ( B > S ), and a higher extinction rate of S patches than of B patches ( S > B ). Biology also dictates that all parameters must be non-negative, but there are no upper limits except as imposed by the aforementioned constraints. With these conditions, the system has two relevant equilibria: an equilibrium in which S is absent.
and an internal equilibrium of Stability conditions imply thatŜ 2 > 0 is required for S to invade. Approximately 15% of parameter space satisfiesŜ 2 > 0 when the aforementioned constraints are satisfied (based on a systematic search of parameter space across the range of (0.01, 2.02) in increments of 0.03 for all five parameters). This 15% value has no quantitative meaning biologically, as we do not know what parameter values would apply in nature; it merely indicates that there is a plausible parameter space in which coexistence is possible. Even when S invades, it cannot drive B extinct with these deterministic processes because there is always a threshold value of S below which B is so unaffected by S (by low colonization) that it persists. However, a choice of SB sufficiently high can push the equilibrium value of B to such a low value that B would not persist in any habitat with a finite number of patches. other-e.g., a B patch can become an S patch, but not the reverse. The rates at which one patch type is converted to another are given by the terms on the arrows, corresponding to Equation (1) Evolution, Medicine, and Public Health Although the haystack model was originally used for insight about group selection, our version here is probably better thought of in the usual metapopulation sense as a model of implicit spatial structure. The resulting spatial segregation is enough to allow coexistence of B and S (under appropriate conditions on the parameters). The model is so highly simplified as to be useful chiefly for intuition and insight to broad classes of outcomes. Its main purpose is to discover whether demographic structure alone can limit gene drive 'success' in extinguishing populations. In this respect, our results mirror those of North et al. [24], who simulated an extinction-causing gene drive in a spatial population with many details specific to mosquito biology. They found that small patches of mosquitoes could escape the extinction wave, suggestive of the patch model here.
To complement this demographic approach, we turn to population genetic models of resistance evolution. These models allow a specific type of structure to evolve, sib mating. In contrast to our haystack model, the sib-mating models operate without extinction-evolution is driven by fitness effects at the individual or family level.

SIB MATING
The above simple metapopulation analysis showed that implicit structure, which segregates a population into smaller groups of interacting individuals, can suppress the expansion of a gene drive, either limiting its spread to a fraction of the population or eliminating it altogether. But groups per se aren't the only structures that have the potential to thwart an extinguishing gene drive. Structure may exist at the family level in the form of inbreeding. It is not that inbreeding need be abundant in a species at the time of a gene-drive release; rather it may exist at low levels across a population or be more prevalent in a population in some conditions than in others. The question is whether inbreeding might evolve to higher levels in response to the gene drive.
A previous theoretical study found that selfing, the most extreme form of inbreeding and possible only in hermaphrodites, can be selected in response to a recessive lethal gene drive and that selfing limits the potential for extinction of the population targeted by the drive [18]. This result is worrying for gene-drive implementations, but there were two hopeful outcomes from that work. First, although evolution of even partial selfing could prevent gene-drive fixation, mean fitness of the targeted population was limited by the magnitude of inbreeding depression. Thus, mean fitness remained low if inbreeding depression was high, preserving much of the intended effect of the gene drive. Second, selfing could sometimes only evolve by major mutations, not small ones. That study considered only recessive lethal drive, and the latter result was evaluated only for the case of drive in one sex only.
We expand upon the selfing work here to address two questions for a different form of inbreeding-sib mating: (i) Does sib mating also evolve as a block to extinguishing gene drives, and can it evolve in small steps? (ii) Do recessive lethal drives and Y-drives equally favor sib mating?
For simplicity, we make the following assumptions. All models assume a life cycle with sexual haploids: male and female parents mate and produce a brief diploid phase, which then undergoes meiosis to produce haploid offspring. Gene drive operates in the diploid phase (regardless of which parent contributes the drive allele in the recessive-lethal model). Drive is always complete, with 100% of the progeny receiving the drive allele. The locus controlling sib mating has two alleles and is unlinked to the drive locus.

Genotypes and phenotypes
For models with genetic control of sib mating, the family's level of sib mating is controlled by the mother's genotype at the A/a locus ( Table 2). This biology simplifies the mathematics, and we do not suspect it has a qualitative effect on the outcomes. There are in fact many ways in which the mother could influence mating among her progeny, such as by ovipositing eggs in clusters, delivering molecules to the eggs that facilitate synchronous hatching, or choosing remote oviposition sites that limit opportunities for encountering non-sibs. Sib mating is limited to families that produce both sexes; with Y-drive, some families are all sons with no sisters to mate. In all trials illustrated here, sib mating for allele 'a' was set to 0-purely outcrossing-and 'A' provided some sib mating. A non-genetic ('ecological') class of models tested here allows the level of sib mating to change dynamically with mean fitness (w 1). The biological justification is that, as mean fitness declines, so will population density, and siblings may increasingly provide the only potential mates. Here, there is no genetic variation for sib mating (all genotypes are 'a' or, equivalently, s a = s A ), and we used an exponential sib-mating function with a single parameter (c) to control steepness: Drive genetics Two models of drive are studied. Recessive-lethal drive (alleles D/ d) operates so that if either (but only one) parent carries D, all offspring inherit D. If both parents carry D, no offspring are produced. That is, D is a recessive lethal only in the brief diploid phase. For the model of Y drive, all females carry X, but males carry either Y or Z chromosome. Of a male that carries Y, half his offspring inherit the Y and half inherit the mother's X. Of a male that carries Z, all his offspring inherit the Z and are thus sons in a brood with no sisters, and hence no possibility of sib mating.

Inbreeding depression
The models assume that the relative brood size of parents who are sibs is , typically lower than that of outcrossed offspring ( < 1).
(Inbreeding depression would be represented as the decrement in fitness, ¼ 1 À .) Inbreeding depression is assumed to be invariant throughout the evolutionary process. In real systems, inbreeding depression is often partially purged upon extended inbreeding, but allowing inbreeding depression to be static is a reasonable starting point and, if anything, provides a conservative measure of the vulnerability of gene-drive systems to be suppressed by inbreeding.

Male reproductive versatility
The net reproductive output of sons from a family with sib mating can be modeled in different ways, each of which may be observed in nature. At one extreme, sons who mate their sisters may then go on to join the random mating pool with no adverse consequence to their abilities in the outcrossing pool. At the other extreme, sons who mate their sisters are forever lost to the outcrossing pool, as if there is a brief time in which all mating occurs and an individual can be at only one place during that time. This latter process, of sons being 'discounted', is conveniently represented by assuming the extreme case that the fraction of a family's sons lost to the outcrossing pool is the same as the fraction of daughters mated by sons. We use K to represent the probability that a sister-mating male also participates in the random-mating pool.

Four models
The preceding account has identified two fundamental biological differences that require specific models: (i) The allele with drive is an autosomal recessive lethal or a Y (Z) chromosome (ii) Sib mating is controlled genetically or ecologically.
Addressing these variables in all combinations, there are four models to study. Equations are given in Appendix 1; all models consist of difference equations that assume discrete generations. Drive was assumed to be complete in all cases.

Genetic control of sib mating
There are minimally 16 viable family types that must be counted: 4 initiated from sib mating, 12 from outcrossing (Appendix 1; formally, there are 20 mating types, but 4 produce no progeny). The drive allele (D) is present in 8 of the outcrossed family types, and for these families, all progeny carry D so any sib mating is nonproductive. Sib-mating rates depend on the mother's genotype at the a/A locus, and the rate of sib mating induced by 'A' was varied systematically in different trials. Both the drive allele D and sibmating allele A were introduced at low frequency at the beginning of each trial. Sib mating rates were unaffected by the presence of D in the brood, so D led to offspring death from sib mating.
Equilibrium. We describe equilibrium outcomes based on mean fitness, which in our case is the average number of daughters per mother within a generation-all daughters are assumed to be mated. (For the equations in Appendix 1, this calculation is wðtÞ when b = 2.) These numerical studies revealed that the sib mating allele always evolved, even when it effected only a small level of sib mating (e.g. s A ¼ 0:01). However, average fitness at equilibrium was strongly dependent on the magnitude of sib mating encoded by allele 'A' up to a value of s A ¼ 0:5 (Fig. 3). With s A < 0:5, allele 'A' fixed and mean fitness remained below . If instead, s A > 0:5, allele 'A' remained polymorphic, and mean fitness equaled . In all cases, mean fitness was bounded by , the fraction of maximum brood size attained with parents who were sibs. The foregoing results apply to complete male discounting (K = 0-males who mate their sisters are lost to the random pool). In the absence of male discounting, mean fitness was observed to exceed when s A > 0:5, but the largest mean fitness observed was 1.3 , and the effect diminished as increased above 0.5. In comparison to the case of complete male discounting, a higher mean fitness with no male discounting is understandable because of the male fitness gained when sib-mating males later join the randommating pool.
The main qualitative result is that, although sib mating evolves in response to a lethal gene drive, mean fitness is approximately bounded by the fitness consequences of sib mating () and also somewhat bounded by the magnitude of sib mating allowed by the genetics.

Ecological adjustment of sib mating
In this model, all evolution is limited to the drive locus because genetic variation in sib mating is not required when the level of sib The female sib mating allele (not depicted) determines the fraction of the brood (here s A ) that is sib mated and the fraction 1 À s A that go to the outcross pool along with gametes from other broods that are available to outcross. A fraction K of males that sib mate a sister also join the outcross pool to possibly mate some more. Contributions to and consequences of the outcrossing pool are shown on the right. (Not all possible family types are illustrated.) Bottom: The mating of a non-drive female and male results in all offspring lacking the drive; the mating of a drive female and a drive male results in no offspring; a mating between one drive parent and one non-drive parent results in all offspring carrying the drive. The influence of is not shown mating increases 'ecologically'. We assumed that sib mating increases as mean fitness declines (Fig. 4, left). There is thus a parallel here with our haystack model (and a difference from the genetic control model) in that gene-drive evolution is being affected by purely demographic consequences of the gene drive's effect on population size. The two types of demographic models do not parallel each other with respect to the type of demographic effect; however-one operates via extinction, the other operates by affecting the level of sib mating.
The evolutionary dynamics are intuitive: the drive allele spreads and depresses mean fitness, thereby increasing sib mating. The increase in sib mating affects spread of the drive allele, thereby limiting further drops in mean fitness or possibly increasing mean fitness. A balance may be reached-dynamic equilibrium-or oscillations may result (Fig. 4). The shape of the sib mating function is critical both to the equilibrium as well as to the outcome of oscillations versus static equilibrium. Mean fitness depends heavily on the shape of the sib-mating function and again on -the fitness of inbred progeny.

Summary
For the genetic model with a recessive lethal drive, the results closely follow those with selfing [18]: sib mating evolves and partly limits gene-drive frequency, possibly avoiding extinction. Sib mating does not restore fitness to wild-type levels, and indeed, the equilibrium fitness cannot be (much) higher than the fitness of a purely sib-mated population.

Y drive
The case of Y drive presents some interesting contrasts to the case of recessive lethal drive. First, the drive allele does not kill any  following more than 1000 initial generations. Note that, in contrast to the genetic models, w= remains well below 1 (this ratio must be inferred from the graph); the ratio depends on and on the shape of the sib-mating function, so it may not obey any simple rule. Full male discounting was assumed (K = 0) individuals; it merely changes the sex ratio so that would-be daughters are instead sons. Second, sib mating is impossible in families with all sons, so a choice must be made whether such males face a penalty as outcrossers if their mother would otherwise have enforced sib mating. Our main interest is whether the evolution of sib mating will be robust to these differences from the recessive lethal case.

Genetic control of sib mating
The model setup is similar in many ways to that of recessive lethal drive (Fig. 5). For convenience, the Y with complete drive is denoted here as Z. There are four types of sib mated families and 8 types of outcrossed families (Appendix 1). Z is present in four types of outcrossed families, and as those families produce only sons, there is no opportunity for sib mating. To provide maximum benefit to the gene drive, we allow all sons carrying Z to join the outcrossing pool. Equilibrium. As with recessive lethal drive, sib mating was found to evolve under Y drive, the details differing somewhat from the case of recessive lethal drive. For example, allele 'A' always evolved to fixation (allele 'a' was for strict outcrossing). Mean fitness relative to closely paralleled the sib mating level (Fig. 6). (Mean fitness was calculated the same as for recessive lethal drive.)

Ecological control of sib mating
Qualitatively similar results were obtained for ecological control of sib mating when assuming Y drive as when assuming recessive lethal drive. Oscillations appeared to require even more extreme deviations from linearity in the sib mating function under Y drive

Summary
Despite some interesting quantitative differences between a recessive lethal drive and a Y drive, the qualitative results are similar in both cases: sib mating evolves, and the fitness recovery is closely tied to the fitness resulting from sib mating.

Dynamics of all cases
Results described above are for long-term, equilibrium behavior. Short-term dynamics are of interest to understand how and how quickly equilibrium is attained. Figure 7 shows representative dynamics from a single trial of each of the four classes of models. The most significant result is that genetic evolution of resistance experiences a large crash in mean fitness when the drive initially sweeps. This crash is due to the low initial frequency of the sibmating 'A' allele, and the rebound in mean fitness is rapid. Nonetheless, the nadir in mean fitness is a vulnerability to extinction. Populations with the environmental inbreeding function do not experience this crash, but mean fitness remains low relative to recoveries under genetic control of sib mating.

DISCUSSION
The intentional engineering of gene drives has become so feasible that this intervention can be entertained for nearly any sexual plant or animal species. The two most basic possible uses of a genedrive system are to drag a genetic cargo through a population (with potentially little fitness consequence) or to suppress population reproduction, possibly to extinction. Extinction is the most profound and far-reaching of these applications, and it is also the most likely to select resistance.
Our interest was in specific ways an extinction-causing gene drive might fail. The widespread enthusiasm for gene-drive implementations to suppress populations (see the Introduction) suggests that there is not a general appreciation for the possibility of failure. Our perspective here was the specific one of population structure and evolution. We analysed mathematical models of evolutionary and demographic processes that deviate from random mating to understand how sensitive might be the outcome of extinction to violations of the commonly assumed random mating of populations. Two broad classes of models were analysed. In one type, the population experienced only demographic effects of the gene drive, and gene-drive evolution was affected by those demographic effects. In the other type, genetic variation for a particular type of resistance was introduced-sib mating. The models complement each other: all address the evolution of gene drives under various forms of structured populations. The biological differences among the processes modeled are large, but all outcomes support a common conclusion that structure thwarts the evolution of an extinction-causing gene drive.
Resistance evolution to block an intentional gene-drive release will nearly always be undesirable, except perhaps in limiting the drive's spread beyond the intended species. In contrast to demographic effects of a lethal gene drive, resistance evolution leaves a permanent genetic mark on the population. Even so, some types of resistance evolution will be more undesirable than others, as some types of resistance will block future efforts that use the same or related technology.
Three classes of resistance evolution can be anticipated: (R1) Resistance blocks the mechanism or action of the drive; with CRISPR, this resistance could involve changes in the nuclease target sequence, could block CRISPR expression, or interfere with the CRISPR RNA/protein complex.
(R2) Resistance may not interfere with the drive but merely compensate for its effects. Lyttle [25] observed a change in Drosophila sex determination that tolerated the Y chromosome in both sexes and thus blocked the effect of a driving Y chromosome without blocking the drive. Burt [6] noted that the effects of a recessive lethal drive could be abrogated by compensatory evolution in a different gene that assumed the function of the targeted gene.
(R3) As studied here, resistance may alter the mating structure of the population, protecting subsets of individuals from invasion by the driving element.
Type R1 underlies many mathematical models of resistance evolution, in that the resistant allele segregates opposite the drive allele [6,26,27]. Type R1 resistance seemed like an inevitable outcome of lethal gene-drive efforts using CRISPR [19,28], although the deployment of multiple guide RNAs to several targets at once was a possible bypass [6,28,29]. However, a recent study identified a target sequence in mosquitoes that is both essential and apparently intolerant of change: caged populations of several hundred mosquitoes did not evolve resistance, instead going extinct [14]. Mutations in the target sequence were observed at some life cycle stages, likely a consequence of imperfect repair of DNA lesions, but they were incompatible with fertility and thus did not evolve.
The potential for resistance-and types of resistance-will depend on the design of a drive system, as shown at least by Kyrou et al. [14], Champer et al. [28] and Oberhofer et al. [29]. The choice of target sequences is obviously important. Some expression constructs may be more easily blocked than others. In addition, drives with cargo may be prone to lose the cargo, thereby introducing a secondary drive that competes with the original design but fails to provide the desired properties. Drives may even be engineered in various ways to prey on other drives, creating a type of arms race that will prevent any one from spreading throughout the population [30].
The type of resistance that evolves has consequences beyond the immediate gene-drive implementation. Specifically, there are different degrees to which resistance will be independent of the molecular mechanism of a gene drive. If resistance is a change in target sequence, that resistance will not affect drives that use other target sites. In contrast, resistance that is somewhat independent of the drive implementation may block future gene-drive implementations that target other sites. Thus, resistance in the form of inbreeding and other changes in mating structure will block population-suppressing drives regardless of the technology. Inbreeding that is only partial may remain effective in limiting future population-suppressing drives but still allow the spread of harmless, cargo-carrying drives. The downside of any evolution of resistance that prevents extinction, even one in which mean fitness remains low, is that it provides a population nucleus in which further evolution of resistance may occur. Following evolution of inbreeding that even partially rescues, subsequent evolution could be of reduced inbreeding depression [31][32][33] or could be other classes of resistance. Being non-essential, CRISPR may be especially prone to suppression in the long term. Rapid extinction may be the only hope for persistent suppression [19]. The speed with which gene drives evolve may in fact greatly facilitate their success in extinction, however. For example, it seems unlikely that gradual evolution of inbreeding would occur in response to a gene drive if the initial evolutionary response did not attain an high enough level of inbreeding to prevent the population from declining.
The results here indicate that various types of population structure have the potential to thwart a gene-drive implementation. We first addressed this with a simple version of the haystack model that invoked implicit spatial structure. We then extended this to structure at the family level in the form of inbreeding, and showed that evolution of sib mating is a threat to an extinction-causing gene drive. Evolution of inbreeding is not assured, and its evolution is far more complicated than is evolution of a change in the target sequence. Even if inbreeding does evolve, the recovery of mean fitness is approximately limited by the magnitude of inbreeding depression. However, other forms of inbreeding are theoretically possible [34], and some of these may be less affected by inbreeding depression than is sib mating.
Evolution of sib mating was robust to the two types of drivebased extinction processes we considered: recessive lethal and Y chromosome. Is sib mating, if it can evolve, also likely to block other types of drive-based extinction? In particular, Kyrou et al. [14] developed a system in which drive operated in both sexes, but only females were rendered infertile; males were normal. From the perspective of sib mating, their system could have parallels to either recessive lethal drive or Y drive. If sterile females behave normally, the system would be similar to recessive lethal drive in that sib mating in a brood with sterile sisters would lead to no progeny. If sterile females were unattractive (abnormal sex phenotype), then the males would not be prone to mate their (sterile) sisters, similar to Y drive. In either case, we expect that sib mating would be favored: it creates lineages that are protected from invasion by the gene drive.

Anomalies in the theory
A previous theoretical analysis of the evolution of selfing in response to a lethal gene drive observed that selfing was favored only if the selfing allele enacted a sufficiently large degree of selfing [18]. That was an encouraging result in suggesting that selfing could not evolve gradually. But that result was demonstrated for models of a recessive lethal drive in which drive was limited to one sex (and was not evaluated for models of drive in both sexes). None of the results here support that outcome.
Further numerical studies of those selfing models (conducted for comparison with the current study; data not shown) suggest that the block to gradual evolution of selfing (the most extreme form of inbreeding) is due to a combination of (i) recessive lethal drive, and (ii) incomplete drive (e.g. drive operating only in one sex, or the segregation distortion of the drive being less than 100%). An intuitive argument for the contrast between the previously observed block to gradual evolution of selfing and for the permitted gradual evolution of sib mating observed here is as follows: when a drive distortion is complete, any mating in which one parent carries the drive allele produces a family in which every offspring carries the drive allele. The alleles for inbreeding in this family never return to the pool of genotypes that lack drive. But when drive is incomplete, a mating in which one parent carries the drive allele sometimes produces descendants that do not carry the drive allele, and they or their descendants can then return to the drive-free pool and influence the frequency of alleles for inbreeding. In such families, an allele that imposes a low level of inbreeding will return more progeny to the drive-free pool than does an allele that imposes a high level of inbreeding-because the drive allele is a recessive lethal. Thus, an incomplete drive partly selects against inbreeding by disproportionately returning low-inbreeding alleles back to the drive-free pool of genotypes.
We conjecture that similar arguments apply to some of the differences in evolutionary dynamics observed here between Y drive and recessive lethal drive (e.g. comparing Figs 3 and 6). In particular, the mating between a parent carrying a recessive lethal drive allele and an allele for high inbreeding kills grandchildren from the inbreeding, whereas there is no such killing of grandchildren with Y drive. This asymmetry may explain why the sib-mating allele always fixes with Y drive but not always with recessive-lethal drive.

Extinction
The models in this study are strictly of population genetics and do not account for population size. Yet the motivation for this work is to understand resistance that blocks extinction. Inferring extinction from mean fitness requires insight to ecology-fecundity, the nature of population regulation, Allee effects at low density, and so on. Predicting extinction will thus face many challenges beyond merely predicting resistance evolution. These difficulties are foreshadowed by a 1977 symposium on the evolution of resistance to the US implementation of the sterile insect technique against the screw worm [35]. The program had suffered a recent rebound in screw worm cases, and various forms of resistance evolution were entertained to anticipate a long term failure of the program. In the final analysis, the rebound was apparently due chiefly to factory evolution of the strains used for release; replacement of those strains in rearing facilities was adequate to restore program efficacy and led to ultimate eradication from North and Central America. Resistance evolution did not prevent program success despite the plausibility of various possible avenues of failure.
Yet even if predicting resistance evolution proves elusive, an understanding of how resistance might evolve makes it possible to take steps to identify the best candidate species for release and to ensure that failures of the first efforts do not thwart later efforts. Experience with the sterile insect technique led to the realization that some species characteristics were more prone to success than others [36,37]. We can likewise suggest a few characteristics that should facilitate extinction by gene drive: (i) High inbreeding depression. If sib mating evolves, the fitness recovery is limited by inbreeding depression. Species with high inbreeding depression, and also for which inbreeding depression is slow to reverse on extended inbreeding, should face difficulty in escaping extinction through inbreeding evolution. Magnitudes of inbreeding depression are easily studied experimentally with any species that can be housed artificially, so species anticipated as targets of extinction-causing gene drives might be screened in advance to decide on the feasibility of inbreeding evolution. (ii) Low fecundity. Declines in mean fitness have a greater impact in suppressing numbers of individuals when females produce few offspring than when they produce many. The nature of population regulation also enters into this effect. The tse-tse fly appears to be an ideal candidate in this respect [38]. Small population size or low density. Opportunities for resistance mutations should be fewer in small populations, and mating opportunities of any resistant individuals should also decline with population density. This principle was an important one behind success of the sterile insect technique. (iv) Intrinsic outbreeding. Some life histories may be intrinsically disposed to outbreeding and thus face difficulty in evolving inbreeding.
Whether extinction-causing gene drives will commonly avoid resistance evolution remains to be seen. Nor is it clear that experience with the sterile insect technique will translate to gene-drive extinctions: overwhelming a wild population with sterile individuals will have different consequences for the evolution of resistance than will suppression of population densities that creates a paucity of opportunities for mating. From this perspective, a Y-chromosome drive may create a demography of extinction more similar to the sterile insect technique than does a recessive lethal drive. However, and as foreshadowed by the sterile insect technique, it is expected that some gene-drive extinction efforts will succeed and others will fail. We may at least hope that we can develop a sense for the difference and understand how to improve the chances of success.
Finally, we comment about the fact that our model does not include population demography. Any version of the model with demography must include some sort of density regulation, effectively making brood size a dynamic quantity. This would take us well beyond the scope of this article, but is certainly relevant for a more complete investigation of extinction. Gene-drive dynamics will likely be very sensitive to the type of density regulation chosen, and so it is difficult to infer general principles of the type we seek here. Female haplotypes: ad; Ad; aD; AD; Male haplotypes: ad; Ad; aD; AD (Any mating between a D male and a D female will be nonviable. All other matings produce offspring.) Maximum brood size per mating: b (Each outcrossed female produces a brood of b offspring, half of which are female. Each sibmated female produces a brood of b offspring, half of which are female.).
Mating frequencies for parents in generation t: When denoting mating pairs below, the first subscripted haplotype is for the female and the second is for the male. Quantities (appropriately subscripted) of the form S(t) and O(t) are fractions corresponding to (viable) mating pairs formed by the parents in generation t. These fractions add to 1. The offspring from these generation t matings have their future (generation t + 1) mating numbers encapsulated in (unnormalized) quantities of the form S 0 ðtÞ and O 0 ðtÞ. The sum of all these 'primed' quantities is the total, T(t), number of female offspring from generation t (the potential mothers of generation t + 1). Not all of these matings produce viable offspring. Dividing the S 0 ðtÞ and O 0 ðtÞ terms by T(t) produces the normalized quantities Sðt þ 1Þ and Oðt þ 1Þ of the next generation. The mean fitness in generation t is defined to be wðtÞ ¼ TðtÞ. For example, a fraction S Ad;ad ðtÞ of the adult females in generation t have haplotype Ad and mate with an adult male (brother) with haplotype ad. These two first appeared as offspring in generation t À 1. Note that we do not list potential mating pairs with both parents carrying the D allele (e.g. S aD;aD ðtÞ and O AD;aD ðtÞ) since they produce no offspring. Note that we track actual numbers (densities) of female types and only fractions of male types. This is because each female is allowed to reproduce only once per generation but males can reproduce multiple times within a generation. Since offspring numbers are limited by the number of females (and not by the number of males), it is natural to track mean fitness in terms of females.
The offspring produced in generation t will be the adults of generation t + 1. Offspring and parents (as well as any unmated males) are assumed to coexist for a brief time in a given generation, but only offspring transition to the next generation (becoming the new adults); i.e. all individuals have a life span of one generation. To compute the mating quantities for the next generation, we first record the numbers of female offspring haplotypes and their anticipated mating types.
Future sib mated families (one per female; before normalization): Note that we could have included the nonviable mating pairs that can form, as long as we used brood size 0 Á b instead of b (for viable outcrossed matings) or Á b (for viable sib matings). Remark: The above equations are retrospective in the sense that quantifying the new mating types in Equations (5) and (6) involves frequencies of mating types for females and males in the parent generation. To see where the different terms come from, it is instructive to think prospectively for a moment. For example, a fraction S Ad;ad ðtÞ of adult females in generation t have haplotype Ad and sib mate with males (brothers) having haplotype ad. Since these are sib matings, they each produce a brood of size b, half of which are female. Of these female offspring, all have allele d at the d/D locus (since both parents carried the d allele), and half have allele a at the a/A locus. Similar fractions apply for the male offspring. Of the b=2 female offspring from each such mating, a fraction s A will mate with a brother and a fraction 1 À s A will outcross. (The mother's haplotype at the a/A locus determines the fraction of her offspring that sib mate.) So, e.g. each of the featured matings will produce b=2 Â s A Â ð1=2Þ Â ð1=2Þ ¼ bs A =8 future matings to S 0 ad;ad ðtÞ. The total serves two purposes. Since parental mating pair frequencies are normalized, T(t) is the mean number of female offspring per female adult in generation t. This is mean fitness in generation t: wðtÞ ¼ TðtÞ: T(t) also serves as the normalizing quantity used to turn the raw counts of (anticipated) matings by generation t offspring into the normalized family frequencies at time t+1.
Updating quantities for generation t + 1 Normalized family frequencies at time t + 1: There are four equations of each of the following types Notice that the maximum brood size, b, cancels in the above ratios. It does appear in mean fitness, though our plots of mean fitness were generated with b = 1; i.e. plots of mean fitness are normalized by brood size. In fact, the effect of brood size is transient in this model since we normalize down to mating frequencies each generation. Y DRIVE EQUATIONS Figure 5 provides a schematic to suggest some (but not all) model components.
Female haplotypes: aX, AX; Male haplotypes: aY ; AY ; aZ ; AZ (Z is the Y drive allele, carried only in males. Any mating involving a Z male results in all Z male offspring; no daughters. All matings involving a Y male results in a brood of half daughters and half Y sons. There can be no sib mating involving Z males since they have no sisters.).
Brood size is as in the recessive lethal model. Mating frequencies for parents in generation t: When denoting mating pairs below, the first subscripted haplotype is for the female and the second is for the male. Other notations parallel those in the recessive lethal drive model description.
Fractions of sib mating pairs: