Mathematical modeling of genetic pest management through female-specific lethality: Is one locus better than two?

Many novel genetic approaches are under development to combat insect pests. One genetic strategy aims to suppress or locally eliminate a species through large, repeated releases of genetically engineered strains that render female offspring unviable under field conditions. Strains with this female-killing characteristic have been developed either with all of the molecular components in a single construct or with the components in two constructs inserted at independently assorting loci. Strains with two constructs are typically considered to be only of value as research tools and for producing solely male offspring in rearing factories which are subsequently sterilized by radiation before release. A concern with the two-construct strains is that once released, the two constructs would become separated and therefore non-functional. The only female-killing strains that have been released in the field without sterilization are single-construct strains. Here, we use a population genetics model with density dependence to evaluate the relative effectiveness of female-killing approaches based on single- and two-construct arrangements. We find that, in general, the single-construct arrangement results in slightly faster population suppression, but the two-construct arrangement can eventually cause stronger suppression and cause local elimination with a smaller release size. Based on our results, there is no a priori reason that males carrying two independently segregating constructs need to be sterilized prior to release. In some cases, a fertile release would be more efficient for population suppression.

INTRODUCTION

Insect pests remain a burden to human health and agriculture (Deutsch et al., 2018; World Health Organization, 2017). Genetic pest management aims to reduce this burden by releasing engineered insects that either introduce a desired trait into a natural population or reduce the size of the population. There have historically been several large area‐wide inundative releases of male insects that were rendered sterile by exposure to radiation (Gould & Schliekelman, 2004). In these releases, local elimination of the target species was achieved as females increasingly mated with the sterile males rather than the wild‐type males with whom they would produce viable offspring. Instead of using radiation to cause sterility, a contemporary alternative is to genetically engineer strains in which the males cause all their offspring or exclusively their daughters to die or to have low fitness (Alphey, 2002). Genetic pest management through such approaches is an area of active research, and genetically engineered strains in a number of species have been cage‐ or field‐tested (Ant et al., 2012; Carvalho et al., 2015; Gorman et al., 2016; Harris et al., 2012; Harvey‐Samuel et al., 2015; Lacroix et al., 2012; Leftwich et al., 2014; Wise de Valdez et al., 2011).

One approach to developing these functionally sterile strains involves inserting a repressible, dominant lethal trait, which can be active in both sexes or in females only (Heinrich & Scott, 2000; Thomas et al., 2000). For either the female killing (FK, also sometimes referred to as fsRIDL, or female‐specific release of insects carrying dominant lethals) or bisex killing (BK), in order to rear the transgenic strain in the generations prior to release, it must be possible to inactivate the dominant lethal gene. Often this is achieved through a Tet‐off system where tetracycline in the diet represses the activator for a lethal gene (Gossen & Bujard, 1992). For an FK strain, the release generation is reared on a diet not containing tetracycline. This results in only males surviving. Further, as the offspring of released FK or BK males would feed on a tetracycline‐free diet under field conditions, the lethal gene is turned on and death ensues.

Intuitively, modeling studies have found that FK can be advantageous over BK because it kills females while allowing the transgene to propagate through multiple generations in heterozygote males (Schliekelman & Gould, 2000; Thomas et al., 2000). This would seem especially useful when females but not males transmit pathogens. However, heterozygous males can also serve as a reservoir for wild‐type alleles, which can make FK less effective than BK under some conditions (Foster et al., 1988; Gentile et al., 2015). It should be noted that BK strains for mosquito disease vectors typically require sex sorting because release of females would be considered unacceptable. It can also be advantageous to release only males as females do not contribute to genetic suppression and tend to mate with the released males and thus reduce their efficiency (Rendón et al., 2004) except in some situations where there is age structuring in the population (Huang et al., 2009).

The full molecular design involves two components: (1) the tetracycline‐repressible transactivator (tTA) with a promoter; and (2) a lethal gene with an enhancer/promoter consisting of multiple tTA binding sites (tetO) and a core promoter. In the initial two‐component systems, tTA was expressed in females by using a female‐specific promoter (Heinrich & Scott, 2000; Thomas et al., 2000). The second component was a lethal gene (e.g., proapoptotic) driven by a tetO enhancer–promoter. The two‐molecular components were built in separate constructs that were inserted independently. Subsequently, a simpler, two‐component system was developed in which tTA acts as both the activator and lethal gene. Here, a single construct includes a tTA coding sequence driven by a tetO enhancer–promoter. In this autoregulated system, high levels of the tTA activator cause lethality in late‐stage larvae or in pupae. The mortality is possibly due to a general interference with transcription (Gong et al., 2005). FK single‐construct strains have included a sex‐specifically spliced intron from the transformer or doublesex genes inserted within the tTA gene (Fu et al., 2007). In these FK strains, only the female tTA transcript encodes a functional protein. A different, single‐construct approach for FK with Aedes aegypti and Aedes albopictus uses a female‐specific indirect flight muscle promoter from the Actin‐4 gene (Fu et al., 2010; Labbé et al., 2012). All field trials with transgenic FK or BK strains have been with single‐construct strains.

More recently, two‐construct FK strains have been made with an early embryo promoter driving tTA expression and a tTA‐regulated lethal gene that contains a sex‐specifically spliced intron (Yan et al., 2020). An advantage of these strains is that female lethality occurs at the embryo or early larval stages, which produces considerable savings in larval diet costs in a mass‐rearing facility, though there can also be a disadvantage in the field compared with lethality at the late larval or pupal stages because transgenic juveniles die without providing competition to viable juveniles. Increased competition with late‐acting lethality results in fewer viable juveniles surviving to adulthood through density‐dependent mortality and thus better suppression (Gentile et al., 2015; Phuc et al., 2007). Although it should be possible to develop any two‐component system as a single construct (Yan & Scott, 2015), they are typically developed as independently segregating constructs. Germline transformation in insects is often achieved through the use of transposable elements such as piggyBac, and due to the randomness of the insertion process, a large number of injections can be required to obtain transgenic strains (Gregory et al., 2016). Furthermore, there are often multiple potential choices for one or both of the components. For this reason, it can be advantageous to separately produce strains with different promoters and lethal genes, then produce individuals bearing both components by crossing to test effectiveness of different combinations. The final transgenic insects have the two components located at two, separate loci (Ogaugwu et al., 2013; Schetelig & Handler, 2012; Schetelig et al., 2016; Scott, 2014; Yan et al., 2020).

Female killing strains with two constructs are generally thought of as useful research tools with potential to be used in rearing facilities so that the final generation before release would only produce males (Ogaugwu et al., 2013; Schetelig & Handler, 2012; Yan & Scott, 2015). It has been suggested that independent inheritance of the components would cause a breakdown in the female killing in the second generation after release (Ogaugwu et al., 2013; Yan & Scott, 2015). However, previous theoretical studies of FK systems have only modeled the components as being inserted together on a single locus (Alphey et al., 2011; Gentile et al., 2015; Schliekelman & Gould, 2000; Thomas et al., 2000) so a comparison of the two approaches has been lacking.

Here, we evaluate the effectiveness of 1‐ and 2‐locus FK, along with BK for comparison. We use a mathematical model parameterized for the Ae. aegypti mosquito that is a vector for several human pathogens. We explore the release of strains with killing in either juveniles or adults. We show that under reasonable assumptions about fitness costs of the insertions, there is not a substantial difference between the 1‐ and 2‐locus FK approaches, particularly when compared to the differences between FK and BK. These results demonstrate the release potential of recently developed 2‐locus FK constructs.

METHODS

Our mathematical model implements the genetics of FK and BS by separately tracking the number of individuals in the population of each genotype, with genotype denoted by subscript i. For the single‐locus system (Table 1), we let the transgenic allele be represented by K and the wild‐type allele at that locus be represented by k, with a total of N = 3 possible genotypes. For the 2‐locus system (Table 2), we let A and B represent the transgenic alleles (i.e., tTA and lethal gene) inserted at two separate loci with wild‐type alleles a and b, respectively, for a total of N = 9 possible diploid genotypes.

TABLE 1

1‐locus genotypes, with associated viabilities and fitnesses

i	Genotype	γia,S	wix
1	kk	1	1
2	Kk	0	1‐h·s^x
3	KK	0	1‐s^x

Viability of genotype i, , takes the value listed when the approach causes loss of viability in sex S (for female killing, only when S = F; for bisex, when S = F or S = M) with timing a (for early approaches, when a = E; for late approaches, when a = L), and are 1 otherwise. Fitnesses apply for both hatching (x = H) and male mating competitiveness (x = M).

TABLE 2

2‐locus genotypes, with associated viabilities and fitnesses

i	Genotype	γia,S	wix
1	Aabb	1	1
2	aaBb	1	1‐h·s^x1‐c_A
3	aaBB	1	1‐s^x(1‐c_A)
4	Aabb	1	1‐h·s^xc_A
5	AaBb	0	1‐h·s^x
6	AaBB	0	1‐h·s^xc_A‐s^x(1‐c_A)
7	AAbb	1	1‐s^xc_A
8	AABb	0	1‐s^xc_A‐h·s^x1‐c_A
9	AABB	0	1‐s^x

Viability of genotype i, , takes the value listed when the approach causes loss of viability in sex S (for female killing, only when S = F; for bisex, when S = F or S = M) with timing a (for early approaches, when a = E; for late approaches, when a = L), and are 1 otherwise. Fitnesses apply for both hatching (x = H) and male mating competitiveness (x = M).

We assume complete effectiveness of the constructs, so when there is no gene repression via tetracycline, all individuals bearing the functional BK system and all females with the functional FK system die, with a genotype viability of zero (Fu et al., 2007; Ogaugwu et al., 2013; Yan et al., 2020). One copy of K is assumed to be sufficient to induce lethality in the 1‐locus system, and only one copy each of A and B is required in the 2‐locus system. We consider lethality acting at different points in the life cycle. In insects that experience strong resource competition during larval stages, having the transgene‐induced mortality occur during or shortly after the pupal stage, instead of during the egg or larval stages, can yield stronger population suppression. This is because the transgenic juveniles consume resources and therefore increase wild‐type juvenile mortality. We model early mortality (E) as occurring in the embryo and late mortality (L) as occurring in pupal stages or in adults before mating, and we assume these differentiate whether the individual contributes toward density‐dependent mortality of all individuals in the population. We let and represent the early (embryonic) and late (adult) expected viabilities for individuals of sex S and genotype i. Tables 1 and 2 give expected viabilities for individuals with each construct and genotype.

We classify constructs into four different approaches depending on when the dominant lethal gene is active, similar to Gentile et al. (2015): early bisex killing (E‐BK), late bisex killing (L‐BK), early female killing (E‐FK), and late female killing (L‐FK). We assume male transgenic homozygotes are released, so mating with wild‐type females will produce offspring that are entirely heterozygous, with a copy of each transgene. If the construct(s) affects both sexes (BK), none of these offspring will survive to mate and pass on their genes, making bisex 1‐locus and 2‐locus equivalent in terms of both population genetics and population dynamics. Female‐specific approaches (FK) allow males to continue to propagate the transgenes, and thus, inheritance differs between 1‐locus and 2‐locus approaches. In all, we consider the following six approaches: E‐BK, L‐BK, 1‐locus E‐FK (E‐FK1), 2‐locus E‐FK (E‐FK2), 1‐locus L‐FK (L‐FK1), and 2‐locus L‐FK (L‐FK2).

Separate from the transgenic, toxin‐induced lethality, we account for potential fitness costs caused by the genetic insertion itself. We allow the fitness costs of inserting a novel genetic element to manifest at an early stage as a reduction in the ability of a zygote to survive beyond the egg stage, that is, the fraction of eggs of that genotype, which survive and hatch into larvae. We let the genotype's hatching fitness, , equal the probability of successfully entering the larval stage, with wild‐type hatching fitness . We also allow for transgenic fitness costs to males in the form of reduced mating competitiveness, , as defined below, with wild‐type mating competitiveness .

We generally assume that the fitness costs are equal for the homozygotes in the 1‐locus and 2‐locus systems to facilitate a direct comparison between the two systems. The 2‐locus system has the same components as the 1‐locus system, which makes equal fitness costs a reasonable base assumption for the purposes of this work (this assumption is relaxed in Figures 3 and S1). We let and be the hatching and mating competitiveness fitness costs, respectively, to the homozygotes KK and AABB, and we allow the two types of costs to vary independently. For simplicity, we assume the degree of dominance for the fitness costs, h, is equal for hatching and mating competitiveness. Unless otherwise noted, we assume costs are additive, with h = 0.5, such that each copy of the K allele alone contributes a fitness cost of for the 1‐locus system (x here indicates that the fitness cost can either be hatching or mating). For the 2‐locus system, we allow for unequal fitness costs between each of the insertions. We let two copies of the A allele contribute a fitness cost of , where is the proportion of the total 2‐locus fitness cost accounted for by the A allele, and one copy of the A allele contributes a fitness cost of . A single B allele contributes a fitness cost of , while being homozygous for B induces a cost of . Note that A can be arbitrarily defined as whichever component is more costly because the two alleles only matter in combination. Resulting fitness expressions for all genotypes are listed in Tables 1 and 2.

FIGURE 3

Effects of fitness cost variation on release efficacy for weekly release ratio r = 7 and β = 3. Each column of panels shows the results for a different genetic approach, while each row of panels depicts a different degree of dominance, h. Within each individual panel, the hatching fitness cost, s, increases from 0 to 1 along the x‐axis, and the cost to male mating competitiveness, s, increases from 0 to 1 along the y‐axis. For every point, a deterministic simulation was run with a unique combination of genetic approach and fitness parameters, and color indicates the number of days until the number of viable adult females is under 0.05% of equilibrium. Darker colors show faster times, with a minimum time of 73 days, and lighter colors show slower times up to 500 days (chosen as a threshold to improve the ability to visually differentiate times below 500 days). White areas indicate that the number of females did not fall below the threshold within 500 days. The colored points in the middle row correspond to the times in Figure 2B at r = 7

We model genotype counts over time using a system of ordinary differential equations adapted from Robert et al. (2013). We let and be the number of juvenile (larvae and pupae) males and females, respectively, of genotype i at time t, and and be the number of viable adult male and adult female mosquitoes in the population, respectively, of genotype i at time t. This gives a maximum of 12 classes of individuals to track (each with different combinations of the three genotypes, two sexes, and two age classes) for the 1‐locus system and 36 classes for the 2‐locus system, though lethality from the genetic construct prevents survival of certain classes. For instance, E‐BK only has five nonzero classes (wild‐type male and female juveniles and adults, and male adult homozygote transgenic, which are released). We also assume that in the absence of early sex‐dependent lethality due to transgenes, males and females have equal birth rates and equal hatching fitness costs, allowing a further reduction in the number of unique classes for late‐acting approaches because for all i and t. This results in seven classes for L‐FK1 (after removing three juvenile classes and two non‐viable adult female classes) and 23 classes for L‐FK2 (after removing nine juvenile classes and four non‐viable adult female classes). These dimensionality reductions can be useful when finding analytical solutions, but for simplicity, we computationally simulated all 12 (1‐locus) or 36 (2‐locus) equations.

Accounting for fitness costs, adult females produce juveniles of genotype at time t at ratewhere is the per‐capita birth rate and is the probability that a juvenile produced from a mating between a female and male of genotypes and , respectively, will be of genotype . The fraction gives the probability that a randomly chosen male adult is of genotype n, weighted by mating competitiveness . The offspring genotype probabilities are calculated assuming Mendelian inheritance, and for the 2‐locus case, independent segregation of genes at each locus.

Juveniles of each genotype and sex emerge to adulthood at per‐capita rate . We assume juveniles, adult males, and adult females have per‐capita density‐independent mortality rates of , , and , respectively. Juveniles also undergo density‐dependent mortality at a per‐capita rate , where is the total number of juveniles at time , and α and β are parameters related to the strength of density‐dependent mortality. The strength of density dependence is adjusted by varying β, with higher β resulting in a faster return to equilibrium population size after a small perturbation. A value of β = 2 gives the logistic model for population dynamics. Given the many complexities and limited understanding of density‐dependent effects in field Ae. aegypti systems (Than et al., 2020), we consider a range of possible β values. By default, we let β = 3 to model an environment that would be more difficult for successful suppression compared with the logistic case (Hibbard et al., 2010). The equilibrium size of an entirely wild‐type population varies with α, and to keep simulations with different values of β comparable, we choose the value of α so that the equilibrium number of wild‐type females remains the same (Robert et al., 2013).

We assume a continuous release of homozygote engineered males at (daily) rate , where is the weekly release ratio (engineered:wild‐type) based on the equilibrium number of males prior to the release. By maintaining a constant number of released males, the effective release ratio increases as the population size decreases. The release genotype is KK for 1‐locus (i = 3) and AABB for 2‐locus (i = 9), and because we assume that no females are released, for all i.

The resulting system of ordinary differential equations (with time dependence of , , , and omitted for simplicity of notation) is for i = 1 N and S = F or M. All model parameters are listed in Table 3 and are based on the empirical estimates for Ae. aegypti used by Robert et al. (2013). While the rates of mortality, larval production, and emergence to adulthood apply to Ae. aegypti, the resulting population dynamics, simulated with different strengths of density dependence, would likely be similar to many other species. All numerical simulations of the ordinary differential equations have initial conditions at the wild‐type equilibrium.

TABLE 3

Model parameters

Parameter	Description	Value/range	Reference
μ_J	Density‐independent per‐capita juvenile mortality rate	0.03 day−1	Rueda et al. (1990)
μ_M	Adult male per‐capita mortality rate	0.28 day−1	Muir and Kay (1998) and Fouque et al. (2006)
μ_F	Adult female per‐capita mortality rate	0.10 day−1	Muir and Kay (1998) and Fouque et al. (2006)
λ	Female per‐capita larval production rate	8 day−1	Harrington et al. (2001) and Styer et al. (2007)
v	Per‐capita rate of emergence to adulthood	0.14 day−1	Muir and Kay (1998)
β	Strength of density dependence (see text)	2 to 4
A1F(0)	Equilibrium number of wild‐type females	2000	Saarman et al. (2017)
r	Weekly release ratio of transgenic: equilibrium wild‐type males	0 to 12
α	Density dependence parameter chosen based on β and A1F(0)	Dependent on β and A1F(0)
wiH, wiM	Egg hatching and male mating competitiveness fitnesses of genotype i	See Tables 1 and 2
s^H, s^M	Egg hatching and male mating competitiveness fitness costs to homozygotes	0 to 1 (0.2 and 0.1 default, respectively)
h	Fitness cost degree of dominance (proportion of s^x incurred in heterozygotes)	0 to 1 (0.5 default)
c_A	Proportion of s contributed by A (2‐locus)	0.5 to 1 (0.55 default)
γiE,S, γiL,S	Viabilities of genotype i for early (E; embryonic) and late (L; adult) approaches in sex S	See Tables 1 and 2

Parameters are based on empirical estimates for Ae. aegypti used by Robert et al. (2013).

In order to explore the effects of demographic stochasticity and genetic drift, we also run simulations using an analogous continuous‐time Markov chain model (for details, see Appendix S1).

RESULTS

Each of the FK (female killing) and BK (bisex killing) genetic strategies has the goal of causing the population to decline by reducing the number of reproductive adult females. The strength of density‐dependent mortality moderates the reductions in population size because stronger density dependence (higher β) causes the juvenile mortality rates to decrease more quickly as population size decreases from equilibrium. In a system with strong density dependence, the weekly release ratio (r) must be larger to achieve the same amount of population suppression as in a system with weak density dependence. Large r can result in target population extinction, which we define as a time point when there are no females left in the population, even though males are still being released. In contrast, small r results in a new, lower equilibrium population density, where the proportion of individuals that die due to bearing the transgene is not high enough to outweigh the increased survival of juveniles due to decreased density‐dependent mortality in the smaller population.

Figure 1 demonstrates the outcome of release for each genetic approach at r = 1 under the deterministic model (see Figure S1 for time series with other release ratios). At this release ratio, approaches that are late‐acting (i.e., mortality in pre‐mated adults, indicated by “L‐”) reduce the number of viable females to a lower number than approaches that are early‐acting (i.e., mortality in the embryonic life stage, indicated by “E‐”). BK approaches reduce the number of females faster and to lower levels than female‐specific approaches with late‐acting mortality, but the opposite is true for early‐acting mortality. Among FK approaches, 1‐locus reduces the number of females more quickly initially, but 2‐locus eventually suppresses the population slightly more than 1‐locus (which we explore further below). Overall, this suggests that L‐BK is most effective, followed by L‐FK, E‐FK, and E‐BK, with little difference between 1‐ and 2‐locus FK. Under stronger density dependence (higher values of β), releases result in weaker suppression of the population, with none of the approaches causing extinction of the population if β = 4.

FIGURE 1

Effect of transgenic releases on population size over time for various strengths of density dependence (β). The number of viable adult females (relative to pre‐release equilibrium) over time is plotted for deterministic simulations with adults for each genetic approach released at a continual weekly release ratio of 1:1 transgenic males to the pre‐release equilibrium wild‐type males (r = ). Line type differentiates embryo (early‐acting, dashed line) and adult (late‐acting, solid line) mortality, and line color differentiates bisex (red), 1‐locus female killing (blue), and 2‐locus female killing (green) constructs. Releases are less effective as strength of density dependence (β) increases across panels from left to right. Simulations use fitness parameters s.2, s =0.1, and c =0.55 and the remaining parameters as listed in Table 3

In general, large release ratios result in extinction (the population goes to an equilibrium size of zero), and small release ratios result in a suppressed but nonzero equilibrium population size. For most sets of parameters, there is a critical release ratio, r, above which the release is large enough to cause the population to go extinct. With such a large number of released males, population extinction is the only stable equilibrium, meaning the release will cause extinction regardless of initial population size. For ongoing release at release ratios below r, there is a nonzero stable equilibrium for the number of viable adult females, meaning that release will not push a wild‐type population to extinction. A population size of zero is also stable, but the system will approach the nonzero equilibrium unless starting from very low population sizes; that is, it is a bistable (or multistable) system with a low unstable equilibrium (as shown in Figure S2 with time series starting from multiple initial conditions). This also means that in an area without any wild‐type individuals, ongoing release below r could prevent (small‐scale) immigration of wild‐type from re‐establishing a population. Mathematical details on the shift in qualitative behavior (i.e., a bifurcation) at r, including analysis using Mathematica (Wolfram Research, 2019) and bifurcation diagrams using MatCont (Dhooge et al., 2008), can be found in Appendix S2, Figures S3,S4.

Figure 2a shows how release ratios affect stable equilibria when β = 3. As a result of the bifurcation, each of the lines is discontinuous, jumping from a nonzero equilibrium to a zero equilibrium at r. The smallest release size required to cause extinction for each approach is  = 0.86 <  = 2.52 <  = 2.86 <  = 4.74 <  = 5.40 <  = 5.97. This order of effectiveness matches that of Figure 1 and mirrors the results for single‐locus constructs in Gentile et al., 2015. However, the order of effectiveness is different at small release ratios (as indicated by crossing lines in Figure 2a); L‐BK has a higher equilibrium population size than either of the L‐FK methods. We explore this result in detail below. There is a small relative difference between the different FK approaches, with similar equilibrium sizes when r is below the critical release ratios and 1‐locus FK requiring a release ratio <15% larger than 2‐locus FK to cause extinction.

FIGURE 2

Release outcomes across different release ratios. (a) Long‐term, stable equilibria for number of viable adult females (relative to pre‐release equilibrium) for different r, found by simulating the system of differential equations until at steady state. The asterisks indicate r = 1, for which the equilibrium values correspond to the middle panel of Figure 1. Each genetic approach exhibits a bifurcation at a critical release ratio, r, indicated by a solid circle. Above the r of each approach, that approach only has an equilibrium at zero. This makes each line have a discontinuity, most clearly visible with the red line (L‐BK). Above r = 2.52, two or more approaches lead to extinction of the population and hence have equilibria at zero: This is indicated using a gray line. Note that below r, each system exhibits bistability, with both nonzero and zero count stable equilibria. Only the nonzero equilibrium is plotted, which is the equilibrium reached when starting simulations from the wild‐type equilibrium. (b) Time until the number of viable adult females is under 0.05% of equilibrium (less than 1 adult female from an equilibrium of 2000) in deterministic simulations for different r. Color and line type match that of Figure 1. Parameters are β = 3, s0.2, s = 0.1, and c = 0.55, as in Figure 1

In settings where release causes the population to go extinct, we can consider the time it takes to reach extinction after starting release (Figure 2b). Given that deterministic simulations will only approach extinction asymptotically, we use the time it takes for the number of females to reach <0.05% of the pre‐control equilibrium, which is suppression to below one adult female when starting from a pre‐release equilibrium of 2000. The times using this threshold are comparable to the average time to extinction in stochastic simulations (Figure S5). Once release ratios are high, L‐BK drops the population under 0.05% of the equilibrium faster than L‐FK methods. Also, both of the 1‐locus FK approaches are slightly faster than the respective 2‐locus FK approaches, though the differences are small for practical purposes. At r = 12, for example, the times are 86 (L‐BK), 114 (L‐FK1), 119 (L‐FK2), 133 (E‐BK), 149 (E‐FK1), and 153 (E‐FK2) days. Note, however, that at lower release ratios, 2‐locus FK systems have a lower time to extinction than 1‐locus FK systems. The results are similar for β = 2 (Figure S6).

Overall, FK1 and FK2 are quite similar. If either of the A or B alleles in 2‐locus FK becomes fixed in the population, the 2‐locus approach becomes nearly identical to 1‐locus FK, where one copy of the unfixed allele causes mortality in females. For example, if all individuals in the population already have the B allele, only one copy of the A component is additionally necessary, just as a single copy of the K allele causes mortality. In this case, the long‐term equilibrium can be identical to 1‐locus FK, though a fitness cost to the fixed allele decreases the average fitness of the entire population and makes the population size lower for 2‐locus FK than 1‐locus FK. Whether one of the 2‐locus FK alleles becomes fixed depends on the fitness costs and release ratios of the system (see Appendices S2,S3 and Figures S4,S7,S8).

An observation from Figure 2a is that L‐BK, L‐FK2, and L‐FK1 cause a similar level of suppression when the release ratio is near r = 0.8. Previous work has suggested that L‐FK allows the wild‐type allele to propagate in heterozygous males, making it less effective than L‐BK (Gentile et al., 2015). This is true at high release ratios, but not at all release ratios. Without fitness costs, the strategies are equally effective if the number of released males is equal to the number of wild‐type males in the population at equilibrium: heterozygotes carry both wild‐type and transgenic alleles, and therefore, the survival of heterozygous males does not affect allele frequency. At small release ratios, when L‐BK release results in low transgenic frequency and thus a large equilibrium population size, survival of heterozygous males would allow the transgenic allele to propagate further and increase in frequency, explaining why L‐FK has a lower equilibrium than L‐BK in this narrow window of small releases.

The main difference exhibited between FK1 and FK2 can also be explained by their propagation of the transgenic and wild‐type alleles. When the components are separated across two loci, the A and B alleles become unlinked, with some individuals only inheriting one allele or the other, while having linked components guarantees inheritance of the transgenic allele and reduces the population size more quickly initially. This also explains the slightly faster time to extinction for FK‐1 than FK‐2 at high release ratios (Figure 2b). Eventually, however, the accumulation of transgenic alleles in the 2‐locus system causes production of a higher proportion of unviable genotypes and greater population suppression (Figure 2).

Results over a wide range of fitness parameters highlight the minimal differences between FK1 and FK2. We use the time it takes to suppress the number of viable females to under 0.05% of the equilibrium when there is a high release ratio of r = 7 as a way to measure effectiveness (Figure 3). For each approach and degree of dominance (h), there is a region with low fitness costs where there is fairly little variation between the time it takes to suppress the population below the threshold. As one or both of the fitness costs increase, there is a margin of longer times separating the successful suppression region from the unsuccessful region. As evident from Figure 2, long times to reduce the population indicate that the release ratio of r = 7 is only slightly greater than r for those parameter values, and a system with even higher fitness costs has r > 7, so suppression below 0.05% of equilibrium is not achieved.

There is little difference between FK1 and FK2, particularly for L‐FK, since the release ratio is much higher than r in most of the region with successful release. Even with different fitness costs, FK1 and FK2 do not differ drastically in time until extinction; FK2 would be substantially less effective only if the costs become large enough for the release to become too small to cause extinction (i.e., r becomes smaller than r because of the higher costs). Unpacking the differences between all approaches, for the E‐BK approach, neither h nor s affect the length of time for reduction; only released males experience costs so the degree of dominance does not affect results, and all offspring that inherit a transgene are inviable so hatching costs do not affect results. For E‐BK, increasing the cost to male mating competitiveness drastically decreases release efficacy because doing so effectively reduces the release ratio (fewer of the released males successfully mate and contribute to population reduction), and with no surviving transgenic offspring, the release ratio is the main contributing factor to E‐BK efficacy. When s = 1 and h = 1, none of the offspring survive for the E‐FK1 and E‐FK2 approaches, making the times shown identical to E‐BK. Decreasing these parameters makes a difference for FK approaches because the males are subject to fitness costs. Dominant fitness costs (bottom row) actually enable successful population reduction for higher male mating costs. With large, dominant mating costs, FK approaches have few mating adult males, becoming effectively similar to the highly effective L‐BK, which has no adult males. The full explanation for increased effectiveness with higher degrees of dominance relates to the propagation of wild‐type alleles, similar to the previous description: homozygote males that are being released are affected regardless of h, and when h > 0, the fitness cost also reduces the spread of wild‐type alleles in heterozygotes. This effect is also evident in other transgenic systems, such as two‐locus underdominance, where dominant transgenic allele fitness costs prevent the wild‐type from being maintained in the population at small frequencies (Dhole et al., 2018).

DISCUSSION

The recent literature on FK systems makes the assumption that strains built with constructs inserted at two independent loci will not be as useful for field releases as those built with a single construct (Ogaugwu et al., 2013; Yan & Scott, 2015). The assumption is that the two constructs will separate from each other in the second generation after a release and will become non‐functional. Our modeling results demonstrate that a 2‐locus FK (FK2) should behave similar to a 1‐locus FK (FK1) and would not present any significant disadvantages in its ability to suppress a population. We generally made the assumption that the 2‐locus and 1‐locus approaches would have similar total fitness costs because they have the same components. If the total cost of either the 1‐locus or the 2‐locus approach was lower than the other, that approach would likely be preferred. Importantly, based on our results, there is no a priori, general reason for genetic engineers to favor a 1‐locus system. The choice will likely depend on specific biological and genetic characteristics of the target species.

Assuming equal costs, FK1 is slightly faster at initial population reduction, but FK2 can eventually suppress the population to lower numbers. FK2 also has a slightly lower critical release ratio than FK1, meaning a smaller release size is necessary to guarantee extinction. For many combinations of fitness costs and release ratios, one of the FK2 alleles would be driven to fixation, resulting in a genetic system similar to FK1. The differences between FK1 and FK2 are much smaller than between FK and BK approaches. Comparing FK and BK approaches, our results are generally similar to previous work (Gentile et al., 2015). Late‐acting approaches cause extinction with a lower release ratio than early‐acting approaches, with L‐BK causing extinction with a lower release ratio than L‐FK, and E‐FK causing extinction with a lower release ratio than E‐BK.

While our modeling results indicate that L‐BK outperforms the other methods, there are other considerations that affect which approach may be best suited for a given scenario. In our model, parameterized for mosquitoes, density‐dependent mortality during early life stages was an important factor and caused early‐acting approaches to result in less population reduction than late‐acting approaches. In species with little density‐dependent dynamics in juveniles, the difference in effectiveness between early and late acting would be minor, though this is not the case for many pest species. Also, as we demonstrate, the impact of partial suppression depends on the nature and strength of density dependence. With overcompensatory density dependence, this could lead to the population size exceeding the pre‐control equilibrium (Alphey & Bonsall, 2014; Rajagopalan et al., 1977).

Beyond population dynamics, there are economic and social factors that differ between approaches. For some systems, it will be necessary to engineer constructs into laboratory strains and then backcross the construct or constructs into a strain that have a genetic makeup similar to the targeted population. In general, it should be easier to do the backcrossing with a one‐locus system. Rearing costs are also expected to vary between approaches. With E‐FK, juvenile females experience mortality before consuming food, whereas E‐BK, L‐BK, and L‐FK require rearing of juveniles of both sexes. Furthermore, BK approaches typically require sexing to remove females prior to release, which increases the total rearing costs and is often difficult to do with complete accuracy. When releasing a species that is a disease vector, sexing accuracy is meaningful from a social perspective as release of females could contribute to disease transmission. In species with little density‐dependent dynamics in juveniles, the difference in effectiveness between early and late acting would be minor, though this is not the case for many pest species.

Apart from engineering and rearing, for most agricultural pests, the juvenile stages of males and females cause damage to crops and livestock. In the first generations of transgenic pest releases, the late‐acting approaches will leave feeding immatures in the environment, and E‐FK will result in male immatures that still cause damage. This may not be favored by farmers even though the overall population could be decreasing rapidly, and E‐BK could be preferred. Finally, even if late‐acting mortality may be ideal for a given scenario, controlling the timing of mortality at the intended life stage may not always be feasible, for example due to leaky expression of the lethal gene.

The model used here has several limitations. An important factor that could affect population genetics is spatial heterogeneity. For example, in a spatial model of FK2, it would be possible, particularly in small populations, for different patches to have different transgenic alleles reach fixation. A spatial model would also be useful to determine whether FK2 has any differences in resilience to wild‐type reinvasion. The details of such a spatial model, including rates of release, would depend on species. A species‐specific model could also implement different forms of density dependence, age structure, mating parameters, and release patterns based on feasibility (instead of assuming continuous release as we do here). Additionally, our model does not account for an Allee effect, where there is a critical population size below which the population is unable to maintain itself. With an Allee effect, suppression does not need to be complete to lead to extinction of the population; this may benefit 1‐locus and 2‐locus approaches differently. Finally, given its generality, our model does not account for any potential mechanisms for resistance development. Depending on the mechanism of lethality, there may be advantages for having both components for lethality inserted together. While these areas require further investigation, our results indicate that overall, there is little difference in the pest population suppression efficacies of 1‐locus FK and 2‐locus FK.

CONFLICT OF INTEREST

None declared.

Supplementary Material

Click here for additional data file.