Spatially structured eco-evolutionary dynamics in a host-pathogen interaction render isolated populations vulnerable to disease

The pathosystem

Plantago lanceolata L. is a perennial monoecious ribwort plantain that reproduces both clonally via the production side rosettes, and sexually via wind pollination. Seeds drop close to the mother plant and usually form a long-term seed bank⁴⁷. Podospharea plantaginis (Castagne; U. Braun and S. Takamatsu) (Erysiphales, Ascomycota) is an obligate biotrophic powdery mildew that infects only P. lanceolata and requires living host tissue through its life cycle⁴⁸. It completes its life cycle as localized lesions on host leaves, only the haustorial feeding roots penetrating the leaf tissue to feed nutrients from its host. Infection causes significant stress for host plant and may increase the host mortality³¹. The interaction between P. lanceolata and P. plantaginis is strain-specific, whereby the same host genotype may be susceptible to some pathogen genotypes while being resistant to others⁴⁹. The putative resistance mechanism includes two steps. First, resistance occurs when the host plant first recognizes the attacking pathogen and blocks its growth. When the first step fails and infection takes place, the host may mitigate infection development. Both resistance traits vary among host genotypes⁴⁹.

Approximately 4000 P. lanceolata populations form a network covering an area of 50 × 70 km in the Åland Islands, SW of Finland. Disease incidence (0/1) in these populations has been recorded systematically every year in early September since 2001 by approximately 40 field assistants, who record the occurrence of the fungus P. plantaginis in the local P. lanceolata populations³⁰. At this time, disease symptoms are conspicuous as infected plants are covered by white mycelia and conidia. The coverage (m²) of P. lanceolata in the meadows was recorded between 2001 and 2008 and is used as an estimate of host population size. In the field survey two technicians estimate Plantago population size by visually estimating how much ground/other vegetation P. lanceolata foliage covers (m²) in each meadow. The proportion of P. lanceolata plants in each population suffering from drought is also estimated annually in the survey. Data on average rainfall (mm) in July and August was estimated separately for each population using detailed radar-measured rainfall (obtained by Finnish Meteorological Institute) and it was available for years 2001–2008.

Host population connectivity (S^H)²⁷ for each local population i was computed with the formula that takes into account the area of host coverage (m²) of all host populations surveyed, denoted with (A_j), and their spatial location compared to other host populations. We assume that the distribution of dispersal distances from a location are described by negative exponential distribution. Under this assumption, the following formula (1) quantifies for a focal population i, the effect of all other host populations, taking into account their population sizes and how strongly they are connected through immigration to it:

here, d_ij is the Euclidian distance between populations i and j and 1/α equals the mean dispersal distance, which was set to be two kilometres based on results from a previous study¹⁶.

The annual survey data has demonstrated that P. plantaginis infects annually 2–16% of all host populations and persists as a highly dynamic metapopulation through extinctions and re-colonizations of local populations¹⁶. The number of host populations has remained relatively stable over the study period⁴⁹. The first visible symptoms of P. plantaginis infection appear in late June as white-greyish lesions consisting of mycelium supporting the dispersal spores (conidia) that are carried by wind to the same or new host plants. Six to eight clonally produced generations follow one another in rapid succession, often leading to local epidemic with substantial proportion of the infected hosts by late summer within the host local population. Podosphaera plantaginis produces resting structures, chasmothecia, that appear towards the end of growing season in August–September³¹. Between 20% and 90% of the local pathogen populations go extinct during the winter, and thus the recolonization events play an important role in the persistence of the pathogen regionally¹⁶.

Inoculation assay: Effect of connectivity and disease history on phenotypic disease resistance

Host and pathogen material for the experiment

To examine whether the diversity and level of resistance vary among host populations depending on their degree of connectivity (S^H) and disease history, we selected 20 P. lanceolata populations for an inoculation assay. These populations occur in different locations in the host network, and were selected based on their connectivity values (S ^H of selected populations was 37–110 in isolated and 237–336 in highly connected category, Fig. 1). We did not include host populations in the intermediate connectivity category that was used in the population dynamic analyses in the inoculation assay due to logistic constraints. Podosphaera plantaginis is an obligate biotrophic pathogen that requires living host tissue throughout its life cycle, and obtaining sufficient inoculum for experiments is extremely time and space consuming. In both isolated and highly connected categories, half of the populations (IDs 193, 260, 311, 313, 337, 507, 1821, 1999, 2818 and 5206) were healthy during the study years 2001–2014, while half of the populations (IDs 271, 294, 309, 321, 490, 609, 1553, 1556, 1676 and 1847) were infected by P. plantaginis for several years during the same period. We collected P. lanceolata seeds from randomly selected ten individual plants around the patch area from each host population in August 2014.

To acquire inoculum for the assay, we collected the pathogen strains as infected leaves, one leaf from ten plant individuals from four additional host populations (IDs 3301, 4684, 1784, and 3108) in August 2014. None of the pathogen populations were same as the sampled host populations and hence, the strains used in the assay all represent allopatric combinations. Both host and pathogen populations selected for the study were separated by at least two kilometres. The collected leaves supporting infection were placed in Petri dishes on moist filter paper and stored at room temperature until later use.

Seeds from ten mother plants from each population were sown in 2:1 mixture of potting soil and sand, and grown in greenhouse conditions at 20 ± 2 °C (day) and 16 ± 2 °C (night) with 16:8 L:D photoperiod. Due to the low germination rate of collected seeds, population 260 (isolated and healthy population) was excluded from the study. Seedlings of ten different mother plants were randomly selected among the germinated plants for each population (n = 190), and grown in individual pots until the plants were eight weeks old.

The pathogen strains were purified through three cycles of single colony inoculations and maintained on live, susceptible leaves on Petri dishes in a growth chamber 20 ± 2 °C with 16:8 L:D photoperiod. Every two weeks, the strains were transferred to fresh P. lanceolata leaves. Purified powdery mildew strains (M1–M4), one representing each allopatric population (3301, 4684, 1784 and 3108), were used for the inoculation assay. To produce enough sporulating fungal material, repeated cycles of inoculations were performed before the assay.

Inoculation assay quantifying host resistance phenotypes

In order to study how the phenotypic resistance of hosts varies depending on population connectivity and infection history, we scored the resistance of 190 host genotypes, ten individuals from each study populations (n = 19), in an inoculation assay. Here, one detached leaf from each plant was exposed to a single pathogen strain (M1–M4) by brushing spores gently with a fine paintbrush onto the leaf. Leaves were placed on moist filter paper in Petri dishes and kept in a growth chamber at 20 ± 2 with a 16/8D photoperiod. All the inoculations were repeated on two individual Petri plates, leading to 760 host genotype—pathogen genotype combinations and a total of 1520 inoculations (19 populations * 10 plant genotypes * 4 pathogen strains * 2 replicates). We then observed and scored the pathogen infection on day 12 post inoculation, under dissecting microscope. The resulting plant phenotypic response was scored as 0 = susceptible (infection) when mycelium and conidia were observed on the leaf surface, and as 1 = resistance (no infection), when no developing lesions could be detected under a dissecting microscope. A genotype was defined resistant only if both inoculated replicates showed similar response (1), and susceptible if one or both replicates became infected (0).

Statistical analyses

Bayesian spatio-temporal INLA model of the changes in host population size

To study how the pathogen infection influences on host population growth, we analyzed the relative change in host population size (m2) (defined as population size (t) − population size (t−1))/population size (t−1)) between consecutive years utilizing data from 2001 to 2008 in response to pathogen presence-absence status at t−1 (Supplementary Table 2). To assess whether this depends on host population connectivity, we estimated the separate effects of pathogen presence/absence in the previous year for connectivity categories—high-, low, and intermediate—that were based on the 0.2 and 0.8 quantiles of the host-connectivity values (Fig. 1A and Supplementary Figs. 1, 2). This allowed us to directly assess and compare the effect of the pathogen on host population growth in the extreme categories between isolated and highly connected host populations which were represented in the sampling for the inoculation study (Fig. 2).

As covariates, we included the proportion (0–100%) of dry host plants measured each year within each local population as well as data on the amount of rainfall at the summer months (June, July, and August) obtained from the satellite images, as these were suggested be relevant for this pathosystem in an earlier analysis¹⁶. Observations where the change in host population size, or the host population coverage had absolute values larger than their 0.99 quantiles in the whole data, were regarded as outliers and omitted from the analysis. Before the analyses, all the continuous covariates were scaled and centred, and the categorical variables were transformed into binary variables.

The relative changes in local host population size between consecutive years was analyzed by a Bayesian spatio-temporal statistical model that simultaneously considers the effects of a set of biologically meaningful predictors. The linear predictor thus consists of two parts (2,3):

where (beta) represents the correlation coefficients corresponding to the effects of environmental covariates, ({z}_{t}) corresponds to the spatiotemporal random effect, and ({X}_{t}) and ({A}_{t}) project these to the observation locations. For ({z}_{t}) we assume that the observations from a location in consecutive time points (t−1) and t are described by 1st order autoregressive process:

where ({w}_{t}) corresponds to spatially structured zero-mean random noise, for which a Matern covariance function is assumed. Statistical inference then targets jointly the covariate effects (beta), the temporal autocorrelation (varphi), and the hyperparameters describing the spatial autocorrelation in w_t. From these the overall variance, as well as spatial range—a distance after which spatial autocorrelation ceases to be significant—can be inferred (Supplementary Fig. 3). For more detailed description of the structure of the statistical model and how to do efficient inference with it using R-INLA, we refer to refs. 16,50.

Identification of resistance phenotypes

The phenotype composition of each study population was defined by individual plant responses to the four pathogen strains, where each response could be “susceptible = 0” or “resistant = 1”. For example, a phenotype “1111” refers to a plant resistant to all four pathogen strains. The diversity of distinct resistance phenotypes within populations was estimated using the Shannon diversity index as implemented in the vegan software package⁵¹. The Shannon diversity index for all four study groups was then analyzed using a linear model with class predictors population type (well-connected or isolated), infection history (healthy or infected), and their interaction.

Analysis of population connectivity and infection history effects on host resistance

To test whether host population resistance varied depending on connectivity (S^H) and infection history, we analyzed the inoculation responses (0 = susceptible, 1=resistant) of each host-pathogen combination by using a logit mixed-effect model in the lme4 package⁵². The model included the binomial dependent variable (resistance-susceptible; 1/0), and class predictors population type (well-connected or isolated), infection history (healthy or infected), mildew strain (M1, M2, M3, and M4) and their interactions. Plant individual and population were defined as random effects, with plant genotype (sample) hierarchically nested under population. Model fit was assessed using chi-square tests on the log-likelihood values to compare different models and significant interactions, and the best model was selected based on AIC-values. P-values for regression coefficients were obtained by using the car package⁵³. We ran all the analyses in R software⁵⁴.

The metapopulation model

We model the ecological and co-evolutionary dynamics of host and pathogen metapopulations to understand key features of the experimental system that impact on the qualitative patterns observed. The structure and parameters in our model are therefore not estimated using experimental data, but rather are chosen to cover a range of possibilities (e.g., low vs high transmission rates, variation in trade-off shapes for fitness costs). We construct the metapopulations in two stages to account for relatively well and poorly connected demes. All demes are identical in quality (i.e., no differences in intrinsic birth or death rates between demes) and only differ in their connectivity. Our metapopulation consists of an outer network of 20 demes, equally spaced around the unit square (0.2 units apart), and a 7×7 inner lattice of demes at a minimum distance of 0.2 units from the outer network (Fig. 3A), giving a total of 69 demes. Demes that are separated by a Euclidean distance of at most 0.2 are then connected to each other. This means that populations near the centre of the metapopulation are highly connected, while those on the boundary of the metapopulation are poorly connected. This also has the effect of making connections between well and poorly connected demes assortative (i.e., well/poorly connected demes tend to be connected to well/poorly connected demes). We relax the assumption of assortativity in a second type of network by randomly reassigning connections between demes, while maintaining the same degree distribution. (i.e., the probability of two demes being connected is proportionate to their degree). While well connected demes still have more connections to other well connected demes than to poorly connected demes, they are not more likely to be connected to a well connected deme than by chance based on the degree distribution. In both types of network structure, we classify a deme as well-connected if it is in the top 20% of the degree distribution and poorly connected if it is in the bottom 20%.

We model the genetics using a multilocus gene-for-gene framework with haploid host and pathogen genotypes characterized by (L) biallelic loci, where 0 and 1 represent the presence and absence, respectively, of resistance and infectivity alleles. Host genotype (i) and pathogen genotype (j) are represented by binary strings: ({x}_{i}^{1}{x}_{i}^{2}ldots {x}_{i}^{L}) and ({y}_{j}^{1}{y}_{j}^{2}ldots {y}_{j}^{L}). Resistance acts multiplicatively such that the probability of host (i) being infected when challenged by pathogen (j) is ({Q}_{{ij}}={sigma }^{{d}_{{ij}}}), where (sigma) is the reduction in infectivity per effective resistance allele and ({d}_{{ij}}={sum }_{k=1}^{L}{x}_{i}^{k}big(1-{y}_{j}^{k}big)) is the number of effective resistance alleles (i.e., the number of loci where hosts have a resistance allele but pathogens do not have a corresponding infectivity allele). Hosts and pathogens with more resistance or infectivity alleles are assumed to pay higher fitness costs, ({c}_{H}left(iright)) eq. (4) and ({c}_{P}left(jright)) eq. (5) with:

and

where (0 , < , {c}_{H}^{1},; {c}_{P}^{1},le, 1) control the overall strength of the costs (i.e., the maximum proportional reduction in reproduction (hosts) or transmission rate (pathogens)) and ({c}_{H}^{2},; {c}_{P}^{2}in {{mathbb{R}}}_{ne 0}) control the shape of the trade-off. When ({c}_{H}^{2},; {c}_{P}^{2}, < , 0) the costs decelerate (increasing returns) and when ({c}_{H}^{2},; {c}_{P}^{2}, > , 0) the costs accelerate the costs accelerate (decreasing returns) (Supplementary Fig. 4). This formulation, therefore, allows for a wide-range of trade-off shapes that may occur in nature.

The dynamics of the (finite) host and pathogen populations are modelled stochastically using the tau-leap method with a fixed step size of (tau=1). For population (p), the mean host birth rate at time (t) for host (i) (6) is

where (a) is the maximum per-capita birth rate, (q) is the strength of density-dependent competition on births, ({N}_{p}left(tright)={S}_{i}^{p}left(tright)+{I}_{icirc }^{p}left(tright)) is the local host population size, ({S}_{i}^{p}left(tright)) and ({I}_{icirc }^{p}left(tright)={sum }_{j=1}^{n}{I}_{{ij}}^{p}left(tright)) are the local sizes of susceptible and infected individuals of genotype (i), and ({I}_{{ij}}^{p}left(tright)) is the local size of hosts of genotype (i) infected by pathogen (j). Host mutations occur at an average rate of ({mu }_{H}) per loci (limited to at most one mutation per time step), so that the mean number of mutations from host type (i) to ({i}^{{prime} }) is ({mu }_{H}{m}_{i{i}^{{prime} }}{B}_{i}^{p}left(tright)), where ({m}_{i{i}^{{prime} }}=1) if genotypes (i) and ({i}^{{prime} }) differ at exactly one locus, and is 0 otherwise.

The mean local mortalities for susceptible and infected individuals are (b{S}_{i}^{p}left(tright)) and (left(b+alpha right){I}_{{ij}}^{p}left(tright)), respectively, where (b) is the natural mortality rate and (alpha) is the disease-associated mortality rate. The average number of infected hosts that recover is (gamma {I}_{{ij}}^{p}left(tright)), where (gamma) is the recovery rate.

The mean number of new local infections of susceptible host type (i) by pathogen (j) eq. (7) is:

where (beta) is the baseline transmission rate and ({Y}_{j}^{p}left(tright)) is the local number of pathogen propagules following mutation and dispersal. Pathogen mutations occur in a similar manner to host mutations, with mutations from type (j) to ({j}^{{prime} }) occurring at rate ({mu }_{P}{m}_{j{j}^{{prime} }}{I}_{circ j}^{p}left(tright)) where ({mu }_{P}) is the mutation rate per loci (limited to at most one mutation per timestep) and ({I}_{circ j}^{p}left(tright)={sum }_{i=1}^{n}{I}_{{ij}}^{p}left(tright)) is the local number of pathogen (j.) Following mutation, the local number of pathogens of type (j) eq. (8) is:

Pathogen dispersal occurs following mutation at a rate of (rho) between connected demes, given by the adjacency matrix ({G}_{{pr}}), with ({G}_{varSigma p}) the total number of connections for deme (p). The mean local number of pathogen propagules following mutation and dispersal eq. (9) is therefore:

We focus our parameter sweep on: (i) the structure of the network (assortative or random connections); (ii) the strength (left({c}_{H}^{1},; {c}_{P}^{1}right)) and shape (left({c}_{H}^{2},; {c}_{P}^{2}right)) of the trade-offs; (iii) the transmission rate (left(beta right)); and (iv) the dispersal rate (left(rho right)), fixing the remaining parameters as described in Supplementary Table 1 (preliminary investigations suggested they had less of an impact on the qualitative outcome) and conducting 100 simulations per parameter set. For each simulation we initially seed all populations with the most susceptible host type and place the least infective pathogen type in one of the well-connected populations to minimize the risk of early extinction. We then solve the dynamics for 10,000 time steps (preliminary investigations indicated this was a sufficient period for the metapopulations to reach a quasi-equilibrium in terms of overall resistance). We calculate the average level of resistance (proportion of loci with a resistance allele) between time steps 4001 and 5000 (transient dynamics) and over the final 1000 time steps (long-term dynamics) for well and poorly connected demes, categorized according to whether the disease is present in (infected) or absent from (uninfected) the local population at a given time point and discarding simulations where the pathogen is driven globally extinct.

We compare the mean level of resistance in infected/uninfected poorly/well-connected populations across all simulations to the empirical results. We say that a simulation is a qualitative ‘match’ for the empirical findings if: (i) in poorly connected demes, the infected populations are on average at least 5% more resistant than uninfected populations; and (ii) in well-connected demes, the uninfected populations are on average at least 5% more resistant than infected populations. In other words, if ({R}_{{CS}}) is the mean resistance for a population with connectivity (C) ((C=W) and (C=P) for well and poorly connected demes, respectively) and infection status (S) ((S=U) and (S=I) for uninfected and infected populations, respectively), then a parameter set is a qualitative ‘match’ for the empirical findings if ({R}_{{WU}} > 1.05{R}_{{WI}}) and (1.05{R}_{{PI}}, > , 1.05{R}_{{PU}}). If these criteria are not met, then the parameter set is a qualitative ‘mismatch’ for the empirical findings. The model is not intended to be a replica of an empirical metapopulation, but rather is used to reveal the key factors which lead to qualitatively similar distributions of resistance and disease incidences observed in the study of the Åland islands. Hence, the purpose of the model is to determine which biological factors are likely to be crucial to the patterns observed herein.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com