Introduction

Insects are often involved in symbiotic relationships with endosymbiotic microorganisms. This symbiosis is often a mutualistic interaction, with the host providing a habitat for the endosymbiont and the endosymbiont providing a fitness benefit to the host. Examples of the benefit include increased thermal tolerance, nutritional benefit and fecundity increases (Bennett and Moran 2015; Hurst 1997; Fast et al. 2011; Baumann 2005). Insect endosymbionts are often primarily or exclusively vertically transmitted, via infected females. Vertical transmission is typically associated with a mutualistic association (Ewald 1987; Bull et al. 1991; Bull and Molineux 1992; Herre 1995; Fisher et al. 2017; Messenger et al. 1999; Sachs and Wilcox 2006). However, vertical transmission is not limited to mutualistic interactions. Parasitic insect-endosymbiont symbioses occur, notably causing son-killing or sexual incompatibility (Hurst 1997; Engelstädter and Hurst 2009). Whether endosymbiotic relationships are ever truly mutualistic is also a point of debate, as endosymbionts could be considered to be trapped and exploited by their hosts, rather than being engaged in a truly mutualistic symbiosis (Keeling and McCutcheon 2017).

There are several examples where a bacterial species living as an insect endosymbiont has evolved to become an insect-vectored plant pathogen. Examples include: various Candidatus Liberibacter spp., the causal organisms of citrus greening disease in citrus, vectored by the psyllid Diaphorina citri (Gottwald 2010; Gottwald and Mccollum 2017); Arsenophonus phytopathogenicus, the causal organism of the syndrome “basses richesses” in sugar beet, vectored by Pentastiridius leporinus (Bressan 2014); and Candidatus Phlomobacter fragariae, the causal organism of marginal chlorosis disease of strawberry plants, vectored by Cixius wagneri (Danet et al. 2003).

Insects often act as vectors for diseases in animals and plants, in both unmanaged and agricultural systems (Burnet 1972). Notable examples of diseases often vectored by insects include phytoplasmas, phloem-limited bacterial pathogens of plants (Weintraub and Beanland 2006), bacteria (Bressan 2014; Gottwald and Mccollum 2017), viruses (Dietzgen et al. 2016), and fungi (Harwood et al. 2011). The insect-vectored plant pathogen lifestyle is a highly specialised one. The pathogen must be able to survive in both insect and plant host and must be able to successfully transmit between the two.

Arthropod-vectored viral pathogens can be categorised as non-persistent, semi-persistent or persistent. The differences between categories are found in the length of the transmission window and the area in the insect within which the viral infection is contained. Persistent viruses are further subcategorised into circulative propagative and circulative non-propagative, depending on whether the virus can replicate in their host or not (Dietzgen et al. 2016). A reduction in the genome size of the pathogen is often seen in insect-vectored plant pathogens, leading to host dependence (Perilla-Henao and Casteel 2016). Bacterial insect-vectored plant pathogens can similarly be classified into circulative propagative and non-circulative propagative (Orlovskis et al. 2015).

This transition in lifestyle from an insect endosymbiont to a plant pathogen is curious. It begs the question of why insect endosymbionts would evolve to give up a more certain and long-established transmission method. The transmission route of an endosymbiont is thought to determine the nature of its relationships with its hosts (Ewald 1987; Bull et al. 1991; Bull and Molineux 1992; Herre 1995; Fisher et al. 2017; Messenger et al. 1999; Sachs and Wilcox 2006), although this has been disputed (Ebert 2013). The argument behind this is based on the following reasoning: if an endosymbiont’s primary method of transmission is vertical, then it is heavily dependent on its host survival for its own survival and that of its future generations. Following this, it would be advantageous for the endosymbiont to positively influence the hosts’ fitness, indirectly increasing its own fitness. Contrary to this, if the endosymbiont is horizontally transmitted, then the need to keep an individual host fit is not as strong: with other potential hosts, coupled with the means to access them, the fate of an individual host could be inconsequential to the endosymbiont and its future generations, depending on the population dynamics and spatial distribution.

In this paper, we:

  1. 1.

    Introduce a population dynamics model describing a vertically transmitted insect-endosymbiont

  2. 2.

    Describe a previously established insect-vectored plant pathogen model

  3. 3.

    Combine the two, assuming there is a trade-off between vertical and horizontal transmission

  4. 4.

    Use adaptive dynamics methods (Geritz et al. 1998; Dieckmann et al. 2002), in order to find parameter ranges and conditions leading to particular evolutionarily stable states for the combined system.

This allows us to investigate the evolution of transmission method in an ecological context.

By varying the values of parameters in the model, we investigate how ecological changes influence the evolution of insect endosymbionts and aim to address the question:

What ecological changes influence evolution from an insect mutualist to a plant pathogen?

Materials and methods

Endosymbiont model

We model the spread of infection within the insect population, altering the structure and assumptions of a previous insect endosymbiont model (Bernhauerová et al. 2015). We begin by splitting the insect population into uninfected (X) and infected (Z) male and female compartments, with males being denoted by subscript m and females by subscript f. We assume that insects are introduced into the system through sexual reproduction, with mating being modelled with the commonly used harmonic mean function, \({\rm M}\left( {f,m} \right)\) (Caswell and Weeks 1986). The harmonic mean function M(f, m) produces the expected births in the population by taking into account all possible matings that occur in the insect population and the per mating birth rate. The harmonic mean function also allows us to account for different mating systems with parameter -h, such as polygamy (h > 1), polyandry (h < 1) and monogamy (h = 1). The harmonic mean function used in the model is given in Eq. (1), where \(f = X_{f} + Z_{f}\), \(m = X_{m} + Z_{m}\) and c is the insect birth rate.

$${\rm M}\left( {f,m} \right) = 2c\frac{fm}{f/h + m}.$$
(1)

After accounting for all possible matings that occur in the insect population with the harmonic mean function, a proportion of all births will be assigned to the uninfected and infected insect populations, the number of which will be determined by how successful vertical transmission of the pathogen is. For many insect endosymbionts, vertical transmission through female progeny is the sole means of transmission from one insect to another, and so insect reproduction is the means in which an insect endosymbiont can maintain itself in the insect population. However, vertical transmission alone cannot maintain insect endosymbionts in the population (Lipsitch et al. 1995).

Mutualistic relationships between insects and their endosymbionts are commonly observed. Many insect endosymbionts provide a fitness benefit to their insect hosts such as aiding insect digestion of food (Akman Gündüz and Douglas 2012) or increasing the number of eggs laid by the mother (Fast et al. 2011; Pelz-Stelinski and Killiny 2016; Ren et al. 2016). Fitness benefits are however not the only way in which endosymbionts can maintain themselves in the population. For example, Wolbachia spp. maintain themselves in their insect hosts by killing male offspring (Hurst 1997; Engelstädter and Hurst 2009). In our model, we assume that infected female insects lay more eggs than uninfected female insects, modelled with parameter R > 1. Vertical transmission of insect endosymbionts is not always successful (Bressan et al. 2009) so we assume that vertical transmission occurs with efficiency α.

The proportion of matings that result in uninfected insects is:

$$\frac{{\left( {X_{f} + R\left( {1 - \alpha } \right)Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}},$$
(2)

and conversely, the proportion of matings that result in infected insects is:

$$\frac{{R\alpha Z_{f} \left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}}.$$
(3)

The insect population is subject to both density independent and density dependent mortality rate, modelled with parameters \(\varXi\) and b respectively. Assuming that the proportion of male and female insects born is equal, the insect dynamics of the system are thus:

$$\begin{aligned} \frac{{dX_{f} }}{dt} & = \frac{1}{2}{\rm M}\left( {f,m} \right)\frac{{\left( {X_{f} + R\left( {1 - \alpha } \right)Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}} \\ & \quad - \left( {\varXi + b\left( {X_{m} + Z_{m} + X_{f} + Z_{f} } \right)} \right)X_{f} , \\ \end{aligned}$$
(4)
$$\begin{aligned} \frac{{dX_{m} }}{dt} & = \frac{1}{2}{\rm M}\left( {f,m} \right)\frac{{\left( {X_{f} + R\left( {1 - \alpha } \right)Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}} \\ & \quad - \left( {\varXi + b\left( {X_{m} + Z_{m} + X_{f} + Z_{f} } \right)} \right)X_{m} , \\ \end{aligned}$$
(5)
$$\frac{{dZ_{f} }}{dt} = \frac{1}{2}{\rm M}\left( {f,m} \right)\frac{{R\alpha Z_{f} \left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}} - \left( {\varXi + b\left( {X_{m} + Z_{m} + X_{f} + Z_{f} } \right)} \right)Z_{f} ,$$
(6)
$$\frac{{dZ_{m} }}{dt} = \frac{1}{2}{\rm M}\left( {f,m} \right)\frac{{R\alpha Z_{f} \left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}} - \left( {\varXi + b\left( {X_{m} + Z_{m} + X_{f} + Z_{f} } \right)} \right)Z_{m} .$$
(7)

As we intend to extend the insect endosymbiont model to incorporate plant dynamics, we simplify the above insect dynamics by letting \(X = X_{f} + X_{m}\) and \(Z = Z_{f} + Z_{m}\), which implies that dX/dt = dXf/dt + dXm/dt and dZ/dt = dZf/dt + dZm/dt. The insect dynamics are now:

$$\begin{aligned} \frac{dX}{dt} & = 2c\frac{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)/h + (X_{m} + Z_{m} )}}\frac{{\left( {X_{f} + R\left( {1 - \alpha } \right)Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}} \\ & \quad - \left( {\varXi + b\left( {X_{m} + Z_{m} + X_{f} + Z_{f} } \right)} \right)\left( {X_{m} + X_{f} } \right), \\ \end{aligned}$$
(8)
$$\begin{aligned} \frac{dZ}{dt} &= 2c\frac{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)/h + (X_{m} + Z_{m} )}}\frac{{R\alpha Z_{f} \left( {X_{m} + Z_{m} } \right)}}{{\left( {X_{f} + Z_{f} } \right)\left( {X_{m} + Z_{m} } \right)}} \hfill \\ &\quad - \left( {\varXi + b\left( {X_{m} + Z_{m} + X_{f} + Z_{f} } \right)} \right)\left( {Z_{m} + Z_{f} } \right). \hfill \\ \end{aligned}$$
(9)

We assume that a proportion, ρ, of the insect population is male and the remaining proportion, (1 − ρ), is female, hence \(X_{m} = \rho X\), Xf = (1 − ρ)X, \(Z_{m} = \rho Z\) and Zf = (1 − ρ)Z. Equations (8) and (9) become:

$$\frac{dX}{dt} = \frac{{2ch\rho \left( {1 - \rho } \right)}}{1 - \rho + \rho h}\left( {X + R\left( {1 - \alpha } \right)Z} \right) - \left( {\varXi + b\left( {X + Z} \right)} \right)X,$$
(10)
$$\frac{dZ}{dt} = \frac{{2ch\rho \left( {1 - \rho } \right)}}{1 - \rho + \rho h}\left( {R\alpha Z} \right) - \left( {\varXi + b\left( {X + Z} \right)} \right)Z.$$
(11)

In the documented events of insect endosymbiont evolving to become plant pathogens, the sex ratios of the insect populations, P. leporinus for A. phytopathogenicus and D. citri for Candidatus Liberibacter spp. have been equal (Sémétey 2006; Gatineau 2002; Nava et al. 2007). As it is our intention to develop a model to investigate the evolutionary emergence of plant pathogens from insect endosymbionts in such systems, we assume an equal sex ratio, modelled by \(\rho = 1/2\), which simplifies the insect dynamics:

$$\frac{dX}{dt} = c\frac{h}{h + 1}\left( {X + R\left( {1 - \alpha } \right)Z} \right) - \left( {\varXi + b\left( {X + Z} \right)} \right)X,$$
(12)
$$\frac{dZ}{dt} = c\frac{h}{h + 1}\left( {R\alpha Z} \right) - \left( {\varXi + b\left( {X + Z} \right)} \right)Z.$$
(13)

The insect endosymbiont model can produce a stable non-trivial steady state:

$$\left( {\bar{X}, \bar{Z}} \right) = \left( { \frac{{R\left( {1 - \alpha } \right)(Rch\alpha - \varXi \left( {h + 1} \right)}}{{b\left( {R - 1} \right)\left( {h + 1} \right)}},\frac{{\left( {R\alpha - 1} \right)(Rch\alpha - \varXi \left( {h + 1} \right)}}{{b\left( {R - 1} \right)\left( {h + 1} \right)}}} \right).$$

The stability analysis of the insect endosymbiont model with equal sex ratio is given in the supplementary material. We find that the basic reproduction number of the insect endosymbiont system is Rα and that when the steady state exists (\(\bar{Z} > 0\)) it is stable.

Plant host dynamics

We take the plant host dynamics from a previous model (van den Bosch et al. 2006). We model the transition from healthy plants to infected plants using variables H and I to represent the density of healthy and infected plant hosts respectively. Plant density is measured as plants per unit area. In the plant population, planting occurs at a constant rate per unit area, represented by parameter σ. Harvest occurs at rate η, relative to the density of the plant populations. Roguing of infected plants occurs at rate ϕ, relative to the density of the infected plant population. Note that ϕ can also be interpreted as an increase in mortality due to infection in the plant population. Horizontal transmission from uninfected to infected states occurs via the insect vector at rate β1 with mass action dynamics; for this transmission to occur an uninfected insect must land and feed on an infected host.

The system dynamics take the following form:

$$\frac{dH}{dt} = \sigma - \eta H - \beta_{1} ZH$$
(14)
$$\frac{dI}{dt} = - \phi I - \eta I + \beta_{1} ZH$$
(15)

A full list of the parameters used is given in Table 1.

Table 1 A description of the modelling parameters used and the dimensions of parameters
Full size table

Combining the models

We combine the insect endosymbiont model, Eqs. (12) and (13), with the plant dynamics of the insect-vectored plant pathogen model (van den Bosch et al. 2006), Eqs. (14) and (15) and then incorporate horizontal transmission of infection in the insect population, via insects feeding on infected plants.

A trade-off between vertical and horizontal transmission of infection in the insect population is assumed. The assumption of a trade-off between vertical and horizontal transmission is common and backed by empirical evidence (Ewald 1987; Bull et al. 1991; Magalon et al. 2010; Dusi et al. 2015).

Since insect endosymbionts can be transmitted vertically and/or horizontally, we incorporate both transmission methods into the model, Fig. 1. In the insect population, insects are born into the system and a portion of the insect births joins the uninfected insect population, X, with the remaining births joining the infected insect population Z. This portion is determined by the proportion of sexual couplings which may produce a vertically transmitted infection (Eq. 3). As before, both infected and uninfected insects die with an intrinsic mortality rate Ξ and a density dependent mortality rate b. Uninfected insects may transition to the infected state by feeding on infected plants. The horizontal transmission of infection from plants to insect occurs at rate β2 and is assumed to follow with mass action dynamics. For this transmission to occur an infected insect must land and feed on an uninfected host.

Fig. 1
figure 1

Diagram of the insect-vectored plant pathogen dynamics, the boxes represent uninfected and infected states for the insect host (X and Z respectively) and plant host (H and I respectively). In the plant host, transmission of infection occurs horizontally via infected insects at rate β1; harvest occurs at rate η; roguing occurs at rate ϕ and planting rate is represented by σ. In the insect host, transmission occurs both vertically and horizontally; vertical transmission is represented in the birth-rate term c and horizontal transmission occurs via feeding on infected plant hosts at rate β2; insect death occurs with a density dependent mortality rate, b, and an intrinsic mortality rate Ξ. Full descriptions of these parameters including their dimensions are given in Table 1

Full size image

Insect endosymbionts can be split into two categories, obligate and facultative. Obligate endosymbionts are solely vertically transmitted, facultative endosymbionts are predominantly vertically transmitted but also can be horizontally transmitted. Insect vectored-plant pathogens are predominantly or only transmitted horizontally. The trade-off between horizontal and vertical transmission, which we take as a proxy for an endosymbiont/pathogen lifestyle, is based on the location where the insect endosymbiont resides within the insect, which then corresponds to how the endosymbiont is transmitted, i.e. vertically when the insect endosymbiont resides near and then within the insect mother’s eggs or horizontally when the insect endosymbiont resides in the insect’s mouthparts (Bressan 2014).

In the model, the trade-off is represented with three functions α(θ), β1(θ) and β2(θ). Parameter θ represents the proportion of the endosymbionts which reside outside the insect ovaries. This definition is chosen to try and differentiate between the two transmission routes, vertical transmission through the ovaries and horizontal transmission through the insect’s mouth. This definition provides a clear way to model a single trade-off between a vertically transmitted lifestyle and various categories of horizontally transmitted plant pathogen lifestyles (i.e. circulative propagative and circulative non-propagative). If endosymbionts reside outside the ovaries of the insect (e.g. in salivary glands), we assume they can transmit through the mouth parts of the insect. As θ is a proportion, it is bounded between 0 and 1. Given our choice of trade-off curves, higher levels of θ indicate higher levels of horizontal transmission; lower levels of θ represent a higher level of vertical transmission.

The system dynamics are now as follows:

$$\frac{dH}{dt} = \sigma - \eta H - \beta_{1} \left( \theta \right)ZH,$$
(16)
$$\frac{dI}{dt} = - \phi I - \eta I + \beta_{1} \left( \theta \right)ZH,$$
(17)
$$\frac{dX}{dt} = c\frac{h}{h + 1}\left( {X + R\left( {1 - \alpha \left( \theta \right)} \right)Z} \right) - \left( {\varXi + b\left( {X + Z} \right)} \right)X - \beta_{2} \left( \theta \right)XI,$$
(18)
$$\frac{dZ}{dt} = c\frac{h}{h + 1}\left( {R\alpha \left( \theta \right)Z} \right) - \left( {\varXi + b\left( {X + Z} \right)} \right)Z + \beta_{2} \left( \theta \right)XI.$$
(19)

We aim to model the evolutionary transition from a primarily vertically transmitted insect endosymbiont to a primarily horizontally transmitted plant pathogen. It would, therefore, be justified to choose trade-off curves which would favour the original vertically transmitted insect endosymbiont lifestyle.

One way of doing this is choosing a curve such that the majority of values of θ vertical transmission is greater than horizontal transmission, meaning that more endosymbionts are transmitted through the ovaries than through the mouth parts of the insect.

Initially, a concave trade-off curve for vertical transmission is chosen: α(θ) = αmaxA(θ − 1)/(θ − A). Parameter A > 1 changes the degree of concavity of the curve, with values closer to 1 increasing the degree of concavity. As vertical transmission is linked to the proportion of insect endosymbiont residing in the ovaries, θ, the vertical transmission function would necessarily be a decreasing function of θ, a property achieved by the concave trade-off curve. Horizontal transmission via the mouthparts would increase as the proportion of insect endosymbiont residing in the insect ovaries decreased (modelled by increasing θ). We therefore model the horizontal transmission trade-off curves as increasing functions of θ. We assume the horizontal transmission trade-off curves are linear, \(\beta_{1} \left( \theta \right) = \beta_{{1_{max} }} \theta , \beta_{2} \left( \theta \right) = \beta_{{2_{max} }} \theta\), see Fig. 2. Given our chosen trade-off curve, when θ = 0, there is only vertical transmission \(\left( {\beta_{1} \left( 0 \right) = \beta_{2} \left( 0 \right) = 0, \alpha \left( 0 \right) = \alpha_{max} } \right)\), representing a solely vertically transmitted insect endosymbiont, incapable of causing infection in the plant. When θ = 1, then there is only horizontal transmission \(\left( {\beta_{1} \left( 1 \right) = \beta_{{1_{max} }} , \beta_{2} \left( 1 \right) = \beta_{{2_{max} }} , \alpha \left( 1 \right) = 0} \right)\), and the symbiont is a plant pathogen, spread purely from insect to plant and from plant to insect, without insect to insect transmission. Intermediate values of θ represent mixed mode transmission, where plant hosts can be infected by the symbiont and the symbiont can be transmitted vertically.

Fig. 2
figure 2

Examples of the vertical and horizontal transmission trade-off function used in the model, the concave function describes vertical transmission efficiency, \(\alpha \left( \theta \right) = \alpha_{max} \left( {\theta - 1} \right)/\left( {\theta - A} \right)\) and the linear functions \(\beta_{1} \left( \theta \right), \beta_{2} \left( \theta \right)\) describe horizontal transmissions between insect and plant hosts. As shown an increase in the concavity of the vertical transmission from decreasing the value of A result in an increase in the value of θ for the intercept between the vertical transmission term, α(θ) and the horizontal transmission terms \(\beta_{1} \left( \theta \right), \beta_{2} \left( \theta \right)\)

Full size image

Currently, not enough data have been gathered on where endosymbionts reside in an insect host to be used to derive a trade-off curve. Therefore, a theoretical trade-off curve is currently our best option in modelling this system. Other trade-off curves are explored in the supplementary material.

In our model, the trade-off curve shaping parameters influence the amount of transmission possible between hosts. The vertical transmission term is shaped by parameters αmax and A, the horizontal transmission terms are shaped by the parameters \(\beta_{{1_{max} }}\) and \(\beta_{{2_{max} }}\). It is likely that these parameters will influence the evolutionary change in transmission route, but more the more interesting question is which of the other modelling parameters result in an evolutionary change in transmission route.

Model dynamics

Stability analysis of the combined model

In the combined model, there are three sets of steady states, \(\left( {\bar{X},\bar{Z},\bar{H},\bar{I}} \right)\). These are a disease-free state:

$$\left( {\frac{{ch - \left( {h + 1} \right)\varXi }}{{b\left( {h + 1} \right)}},0,\frac{\sigma }{\eta },0} \right),$$

an insect-free and therefore plant disease-free state:

$$\left( {0,0,\frac{\sigma }{\eta },0} \right),$$

and a non-trivial steady state. The stability analysis of the disease free-state and insect-free state are given in the supplementary material. The existence and stability of the non-trivial equilibrium point was determined through numerically solving the system governed by Eqs. 16–19 using the programming language Python and the ‘odeint’ function in the package scipy until the system was at rest. Figure 3 displays an example of the combined model solutions converging towards the stable non-trivial equilibrium point.

Fig. 3
figure 3

The full system can produce a stable non-trivial equilibrium point as shown in this figure. Here the system is simulated using the odeint function from the Python library, scipy, until the equilibrium points are at rest at the non-trivial steady states \(\left( {\bar{X},\bar{Z},\bar{H},\bar{I}} \right)\), denoted with the dotted line. In the making of this plot the following parameters were used: \(c = 3, h = 1, R = 1.5, \varXi = 0.1, b = 0.1, \sigma = 1, \eta = 0.1\), Φ = 0.05 and the trade-off curve shaping parameters \(A = 1.2,\alpha_{max} = 0.95, \theta = 0.5, \beta_{{1_{max} }} = \beta_{{2_{max} }} = 0.3\)

Full size image

Adaptive dynamics

Adaptive dynamics allows us to model how a particular trait will evolve in response to its environment (Geritz et al. 1998). We choose our trait of interest and find its associated parameter(s) within the model. The environment is represented by the other model parameters. Here, we investigate the evolution of the proportion of the endosymbionts which reside outside the insect ovaries, θ. This proportion corresponds to our trait of interest, the transmission method of a pathogen/endosymbiont. The remaining model parameters in the system are used to determine the environment. The environmental insect dynamic parameters are the birth rate, c, per capita death rate, Ξ, density dependent death rate, b, fitness benefit, R and the insect mating system parameter, h. In the plant dynamics, the environmental parameters are the planting rate, σ, the harvesting rate, η, and the roguing rate, Φ.

Using adaptive dynamics, we can see how a parameter within a model can evolve to an optimum point, referred to as a convergent stable strategy (CSS). A CSS is a strategy that is both convergence stable, meaning that the strategy is selected through evolution, and is an evolutionary stable strategy (ESS). An ESS in this model’s context is a parameter value of θ which cannot be invaded, outcompeted and subsequently replaced by any other rare invading strategy (a rare invader with an alternate value of θ). One of the central assumptions of the adaptive dynamics’ framework is that the invading mutant phenotype, represented here with parameter θi, is always assumed to be rare (Brännström and Festenberg 2007; Geritz et al. 1998). Another relevant category of CSS is the locally stable CSS, which is a CSS point of θ where rare neighbouring values of θ cannot invade the local CSS value of θ, but there exists values of θ in the parameter space which could outcompete the local CSS value of θ (McGill and Brown 2007). Often the CSS is dependent on the other environmental parameters within the model. Changes in the value of the CSS represent evolutionary changes in the trait of interest. Finding the parameters which correspond to evolutionary changes in the trait of interest can help us identify environmental sources of selection.

To use the adaptive dynamics methods, we need to further expand this system to introduce the resident phenotype and invading phenotype. To do this we assume that the system is at equilibrium with a resident infection with trait θr. This will determine the biotic environment experienced by any initially rare invading mutant strain. The invader is assumed to be an evolutionary descendent of the resident infection with a slightly different trait value θi. The infection takes place in both infected insect and plant hosts; thus, our system becomes expanded to include equations for Zr and Ir, the insect and plant population host to the resident infection with trait θr; and Zi and Ii, the insect and plant population host to the invading infection with trait θi.

Initially, we consider only the resident system (Zi = Ii = 0) as we want to determine the environment the invading mutant encounters. This model is equivalent to the ‘combined’ model above. As such, we have already shown the steady states of this system and determined them to be biologically realistic, so with this we can explore the invasion dynamics. The invader system initially encounters the resident system when the resident system is at equilibrium, the invading phenotype in the presence of the resident at its equilibrium population is determined by:

$$\frac{{dZ_{i} }}{dt} = c\frac{h}{h + 1}\left( {R\alpha \left( {\theta_{i} } \right)Z_{i} } \right) - \left( {\varXi + b\left( {\bar{X} + \overline{{Z_{r} }} + Z_{i} } \right)} \right)Z_{i} + \beta_{2} \left( {\theta_{i} } \right)\bar{X}I_{i} ,$$
(20)
$$\frac{{dI_{i} }}{dt} = - \phi I_{i} - \eta I_{i} + \beta_{1} \left( {\theta_{i} } \right)Z_{i} \bar{H},$$
(21)

where \(\bar{X}, \bar{H}\) and \(\overline{{Z_{r} }}\) have been numerically estimated by taking the endpoints of the simulated resident system using the odeint function in the Python’s scipy library.

When the equilibrium point \(\left( {\bar{X},\overline{{Z_{r} }} ,\bar{H},\overline{{I_{r} }} ,0,0} \right)\) is unstable, then the invading mutant can successfully invade and take over the system. A simple way of finding this criterion is to evaluate the Jacobian matrix of the invader system at (0, 0):

$$J\left( {0,0} \right) = \left( {\begin{array}{*{20}l} {\frac{{R\alpha_{i} ch}}{h + 1} - \left( {\varXi + b\left( {\bar{X} + \overline{{Z_{r} }} } \right)} \right)} \hfill & {\beta_{2} \left( {\theta_{i} } \right)\bar{X}} \hfill \\ {\beta_{1} \left( {\theta_{i} } \right)\bar{H}} \hfill & { - \left( {\phi + \eta } \right)} \hfill \\ \end{array} } \right)$$

if there is at least one eigenvalue of J(0, 0) with a positive real part, then the equilibrium point is unstable (Britton 2003). In terms of the biology the system is representing, this means that the invader population must increase from zero.

To find an evolutionary endpoint in this system, we numerically evaluated the invasion criterion with a fairly dense sampling of possible pairs of (θr, θi) where 0 ≤ θr, θi ≤ 1, then plot on a graph where an invader with trait value θi is successful in invading a resident with trait value θr with a black dot. Unsuccessful invasion pairs are left blank; this type of plot is called a ‘Pairwise Invasibility Plot’ (PIP) (Kisdi et al. 1993). An annotated example of a PIP is shown in Fig. 4.

Fig. 4
figure 4

An example of an evolutionary singular strategy with convergence stability, otherwise known as a continuously stable strategy (CSS) shown on a pair-wise invasibility plot (PIP). The invasion process begins when a mutation occurs in trait θr,the mutated trait value is denoted as θi. We assume that prior to this mutation our population is homogenous in trait θr. This is shown on the PIP by following the line θr = θi, highlighted in red; a mutation in trait θr will result in a value of θi which is either bigger or smaller than θr, meaning that the trait value \(\theta_{i}\) will be directly above or below a point on the red line. This trait value, θi, determines whether the invasion of the resident system, defined by parameter θr, is successful or not. Successful invasions are denoted by a black dot, unsuccessful invasions are left white. When an invasion is successful, the value of θi which successfully invaded the resident trait value, θr, becomes the new resident trait value. Until a CSS is reached, evolution will go through a stage of repeatedly substituting various values of θ, the red arrows on the PIP indicate this process happening; the vertical arrows indicate where successful invasion takes place, the horizontal arrows indicate the new value of θi establishing itself as the new resident, eventually evolving towards the CSS. The CSS is shown on the PIP as a point along the line θr = θi with a white area directly above and below, this means that we have found a trait value of θr where neighbouring values of θi cannot invade the resident value

Full size image

We can then construct PIPs across the parameter space, to see which parameters change the evolutionary outcome. In the present case, we can determine which factors are likely to influence the evolution of the transmission method of an insect endosymbiont/plant pathogen.

Finding a suitable range for each parameter is difficult. Some have natural bounds, for example, vertical transmission efficiency is bounded between 0 and 1, but others are more ambiguous. We now apply units to our parameters. The hosts are measured as the density of hosts per unit area over time, we take our unit area to be hectares and time as hours. For the parameters we could not find an estimate of, we will investigate the evolutionary dynamics over a wider range of plausible parameter interactions, doing so until the consequences of changing the parameter’s value are clear. This method only takes a subset of the infinite number of possible parameter combinations, but with this method we are still able to see how the model’s parameters influence the evolution of transmission method. The range of values tested in our evolutionary simulations is given in Table 1.

We found an estimated range for σ by looking at the recommended density for sugar beet hosts per hectare, the British Beet Research Organisation recommends a plant density of 100,000 plants per hectare (BBRO 2018). We use the average introduced plant density per hour as the estimate for σ. Sugar beet is an annual crop so this is established with by the following: \(100,000/\left( {24 \times 365} \right) \approx 12\), this is taken as the upper bound for \(\sigma\).

We tested a large subset of the parameter interactions within the chosen parameter range, calculating the evolutionary endpoints of multiple parameter combinations spanning the ranges stated in Table 1 (see supplementary material). We state the general effect of increasing or decreasing the model’s parameters on the value of the global/local CSS of θ in the results section.

Results

Through adaptive dynamics, we are able to find which parameters could cause an evolutionary shift in transmission method. An evolutionary shift in lifestyle from mutualistic endosymbiont to an insect-vectored plant pathogen due to changing environments (determined by the modelling parameters) would be represented by increasing the value of the local or global CSS for parameter θ. Conversely, an evolutionary shift in lifestyle from insect-vectored plant pathogen to insect endosymbiont would be represented by decreasing the value of the local or global CSS for parameter θ. For an illustration of the interaction between the environmental modelling parameters, the pairwise invasibility plots and the evolution of transmission method see Fig. 5.

Fig. 5
figure 5

We see here that changes in the environmental parameter set in step 1, (in this example the planting rate, σ) in part determines the environment set by the resident system in step 2 (for example simulated here with θ = 0.5). In step 3, we simulate the evolutionary dynamics and show that the changes to the resident environment from differing σ produce significantly different evolutionary dynamics. This can be seen by comparing the differing the evolutionary endpoint of the transmission method θ (marked by a red dot) in rows a and b. How the evolutionary endpoint of θ relates to the endosymbiont’s transmission method is shown in step 4, the red dashed vertical line corresponds to the evolutionary endpoint of transmission method, θ*, calculated in step 3. The changes to the environmental parameter σ in rows a and b resulted in the CSS θ* = 0 when σ = 0.1 (row a) and the CSS at θ* ≈ 0.68 when σ = 10 (row b). The CSS θ* = 0 corresponds to vertical transmission only as \(\left( {\beta_{1} \left( 0 \right) = \beta_{2} \left( 0 \right) = 0,\alpha \left( 0 \right) = \alpha_{max} } \right)\), whereas the CSS at θ* ≈ 0.68 corresponds to mixed mode transmission as \(\beta_{1} \left( {0.68} \right),\beta_{2} \left( {0.68} \right),\alpha \left( {0.68} \right)\) are all non-zero. As the only parameter to change between rows a and b was the planting rate σ, thus we can conclude that increased planting rate results in the evolution of mixed mode transmission and therefore a plant-pathogen lifestyle

Full size image

In the plant dynamics, Eqs. (14) and (15), the environmental parameters tested were the planting rate (σ), the harvest rate (η) and the roguing rate (ϕ), Fig. 6. With lower values of σ, there existed a local CSS at θ = 0, meaning that the endosymbiont is transmitted entirely through vertical transmission, however, increases in planting rate (σ, row a) changed the evolutionarily dynamics such that a local CSS appeared for an intermediate value of θ, meaning that a combination of horizontal and vertical transmission could be maintained in the population. Increasing the harvest rate (η, row b) and the roguing rate (ϕ, row c) changed the evolutionary dynamics such that a sole locally stable CSS existed corresponding with only vertical transmission, with roguing appearing to be less sensitive. Simply put these results imply that if we introduce more plant hosts in the system by increasing the planting rate (σ), horizontal transmission becomes an evolutionary viable strategy as a local CSS appears. Conversely, if we reduce the density of plant hosts in the system by increasing harvesting (η) and/or roguing (Φ), then the effects of increasing planting rate disappear and the only viable transmission strategy is vertical transmission.

Fig. 6
figure 6

Sensitive parameters in the plant dynamics, showing shifts in the predominant transmission method, and thus lifestyle. Mathematically these shifts are represented by changes in the local evolutionary stable strategy, shown here with red dots. Increases in planting rate of uninfected hosts (σ, row a) will result in favoured horizontal transmission and a plant pathogen lifestyle. Increasing harvesting rate (η, row b) and roguing rate (ϕ, row c) will reverse these changes, resulting in favoured vertical transmission. In all rows the following parameters were used unless otherwise stated in the subplot’s title: \(c = 3, h = 1, R = 1.5, \varXi = 0.1, \sigma = 10, \eta = 0.1, \varPhi = 0.05, A = 1.2, \alpha_{max} = 1,\) and \(\beta_{{1_{max} }} = \beta_{{2_{max} }} = 0.3\)

Full size image

Within the insect dynamics, Eqs. (18) and (19), the parameters representing the environment were: the fecundity benefit ratio in the infected insect population (R); the insect growth rate (c); the insect mating system (h); the per capita insect mortality rate (Ξ); and the density dependent insect mortality rate (b). The role these parameters played in evolutionary change of the CSS value for transmission method was investigated by changing the values for these parameters, and then observing the resulting change, if any, of the CSS value. Sample results are shown in Fig. 7, in which the pairwise invasibility plots demonstrate changing evolutionary endpoints in transmission method according to altered parameter values. Increases to the insect growth rate, c, density dependent insect mortality rate, b, and the fitness benefit provided by the insect endosymbiont, R, resulted in favoured vertical transmission. The insect mating system also proved significant with monogamy (h = 1) and polyandry (h < 1) resulting in mixed mode transmission. The intrinsic insect mortality rate, Ξ, did not qualitatively change the optimum transmission mode.

Fig. 7
figure 7

Parameters in the insect dynamics, showing shifts in the predominant transmission method, and thus lifestyle. Mathematically these shifts are represented by changes in the locally stable evolutionary stable strategies. Increases in insect birth rate (c, row a) will result in favoured vertical transmission a shift to an insect endosymbiont lifestyle. The mating system (h, row b) proved to be evolutionarily significant, with polygamy (h > 1) resulting in evolution towards a favoured insect endosymbiont lifestyle and polyandry (h < 1) resulting in the possibility of mixed mode transmission. Changes resulting from increases in insect mortality differ depending on whether they are intrinsic (Ξ,  row c) or density dependent (b row d); the intrinsic mortality rate did not significantly alter the evolutionary dynamics, but increases density dependent mortality led to favoured vertical transmission and an insect endosymbiont lifestyle. Finally, increasing the fitness benefit (R, row d) of infection led to favoured vertical transmission and an insect endosymbiont lifestyle. Unless otherwise stated, the parameters used in the making of these PIPs in rows a and b were: \(c = 3, h = 1, R = 1.5, \varXi = 0.1, b = 0.1, \sigma = 1, \eta = 0.1, \varPhi = 0.05, A = 1.2, \alpha_{max} = 0.95\) and \(\beta_{{1_{max} }} = \beta_{{2_{max} }} = 0.2\). The same parameters were used in rows c to d, with the exception of c = 10 and σ = 3.5

Full size image

We also investigated the parameters used to shape the trade-off functions, α(θ), β1(θ) and β2(θ): \(A, \alpha_{max}\), \(\beta_{{1_{max} }}\) and \(\beta_{{2_{max} }}\) Fig. 8. Increasing the shaping parameter A, which reduces the concavity of the trade-off curve and the degree in vertical transmission is larger than the horizontal transmission terms, results in an increase in the CSS value of θ (row a). Increases in αmax result in favoured vertical transmission (row b) and increases in \(\beta_{{1_{max} }}\) or \(\beta_{{2_{max} }}\) (row c and d respectively) result in favoured horizontal transmission. These results are unsurprising, as these parameters directly determine how effective each transmission route can be, so increasing the effectiveness of one transmission route should increase selection towards that transmission route.

Fig. 8
figure 8

The trade-off curve shaping parameters alter the evolutionary dynamics. Increasing A, which decreases the degree of concavity of the trade-off curve results in an increase in the CSS value of θ, corresponding to an increase in horizontal transmission. Increasing αmax reduces the local CSS and increases vertical transmission. Increasing \(\beta_{{1_{max} }}\), \(\beta_{{2_{max} }}\) results in an increase in the prevalence of horizontal transmission in mixed mode transmission, however, the increase of the CSS value is more notable for increasing \(\beta_{{1_{max} }}\)

Full size image

The figures shown here are a select few taken to illustrate the general trend seen in our analysis.

The shape of the trade-off function also affected the evolutionarily stable transmission method. Different trade-off curves shared the same qualitative results as described above, but some notable differences in the shape of the PIPs were evident. Evolutionary bi-stability was far more common when other trade-off curves were used, compared to the concave trade-off curve (see supplementary material).

Central to our model derivation is the assumption of an equal insect sex ratio, based on the roughly equal sex ratio reported in the vector of A. phytopathogenicus (Sémétey 2006; Gatineau 2002; Nava et al. 2007). The robustness of our results given this assumption are considered further in the supplementary material, but we find that altering the insect sex ratio in either direction (male majority or female majority) had a reductive effect on the evolutionary dynamics. Specifically, whether the CSS increased or decreased due to environmental parameter value changes did not change, but the magnitude in which the CSS changed due to an altered environmental was reduced, meaning the results produced here are robust when considering different sex structures in insect populations (see supplementary material).

Discussion

With this model, we begin to investigate why an insect endosymbiont may evolve to become a plant pathogen. We developed a model, based on previously established models of the component parts (Jeger et al. 2004; Madden et al.2007; van den Bosch et al. 2006; Bernhauerová et al. 2015), which allows for both horizontal and vertical transmission of a symbiont, the trade-off being based on the position of the endosymbiont within the insect host. Using adaptive dynamics, we can look for parameters in the model which, when changed, alter the evolutionary endpoint of the endosymbionts’ transmission method. The transmission method of the endosymbiont was used as a proxy for its lifestyle, with horizontal transmission representing a plant pathogen lifestyle and vertical transmission representing a mutualistic insect endosymbiont lifestyle.

Previous attempts to understand evolutionary changes in symbiosis have used a wide variety of modelling techniques, each with their own strengths and weaknesses (Hoeksema and Bruna 2000). The strength of using an adaptive dynamic modelling approach is in the ecological background it provides, which makes it easier to validate this model and the ecological changes it reveals with empirical methods than game theory models, where evolutionary changes are explained in terms of trades of resources.

Planting rate—σ

Increasing the parameter σ leads to an increased healthy plant density per unit area. Increased host density is known to have epidemiological consequences for both vectored and non-vectored plant disease. For example, increasing host density alters the logistical ease of spread of the disease, by decreasing the distance the vector has to travel between hosts. Increasing host density can also affect the plant host, changing the growth rate, size shape and nutrient status (Burdon and Chilvers 1982; Rodelo-Urrego et al. 2013). Previous models researching the horizontal-vertical transmission trade-off have also indicated that an increase in the density of susceptible hosts will affect the optimum balance between horizontal and vertical transmission (Turner et al. 1998).

Our model suggests that an increase in the planting rate may lead to the transition of an endosymbiont towards a plant-pathogen. This can be interpreted as more intensive agricultural practices. Here we present two real-life examples of how this may have happened in the past.

A parallel can be drawn between the production of sugar beet and strawberries. The modern varieties of each crop were bred in the 18th century, which is recent in terms of evolutionary time. In the case of the evolution of A. phytopathogenicus, a pathogen of sugar beet, and Candidatus Phlomobacter fragariae, a pathogen of strawberries, our model suggests that the increase in production would select for an insect endosymbiont to evolve to a plant pathogen. An increase in the production of these crops can be represented in the model by an increase in σ, leading to an evolution towards higher amounts of horizontal transmission, Fig. 6.

Harvesting rate \(\left( \eta \right)\) and roguing rate (ϕ)

An increase in the harvesting and roguing rates results in selection towards an insect endosymbiont lifestyle, which can effectively reverse the effects of increased planting rate, Fig. 6. This result isn’t surprising as an increase in η decreases the overall amount of plant hosts in the resident system. Similarly, an increase in the roguing rate selects for an insect endosymbiont lifestyle. This is significant as it provides evidence to suggest that roguing can be an effective counter measure in preventing the evolutionary emergence of insect vectored plant pathogens, evolving from insect endosymbionts. However, the logistics of large scale roguing can be problematic.

Insect birth rate—c

Cixiid planthoppers have been known to shift host, followed by explosive population growth and, by extension, an increase in insect population density. Mathematically, this can be represented by increasing the insect’s birth rate. An example is the planthopper Hyalesthes obsoletus. This planthopper changed primary hosts from bindweed (Darimont and Maixner 2001) to nettle (Johannesen et al. 2008). Since this host change, the population size has dramatically increased. Similarly, P. leporinus, the host and vector of A. phytopathogenicus was previously observed using the common reed, Phragmites australis L. as its primary host (Bressan et al. 2009). More recently P. leporinus has been found feeding and reproducing within sugar beet-wheat systems, indicating a shift in primary host (Bressan et al. 2010). Our model predicts that an increase in insect birth rate results in favoured vertical transmission Fig. 7. This result is interesting as if this host shift in P. leporinus increased insect population density (represented in the model as an increase in insect birth rate), as it did in the plant hopper H. obsoletus, then the model would suggest that the increase in insect growth rate would not explain why A. phytopathogenicus evolved from an insect endosymbiont to a plant pathogen and could suggest that even despite the potentially increased birth rate experienced in the insect population resulting from a plant host shift, an insect vectored plant pathogen lifestyle would still be selected for.

Density dependent insect mortality rate—b

An increase in the density dependent insect mortality rate selects for an insect endosymbiont lifestyle, Fig. 7. Interestingly the intrinsic insect mortality rate, Ξ, did not alter the evolutionary dynamics. This result suggests that increased insect mortality due to competition or other density-dependent mortality factors may prevent the evolutionary emergence of insect vectored plant pathogens, evolving from insect endosymbionts and a reduction in density dependent mortality may result contribute to the evolutionary emergence of new plant pathogens from insect endosymbionts. A lower density dependent mortality rate would be experienced by a population establishing a new environment, meaning that if an insect’s plant host changed and therefore the insect’s environment changed to that of a managed plant population, the resulting reduction in density dependent mortality could alter the evolutionary optimum of transmission method so that horizontal transmission and therefore a plant pathogen lifestyle may be viable.

Fecundity benefit ratio—R

A reduction in the fecundity benefit ratio selects for higher levels of horizontal transmission and a plant pathogen lifestyle, Fig. 7. Several external factors could affect the fecundity benefit ratio. One of the most obvious is heat stress. An example comes from a study of endosymbionts in Aphids. Heat stress reduces the number of bacteriocytes in which obligate symbionts reside in aphids, thereby reducing endosymbiont number (Montllor et al. 2002). If the fecundity benefit of an endosymbiont was dependent on the endosymbiont population size within an insect host, which seems very likely then a reduction in endosymbiont population size would reduce fecundity.

Insect mating type—h

Whether the insect mating system was polygamous, monogamous or polyandrous altered the evolutionary optimum value of transmission method, with monogamous and polyandrous mating systems causing the possibility of local CSSs where horizontal transmission is viable, whereas polygamous mating systems favoured sole vertical transmission and polyandry mixed mode transmission Fig. 7. This may be due to the relationship between the fact that vertical transmission occurs only through the female and the increased availability of females in polygamous mating systems. If more female insects are available for mating per male as in a polygamous system, then the benefit of vertical transmission via infected females to the insect endosymbiont population would be increased as this mating system offers increased opportunities of vertical transmission to the endosymbiont. Interestingly, the vector of citrus greening disease, D. citri, has a polyandrous mating system (Wenninger and Hall 2008). After reviewing the literature, we could not find any published records of the mating systems in P. leporinus or C. wagneri, vectors of the pathogens A. pyhopathogenicus and Candidatus Phlomobacter fragariae.

The role of mating system in the evolutionary dynamics presented here differs from another modelling study, focusing instead on the evolution of male-killing, which showed greater levels of polygamy resulting in favoured horizontal transmission (Bernhauerová et al. 2015). The opposite selection in transmission route between Bernhauerová et al. and our model is likely due to the ecological differences in the systems being modelled but both models highlight the importance of insect mating system on insect endosymbiont evolution.

Trade-off curve shaping parameters, \(\alpha_{\hbox{max} } , \beta_{1\hbox{max} } , \beta_{2\hbox{max} }\) and A

The parameters αmax, β1max and β2max reduce the maximum values the trade-off curves α(θ), β1(θ) and β2(θ), whereas A alters the degree of concavity of α(θ). When either αmax, β1max or β2max is increased, the evolutionary endpoint shifts such that vertical transmission is increased (increasing αmax) or horizontal transmission is increased (increasing \(\beta_{1max} , \beta_{2max}\)). The change in CSS is more notable for β1max than in β2max, the maximum value for transmission from infected insect to uninfected plants, Fig. 8. An increase in the concavity of the function α(θ) results in a larger range of values of θ in which α(θ) is greater than β1(θ) and β2(θ) (e.g. see how the increase in concavity of α(θ) increases the θ intercept value between the vertical and horizontal transmission the trade-off curves Fig. 2). The evolutionary changes that result from different values of A reflect this, with horizontal transmission becoming more prevalent when A is increased, effectively reducing the benefit of vertical transmission over horizontal transmission.

There are a number of ecological factors which could be linked to these parameters. One example is climate change. Changes in temperature and drought are linked to climate change (Arnell 2008). Vertical transmission of endosymbionts is affected by temperature in a range of insects: in one review, insects exposed to both hot and cold temperature treatments showed reduced vertical transmission of endosymbionts (Corbin et al. 2016; Chen et al. 2009). Another example comes from heat stress in insects. Heat stress reduces the population size of insect endosymbionts (Wernegreen 2012). With fewer endosymbionts to transmit, the amount of vertically transmitted endosymbionts will likely be reduced. Both processes can be modelled by reducing αmax, the evolutionary consequences of reducing αmax would result in horizontal transmission becoming a viable strategy and could explain the emergence of novel insect vectored plant pathogens. If the model focused on the scale of the endosymbiont population, this result would not be applicable as heat stress would induce extra mortality, rather than reducing the efficiency of vertical transmission. Since the model was built on the scale of host density, heat stress would effectively reduce vertical transmission efficiency. The role of temperature on the evolutionary emergence of plant pathogens could be investigated formally by expressing the model’s parameters as functions of temperature, as has been done with other modelling studies (Berec 2019). The intention of this study was to generate hypotheses as to why this shift has occurred, rather than to investigate the role of any one mechanism, such as temperature, in causing these shifts. As such, the formal investigation of temperature in the shift has not been done here, but would make for interesting future research.

Changes to the shaping parameter \(\beta_{{1_{max} }}\) could be explained by the symptoms of infection in a plant host. Insect-vectored plant viruses can ‘entice’ insects to land on infected plant hosts (Mauck et al. 2012; Mwando et al. 2018). This mechanism would indirectly increase the transmission of infection from plant to insect. Citrus greening disease in citrus, syndrome basses richesses in sugar beet and marginal chlorosis in strawberries all alter the appearance of their host plants. No study has investigated vector feeding preference for infected or uninfected hosts in syndrome basses richesses in sugar beet or marginal chlorosis in strawberry systems. However, there is evidence that the altered appearance of hosts infected by citrus greening disease is more attractive to the psyllid D. citri (Mauck et al. 2012). Another study on the feeding response of the citrus greening disease vector D. citri to plant hosts infected by Candidatus liberibacter asiaticus and non-infected hosts focused on the volatiles released by the infected plant host (Mann et al. 2012), a significant difference was found in host preference, with psyllids preferring to feed on infected hosts. The change in volatiles produced by infected hosts to increase attractiveness to insects is seen in other insect-vectored plant pathogen systems (Shapiro and Turner 2014; Mas et al. 2014).

Feeding studies of the Asiatic citrus psyllid, another of the vectors of citrus greening disease, indicate that infected citrus hosts are a less favourable food source to this insect. However, the duration of sustained ingestion from phloem, a process required for the acquisition and transmission of the pathogen in the insect population (Bonani et al. 2010), increases with the symptom severity in infected citrus hosts (Cen et al. 2012). This implies that an increase in the transmission between both insect to plant host, and plant to insect is present in this system. This means that the direction of evolution is unclear, we do not know if this increase in \(\beta_{{1_{max} }}\) sparked an evolutionary transition in Candidatus Liberibacter spp., from an insect endosymbiont to a plant pathogen, or Candidatus Liberibacter spp. emerged as a plant pathogen and then an increase in transmission efficiencies was selected for.

As reflected in our model Fig. 8, these mechanisms would cause evolution to favour horizontal transmission; if an insect endosymbiont could evolve the ability to attract insects to the infected plant hosts by some means, our model shows that this could make it more likely that the endosymbiont would evolve to be a plant pathogen.

Insect sex ratio—ρ

In the supplementary material, we explore what effect an altered insect sex ratio had on the evolutionary dynamics, finding that a male dominated or female dominated population had a dampening effect on the evolutionary dynamics. Specifically, the increase or decrease in the CSS values resulting from changes in the environmental model parameters remained the same, but the magnitude of change in the CSS value was reduced (see supplementary material). This means the evolutionary dynamics remain the same when we relax the assumption of an equal insect sex ratio.

Why would insect endosymbionts evolve away from established vertical transmission to horizontal transmission?

We believe that it is competition amongst endosymbionts that explains this. Say the invading mutant phenotype evolves some ability to transmit horizontally and encounters a resident system such that horizontal transmission via the plant host provides a fitness benefit (for example one with a high density of plant hosts). Further, assume that the resident phenotype is entirely vertically transmitted. Despite vertical transmission being a certain way of ensuring the endosymbiont’s survival, having the ability to transmit infection horizontally (even at the expense of the more certain vertical transmission route), could give the invading pathogen a sufficient fitness benefit such that the invading pathogen would outcompete the resident phenotype. The invader phenotype would then establish itself in the population and becoming the new resident phenotype. This therefore implies that subsequent generations of the pathogen would need some level of horizontal transmission to remain competitive in this new environment, thus explaining why the vertical transmission route could become less prevalent.

Conclusion

From our modelling approach, we propose that the following factors can cause an evolutionary transition from mutualistic endosymbiont to insect-vectored plant pathogen lifestyles:

  1. 1.

    Increases in temperature, possibly due to climate change, which may reduce the vertical transmission efficiency of infection in the insect population or the fecundity benefit given to infected insects.

  2. 2.

    A primary host shift in the insect, from an uncultivated crop to a more abundant cultivated one, resulting in decreased density dependent mortality in the insect population.

  3. 3.

    An increase in population density of the host plant per hectare.

  4. 4.

    An increase in the capability to transmit infection horizontally.

We also find a method of control against the evolutionary emergence of an insect-vectored plant pathogen, from an insect endosymbiont lifestyle, namely roguing of infected host plants.

The information available concerning former insect endosymbionts which evolved to become insect-vectored plant diseases—including citrus greening disease in citrus (Gottwald 2010; Gottwald and Mccollum 2017), syndrome basses richesses in sugar beet and marginal chlorosis in strawberries (Bressan et al. 2012)—is consistent with these driving changes.