Investigating the consequences of the mating system for drug resistance evolution in Caenorhabditis elegans

TL;DR

This study explores how the mating system of a nematode affects the evolution of drug resistance, showing that androdioecious populations adapt quickly and efficiently.

Contribution

The study introduces a polygenic population genetic model combined with pharmacodynamic approaches to analyze drug resistance evolution in androdioecious populations.

Findings

Androdioecious populations show rapid initial adaptation and high drug concentration tolerance.

These populations exhibit the highest diversity and fastest fixation of beneficial alleles.

The findings suggest androdioecious systems optimize reproductive strategies under drug selection.

Abstract

The rise of anthelmintic-resistant strains in livestock threatens both animal and human health. Understanding the factors influencing anthelmintic resistance is crucial to mitigate the threat posed by these parasites. Due to difficulties in studying parasitic worms in the laboratory, the non-parasitic nematode Caenorhabditis elegans is used as a model organism to investigate anthelmintic resistance evolution. However, the suitability of this free-living nematode as a model for parasitic worms is debatable due to its rare androdioecious reproductive system, raising questions about the generalizability of findings from evolutionary experiments in C. elegans to other species. In this study, we developed a polygenic, population genetic model combined with pharmacodynamic approaches to investigate the effects of reproductive strategy and other aspects, such as dominance, mutational effects, the number of loci and population size, on determining the dynamics and outcome of evolutionary processes. We found that androdioecious populations showed both rapid initial adaptation typical for hermaphrodites and tolerance to high drug concentrations observed in dioecious populations. They also exhibited the highest diversity and shortest time for the fixation of the beneficial allele. These results suggest that androdioecious populations can harness the advantages of both selfing and outcrossing, optimizing their reproductive strategy in response to drug selection.

Figures

Modelling evolution of resistance in diploid organisms . (A) Diploid genomes with four loci, each with two alleles (0 or 1), whereby 1 represents a mutation. (B) A serial passage experiment with increasing drug concentration. (C) Drug resistance for wild-type and mutant strains, indicated through the relationship between drug concentration (X axis) and fitness (Y axis). (D) Different reproduction alternatives between a female/hermaphrodite (in purple) and a male (blue). While selfing can generate offspring of various genotypes, all are hermaphrodites. Mating produces various genotypes equally split between males and hermaphrodites. (E) Fitness function of different strains. In this illustrated example, resistance is encoded by two independent recessive mutations. Sensitive strain (cyan) has the highest fitness in the absence of the drug, but the lowest EC50 (shown as a dashed line). Two different single mutants (with homozygous mutations at one of the two loci, shown in blue and purple) have lower fitness but higher EC50. The double mutant (homozygous for mutations at both loci, shown in purple) has the lowest fitness of all the strains in the absence of the drug but the highest EC50. Mutations are recessive; thus, heterozygous loci do not contribute to resistance, nor do they reduce fitness.

Adaptation to increasing drug concentration, assuming recessive mutations and six loci. Mean and s.d. of 100 simulation trials. In androdioecious populations, 50% of hermaphrodites reproduce by mating and 50% by selfing. (A–C) First 50 generations. (D–F) All 500 generations.

Number of homozygous loci over time, assuming recessive mutations and six loci. Mean and s.d. of 100 simulation trials. In androdioecious populations, 50% of hermaphrodites reproduce by mating and 50% by selfing. (A–C) First 50 generations. (D–F) All 500 generations.

Time to fixation of the first beneficial allele for various fractions of mating hermaphrodites ξ . Resistance is encoded by six loci. Note the various scales on the Y axis. Box plots show the median, interquartile range and whiskers extending to 1.5× the interquartile range from the box edges.

Biodiversity changes over time. Recessive and dominant mutations, six loci encoding resistance. (A–C) Recessive mutations, and (D–F) dominant mutations. In androdioecious populations, 50% of hermaphrodites reproduce by mating and 50% by selfing.

Drug concentration reached by populations at different time points, recessive mutations. Blue represents hermaphrodite populations, red represents dioecious ones and green represents various levels of androdioecy. Box plots show the median, interquartile range and whiskers extending to 1.5× the interquartile range from the box edges. Red line connects the means.

Equations

Full text

Introduction

Helminths, particularly parasitic nematodes, are diverse parasites of plants, animals and humans. Despite advances in sanitation, helminth infections are among the most common infections worldwide, affecting large human populations, mainly in low-income countries and in the poorest and most deprived communities. They cause a significant burden on the population, drastically decreasing the quality of life and life expectancy, especially when combined with malnutrition.

Anthelmintic drugs are generally very efficient, clearing infections within days. However, extensive use of these drugs in grazing animals led to the evolution of drug-resistant strains [1–4]. Their prevalence has increased globally in recent years, representing a severe problem in various livestock sectors [5,6] becoming an increasingly severe threat to agriculture and human health [1,7–9]. Indeed, the reduced sensitivity of Ascaris populations to benzimidazole has already been detected in children in Rwanda [ 10]. As many parasitic worms are zoonotic, anthelmintic resistance in worms is a significant focus of One Health initiatives.

Understanding the factors that influence the evolution of anthelmintic resistance in parasitic worms is essential to mitigate their threat. Unfortunately, conducting laboratory experiments involving parasitic worms poses challenges due to the requirement of an experimentally accessible host organism, the lack of genetic tools for analysis of the parasites and the long life cycle of many parasites. Consequently, these studies are relatively rare. Therefore, the intensively studied nematode Caenorhabditis elegans serves as a model organism for parasitic nematode infections despite not being a parasite, mainly due to its short life span and easy maintenance in the laboratory setting, well-annotated genome and similar response to anthelmintics [11,12]. It has been used successfully to gain important insights into the evolution and genetics of resistance against anthelminthic drugs [13–18].

However, the suitability of the free-living C. elegans as a model organism for parasitic worms has also been challenged [12]. Beyond its free-living lifestyle, one of the key distinctions lies in its reproductive mode. Caenorhabditis elegans is androdioecious, with populations consisting of self-fertilizing hermaphrodites and a small proportion of males. While nematode reproductive modes are remarkably diverse [19], most of the parasitic worms reproduce by mating between males and females (dioecy). Overall, androdioecy is relatively rare, with only a few examples [20,21].

This diversity in reproductive strategies raises crucial questions for understanding adaptation. Theoretical studies in population genetics have emphasized the potential impact of reproductive mode on the rate and trajectory of evolution [22–30]. However, these studies often focus on the broad divide between sexual and asexual reproduction [31], rather than the finer distinctions among mating systems like dioecy, hermaphroditism and the complex dynamics of androdioecy. As a rare exception, [32] developed a one-locus model to explore the dynamics of androdioecious species, finding that a slight reproductive advantage of males can maintain high male frequencies. Weeks et al. [33] conducted a comparative analysis of androdioecy in animals, suggesting that it may have evolved as a response to reproductive assurance in species with episodic low densities.

Due to the limited number of studies on androdioecy in animals (see [34] for exceptions), our understanding of the androdioecious systems and their consequences on evolutionary adaptation comes primarily from plant research, where androdioecy is more common. Some researchers suggest that androdioecy is an intermediate state between hermaphroditism and dioecy (see [21] for a detailed review). It has been shown that selfing by hermaphrodites provides reproductive assurance, allowing populations to persist through periods of low density when mates are scarce [35]. In addition, selfing increases homozygosity, exposing new mutations (both beneficial and deleterious) to natural selection, thereby increasing its efficacy [36]. However, selfing may reduce genetic variation over time, limiting the potential to adapt to novel conditions. In contrast, when hermaphrodites outcross with males, new genetic variation can be introduced into the population, which may be especially suitable for adaptation to changing environments [34,36]. Furthermore, outcrossing is expected to compensate for lower fitness caused by inbreeding depression [37].

The evolution of androdioecy in C. elegans has already been studied by experimental and modelling approaches, taking into consideration mating ability and spontaneous production of males [38]. Chasnov & Chow [39] suggest that dioecy can transition to androdioecy when the benefits of reproductive assurance outweigh the costs of selfing. Using experimental evolution, Anderson et al. [40] show that males can be maintained at high frequencies when populations have genetic diversity or are exposed to novel environments. Teotonio et al. [41] suggest that in androdioecious populations with standing genetic diversity, outcrossing rates were maintained at an optimal level of around 50% through stabilizing selection on male frequencies. On the other hand, males appear to reproduce at a rate just below that necessary for them to be maintained [42]. However, the previous models usually assumed incomplete mating abilities of the males, which is particularly true for the lab-adapted C. elegans strain N2, but not for most natural isolates of C. elegans [43,44]. Therefore, Cutter et al. [45] state that males in Caenorhabditis are understudied compared to hermaphrodites while playing an important role in sexual selection, sexual conflict and reproductive mode transitions. In summary, the androdioecious reproductive mode presents a complex landscape of potential evolutionary advantages and disadvantages. However, due to the limited number of studies, the evolutionary consequences of androdioecy remain poorly understood, and the generalizability of findings from C. elegans to other nematode species with diverse reproductive modes remains an open question.

To address this gap, this study investigates the suitability of C. elegans as a model organism for studying the evolution of antihelmintic resistance in parasitic nematodes. We aim to improve our understanding of the evolutionary consequences of androdioecy using a mathematical modelling approach. For this, we developed a polygenic population genetic model of a nematode population with dioecious (mating), hermaphroditic (selfing) or androdioecious (combined selfing and mating) reproduction. We investigate how these reproductive modes influence the rate and outcome of adaptation to an anthelmintic-containing environment. Moreover, we investigate how other factors, such as dominance, population size, the number of loci and the size of their effects, influence the consequences of the mating system. To incorporate the selective pressures of this environment, we utilize a modified pharmacodynamic framework from Regoes et al. [46], which captures the varying effects of drug concentrations on populations with different mutations, thus enabling us to model anthelmintic resistance in diploid organisms. This approach, combined with the population genetics [47], promises novel insights into the role of reproductive strategy in drug resistance evolution among parasitic worms.

The model

Below, we describe our compartmental, polygenic, population genetic model. Inputs of the model are various parameters determining the genetics (e.g. number of loci), the translation to phenotype (e.g. fitness effects of individual mutations) and the details of the experiment (e.g. length of the treatment, drug concentration; see electronic supplementary material, table S1 for details). The model compartments represent a particular class of individuals defined by their sex and genotype. We assume that the organisms are diploid.

How each genotype manifests depends on whether the mutations are dominant or recessive. Each genotype is associated with a basal fitness (in the absence of the drug), and it has an $[eqn]$ value (concentration of a substance that produces 50% of its maximum possible effect), as explained in detail below. Hermaphrodites from each compartment can self-fertilize or mate, laying a number of eggs according to the origin of sperm, their fitness and genotype. These eggs may have various genotypes, assuming that genetic loci are not linked and are carried on different chromosomes. Additionally, mutations may occur. Finally, offspring generation replaces the parental one and may be reduced to a required size depending on the modelled scenario. All of these processes are explained in detail below. See figure 1 for a schematic representation of the overall simulation process. See electronic supplementary material for details of the simulation process.

Modelling evolution of resistance in diploid organisms . (A) Diploid genomes with four loci, each with two alleles (0 or 1), whereby 1 represents a mutation. (B) A serial passage experiment with increasing drug concentration. (C) Drug resistance for wild-type and mutant strains, indicated through the relationship between drug concentration (X axis) and fitness (Y axis). (D) Different reproduction alternatives between a female/hermaphrodite (in purple) and a male (blue). While selfing can generate offspring of various genotypes, all are hermaphrodites. Mating produces various genotypes equally split between males and hermaphrodites. (E) Fitness function of different strains. In this illustrated example, resistance is encoded by two independent recessive mutations. Sensitive strain (cyan) has the highest fitness in the absence of the drug, but the lowest EC50 (shown as a dashed line). Two different single mutants (with homozygous mutations at one of the two loci, shown in blue and purple) have lower fitness but higher EC50. The double mutant (homozygous for mutations at both loci, shown in purple) has the lowest fitness of all the strains in the absence of the drug but the highest EC50. Mutations are recessive; thus, heterozygous loci do not contribute to resistance, nor do they reduce fitness.

The model’s output is a complete record of the population sizes of all compartments over time. Therefore, it is possible to observe population dynamics (e.g. changes in population sizes, allele frequencies or accumulation of the mutations), the timing of events (e.g. when a mutant appears or reaches a certain frequency) and to calculate the emerging population properties (e.g. diversity) and observe their changes over time. The model is implemented in Python. The simulations are stochastic, allowing us to determine the probability distributions of the outcomes.

Genotypes: Parasitic worms, as well as C. elegans, are diploid organisms, in which the drug resistance is known to be encoded by multiple loci [13,17,48–50]. Therefore, we represent their genotype, consisting of $[eqn]$ unlinked loci relevant for the modelled adaptation process, as a string of length $[eqn]$ . Each position describes the number of mutated alleles at that particular locus: 0 represents a sensitive homozygote (at the particular locus), 1 represents a heterozygote and 2 represents a homozygote with two mutated alleles. With $[eqn]$ loci, there are $[eqn]$ possible genotypes. These represent sub-populations (classes) with possibly distinct pharmacodynamic properties.

As we assume a diploid genome, the relationship between the alleles of the same locus must also be defined. Vector $[eqn]$ (length $[eqn]$ ) determines the dominance of mutation at individual loci. If $[eqn]$ , the mutated allele at locus $[eqn]$ is dominant, and genotypes that have one or two mutated alleles at this locus manifest in the same way. If $[eqn]$ , the mutated allele at locus $[eqn]$ is recessive, and only homozygous mutations lead to a fitness effect. Using $[eqn]$ , for each genotype $[eqn]$ , we determine phenotype $[eqn]$ (length $[eqn]$ ), with its elements capturing whether the effect of each mutated locus should be considered or not.

Sexes: From each genotype, it is possible to have individuals who are capable of contributing their genetic information to the next generation but do not lay eggs, further referred to as males. In addition, we also include individuals who contribute their genetic information and are capable of laying eggs. Depending on the simulated scenario, these individuals may be able to self (hermaphrodites) or require mating in order to lay eggs (representing females). However, due to simplicity, we will refer to all egg-laying individuals as hermaphrodites throughout this manuscript, and the ability to lay eggs with or without male presence will be modulated by a parameter $[eqn]$ , as explained below. Individuals of both sexes may have any possible genotype, altogether creating $[eqn]$ of sub-populations (compartments) that are simulated.

Pharmacodynamics and fitness: The effects of anthelmintic drugs on worms can vary at different stages of their life cycle. These drugs can kill adult worms, paralyze them, prevent them from reproducing or block the hatching of their eggs. In our model, we simplified the impact of anthelmintics on worm fitness by focusing solely on a reduced number of hatched eggs and thus the number of offspring a particular genotype can produce, regardless of the underlying causes. Moreover, we define fitness as the number of offspring that survive until adulthood and are able to reproduce. In this sense, fitness is a result of various processes—from decreased survival of the genotype in focus, with subsequently reduced offspring production and decreased mobility, thus affecting the ability to feed/find a mating partner(s), or just having fewer eggs per individual.

In the presence of an anthelmintic, the fitness of all individuals is reduced. The concentration of the anthelmintic required to produce a specific effect (such as killing, paralyzing the worms or preventing reproduction) is higher for the mutated worms, indicating an increase in the $[eqn]$ compared to the sensitive strain. However, we also assume that mutations conferring drug resistance come with a cost in the absence of the drug, leading to a reduction in fitness.

Therefore, in our model, each resistant allele is associated with a benefit increasing the carrier’s $[eqn]$ in the presence of drugs and a cost-reducing fitness in its absence. Costs and benefits are given as vectors c and b of length $[eqn]$ , respectively, with each element corresponding to the respective locus. Mutational effects of resistance-conferring mutations can be of various sizes and can interact together [51]. To avoid introducing epistatic effects, we assume that both costs and benefits combine in an additive manner [47,52–55].

The total cost $[eqn]$ associated with phenotype $[eqn]$ is $[eqn]$ And the total benefit $[eqn]$ associated with phenotype $[eqn]$ is $[eqn]$

The total cost associated with a particular genotype (strain) determines its fitness $[eqn]$ in the absence of the drug $[eqn]$ , where the fitness of the fully sensitive strain is 1. The total benefit captures the increase of the $[eqn]$ of this strain as $[eqn]$ , where $[eqn]$ is the $[eqn]$ of the sensitive wild-type. The fitness in the absence of the drug and $[eqn]$ define pharmacodynamic curves that allow us to determine the expected fitness of any strain $[eqn]$ (figure 1)

$$[eqn]$$

where $[eqn]$ is the drug concentration and $[eqn]$ defines the steepness of the curve.

Egg count: Apart from fitness, the number of eggs the hermaphrodite can lay is also influenced by the origin of sperm: whether it comes from mating or selfing. In line with the experimental observations in C. elegans, we will assume that mated, wild-type individuals can lay, on average, 1000 eggs in a drug-free environment when out crossed ( $[eqn]$ ) but only 300 when selfing ( $[eqn]$ ) [56]. Therefore, the number of eggs that the individual with genotype $[eqn]$ lays is $[eqn]$ if the individual reproduces by selfing, or $[eqn]$ if outcrossed.

Reproduction: The fraction of males in the population is denoted $[eqn]$ , the rest are hermaphrodites. Depending on the simulated scenario, a fraction $[eqn]$ of hermaphrodites in the population engage in mating with males, while the remainder reproduces through self-fertilization. Mating between hermaphrodites and males is random with respect to male genotype and male fraction $[eqn]$ , provided that there are at least some males. Thus, the frequency of mating between a hermaphrodite and a male of a given genotype is proportional to the frequency of that male genotype in the population. Offspring are generated with probabilities calculated using Punnet squares. Half of the offspring (see the definition of egg count) in each genotype class are males, while the other half are hermaphrodites, merged with the populations calculated in the selfing step outlined below.

Due to recombination and/or chromosomal shuffling, selfing individuals can also produce offspring with different genotypes, but these are always hermaphrodites. As the formation of males during selfing is rare, in our model, we assume that males are only generated by mating. A probability of a hermaphrodite that belongs to the genotype $[eqn]$ to generate eggs/offspring that belongs to the sub-population with genotype $[eqn]$ is calculated based on the Mendelian rules using Punnet squares for $[eqn]$ loci. The linkage between loci is not currently considered; any allele can be inherited independently from another allele at a different locus.

Mutations: During reproduction, mutations may occur. We assume that one mutation in the genome occurs directly after fertilization (either by mating or selfing) with probability $[eqn]$ , changing one locus in the whole genome from 0 to 1, 1 to 0 or 2 or from 2 to 1. Only a single mutation is allowed per individual, so approximately (sampling from Poisson distribution) $[eqn]$ individuals from genotype class $[eqn]$ can mutate into any genotype class in the mutational distance 1. If there are multiple target genotype classes, $[eqn]$ individuals are divided into all of them in equal proportions.

Environment: The experiment, or the treatment, consists of any number of cycles (length of the cycle corresponds to a generation time), starting drug concentration, changes in drug concentrations and their schedule, initial population sizes, the drug degradation rate, conditions for the environmental change (e.g. increase of concentration) and the dilution of the final population before it enters the next cycle. By setting these parameters, many different treatment regimens can be modelled: from a single treatment through serial passages commonly used in evolutionary experiments to a chemostat environment.

Simulated scenarios

We take advantage of this study’s computational nature and investigate a wide range of parameters and aspects that would not be feasible in a laboratory setting (see electronic supplementary material, table S3 for the summary of simulated scenarios). For each scenario, we conducted 100 repeats. Moreover, we followed the populations for 500 generations, which would correspond to more than four years of continuous evolution experiments with C. elegans. Below, we explain the different aspects we focused on and the parameters that we used.

Mating systems

(a)

In our simulations, the fraction of hermaphrodites that mate ( $[eqn]$ ) is independent of the fraction of males in the population. This corresponds to a biological scenario where the availability of males is not a limiting factor, but rather, the mating fraction is determined by the avoidance or inclination of hermaphrodites to mate. This means that if any males are present in the population, the fraction of hermaphrodites $[eqn]$ will mate, with the frequency of matings contributed by males of different genotypes proportional to the frequency of those genotypes in the male population. By varying the parameter $[eqn]$ , we can simulate different mating systems, from hermaphroditic through androdioecious to dioecious.

As males are only generated during mating, their fraction is directly influenced by the fraction of mating hermaphrodites

The number of male offspring $[eqn]$ and $[eqn]$ is the number of hermaphrodites in the parental generation. The number of hermaphrodite offspring is $[eqn]$ , where the first term is generated by mating and the second term by selfing. Therefore, the fraction of males is

$$[eqn]$$

and it is independent of $[eqn]$ , as well as the initial fraction of males (provided there were some males). It only depends on the mating fraction of hermaphrodites $[eqn]$ , and the fecundity of hermaphrodites when selfing and mating. For simplicity, we initialize all simulations with a fraction of males of 0.5.

To simulate the hermaphroditic mating system, we set the fraction of hermaphrodites that mate to $[eqn]$ . This means that males are purged from the population after a single generation, and the population remains hermaphroditic.

On the other hand, to simulate the dioecious mating system, we set the fraction of mating hermaphrodites to $[eqn]$ . This means that all the egg-laying individuals reproduced only by mating, essentially representing females. Due to the stochasticity of the system, the fraction of males fluctuates around $[eqn]$ .

When $[eqn]$ , some hermaphrodites reproduce by mating, while others by self-fertilization, representing an androdioecious mating system. We simulated a range of values of $[eqn]$ in increments of 0.1, where $[eqn]$ . The fraction of males fluctuates around $[eqn]$ given by equation (3.1).

Genetics, population size and time scale

(b)

First, we assume two loci with relatively small effects (cost 5 %, benefit 0.5 nM), resulting in nine possible genotypes. Then, to investigate evolutionary dynamics in more complex scenarios, we assume six independent loci, three with small effects (cost 5 %, benefit 0.5 nM) and three with large effects (cost 10 %, benefit 3 nM). It allows the generation of 729 different genotypes, reflecting the complexity found in nature. Note that the fitness (in the absence of a drug) of the sensitive strain (with no mutation) is 1, and its $[eqn]$ is 1 nM. The mutation rate is $[eqn]$ mutations per genome. We investigate adaptation when all loci are dominant or all loci are recessive. Size of the population is known to play an important role in the rate of adaptation. Therefore, we performed simulations with 200, 2000, 20 000 and 200 000 individuals. When evaluating the evolutionary consequences of the mating strategy, it is important to consider time scales. While some mating strategies may encourage short-term survival and rapid adaptation, others enable better survival in the long run. Therefore, we performed 100 simulations of 500 generations for each parameter set. Having complete information about the dynamics (population sizes of all compartments at all times), we compared observations made at different time points and discussed conclusions that would be made if only the particular time point was investigated. In particular, we looked at populations after 10, 25, 100, 250 and 500 generations.

Environment

(c)

Inspired by in vitro evolution experiments performed in C. elegans to investigate drug resistance evolution, we simulated a gradually increasing drug concentration (e.g. ivermectin), with concentration increases following a pre-defined sequence shown in electronic supplementary material, table S2. The population starts in an environment with no drug concentration, which increases over time when the population is deemed sufficiently adapted (having an average egg count of at least 20 per individual). Populations are passed from one generation to another, increasing the drug concentration but allowing constant growth (no limits imposed by carrying capacity). This approach has been shown to enable the evolution of high resistance while limiting the probability of extinction [15–17], which may occur if the concentration is increased too steeply. Electronic supplementary material, figure S1A, shows preliminary simulations of populations exposed to 10 nM drug concentration without prior exposure to lower drug concentrations. The majority of these populations go extinct within the first 20 generations, especially in smaller populations. Moreover, gradual exposure to drugs leads to the survival of all populations regardless of their size (electronic supplementary material, figure S1B) and their adaptation to 10 nM concentration or even higher (electronic supplementary material, figure S1C).

Results

In this study, we simulated the evolution of worm populations in an environment with increasing concentrations of ivermectin as a model to study nematode resistance evolution. The concentration was raised when the population was deemed sufficiently adapted to its current conditions. Thus, the rate of increase in drug concentration acts as a proxy for the rate of adaptation. We examined how various genetic and environmental factors, along with different mating strategies, influenced the rate of adaptation and the maximum concentration that the populations reached over a specified period.

Mating facilitates sustained adaptation in large populations with recessive mutations

(a)

Our findings demonstrate that a gradual drug concentration increase enables populations to adapt to relatively high drug concentrations (figure 2). Population size is a critical factor, with larger populations exhibiting, in general, faster adaptation. As simulations are initiated with wild-type populations (no mutations are present), it suggests that mutational supply limits adaptation in smaller populations, and mutations are less frequent in smaller populations.

Adaptation to increasing drug concentration, assuming recessive mutations and six loci. Mean and s.d. of 100 simulation trials. In androdioecious populations, 50% of hermaphrodites reproduce by mating and 50% by selfing. (A–C) First 50 generations. (D–F) All 500 generations.

However, the interplay between population size and adaptation is more complex, as is particularly evident in figure 2C, for mating populations. Here, larger populations (200 000 individuals) initially lag behind smaller ones (2000 and 20 000 individuals) in progressing to higher concentrations. A similar effect is also seen in androdioecious populations ( $[eqn]$ ), where populations of 2000 individuals are the first to overcome the 24 nM threshold (figure 2E). This indicates that while the mutation supply limits small populations, once the mutation appears, it takes longer for the mutant to outcompete the others in larger populations.

Furthermore, figure 2 highlights distinct adaptation patterns across the three mating strategies. Strictly selfing populations (hermaphrodites only) rapidly progress to 24 nM but fail to exceed this threshold within the simulation timeframe (figure 2D). In contrast, strictly mating populations exhibit slower initial adaptation but then continue to progress to 28 nM (figure 2F). Androdioecious populations display a rapid initial adaptation, similar to that of hermaphrodites, followed by a slower progression to higher concentrations. Nonetheless, they ultimately adapt to a greater extent than hermaphrodites (figure 2E). Moreover, adaptation of the smallest androdioecious populations in the given timeframe exceeded that of the smallest dioecious populations. This suggests that androdioecious populations may be able to reap the benefits of both mating and selfing strategies. Therefore, evolution experiments with C. elegans might overestimate the adaptation to anthelmintics, especially at low population sizes.

Rate of adaptation to higher concentration does not directly reflect genetic changes

(b)

The increasing ability to survive in higher concentrations is enabled by the gradual accumulation of resistance mutations. However, the rate of adaptation does not necessarily reflect the rate of mutation accumulation, as other differences between reproductive strategies are also present.

Figure 3 illustrates the average number of homozygous loci within the population over time, considering that resistance mutations are recessive, with six loci potentially contributing to resistance. Hermaphroditic populations exhibit a faster accumulation of homozygous mutations compared to strictly mating populations, which demonstrate the slowest adaptation. Furthermore, hermaphroditic populations display a stepwise increase in homozygous loci, whereas other mating systems exhibit a more gradual progression.

Number of homozygous loci over time, assuming recessive mutations and six loci. Mean and s.d. of 100 simulation trials. In androdioecious populations, 50% of hermaphrodites reproduce by mating and 50% by selfing. (A–C) First 50 generations. (D–F) All 500 generations.

While figure 3 shows that hermaphrodites reach the maximum of six homozygous loci within the simulated period, this is not observed in dioecious populations. Only three (or fewer) loci with large effects are homozygous for the mutated allele in strictly mating populations by the end of the simulations. Androdioecious populations display a similar pattern, remaining stagnant at three homozygous mutated loci for an extended duration. However, figure 3B shows that this threshold is eventually crossed, and these populations continue to accumulate mutated loci.

Initially, this appears to contradict the earlier observation of dioecious populations achieving resistance to the highest concentrations. Therefore, it is important to note here that in the same conditions, dioecious populations produce, on average, a greater number of offspring (1000 per hermaphrodite, so 500 per individual) than selfing hermaphrodites (300 per individual; see §2e) in a drug-free environment [56]. This difference in average egg count persists across genotypes and drug concentrations, with hermaphrodites always laying more eggs after mating compared to self-fertilization. It allows dioecious populations to have satisfactory fertility to meet the requirements for increasing drug concentration. Consequently, the selection pressure on dioecious populations for further evolution is reduced, decelerating the accumulation of mutations. At the same time, hermaphrodites, while having all six loci adapted, do not progress to higher concentrations. This is caused by the small effects of three of the loci. Even their combined effect is not sufficient to allow them to reach the egg count threshold required to proceed to higher concentrations (see the section on Genetics and sections below Sections 2.a, 2e, and 4e.).

The androdioecious populations also slow down mutation accumulation after the loci of large effect are fixed, but eventually accumulate mutations with small effects faster than dioecious populations, as selection pressure is higher. As we see in figure 2, large populations reach some kind of adaptation threshold after approximately 30 generations. This is because three loci with large effects are already mutated, and these three mutations are now fixed in the population. The populations progressed quickly into high concentration, where they were surviving but not doing well enough to proceed into higher concentrations (not producing an average egg count of at least 20). The remaining smaller mutations are not yet fixed, which we can see in figure 3B (the average number of mutated loci is just three). However, these mutations with small effects are present and slowly (because of their small effect) increase in frequency. After about 400 generations, their combined effect is strong enough to allow concentrations to increase (in some of the simulations), which, in turn, increases selection pressure, and the mutation accumulation of these three remaining mutations with small effects accelerates (as seen in figures 3E and 2E).

These observations have implications for the relevance of evolution experiments using C. elegans as a model organism for dioecious parasitic worms. Our results suggest that these experiments might overestimate the number of resistance mutations that would evolve in natural parasitic populations.

Androdioecy accelerates the fixation of beneficial mutations

(c)

If the environment changes rapidly, the beneficial allele can increase in frequency quickly. Figure 4 depicts the time to fixation (defined as reaching 95% frequency) of the first homozygous mutation across various mating systems, with the proportion of hermaphrodites engaging in mating represented by $[eqn]$ . This ranges from exclusive selfing in hermaphrodites ( $[eqn]$ , shown in blue), through androdioecious populations ( $[eqn]$ , green), to dioecious populations ( $[eqn]$ , red). Multiple panels within the figure illustrate results for different population sizes.

Time to fixation of the first beneficial allele for various fractions of mating hermaphrodites ξ . Resistance is encoded by six loci. Note the various scales on the Y axis. Box plots show the median, interquartile range and whiskers extending to 1.5× the interquartile range from the box edges.

Figure 4 demonstrates that while selfing results in the fastest fixation of the mutated allele in small populations, a combination of mating and selfing proves advantageous in larger populations. The proportion of mating that appears beneficial increases with population size. Furthermore, this relationship is influenced by dominance. For equivalent population sizes, a considerably lower mating frequency led to the shortest fixation time when mutations are dominant. In populations of 20 000 and 200 000 individuals, the shortest time to fixation was observed under different conditions. When mutations were dominant, the shortest time occurred when only 10% of hermaphrodites outcrossed. In contrast, when mutations were recessive, the shortest time occurred when 60 or 70% of hermaphrodites outcrossed, respectively. Again, this highlights the differences between C. elegans, which maintains males at relatively low frequencies, and dioecious parasitic worms, particularly when resistance mutations are dominant.

Dominant mutations lead to faster adaptation to higher drug concentration

(d)

Comparing figure 2 and the electronic supplementary material, figure S2 demonstrates that adaptation to higher concentrations is faster when mutations are dominant. This is expected, as dominant mutations are exposed to selection even in the heterozygous state, conferring an immediate fitness advantage.

The time to fixation of the first mutation (figure 4) is also generally longer for recessive mutations, especially in dioecious populations. This difference is less pronounced in androdioecious and hermaphrodite populations. We observed that the fastest adaptation always occurred for purely selfing populations or for populations with a high selfing rate. Figure 4 further reveals an unexpected interaction between population size and dominance. While recessive mutations require more time to fix in smaller populations, this is not the case for dominant mutations, where the differences between different population sizes were much smaller. This suggests that in smaller populations, mutations take longer to appear due to reduced mutation supply. However, once a beneficial mutation arises, it can reach prominence more rapidly due to stronger genetic drift in smaller populations. This effect is particularly pronounced for dominant mutations, which are immediately visible to selection, unlike recessive mutations that require homozygosity for phenotypic manifestation.

Mating hinders the fixation of mutations with small effects, but increases diversity

(e)

As mentioned previously, dioecious populations exhibited limited fixation of recessive mutations, with an average of only three loci homozygous for the recessive allele by the simulation’s end. To further investigate this phenomenon, we analysed the identity of the locus where the mutated allele first reached fixation (95% frequency). Electronic supplementary material, figures S3 and S4, reveal that in small, hermaphroditic populations, any of the six loci can be the first to reach fixation, suggesting a stochastic process driven by the random appearance of mutations. However, in larger populations, fixation consistently occurs first at one of the loci with a large effect. This indicates that in larger populations, the evolutionary trajectory is shaped primarily by competition between mutants, where those with larger effects have a selective advantage. Furthermore, both dominance and the presence of mating increased the likelihood of large-effect loci reaching fixation first, highlighting the influence of these factors on the adaptive process. Mating is known to have significant impacts on genetic diversity, with selfing leading to much lower diversity than mating [57,58]. Diversity is also expected to be lower in smaller populations. Figure 5 depicts the diversity, calculated as the Shannon index, over time (see electronic supplementary material, figure S5 for a zoomed-in version of the first 50 generations). In agreement with expectations, it shows that diversity is the lowest in the strictly selfing populations and the smallest populations. However, it shows that in androdioecious populations, diversity is as high as in strictly mating populations. In simulations where only two loci encode variance, the diversity is even higher in androdioecious populations than in dioecious ones (electronic supplementary material, figures S6 and S7). In all scenarios, we observed an initial steep increase in genetic diversity as new mutations appeared and increased in frequency. This was followed by a rapid decline as sensitive alleles were driven to extinction. The first peak in diversity can be attributed to the emergence and fixation of mutations with large effects. In contrast, the second peak arises from the increasing prevalence of mutations with smaller effects, which eventually also become fixed in the populations (except for dioecious ones). Since the selection advantage provided by these small-effect loci is weak, the adaptive process is slower, resulting in a broader peak. Notably, we did not observe a second peak in dioecious populations with recessive mutations, which aligns with our earlier finding that only three loci became fixed in this scenario.

Biodiversity changes over time. Recessive and dominant mutations, six loci encoding resistance. (A–C) Recessive mutations, and (D–F) dominant mutations. In androdioecious populations, 50% of hermaphrodites reproduce by mating and 50% by selfing.

Short-term observations can contradict the long-term ones

(g)

Simulations were conducted to investigate population dynamics over 500 generations, equivalent to more than four years of continuous in vitro experimentation and at least 80 years of continued in vivo selection of parasitic nematodes. In addition to being time-consuming, such experiments are likely to be extremely costly and unlikely to be conducted often. In practice, experiments containing approximately 50 or fewer generations are more feasible [15]. However, our findings demonstrate that observations from shorter timeframes may yield conclusions that differ significantly from those derived from longer-term investigations.

Figure 6 shows the drug concentration that was reached in the experiments with populations with a range of mating systems (from hermaphroditism, $[eqn]$ , to dioecy, $[eqn]$ ) after 10, 25, 100, 250 and 500 generations. ‘Drug concentration reached’ refers to a concentration at which the worms were able to survive in a particular generation but did not produce a sufficient egg count (20 on average), allowing us to increase the drug concentration further. The adaptive advantage conferred by various mating systems varies considerably across different time points and population sizes (see electronic supplementary material, figure S8 for bigger populations). Specifically, at early time points, hermaphroditic populations exhibit the most rapid adaptation, with sexual reproduction seemingly conferring a disadvantage across all population sizes. This trend, however, reverses over time as the benefits of outcrossing become increasingly evident, particularly in populations with recessive mutations (see electronic supplementary material, figure S9 for dominant mutations). These results highlight the critical importance of incorporating a temporal dimension into evolutionary studies, especially when considering the interplay between reproductive strategies and multi-locus evolution.

Drug concentration reached by populations at different time points, recessive mutations. Blue represents hermaphrodite populations, red represents dioecious ones and green represents various levels of androdioecy. Box plots show the median, interquartile range and whiskers extending to 1.5× the interquartile range from the box edges. Red line connects the means.

Discussion

Our study investigated how different mating systems influence the evolution of drug resistance in worm populations. We examined a range of mating systems, including dioecy, hermaphroditism and androdioecy as an intermediate state. Through simulations based on the biology of C. elegans and experiments conducted over 500 generations, we discovered that mating systems significantly impact evolutionary trajectories.

Hermaphrodites, capable of self-fertilization, exhibited faster adaptation under recessive mutations. Selfing promotes homozygosity, effectively exposing and fixing new mutations [36]. This is due to the increased probability of generating homozygous offspring from a single mutant individual, allowing for rapid fixation of the beneficial allele. In contrast, dioecious populations, reliant on outcrossing, faced a greater challenge in establishing recessive mutations, as these require the rare occurrence of mating between two carriers of the allele.

It is worth noting that the classic ‘cost of males’ does not apply in our model, as the presence of males actually leads to increased average egg production per individual. This increase in fecundity with the mating fraction $[eqn]$ arises from the fact that mated hermaphrodites lay more than thrice the number of eggs compared to those that reproduce through self-fertilization [56]. Consequently, the average fecundity of a population—or metapopulation, when focusing on a specific strain—depends on the proportion of mated hermaphrodites. This reflects the maximum fecundity in the absence of drugs, where the actual egg count is influenced by the fitness of the strain, which is determined by both genetic and environmental factors.

This deviation from traditional theoretical models results in altered selective pressures, whereby dioecious populations generate more offspring than hermaphrodites. This increased reproductive output enables dioecious populations to achieve higher drug concentrations despite initially adapting and accumulating mutations at a slower rate. Thus, over time, dioecious populations manage to survive and reproduce in environments with higher drug concentrations.

Androdioecious populations, characterized by a combination of selfing and outcrossing, appear to be able to access the advantages of both reproductive modes. Androdioecious populations adapted to increasing drug concentrations at a speed initially comparable to hermaphrodites and exceeded the speed of dioecious populations. While hermaphrodites lagged in their ability to reach higher concentrations, androdioecious populations, after a delay, successfully adapted to it. This pattern suggests that the flexibility afforded by androdioecy allows populations to rapidly exploit beneficial mutations through selfing while maintaining sufficient genetic diversity through outcrossing to overcome more long-term selective pressures. These results align with previous empirical observations of lower extinction rates in androdioecious lineages of C. elegans during experimental evolution in a new environment, supporting the adaptive benefits of combining selfing and outcrossing [59–61].

Our study allowed us to investigate a range of reproductive strategies, ranging from pure selfing to pure mating, with intermediate stages. We showed that the optimal balance between mating and selfing depends both on genetics—whether the mutations are dominant or recessive, as well as on the population size. A small fraction of males and mating may lead to the fastest adaptation. This is consistent with three independent studies that identified an advantage of outcrossing of androdioecious C. elegans populations during coevolution with pathogenic bacteria. Two of these studies included coevolving Bacillus thuringiensis as an antagonist and demonstrated that the favoured outcrossing rate never reached 100%, but remained at an intermediate level of 20−60% [60,62]. Outcrossing with males was directly shown to increase offspring resistance against the pathogen [60]. The third coevolution experiment, which included a C. elegans population with less genetic diversity than in the above studies and also a different pathogen, Serratia marcescens, found a steady increase in outcrossing rates from 20% to above 70% [61]. Overall, these findings strongly suggest an advantage of the intermediate outcrossing rates facilitated by the androdioecious reproductive system.

Our findings further highlight the role of population size in shaping adaptive outcomes. We observed that in small populations, any locus could be the first to reach fixation, indicating a stochastic process driven by the chance emergence of beneficial alleles. Conversely, in larger populations, adaptation becomes more predictable, with loci harbouring large-effect mutations consistently reaching fixation first. This shift suggests that competition between mutants, rather than mere chance, dictates the evolutionary trajectory in larger populations. Furthermore, if natural selection favours the recombination of advantageous alleles at two loci or if it favours heterozygosity (not investigated in this model), then adaptation should be faster in large populations. One example of such a scenario comes from previous evolution experiments, in which C. elegans coevolved with its pathogen B. thuringiensis and where reciprocal co-adaptations were favoured in larger host populations [63].

Most importantly, our study reveals the strong influence of that mating strategy on evolutionary dynamics and its complex interaction with other factors, such as population size and genetic structure. This has important implications for the use of C. elegans as a model species for the study of drug resistance in parasitic nematodes that are, to our knowledge, almost entirely dioecious. It suggests that androdioecious C. elegans benefits from the high selfing rates that lead to fast adaptation, especially in relatively short timescales typically used in evolutionary experiments. By generalizing the finding for C. elegans to other species, we may overestimate their ability to evolve resistance and, in turn, suggest non-optimal drug treatments with possible detrimental effects. At the same time, our findings indicate that androdiecious pathogenic helminths, should they emerge, will have a higher risk of evolving resistance against anthelmintics.

Our modelling study further highlights the necessity of theoretical approaches complementary to experimental observations, such as our model, in order to capture the full complexity of evolutionary dynamics. Mathematical models and simulations are especially valuable when conducting laboratory experiments that are impractical or unfeasible, such as when long time scales, large population sizes or high replication numbers are needed. In particular, our long-term simulations underscore the limitations of empirical studies when investigating the evolutionary processes across long time periods. While practical constraints often limit in vitro and in vivo experiments to a feasible timeframe of 50 or fewer generations, our model reveals that different outcomes can be observed over periods of many hundreds of generations. This is particularly evident in the context of reproductive mode, where the apparent short-term advantage of hermaphroditism can give way to the long-term benefits of outcrossing.

It is important to acknowledge certain limitations of our model. We assumed unlinked autosomal genes, which may not fully reflect the complexities of real-world genomes and considering the linkage between the genes is an important expansion of our model. Furthermore, we assumed that both the benefits and costs of individual mutations combine in an additive manner. While this is a common assumption in many polygenic models, several recent studies suggest that mutations may interact in a multiplicative way instead. Using a multiplicative approach can be beneficial, especially for cost, because it prevents negative fitness values, even when multiple costly mutations are combined, thus allowing for the simulation of such mutations. However, we used an additive model to avoid introducing epistasis into our analysis. Epistasis occurs when the effect of one mutation depends on the presence of another, which would complicate the interpretation of our results and obscure the direct effects of the mutations themselves. While we recognize the significance of epistasis in evolutionary processes, our primary goal for this research was to concentrate on other effects, and the consideration of epistasis is a promising extension for future research.

The influence of drug concentration on egg count was modelled solely through hermaphrodites, neglecting any potential role of male fitness. Specifically, we simplified the impact of anthelmintics on worm fitness by focusing solely on reduced egg count, regardless of the underlying causes—such as a shorter lifespan of the hermaphrodite, lower fecundity or decreased likelihood of eggs hatching. We did not consider other factors, like a reduced chance of finding a mating partner. While this simplification was necessary, the explicit consideration of a larger number of factors represents a worthwhile avenue for future work. In particular, future studies could address these limitations by incorporating linked genes, exploring the impact of male fitness on adaptation and extending the model to encompass other reproductive modes, such as asexual reproduction observed in some parasitic worm species.

Finally, we note that the gradually increasing drug concentrations modelled here are more representative of the laboratory experiments than of the anthelmintic use in real populations. However, as our main objective was to discuss the suitability of C. elegans for studies of drug resistance evolution in the laboratory setting, we chose to gradually increase drug concentration as it agrees well with methods employed by most researchers in order to select resistant parasite strains in animal experiments. Moreover, we believe that this study is also relevant for natural populations. Pathogens are often exposed to lower effective concentrations than those that would kill them. For instance, reinfection can occur when drug concentrations have already dropped to levels that do not completely prevent further parasite development but still select for certain parasite genotypes. In addition, drug concentration is not constant throughout the human and animal body, and some worms may experience higher concentrations, while some lower. The differences also occur during the course of drug metabolism and between hosts who differ in their metabolism and uptake of the drugs.

In conclusion, our study demonstrates the complex interplay between mating systems, population size and genetic architecture in driving the evolution of drug resistance in worms. While hermaphroditism initially facilitates rapid adaptation through selfing, dioecy and androdioecy ultimately prevail in high drug concentrations due to increased genetic diversity and offspring production. Notably, androdioecy, by combining selfing and outcrossing, appears to offer the most robust adaptive strategy. These findings underscore the importance of long-term evolutionary perspectives and integrated theoretical approaches for understanding and predicting drug resistance dynamics in parasitic worms, with implications for effective control strategies. Our model, in combination with the experimental system, could be used for the testing of new drugs and the general potential of nematodes to evolve resistance under a variety of conditions. This approach should be expanded in the future to include nematodes with different reproductive strategies, especially those of parasitic worms, and population structures that permit migration and gene flow, reflecting natural complexity.