Synthesizing selection mosaic theory and host-pathogen theory to explain large-scale pathogen coexistence

Domain: chemistry
Published: 2025-12-26
Editorial outcome: minor_revision

Abstract
Selection mosaic theory explains observations of polymorphism in host-pathogen interactions in terms of spatially variable natural selection but does not account for population dynamics. In contrast, classical host-pathogen theory easily explains observations of population cycles, but does not explain the persistence of pathogen polymorphism. Here, we synthesize these two frameworks to understand the effects of population cycles on pathogen polymorphism. We show that geographic variation in the frequency of two morphotypes of a baculovirus that infects the Douglas-fir tussock moth (Orgyia pseudotsugata) depends on the frequency of Douglas-fir (Pseudotsuga menziesii), an important tussock moth host tree. The morphotype frequency data are best explained by host-pathogen models that combine a selection mosaic with population cycles. In our model, population cycles intensify pathogen competition across a selection mosaic, leading to a strong effect of Douglas-fir frequency on morphotype frequency that matches the data. Models without host-pathogen cycles or a selection mosaic project only weak effects of varying Douglas-fir frequency. Our model further projects that a biopesticide made up of both viral morphotypes would be more effective than the current single-morphotype biopesticide, demonstrating that our synthesis of selection mosaic theory and host-pathogen theory provides useful insights into pest management.

Introduction
A central principle from the theory of host-pathogen interactions of host-pathogen ecology and host-pathogen evolution is based largely on theory rather than on data. The many host-pathogen systems that show high levels of pathogen polymorphism suggest that evolutionary dynamics are indeed likely to affect host-pathogen ecology, and have motivated efforts to incorporate evolutionary theory into wildlife disease management. The theory also states that ecological dynamics should, in turn, have strong effects on host and pathogen evolutionary dynamics, but empirical evidence for such effects is extremely rare. As a result, our understanding of the interacting effects of host-pathogen ecology and host-pathogen evolution is based largely on theory rather than on data. Another basic principle of host-pathogen theory is that pathogens can easily regulate populations of their hosts, either by maintaining the host’s population at a stable point equilibrium, or by causing the host’s population to undergo cycles or chaotic fluctuations. Field studies of spillover events and emerging diseases have indeed shown that pathogens often suppress their host’s population, while longer-term field studies have shown that some parasites and pathogens do indeed cause fluctuations in their host’s population. Addressing the knowledge gap between evolutionary and ecological and evolutionary and ecological studies of hosts and pathogen thus requires a consideration of the complex fluctuations in density that may result from interactions between hosts and their pathogens.

A lack of consideration of ecological dynamics is particularly a problem in studies of pathogen coexistence. Although some pathogen coexistence theory considers ecological mechanisms, by far the most widely cited theory invokes genotype-by-genotype (G × G) interactions. In G × G interactions, pathogens are specialized to attack particular host strains, while hosts are specialized to defend against particular pathogen strains. A key extension of this work allows for environmental variation across landscapes, thus yielding G × G × E theory. G × G × E theory has shown that pathogen coexistence can result from spatially heterogeneous selection pressure, a type of spatial structure known as a “selection mosaic”.

Although selection mosaics have been documented in many different organisms (Supplementary Note 1), to the best of our knowledge, the empirical evidence for effects of selection mosaics and G × G × E interactions on host-pathogen systems has come entirely from studies of plant pathogens and not animal pathogens. The immobility of plants suggests that spatial structure is indeed likely to affect plant-pathogen interactions, whereas the high mobility of many animals suggests that dispersal may eliminate the effects of spatial structure on animal-pathogen interactions, so that G × G × E theory may be irrelevant for animal-pathogen interactions. Moreover, G × G × E theory almost invariably assumes that host and pathogen densities are constant. The intense spatial competition that plants experience supports this assumption (but see Papaïx et al.), but as we have described, animal hosts of many infectious pathogens undergo fluctuations in density.

Whether selection mosaics and G × G × E interactions play a role in animal-pathogen interactions is therefore unclear. To understand how ecological dynamics affect pathogen polymorphisms, and to test the relevance of selection mosaics for hosts with pathogen-driven population cycles, we studied a widespread polymorphism in a fatal, specialist baculovirus of a forest pest insect, the Douglas-fir tussock moth (Orgyia pseudotsugata). Populations of the Douglas-fir tussock moth periodically increase from densities at which the insect is practically undetectable to densities at which forest damage is widespread and severe. Because peak insect populations are regularly decimated by baculovirus epizootics (Mihaljevic et al., epizootics are epidemics in animals), and because host-pathogen models can reproduce the 7–10 year period and 3–5 order of magnitude amplitude of the insect’s population cycle, the baculovirus appears to be causing predator-prey type cycles in the insect.

The Douglas-fir tussock moth’s population cycles may also be affected by genetic variation in the baculovirus. The Douglas-fir tussock moth baculovirus consists of two morphotypes, a multi-capsid morphotype in which the virions or “nucleocapsids” are arranged in bundles within infectious particles or “occlusion bodies”, and a single-capsid morphotype, in which the nucleocapsids are arranged singly within occlusion bodies. The two morphotypes are more closely related to other Lepidopteran baculoviruses than to each other, and have long co-occurred across much of the range of the Douglas-fir tussock moth in North America, including the province of British Columbia, Canada, and the states of Washington, Oregon, Idaho, Montana, Colorado, California, Nevada, Arizona, and New Mexico in the USA. This long-term record of co-occurrence suggests that the two morphotypes coexist, but the mechanisms that allow for this coexistence are unknown.

Because many host-pathogen models predict that pathogen coexistence is unlikely, it is not immediately obvious why the morphotypes are able to coexist, but G × G × E theory may provide the answer. Like most insect baculoviruses, the Douglas-fir tussock moth baculovirus is transmitted when host larvae consume foliage contaminated with occlusion bodies released from infectious cadavers of other larvae. This is important because larvae feed on both Douglas-fir and Abies firs, and because an insect’s risk of baculovirus infection can be affected by variation in foliage quality, structure, nutrient content, and secondary metabolites, or by effects of foliage variation on an insect’s mobility, health, and feeding rate, anyone of which can affect infection risk. The complex differences in foliage chemistry between Douglas-fir and Abies firs, therefore, suggest that a tussock moth’s infection risk, and thus the fitness of the two morphotypes, might differ between Douglas-fir and Abies firs.

Here, we ask, can population cycles and selection mosaics explain the pattern of coexistence of the two morphotypes of the Douglas-fir tussock moth baculovirus? Our approach to answering this question combines methods from studies of plant pathogens with methods from studies of animal pathogens. We follow studies of plant pathogens in using a field experiment to quantify the fitness of the two morphotypes on two important Douglas-fir tussock moth host trees, Douglas-fir, Pseudotsuga menziesii, and grand fir, Abies grandis. Our field experiment allows us to directly test whether a selection mosaic driven by variation in host diet has the potential to affect morphotype coexistence. We follow studies of animal pathogens in using observational field data to choose between host-pathogen models with and without population cycles and with and without selection mosaics. Our model allows us to test whether morphotype coexistence can be explained by a combination of population cycles and a selection mosaic (Fig. 1).

Our model also allows us to consider how Douglas-fir tussock moth populations are affected by the use of the baculovirus in pest control. Because the Douglas-fir tussock moth damages valuable timber, and because its baculovirus infects only a few insect species, between 1975 and 1992 the U.S. Forest Service produced large quantities of the baculovirus for use as a microbial biopesticide to control Douglas-fir tussock moth populations. This biopesticide, known as “Tussock Moth BioControl-1” (TMB-1), consists only of the multi-capsid morphotype, which preliminary research suggested is slightly more virulent than the single-capsid morphotype. Given that the two morphotypes coexist in nature, an unanswered question is, how are tussock moth population dynamics affected by the use of a biopesticide that consists only of a single baculovirus morphotype?

Our models track competition between morphotypes across forests with varying relative abundance of Douglas-fir versus other tree species. The models include genetic variation in both the host and the pathogen, and allow the level of host variation to differ between pathogen strains, thereby incorporating a G × G interaction. Because the models allow the G × G interaction to depend on forest composition, and because forest composition varies spatially, the models incorporate G × G × E interactions.

Host and pathogen variation across morphotypes and across forests of different tree species composition is incorporated in terms of differences in average pathogen infectiousness, represented by the average per-capita transmission rate ν, and in terms of differences in variation in pathogen infectiousness, represented by the coefficient of variation of the transmission rate C. The average transmission rate ν can equivalently represent either pathogen infectiousness or host infection risk, averaged across hosts, while variation in infectiousness C can equivalently represent either variation in pathogen infectiousness or variation in host infection risk. Higher ν always leads to higher host infection risk and thus higher pathogen fitness, but higher C leads to lower pathogen fitness. This is because increasing C while holding ν constant is equivalent to adding more-susceptible hosts and more-resistant hosts in such a way that the mean is unaffected even though the C.V. increases; more-resistant hosts have a disproportionate impact on infection risk, however, and so pathogen strains with higher values of C have lower fitness.

Although here we do not decompose our estimates of host variation into genetic and environmental components, insect resistance to baculoviruses has long been known to be heritable, and in the baculovirus of the closely related spongy moth, Lymantria dispar, variation in infectiousness C has been shown to be heritable. We therefore allow for heritable variation in our models.

Importantly, variation in infectiousness plays a key role not just in the selection mosaic in our models, but also in the host-pathogen population cycles generated by the models. When variation in host infection risk C > 1, as it is known to be in this and other insect-baculovirus systems, and when this variation is entirely due to environmental causes and is thus not heritable, then host and pathogen densities reach a stable point equilibrium. Even in the presence of high environmental stochasticity, models producing a stable point equilibrium produce only short-period, small-amplitude cycles, which fail to capture the long-period, large-amplitude cycles seen in the Douglas-fir tussock moth and other forest defoliators. In contrast, eco-evolutionary models that assume that host variation is at least partly heritable produce realistic population cycles even for realistically high values of variation in host infection risk C.

This effect is important partly because it means that the heritable host variation that is fundamental to host-pathogen selection mosaics is also fundamental to host-pathogen population cycles. In host-pathogen population cycles, when host variation in infection risk is high, a significant fraction of the host population will always be relatively resistant to infection. If host variation is not heritable, the ability of resistant hosts to avoid infection will prevent the pathogen from imposing the high levels of density-dependent mortality that drive population cycles. If host variation is heritable, and if there is a fitness tradeoff between host infection risk and host fecundity, then high pathogen densities will drive the evolution of both higher host resistance and lower host fecundity, strengthening the effects of density-dependence and allowing for population cycles even when host variation is high. As we will show, this effect interacts with selection mosaics in complex ways.

To synthesize selection-mosaic theory and host-pathogen theory, we construct a novel two-pathogen eco-evolutionary model that incorporates natural selection on host resistance and pathogen-driven population cycles. To evaluate the usefulness of the model, we first use field experiments to show that two key transmission parameters differ between the two morphotypes of the Douglas-fir tussock moth baculovirus and that the extent of these differences depends on the host tree species consumed during infection. We then test the model’s ability to explain variation in morphotype coexistence, as quantified by a combination of published data and our own extensive dataset of morphotype frequencies. This model explains the morphotype frequency data better than models without selection on host resistance or without population cycles, demonstrating that there is a selection mosaic determined by geographic variation in forest composition and that this mosaic is modulated by pathogen-driven population cycles. Our model further shows that the morphotype composition of viral biopesticides can influence host-pathogen evolution, such that mixtures containing equal proportions of the two morphotypes provide the most effective insect population control.

Results
Geographic variation in morphotype frequency
To assess the usefulness of G × G × E theory for explaining the long-term polymorphism of the two baculovirus morphotypes, we first quantified their relative frequencies across the range of the Douglas-fir tussock moth. Published data on morphotype frequencies included mostly small sample sizes, and so we collected a large number of new samples. Sampling from the entire range was beyond what we could accomplish, and insects can only be found during outbreaks, so most of our samples were collected during a period of widespread outbreaks in Washington, Oregon, Idaho, and Montana USA (Supplementary Table 3). The combined dataset reveals that the two morphotypes co-occur at high frequencies across much of the insect’s range (Fig. 2). In 61% of populations in which only one morphotype was recorded, most had small sample sizes, such that only 14% (11 total) had a sample size greater than 5 larvae, with an average of 4.27 larvae. In contrast, among the 39% of populations in which both morphotypes were recorded, 46% had a sample size greater than 5 larvae, with an average of 36 larvae. This pattern in the data is consistent with the widespread co-occurrence of the two morphotypes at high frequency.

A comparison of morphotype frequencies across years shows that the polymorphism has persisted in roughly its current form for at least 50 years (Supplementary Table 4). The data thus imply that the two morphotypes have coexisted for many decades.

To test for effects of forest tree species composition on morphotype frequencies, we fit generalized linear models (GLMs) to the data, and we used Akaike information criterion (AIC) to compare the fit of models with and without effects of Douglas-fir frequency. This analysis showed that the model with an effect of Douglas-fir frequency provided a much better explanation for the data, such that ΔAIC = 134.9 (Supplementary Note 2). We thus conclude that variation in the frequency of the multi-capsid morphotype across sites is strongly related to variation in the relative percent Douglas-fir across sites. This statistical result does not identify the biological mechanism that allows for coexistence or for dynamical changes in the relative frequency of the morphotypes. It nevertheless suggests an experimentally testable hypothesis, namely that the fitness of the two morphotypes varies across the range of the Douglas-fir tussock moth because of variation in the frequency of Douglas-fir relative to other tree species.

Variation in morphotype fitness
To estimate fitness differences between the two morphotypes on different host trees, we carried out a field experiment in which we allowed uninfected Douglas-fir tussock moth larvae to feed on branches of either Douglas-fir or grand fir that were contaminated with infectious cadavers. The cadavers were infected with one of five multi-capsid (MNPV) baculovirus isolates or one of five single-capsid (SNPV) baculovirus isolates (Fig. 1). We defined an isolate to be a sample of virus from one infected insect collected in the field. This experiment allowed us to quantify the two key components of pathogen fitness described above, average infectiousness ν and variation in infectiousness C.

To analyze our transmission data, we embedded our mechanistic transmission model (Fig. 1) in a Bayesian hierarchical model that we fit to the experimental data. This approach showed that our estimates of pathogen fitness components varied by viral morphotype and host tree species. When insects were exposed to the virus while feeding on Douglas-fir branches, the multi-capsid morphotype had slightly lower average infectiousness ν than the single-capsid morphotype, but much lower variation in infectiousness C (Fig. 3A). When insects were feeding on grand fir, the multi-capsid morphotype instead had substantially higher ν and slightly lower C than the single-capsid morphotype. To explain the implications of these experimental parameter estimates for morphotype fitness, ultimately, we will use long-term models of morphotype fitness in forests of varying composition. As a simple preliminary indicator of fitness, first, we use the cumulative fraction of hosts infected by each morphotype in a single epizootic in forests composed of a single tree species.

When pathogen transmission is strongly density-dependent, as it is in the case of baculoviruses, pathogen fitness increases sharply with initial host density, and with increases in a pathogen’s average infectiousness. The theory of pathogen variation has further shown that pathogen fitness is sharply reduced by increases in a pathogen’s variation in infectiousness; these negative effects, however, are stronger at higher densities than at lower densities, whereas the positive effects of increases in the mean are similar across densities. If two pathogen strains differ with respect to both their mean infectiousness and their variation in infectiousness, then one pathogen strain may have an advantage at low host densities, while the other has an advantage at high host densities.

These effects are important because our two morphotypes differ in both their mean infectiousness and in their variation in infectiousness, and the direction of these differences in turn differs between host trees. To understand the importance of such effects for our system, we incorporated our experimental parameter estimates into an SEIR model, and we used the model to project changes in morphotype fitness with density. Here we (temporarily) define fitness as the cumulative fraction infected at the end of an epizootic, and we (also temporarily) assume single-pathogen epizootics in single-host-tree forests. Because the multi-capsid morphotype has slightly lower average infectiousness but much lower variation in infectiousness on Douglas-fir, in Douglas-fir-only forests our model projects that the fitness of the multi-capsid morphotype will be slightly lower than the fitness of the single-capsid morphotype at low densities, but slightly higher at high densities (Fig. 3B). Because the multi-capsid morphotype has moderately higher average infectiousness but slightly lower variation in infectiousness on grand fir, in grand fir-only forests the model projects that the fitness of the multi-capsid morphotype will be higher than the fitness of the single-capsid morphotype at all densities, but the difference will be substantially less than in Douglas-fir-only forests.

Our experiments thus provided initial evidence that the fitnesses of the two morphotypes of the Douglas-fir tussock moth baculovirus differ between the two host tree species, suggesting that a selection mosaic may indeed explain variation in morphotype frequency across forests of different tree species composition. The interacting effects of host tree species and insect density on morphotype fitness meant that we could not easily intuit how morphotype frequency would vary with the frequency of Douglas-fir. Our next step was to construct long-term models of pathogen competition to understand the implications of our experiment for variation in morphotype frequency across forest types. Because of the complex effects of host density in Fig. 3B, it was at least possible that allowing for host-pathogen population cycles would alter morphotype fitness, and so our models allowed for the possibility of population cycles.

To test whether variation in morphotype frequency across forests could be explained by variation in morphotype fitness across tree species, we implemented our two-pathogen, eco-evolutionary models on a spatial grid representing a forest. On the grid, the tree species at each point was either Douglas-fir or grand fir (Fig. 1). We used grand fir as a proxy for species other than Douglas-fir because grand fir was the most prevalent Abies host in our tree species dataset, and was thus the most prevalent tussock moth host tree species besides Douglas-fir. We assumed that the model forest encompassed a scale of 8076.6 m2. Because hatching larvae can disperse on the wind for distances of 400 m or more, this scale allowed larvae to easily disperse to any location in the forest. This assumption ensured that pathogen fitness would not be unduly influenced by the clumping of cadavers within a forest; given that our ultimate interest is in whether G × G × E theory is useful for understanding polymorphism in animal pathogens, this was a conservative assumption. As we explain in Supplementary Note 8, the scale of our model forest is appropriate to the scale at which tree species vary across the range of the Douglas-fir tussock moth.

We then tested whether our morphotype frequency data are best explained by models that include selection mosaics or by models that do not include selection mosaics. To test whether host population cycles are necessary to explain the data, we also included a non-evolutionary model without heritable variation. We chose the best model using AIC.

Host-pathogen population cycles and a selection mosaic together explain variation in morphotype frequency
To test whether models with selection mosaics explained the morphotype frequency data better than models without selection mosaics, we compared models in which morphotype variation in infectiousness C differed between the two host tree species to models in which variation in infectiousness did not differ between host tree species. Because variation in host plant foliage quality can have complex effects on insect fecundity, we also allowed for a mosaic in the host’s resistance-fecundity tradeoff by assuming that the tradeoff parameter differed between tree species. Because the tradeoff parameter modulates average morphotype infectiousness ν, allowing the value of s for each morphotype to vary between tree species leads to a selection mosaic in average morphotype infectiousness.

Model selection confirmed that the best-fitting model includes selection mosaics driven by variation in tree species with respect to both variation in morphotype infectiousness and with respect to the resistance-fecundity tradeoff parameter, and thus in average morphotype infectiousness (Table 1). Models that did not include both selection mosaics provided a much worse fit to the data. The model without either selection mosaic attempted to explain the morphotype frequency data by invoking high levels of stochasticity, leading to a much worse fit to the data (Fig. 4).

Notably, a model that included a selection mosaic in the resistance-fecundity tradeoff parameters but not in variation in infectiousness C provided a much better explanation for the data, emphasizing the importance of a selection mosaic in average morphotype infectiousness. This effect occurs because the tradeoff parameters can have a strong effect on the ability of the insect to develop high pathogen resistance, and thus low infection risk ν. In comparison, variation in infection risk C has the more complex effect of modulating how infection risk rises with insect density. Model comparison, however, show that both selection mosaics are necessary to explain the data.

Because the morphotype frequency data do not include any information about host population density, it was not obvious in advance that models with heritable host variation, and thus population cycles, would do a better job of explaining the morphotype frequency data than a model without heritable host variation (Remember that the level of host variation differs between pathogen morphotypes, and therefore depends on both the host and the pathogen). Indeed, like the evolutionary models, the non-evolutionary model includes a selection mosaic in terms of variation in infectiousness C. Because our experimental estimate of variation C > 1 for each morphotype on each tree species, the non-evolutionary model, Model 7 in Table 1, projects that the host population will undergo short-period, small-amplitude fluctuations around the equilibrium of the corresponding deterministic model. These fluctuations do not resemble the long-period, large-amplitude cycles observed in Douglas-fir tussock moth populations (Supplementary Fig. 12). For the model that instead assumes that host variation is (partly) heritable, host and pathogen populations fluctuate with a realistically long period and large amplitude. The extremely poor fit of Model 7 thus suggests that explaining the data requires models that allow for realistic host-pathogen population cycles. In Supplementary Fig. 14, we confirm this point by showing that versions of our eco-evolutionary model that show only short-period, short-amplitude cycles also provide a poor fit to the data.

Understanding the interacting effects of selection mosaics and population cycles
To explain why we need both a selection mosaic and population cycles to reproduce our morphotype frequency data, in Fig. 5 we show how the model dynamics change with changes in forest tree species composition. As Fig. 5B shows, the two pathogens together drive predator-prey type cycles in the density of the host, with a period and an amplitude that match the 8–10 year period and 3–5 order of magnitude fluctuations typical of Douglas-fir tussock moth populations in nature.

In the model, host infection risk ν, which is equivalent to pathogen infectiousness, changes along with host and pathogen densities, due to the balancing (but fluctuating) selection that results from the tradeoff between host infection risk and host fecundity. (Note that host fecundity depends on host infection risk through the risk-fecundity tradeoff, and therefore changes with changes in host infection risk/pathogen infectiousness.) Reassuringly, the average infectiousness of each morphotype in the model fluctuates around values that are close to our experimental estimate of the infectiousness of each morphotype (Fig. 5C), even though we did not fit the model to the experimental data. This latter result suggests that the infectiousness of each morphotype in our model was close to the infectiousness of each morphotype in nature. As the figure further shows, these fluctuations are strongly affected by the frequency of Douglas-fir relative to other tree species.

Crucially, as the relative frequency of Douglas-fir increases, the higher average infectiousness of the multi-capsid morphotype on Douglas-fir causes its relative frequency to change from being high only some of the time to being high most of the time (Fig. 5C, D). Meanwhile, because the morphotype frequency data were collected in populations with at least moderately high infection rates, in comparing the model to the data, we only included model generations in which the overall fraction of infected hosts was above 0.3, which usually only occurs just after peaks in the host population. This is important because the strong density-dependence of the transmission process causes epizootics near population peaks to be dominated by one pathogen morphotype or the other (Fig. 5D, E, F). Increases in the relative frequency of Douglas-fir therefore lead to sharp increases in the frequency of the multi-capsid morphotype, which is why the model provides a good fit to our data (Fig. 4). Notably, the fit of the best model to the data is almost exactly the same if we follow Dwyer et al. in allowing for a generalist predator or parasitoid (Supplementary Note 9). Also, the scales at which forest patches differ in tree species composition are typically far more than 10 km; at such scales, dispersal between dissimilar forest patches has no effect on our results (Supplementary Fig. 18). Our focus on single forest patches is therefore appropriate to our data.

In Supplementary Fig. 15, we show that when the effect of density-dependent transmission is stronger, so that the amplitude of the population cycles is larger, then as the frequency of Douglas-fir increases, the transition between domination by the single-capsid morphotype and domination by the multi-capsid morphotype is much sharper than it is in the data. Conversely, when the effects of density-dependent transmission are weaker, so that the amplitude of the population cycles is smaller, then as the frequency of Douglas-fir increases, the transition between domination by the single-capsid morphotype and domination by the multi-capsid morphotype is much more gradual than it is in the data. Reproducing the data thus requires that the effects of density-dependent transmission and population cycles be in rough balance with the effects of the selection mosaic, so that the transition between single-capsid domination and multi-capsid domination is neither too steep nor too gradual compared to the data.

Explaining the relationship between forest composition and our morphotype frequency data thus requires a combination of host-pathogen population cycles and a G × G × E interaction. The lack of this effect in the non-evolutionary model helps to explain why that model provides a poor fit to the data. Also, we note that the cycles in our data are driven both by the pathogen and by environmental stochasticity, demonstrating that stochasticity plays a key role in the host-pathogen population cycles.

Projecting the effects of biopesticide use on Douglas-fir tussock moth populations
The interacting effects of host population fluctuations and a selection mosaic in our model together suggest that the long-term use of a single morphotype in the TMB-1 baculovirus biopesticide product may alter the population dynamics of the Douglas-fir tussock moth. To understand such effects, we modified our best-fit model to include the use of the baculovirus as a biopesticide under three different scenarios: (1) a multi-capsid-only formulation like TMB-1; (2) a hypothetical single-capsid-only formulation; (3) a hypothetical 50:50 morphotype mixture formulation. In the model, we applied the biopesticide as an input of baculovirus at the beginning of the epizootic period (summer), such that the biopesticide was applied when the insect density N in generation n increased from the previous generation (N n > N n−1) and crossed a defined threshold N , such that N > N . To consider different levels of population control, we varied both the density threshold and the amount of biopesticide being applied. Because the biopesticide was expensive to produce and because supplies are limited, in practice, the biopesticide is not applied every time the insect density reaches a high level, so we also varied the likelihood of deploying the biopesticide in high-density years. We then projected the average insect population over 100 years of biopesticide application, and we used our projections to calculate the percent reduction in Douglas-fir tussock moth densities compared to the baseline no-biopesticide case.

We focus on the effect of varying the morphotype composition of the biopesticide and the probability of biopesticide application by fixing the density threshold N at 0.5 insects/m2 and by fixing the amount of biopesticide applied at the equivalent of 150 cadavers/m2, which were values that yielded low insect densities in a previous model of biopesticide use in a similar insect-virus system, and that yielded low insect densities in our model. Varying these parameters showed that lowering the threshold for biopesticide use and increasing the amount of biopesticide applied produced better control, as one might expect (Supplementary Note 10).

As one might further expect, more frequent biopesticide applications due to a higher probability of application also produced better population control, but the extent of this latter effect varied strongly with forest tree species composition and with the morphotype composition of the biopesticide. In forests with ≤ 80% Douglas-fir, all biopesticide formulations and all levels of biopesticide application reduced the host population, with the extent of the reduction becoming more severe as the probability of biopesticide application increased (Fig. 6).

Beyond this broad trend, however, there were also interesting differences between biopesticides of different compositions. In general, the hypothetical single-capsid-only biopesticide and the hypothetical 50:50 mixture biopesticide produced the best control, in that increasing the probability of applying these formulations produced ever-lower population densities, from a reduction of 30% to a reduction of more than 80% (Fig. 6). The current multi-capsid-only biopesticide in contrast produced only a modest improvement in control as the probability of application increased, with the population reduction only increasing with the probability of application from about 20% to about 30%. As we show in Supplementary Fig. 22, the multi-capsid-only biopesticide has lower efficacy because the host is able to evolve high resistance to the multi-capsid morphotype. More surprisingly, in forests with 95% Douglas-fir, a low probability of biopesticide application produced a worse outcome than no biopesticide at all. This counterintuitive result arises because infrequent biopesticide use sharply increases insect resistance without killing enough insects to significantly lower the insect’s density.

Methods
Geographic morphotype frequency data
Our morphotype frequency dataset combines published data sources and our own field collections (Supplementary Table 4), which together include 1914 samples across 128 locations. The first published source is Hughes, who sampled 65 locations across the range of the Douglas-fir tussock moth, with most locations in Washington, Oregon, and Idaho. To locate the sampled sites from the map included in Hughes, we used ArcGIS to determine the approximate latitudes and longitudes of the locations on the map. Because the Hughes data do not provide sample sizes, we assumed that the sample size at each location was one larva when one morphotype was present and two larvae when both morphotypes were present. This approach allowed us to take into account this dataset without allowing it to strongly influence our estimates of morphotype frequencies.

The second published source in our morphotype frequency dataset is Williams and Otvos, who collected 185 baculovirus samples in five geographic regions of southern British Columbia, Canada, by rearing larvae from field-collected egg masses. Their samples consisted only of the multi-capsid morphotype. The third published source is Williams et al., who collected 165 virus samples from 29 locations across the range of the Douglas-fir tussock moth. These samples came from both field-collected egg masses and larval cadavers.

Our addition to the data came from 6091 larvae that we collected at 22 locations in Washington, Oregon, Idaho, and Montana in the summers of 2019 and 2020. The larvae we collected were reared in the lab until death or pupation, yielding 1715 virus-infected cadavers. The overall fraction infected was 28%, although the fraction infected ranged from 0 to 60% across sites. Four locations did not have any virus infections, so they could not be used in our analyses of morphotype frequencies. The multi-capsid-only biopesticide (TMB-1) is sufficiently expensive and is in sufficiently limited supply that it is only ever applied to very small areas, and there are no records that it was ever used at any of our sample sites. All fieldwork and data collection in this study were carried out with the appropriate permits in coordination with the USDA Forest Service.

To identify the morphotype in each infected larva, we carried out PCR amplification of the polyhedrin gene polh, a conserved protein that occludes virions during environmental transmission (see Supplementary Table 2). We were able to successfully classify 1316 (76.7%) virus-infected cadavers as being infected with either the single-capsid morphotype, the multi-capsid morphotype, or both. Of these, 41.8% were infected with only the multi-capsid morphotype, 47.2% cadavers were infected with only the single-capsid morphotype, and 11.2% were infected with both morphotypes. Multiple lines of evidence indicate that infection with one morphotype did not affect an insect’s chance of subsequent infection with the other morphotype (Supplementary Fig. 1). We therefore recorded each coinfected insect as both a single-capsid and a multi-capsid isolate, so that our data consisted of 1463 morphotype samples.

Prior to the development of our PCR protocol, the corresponding author and colleagues used electron microscopy to morphotype 12 samples that had been collected by forest managers from across the insect’s range. Electron microscopy is much more expensive than PCR, but it makes it possible to directly visualize whether a virus-killed insect has occlusion bodies with virions arranged singly, as in the single-capsid morphotype, or arranged in bundles, as in the multi-capsid morphotype.

For each of our 123 sample sites in the United States, we estimated the relative frequency of Douglas-fir using the National Forest Type Dataset provided by the USDA Forest Service Forest Inventory and Analysis Program & Remote Sensing Applications Center. The dataset was downloaded on August 22, 2022. We used the extract function in the raster v3.6.26 R package to quantify relative Douglas-fir frequency within a 5 km radius around each site. For each of the Canadian sites, we used the Canadian Forest Monitoring MODIS kNN map on Canada’s National Forest Inventory website (https://nfi.nfis.org/en/maps) to confirm that the site was dominated by Douglas-fir. See Supplementary Note 2 for details on the tree species dataset.

Quantifying morphotype fitness with a field experiment
Overall infection risk depends on natural feeding behavior, which cannot be easily mimicked using laboratory experiments. We therefore used a field experiment to estimate average pathogen infectiousness ν and variation in pathogen infectiousness C for each morphotype on each host tree. We carried out this field experiment in the Okanogan-Wenatchee National Forest in Washington State, USA (47.315°N, 120.425°W) in 2021. Our experimental protocol was refined in two initial pilot experiments. Because the first of these pilot experiments had very low infection rates, while the second had high contamination levels in controls, we did not use the data from either pilot experiment.

In the first step in our experiment, we reared uninfected Douglas-fir tussock moth larvae in the laboratory from field-collected egg masses. The egg masses were collected in the spring of 2021 from a natural population in the Caribou-Targhee National Forest in Idaho, USA (43.75°N, 111.29°W). To ensure that the larvae in the experiment were not infected at hatch, we surface-sterilized all egg masses in 10% formalin, which has been shown to effectively eliminate viruses on insect eggs without killing the eggs.

In the second step in the experiment, we produced infected larvae by feeding neonates an artificial diet inoculated with a virus dose that we had determined using an initial dose-response experiment (Supplementary Note 3). We reared the neonates in the lab for 5 days at 25.5°C to ensure that they were near death. We then deployed the infected neonates on branches in the field.

In the third step in the experiment, we placed the infected neonates on small branches of the host trees, and we encased the branches in mesh bags to keep the neonates from escaping. The mesh bags were made of spun-bonded polyester, a translucent material that blocks UV, thereby eliminating the breakdown of the virus while permitting the passage of air and water. The larvae were thus allowed to feed under natural conditions. This type of experiment provides estimates of transmission rates that are close to estimates from naturally occurring epizootics. Because using multiple pathogen densities allows for estimation of variation in pathogen infectiousness, we used 5, 10, and 20 infected cadavers per branch to estimate variation in infectiousness.

After 5 days, a period long enough to ensure that the infected neonates died naturally on the foliage, we added ten fourth instars to each bag. To ensure that the fourth instars were uninfected, we reared them in the lab before deployment in the field. The resulting larval densities were within the range of larval densities observed in nature. We then allowed the uninfected fourth instars to feed on the contaminated foliage for 7 days, a period long enough to ensure that infection rates would reach at least moderate levels, but short enough to ensure that none of the initially uninfected fourth instars would become infected and die of the virus while in a bag. Our experiment thus included only a single round of transmission.

The branches in our experiment were on Douglas-fir and grand fir trees. We used five viral isolates of each morphotype, 9 of which were taken from our morphotype collections (see Supplementary Table 5 for the location from which each isolate was collected). The tenth isolate was the biopesticide TMB-1 that consists only of the multi-capsid morphotype.

After 7 days, we removed the branches to the lab, where we placed the fourth instars in individual 2-ounce cups containing cubes of uncontaminated artificial diet. We reared the larvae for 21 days, a period long enough to ensure that infected larvae would die of the virus. Baculovirus occlusion bodies can be seen at 400X under a light microscope, so all dead larvae were autopsied to determine if they were virus-killed. We estimated total branch area by photographing each branch at the end of the experiment and then by using the open-source image-analysis software ImageJ v1.53a to estimate the area. We then used the area of each branch and the number of infected larvae added to each branch to estimate initial pathogen density in terms of cadavers/m2.

Our field experiment included 92 experimental branches, 48 for the multi-capsid morphotype and 44 for the single-capsid morphotype, and 8 control branches with no infectious cadavers. Our goal in including control branches was to test whether there was any baculovirus on the branches besides the virus that we added, and to test whether any of our uninfected insects had accidentally become infected in the lab. To ensure that the control branches differed from the experimental branches only in having no added infectious cadavers, we enclosed the control branches in mesh bags at the same time that we enclosed the experimental branches (plus infected larvae) in mesh bags. Out of the 73 control larvae, only one became infected after being deployed in the field. This low rate (1.3%) indicates a lack of systematic contamination in either the lab or the field, and we therefore did not include the infected control larva in our analyses. A modest number of the initially uninfected fourth instars were dead when we opened the bags, but we confirmed via autopsy that none were infected with the baculovirus. We used the fraction infected from each experimental branch to fit the transmission parameters ν and C.

To determine the extent to which baculovirus fitness is determined by tree species and/or viral morphotype, we used model selection to compare the ability of different models to explain our experimental data, such that different models made different assumptions about whether transmission in the experiment depended on morphotype, morphotype and tree species, or neither morphotype nor tree species. Our models were based on a simplified version of our Susceptible-Exposed-Infectious-Recovered (SEIR) model. To produce this model, we adapted the SEIR model to our experiments by assuming that there is only one round of virus transmission and that the mesh bags eliminate the breakdown of the virus, as in our experiments. The density of infectious cadavers, therefore, did not change during the experiment, so we can set the rate of change of the pathogen density in our SEIR model to zero. The fraction of hosts that became infected during the experiment is then:

I(T)/S(0) = 1 - (1 + αC^2 ν P(0) T)^(-1/C^2)

Here, I(T) represents the number of infected insects that become infected over T days, the number of days for which the susceptible insects in the experiment were exposed to the baculovirus (in practice T = 7). The symbols S(0) and P(0) represent the initial density of susceptible insects and infectious cadavers, respectively. The fraction infected during the experiment is therefore I(T)/S(0). The average pathogen infectiousness/host infection risk is ν, and the variation in infectiousness/risk is C. Because first instar larvae are much smaller than fourth instar larvae, we include the parameter α to represent the ratio of the number of occlusion bodies produced by first instars to the number produced by fourth instars. Here we use an estimate α = 0.032 from previous work.

Because of genetic similarities, it is reasonable to assume that virus isolates within a morphotype are not independent. One rigorous way to account for this dependence is by using a Bayesian hierarchical model in which isolate is nested within morphotype. When we estimated the morphotype fitness parameters ν and C, we thus assumed that each isolate consisted of a random draw from a morphotype-specific distribution of each parameter.

To test whether baculovirus fitness is determined by tree species and/or viral morphotype, we compared the fit of models in which the fraction infected in the experiment depended on morphotype, morphotype and tree species, or neither morphotype nor tree species. In addition to models in which each isolate consisted of a random draw from a morphotype-specific distribution of parameters, we also considered a model in which each isolate consisted of a random draw from a single distribution that represented both morphotypes at once and that therefore did not include effects of morphotype on pathogen fitness. To choose between models, we used the leave-one-out cross-validation (LOO-CV) model-selection criterion; because the best models had very similar LOO-CV scores, we chose the overall best model using the minimum sum of squared errors (SSE, Supplementary Table 6). We implemented our Bayesian hierarchical models using the Stan statistical computing language v2.32.2, and we analyzed the output using the rstan v2.32.6, loo v2.7.0, bayesplot v1.11.1, and tidyverse v2.0.0 R packages (see Supplementary Note 3 for the fit of the models to the data and other model details). For purposes of comparison, we also considered a model that allowed for fitness differences between morphotypes but not tree species, and a model that included no fitness differences between morphotypes. We then parameterized our mechanistic models of morphotype competition using the parameter estimates from the best model of our experimental data.

To understand how our field estimates of the transmission parameters affect individual morphotype fitness over the course of a single epizootic, we used our SEIR model to quantify the fraction infected at the end of a season over a range of host densities. To show the range of the dynamics, we used 225 draws from the thinned posterior distribution of parameter values for the best model, which were subsetted using the as_draws_matrix and subset_draws functions in the posterior v1.5.0 R package. We simulated the SEIR model in an all Douglas-fir forest and in an all grand fir forest for each morphotype independently and calculated the cumulative fraction infected at the end of the epizootic.

Modeling the effects of variation in forest tree species on morphotype frequency
To test the extent to which our experimental transmission data could help explain patterns in our morphotype frequency data, we fit spatial, long-term models to the morphotype frequency data. In these models, we used an SEIR-type model to describe combined epizootics of the two viral morphotypes (Supplementary Note 4). At the end of the epizootic in the model, surviving hosts reproduce, and the two viral morphotypes experience over-winter mortality. Because we allow for a resistance-fecundity tradeoff, the average infectiousness of each morphotype changes after the epizootic due to host reproduction and the (incomplete) heritability of infection risk. The eco-evolutionary model thus describes changes in the density of the host, in the density of cadavers infected with each morphotype, and in the average infectiousness of each morphotype:

N(n+1) = e^(ϵ_n) N(n) (1 - i_1 - i_2) [r + r (s_1 ν*_1,n + s_2 ν*_2,n)]

Z(1,n+1) = ϕ N(n) i_1 + γ Z(1,n) (Infectiousness_Morph1)

Z(2,n+1) = ϕ N(n) i_2 + γ Z(2,n) (Infectiousness_Morph2)

Here, host density N(n+1) depends on the host density in the preceding generation N(n), the baseline fecundity r, a stochastic term ϵ_n, the fraction surviving the epizootic 1 - i_1 - i_2, and the resistance-fecundity tradeoff parameters s_1 and s_2 for pathogens 1 and 2, respectively. The symbols i_1 and i_2 represent the respective fraction of hosts infected with each pathogen. The variable ϵ_n is a normal random variate with mean 0 and standard deviation σ, representing “environmental stochasticity”, the randomness due to fluctuations in environmental variables such as weather. To allow for a cost of resistance in the host, we assume that fecundity increases with increasing average infectiousness, such that s_1 and s_2, the tradeoff parameters, are the respective rates at which fecundity increases with infection risk with respect to each pathogen. The symbols ν*_1 and ν*_2 represent the respective average infectiousness of each morphotype at the end of the previous epizootic, which may equivalently be interpreted as the infection risk of the host with respect to each morphotype. Note that pathogen-induced mortality always causes host infection risk at the end of an epizootic to be lower, and often much lower, than it was at the beginning of the epizootic. We therefore have ν*_1,n < ν_1,n and ν*_2,n < ν_2,n. Changes in host density are therefore partially driven by balancing selection, such that higher mean infection risks ν_1 and ν_2 lead to increased mortality from each morphotype but also to increased fecundity.

The state variables Z_1,n and Z_2,n represent the densities of the two pathogens in generation n. The symbol ϕ represents the effective overwintering rate of the pathogen, which incorporates both the survival of the pathogen from one generation to the next and the higher susceptibility of neonates to later instars. It is often true that ϕ > 1 because neonates have much higher susceptibility. Longer-term pathogen survival is represented by the survival rate γ.

The average infectiousness ν_1,n+1 for pathogen 1 in generation n+1 is modulated by the heritability of infectiousness b_1 for pathogen 1 and the correlation in the host’s susceptibility to each pathogen ρ. Although theory has shown that a negative correlation ρ < 0, corresponding to negative frequency-dependent selection, can strongly promote pathogen coexistence, in practice, correlations in an insect’s susceptibility to different baculovirus strains are almost invariably positive, and so negative frequency-dependent selection is unlikely to play a role in our system.

The parameter combinations b_1^2 C_1^2 and b_2^2 C_2^2 represent the fraction of overall variation that is due to additive genetic factors, so that higher values of b_1 and b_2 strengthen the effects of selection. Changes in infection risk are thus determined by balancing selection, as in the case of host density, except that, in contrast to host density, the change in infection risk does not depend on the baseline fecundity r.

To represent spatial variation in tree species within the forest in the model, we implemented equations on a hexagonal spatial grid. A key consideration in spatial modeling is how the grid structure interacts with the spatial domain, particularly at the edges. Hexagonal grids offer advantages over square grids because they have a symmetric, orthogonal coordinate system, a low perimeter-to-area ratio, and because the center of each hexagon is equidistant from the centers of its six nearest neighbors. These properties help reduce bias from edge effects. In contrast, square grids often have stronger edge effects because of their higher perimeter-to-area ratio and because they have two different kinds of nearest neighbors; four orthogonal neighbors that share an edge, and four diagonal neighbors that share only a corner. This difference creates asymmetries in local interactions.

We incorporated variation in the fraction of Douglas-firs in the hexagonal grids in our models by varying the fraction of hexagons that represented Douglas-fir versus grand fir. Our sample sites spanned a wide range of Douglas-fir frequencies, so in the model, we varied the number of Douglas-fir cells in the 37-cell grid from 1 (≈2.7%) to 36 (≈97.3%). Although our forest composition dataset includes sites with 0 and 100% Douglas-fir, we excluded these extreme cases from model fitting because including them led the model fitting routine into regions of parameter space in which the host and pathogen populations would frequently go extinct. Moreover, the extreme cases are probably not ecologically realistic, because forests defoliated by the Douglas-fir tussock moth are not likely to contain exactly 0% or 100% Douglas-fir. This limitation stems from the forest-types dataset, which records only the dominant tree species at a location; because not all species were recorded at each location, values of 0% or 100% Douglas-fir likely reflect sampling error.

Because it was possible that Douglas-fir tussock moth dispersal could eliminate the effects of the spatially varying selection that is fundamental to selection mosaics, it was important to allow for dispersal in our model. The only meaningful dispersal of larvae within a forest occurs when first instar larvae balloon on silk strands, which rarely leads to dispersal distances of more than a few hundred meters. Also, the vast majority of ballooning takes place before larvae begin feeding, so we do not expect host tree effects on ballooning. Ballooning distances in our model therefore follow an exponential distribution, so that most larvae disperse to their nearest neighbors, but some larvae disperse much farther, and dispersal is identical on the two host-tree species (Supplementary Note 5). We set the dispersal parameter in our model so that the average distance that larvae dispersed was 2.69 grid units. If we assume that the average distance that larvae disperse in the field is approximately 50 m, then the distance between the centers of two adjacent grid cells is approximately 18.59 m and the distance between two grid points at opposite edges is approximately 111.5 m. Given the geometry of our spatial grid, the total area of our forest is 8076.6 m2.

To test whether our morphotype frequency data can be best explained by variation in host tree species and/or population cycles, we constructed models that allowed for selection mosaics in the variation in infectiousness C and the fecundity-resistance tradeoff parameter s, and we compared them to models that had no selection mosaic in either C or s. For models that allowed for a selection mosaic in C (Models 1 and 5 in Table 1), we used estimates of C from the best transmission model that allowed for host tree differences. For the models that did not allow for a selection mosaic in C or s (Model 4 in Table 1), we used parameter estimates from the best transmission model that did not allow for effect of tree species on C. For the model that allowed a selection mosaic for s but not C (Model 3 in Table 1), we used estimates of C from the transmission model that did not allow for effects of tree species on C. For the model that included a selection mosaic for s, but no selection mosaic or morphotype effects on C (Model 2 in Table 1), we used estimates from the transmission model that did not allow for either host tree or morphotype effect on C. For the model that did not include host evolution, but did allow C and ν to vary spatially (Model 7 in Table 1), we used estimates of C and ν from the transmission model in which C and ν both varied between morphotypes and between host tree species.

To simplify the fitting process, we used estimates of the remaining parameters of the SEIR model from previous work with colleagues in which SEIR models were fit to data from baculovirus epizootics in insect populations (see Supplementary Table 9). We then estimated five parameters from the morphotype frequency data: the virus overwintering parameter ϕ, the correlation between host infection risk to the two pathogens ρ, the fecundity-resistance tradeoff parameters, and the stochasticity parameter σ, which is the standard deviation in the stochasticity term ϵ_n in equations (2). We used maximum likelihood to estimate these parameters from our morphotype frequency data.

Our model fitting approach involved generating stochastic realizations of the model and using the model output to calculate the average fraction of larvae that were infected by the multi-capsid morphotype in years in which the overall fraction infected was greater than or equal to 0.3. We then computed the likelihood of observing the morphotype frequency data given the average fraction infected in the model. To align the morphotype frequency data with the model’s spatial representation of forest structure, each field site was assigned to the hexagonal grid with the most similar forest composition. Because the morphotype frequency data consist of the proportion of each morphotype observed at each location, the data are fundamentally binomial, and so we used a binomial likelihood function. To account for model stochasticity, we averaged the likelihoods across many model realizations, yielding an integrated likelihood estimate.

Because generating model realizations is highly computationally intensive, we used a simple likelihood maximization approach known as “line search”. The line search algorithm works by iteratively optimizing one parameter at a time across a predetermined range of values. The target parameter is then set to its optimal value across its range, and the process is repeated for the next parameter. The algorithm cycles through the parameters until each parameter has been varied a specified number of times.

We first implemented our line search algorithm by averaging likelihoods across 96 model realizations, while fixing the stochasticity parameter σ = 0.5. Each parameter was sampled independently in the line search routine with no nesting for tree species or morphotype. We then optimized σ for each of the 30 best parameter sets from the line search step. In this second optimization step, we improved the accuracy of our likelihood estimates by averaging likelihoods across 2000 realizations. We then chose the best model using the AIC model-selection criterion. All model output in the figures is based on the best-fitting parameter set from each model (Supplementary Note 6).

To understand how the eco-evolutionary dynamics change as a function of forest composition, we simulated the model over 2000 generations using the best-fit parameters. To consider the effects of more complex ecological interactions, we also considered a model with a generalist predator (Supplementary Note 9). To understand the extent to which our results are affected by dispersal, we extended our best model to allow for two hexagonal grids with dispersal between grids (Supplementary Note 8).

Projecting the effects of long-term biopesticide use
To implement the biopesticide in our model, we followed the general approach of Reilly and Elderd, who tested the effects of a baculovirus biopesticide on outbreak cycles in the spongy moth. Reilly and Elderd found that the threshold density at which the biopesticide is applied and the amount applied have important effects on long-term pest dynamics, but that the timing of the biopesticide application within a season has little effect. We therefore applied biopesticides after the dispersal of first instars, which is typically when forest managers apply biopesticides in forest insect populations. The biopesticide was then applied in the model as an input of the baculovirus when the insect density N in generation n increased from the previous year N_n > N_n-1 and when the population crossed a threshold insect density N_T, so that N > N_T.

We varied four parameters that affect biopesticide use. First, we varied the threshold insect density that triggers an application, using values between 0.2 and 2.0 insects/m2. Second, we varied the amount of biopesticide applied, using values between 1 and 200 cadavers/m2. Third, we varied the relative morphotype compositions that make up the biopesticide, using either only the multi-capsid morphotype, only the single-capsid morphotype, or a mixture of the two. Fourth, we varied the probability of applying the biopesticide when the threshold and density criteria were met, using values between 0 and 1. For this, we fixed the threshold at 0.5 insects/m2 and the biopesticide amount at 150 cadavers/m2.

We simulated each scenario for 350 years/tussock moth generations, such that the biopesticide was only ever applied between years 150–250 and thus not before year 150 or after year 250. Eliminating biopesticide application before year 150 and after year 250 allowed us to compare the effects on the host population and average infectiousness before, during, and after consistent biopesticide treatment. We then calculated average host population densities across 20 model realizations.

All simulations and morphotype frequency model fitting routines were performed using the DifferentialEquations v7.4.0, Distributions v0.25.73, and MPI v0.20.8 packages in the Julia programming language v1.7.3 and OpenMPI v4.1.0 library project for the message passing implementation.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an exciting demonstration of how selection mosaics can maintain genetic diversity in pathogens through GxGxE interactions. The authors find that a well characterized pathogen polymorphism is best explained by variation in host habitat, which then influences pathogen fitness in a density-dependent manner. While such dynamics have been shown in plants and fungi, I am unaware of any other studies which have demonstrated a similar interaction in animal pathogens. This result should be of interest to disease ecologists, and evolutionary biologists more broadly, as how dynamical forces can maintain such polymorphisms is a major open question. However, I feel like these results are let down by the writing, which can be somewhat disorganized. Particularly, this manuscript would be more accessible if the authors frontload some of the exposition (especially with regards to infectiousness variation) currently in the results to the introduction. Furthermore, while the details necessary to replicate this work are available in the supplemental information, the model descriptions in the main text left me entirely confused about what was happening under the hood. If the authors can revise this manuscript so that the methods are more transparent, and some of the background is presented more clearly in the introduction, I believe this study would make a valuable contribution to Nature Communications.\nMajor Comments\nI found it difficult to understand the broad structure of the predictive model, it would be helpful to give more detail at the end of the introduction or the beginning of the results to explain the overall framework.\nToo much of the results section is spent explaining concepts which should be addressed in the introduction. Specifically, I believe it would make the paper more readable to address the effects of variation in infectiousness of pathogens in the introduction. While the effects of infectiousness heterogeneity as described in this manuscript are well supported by Fleming-Davies et al., 2015, I don’t feel this is sufficiently clear from the manuscript alone. To this point, the authors should better explain their equation for the proportion of infected individuals over time (I(T)/S(0) in figure 1). While the authors have derived and used this expression in previous papers, it is not common knowledge and should be further detailed.\nThe description of the eco-evolutionary model in the main text is entirely insufficient. While I understand the complexity of the model means some details must be tucked away in the supplemental information, the main text description offers very little specificity about their implementation. When the authors state that they include variation in both the host and pathogen (L182), what specifically are they varying other than v and C? Is there another host trait included in the model that affects transmission?\nMinor Comments\nL25-27: The challenge of limiting the use of “host-pathogen” in disease ecology manuscripts is a persistent, yet important for readability.\nL42: To say “A lack of consideration of host-pathogen ecology in studies of host-pathogen evolution is thus a major oversight” is a little strong, the disease ecology and evolution community has long recognized the importance of ecological variation to host-pathogen evolution.\nL103-104: This statement is too broad, standard models such as the gene-for-gene model and the matching alleles model often predict coexistence. While these models are often used for plant resistance, they have also seen applications for animal immune systems.\nL202: How is this tradeoff implemented? Does the model assume continuous variability in fecundity / resistance, and thus requiring a particular tradeoff function?\nL:298-299: I don’t like the terminology of “allowing host-pathogen population cycles” rather than incorporating aspects which lead to population cycles\nFig 5: The terminology “% of Fraction Infected” or “Fraction of fraction infected” is confusing. I also find the panels in G not particularly informative, is the vertical axis host population density?\nSupplemental Materials\nIn equations 4-5 are the k’s and deltas supposed to have different subscripts?\nL221: Typo: “an SEIR models” - > “an SEIR model”\nL261: Typo in “vriation”\nL305-311: Is the direction sampled from a uniform distribution?\nL321: Equations ?? should be clarified\nL350: Figure ?? should be clarified", "reviewer_2": "This paper contributes significantly to understanding host-pathogen interactions by integrating selection-mosaic theory with population cycles. The study provides empirical evidence on the fitness differences between two morphotypes of the Douglas-fir tussock moth baculovirus. It explores the effects of host tree species composition on their coexistence. The authors employ a combination of field experiments, statistical modelling, and theoretical frameworks to support their conclusions. The study is innovative in its attempt to bridge theories from plant-pathogen and animal-pathogen systems. However, certain key assumptions require further clarification, and additional discussion on alternative mechanisms and ecological variables would enhance the robustness of the study.\n1. The assumption that heritable host variation is necessary for population cycles to occur could benefit from further discussion. While this assumption aligns with previous eco-evolutionary models, it would be helpful to consider whether population cycles might also arise from purely ecological factors such as density-dependent effects or environmental fluctuations. Discussing cases where population cycles occur without heritable variation could strengthen the justification.\n2. The application of selection-mosaic theory to an insect-virus system is an interesting approach, but the potential impact of host mobility on selection mosaics may require further exploration. Unlike plants, insects can disperse across different environments, potentially weakening spatially structured selection pressures. Clarifying how this mobility influences the selection mosaic framework in this context would enhance the theoretical foundation.\n3. The role of occasional coinfections in influencing morphotype coexistence is briefly mentioned. While coinfections may be rare, they could introduce additional selection pressures that influence transmission dynamics and morphotype balance. If the authors have data on coinfection rates, incorporating this information would clarify whether coinfections significantly impact the system. Even if infrequent, a brief discussion of their possible evolutionary consequences would be valuable.\n4. The model does not explicitly account for ecological factors such as predators, parasitoids, and weather conditions, all of which can significantly influence host-pathogen interactions. A sensitivity analysis examining how these variables might affect model outcomes would enhance the study’s robustness. If previous studies on this insect-virus system have documented strong effects of these ecological variables, including a discussion of their potential impact on the model’s conclusions would be helpful.\n5. The paper discusses environmental variation in infection risk but does not fully explore how it interacts with evolutionary trade-offs between host susceptibility and pathogen infectivity. A more detailed discussion of how fluctuating environmental conditions might shift the balance of these trade-offs would be beneficial. It would also be useful to clarify whether the model assumes these traits to be static or if they vary dynamically over time.\n6. The finding that transmission risk decreases with higher Douglas-fir density is interesting but may require further discussion. Exploring possible behavioural or ecological mechanisms that might explain this pattern, such as changes in larval movement or feeding preferences, could strengthen this result.\n7. The study suggests that morphotype coexistence is driven by selection mosaics and population cycles. However, alternative explanations such as frequency-dependent selection, bet-hedging strategies, or trade-offs in dispersal ability could also play a role. If there are data or theoretical arguments that rule out these mechanisms, incorporating this discussion would provide a more comprehensive view of polymorphism maintenance.\n8. The paper states that surface sterilization with formalin is effective for virus removal. To ensure methodological rigour, control tests confirming that no residual contamination occurred would be beneficial. Additionally, the baculovirus doses used to inoculate neonates should ideally be validated through dose-response trials to confirm that they reflect naturally occurring infection levels. If similar dose-response validations have been conducted in prior studies of this virus, referencing them would strengthen the methodology.\n9. The use of mesh bags to enclose branches for field experiments is a practical approach, but it may have influenced microenvironmental conditions, potentially affecting viral persistence or larval behaviour. If previous studies using similar enclosures have addressed whether they impact host-pathogen interactions, referencing those findings would provide useful context. Additionally, considering the possibility of undetected low-level secondary transmission within the bags could add further transparency.\n10. The potential effects of host tree chemistry, such as differences in secondary metabolites or nutrient content, on infection dynamics should be acknowledged. If differences in foliage chemistry between Douglas-fir and grand-fir influence larval feeding behaviour or viral transmission, incorporating this information would strengthen the discussion. If direct data are unavailable, acknowledging the potential impact and suggesting future research directions would be helpful.\n11. The authors mention that selection mosaics have been extensively studied in plant-pathogen systems but not in animal-pathogen systems. If there are any studies on selection mosaics in animals, referencing them would provide a broader context and strengthen the study’s position within the literature.\n12. The model assumes that host infection risk variation is heritable. If this assumption has not been explicitly validated for this system, discussing its implications and potential limitations would be useful. If relevant, the authors might also consider referencing studies that have tested the heritability of similar traits in related systems.\n13. The policy implications for Douglas-fir tussock moth management are particularly interesting, especially the recommendation to use a single morphotype for biocontrol. However, could using only one morphotype lead to unintended evolutionary consequences, such as increased resistance or compensatory adaptation by the host? Addressing potential long-term effects would enhance the applied relevance of the study.\n14. In the supplementary material, Eq. 6 should be written as E_k1. In line 230, the mean should be 1/delta_1, and the variance should be 1/k*delta_1^2. Equation 7 needs a bracket around [r+r(s1v1+s2v2)]. In section 3.3 (line 323), the notations I1 and I2 should likely be lowercase. Please check the other equations for consistency, ensuring that notations are properly continued and all terms are clearly defined. The equation number is missing in line 321, and the figure number is missing in line 350.\n15. The study could benefit from a discussion on how heterogeneous dispersal patterns might influence pathogen spread and morphotype coexistence. If individuals disperse non-randomly based on host tree composition or pathogen load, how would this affect the selection-mosaic framework? Exploring whether heterogeneous dispersal could amplify or dampen selection mosaics would add depth to the analysis.\n16. The influence of heterogeneous node properties in the spatial network should be considered. If certain locations in the forest act as transmission hubs due to higher host density or favourable environmental conditions, how might this affect pathogen coexistence? Addressing whether the spatial model incorporates such heterogeneities would provide a clearer picture of the system’s dynamics."}, "author_response": "Response Letter, NCOMMS-24-74983-T, Synthesizing Selection-\nMosaic Theory and Host-Pathogen Theory to Explain Large-Scale\nPathogen Coexistence in an Insect-Virus System\nGeneral Comment\nIn an effort to thoroughly respond to the reviewer comments, we deleted, substituted or moved more than\n95% of the text in both the main text and in the Supplemental Information. We therefore did not indicate\nthe new text in either document. Instead, we now indicate key changes by referring to the relevant lines in\nbold face.\nReviewer 1\nOverall Comments\nThis is an exciting demonstration of how selection mosaics can maintain genetic diversity in\npathogens through GxGxE interactions. The authors find that a well characterized pathogen\npolymorphism is best explained by variation in host habitat, which then influences pathogen\nfitness in a density-dependent manner. While such dynamics have been shown in plants and\nfungi, I am unaware of any other studies which have demonstrated a similar interaction in\nanimal pathogens. This result should be of interest to disease ecologists, and evolutionary\nbiologists more broadly, as how dynamical forces can maintain such polymorphisms is a major\nopen question. However, I feel like these results are let down by the writing, which can be\nsomewhat disorganized. Particularly, this manuscript would be more accessible if the authors\nfrontload some of the exposition (especially with regards to infectiousness variation) currently\nin the results to the introduction. Furthermore, while the details necessary to replicate this\nwork are available in the supplemental information, the model descriptions in the main text\nleft me entirely confused about what was happening under the hood. If the authors can revise\nthis manuscript so that the methods are more transparent, and some of the background is\npresented more clearly in the introduction, I believe this study would make a valuable contri-\nbution to Nature Communications.\nWe have completely rewritten the manuscript. As part of this rewrite, we have moved key explanatory\nmaterial fromthe Results to theIntroduction at lines 121-154. Wealso movedthe eco-evolutionary model\nfrom the SI to the Methods in the main text at lines 633-714, with the model in equations (2)-(5). We\nkept the SEIR model in the SI, partly in the interests of brevity and partly because the SEIR model is less\ncentral to our overall argument. We also now explain the model at lines 116-130.\nMajor Comments\nI found it difficult to understand the broad structure of the predictive model, it would be\nhelpful to give more detail at the end of the introduction or the beginning of the results to\nexplain the overall framework.\nWe now explain the predictive model in detail in the Introduction at lines 121-154. In the Results, we\nnow explain the features of the model that strongly affect the model projections at lines 285-327 .\nToo much of the results section is spent explaining concepts which should be addressed in\nthe introduction. Specifically, I believe it would make the paper more readable to address\nthe effects of variation in infectiousness of pathogens in the introduction. While the effects of\ninfectiousness heterogeneity as described in this manuscript are well supported by Fleming-\nDavies et al., 2015, I don’t feel this is sufficiently clear from the manuscript alone. To this\npoint, the authors should better explain their equation for the proportion of infected individu-\nals over time (I(T)/S(0) in figure 1). While the authors have derived and used this expression\n1\nin previous papers, it is not common knowledge and should be further detailed.\nWe now explain the effects of pathogen variation in infection risk in the Introduction at lines 121-154.\nIn the Methods in the main text, at lines 590-604, we now explain how we derived the expression for\nI(T)/S(0) from our SEIR model and we include the expression itself as equation (1).\nThe description of the eco-evolutionary model in the main text is entirely insufficient. While\nI understand the complexity of the model means some details must be tucked away in the\nsupplemental information, the main text description offers very little specificity about their\nimplementation. When the authors state that they include variation in both the host and\npathogen (L182), what specifically are they varying other than v and C? Is there another host\ntrait included in the model that affects transmission?\nThis is an excellent point. At lines 116-130 in the main text, we now explain that our models track\ncompetition between morphotypes across forests that vary in their relative abundance of Douglas-fir versus\nother tree species. We further explain that our models allow for variation in the host and the pathogen in\nterms of average transmission ν¯ and variation in transmission C, as the reviewer states, and that average\ntransmission can represent either host infection risk or pathogen infectiousness. Variation in transmission\ncorrespondinglyrepresentsvariationineitherthehostorthepathogen. Wedidnotincludeotherhosttraits\nthat affect transmission, largely because to the best of our knowledge there are no other such traits.\nMinor Comments\nL25-27: The challenge of limiting the use of “host-pathogen” in disease ecology manuscripts\nis a persistent, yet important for readability.\nWe appreciate the author’s concern, and we have therefore deleted many uses of “host-pathogen”. We\nemphasize, however, that in some cases “host-pathogen” is clearer and briefer than the alternatives.\nL42: To say “A lack of consideration of host-pathogen ecology in studies of host-pathogen\nevolution is thus a major oversight” is a little strong, the disease ecology and evolution com-\nmunity has long recognized the importance of ecological variation to host-pathogen evolution.\nWe have rewritten this text to instead argue that our understanding of the effects of ecological change\non host-pathogen evolution is based largely on theory rather than data. The relevant sentence is at line 34\nin the Introduction.\nL103-104: This statement is too broad, standard models such as the gene-for-gene model and\nthe matching alleles model often predict coexistence. While these models are often used for\nplant resistance, they have also seen applications for animal immune systems.\nAt line 46 in the Introduction we now instead say that some models show how ecological factors can\ndrive coexistence, but that most models, including gene-for-gene models, focus on G×G interactions.\nL202: How is this tradeoff implemented? Does the model assume continuous variability in\nfecundity / resistance, and thus requiring a particular tradeoff function?\nThe tradeoff is implemented in terms of increasing fecundity with increasing host infection risk. We now\nclarify this point at line 150 in the Introduction.\nL:298-299: I don’t like the terminology of “allowing host-pathogen population cycles” rather\nthan incorporating aspects which lead to population cycles.\nFig5: Theterminology“%ofFractionInfected”or“Fractionoffractioninfected”isconfusing.\nI also find the panels in G not particularly informative, is the vertical axis host population\ndensity?\nIn fig. 5, We have replaced “% of Fraction Infected” with “Relative Frequency”, and we have eliminated\nthe panels in G.\nSupplemental Materials\nIn equations 4-5 are the k’s and deltas supposed to have different subscripts?\nYes. We have fixed this mistake in equations (3)-(5) at line 254 in the SI. We appreciate the Reviewer’s\nattention to detail.\nL221: Typo: “an SEIR models” → “an SEIR model”.\nAs part of our revision, this phrase was modified, and so this problem is no longer an issue.\nL261: Typo in “vriation”\nCorrected. We appreciate the reviewer’s attention to detail.\nL305-311: Is the direction sampled from a uniform distribution?\nBecause we are modeling the density of insects rather than individual insects, dispersal in our model is\nnot stochastic. Instead, the density of insects moving to other points is based on distance and not direction,\nand so dispersal is equal in all directions.\nL321: Equations ?? should be clarified.\nTo avoid this problem, we now simply refer to the evolutionary model without listing the equation num-\nbers.\nL350: Figure ?? should be clarified.\nWe replaced this sentence as part of the revision.\nReviewer 2\nOverall Comments\nThis paper contributes significantly to understanding host-pathogen interactions by integrat-\ning selection-mosaic theory with population cycles. The study provides empirical evidence on\nthe fitness differences between two morphotypes of the Douglas-fir tussock moth baculovirus.\nIt explores the effects of host tree species composition on their coexistence. The authors\nemploy a combination of field experiments, statistical modelling, and theoretical frameworks\nto support their conclusions. The study is innovative in its attempt to bridge theories from\nplant-pathogen and animal-pathogen systems. However, certain key assumptions require fur-\nther clarification, and additional discussion on alternative mechanisms and ecological variables\nwould enhance the robustness of the study.\nand to discuss alternative mechanisms and ecological variables.\n1. The assumption that heritable host variation is necessary for population cycles to oc-\ncur could benefit from further discussion. While this assumption aligns with previous eco-\nevolutionary models, it would be helpful to consider whether population cycles might also\narise from purely ecological factors such as density-dependent effects or environmental fluc-\ntuations. Discussing cases where population cycles occur without heritable variation could\nstrengthen the justification.\nAtlines 398-402intheDiscussion,wenowexplainthatDouglas-firtussockmothoutbreakcyclescould\nin theory be driven by specialist parasitoids or food limitation. We also emphasize, however, that the ev-\nidence that these alternatives play a key role in Douglas-fir tussock moth cycles is much weaker than the\nevidence that the baculovirus plays a key role.\n2. The application of selection-mosaic theory to an insect-virus system is an interesting ap-\nproach, but the potential impact of host mobility on selection mosaics may require further\nexploration. Unlike plants, insects can disperse across different environments, potentially\nweakening spatially structured selection pressures. Clarifying how this mobility influences the\nselection mosaic framework in this context would enhance the theoretical foundation.\nThis is such an important point that we now include a new subsection in the SI, subsection 4.4 at lines\n326-351 in the SI, in which we consider how our results would be affected by host dispersal between forest\npatches. Todothis,weextendedourmodeltoallowfortwopatches,wevariedthedistancebetweenthetwo\npatches, and we compared the results to the results from our single-patch model. To increase the chances\nthat we would see a strong effect, we set the frequency of Douglas-fir in one patch to 5%, and the frequency\nof Douglas-fir in the other patch to 95%.\nThe results of this model differ somewhat from the results of our original model, but only if the two\npatches are less than 400m apart. The lack of an effect at distances beyond 400m is important because\nof how the percent of Douglas-fir between forests in our data set varied with distance between forests. To\nshowthis, wenowincludeaplotthatexpressesthepairwisedifferenceinthefractionofDouglas-firbetween\nlocations as a function of the distance between the locations. The few locations that are less than 10km\\apartdifferalmostnotatallintheirpercentDouglas-fir, whilelocationsthatdifferstronglyintheirpercent\nDouglas-fir are all at least 10km apart, and are often more than 100km apart. We conclude that, in our\nsystem, dispersalisunlikelytohaveplayedanimportantroleindeterminingtheeffectsofforesttree-species\ncomposition on the relative frequency of the two baculovirus morphotypes.\n3. The role of occasional coinfections in influencing morphotype coexistence is briefly men-\ntioned. While coinfections may be rare, they could introduce additional selection pressures\nthat influence transmission dynamics and morphotype balance. If the authors have data on\ncoinfection rates, incorporating this information would clarify whether coinfections signifi-\ncantly impact the system. Even if infrequent, a brief discussion of their possible evolutionary\nconsequences would be valuable.\nIn the main text, we now emphasize more strongly that we present data on coinfections in the SI, in\nsubsection 2.1 at lines 62-79. These data strongly suggest that coinfections are not a key issue. Although\nthe reviewer is undoubtedly correct that a discussion of the evolutionary consequences of coinfections would\nbe interesting, the coinfection effects are so weak that we would not know what to say. In the interests of\nbrevity we therefore consider coinfections only briefly.\n4. The model does not explicitly account for ecological factors such as predators, parasitoids,\nand weather conditions, all of which can significantly influence host-pathogen interactions.\nA sensitivity analysis examining how these variables might affect model outcomes would en-\n4\nhance the study’s robustness. If previous studies on this insect-virus system have documented\nstrong effects of these ecological variables, including a discussion of their potential impact on\nthe model’s conclusions would be helpful.\nThis is a sufficiently important point that we have added a model that includes the effects of generalist\npredators, or equivalently of generalist parasitoids, in the SI at subsection 4.8, at lines 428-449. Previous\nwork by the last author and colleagues showed that generalist predators or parasitoids can increase the\nvariance in the outbreak period, but are unlikely to have much effect on the mean of of the outbreak period,\norontheamplitude. Totestwhetherthisgeneralizationholdsinourcase,wecreatedaversionofourmodel\nthat includes a generalist; as we now show in subsection 4.8 of the SI, the generalist does indeed have very\nlittle effect on our model’s results.\nIn the first author’s dissertation, of which the current paper was a chapter, we explored the effects of\nweather on the insect’s range, and we now cite the dissertation (the associated paper is in prep.) at lines\n432-436 in the Discussion. Meanwhile, systematic weather effects on insect outbreak cycles are not well\nunderstood, and there is no obvious way to allow for weather in our model, other than as the environmental\nstochasticity that we already include. We now say that this is an important area of future research, at line\n407 in the Discussion.\n5. The paper discusses environmental variation in infection risk but does not fully explore\nhow it interacts with evolutionary trade-offs between host susceptibility and pathogen infec-\ntivity. A more detailed discussion of how fluctuating environmental conditions might shift the\nbalance of these trade-offs would be beneficial. It would also be useful to clarify whether the\nmodel assumes these traits to be static or if they vary dynamically over time.\nWenowemphasizemorestronglythathostfecunditychangesovertime,atlines292-293andlines651-\n653. Wealsonowemphasizemorestronglythatourmodelincludesenvironmentalfluctuationsintheformof\nstochasticfluctuationsinfecundity,atlines 325-327andatlines 717-720. Atlines 737-738,wealsonow\nclarifythatwefitthevalueofstochasticitytothedata,andthusthatstochasticityisimportantinourmodel.\n6. The finding that transmission risk decreases with higher Douglas-fir density is interesting\nbut may require further discussion. Exploring possible behavioural or ecological mechanisms\nthat might explain this pattern, such as changes in larval movement or feeding preferences,\ncould strengthen this result.\nIn work that is in revision at The American Naturalist, we showed that the mechanisms underlying the\ndifferences in the two morphotypes in our work include both variation in the innate infectiousness of the\ntwo morphotypes, and variation in the host’s ability to avoid the two morphotypes. We now briefly summa-\nrize this work in the Discussion section at lines 423-430, and we include a citation to a draft on bioRxiv\n(Freedman et al. 2025).\n7. The study suggests that morphotype coexistence is driven by selection mosaics and population\ncycles. However, alternative explanations such as frequency-dependent selection, bet-\nhedging strategies, or trade-offs in dispersal ability could also play a role. If there are data\nor theoretical arguments that rule out these mechanisms, incorporating this discussion would\nprovide a more comprehensive view of polymorphism maintenance.\nIn the revision, we now explain that our model includes frequency-dependent selection. An important\npoint, however, is that, to allow for pathogen coexistence, frequency-dependent selection requires negative\ncorrelations in a host’s risk of infection to two or more pathogen strains. In the case of insect baculoviruses,\nprevious work has shown that correlations in host infection risk between morphotypes tend to be positive\nrather than negative. It therefore appears that frequency-dependent selection is unlikely to play a role in\nour system. We explain these points at lines 661-664.\nhave begun to feed, and thus before the baculovirus has begun to replicate. It is therefore also unlikely that\ndifferences in dispersal play a role in our system.\nIn the interests of brevity, we do not include a discussion of bet-hedging.\n8. The paper states that surface sterilization with formalin is effective for virus removal. To\nensuremethodologicalrigour, controltestsconfirmingthatnoresidualcontaminationoccurred\nwould be beneficial. Additionally, the baculovirus doses used to inoculate neonates should ide-\nally be validated through dose-response trials to confirm that they reflect naturally occurring\ninfection levels. If similar dose-response validations have been conducted in prior studies of\nthis virus, referencing them would strengthen the methodology.\nAt lines 582-584, We now emphasize more strongly that the very low infection rate in our control\ntreatments (1 out of 73 larvae) demonstrates that surface-sterilization was effective. Second, we did carry\nout a dose-response trial to quantify the effects of our doses on neonates. We now include the dose-response\ndata in a new subsection of the SI, subsection 3.1, pages 12-15, lines 142-178, with the results in SI figs. 7\nand 8.\n9. The use of mesh bags to enclose branches for field experiments is a practical approach,\nbut it may have influenced microenvironmental conditions, potentially affecting viral persis-\ntence or larval behaviour. If previous studies using similar enclosures have addressed whether\nthey impact host-pathogen interactions, referencing those findings would provide useful con-\ntext. Additionally, considering the possibility of undetected low-level secondary transmission\nwithin the bags could add further transparency.\nFirst, as we described in our response to the previous comment, previous work has provided strong that\nour experimental design produces realistic transmission rates. We now emphasize this point more strongly\nat lines 552-553 . Second, we have strong evidence that low-level secondary transmission did not occur\nbecause, among the few deaths of the initially uninfected larvae that occurred in our experiment before\nthe branches were removed to the lab, there were no virus infections. We now mention this point at lines\n584-586. Third, we now explicitly state at lines 552-553 in the main text that in previous work our\nexperimental design provided estimates of baculovirus transmission rates that matched estimates from bac\nulovirus epizootics in naturally occurring populations. It therefore seems unlikely that the conditions of our\nexperiments biased our estimates of transmission rates.\n10. The potential effects of host tree chemistry, such as differences in secondary metabolites or\nnutrientcontent, oninfectiondynamicsshouldbeacknowledged. Ifdifferencesinfoliagechem-\nistrybetweenDouglas-firandgrand-firinfluencelarvalfeedingbehaviourorviraltransmission,\nincorporating this information would strengthen the discussion. If direct data are unavailable,\nacknowledgingthepotentialimpactandsuggestingfutureresearchdirectionswouldbehelpful.\nAt lines 93-95 we now re-emphasize that there are likely to be strong differences in foliage chemistry\nbetween Douglas-fir and grand fir, and that these effects are likely to affect transmission. Also, as we men-\ntioned in a response to an earlier comment, in work that is in revision in The American Naturalist we show\nthat there are indeed foliage effects on infection risk, and that insect behavior plays a key role in deter-\nmining transmission. At lines 89-93, we now refer to this and other relevant citations on factors affecting\nbaculovirus transmission.\n11. The authors mention that selection mosaics have been extensively studied in plant-\npathogen systems but not in animal-pathogen systems. If there are any studies on selection\nmosaics in animals, referencing them would provide a broader context and strengthen the\nstudy’s position within the literature.\nto studies documenting selection mosaics involving animals.\n12. The model assumes that host infection risk variation is heritable. If this assumption has\nnot been explicitly validated for this system, discussing its implications and potential limita-\ntions would be useful. If relevant, the authors might also consider referencing studies that\nhave tested the heritability of similar traits in related systems.\nAt lines 465-472 in the revised text we now point out that heritability has been documented in many\nother insect-pathogen systems, including in the closely related spongy moth, Lymantria dispar.\n13. The policy implications for Douglas-fir tussock moth management are particularly inter-\nesting, especially the recommendation to use a single morphotype for biocontrol. However,\ncould using only one morphotype lead to unintended evolutionary consequences, such as in-\ncreased resistance or compensatory adaptation by the host? Addressing potential long-term\neffects would enhance the applied relevance of the study.\nThis is such an excellent suggestion that we now include a new figure in the main text showing the\neffects of using different formulations for biocontrol (which we now refer to as a “biopesticide”, following\nthe usage of the US Environmental Protection Agency). As the reviewer suspected, there are unexpected\nconsequences, inthatourmodelprojectsthatthecurrentbiopesticide, whichiscomposedonlyofthemulti-\ncapsid morphotype, will be much less effective than would biopesticides composed of both morphotypes, or\nofonlythesingle-capsidmorphotypes. Theseeffectsareindeedduetothehost’sabilitytoevolveresistance,\nwhichislikelytobemuchstrongerinthecaseofthemulti-capsidmorphotype. Wenowdiscussthesepoints\nin detail. The new figure is fig. 6, the text describing the figure is at lines 329-367, additional discussion\nof biopesticide effects is at lines 437-459, and in the SI subsection 4.9 at lines 450-491 we describe how\nwe modeled the biopesticide.\n14. In the supplementary material, Eq. 6 should be written as E . In line 230, the mean\nk,1\nshould be 1/δ , and the variance should be 1/k ∗ δ2. Equation 7 needs a bracket around\n1 1 1\n[r+r(s ν¯ +s ν¯ )]. In section 3.3 (line 323), the notations I and I should likely be lowercase.\n1 1 2 2 1 2\nPlease check the other equations for consistency, ensuring that notations are properly contin-\nued and all terms are clearly defined. The equation number is missing in line 321, and the\nfigure number is missing in line 350.\nWe appreciate the Reviewer’s attention to detail. In the SI, at line 254, equations (1)-(9), and at lines\n260-262, we have fixed these problems. We also now include the equations for both morphotypes. In the\ncaseofr+r(s ν¯ +s ν¯),whichisnowinequation(2)atline 640inthemaintext,wenowusetwodifferent\n1 1 2 2\nsizes of parentheses in an effort to achieve the same goal.\n15. The study could benefit from a discussion on how heterogeneous dispersal patterns might\ninfluence pathogen spread and morphotype coexistence. If individuals disperse non-randomly\nbased on host tree composition or pathogen load, how would this affect the selection-mosaic\nframework? Exploring whether heterogeneous dispersal could amplify or dampen selection\nmosaics would add depth to the analysis.\nAs we explained in our response to an earlier reviewer comment, there are good reasons to believe that\ndispersal is entirely random, rather than depending on infection or on host tree. In the interests of brevity\nwe therefore do not consider non-random dispersal.\n16. The influence of heterogeneous node properties in the spatial network should be consid-\nered. If certain locations in the forest act as transmission hubs due to higher host density or\nfavourable environmental conditions, how might this affect pathogen coexistence? Addressing\nwhether the spatial model incorporates such heterogeneities would provide a clearer picture\nof the system’s dynamics.\nThis is an interesting point, but there is no evidence that there are areas of higher host density or favor-\nable environmental conditions in forests infested with Douglas-fir tussock moths. Indeed, at the scale of our\nmodel forests, and at the scale of most forest patches in nature, the dispersal of the insect is high enough\nthat variation in host density is unlikely to play much of a role, as we now explain in the SI in subsection\n4.4 at lines 326-351. We now clarify these points.\nWe hope our revisions are sufficient.\nKatie Dixon and Greg Dwyer\nMosaic Theory and Host-Pathogen Theory to Explain Large-Scale\nPathogen Coexistence in an Insect-Virus System\nREVIEWERS’ COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have substantially changed their manuscript, and I feel that these changes have\nsignificantly improved the overall quality and readability of their work. While there are a few\naspects of their model I would like to see clarified, overall, I feel that this manuscript is ready\nfor publication in Nature Communications after minor revisions.\nWe appreciate the reviewer’s kind words. In what follows we address each of the reviewer’s comments.\nL271: Here are the authors referring to genetic variation (C) in the mean infectiousness as\nhost genetic variation? While infection risk is both a host and pathogen trait, by assuming\ninfection risk increases proportionally to C in the eco-evolutionary model, the authors are\ntreating this as genetic variation for the pathogen (unless I am misinterpreting the equations\nin the eco-evolutionary model).\nBecause C varies between pathogen strains, in a larger sense the reviewer is correct that C is a pathogen\ntrait. The tricky part is that C can only be larger if there is more variation between hosts, rather than be-\ntweenpathogens. ItthereforeseemsclearesttoustorefertoC asahosttraitthatdiffersbetweenpathogen\nstrains rather than as a pathogen trait.\nWe also note that we are not assuming that infection risk increases proportionally to C in the model.\nInstead, infection risk falls with C, and that effect is not because of an assumption. Instead, it is an effect\nthat arises out of the model, and is therefore a consequence of the model’s assumptions.\nL667: For the evolution of increased infectiousness, it seems that as long as heritability and\ninfectiousness variation are greater than zero, the average infectiousness will always increase,\nsince b2C2+1>1 and ρb C b C +1>1. Clearly this is not the case, as the average infectious-\n1 1 1 1 2 2\nness oscillates in Fig. 5, but I am missing how this is reflected in the equations. Can s take\non negative values? Furthermore, it seems that in this model, pathogen infectiousness can\nrespond to selection without any decrease in genetic variation.\nNo, s and s cannot take on negative values. We believe that we can dispel the reviewer’s confusion by\n1 2\nexplaining that infection risk at the end of the epizootic ν∗ is almost certainly less than infection risk at\n1,n\nthe beginning of the epizootic ν , so that ν∗ <ν . Similarly, ν∗ <ν . Indeed, it is almost certainly\n1,n 1,n 1,n 2,n 2,n\ntrue that infection risk at the end of the epizootic is vastly less than it is at the beginning of the epizootic.\nIt is therefore usually the case that after an outbreak peak we have ν <<ν and ν <<ν .\n1,n+1 1,n 2,n+1 2,n\nWe have therefore added the following text on lines 639-641:\nNotethatpathogen-inducedmortalityalwayscauseshostinfectionriskattheendofanepizootic\nto be lower, and often much lower, than it was at the beginning of the epizootic. We therefore\nhave ν∗ <ν and ν∗ <ν .\n1,n 1,n 2,n 2,n\nMinor Comments:\nL254: How are the values of s sampled? Is there any nesting for pathogen morphotype and\nhost tree species?\nwhich parameters are systematically varied one at a time. We only used a hierarchical Bayesian routine in\nthe case of our transmission data, for which the models we were fitting to the data were quite simple and\ntherefore could be used with a somewhat slow fitting algorithm like the Hamiltonian Markov-chain Monte\nCarlo routine that we used. This distinction between the two types of fitting routine is important because,\nforallthatthevaluesofsinsomeofourmodelsweremorphotypeand/ortree-speciesspecific,ourline-search\nroutine did not allow for the possibility of a Bayesian hierarchical model for the values of s. Accordingly, to\nthe extent that our line-search routine can be said to sample values of s, all the values of s were sampled\nindependently, as were the samples of all the other parameters. We have clarified these points in the text\n(lines 727-728).\nL129: It may be worth adding a sentence explaining why more resistant hosts have a greater\nimpact on overall infection than less resistant hosts.\nThis is a good point. We therefore made the following change on lines 119-122. The previous text here\nwas:\nHigher ν¯ always leads to higher host infection risk and thus higher pathogen fitness, but because\nmore resistant hosts have a stronger effect on infection rates than less-resistant hosts, higher C\nleads to lower host infection risk and thus lower pathogen fitness. Pathogen strains with higher\nvariation in host infection risk therefore have lower fitness (Fleming-Davies et al., 2015).\nThis text now says:\nHigher ν¯ always leads to higher host infection risk and thus higher pathogen fitness, but higher\nC leads to lower pathogen fitness. This is because increasing C while holding ν¯ constant is\nequivalent to adding more-susceptible hosts and more-resistant hosts in such a way that the\nmean is unaffected even though the C.V. increases; more-resistant hosts have a disproportionate\nimpactoninfectionrisk(Dwyeretal.,2000),however,andsopathogenstrainswithhighervalues\nof C have lower fitness (Fleming-Davies et al., 2015).\nNitpicks\nL45: instead of “ecological dynamics”, it might be more accurate to say a lack of considera-\ntion of environmental heterogeneity, or special dynamics, since frequency-dependent selection\nis often considered an ecological feedback.\nIn the context of this section of text we believe that the meaning of “ecological dynamics” is clear. We\ntherefore did not change this sentence.\nEquation 6: It looks like the last term in the numerator is missing s before ν∗ . It also seems\n1 1,n\nthat the text refers to the average infectiveness at the end of the epizootic as ν∗ while the\n1,n\nequations denote this as ν∗ .\n1,n\nEpically awesome catches. Both have been fixed.\nL665: Is this supposed to read b C2?\n1 1\nYes. Fixed. Again an excellent catch. We thank the reviewer for their attention to detail.\nReviewer #2 (Remarks to the Author):\nI sincerely thank the authors for their careful and thorough revision. They have addressed all\nmy concerns with great effort and clarity. The manuscript is now much improved, and I am\n2\nsatisfied with the revisions.\nWe thank the reviewer for their hard work, and we appreciate their kind words.\n3"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have substantially changed their manuscript, and I feel that these changes have significantly improved the overall quality and readability of their work. While there are a few aspects of their model I would like to see clarified, overall, I feel that this manuscript is ready for publication in Nature Communications after minor revisions.\nL271: Here are the authors referring to genetic variation (C) in the mean infectiousness as host genetic variation? While infection risk is both a host and pathogen trait, by assuming infection risk increases proportionally to C in the eco-evolutionary model, the authors are treating this as genetic variation for the pathogen (unless I am misinterpreting the equations in the eco-evolutionary model).\nL667: For the evolution of increased infectiousness, it seems that as long as heritability and infectiousness variation are greater than zero, the average infectiousness will always increase, since b_1^2C_1^2+1>1 and \\rho{b_1C}_1b_2C_2+1>1. Clearly this is not the case, as the average infectiousness oscillates in Fig. 5, but I am missing how this is reflected in the equations. Can s take on negative values? Furthermore, it seems that in this model, pathogen infectiousness can respond to selection without any decrease in genetic variation.\nMinor Comments:\nL254: How are the values of s sampled? Is there any nesting for pathogen morphotype and host tree species?\nL129: It may be worth adding a sentence explaining why more resistant hosts have a greater impact on overall infection than less resistant hosts.\nNitpicks\nL45: instead of “ecological dynamics”, it might be more accurate to say a lack of consideration of environmental heterogeneity, or special dynamics, since frequency-dependent selection is often considered an ecological feedback.\nEquation 6: It looks like the last term in the numerator is missing s_1 before v_{1,n}^\\ast. It also seems that the text refers to the average infectiveness at the end of the epizootic as {\\overline{v}}_{1,n}^\\ast while the equations denote this as\nv_{1,n}^\\ast.\nL665: Is this supposed to read b_1C_1^2?\n(Remarks on code availability)", "reviewer_2": "I sincerely thank the authors for their careful and thorough revision. They have addressed all my concerns with great effort and clarity. The manuscript is now much improved, and I am satisfied with the revisions.\n(Remarks on code availability)"}, "author_response": ""}]

Decision letter