Phenotypic plasticity explains apparent reverse evolution of fat synthesis in parasitic wasps

Experimental study and protein domain analysis

Insects

Hosts and parasitoids were maintained as previously described²⁵. Five Leptopilina heterotoma (Hymenoptera: Figitidae) populations were used for experiments: a population from Japan (Sapporo), two populations from the United Kingdom (1: Whittlesford; 2: Great Shelford) and two populations from Belgium (1: Wilsele; 2: Eupen). Information on collection sites, including GPS coordinates, can be found in²⁵.

Determination of host fat content

D. simulans and D. melanogaster hosts were allowed to lay eggs during 24 h in glass flasks containing ~ 50 mL standard medium²⁵. After two days, developing larvae were sieved and ~ 200 were larvae placed in a Drosophila tube containing ~ 10 mL medium. Seven days after egg laying, newly formed pupae were frozen at – 18 °C, after which fat content was determined as described in²⁵, where dry weight before and after neutral fat extraction was used to calculate absolute fat amount (in μg) for each host. The host pupal stage was chosen for estimating fat content, because at this point the host ceases to feed, while the parasitoid starts consuming the entire host³⁶. All data were analysed using R Project version 3.4.3⁶⁰. Fat content of hosts was compared using a one-way ANOVA with host species as fixed factor.

Manipulation of host fat content

To generate leaner D. melanogaster hosts, we adapted our standard food medium²⁵ to contain 100 times less (0.5 g) sugar per litre water. Manipulating sugar content did not alter the structure of the food medium, thus maintaining similar rearing conditions, with the exception of sugar content. Fat content of leaner and fatter D. melanogaster hosts was determined and analysed as described above.

Fat synthesis quantification of wasp populations

Mated female L. heterotoma were allowed to lay eggs on host fly larvae collected as described above with ad libitum access to honey as a food source until death. Honey consists of sugars and other carbohydrates that readily induce fat synthesis. After three weeks, adult offspring emergence was monitored daily and females were haphazardly placed in experimental treatments. Females were either killed at emergence (to measure teneral lipid reserves) or after feeding for 7 days on honey. Wasps were frozen at − 18 °C after completion of experiments. Fat content was determined as described above for hosts. The ability for fat synthesis was then determined by comparing mean fat levels of recently emerged compared to fed individuals, similar to procedures described in^15,25,28. An increase in fat levels after feeding is indicative of active fat synthesis; equal or lower fat levels suggest fat synthesis did not take place. Each population tested on D. melanogaster or D. simulans represented an independent dataset that was analysed separately, as in Visser et al. 2018²⁵, because we are interested in the response of each population on each host species. We used T-tests when data was normally distributed and variances equal, log-transformed data for non-normal data, and a Welch’s t-test when variances were unequal. We corrected for multiple testing using Benjamini and Hochberg’s False Discovery Rate⁶¹.

Fat synthesis quantification using a familial design and GC–MS analyses

To tease apart the effect of wasp genotype and host environment, we used a split-brood design where the offspring of each mother developed on lean D. simulans or fat D. melanogaster hosts in two replicated experiments (experiment 1 and 2). In both experiments, mothers were allowed to lay eggs in ~ 200 2nd to 3rd instar host larvae of one species for four days, after which ~ 200 host larvae of the other species were offered during four days. The order in which host larvae were presented was randomized across families. Following offspring emergence, daughters were allocated into two treatment groups: a control where females were fed a mixture of honey and water (1:2 w/w) or a treatment group fed a mixture of honey and deuterated water (Sigma Aldrich) (1:2 w/w; stable isotope treatment) for 7 days. Samples were prepared for GC–MS as described in ²⁸. Incorporation of up to three deuterium atoms can be detected, but percent incorporation is highest when only 1 deuterium atom is incorporated. As incorporation of a single atom unequivocally demonstrates active fat synthesis, we only analysed percent incorporation (in relation to the parent ion) for the abundance of the m + 1 ion. Percent incorporation was determined for five fatty acids, C16:1 (palmitoleic acid), C16:0 (palmitate), C18:2 (linoleic acid), C18:1 (oleic acid), and C18:0 (stearic acid), and the internal standard C17:0 (margaric acid). Average percent incorporation for C17:0 was 19.4 (i.e. baseline incorporation of naturally occurring deuterium) and all values of the internal standard remained within 3 standard deviations of the mean (i.e. 1.6). Percent incorporation of control samples was subtracted from treatment sample values to correct for background levels of deuterium (i.e. only when more deuterium is incorporated in treatment compared to controls fatty acids are actively being synthesized). For statistical analyses, percent incorporation was first summed for C16:1, C16:0, C18:2, C18:1 and C18:0 to obtain overall incorporation levels, as saturated C16 and C18 fatty acids are direct products of the fatty acid synthesis pathway (that can subsequently be desaturated).

Data (presented in Fig. 1) was analysed by means of a linear mixed effects model (GLMM, lme4 package) with host (lean D. simulans and fat D. melanogaster) and experiment (conducted twice) as fixed effect, family nested within population (Japan, United Kingdom 1 and 2, Belgium 1 and 2) as random factor, and percentage of incorporation of stable isotopes as dependent variable (log transformed; n = 138). Non-significant terms (i.e., experiment) were sequentially removed from the model to obtain the minimal adequate model as reported in Table 2. When referring to “families,” we are referring to the comparison of daughters of singly inseminated females, which (in these haplodiploid insects) share 75% of their genome.

Identification of functional acc and fas genes in distinct parasitoid species

To obtain acc and fas nucleotide sequences for L. clavipes, G. legneri, P. maculata and A. bilineata, we used D. melanogaster mRNA ACC transcript variant A (NM_136498.3 in Genbank) and FASN1-RA (FBtr0077659 in FlyBase) and blasted both sequences against transcripts of each parasitoid (using the blast function available at http://www.parasitoids.labs.vu.nl^62,63). Each nucleotide sequence was then entered in the NCBI Conserved Domain database⁶⁴ to determine the presence of all functional protein domains. All sequences were then translated using the Expasy translate tool (https://web.expasy.org/translate/), where the largest open reading frame was selected for further use and confirming no stop codons were present. Protein sequences were then aligned using MAFFT v. 7 to compare functional amino acid sequences between all species (Supplementary files 1 and 2)⁶⁵.

Simulation study

We consider the general situation where phenotypic plasticity is only sporadically adaptive and ask the question whether and under what circumstances plasticity can remain functional over long evolutionary time periods when the regulatory processes underlying plasticity are gradually broken down by mutations. We consider a regulatory mechanism that switches on or off a pathway (like fat synthesis) in response to environmental conditions (e.g., host fat content).

Fitness considerations

We assume that the local environment of an individual is characterized by two factors: fat content F and nutrient content N, where nutrients represent sugars and other carbohydrates that can be used to synthesize fat. Nutrients are measured in units corresponding to the amount of fat that can be synthesized from them. We assume that fitness (viability and/or fecundity) is directly proportional to the amount of fat stored by the individual. When fat synthesis is switched off, this amount is equal to F, the amount of fat in the environment. When fat synthesis is switched on, the amount of fat stored is assumed to be (N – c + (1 – k)F). This expression reflects the following assumptions: (i) fat is synthesized from the available nutrients, but this comes at a fitness cost c; (ii) fat can still be absorbed from the environment, but at a reduced rate ((1 – k)). It is adaptive to switch on fat synthesis if (N – c + (1 – k)F) is larger than F, or equivalently if (F < tfrac{1}{k}(N – c)).

The right-hand side of this inequality is a straight line, which is illustrated by the blue line in Fig. 4. The three boxes in Fig. 4 illustrate three types of environmental conditions.

• Red box low-fat environments. Here, (F < tfrac{1}{k}(N – c)) is always satisfied, implying that fat synthesis should be switched on constitutively.

• Yellow box high-fat environments. Here, (F > tfrac{1}{k}(N – c)), implying that fat synthesis should be switched off constitutively.

• Orange box intermediate-fat environments. Here, fat synthesis should be plastic and switched on if for the given environment (N, F) the fat content is below the blue line and switched off otherwise.

The simulations reported here were all run for the parameters (k = tfrac{1}{2}{text{ and }}c = tfrac{1}{4}). We also investigated many other combinations of these parameters; in all cases, the results were very similar to those reported in Fig. 3.

Gene regulatory networks (GRN)

In our model, the switching device was implemented by an evolving gene regulatory network (as in van Gestel and Weissing⁶⁶). The simulations shown in Fig. 3 of the main text are based on the simplest possible network that consists of two receptor nodes (sensing the fat and the nutrient content in the local environment, respectively) and an effector node that switches on fat synthesis if the combined weighted input of the two receptor nodes exceeds a threshold value T and switches it off otherwise. Hence, fat synthesis is switched on if (w_{F} F + w_{N} N > T) (and off otherwise). The GRN is characterized by the weighing factors (w_{F} {text{ and }}w_{N}) and the threshold T. These parameters are transmitted from parents to offspring, and they evolve subject to mutation and selection. We also considered alternative network structures (all with two receptor nodes and one effector node, but with a larger number of evolvable weighing factors⁶⁷, and obtained very similar results, see below).

For the simple GRN described above, the switching device is 100% adaptive when the switch is on (i.e., (w_{F} F + w_{N} N > T)) if (F < tfrac{1}{k}(N – c)) and off otherwise. A simple calculation yields that this is the case if: (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N}).

Evolution of the GRN

For simplicity, we consider an asexual haploid population with discrete, non-overlapping generations and fixed population size (N = 10,000). Each individual has several gene loci, each locus encoding one parameter of the GRN. In case of the simple network described above, there are three gene loci, each with infinitely many alleles. Each individual harbours three alleles, which correspond to the GRN parameters (w_{F} {, }w_{N} {text{ and }}T), and hence determine the functioning of the genetic switch. In the simulations, each individual encounters a randomly chosen environment ((N{, }F)). Based on its (genetically encoded) GRN, the individual decides on whether to switch on or off fat synthesis. If synthesis is switched on, the individual’s fitness is given by (N – c + (1 – k)F); otherwise its fitness is given by F. Subsequently, the individuals produce offspring, where the number of offspring produced is proportional to the amount of fat stored by an individual. Each offspring inherits the genetic parameters of its parent, subject to mutation. With probability μ (per locus) a mutation occurs. In such a case the parental value (in case of a simple network: the parent’s allelic value (w_{F} {, }w_{N} {text{ or }}T)) is changed to a mutated value ((w_{F} { + }delta {, }w_{N} { + }delta {text{ or }}T + delta)), where the mutational step size δ is drawn from a normal distribution with mean zero and standard deviation σ. In the reported simulations, we chose (mu = 0.001) and (sigma = 0.1). The speed of evolution is proportional to (mu cdot sigma^{2}), implying that the rate of change in Fig. 3 (both the decay of plasticity and the rate of regaining adaptive plasticity) are positively related to μ and σ.

Preadaptation of the GRNs

Starting with a population with randomly initialized alleles for the GRN parameters, we first let the population evolve for 10,000 generations in the intermediate-fat environment (the orange box in Fig. 4). In all replicate simulations, a “perfectly adapted switch” (corresponding to (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N})) evolved, typically within 1,000 generations. Still, the evolved GRNs differed across replicates, as they evolved different values of (w_{N} > 0). These evolved networks were used to seed the populations in the subsequent “decay” simulations.

Evolutionary decay of the GRNs

For the decay experiments reported in Fig. 3 of the main text, we initiated a large number of monomorphic replicate populations with one of the perfectly adapted GRNs from the preadaptation phase. These populations were exposed for an extended period of time (1,000,000 generations) to a high-fat environment (the yellow box in Fig. 4), where all preadapted GRNs switched off fat synthesis. However, in some scenarios, the environmental conditions changed back sporadically (with probability q) to the intermediate-fat environment (the orange box in Fig. 4), where it is adaptive to switch on fat metabolism in 50% of the environmental conditions (when (N, F) is below the blue line in Fig. 4). In Fig. 3, we report on the changing rates (q = 0.0) (no changing back; red), (q = 0.001) (changing back once every 1,000 generations; purple), and (q = 0.01) (changing back once every 100 generations; pink). When such a change occurred, the population was exposed to the intermediate-fat environment for t generations (Fig. 3 is based on t = 3).

Throughout the simulation, the performance of the network was monitored every 100 generations as follows: 100 GRNs were chosen at random from the population, and each of these GRNs was exposed to 100 randomly chosen environmental conditions from the intermediate-fat environment (orange box in Fig. 4). From this, we could determine the average percentage of “correct” decisions (where the network should be switched on if and only if (F < tfrac{1}{k}(N – c)). 1.0 means that the GRN is still making 100% adaptive decisions; 0.5 means that the GRN only makes 50% adaptive decision, as would be expected by a random GRN or a GRN that switches the pathway constitutively on or off. This measure for performance in the “old” intermediate-fat environment was determined for 100 replicate simulations per scenario and plotted in Fig. 3 (mean ± standard deviation).

Evolving robustness of the GRNs

The simulations in Fig. 3 are representative for all networks and parameters considered. Whenever (q = 0.0), the performance of the regulatory switch eroded in evolutionary time, but typically at a much lower rate in case of the more complex GRNs. Whenever (q = 0.01), the performance of the switch went back to levels above 90% and even above 95% for the more complex GRNs. Even for (q = 0.001), a sustained performance level above 75% was obtained in all cases.

Intriguingly, in the last two scenarios the performance level first drops rapidly (from 1.0 to a much lower level, although this drop is less pronounced in the more complex GRNs) and subsequently recovers to reach high levels again. Apparently, the GRNs have evolved a higher level of robustness, a property that seems to be typical for evolving networks⁸. For the simple GRN studied in Fig. 3, this outcome can be explained as follows. The initial network was characterized by the genetic parameters (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N}) (see above), where (w_{N}) was typically a small positive number. In the course of evolutionary time, the relation between the three evolving parameters remained approximately the same, but (w_{N}) (and with it the other parameters) evolved to much larger values. This automatically resulted in an increasingly robust network, since mutations with a given step size distribution affect the performance of a network much less when the corresponding parameter is large in absolute value.

Costs of plasticity

Phenotypically plastic organisms can incur different types of costs⁶⁸. In our simple model, we only consider the cost of phenotype-environment mismatching, that is, the costs of expressing the ‘wrong’ phenotype in a given environment. When placed in a high-fat environment, the preadapted GRNs in our simulations take the ‘right’ decision to switch off fat metabolism. Accordingly, they do not face any costs of mismatching. Yet, the genetic switch rapidly decays (as indicated in Fig. 3 by the rapid drop in performance when tested in an intermediate-fat environment), due to the accumulation of mutations.

It is not unlikely that there are additional fitness costs of plasticity, such as the costs for the production and maintenance of the machinery underlying plasticity⁶⁸. In the presence of such constitutive costs, plasticity will be selected against when organisms are living in an environment where only one phenotype is optimal (as in the high- and low-fat environments in Fig. 4). This would obviously affect the evolutionary dynamics in Fig. 3, but the size of the effect is difficult to judge, as the constitutive costs of plasticity are notoriously difficult to quantify. In case of the simple switching device considered in our model, we consider the constitutive costs of plasticity as marginal, but these costs might be substantial in other scenarios.

Source: Ecology - nature.com