<figure><img src="image_1.png" /></figure>

HAL
open science

# The impact of environmental fluctuations, sexual dimorphism, dominance reversal and plasticity on the pigmentation-related genetic and phenotypic variation in D. melanogaster populations - A modelling study

Laurent Freoa, Jean-Michel Gibert, Amandine Véber

► To cite this version:

Laurent Freoa, Jean-Michel Gibert, Amandine Véber. The impact of environmental fluctuations, sexual dimorphism, dominance reversal and plasticity on the pigmentation-related genetic and phenotypic variation in D. melanogaster populations - A modelling study. 2024. hal-04597232

HAL Id: hal-04597232

https://hal.science/hal-04597232

Preprint submitted on 2 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The impact of environmental fluctuations, sexual dimorphism, dominance reversal and plasticity on the pigmentation-related genetic and phenotypic variation in D. melanogaster populations - A modelling study

Laurent Freoa, Jean-Michel Gibert, Amandine Véber

► To cite this version:

Laurent Freoa, Jean-Michel Gibert, Amandine Véber. The impact of environmental fluctuations, sexual dimorphism, dominance reversal and plasticity on the pigmentation-related genetic and phenotypic variation in D. melanogaster populations - A modelling study. 2024. hal-04597232

HAL Id: hal-04597232

https://hal.science/hal-04597232v1

Preprint submitted on 2 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The impact of environmental fluctuations, sexual dimorphism, dominance reversal and plasticity on the pigmentation-related genetic and phenotypic variation in *D. melanogaster* populations – A modelling study.

Laurent Freoa<fnref n="1" />,<fnref n="2" />, Jean-Michel Gibert<fnref n="2" />, and Amandine Véber*,<fnref n="1" />

June 2, 2024

## Abstract

In this article, we consider a phenotypic trait which is variable only in females and depends both on the individual's allelic state at a given locus and on environmental conditions at birth. Whatever their genotype and environment, males express the same phenotype. The key example of such traits that we study here is the pigmentation of the posterior abdomen in *Drosophila melanogaster*, which, in females, is influenced by the genotype of the individual and by the environmental temperature during development (and is approximately constant in males). In the absence of sexual dimorphism and with selection acting on a trait, it has been shown that periodically fluctuating environments can maintain genetic variation at the underlying locus. Here, we introduce a more complex model integrating several important features that can also influence the genetic and phenotypic composition of the population, in order to partially disentangle their effects on the maintenance of variation at the genetic and phenotypic levels. In this model, selection acts on the trait, the population is sexually dimorphic, individuals are diploid, the trait considered is plastic and dominance reversal renders the contribution of the genotype to the phenotype dependent on the environment (i.e., an individual heterozygous at the locus of interest tends to have a trait value which is favourable in the environment in which it was born). Because drosophila populations experience massive seasonal fluctuations in size due to variations in resource availability and environmental conditions, we model a spatially structured population with a local regulation of population sizes and whose demography exhibits boom and bust dynamics.

Using simulations of a rather caricatural model of thermal adaptation based on abdominal pigmentation in *D. melanogaster* that constitute our case study, we find that, in periodically fluctuating environments, the combination of plasticity, dominance reversal and sexual dimorphism leads to oscillations of the mean female phenotype but not of allele frequencies, which tend to stabilise around some random value away from 0 and 1 (that is, genetic variability within the population is conserved). Plasticity and dominance reversal allow gradual but fast phenotypic adaptation in the female population, which is all the more efficient as the period of environment oscillations is long. However, the absence of sexual dimorphism leads to higher extinction probabilities, as not enough individuals exhibit the best adapted phenotype during the periods when surviving tough environmental conditions is crucial. Conversely, sexual dimorphism without plasticity and dominance reversal leads to phenotypic variability but no oscillations of the female population mean phenotype. The male population thus serves as a reserve of individuals with optimal phenotypes to survive harsh conditions and restart the population when better times are coming; because male phenotypes are independent of their genotypes, this reserve also shelters alleles associated with less adapted phenotypes, thereby helping to maintain genetic variability in the whole population.

**Keywords:** birth-death process, locally regulated population sizes, fluctuating environment, dominance reversal, plasticity.

<footnote n="1">¹Université Paris Cité, CNRS, MAP5, 45 rue des Saints-Pères 75006 Paris, France</footnote>

<footnote n="2">²Laboratoire de Biologie du Développement, UMR 7622, CNRS, Institut de Biologie Paris-Seine (IBPS), Sorbonne Université, 9 Quai St-Bernard, 75005, Paris, France.</footnote>

<footnote>*Corresponding author: amandine.veber@parisdescartes.fr</footnote>
# The impact of environmental fluctuations, sexual dimorphism, dominance reversal and plasticity on the pigmentation-related genetic and phenotypic variation in *D. melanogaster* populations – A modelling study.

Laurent Freoa<fnref n="1" />,<fnref n="2" />, Jean-Michel Gibert<fnref n="2" />, and Amandine Véber*,<fnref n="1" />

June 2, 2024

## Abstract

In this article, we consider a phenotypic trait which is variable only in females and depends both on the individual's allelic state at a given locus and on environmental conditions at birth. Whatever their genotype and environment, males express the same phenotype. The key example of such traits that we study here is the pigmentation of the posterior abdomen in *Drosophila melanogaster*, which, in females, is influenced by the genotype of the individual and by the environmental temperature during development (and is approximately constant in males). In the absence of sexual dimorphism and with selection acting on a trait, it has been shown that periodically fluctuating environments can maintain genetic variation at the underlying locus. Here, we introduce a more complex model integrating several important features that can also influence the genetic and phenotypic composition of the population, in order to partially disentangle their effects on the maintenance of variation at the genetic and phenotypic levels. In this model, selection acts on the trait, the population is sexually dimorphic, individuals are diploid, the trait considered is plastic and dominance reversal renders the contribution of the genotype to the phenotype dependent on the environment (i.e., an individual heterozygous at the locus of interest tends to have a trait value which is favourable in the environment in which it was born). Because drosophila populations experience massive seasonal fluctuations in size due to variations in resource availability and environmental conditions, we model a spatially structured population with a local regulation of population sizes and whose demography exhibits boom and bust dynamics.

Using simulations of a rather caricatural model of thermal adaptation based on abdominal pigmentation in *D. melanogaster* that constitute our case study, we find that, in periodically fluctuating environments, the combination of plasticity, dominance reversal and sexual dimorphism leads to oscillations of the mean female phenotype but not of allele frequencies, which tend to stabilise around some random value away from 0 and 1 (that is, genetic variability within the population is conserved). Plasticity and dominance reversal allow gradual but fast phenotypic adaptation in the female population, which is all the more efficient as the period of environment oscillations is long. However, the absence of sexual dimorphism leads to higher extinction probabilities, as not enough individuals exhibit the best adapted phenotype during the periods when surviving tough environmental conditions is crucial. Conversely, sexual dimorphism without plasticity and dominance reversal leads to phenotypic variability but no oscillations of the female population mean phenotype. The male population thus serves as a reserve of individuals with optimal phenotypes to survive harsh conditions and restart the population when better times are coming; because male phenotypes are independent of their genotypes, this reserve also shelters alleles associated with less adapted phenotypes, thereby helping to maintain genetic variability in the whole population.

**Keywords:** birth-death process, locally regulated population sizes, fluctuating environment, dominance reversal, plasticity.

<footnote n="1">¹Université Paris Cité, CNRS, MAP5, 45 rue des Saints-Pères 75006 Paris, France</footnote>

<footnote n="2">²Laboratoire de Biologie du Développement, UMR 7622, CNRS, Institut de Biologie Paris-Seine (IBPS), Sorbonne Université, 9 Quai St-Bernard, 75005, Paris, France.</footnote>

<footnote>*Corresponding author: amandine.veber@parisdescartes.fr</footnote>
# 1 Introduction

In many animal species, body pigmentation exhibits high variability, stemming either from genetic diversity or from phenotypic plasticity. Several mechanisms associate a selective advantage to color variation: mate recognition, involving sexual dimorphism [25, 41, 42, 48, 61], physiological advantage in ectotherms linked to thermal balance and light radiation absorption [34, 37, 45, 67], longevity, fecundity, production of cuticular hydrocarbons, and resistance against pathogens, parasites, UV or desiccation [4, 24, 46, 52, 55, 58, 69].

Insect melanism has long served as a case study in evolutionary biology, offering examples of natural selection [44], genetic regulation [70], and pleiotropic effects [69]. Melanism stands out as one of the simplest and most prevalent instances of biodiversity in nature [64], manifesting both as intraspecific polymorphisms and fixed differences among closely related species, significantly influencing various aspects of insect biology. *Drosophila melanogaster* exhibits sexual dimorphism in pigmentation of the posterior abdomen (tergites 5 and 6): males possess a black posterior abdomen, while females display yellow tergites with a black stripe at the posterior margin. The extension of black pigmentation on female tergites is highly variable and depends on the genetic variation of a small number of genes and on the environmental temperature during pupal development [3, 21, 30]. The relation between body pigmentation and adaptation to environmental temperature and light in ectotherms as small as drosophilae has long been questioned. In [68], it was argued that pigmentation did not significantly affect the body temperature of small insects exposed to sunlight. This conclusion was challenged in [29], where a thermal camera was used to show that mutants or species of *Drosophila* with distinct pigmentation exhibited statistically significant differences in the rate and level at which individuals warmed when they were exposed to a light source mimicking the sunlight (with warming differences lower than one Celsius degree, hence the difficulty to measure them with older calorimetry tools). Consequently, even in small insects, pigmentation appears to be involved in thermoregulation, with fully pigmented flies reaching a slightly higher body temperature under sunlight compared to less pigmented flies [29].

In temperate areas, drosophila populations live in seasonally fluctuating environments, with abundant food and warm temperatures leading to a large increase in local population sizes around summer, followed by more severe living conditions causing a sharp decline in population size around winter. In this work, we aim at disentangling the roles of the genetic and environmental components of body pigmentation in the maintenance of phenotypic diversity in *D. melanogaster*, taking into account the species ecology as well as specific genotype-phenotype features impacting the distribution of phenotypes in such populations, namely sexual dimorphism, plasticity and dominance reversal (we recall the definition of these terms thereafter). We also investigate the patterns of genetic diversity, and in particular the presence or absence of fluctuations in allele frequencies, in this situation where different forces go in opposite directions. Indeed, selection on an environmentally sensitive trait, coupled with environmental fluctuations, can lead to oscillations in allele frequencies at the loci underlying the trait under selection. On the other hand, plasticity and dominance reversal can limit individual maladaptation by producing phenotypes that are favourable despite the individual carrying alleles that are, in principle, associated with deleterious phenotypes. This may result in reduced fluctuations in (or even stabilisation of) allele frequencies without leading to the extinction of some alleles or of the population itself. Because we do not assume selection due to mating preferences, the main effect that may be expected from the particular form of sexual dimorphism considered here (and all the more so when it is combined with dominance reversal) is that alleles associated with currently deleterious phenotypes can remain present in the male population at low cost, until the environment comes back to a state in which the phenotypes they produce are favourable. This phenomenon is sometimes referred to as a *genomic storage effect* [40]. In this case, the male population acts as a genetic reserve that buffers the allele frequency oscillations generated by the female population, which is subject to fluctuating selection due to alternating environments. Which of these effects has the largest impact on the genetic diversity seen in the whole population is one of the key questions we pose here.

Cyclic changes in allele frequencies due to strong selection in a seasonally fluctuating environment have been extensively studied from an empirical and theoretical point of view, starting from the seminal work of Timofeef-Ressovsky on the beetle *Adalia bipunctata* [63] and of Dobzhansky on
# 1 Introduction

In many animal species, body pigmentation exhibits high variability, stemming either from genetic diversity or from phenotypic plasticity. Several mechanisms associate a selective advantage to color variation: mate recognition, involving sexual dimorphism [25, 41, 42, 48, 61], physiological advantage in ectotherms linked to thermal balance and light radiation absorption [34, 37, 45, 67], longevity, fecundity, production of cuticular hydrocarbons, and resistance against pathogens, parasites, UV or desiccation [4, 24, 46, 52, 55, 58, 69].

Insect melanism has long served as a case study in evolutionary biology, offering examples of natural selection [44], genetic regulation [70], and pleiotropic effects [69]. Melanism stands out as one of the simplest and most prevalent instances of biodiversity in nature [64], manifesting both as intraspecific polymorphisms and fixed differences among closely related species, significantly influencing various aspects of insect biology. *Drosophila melanogaster* exhibits sexual dimorphism in pigmentation of the posterior abdomen (tergites 5 and 6): males possess a black posterior abdomen, while females display yellow tergites with a black stripe at the posterior margin. The extension of black pigmentation on female tergites is highly variable and depends on the genetic variation of a small number of genes and on the environmental temperature during pupal development [3, 21, 30]. The relation between body pigmentation and adaptation to environmental temperature and light in ectotherms as small as drosophilae has long been questioned. In [68], it was argued that pigmentation did not significantly affect the body temperature of small insects exposed to sunlight. This conclusion was challenged in [29], where a thermal camera was used to show that mutants or species of *Drosophila* with distinct pigmentation exhibited statistically significant differences in the rate and level at which individuals warmed when they were exposed to a light source mimicking the sunlight (with warming differences lower than one Celsius degree, hence the difficulty to measure them with older calorimetry tools). Consequently, even in small insects, pigmentation appears to be involved in thermoregulation, with fully pigmented flies reaching a slightly higher body temperature under sunlight compared to less pigmented flies [29].

In temperate areas, drosophila populations live in seasonally fluctuating environments, with abundant food and warm temperatures leading to a large increase in local population sizes around summer, followed by more severe living conditions causing a sharp decline in population size around winter. In this work, we aim at disentangling the roles of the genetic and environmental components of body pigmentation in the maintenance of phenotypic diversity in *D. melanogaster*, taking into account the species ecology as well as specific genotype-phenotype features impacting the distribution of phenotypes in such populations, namely sexual dimorphism, plasticity and dominance reversal (we recall the definition of these terms thereafter). We also investigate the patterns of genetic diversity, and in particular the presence or absence of fluctuations in allele frequencies, in this situation where different forces go in opposite directions. Indeed, selection on an environmentally sensitive trait, coupled with environmental fluctuations, can lead to oscillations in allele frequencies at the loci underlying the trait under selection. On the other hand, plasticity and dominance reversal can limit individual maladaptation by producing phenotypes that are favourable despite the individual carrying alleles that are, in principle, associated with deleterious phenotypes. This may result in reduced fluctuations in (or even stabilisation of) allele frequencies without leading to the extinction of some alleles or of the population itself. Because we do not assume selection due to mating preferences, the main effect that may be expected from the particular form of sexual dimorphism considered here (and all the more so when it is combined with dominance reversal) is that alleles associated with currently deleterious phenotypes can remain present in the male population at low cost, until the environment comes back to a state in which the phenotypes they produce are favourable. This phenomenon is sometimes referred to as a *genomic storage effect* [40]. In this case, the male population acts as a genetic reserve that buffers the allele frequency oscillations generated by the female population, which is subject to fluctuating selection due to alternating environments. Which of these effects has the largest impact on the genetic diversity seen in the whole population is one of the key questions we pose here.

Cyclic changes in allele frequencies due to strong selection in a seasonally fluctuating environment have been extensively studied from an empirical and theoretical point of view, starting from the seminal work of Timofeef-Ressovsky on the beetle *Adalia bipunctata* [63] and of Dobzhansky on
*Drosophila pseudoobscura* [23, 72]. Advances in fast genome sequencing led to the uncovering of hundreds of polymorphisms in *Drosophila* populations whose frequencies oscillate in a seasonal way [5, 50], some of them underlying fitness-related traits, and recent efforts to produce bioinformatics pipelines to integrate different types of data sets [43] will probably lead to a better understanding of the spatio-temporal genetic patterns and evolutionary dynamics of *D. melanogaster* populations. At the phenotypic level, cyclic fluctuations in life-history traits or in traits related to tolerance to different forms of stress have also been reported (see [1] and references therein). From the theoretical point of view, the long term maintenance of genetic and phenotypic diversity in multivoltine species like *Drosophila* species can be explained by different combinations of evolutionary forces such as balancing selection, polygenic adaptation, gene regulation, sexual antagonism, or temporally varying selection, as reviewed in [27, 36]. In particular, in the last decade, the role of temporally fluctuating environments in the maintenance of diversity (at the population and species scales, and for different organisms) has been the object of renewed attention [73]. Of particular interest are the impact of the period and amplitude of the environmental fluctuations on selection gradients [66], the conditions under which the interplay between selection and environmental variations generate stable, or oscillating, polymorphisms [12], and the impact of the structuring of the population into different classes with distinct responses to selection [49]. To our knowledge, only asexual populations are considered in these models. Empirical observations of *adaptive tracking*, that is, the continuous adaptation in response to rapid environmental change, have recently been reported in *D. melanogaster* [50, 60] and in other diptera species [56].

Phenotypic plasticity, defined as the property of a given genotype to produce different phenotypes in response to distinct environmental conditions [57], tends to partly decorrelate genetic and phenotypic diversity by allowing genotypes *a priori* associated with less favourable trait values to give rise to phenotypes which are adapted, or at least less maladapted (and conversely). In drosophilae, body pigmentation plasticity is primarily associated with developmental temperature: lower temperatures during pupal development result in increasingly darker flies, aligning with the thermal budget adaptive hypothesis [18, 32, 67]. Plasticity itself is a heritable and adaptive trait [47]; it may evolve in variable environments due to selection pressures, for instance related to the variability and reliability of environmental cues [8, 40]. In models and empirical investigations, a standard proxy for plasticity is the maximum slope of the reaction norm, and this is also our choice for the model we propose later. Another important question we want to address with this model is to evaluate the relative contributions of phenotypic plasticity and of allele frequency variations in the pigmentation distribution in a *D. melanogaster* population living in a periodically fluctuating environment. Because the number of mechanisms that we integrate in the model is already large, we do not allow plasticity to evolve and we leave this interesting question for future investigations.

Dominance reversal is another mechanism which was identified as being potentially important in the maintenance of genetic diversity in fluctuating environments, at least in theory [6, 15, 38, 71]. It is defined as the property that alleles at a given locus may be dominant in selective contexts in which they are favourable, and recessive in contexts in which they are deleterious. In the case of body pigmentation in *D. melanogaster*, dominance reversal corresponds to heterozygotes at the loci involved in pigmentation being darker when environmental temperatures during pupal development are lower, and lighter when developing at warmer temperatures. Theoretical exploration of the consequences of dominance reversals on genetic diversity led to the concept of *segregation lift* [71] (as opposed to *segregation load*, corresponding to the cost of segregation at selectively overdominant loci), since segregation at the loci of interest may be advantageous at the population level due to an increased fraction of the population displaying favourable phenotypes. Empirical evidence of dominance reversal in natural populations is scarce and difficult to obtain [38], but temperature dependent dominance was identified in *D. melanogaster* [10].

Finally, the role of a spatially heterogeneous and fluctuating environment in the maintenance of local genetic and phenotypic diversity has also been theoretically investigated [11, 39]. Indeed, cyclic selection that spatially varies in magnitude, coupled with migration, can lead to the creation of local reservoirs of alleles giving rise to phenotypes which are not presently favourable in neighbouring areas, but will be in a close future (either due to the fact that genetic drift acts less efficiently in densely
*Drosophila pseudoobscura* [23, 72]. Advances in fast genome sequencing led to the uncovering of hundreds of polymorphisms in *Drosophila* populations whose frequencies oscillate in a seasonal way [5, 50], some of them underlying fitness-related traits, and recent efforts to produce bioinformatics pipelines to integrate different types of data sets [43] will probably lead to a better understanding of the spatio-temporal genetic patterns and evolutionary dynamics of *D. melanogaster* populations. At the phenotypic level, cyclic fluctuations in life-history traits or in traits related to tolerance to different forms of stress have also been reported (see [1] and references therein). From the theoretical point of view, the long term maintenance of genetic and phenotypic diversity in multivoltine species like *Drosophila* species can be explained by different combinations of evolutionary forces such as balancing selection, polygenic adaptation, gene regulation, sexual antagonism, or temporally varying selection, as reviewed in [27, 36]. In particular, in the last decade, the role of temporally fluctuating environments in the maintenance of diversity (at the population and species scales, and for different organisms) has been the object of renewed attention [73]. Of particular interest are the impact of the period and amplitude of the environmental fluctuations on selection gradients [66], the conditions under which the interplay between selection and environmental variations generate stable, or oscillating, polymorphisms [12], and the impact of the structuring of the population into different classes with distinct responses to selection [49]. To our knowledge, only asexual populations are considered in these models. Empirical observations of *adaptive tracking*, that is, the continuous adaptation in response to rapid environmental change, have recently been reported in *D. melanogaster* [50, 60] and in other diptera species [56].

Phenotypic plasticity, defined as the property of a given genotype to produce different phenotypes in response to distinct environmental conditions [57], tends to partly decorrelate genetic and phenotypic diversity by allowing genotypes *a priori* associated with less favourable trait values to give rise to phenotypes which are adapted, or at least less maladapted (and conversely). In drosophilae, body pigmentation plasticity is primarily associated with developmental temperature: lower temperatures during pupal development result in increasingly darker flies, aligning with the thermal budget adaptive hypothesis [18, 32, 67]. Plasticity itself is a heritable and adaptive trait [47]; it may evolve in variable environments due to selection pressures, for instance related to the variability and reliability of environmental cues [8, 40]. In models and empirical investigations, a standard proxy for plasticity is the maximum slope of the reaction norm, and this is also our choice for the model we propose later. Another important question we want to address with this model is to evaluate the relative contributions of phenotypic plasticity and of allele frequency variations in the pigmentation distribution in a *D. melanogaster* population living in a periodically fluctuating environment. Because the number of mechanisms that we integrate in the model is already large, we do not allow plasticity to evolve and we leave this interesting question for future investigations.

Dominance reversal is another mechanism which was identified as being potentially important in the maintenance of genetic diversity in fluctuating environments, at least in theory [6, 15, 38, 71]. It is defined as the property that alleles at a given locus may be dominant in selective contexts in which they are favourable, and recessive in contexts in which they are deleterious. In the case of body pigmentation in *D. melanogaster*, dominance reversal corresponds to heterozygotes at the loci involved in pigmentation being darker when environmental temperatures during pupal development are lower, and lighter when developing at warmer temperatures. Theoretical exploration of the consequences of dominance reversals on genetic diversity led to the concept of *segregation lift* [71] (as opposed to *segregation load*, corresponding to the cost of segregation at selectively overdominant loci), since segregation at the loci of interest may be advantageous at the population level due to an increased fraction of the population displaying favourable phenotypes. Empirical evidence of dominance reversal in natural populations is scarce and difficult to obtain [38], but temperature dependent dominance was identified in *D. melanogaster* [10].

Finally, the role of a spatially heterogeneous and fluctuating environment in the maintenance of local genetic and phenotypic diversity has also been theoretically investigated [11, 39]. Indeed, cyclic selection that spatially varies in magnitude, coupled with migration, can lead to the creation of local reservoirs of alleles giving rise to phenotypes which are not presently favourable in neighbouring areas, but will be in a close future (either due to the fact that genetic drift acts less efficiently in densely
populated regions, and therefore locally maladapted phenotypes and genotypes are not totally wiped
out before the environment changes and they become adapted again; or due to the fact that selection
does not have the same direction in all regions of space and maladapted individuals in one area may
be adapted in other areas). The effect is called the *spatial storage effect* [11, 39]. Note that *ecological
storage effects* were also identified in models for populations with *boom and bust* dynamics [6], such as
local populations of *D. melanogaster*, and the two effects may combine to promote genetic diversity in
drosophila populations considered over large, inhomogeneous, regions of space. Massive genomic data
have been collected in recent years to investigate these questions in real populations [13, 43, 59].

To incorporate all these features into a model, we use the standard approach of defining a stochastic
(individual-based) birth-death model, in which each individual is represented by a spatial location,
a genotype and a phenotype and is male or female. We represent the state of the whole population
at a given time *t* by a pair ($v_t^m, v_t^f$) corresponding to the empirical distributions over the set of all
possible spatial locations, genotypes, and phenotypes, of the subpopulations of males and females
considered separately. This approach finds its origins in work on spatially structured populations
such as [7, 22], and was made rigorous in a probabilistic setting in [28]. It is now extensively used
for the modelling and analysis of structured populations, as it constitutes a very flexible framework
to encode and analyse dynamics in which only the way individuals are distributed over the space of
characteristics matters (and not the precise labelling of individuals), and allows to take the ecology
as well as the stochasticity inherent to finite populations into account. Furthermore, it comes with
a well-developed toolbox for mathematical analysis, for the derivation of limiting dynamics (as the
population size tends to infinity, in particular) [53] and for simulation [9, 35]. Originally used to
model asexual populations, these measure-valued processes were then generalised to incorporate the
case of sexual reproduction, in both haploid [17, 62] and diploid populations [14, 16, 54]. Our choice
of birth-death dynamics departs from the classical population genetics (*Wright-Fisher*-type) approach
that most of the previously cited work use, in which population size is held constant by assuming that
every birth is compensated by the death of another individual. In our model, birth and death rates
are modulated by the environment, and a local regulation of population sizes is enforced through a
density-dependent competition term.

Another biological field and mathematical framework in which the relative contributions of the
genetic and environmental components of a phenotypic trait have been extensively studied is quanti-
tative genetics. In this family of models, population sizes are allowed to vary but one usually assumes
that the trait value is the sum of a genetic part, which results from the contribution of many loci
[26, 65], and of an environmental noise encompassing all factors that are not heritable. In the extreme
case of Fisher's *infinitesimal model* [26], where an 'infinite' number of loci each have an 'infinitesimal'
contribution to the trait, the genetic part of the value of the offspring trait knowing the parental
traits follows a Gaussian distribution [2, 26]. Based on this approximation, the impact of sexual re-
production and selection acting on the trait in an inhomogeneous environment has been the object of
recent theoretical work [20]. Here we consider the other extreme in which the trait value depends on a
single locus, caricaturing the fact that only few genes are involved in female abdominal pigmentation
variation in *D. melanogaster*.

The model is introduced in Section 2. In Section 3, using simulations we investigate, separately and
in different combinations, the impact of sexual dimorphism, dominance reversal and plasticity on the
phenotypic and genetic diversity when the population of interest experiences periodic environmental
fluctuations and selection acts on the phenotypic trait. We do not consider spatially heterogeneous
selection, as many scenarios would then have to be explored and this requires a full study on its own.
Finally, these results are discussed in Section 4.

## 2 The stochastic individual-based model

Let us first set the main notation. Each individual in our population has a spatial location *x* taking
values in a bounded domain *X* of R<sup>2</sup>, a genotype *g* taking values in a finite set *G* = {*g*<sub>1</sub>, ..., *g*<sub>*K*</sub>}, and
a phenotype *u*<sub>*p*</sub> taking values in a compact set U<sub>*p*</sub> of R<sup>q</sup>. When convenient, we shall abbreviate the
populated regions, and therefore locally maladapted phenotypes and genotypes are not totally wiped
out before the environment changes and they become adapted again; or due to the fact that selection
does not have the same direction in all regions of space and maladapted individuals in one area may
be adapted in other areas). The effect is called the *spatial storage effect* [11, 39]. Note that *ecological
storage effects* were also identified in models for populations with *boom and bust* dynamics [6], such as
local populations of *D. melanogaster*, and the two effects may combine to promote genetic diversity in
drosophila populations considered over large, inhomogeneous, regions of space. Massive genomic data
have been collected in recent years to investigate these questions in real populations [13, 43, 59].

To incorporate all these features into a model, we use the standard approach of defining a stochastic
(individual-based) birth-death model, in which each individual is represented by a spatial location,
a genotype and a phenotype and is male or female. We represent the state of the whole population
at a given time *t* by a pair ($v_t^m, v_t^f$) corresponding to the empirical distributions over the set of all
possible spatial locations, genotypes, and phenotypes, of the subpopulations of males and females
considered separately. This approach finds its origins in work on spatially structured populations
such as [7, 22], and was made rigorous in a probabilistic setting in [28]. It is now extensively used
for the modelling and analysis of structured populations, as it constitutes a very flexible framework
to encode and analyse dynamics in which only the way individuals are distributed over the space of
characteristics matters (and not the precise labelling of individuals), and allows to take the ecology
as well as the stochasticity inherent to finite populations into account. Furthermore, it comes with
a well-developed toolbox for mathematical analysis, for the derivation of limiting dynamics (as the
population size tends to infinity, in particular) [53] and for simulation [9, 35]. Originally used to
model asexual populations, these measure-valued processes were then generalised to incorporate the
case of sexual reproduction, in both haploid [17, 62] and diploid populations [14, 16, 54]. Our choice
of birth-death dynamics departs from the classical population genetics (*Wright-Fisher*-type) approach
that most of the previously cited work use, in which population size is held constant by assuming that
every birth is compensated by the death of another individual. In our model, birth and death rates
are modulated by the environment, and a local regulation of population sizes is enforced through a
density-dependent competition term.

Another biological field and mathematical framework in which the relative contributions of the
genetic and environmental components of a phenotypic trait have been extensively studied is quanti-
tative genetics. In this family of models, population sizes are allowed to vary but one usually assumes
that the trait value is the sum of a genetic part, which results from the contribution of many loci
[26, 65], and of an environmental noise encompassing all factors that are not heritable. In the extreme
case of Fisher's *infinitesimal model* [26], where an 'infinite' number of loci each have an 'infinitesimal'
contribution to the trait, the genetic part of the value of the offspring trait knowing the parental
traits follows a Gaussian distribution [2, 26]. Based on this approximation, the impact of sexual re-
production and selection acting on the trait in an inhomogeneous environment has been the object of
recent theoretical work [20]. Here we consider the other extreme in which the trait value depends on a
single locus, caricaturing the fact that only few genes are involved in female abdominal pigmentation
variation in *D. melanogaster*.

The model is introduced in Section 2. In Section 3, using simulations we investigate, separately and
in different combinations, the impact of sexual dimorphism, dominance reversal and plasticity on the
phenotypic and genetic diversity when the population of interest experiences periodic environmental
fluctuations and selection acts on the phenotypic trait. We do not consider spatially heterogeneous
selection, as many scenarios would then have to be explored and this requires a full study on its own.
Finally, these results are discussed in Section 4.

## 2 The stochastic individual-based model

Let us first set the main notation. Each individual in our population has a spatial location *x* taking
values in a bounded domain *X* of R<sup>2</sup>, a genotype *g* taking values in a finite set *G* = {g<sub>1</sub>, ..., g<sub>K</sub>}, and
a phenotype *u*<sub>*p*</sub> taking values in a compact set U<sub>*p*</sub> of R<sup>q</sup>. When convenient, we shall abbreviate the
set of characteristics of an individual as $c = (x, g, u_p)$, this variable thus belonging to

$$ \mathcal{C} := \mathcal{X} \times \mathcal{G} \times \mathcal{U}_p. \qquad (1) $$

In the case study presented in Section 3, we take $\mathcal{X} = (-1000, 1000)^2$, $\mathcal{G} = \{AA, Aa, aa\}$ (corresponding to a single biallelic locus in a diploid population), and the phenotype is a one-dimensional variable taking values in $\mathcal{U}_p = [0, 10]$ and describing the intensity of abdominal pigmentation.

Because we want to explicitly model density dependence on both sexes for reproduction, as well as sexual dimorphism, it is more convenient to consider the male and female subpopulations separately. Mathematically speaking, at any time $t \ge 0$, each of these subpopulations is represented by its empirical distribution, that is, the counting measure on $\mathcal{C}$ giving mass 1 to the triplet of characteristics $c = (x, g, u_p)$ of each individual alive at time $t$:

$$ \nu_t^m = \sum_{i \in V_t^m} \delta_{(x_t^i, g^i, u_p^i)} \quad \text{and} \quad \nu_t^f = \sum_{i \in V_t^f} \delta_{(x_t^i, g^i, u_p^i)}, \qquad (2) $$

where $V_t^m$ and $V_t^f$ are respectively the set of labels of males and females present in the population at time $t$, $\delta_c$ denotes a Dirac mass at $c$, $x_t^i$ is the position at time $t$ of individual $i$ (which will change through time, hence the dependence on $t$), and $g^i$, $u_p^i$ are the genotype and phenotype of individual $i$ (supposed to be fixed throughout its life). We write $\mathcal{M}$ for the set of all finite counting measures on $\mathcal{C}$. The state of the whole population at time $t$ is then described by the pair $(\nu_t^m, \nu_t^f) \in \mathcal{M}^2$, and we are interested in the dynamics of the process $(\nu_t^m, \nu_t^f)_{t \ge 0}$. Note that this formalism is totally equivalent to following the collection of characteristics $\{(x_t^i, g^i, u_p^i), i \in V_t^m \cup V_t^f\}$ (which is what is done in simulations of this process), but it will be more convenient to use when formulating the population dynamics below.

For a bounded function $f: \mathcal{C} \to \mathbb{R}$ and a counting measure $\nu \in \mathcal{M}$, we have

$$ \int_{\mathcal{C}} f(c) \nu/dc = \sum_{i} f(x^i, g^i, u_p^i), \qquad (3) $$

with the sum on the right-hand side being taken over all atoms in $\nu$. For convenience, in what follows we shall mostly use the integral form on the left-hand side of (3) when referring to such quantities. In terms of our application, we are simply summing the value of the function $f$ applied to each individual's characteristic, over all individuals in the population. For example, if the function $f$ is the indicator function of a set $E_1 \times E_2 \times E_3 \subset \mathcal{X} \times \mathcal{G} \times \mathcal{U}_p$, then $\int f(x, g, u_p)\nu(dx, dg, du_p)$ counts how many individuals in the population encoded by $\nu$ have a spatial location which is in $E_1$, a genotype in $E_2$ and a phenotype in $E_3$.

We start from an initial empirical distribution of males and females $(\nu_0^m, \nu_0^f)$ at time 0.

## Movement

We suppose that individuals move in space during their lifetimes, independently of each other and according to Brownian motion normally reflected at the boundary of the domain $\mathcal{X}$. In this work, we assume that the diffusion coefficient $\sigma^2$ is the same for all individuals. The generalisation to a trait-dependent diffusion coefficient is briefly discussed in Section 4.

## Reproduction

The reproductive dynamics of the population can be summarised as follows. Each female reproduces at a given rate, which depends on the availability (and amount) of males in a neighbourhood around her. To ease the formulation of the dynamics, the total reproduction rate of a female is split into the different rates at which it may reproduce with males of given genotypes. The genotypes of the offspring then follow Mendelian inheritance (likewise, they have probability 1/2 to be males, independently of each other), and offspring are born at the location of their mother. Their phenotypes depend on their sexes, genotypes and environment at birth (characterised by the time and spatial position of their birth).
set of characteristics of an individual as $c = (x, g, u_p)$, this variable thus belonging to

$$ \mathcal{C} := \mathcal{X} \times \mathcal{G} \times \mathcal{U}_p. \qquad (1) $$

In the case study presented in Section 3, we take $\mathcal{X} = (-1000, 1000)^2$, $\mathcal{G} = \{AA, Aa, aa\}$ (corresponding to a single biallelic locus in a diploid population), and the phenotype is a one-dimensional variable taking values in $\mathcal{U}_p = [0, 10]$ and describing the intensity of abdominal pigmentation.

Because we want to explicitly model density dependence on both sexes for reproduction, as well as sexual dimorphism, it is more convenient to consider the male and female subpopulations separately. Mathematically speaking, at any time $t \ge 0$, each of these subpopulations is represented by its empirical distribution, that is, the counting measure on $\mathcal{C}$ giving mass 1 to the triplet of characteristics $c = (x, g, u_p)$ of each individual alive at time $t$:

$$ \nu_t^m = \sum_{i \in V_t^m} \delta_{(x_t^i, g^i, u_p^i)} \quad \text{and} \quad \nu_t^f = \sum_{i \in V_t^f} \delta_{(x_t^i, g^i, u_p^i)}, \qquad (2) $$

where $V_t^m$ and $V_t^f$ are respectively the set of labels of males and females present in the population at time $t$, $\delta_c$ denotes a Dirac mass at $c$, $x_t^i$ is the position at time $t$ of individual $i$ (which will change through time, hence the dependence on $t$), and $g^i, u_p^i$ are the genotype and phenotype of individual $i$ (supposed to be fixed throughout its life). We write $\mathcal{M}$ for the set of all finite counting measures on $\mathcal{C}$. The state of the whole population at time $t$ is then described by the pair $(\nu_t^m, \nu_t^f) \in \mathcal{M}^2$, and we are interested in the dynamics of the process $(\nu_t^m, \nu_t^f)_{t \ge 0}$. Note that this formalism is totally equivalent to following the collection of characteristics $\{(x_t^i, g^i, u_p^i), i \in V_t^m \cup V_t^f\}$ (which is what is done in simulations of this process), but it will be more convenient to use when formulating the population dynamics below.

For a bounded function $f: \mathcal{C} \to \mathbb{R}$ and a counting measure $\nu \in \mathcal{M}$, we have

$$ \int_{\mathcal{C}} f(c) \nu/dc = \sum_{i} f(x^{i}, g^{i}, u_{p}^{i}), \qquad (3) $$

with the sum on the right-hand side being taken over all atoms in $\nu$. For convenience, in what follows we shall mostly use the integral form on the left-hand side of (3) when referring to such quantities. In terms of our application, we are simply summing the value of the function $f$ applied to each individual's characteristic, over all individuals in the population. For example, if the function $f$ is the indicator function of a set $E_1 \times E_2 \times E_3 \subset \mathcal{X} \times \mathcal{G} \times \mathcal{U}_p$, then $\int f(x, g, u_p)\nu(dx, dg, du_p)$ counts how many individuals in the population encoded by $\nu$ have a spatial location which is in $E_1$, a genotype in $E_2$ and a phenotype in $E_3$.

We start from an initial empirical distribution of males and females $(\nu_0^m, \nu_0^f)$ at time 0.

## Movement

We suppose that individuals move in space during their lifetimes, independently of each other and according to Brownian motion normally reflected at the boundary of the domain $\mathcal{X}$. In this work, we assume that the diffusion coefficient $\sigma^2$ is the same for all individuals. The generalisation to a trait-dependent diffusion coefficient is briefly discussed in Section 4.

## Reproduction

The reproductive dynamics of the population can be summarised as follows. Each female reproduces at a given rate, which depends on the availability (and amount) of males in a neighbourhood around her. To ease the formulation of the dynamics, the total reproduction rate of a female is split into the different rates at which it may reproduce with males of given genotypes. The genotypes of the offspring then follow Mendelian inheritance (likewise, they have probability 1/2 to be males, independently of each other), and offspring are born at the location of their mother. Their phenotypes depend on their sexes, genotypes and environment at birth (characterised by the time and spatial position of their birth).
In view of our application to abdominal pigmentation in *D. melanogaster*, we suppose that males have a fixed phenotype $u_p^m$. Female phenotypes are determined by functions of their genotype and environment at birth that take plasticity and dominance reversal effects into account. All individuals (males and females) can also die independently of each other at a given rate (corresponding to natural death), or through local density-dependent competition for resources.

Let us now introduce this dynamics in more details. Recall that $K$ is the number of possible genotypes (i.e., the cardinality of $\mathcal{G}$). For each genotype $g_l \in \mathcal{G}$, let $\lambda_{g_l} : (x, g, u_p, i, t) \mapsto \lambda_{g_l}(x, g, u_p, i, t)$ be a nonnegative function defined on $\mathcal{C} \times (\mathbb{R}_+)^2$. These functions will be used to define the instantaneous reproduction rates of females with characteristics $(x, g, u_p)$ at time $t$, and the $i$ coordinate will be used to encode the dependence of these reproduction rates on the current state of the male population. We also define a spatial interaction kernel $I : \mathbb{R}^2 \to \mathbb{R}_+$, which will encode the intensity of the interaction between two individuals at separation $z \in \mathbb{R}^2$. For instance, $I$ may be a truncated Gaussian kernel, or the indicator function of a ball of radius $R$ around 0 (i.e., individuals need to be within distance $R$ of each other to interact). For simplicity, we shall use the same kernel $I$ for births and deaths.

**Remark 2.1.** Here the interaction kernel only depends on the spatial separation between pairs of individuals. We may use a more general interaction kernel of the form $I(x - x')W((g, u_p), (g', u'_p))$, with $W : (\mathcal{G} \times \mathcal{U}_p)^2 \to \mathbb{R}_+$ describing how the intensity of interaction depends on the individuals' phenotypes and genotypes. This could allows us to incorporate mating preferences based on phenotypes, or genetic incompatibilities.

We start from the initial distributions $(\nu_0^m, \nu_0^f)$ for the male and female subpopulations, and then proceed as follows. Each female present in the population, say with characteristics $(x, g, u_p)$ at the time $t$ we consider, reproduces with a male of genotype $g_l$ at instantaneous rate

$$ \Lambda_{g_l}(x, g, u_p, \nu_t^m, t) := \lambda_{g_l} \left( x, g, u_p, \int_{\mathcal{X} \times \{g_l\} \times \mathcal{U}_p} I(x' - x) \nu_t^m(dx', dg', du'_p), t \right). \quad (4) $$

In words, the fourth coordinate of $\lambda_{g_l}$ in (4) counts the total intensity of interaction of the female located at $x$ with males of genotype $g_l$ around her.

This holds true independently for every $g_l \in \mathcal{G}$, and so the total instantaneous rate at time $t$ at which a given female with characteristics $(x, g, u_p)$ reproduces is

$$ \Lambda(x, g, u_p, \nu_t^m, t) := \sum_{l=1}^{K} \Lambda_{g_l}(x, g, u_p, \nu_t^m, t). \quad (5) $$

Females are supposed to reproduce independently of each other, and all of them produce an average number of offspring equal to $\lambda$, where $\lambda > 0$ is fixed. When the genotype of the father is $g_l$ (and $g$ is the genotype of the mother, as above), the reproduction event gives rise to a random number of offspring of different sexes and genotypes, and their phenotypes is decided in a second step. The number of male (resp., female) offspring with genotype $g_l'$ are all independent of each other and follow a Poisson distribution with parameter $(\lambda/2)p_{g_l'}(g, g_l)$, where $p_{g_l'}(g, g_l)$ is the probability that an offspring of parents with genotypes $g$ and $g_l$ has genotype $g_l'$ under Mendelian inheritance. These variables are respectively denoted by $N_{g_l'}^m$ and $N_{g_l'}^f$, and we write $N$ for the total number of offspring.

**Remark 2.2.** By standard properties of independent Poisson random variables, the total number of offspring

$$ N = \sum_{l'=1}^{K} N_{g_l'}^m + \sum_{l'=1}^{K} N_{g_l'}^f $$

follows a Poisson distribution with parameter $\lambda$. Conditionally on $N = n$, the full vector of offspring numbers $(N_1^m, \dots, N_K^m, N_1^f, \dots, N_K^f)$ follows a multinomial distribution $\text{Multinomial}(n, \pi)$, with

$$ \pi := \left( \frac{p_{g_1}(g, g_l)}{2}, \dots, \frac{p_{g_K}(g, g_l)}{2}, \frac{p_{g_1}(g, g_l)}{2}, \dots, \frac{p_{g_K}(g, g_l)}{2} \right). $$
In view of our application to abdominal pigmentation in *D. melanogaster*, we suppose that males have a fixed phenotype $u_p^m$. Female phenotypes are determined by functions of their genotype and environment at birth that take plasticity and dominance reversal effects into account. All individuals (males and females) can also die independently of each other at a given rate (corresponding to natural death), or through local density-dependent competition for resources.

Let us now introduce this dynamics in more details. Recall that $K$ is the number of possible genotypes (i.e., the cardinality of $\mathcal{G}$). For each genotype $g_l \in \mathcal{G}$, let $\lambda_{g_l} : (x, g, u_p, i, t) \mapsto \lambda_{g_l}(x, g, u_p, i, t)$ be a nonnegative function defined on $\mathcal{C} \times (\mathbb{R}_+)^2$. These functions will be used to define the instantaneous reproduction rates of females with characteristics $(x, g, u_p)$ at time $t$, and the $i$ coordinate will be used to encode the dependence of these reproduction rates on the current state of the male population. We also define a spatial interaction kernel $I: \mathbb{R}^2 \to \mathbb{R}_+$, which will encode the intensity of the interaction between two individuals at separation $z \in \mathbb{R}^2$. For instance, $I$ may be a truncated Gaussian kernel, or the indicator function of a ball of radius $R$ around 0 (i.e., individuals need to be within distance $R$ of each other to interact). For simplicity, we shall use the same kernel $I$ for births and deaths.

**Remark 2.1.** Here the interaction kernel only depends on the spatial separation between pairs of individuals. We may use a more general interaction kernel of the form $I(x - x')W((g, u_p), (g', u'_p))$, with $W: (\mathcal{G} \times \mathcal{U}_p)^2 \to \mathbb{R}_+$ describing how the intensity of interaction depends on the individuals' phenotypes and genotypes. This could allows us to incorporate mating preferences based on phenotypes, or genetic incompatibilities.

We start from the initial distributions $(\nu_0^m, \nu_0^f)$ for the male and female subpopulations, and then proceed as follows. Each female present in the population, say with characteristics $(x, g, u_p)$ at the time $t$ we consider, reproduces with a male of genotype $g_l$ at instantaneous rate

$$ \Lambda_{g_l}(x, g, u_p, \nu_t^m, t) := \lambda_{g_l} \left( x, g, u_p, \int_{\mathcal{X} \times \{g_l\} \times \mathcal{U}_p} I(x' - x) \nu_t^m(dx', dg', du'_p), t \right). \quad (4) $$

In words, the fourth coordinate of $\lambda_{g_l}$ in (4) counts the total intensity of interaction of the female located at $x$ with males of genotype $g_l$ around her.

This holds true independently for every $g_l \in \mathcal{G}$, and so the total instantaneous rate at time $t$ at which a given female with characteristics $(x, g, u_p)$ reproduces is

$$ \Lambda(x, g, u_p, \nu_t^m, t) := \sum_{l=1}^{K} \Lambda_{g_l}(x, g, u_p, \nu_t^m, t). \quad (5) $$

Females are supposed to reproduce independently of each other, and all of them produce an average number of offspring equal to $\lambda$, where $\lambda > 0$ is fixed. When the genotype of the father is $g_l$ (and $g$ is the genotype of the mother, as above), the reproduction event gives rise to a random number of offspring of different sexes and genotypes, and their phenotypes is decided in a second step. The number of male (resp., female) offspring with genotype $g_l'$ are all independent of each other and follow a Poisson distribution with parameter $(\lambda/2)p_{g_l'}(g, g_l)$, where $p_{g_l'}(g, g_l)$ is the probability that an offspring of parents with genotypes $g$ and $g_l$ has genotype $g_l'$ under Mendelian inheritance. These variables are respectively denoted by $N_{g_l'}^m$ and $N_{g_l'}^f$, and we write $N$ for the total number of offspring.

**Remark 2.2.** By standard properties of independent Poisson random variables, the total number of offspring

$$ N = \sum_{l'=1}^{K} N_{g_l'}^m + \sum_{l'=1}^{K} N_{g_l'}^f $$

follows a Poisson distribution with parameter $\lambda$. Conditionally on $N = n$, the full vector of offspring numbers $(N_1^m, \dots, N_K^m, N_1^f, \dots, N_K^f)$ follows a multinomial distribution $\text{Multinomial}(n, \pi)$, with

$$ \pi := \left( \frac{p_{g_1}(g, g_l)}{2}, \dots, \frac{p_{g_K}(g, g_l)}{2}, \frac{p_{g_1}(g, g_l)}{2}, \dots, \frac{p_{g_K}(g, g_l)}{2} \right). $$
Finally, let us allocate phenotypes to these offspring. Male offspring always have the same phenotype $u_p^m$, irrespective of their genotype and environment at birth. The phenotypes of female offspring are given by a function $\Psi_f : \mathcal{X} \times \mathcal{G} \times \mathbb{R}_+ \to \mathcal{U}_p$. That is, a female offspring born at location $x$, with genotype $g$ and at time $t$ will have phenotype $\Psi_f(x, g, t) \in \mathcal{U}_p$. Phenotypes are assumed to remain constant throughout the individual's life.

**Remark 2.3.** Here we do not explicitly model developmental stages that may take place between birth and adulthood (such as the pupal stages in D. melanogaster). Instead, we assume that newborns appear in the adult stage and can start reproducing directly after being born.

## Death

Recall the interaction kernel $I$ defined for reproduction events, and let $\mu : (x, g, u_p, i, t) \mapsto \mu(x, g, u_p, i, t)$ be a nonnegative function defined on $\mathcal{C} \times (\mathbb{R}_+)^2$. Each individual alive at time $t$ in the population, male or female and with characteristics $(x, g, u_p)$, dies at an instantaneous rate

$$D(x, g, u_p, \nu_t^m, \nu_t^f, t) := \mu(x, g, u_p, \int_C I(x' - x)(\nu_t^m(dx', dg', du_p') + \nu_t^f(dx', dg', du_p')), t). \quad (6)$$

In words, the death rate of an individual is a function of its characteristics, of the intensity of its interactions with the other individuals (males and females) around it and of the time considered.

## Simulation of the process

An algorithmic construction of such processes is standard and can be found in [9]. Due to the continuous change in positions of individuals generated by reflected Brownian motion and the dependency of the reproduction and death rates on these positions, we could not use the efficient simulation tools developed in [35] for our case study and, instead, we wrote a simulation code in Python which is available on GitHub. Time and memory constraints limited our ability to generate large numbers of replicates for each of the scenarios we explored, which explains the rather small number of realisations of the stochastic process displayed in the figures presented in the following section.

# 3 Pigmentation-related genetic and phenotypic diversity in *Drosophila melanogaster* populations

We now turn to our case study on *D. melanogaster* body pigmentation. Figure 1 shows the significant variation in pigmentation of the posterior abdomen in females of different genetic variants of *D. melanogaster*, and in particular in the A5, A6 and A7 cuticular segments. Male always display a dark pigmentation (not shown).

Because our aim is to highlight principles, and not to provide quantitative predictions of allele frequencies and phenotypic variability, we study a caricature of the interplay between the genetic and the environmental components of pigmentation. The phenotype is encoded as a number in [0, 10], with 0 corresponding to the lightest pigmentation and 10 to the darkest. Males always have pigmentation intensity 10. The environmental component of pigmentation in females is driven by the environmental temperature at birth, which is a function of the spatial location and time of the birth event denoted by $\theta(x, t)$. To simplify the analysis, we suppose that a single biallelic gene influences the pigmentation intensity, with an allele $A$ giving rise to darker phenotypes (favoured at low temperatures, as dark flies warm up slightly better) and an allele $a$ giving rise to lighter phenotypes (favoured at higher temperatures, as lightly pigmented flies then warm up slightly less). Moreover, we suppose that pigmentation in females is plastic (the same genotype possibly leading to different phenotypes) and subject to dominance reversal: a female with genotype $Aa$ will have a pigmentation similar to that of a female homozygous for the $A$ allele when the temperature at birth is low, while it will be similar to the pigmentation of a female homozygous for the $a$ allele when the temperature at birth is in the highest part of the range. The mapping relating genotype, environment and phenotype is a function $\Psi_f$ defined thereafter. Let us emphasise again that female phenotypes are decided at their birth and
Finally, let us allocate phenotypes to these offspring. Male offspring always have the same phenotype $u_p^m$, irrespective of their genotype and environment at birth. The phenotypes of female offspring are given by a function $\Psi_f : \mathcal{X} \times \mathcal{G} \times \mathbb{R}_+ \to \mathcal{U}_p$. That is, a female offspring born at location $x$, with genotype $g$ and at time $t$ will have phenotype $\Psi_f(x, g, t) \in \mathcal{U}_p$. Phenotypes are assumed to remain constant throughout the individual's life.

**Remark 2.3.** Here we do not explicitly model developmental stages that may take place between birth and adulthood (such as the pupal stages in D. melanogaster). Instead, we assume that newborns appear in the adult stage and can start reproducing directly after being born.

## Death

Recall the interaction kernel $I$ defined for reproduction events, and let $\mu : (x, g, u_p, i, t) \mapsto \mu(x, g, u_p, i, t)$ be a nonnegative function defined on $\mathcal{C} \times (\mathbb{R}_+)^2$. Each individual alive at time $t$ in the population, male or female and with characteristics $(x, g, u_p)$, dies at an instantaneous rate

$$D(x, g, u_p, \nu_t^m, \nu_t^f, t) := \mu\left(x, g, u_p, \int_C I(x' - x)(\nu_t^m(dx', dg', du_p') + \nu_t^f(dx', dg', du_p')) , t\right). \quad (6)$$

In words, the death rate of an individual is a function of its characteristics, of the intensity of its interactions with the other individuals (males and females) around it and of the time considered.

## Simulation of the process

An algorithmic construction of such processes is standard and can be found in [9]. Due to the continuous change in positions of individuals generated by reflected Brownian motion and the dependency of the reproduction and death rates on these positions, we could not use the efficient simulation tools developed in [35] for our case study and, instead, we wrote a simulation code in Python which is available on GitHub. Time and memory constraints limited our ability to generate large numbers of replicates for each of the scenarios we explored, which explains the rather small number of realisations of the stochastic process displayed in the figures presented in the following section.

# 3 Pigmentation-related genetic and phenotypic diversity in *Drosophila melanogaster* populations

We now turn to our case study on *D. melanogaster* body pigmentation. Figure 1 shows the significant variation in pigmentation of the posterior abdomen in females of different genetic variants of *D. melanogaster*, and in particular in the A5, A6 and A7 cuticular segments. Male always display a dark pigmentation (not shown).

Because our aim is to highlight principles, and not to provide quantitative predictions of allele frequencies and phenotypic variability, we study a caricature of the interplay between the genetic and the environmental components of pigmentation. The phenotype is encoded as a number in [0, 10], with 0 corresponding to the lightest pigmentation and 10 to the darkest. Males always have pigmentation intensity 10. The environmental component of pigmentation in females is driven by the environmental temperature at birth, which is a function of the spatial location and time of the birth event denoted by $\theta(x, t)$. To simplify the analysis, we suppose that a single biallelic gene influences the pigmentation intensity, with an allele A giving rise to darker phenotypes (favoured at low temperatures, as dark flies warm up slightly better) and an allele a giving rise to lighter phenotypes (favoured at higher temperatures, as lightly pigmented flies then warm up slightly less). Moreover, we suppose that pigmentation in females is plastic (the same genotype possibly leading to different phenotypes) and subject to dominance reversal: a female with genotype $Aa$ will have a pigmentation similar to that of a female homozygous for the A allele when the temperature at birth is low, while it will be similar to the pigmentation of a female homozygous for the a allele when the temperature at birth is in the highest part of the range. The mapping relating genotype, environment and phenotype is a function $\Psi_f$ defined thereafter. Let us emphasise again that female phenotypes are decided at their birth and
<figure><img src="image_3.png" /><figcaption>Figure 1: Abdominal pigmentation of female *Drosophila melanogaster* flies from the isogenic lines *Dark* and *Pale*, and F1 crosses of *Dark* and *Pale* individuals grown at different environmental temperatures (18°C, 25°C and 29°C). For all genotypes, most of the phenotypic variation occurs on the A6 and A7 segments, which are the segments that we use to calibrate our model parameters. At all developmental temperatures, the pigmentation of the *Dark* variant is darker than the pigmentation of the *Pale* variant; in all genetic variants, flies developing at lower temperatures are darker than flies developing at higher temperatures (these differences are quantified by differences in mean grey values measured on imaged cuticles, see [19]). F1 crosses display dominance reversal: their pigmentation is close to the *Dark* pigmentation at low developmental temperatures, and close to the *Pale* pigmentation at high developmental temperatures. See [19, 30] for more details.</figcaption></figure>

293 depend on the environmental temperatures at that moment; if the environment then changes, the 294 individual may become maladapted even though it was adapted in the early moments of its life.

295 **3.1 Model parameters**

296 We take $\mathcal{X} = (-1000, 1000)^2$, $\mathcal{G} = \{AA, Aa, aa\}$ (genotypes), $U_p = [0, 10]$ (phenotypes). Time is 297 considered in days, so that $T = 365$ corresponds to a complete year.

The choice of considering a single bi-allelic locus exhibiting dominance reversal for pigmentation is based on previous empirical work on the *Dark* and *Pale* isogenic lines [19], in which most of the difference in pigmentation is due to a single autosomal locus. The reaction norm of the pigmentation of the A6 segment was studied in [31], where it was found to be well-described by a polynomial of order 3. For the sake of simplicity (and because it is more easily compared with other studies based on reaction norms), we approximate this curve by a sigmoid function with an inflection point at 25°C (see Figure 3, left). To properly define the reaction norm we use, for $\alpha \ge 0$ let $\sigma_\alpha$ be the function

$$ \sigma_{\alpha}(x) = \frac{1}{1 + \exp(-\alpha x)}, \quad \text{for all } x \in \mathbb{R}. \qquad (7) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Abdominal pigmentation of female *Drosophila melanogaster* flies from the isogenic lines *Dark* and *Pale*, and F1 crosses of *Dark* and *Pale* individuals grown at different environmental temperatures (18°C, 25°C and 29°C). For all genotypes, most of the phenotypic variation occurs on the A6 and A7 segments, which are the segments that we use to calibrate our model parameters. At all developmental temperatures, the pigmentation of the *Dark* variant is darker than the pigmentation of the *Pale* variant; in all genetic variants, flies developing at lower temperatures are darker than flies developing at higher temperatures (these differences are quantified by differences in mean grey values measured on imaged cuticles, see [19]). F1 crosses display dominance reversal: their pigmentation is close to the *Dark* pigmentation at low developmental temperatures, and close to the *Pale* pigmentation at high developmental temperatures. See [19, 30] for more details.</figcaption></figure>

293 depend on the environmental temperatures at that moment; if the environment then changes, the
294 individual may become maladapted even though it was adapted in the early moments of its life.

295 **3.1 Model parameters**

296 We take $\mathcal{X} = (-1000, 1000)^2$, $\mathcal{G} = \{AA, Aa, aa\}$ (genotypes), $U_p = [0, 10]$ (phenotypes). Time is
297 considered in days, so that $T = 365$ corresponds to a complete year.

The choice of considering a single bi-allelic locus exhibiting dominance reversal for pigmentation is based on previous empirical work on the *Dark* and *Pale* isogenic lines [19], in which most of the difference in pigmentation is due to a single autosomal locus. The reaction norm of the pigmentation of the A6 segment was studied in [31], where it was found to be well-described by a polynomial of order 3. For the sake of simplicity (and because it is more easily compared with other studies based on reaction norms), we approximate this curve by a sigmoid function with an inflection point at 25°C (see Figure 3, left). To properly define the reaction norm we use, for $\alpha \ge 0$ let $\sigma_\alpha$ be the function

$$ \sigma_{\alpha}(x) = \frac{1}{1 + \exp(-\alpha x)}, \quad \text{for all } x \in \mathbb{R}. \qquad (7) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Representation of the temperature function $\theta$ over 6 years, with $C = 11$ (constant across space) and $T = 365$. Temperatures oscillate between 10°C and 32°C.</figcaption></figure>

## Environmental temperature

The temperature at any point in space and any time is given by the periodic function $\theta$ defined by

$$ \theta(x, t) = C(x) \left( \sin \left( \frac{2\pi t}{T} \right) + 1 \right) + 10, \qquad (8) $$

where $C(x)$ is the local amplitude of the fluctuations and $T$ is their period. In what follows, we do not explore scenarios where the amplitude $C(x)$ depends on the spatial location $x$, and we simply take $C$ to be a constant. Figure 2 shows an example of temperature functions. Our simulations always start by a period of increase in temperature, so that the burn-in period of the simulation is always favourable. This prevents erratic initial crashes of population sizes (followed by extinction) before the population has actually experienced a full period of environmental fluctuations.

The minimal and maximal possible temperatures at site $x$ are respectively

$$ \theta_{\min} = 10 \quad \text{and} \quad \theta_{\max}(x) = 2C(x) + 10. \qquad (9) $$

## Phenotype as a function of the genotype and the environment

All males have phenotype $u_p^m = 10$. The pigmentation of a female with genotype $g$ born at site $x$ at time $t$ (and thus at temperature $\theta(x, t)$) is given by the function $\Psi_f$ defined as follows:

$$ \Psi_f(x, g, t) = \begin{cases} 3.5 \sigma_1(26 - \theta(x, t)) + 6.5 & \text{if } g = AA, \\ 10 \sigma_{0.8}(24.8 - \theta(x, t)) & \text{if } g = Aa, \\ 6.5 \sigma_{0.6}(20 - \theta(x, t)) & \text{if } g = aa. \end{cases} \qquad (10) $$

These functions can be visualised in Figure 3 (left). Observe that females with genotypes $AA$ tend to have a dark pigmentation (i.e., phenotypic values close to 10), while females with genotypes $aa$ tend to have a light pigmentation (phenotypic values closer to 0). Heterozygous females exhibit dominance reversal.

These choices of functions are based on experimental measures obtained in previous work [19, 30]. The reaction norm for heterozygous females was inspired by [31] and calibrated with data from F1 crosses obtained in the lab (see also Figure 1). An estimation of the pigmentation intensity in two *D. melanogaster* lines, *Dark* and *Pale*, for different temperatures at development (18°C, 25°C and 29°C) was obtained in [19] (see Figure 1B in this article). Interpolating between these values to obtain a sigmoid function, and using the reaction norm for *Dark* (*resp.*, *Pale*) females as a proxy for the reaction norm corresponding to *AA* (*resp.*, *aa*) females, we obtain the functions $\Psi_f$ for homozygous females displayed in (10).
<figure><img src="image_3.png" /><figcaption>Figure 2: Representation of the temperature function $\theta$ over 6 years, with $C = 11$ (constant across space) and $T = 365$. Temperatures oscillate between 10°C and 32°C.</figcaption></figure>

## 298 Environmental temperature

The temperature at any point in space and any time is given by the periodic function $\theta$ defined by

$$ \theta(x, t) = C(x) \left( \sin \left( \frac{2\pi t}{T} \right) + 1 \right) + 10, \qquad (8) $$

where $C(x)$ is the local amplitude of the fluctuations and $T$ is their period. In what follows, we do not explore scenarios where the amplitude $C(x)$ depends on the spatial location $x$, and we simply take $C$ to be a constant. Figure 2 shows an example of temperature functions. Our simulations always start by a period of increase in temperature, so that the burn-in period of the simulation is always favourable. This prevents erratic initial crashes of population sizes (followed by extinction) before the population has actually experienced a full period of environmental fluctuations.

The minimal and maximal possible temperatures at site $x$ are respectively

$$ \theta_{\min} = 10 \quad \text{and} \quad \theta_{\max}(x) = 2C(x) + 10. \qquad (9) $$

## 305 Phenotype as a function of the genotype and the environment

All males have phenotype $u_p^m = 10$. The pigmentation of a female with genotype $g$ born at site $x$ at time $t$ (and thus at temperature $\theta(x, t)$) is given by the function $\Psi_f$ defined as follows:

$$ \Psi_f(x, g, t) = \begin{cases} 3.5 \sigma_1(26 - \theta(x, t)) + 6.5 & \text{if } g = AA, \\ 10 \sigma_{0.8}(24.8 - \theta(x, t)) & \text{if } g = Aa, \\ 6.5 \sigma_{0.6}(20 - \theta(x, t)) & \text{if } g = aa. \end{cases} \qquad (10) $$

These functions can be visualised in Figure 3 (left). Observe that females with genotypes $AA$ tend to have a dark pigmentation (i.e., phenotypic values close to 10), while females with genotypes $aa$ tend to have a light pigmentation (phenotypic values closer to 0). Heterozygous females exhibit dominance reversal.

These choices of functions are based on experimental measures obtained in previous work [19, 30]. The reaction norm for heterozygous females was inspired by [31] and calibrated with data from F1 crosses obtained in the lab (see also Figure 1). An estimation of the pigmentation intensity in two *D. melanogaster* lines, *Dark* and *Pale*, for different temperatures at development (18°C, 25°C and 29°C) was obtained in [19] (see Figure 1B in this article). Interpolating between these values to obtain a sigmoid function, and using the reaction norm for *Dark* (*resp.*, *Pale*) females as a proxy for the reaction norm corresponding to *AA* (*resp.*, *aa*) females, we obtain the functions $\Psi_f$ for homozygous females displayed in (10).
<figure><img src="image_5.png" /><figcaption>Figure 3: (Left) Pigmentation reaction norms (i.e., pigmentation intensity as a function of temperature at birth) used for each genotype of male and female *Drosophila melanogaster* flies. The light blue curve depicts male pigmentation, which is independent of temperature and set to 10. The dark blue and red curves give the reaction norms for homozygous females (AA, with darker pigmentation, in red and aa, with lighter pigmentation, in blue). The green curve corresponds to heterozygous females, which exhibit dominance reversal (Aa females born at low temperatures have a pigmentation similar to AA females, while those born at high temperatures have a pigmentation similar to aa females). (Right) Mismatch functions for a dark phenotype (blue) and a light phenotype (red), for $\alpha = \beta = 1$.</figcaption></figure>

318 **Interaction kernel**

We fix a distance $\delta > 0$, below which the intensity of pairwise interaction is constant and above which this intensity is 0 (that is, individuals interact only if they are at distance at most $\delta$ of each other). In formula:

$$I(x) = \frac{1}{\mathrm{Vol}(B(0, \delta))} \mathbf{1}_{\{\|x\| \le \delta\}}, \quad \text{for all } x \in \mathbb{R}^2, \qquad (11)$$

319 where $\mathrm{Vol}(B(0, \delta))$ is the volume of a ball of radius $\delta$.

320 **Reproduction rates**

We fix a constant $c > 0$, and we recall the notation $\theta_{\min}$ for the minimal environmental temperature at site $x$, which was introduced in (9). The instantaneous reproduction rate at time $t$ of a female with characteristics $(x, g, u_p) \in C$ with some male of genotype $g_l \in \{AA, aA, aa\}$ when the current male population is described by the measure $\nu^m$, as specified in (4), is taken to be

$$\Lambda_{g_l}(x, g, u_p, \nu^m, t) := c(\theta(x, t) - \theta_{\min}) \left( \sigma_{0.18} \left( \int_{X \times \{g_l\} \times [0,10]} I(x' - x) \nu^m(dx', dg', du'_p) \right) - 0.5 \right). \quad (12)$$

321 The reproduction rate is proportional to $\theta(x, t) - \theta_{\min}$, which decreases when the temperature decreases (modelling the reduction in available resources, and hence the slowdown in reproduction, as winter approaches) and increases as the temperature increases (modelling the increase in resources as summer approaches). The choice of $\sigma_{0.18}$ here was made so that the reproduction rate starts saturating when about 10 individuals are in the neighbourhood of the focal female.

322 **Death rates**

323 First, let us define two families of mismatch functions, Dark and Light, that will be used to quantify the mismatch between pigmentation and environmental temperature induced respectively by a dark and a light pigmentation. We shall assume that dark flies are more adapted (and therefore die at a slower rate) in a cold environment, while they are less adapted and die faster in a warm environment. The reverse holds true for lightly pigmented drosophilae. This influence of the temperature is fine-tuned by two parameters $\alpha, \beta$, which are meant to allow the impact of pigmentation to depend on temperature differently in winter and in summer. The functions Dark$_{\alpha,\beta}$ and Light$_{\alpha,\beta}$ will be used in the definition
<figure><img src="image_4.png" /><figcaption>Figure 3: (Left) Pigmentation reaction norms (i.e., pigmentation intensity as a function of temperature at birth) used for each genotype of male and female *Drosophila melanogaster* flies. The light blue curve depicts male pigmentation, which is independent of temperature and set to 10. The dark blue and red curves give the reaction norms for homozygous females (AA, with darker pigmentation, in red and aa, with lighter pigmentation, in blue). The green curve corresponds to heterozygous females, which exhibit dominance reversal (Aa females born at low temperatures have a pigmentation similar to AA females, while those born at high temperatures have a pigmentation similar to aa females). (Right) Mismatch functions for a dark phenotype (blue) and a light phenotype (red), for $\alpha = \beta = 1$.</figcaption></figure>

318 **Interaction kernel**

We fix a distance $\delta > 0$, below which the intensity of pairwise interaction is constant and above which this intensity is 0 (that is, individuals interact only if they are at distance at most $\delta$ of each other). In formula:

$$I(x) = \frac{1}{\mathrm{Vol}(B(0, \delta))} \mathbf{1}_{\{\|x\| \le \delta\}}, \quad \text{for all } x \in \mathbb{R}^2, \qquad (11)$$

319 where $\mathrm{Vol}(B(0, \delta))$ is the volume of a ball of radius $\delta$.

320 **Reproduction rates**

We fix a constant $c > 0$, and we recall the notation $\theta_{\min}$ for the minimal environmental temperature at site $x$, which was introduced in (9). The instantaneous reproduction rate at time $t$ of a female with characteristics $(x, g, u_p) \in C$ with some male of genotype $g_l \in \{AA, aA, aa\}$ when the current male population is described by the measure $\nu^m$, as specified in (4), is taken to be

$$\Lambda_{g_l}(x, g, u_p, \nu^m, t) := c(\theta(x, t) - \theta_{\min}) \left( \sigma_{0.18} \left( \int_{\mathcal{X} \times \{g_l\} \times [0,10]} I(x' - x) \nu^m(dx', dg', du'_p) \right) - 0.5 \right). \quad (12)$$

321 The reproduction rate is proportional to $\theta(x, t) - \theta_{\min}$, which decreases when the temperature decreases (modelling the reduction in available resources, and hence the slowdown in reproduction, as winter approaches) and increases as the temperature increases (modelling the increase in resources as summer approaches). The choice of $\sigma_{0.18}$ here was made so that the reproduction rate starts saturating when about 10 individuals are in the neighbourhood of the focal female.

322 **Death rates**

323 First, let us define two families of mismatch functions, Dark and Light, that will be used to quantify the mismatch between pigmentation and environmental temperature induced respectively by a dark and a light pigmentation. We shall assume that dark flies are more adapted (and therefore die at a slower rate) in a cold environment, while they are less adapted and die faster in a warm environment. The reverse holds true for lightly pigmented drosophilae. This influence of the temperature is fine-tuned by two parameters $\alpha, \beta$, which are meant to allow the impact of pigmentation to depend on temperature differently in winter and in summer. The functions Dark$_{\alpha,\beta}$ and Light$_{\alpha,\beta}$ will be used in the definition
of the death rates. Recall the notation $\theta_{\max}(x)$ for the maximal environmental temperature at location $x$ defined in (9).

Let $\alpha, \beta > 0$, and for every temperature value $\theta$, let

$$ \mathrm{Dark}_{\alpha,\beta}(\theta) = \alpha\sigma_{0.6}(\theta - 11) \mathbf{1}_{\{\theta \le 21\}} + (\beta\sigma_{0.6}(\theta - 29) + \alpha) \mathbf{1}_{\{\theta > 21\}}, \quad (13) $$

and

$$ \mathrm{Light}_{\alpha,\beta}(\theta) = (\alpha\sigma_{0.6}(11 - \theta) + \beta) \mathbf{1}_{\{\theta \le 21\}} + \alpha\sigma_{0.6}(29 - \theta) \mathbf{1}_{\{\theta > 21\}}. \quad (14) $$

Figure 3 (right) shows a graph of these functions when $\alpha = 1 = \beta$. The form of these functions is motivated by the biology of *D. melanogaster*: outside the (approximate) range [10°C,32°C], the fly survival probability is low without specific adaptation. For dark flies, the mismatch between pigmentation and temperature increases fast over [10°C,17°C] and over [25°C,32°C], as the advantage of a dark pigmentation is mainly felt when temperatures are low and its disadvantage is particularly felt when temperatures are high. The converse is true for lightly pigmented flies.

Finally, we can write down an expression for the individual death rate. To this end, we fix $\alpha, \beta > 0$ and let $d_n, d_c \ge 0$ be two constants which will be respectively used in the natural death rate, and in the death rate due to competition. The instantaneous death rate at time $t$ of an individual (male or female) with characteristics $(x, g, u_p)$, in a global population described by the pair $(\nu^m, \nu^f)$, as specified in (6), is given by

$$ D(x, g, u_p, \nu^m, \nu^f, t) := \left[ \frac{u_p}{10} \mathrm{Dark}_{\alpha,\beta}(x,t) + \left(1 - \frac{u_p}{10}\right) \mathrm{Light}_{\alpha,\beta}(x,t) \right] (\theta_{\max}(x) - \theta(x,t)) \\ \times \left( d_n + d_c \int_C I(x' - x)(\nu^m(dx', dg', du_p') + \nu^f(dx', dg', du_p')) \right). \quad (15) $$

In this formula, the first part describes the mismatch between the individual pigmentation and the current environmental temperature, obtained by interpolating between the mismatch corresponding to the darkest and the lightest phenotypes. The second part, of the form $\theta_{\max}(x) - \theta(x, t)$, decreases when the environmental temperature $\theta(x, t)$ increases (that is, flies die at a slower rate as summer approaches), and increases when $\theta(x, t)$ decreases (mimicking the global increase in death rates as environmental conditions become harsher). Following (8), this component of the death rate is periodic, with period $T$. Finally, the last part is the sum of a constant death rate $d_n$ modelling natural death, and of a competition term given by the product of the constant $d_c$ and the total density of individuals in a neighbourhood of radius $\delta$ around our focal individual.

## Simulation parameters

All the simulations of our model presented thereafter span 6 years, with a time unit of a day. Unless otherwise specified, the period of the environmental oscillations is taken to be $T = 365$ days (cf. Figure 2), all simulations starting at 21°C and with temperatures initially increasing in order to avoid early population crashes before the dynamics generated by seasonal fluctuations start taking effect. The individual diffusion coefficient is equal to a constant $m$, and the average total number of offspring of a female during a reproduction event is $\lambda$. The time step used in the simulations is written $\Delta t$.

In most of the simulations presented here, we take

$$ m = 0.001, C = 11, T = 365, \lambda = 100, \delta = 50, c = 10, \quad (16) $$

$$ \alpha = 1, \beta = 1, d_n = 10, d_c = 10, \Delta_t = 0.0001. $$

For this choice of parameters, we have $\theta_{\min} = 10$ and $\theta_{\max} = 32$.

## 3.2 Impact of temperature variations on the maintenance of genetic and phenotypic variation

Figure 4 (left) shows a typical realisation of the trajectories of the total population size and the total number of alleles of a given genotype (including males and females – the genetic composition of the
of the death rates. Recall the notation $\theta_{\max}(x)$ for the maximal environmental temperature at location $x$ defined in (9).

Let $\alpha, \beta > 0$, and for every temperature value $\theta$, let

$$ \mathrm{Dark}_{\alpha,\beta}(\theta) = \alpha\sigma_{0.6}(\theta - 11) \mathbf{1}_{\{\theta \le 21\}} + (\beta\sigma_{0.6}(\theta - 29) + \alpha) \mathbf{1}_{\{\theta > 21\}}, \quad (13) $$

and

$$ \mathrm{Light}_{\alpha,\beta}(\theta) = (\alpha\sigma_{0.6}(11 - \theta) + \beta) \mathbf{1}_{\{\theta \le 21\}} + \alpha\sigma_{0.6}(29 - \theta) \mathbf{1}_{\{\theta > 21\}}. \quad (14) $$

Figure 3 (right) shows a graph of these functions when $\alpha = 1 = \beta$. The form of these functions is motivated by the biology of *D. melanogaster*: outside the (approximate) range [10°C,32°C], the fly survival probability is low without specific adaptation. For dark flies, the mismatch between pigmentation and temperature increases fast over [10°C,17°C] and over [25°C,32°C], as the advantage of a dark pigmentation is mainly felt when temperatures are low and its disadvantage is particularly felt when temperatures are high. The converse is true for lightly pigmented flies.

Finally, we can write down an expression for the individual death rate. To this end, we fix $\alpha, \beta > 0$ and let $d_n, d_c \ge 0$ be two constants which will be respectively used in the natural death rate, and in the death rate due to competition. The instantaneous death rate at time $t$ of an individual (male or female) with characteristics $(x, g, u_p)$, in a global population described by the pair $(\nu^m, \nu^f)$, as specified in (6), is given by

$$ D(x, g, u_p, \nu^m, \nu^f, t) := \left[ \frac{u_p}{10} \mathrm{Dark}_{\alpha,\beta}(x,t) + \left(1 - \frac{u_p}{10}\right) \mathrm{Light}_{\alpha,\beta}(x,t) \right] (\theta_{\max}(x) - \theta(x,t)) \\ \times \left( d_n + d_c \int_C I(x' - x)(\nu^m(dx', dg', du_p') + \nu^f(dx', dg', du_p')) \right). \quad (15) $$

In this formula, the first part describes the mismatch between the individual pigmentation and the current environmental temperature, obtained by interpolating between the mismatch corresponding to the darkest and the lightest phenotypes. The second part, of the form $\theta_{\max}(x) - \theta(x, t)$, decreases when the environmental temperature $\theta(x, t)$ increases (that is, flies die at a slower rate as summer approaches), and increases when $\theta(x, t)$ decreases (mimicking the global increase in death rates as environmental conditions become harsher). Following (8), this component of the death rate is periodic, with period $T$. Finally, the last part is the sum of a constant death rate $d_n$ modelling natural death, and of a competition term given by the product of the constant $d_c$ and the total density of individuals in a neighbourhood of radius $\delta$ around our focal individual.

## Simulation parameters

All the simulations of our model presented thereafter span 6 years, with a time unit of a day. Unless otherwise specified, the period of the environmental oscillations is taken to be $T = 365$ days (cf. Figure 2), all simulations starting at 21°C and with temperatures initially increasing in order to avoid early population crashes before the dynamics generated by seasonal fluctuations start taking effect. The individual diffusion coefficient is equal to a constant $m$, and the average total number of offspring of a female during a reproduction event is $\lambda$. The time step used in the simulations is written $\Delta t$.

In most of the simulations presented here, we take

$$ m = 0.001, C = 11, T = 365, \lambda = 100, \delta = 50, c = 10, \quad (16) $$

$$ \alpha = 1, \beta = 1, d_n = 10, d_c = 10, \Delta t = 0.0001. $$

For this choice of parameters, we have $\theta_{\min} = 10$ and $\theta_{\max} = 32$.

## 3.2 Impact of temperature variations on the maintenance of genetic and phenotypic variation

Figure 4 (left) shows a typical realisation of the trajectories of the total population size and the total number of alleles of a given genotype (including males and females – the genetic composition of the
<figure><img src="image_6.png" /><figcaption>Figure 4: (Left) A single realisation of the stochastic population dynamics, for the set of parameter values (16). The dark blue curve shows the change in time of the total number of individuals with genotype *aa* (including males and females), the red curve corresponds to the total number of individuals with genotype *AA* and the green curve to the total number of individuals with genotype *Aa*. The change in time of the total population size is given by the light blue curve. (Right) Fifty independent realisations of the total population size process, using the same parameter values.</figcaption></figure>

subpopulations of each sex is investigated in Figure 7). In this realisation of the stochastic population process, population sizes display marked seasonal oscillations, and heterozygotes always outnumber the *AA* and *Aa* subpopulations. The latter may be expected either from a fast enough alternation of seasons favouring the maintenance of both types of homozygotes, and therefore the constant re-creation of heterozygotes; it could also be due to dominance reversal effects, as female heterozygotes tend to be born with a pigmentation close to that of adapted homozygotes, enabling them to better survive in all environments. When considering multiple realisations (see Figure 4, right, where 50 replicates of the total population size dynamics are shown), the oscillatory pattern seems to be a robust prediction, the other possible outcome being the extinction of the whole population during the first few years. The sizes of populations that survive tend to globally increase over the years; however, this may be an artefact of the fact that we had to restrict ourselves to parameter values giving rise to not too large population sizes, due to memory and time constraints. Indeed, a population whose size remains relatively small during a summer period is more at risk of extinction during the following winter, and therefore surviving several winters is more likely for populations that manage to reach larger sizes during summer (in other words, conditioning on non extinction biases the population size distribution towards larger values). Accordingly, we can observe in these simulations that extinction in a given winter period happens only to populations whose sizes are amongst the smallest during the preceding summer. We may expect this effect to be less strong when population sizes are much larger (as in natural populations of fruit flies, at least in summer).

Due to the supercritical growth of the population size as soon as temperatures are high enough in our model, the overall survival of a population hinges on its ability to survive winters. Consequently, we may expect *AA* females, with darker phenotypes, to have a selective advantage in the long term and therefore to remain in relatively high frequencies in populations surviving during several years (recall that males are supposed to display the darkest phenotype, whatever their genotypes and environments). To measure this potential long term advantage of the *AA* subpopulation (again, without splitting the analysis into males and females for now), Figure 5 plots the difference between the total number of heterozygotes and of individuals with genotype *AA* (left), and the difference between the numbers of *AA* and *aa* individuals (right) both as a function of time, using the same 50 independent realisations of the whole population process as in Figure 4. Figure 5 (left) shows that if heterozygotes are constantly more numerous in a substantial fraction of the replicates, in nearly 1/5 of them *AA* individuals actually outnumber the heterozygote subpopulation. Furthermore, while we would expect *aa* females to have an advantage over *AA* individuals in summer, Figure 5 (right) shows that one of the two homozygous subpopulations seems to take the lead during the first seasons and then remain prevalent over the other from then on. The absolute difference between the numbers of *AA* and *aa* individuals seems to be larger when *AA* individuals are more numerous, but this effect is
<figure><img src="image_5.png" /><figcaption>Figure 4: (Left) A single realisation of the stochastic population dynamics, for the set of parameter values (16). The dark blue curve shows the change in time of the total number of individuals with genotype *aa* (including males and females), the red curve corresponds to the total number of individuals with genotype *AA* and the green curve to the total number of individuals with genotype *Aa*. The change in time of the total population size is given by the light blue curve. (Right) Fifty independent realisations of the total population size process, using the same parameter values.</figcaption></figure>

subpopulations of each sex is investigated in Figure 7). In this realisation of the stochastic population process, population sizes display marked seasonal oscillations, and heterozygotes always outnumber the *AA* and *Aa* subpopulations. The latter may be expected either from a fast enough alternation of seasons favouring the maintenance of both types of homozygotes, and therefore the constant re-creation of heterozygotes; it could also be due to dominance reversal effects, as female heterozygotes tend to be born with a pigmentation close to that of adapted homozygotes, enabling them to better survive in all environments. When considering multiple realisations (see Figure 4, right, where 50 replicates of the total population size dynamics are shown), the oscillatory pattern seems to be a robust prediction, the other possible outcome being the extinction of the whole population during the first few years. The sizes of populations that survive tend to globally increase over the years; however, this may be an artefact of the fact that we had to restrict ourselves to parameter values giving rise to not too large population sizes, due to memory and time constraints. Indeed, a population whose size remains relatively small during a summer period is more at risk of extinction during the following winter, and therefore surviving several winters is more likely for populations that manage to reach larger sizes during summer (in other words, conditioning on non extinction biases the population size distribution towards larger values). Accordingly, we can observe in these simulations that extinction in a given winter period happens only to populations whose sizes are amongst the smallest during the preceding summer. We may expect this effect to be less strong when population sizes are much larger (as in natural populations of fruit flies, at least in summer).

Due to the supercritical growth of the population size as soon as temperatures are high enough in our model, the overall survival of a population hinges on its ability to survive winters. Consequently, we may expect *AA* females, with darker phenotypes, to have a selective advantage in the long term and therefore to remain in relatively high frequencies in populations surviving during several years (recall that males are supposed to display the darkest phenotype, whatever their genotypes and environments). To measure this potential long term advantage of the *AA* subpopulation (again, without splitting the analysis into males and females for now), Figure 5 plots the difference between the total number of heterozygotes and of individuals with genotype *AA* (left), and the difference between the numbers of *AA* and *aa* individuals (right) both as a function of time, using the same 50 independent realisations of the whole population process as in Figure 4. Figure 5 (left) shows that if heterozygotes are constantly more numerous in a substantial fraction of the replicates, in nearly 1/5 of them *AA* individuals actually outnumber the heterozygote subpopulation. Furthermore, while we would expect *aa* females to have an advantage over *AA* individuals in summer, Figure 5 (right) shows that one of the two homozygous subpopulations seems to take the lead during the first seasons and then remain prevalent over the other from then on. The absolute difference between the numbers of *AA* and *aa* individuals seems to be larger when *AA* individuals are more numerous, but this effect is
<figure><img src="image_7.png" /><figcaption>Figure 5: (Left) Difference between the total number of *Aa* and *AA* individuals as a function of time. (Right) Difference between the total number of *AA* and *aa* individuals as a function of time. Both graphs are based on the same 50 independent realisations of the stochastic process describing the whole population, with the same set of parameter values (16) as in Figure 4. Identical colors are used on all plots for trajectories coming from the same realisation of the process. If the population becomes extinct before the end of the simulation, the difference is set to 0 by convention.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: (Left) Frequency of the *A* allele in the total population (including males and females) as a function of time, for the same 50 realisations of the stochastic population process as in the previous figures (and identical color allocations). (Right) Frequency of the *A* allele in the male population as a function of time (same realisations of the population process as earlier). The vertical coloured lines correspond to population extinctions, after which allele frequencies are not defined.</figcaption></figure>

also associated with total population sizes being larger, at least for the realisations with the largest differences between *AA* and *aa* subpopulation sizes (see Figure 4, right, where the same colors are used for trajectories coming from the same realisation of the whole population process).

Turning to allele frequencies, the observed patterns are quite different and do not exhibit obvious seasonal fluctuations. In Figure 6 (left), 50 realisations of the frequency of allele *A* as a function of time are displayed. As suggested by the counts of different genotypes shown in Figure 5, in more than 50% of the simulations, allele *A* prevails over allele *a*. Observed *A* allele frequencies at the end of the simulations (corresponding to 6 periods of environment changes) range from 0.35 to 0.9, with a large variability between different realisations witnessing the stochastic nature of the allele frequency “trend” initiated early in the population history and then conserved through time. Figure 6 (right) shows the dynamics of the *A* allele frequency in the male subpopulation, and Figure 7 (left) shows the dynamics of the *A* allele in the female subpopulation, for the same set of realisations. No significant differences can be observed between the male and female populations in terms of allele frequency variations.

Using the same 50 realisations of the population process, we can now explore the patterns of pigmentation observed in the female subpopulation (males always have the darkest phenotype, $u_p^m = 10$). For each realisation, we compute the mean of the phenotypic values of all females alive at every given time $t$, and plot this mean pigmentation intensity as a function of time. In Figure 7
<figure><img src="image_6.png" /><figcaption>Figure 5: (Left) Difference between the total number of *Aa* and *AA* individuals as a function of time. (Right) Difference between the total number of *AA* and *aa* individuals as a function of time. Both graphs are based on the same 50 independent realisations of the stochastic process describing the whole population, with the same set of parameter values (16) as in Figure 4. Identical colors are used on all plots for trajectories coming from the same realisation of the process. If the population becomes extinct before the end of the simulation, the difference is set to 0 by convention.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: (Left) Frequency of the *A* allele in the total population (including males and females) as a function of time, for the same 50 realisations of the stochastic population process as in the previous figures (and identical color allocations). (Right) Frequency of the *A* allele in the male population as a function of time (same realisations of the population process as earlier). The vertical coloured lines correspond to population extinctions, after which allele frequencies are not defined.</figcaption></figure>

also associated with total population sizes being larger, at least for the realisations with the largest differences between *AA* and *aa* subpopulation sizes (see Figure 4, right, where the same colors are used for trajectories coming from the same realisation of the whole population process).

Turning to allele frequencies, the observed patterns are quite different and do not exhibit obvious seasonal fluctuations. In Figure 6 (left), 50 realisations of the frequency of allele *A* as a function of time are displayed. As suggested by the counts of different genotypes shown in Figure 5, in more than 50% of the simulations, allele *A* prevails over allele *a*. Observed *A* allele frequencies at the end of the simulations (corresponding to 6 periods of environment changes) range from 0.35 to 0.9, with a large variability between different realisations witnessing the stochastic nature of the allele frequency “trend” initiated early in the population history and then conserved through time. Figure 6 (right) shows the dynamics of the *A* allele frequency in the male subpopulation, and Figure 7 (left) shows the dynamics of the *A* allele in the female subpopulation, for the same set of realisations. No significant differences can be observed between the male and female populations in terms of allele frequency variations.

Using the same 50 realisations of the population process, we can now explore the patterns of pigmentation observed in the female subpopulation (males always have the darkest phenotype, $u_p^m = 10$). For each realisation, we compute the mean of the phenotypic values of all females alive at every given time $t$, and plot this mean pigmentation intensity as a function of time. In Figure 7
<figure><img src="image_9.png" /><figcaption>Figure 7: (Left) Frequency of the *A* allele in the female population as a function of time (same realisations as earlier). The vertical coloured lines correspond to population extinctions, after which allele frequencies are not defined. (Right) Mean female pigmentation as a function of time, in the same 50 realisations of the population process as in Figures 4 and 5. When a population goes extinct, its mean pigmentation remains equal to the pigmentation intensity of the last female alive (hence the horizontal coloured lines).</figcaption></figure>

(right), we see that oscillations in temperatures generate synchronised oscillations in the mean female pigmentation in all of the 50 realisations of the population process, with darker average phenotypes (i.e., higher phenotypic values) during periods of lower temperatures and *vice-versa*. The amplitude of the phenotypic oscillations does not cover the whole range of possible phenotypes which, in agreement with the variations in allele frequencies depicted in Figure 7 (left), indicates that polymorphism is maintained at the genetic and phenotypic level at any time in the female population.

The discrepancy between the oscillating mean (female) phenotype and the relatively stable allele frequencies suggests that plasticity and dominance reversal allows the population to accommodate cyclic environmental changes by a rapid phenotypic adaptation. To better understand this phenomenon, we considered several alternative scenarios of environmental variations. Simulating the population process with the same parameter values as earlier except the environmental temperature being held constant and equal to 21°C (the median of the [10°C,32°C] interval of temperatures in the previous setting) led to the extinction of the population in at most 1.5 periods of time in all of the 10 independent replicates. Warmer constant temperatures were also tested, but they caused numerical issues due to the fact that local population sizes would then take very large values that our algorithm was not able to handle. We also explored the impact of the periodicity of environmental fluctuations, using again the same parameter values (16). Figure 8 shows the fluctuations in population sizes for 4 values of the period, *T* = 30 days (i.e., a cycle is completed in the equivalent of a month), *T* = 182.5 (a period covering half a year), *T* = 365 and *T* = 547.5 (a period covering 1.5 years). Oscillations of population sizes occur in all scenarios, and their frequency matches that of environmental changes. However, faster temperature oscillations lead to a higher extinction probability, due to the many more periods of low temperature experienced by the population over the course of a simulation and a shorter amount of time for the female population to adapt by giving birth to darker individuals when winter approaches. The amplitude of population size fluctuations increases with *T*, mostly because population sizes decrease to smaller values in winter without leading to population extinction. Comparing with the extinctions caused by small winter population sizes observed in Figure 4 (right), the larger survival probability during winter observed in Figure 8(d) suggests that, in the latter case, the female population has the time to better adapt by producing darker individuals before the environmental temperature reaches its minimum. The lack of oscillations in allele frequencies observed in the same replicate populations and the fact that both the prevalence of the *A* allele and of the *a* allele lead to the same survival pattern and oscillations in population sizes (see Figure 9) highlight again the role of plasticity and dominance reversal in this adaptation.
<figure><img src="image_8.png" /><figcaption>Figure 7: (Left) Frequency of the *A* allele in the female population as a function of time (same realisations as earlier). The vertical coloured lines correspond to population extinctions, after which allele frequencies are not defined. (Right) Mean female pigmentation as a function of time, in the same 50 realisations of the population process as in Figures 4 and 5. When a population goes extinct, its mean pigmentation remains equal to the pigmentation intensity of the last female alive (hence the horizontal coloured lines).</figcaption></figure>

(right), we see that oscillations in temperatures generate synchronised oscillations in the mean female pigmentation in all of the 50 realisations of the population process, with darker average phenotypes (i.e., higher phenotypic values) during periods of lower temperatures and *vice-versa*. The amplitude of the phenotypic oscillations does not cover the whole range of possible phenotypes which, in agreement with the variations in allele frequencies depicted in Figure 7 (left), indicates that polymorphism is maintained at the genetic and phenotypic level at any time in the female population.

The discrepancy between the oscillating mean (female) phenotype and the relatively stable allele frequencies suggests that plasticity and dominance reversal allows the population to accommodate cyclic environmental changes by a rapid phenotypic adaptation. To better understand this phenomenon, we considered several alternative scenarios of environmental variations. Simulating the population process with the same parameter values as earlier except the environmental temperature being held constant and equal to 21°C (the median of the [10°C,32°C] interval of temperatures in the previous setting) led to the extinction of the population in at most 1.5 periods of time in all of the 10 independent replicates. Warmer constant temperatures were also tested, but they caused numerical issues due to the fact that local population sizes would then take very large values that our algorithm was not able to handle. We also explored the impact of the periodicity of environmental fluctuations, using again the same parameter values (16). Figure 8 shows the fluctuations in population sizes for 4 values of the period, *T* = 30 days (i.e., a cycle is completed in the equivalent of a month), *T* = 182.5 (a period covering half a year), *T* = 365 and *T* = 547.5 (a period covering 1.5 years). Oscillations of population sizes occur in all scenarios, and their frequency matches that of environmental changes. However, faster temperature oscillations lead to a higher extinction probability, due to the many more periods of low temperature experienced by the population over the course of a simulation and a shorter amount of time for the female population to adapt by giving birth to darker individuals when winter approaches. The amplitude of population size fluctuations increases with *T*, mostly because population sizes decrease to smaller values in winter without leading to population extinction. Comparing with the extinctions caused by small winter population sizes observed in Figure 4 (right), the larger survival probability during winter observed in Figure 8(d) suggests that, in the latter case, the female population has the time to better adapt by producing darker individuals before the environmental temperature reaches its minimum. The lack of oscillations in allele frequencies observed in the same replicate populations and the fact that both the prevalence of the *A* allele and of the *a* allele lead to the same survival pattern and oscillations in population sizes (see Figure 9) highlight again the role of plasticity and dominance reversal in this adaptation.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 8: Population size trajectories for different values of the period *T* of environmental fluctuations. All other parameter values are as in (16). For each value of *T*, 10 independent realisations are displayed.</figcaption></figure>

## 3.3 Separate and combined effects of sexual dimorphism, plasticity and dominance reversal

After analysing the contribution of environmental fluctuations in the maintenance and patterns of genetic and phenotypic diversity, we now investigate the contribution of sexual dimorphism, plasticity and dominance reversal.

We start with the impact of sexual dimorphism, by contrasting our previous simulations with simulations where the phenotypes of males and females depend on genotypes and environmental temperatures at birth in the same way. In Figure 10 (left), male phenotypes are no longer constant; instead, we suppose that they are determined by the same function $\Psi_f$ of genotype and temperature at birth as the function determining female phenotypes (see Equation (10)). All other parameter values are set to those stated in (16). Somewhat surprisingly, in all of the ten replicates the population goes extinct early (after at most 3 temperature cycles). By contrast, Figure 10 (right) shows that population survival (after the 6 cycles of environmental changes spanned by the simulations) is possible when the phenotypes of males and females are always equal to 10. This suggests that plasticity and dominance reversal are not sufficient to ensure the survival of the population when there is no guaranteed reserve of males adapted to low temperatures capable of helping the population to restart its growth when temperature starts to rise again. When phenotypes do not depend on genotypes in males, this reserve of more adapted individuals may also play the role of genetic reserve for the *a* allele, allowing the maintenance of both alleles during winter despite the association of the *a* allele with less adapted phenotypes in females in the case of sexual dimorphism.

To further understand the role of plasticity in shaping the allelic and phenotypic diversity in the population, we simulated our population process assuming that there was neither plasticity nor dominance reversal: all males have a phenotype value 10, all *AA* females have phenotype 9, all *Aa* females have phenotype 5 and all *aa* females have phenotype 1. Setting the darkest female phenotype to 9 allows males to still have a small advantage over the most adapted females during winter. As in
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8: Population size trajectories for different values of the period *T* of environmental fluctuations. All other parameter values are as in (16). For each value of *T*, 10 independent realisations are displayed.</figcaption></figure>

## 3.3 Separate and combined effects of sexual dimorphism, plasticity and dominance reversal

After analysing the contribution of environmental fluctuations in the maintenance and patterns of genetic and phenotypic diversity, we now investigate the contribution of sexual dimorphism, plasticity and dominance reversal.

We start with the impact of sexual dimorphism, by contrasting our previous simulations with simulations where the phenotypes of males and females depend on genotypes and environmental temperatures at birth in the same way. In Figure 10 (left), male phenotypes are no longer constant; instead, we suppose that they are determined by the same function $\Psi_f$ of genotype and temperature at birth as the function determining female phenotypes (see Equation (10)). All other parameter values are set to those stated in (16). Somewhat surprisingly, in all of the ten replicates the population goes extinct early (after at most 3 temperature cycles). By contrast, Figure 10 (right) shows that population survival (after the 6 cycles of environmental changes spanned by the simulations) is possible when the phenotypes of males and females are always equal to 10. This suggests that plasticity and dominance reversal are not sufficient to ensure the survival of the population when there is no guaranteed reserve of males adapted to low temperatures capable of helping the population to restart its growth when temperature starts to rise again. When phenotypes do not depend on genotypes in males, this reserve of more adapted individuals may also play the role of genetic reserve for the *a* allele, allowing the maintenance of both alleles during winter despite the association of the *a* allele with less adapted phenotypes in females in the case of sexual dimorphism.

To further understand the role of plasticity in shaping the allelic and phenotypic diversity in the population, we simulated our population process assuming that there was neither plasticity nor dominance reversal: all males have a phenotype value 10, all *AA* females have phenotype 9, all *Aa* females have phenotype 5 and all *aa* females have phenotype 1. Setting the darkest female phenotype to 9 allows males to still have a small advantage over the most adapted females during winter. As in
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Figure 9: Frequencies of the *A* allele as a function of time, in the same replicate populations as in Figure 8 (identical colors are used for trajectories coming from the same realisation of the process). The vertical coloured lines correspond to population extinctions, after which allele frequencies are not defined.</figcaption></figure>

Figure 6, a global trend with no clear seasonal fluctuations can be observed in *A* allele frequencies in the female population (Figure 11, left), but now mean female phenotypes also do not exhibit seasonal fluctuations (Figure 11, right). Instead, they seem to stabilise around a random value ranging from 4 and 7.5 in the simulations displayed in Figure 11 (with nontrivial fluctuations remaining around this value, due at least to the stochasticity of births and deaths).

Finally, we removed plasticity, dominance reversal and sexual dimorphism from the dynamics by assuming that all *AA* individuals (males and females) have phenotype 9, all *Aa* individuals have phenotype 5 and all *aa* individuals have phenotype 1 (independently of the environment in which they were born). Figure 12 displays ten realisations of such dynamics, with the same parameter values as in (16). In line with what was already observed in Figure 10 (left), the absence of sexual dimorphism leads to higher a higher extinction probability for the population. However, the variability in extinction times is higher, and some populations survive more periods of low temperatures than in the case where individuals can gradually adapt their phenotypes in the approach of winter. This speaks again in favour of the necessity of having a sufficiently large reserve of individuals with the best adapted (i.e., the darkest) phenotypes to survive winter. This reserve is more difficult to build in the presence of plasticity and dominance reversal, as adaptation during periods of temperature decrease is gradual and gives rise to a gradient of dark pigmentations within the population, rather than to a subpopulation of individuals all with pigmentation very close to 10. In the few surviving populations shown in Figure 12, the frequency of the *A* allele in the female population increases over time, and seems to be converging to 1. This is in keeping with the conclusions drawn from Figure 10 (right), in which setting all phenotypes to the darkest value 10 allowed the population to survive better over the few periods of environmental fluctuations considered in our simulations.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 9: Frequencies of the *A* allele as a function of time, in the same replicate populations as in Figure 8 (identical colors are used for trajectories coming from the same realisation of the process). The vertical coloured lines correspond to population extinctions, after which allele frequencies are not defined.</figcaption></figure>

Figure 6, a global trend with no clear seasonal fluctuations can be observed in *A* allele frequencies in the female population (Figure 11, left), but now mean female phenotypes also do not exhibit seasonal fluctuations (Figure 11, right). Instead, they seem to stabilise around a random value ranging from 4 and 7.5 in the simulations displayed in Figure 11 (with nontrivial fluctuations remaining around this value, due at least to the stochasticity of births and deaths).

Finally, we removed plasticity, dominance reversal and sexual dimorphism from the dynamics by assuming that all *AA* individuals (males and females) have phenotype 9, all *Aa* individuals have phenotype 5 and all *aa* individuals have phenotype 1 (independently of the environment in which they were born). Figure 12 displays ten realisations of such dynamics, with the same parameter values as in (16). In line with what was already observed in Figure 10 (left), the absence of sexual dimorphism leads to higher a higher extinction probability for the population. However, the variability in extinction times is higher, and some populations survive more periods of low temperatures than in the case where individuals can gradually adapt their phenotypes in the approach of winter. This speaks again in favour of the necessity of having a sufficiently large reserve of individuals with the best adapted (i.e., the darkest) phenotypes to survive winter. This reserve is more difficult to build in the presence of plasticity and dominance reversal, as adaptation during periods of temperature decrease is gradual and gives rise to a gradient of dark pigmentations within the population, rather than to a subpopulation of individuals all with pigmentation very close to 10. In the few surviving populations shown in Figure 12, the frequency of the *A* allele in the female population increases over time, and seems to be converging to 1. This is in keeping with the conclusions drawn from Figure 10 (right), in which setting all phenotypes to the darkest value 10 allowed the population to survive better over the few periods of environmental fluctuations considered in our simulations.
<figure><img src="image_18.png" /></figure>

Figure 10: (Left) Population size trajectories when the same function $\Psi_f$ as for female phenotypes is used to determine male phenotypes at birth. Ten independent realisations are displayed. (Right) Population size trajectories when both female and male phenotypes are constant equal to 10 (the darkest pigmentation possible). Again, ten independent realisations are displayed. There is no clear explanation for the plateau observed during each ‘spring’; it may be due to the fact that the maladaptation of dark phenotypes above 25°C transiently compensates for the increased reproduction rate due to warmer temperatures, until temperatures are sufficiently warm for reproduction to win over individual death and population growth to start again.

# 4 Discussion

We have introduced a general model for spatially structured populations in fluctuating environments, exhibiting sexual dimorphism for a phenotype which is (i) dependent on both the individual genotype and the environment state at its birth, (ii) plastic and (iii) subject to dominance reversal, and with selection acting on this phenotype. To analyse such a multi-factorial system, we have focused on the case study of abdominal pigmentation in *D. melanogaster*, in which natural selection is driven by local environmental temperatures: during summer (i.e., when temperatures are high), lightly pigmented flies are more adapted as they warm up slightly less than darker flies; on the other hand, during winter (i.e., when temperatures are low), flies with a darker pigmentation are more adapted as they warm up slightly more efficiently. Males have a constant phenotype, corresponding to the darkest pigmentation possible, while female pigmentation is variable. The explicit spatial structure of the model allows us to incorporate local competition effects that regulate population sizes in summer (when the population would grow exponentially fast without regulation), as well as the small population effect of winter during which low population densities may prevent individuals to find mates in an reasonably close neighbourhood around them, thereby adding to the already difficult temperature conditions experienced by the individuals. Temperatures (and hence selection pressures) were supposed to be the same everywhere in space and to oscillate in a seasonal manner. Simulations were carried out for a number of steps corresponding to 6 years, with the period of environmental oscillations being fixed to a year (mimicking the annual cycle of environmental temperatures) unless otherwise specified. In particular, our aim was not to investigate the *long term* behaviour of such a population, but to understand its dynamics and diversity over a few periods of environmental fluctuations.

The most striking outcome of our analysis is the fact that, in all populations which persisted over the time span of the simulations, genetic as well as phenotypic variation was maintained. When sexual dimorphism, plasticity and dominance reversal acted in combination, oscillations of the mean female pigmentation were observed but no seasonal patterns were found for the frequency of the *A* allele in the female population. This can be interpreted as the property that plasticity and dominance reversal allow the female population to adapt to variations of its environment by using the genetic diversity presently available; as a consequence, the effect of selection on allelic frequencies is limited. This seems to be at odds with the observed oscillations of allele frequencies at many loci in *Drosophila* species [5, 50]. However, these loci were not proved to be associated with sexually dimorphic phenotypes that are plastic and/or subject to dominance reversal, and therefore our predictions do not apply to them. In addition, in our simulations alleles frequencies did fluctuate in time, even though these fluctuations
<figure><img src="image_17.png" /></figure>

Figure 10: (Left) Population size trajectories when the same function $\Psi_f$ as for female phenotypes is used to determine male phenotypes at birth. Ten independent realisations are displayed. (Right) Population size trajectories when both female and male phenotypes are constant equal to 10 (the darkest pigmentation possible). Again, ten independent realisations are displayed. There is no clear explanation for the plateau observed during each ‘spring’; it may be due to the fact that the maladaptation of dark phenotypes above 25°C transiently compensates for the increased reproduction rate due to warmer temperatures, until temperatures are sufficiently warm for reproduction to win over individual death and population growth to start again.

# 4 Discussion

We have introduced a general model for spatially structured populations in fluctuating environments, exhibiting sexual dimorphism for a phenotype which is (i) dependent on both the individual genotype and the environment state at its birth, (ii) plastic and (iii) subject to dominance reversal, and with selection acting on this phenotype. To analyse such a multi-factorial system, we have focused on the case study of abdominal pigmentation in *D. melanogaster*, in which natural selection is driven by local environmental temperatures: during summer (i.e., when temperatures are high), lightly pigmented flies are more adapted as they warm up slightly less than darker flies; on the other hand, during winter (i.e., when temperatures are low), flies with a darker pigmentation are more adapted as they warm up slightly more efficiently. Males have a constant phenotype, corresponding to the darkest pigmentation possible, while female pigmentation is variable. The explicit spatial structure of the model allows us to incorporate local competition effects that regulate population sizes in summer (when the population would grow exponentially fast without regulation), as well as the small population effect of winter during which low population densities may prevent individuals to find mates in a reasonably close neighbourhood around them, thereby adding to the already difficult temperature conditions experienced by the individuals. Temperatures (and hence selection pressures) were supposed to be the same everywhere in space and to oscillate in a seasonal manner. Simulations were carried out for a number of steps corresponding to 6 years, with the period of environmental oscillations being fixed to a year (mimicking the annual cycle of environmental temperatures) unless otherwise specified. In particular, our aim was not to investigate the *long term* behaviour of such a population, but to understand its dynamics and diversity over a few periods of environmental fluctuations.

The most striking outcome of our analysis is the fact that, in all populations which persisted over the time span of the simulations, genetic as well as phenotypic variation was maintained. When sexual dimorphism, plasticity and dominance reversal acted in combination, oscillations of the mean female pigmentation were observed but no seasonal patterns were found for the frequency of the *A* allele in the female population. This can be interpreted as the property that plasticity and dominance reversal allow the female population to adapt to variations of its environment by using the genetic diversity presently available; as a consequence, the effect of selection on allelic frequencies is limited. This seems to be at odds with the observed oscillations of allele frequencies at many loci in *Drosophila* species [5, 50]. However, these loci were not proved to be associated with sexually dimorphic phenotypes that are plastic and/or subject to dominance reversal, and therefore our predictions do not apply to them. In addition, in our simulations alleles frequencies did fluctuate in time, even though these fluctuations
<figure><img src="image_19.png" /></figure>

Figure 11: (Left) Frequencies of the *A* allele in the female population as a function of time, in the absence of plasticity and dominance reversal. Parameter values are as in (16), except that all females with genotype *AA* have pigmentation intensity 9, *Aa* females have pigmentation intensity 5 and *aa* females have pigmentation intensity 1. Ten independent replicates are shown. (Right) Mean female pigmentation as a function of time, in the same replicate populations (identical color allocation in the two subfigures). Vertical coloured lines correspond to population extinctions, after which allele frequencies and mean pigmentation are not defined.

were of limited amplitude and not cyclic. Identifying small ‘Brownian’-like fluctuations vs. cyclic oscillations of small amplitude from real allele frequency data may not be an easy task.

Adaptation to the fluctuating environment occurs fast enough to generate oscillations of the mean female phenotype that are synchronised with the oscillations of the environment. This adaptive tracking appears to be all the more efficient as the period of environmental oscillations increases. However, in the absence of sexual dimorphism, it seems not to be sufficient to guarantee that winter populations are adapted enough to persist and start growing again when temperatures increase again. This suggests that the variability in phenotypes obtained *via* gradual adaptation leads to a global maladaptation of the population in winter, with too small a fraction of individuals displaying nearly optimal phenotypes. In contrast, when males have a constant phenotype which is optimal at low temperatures, as in our base scenario, the male population constitutes a reserve of adapted individuals that remain available for reproduction when temperatures rises again. Because in males phenotypes are independent of genotypes, this reserve of phenotypically adapted individuals can also shelter alleles associated with maladapted phenotypes, and act as a genetic reserve too, in line with the genomic storage theory developped for asexual populations [40]. To corroborate this explanation, simulations in which males and females have the darkest possible phenotype show a better ability to persist.

Actually, the effect of winter periods deserves a finer investigation. Indeed, it is hypothesized that in natural populations, a (potentially very small) fraction of the individuals surviving the temperature decline in fall find shelters that allows them to survive winter. Nothing is known of the order of magnitude of the size of the overwintering subpopulation, and very strong stochastic effects could be felt if this order of magnitude was very low. Moreover, natural populations are likely to live in inhomogeneous environments, with spatially variable temperature fluctuations due to differences in altitude or landscape. The number of scenarios to consider being very large (ranging from clines to highly heterogeneous and stochastic environments), we have not addressed this question here, although it could be done by simply changing the definition of the function $\theta(x, t)$. Another generalisation which can be accommodated by our model is to make the speed of spatial diffusion depend on the individual's phenotype, assuming for instance that more adapted individuals are more active and, therefore, move around faster.

Interestingly, despite the caricatural aspects of our model, we were able to identify biological relevant strategies for thermal adaptation in *Drosophila* species. Indeed, we observed several conditions allowing survival of the population in a regime of seasonal fluctuations of temperature: (i) thermal plasticity and genetic variation with sexual dimorphism, (ii) genetic variation only, with sexual dimorphism, and (iii) dark pigmentation with no plasticity, no genetic variation and no sexual dimorphism. These conditions correspond to what is observed in different species of drosophilids,
<figure><img src="image_18.png" /><figcaption>Figure 11: (Left) Frequencies of the *A* allele in the female population as a function of time, in the absence of plasticity and dominance reversal. Parameter values are as in (16), except that all females with genotype *AA* have pigmentation intensity 9, *Aa* females have pigmentation intensity 5 and *aa* females have pigmentation intensity 1. Ten independent replicates are shown. (Right) Mean female pigmentation as a function of time, in the same replicate populations (identical color allocation in the two subfigures). Vertical coloured lines correspond to population extinctions, after which allele frequencies and mean pigmentation are not defined.</figcaption></figure>

were of limited amplitude and not cyclic. Identifying small ‘Brownian’-like fluctuations vs. cyclic oscillations of small amplitude from real allele frequency data may not be an easy task.

Adaptation to the fluctuating environment occurs fast enough to generate oscillations of the mean female phenotype that are synchronised with the oscillations of the environment. This adaptive tracking appears to be all the more efficient as the period of environmental oscillations increases. However, in the absence of sexual dimorphism, it seems not to be sufficient to guarantee that winter populations are adapted enough to persist and start growing again when temperatures increase again. This suggests that the variability in phenotypes obtained *via* gradual adaptation leads to a global maladaptation of the population in winter, with too small a fraction of individuals displaying nearly optimal phenotypes. In contrast, when males have a constant phenotype which is optimal at low temperatures, as in our base scenario, the male population constitutes a reserve of adapted individuals that remain available for reproduction when temperatures rises again. Because in males phenotypes are independent of genotypes, this reserve of phenotypically adapted individuals can also shelter alleles associated with maladapted phenotypes, and act as a genetic reserve too, in line with the genomic storage theory developped for asexual populations [40]. To corroborate this explanation, simulations in which males and females have the darkest possible phenotype show a better ability to persist.

Actually, the effect of winter periods deserves a finer investigation. Indeed, it is hypothesized that in natural populations, a (potentially very small) fraction of the individuals surviving the temperature decline in fall find shelters that allows them to survive winter. Nothing is known of the order of magnitude of the size of the overwintering subpopulation, and very strong stochastic effects could be felt if this order of magnitude was very low. Moreover, natural populations are likely to live in inhomogeneous environments, with spatially variable temperature fluctuations due to differences in altitude or landscape. The number of scenarios to consider being very large (ranging from clines to highly heterogeneous and stochastic environments), we have not addressed this question here, although it could be done by simply changing the definition of the function $\theta(x, t)$. Another generalisation which can be accommodated by our model is to make the speed of spatial diffusion depend on the individual's phenotype, assuming for instance that more adapted individuals are more active and, therefore, move around faster.

Interestingly, despite the caricatural aspects of our model, we were able to identify biological relevant strategies for thermal adaptation in *Drosophila* species. Indeed, we observed several conditions allowing survival of the population in a regime of seasonal fluctuations of temperature: (i) thermal plasticity and genetic variation with sexual dimorphism, (ii) genetic variation only, with sexual dimorphism, and (iii) dark pigmentation with no plasticity, no genetic variation and no sexual dimorphism. These conditions correspond to what is observed in different species of drosophilids,
<figure><img src="image_20.png" /><figcaption>Figure 12: (Left) Total population size trajectories in the absence of plasticity, dominance reversal and sexual dimorphism. Parameter values are as in (16), except that all individuals (males and females) with genotype AA have pigmentation intensity 9, Aa individuals have pigmentation intensity 5 and aa individuals have pigmentation intensity 1. Ten independent replicates are shown. (Right) Frequency of the A allele in the female subpopulation as a function of time, in the same replicate populations (identical color allocation in the two subfigures). Vertical coloured lines correspond to population extinctions, after which allele frequencies are not defined.</figcaption></figure>

which seem to follow distinct evolutionary strategies to adapt to fluctuating thermal environment. *Drosophila melanogaster*, the species which has inspired the model, shows thermal plasticity, genetic variation and sexual dimorphism for pigmentation [3, 31]. *Drosophila kikkawai* shows reduced thermal plasticity and strong genetic variation for abdominal pigmentation located at a single locus [33]. Many species of drosophilids have a fixed dark abdominal pigmentation in both sexes, such as *Drosophila obscura* and *Drosophila saltans* [51].

**Code availability:** The Python script which was used to produce all simulations of the stochastic population model is available on GitHub: github.com/amandineveber/PigmentationModel.

**Author contributions:** The conception of the model was carried out by L.F. and A.V., with inputs from J.M.G. Simulations and analyses were carried out by L.F., under the joint supervision of J.M.G. and A.V. The manuscript was written by L.F. and A.V., with inputs from J.M.G.

**Acknowledgements:** The authors thank Luis-Miguel Chevin, Philippe Christol, Amaury Lambert, Sylvie Méléard, Michael Rera and Amir Yassin for interesting discussions and feedback. The project was funded by CNRS Mission for Transversal and Interdisciplinary Initiatives (MITI) through the PhD grant of L.F. and the funding of the collaborative research project *PigmTempAdapt*. The authors acknowledge partial support from the chaire program "Mathematical modeling and biodiversity" (Ecole Polytechnique, Muséum National d'Histoire Naturelle, Veolia Environnement, Fondation X).

**Declaration of interests:** The authors declare no conflicts of interest.

## References

[1] D.D. Aggarwal, S. Rybnikov, S. Sapielkin, E. Rashkovetsky, Z. Frenkel, M. Singh, P. Michalak, and A.B. Korol. Seasonal changes in recombination characteristics in a natural population of *Drosophila melanogaster*. *Heredity*, 127:278-287, 2021.

[2] N.H. Barton, A.M. Etheridge, and A. Véber. The infinitesimal model with dominance. *Genetics*, 225(2):iyad133, 2023.

[3] H. Bastide, A. Betancourt, V. Nolte, R. Tobler, P. Stöbe, A. Futschik, and C. Schlötterer. A genome-wide, fine-scale map of natural pigmentation variation in *Drosophila melanogaster*. *PLOS Genetics*, 9(6):1-8, 2013.
<figure><img src="image_19.png" /><figcaption>Figure 12: (Left) Total population size trajectories in the absence of plasticity, dominance reversal and sexual dimorphism. Parameter values are as in (16), except that all individuals (males and females) with genotype AA have pigmentation intensity 9, Aa individuals have pigmentation intensity 5 and aa individuals have pigmentation intensity 1. Ten independent replicates are shown. (Right) Frequency of the A allele in the female subpopulation as a function of time, in the same replicate populations (identical color allocation in the two subfigures). Vertical coloured lines correspond to population extinctions, after which allele frequencies are not defined.</figcaption></figure>

which seem to follow distinct evolutionary strategies to adapt to fluctuating thermal environment. *Drosophila melanogaster*, the species which has inspired the model, shows thermal plasticity, genetic variation and sexual dimorphism for pigmentation [3, 31]. *Drosophila kikkawai* shows reduced thermal plasticity and strong genetic variation for abdominal pigmentation located at a single locus [33]. Many species of drosophilids have a fixed dark abdominal pigmentation in both sexes, such as *Drosophila obscura* and *Drosophila saltans* [51].

**Code availability:** The Python script which was used to produce all simulations of the stochastic population model is available on GitHub: github.com/amandineveber/PigmentationModel.

**Author contributions:** The conception of the model was carried out by L.F. and A.V., with inputs from J.M.G. Simulations and analyses were carried out by L.F., under the joint supervision of J.M.G. and A.V. The manuscript was written by L.F. and A.V., with inputs from J.M.G.

**Acknowledgements:** The authors thank Luis-Miguel Chevin, Philippe Christol, Amaury Lambert, Sylvie Méléard, Michael Rera and Amir Yassin for interesting discussions and feedback. The project was funded by CNRS Mission for Transversal and Interdisciplinary Initiatives (MITI) through the PhD grant of L.F. and the funding of the collaborative research project *PigmTempAdapt*. The authors acknowledge partial support from the chaire program "Mathematical modeling and biodiversity" (Ecole Polytechnique, Muséum National d'Histoire Naturelle, Veolia Environnement, Fondation X).

**Declaration of interests:** The authors declare no conflicts of interest.

## References

[1] D.D. Aggarwal, S. Rybnikov, S. Sapielkin, E. Rashkovetsky, Z. Frenkel, M. Singh, P. Michalak, and A.B. Korol. Seasonal changes in recombination characteristics in a natural population of *Drosophila melanogaster*. *Heredity*, 127:278-287, 2021.

[2] N.H. Barton, A.M. Etheridge, and A. Véber. The infinitesimal model with dominance. *Genetics*, 225(2):iyad133, 2023.

[3] H. Bastide, A. Betancourt, V. Nolte, R. Tobler, P. Stöbe, A. Futschik, and C. Schlötterer. A genome-wide, fine-scale map of natural pigmentation variation in *Drosophila melanogaster*. *PLOS Genetics*, 9(6):1-8, 2013.
588

589

590

[4] H. Bastide, A. Yassin, E.J. Johanning, and J.E. Pool. Pigmentation in *Drosophila melanogaster* reaches its maximum in Ethiopia and correlates most strongly with ultra-violet radiation in sub-Saharan Africa. *BMC Evolutionary Biology*, 14:179, 2014.

591

592

593

[5] A.O. Bergland, E.L. Behrman, K.R. O'Brien, P.S. Schmidt, and D.A. Petrov. Genomic evidence of rapid and stable adaptive oscillations over seasonal time scales in *Drosophila*. *PLOS Genetics*, 10(11):1-19, 2014.

594

595

[6] J. Bertram and J. Masel. Different mechanisms drive the maintenance of polymorphism at loci subject to strong versus weak fluctuating selection. *Evolution*, 73(5):883-896, 2019.

596

597

[7] B. Bolker and S.W. Pacala. Using moment equations to understand stochastically driven spatial pattern formation in ecological systems. *Theoretical Population Biology*, 52(3):179-197, 1997.

598

599

[8] S. Bonamour, L.-M. Chevin, A. Charmantier, and C. Teplitsky. Phenotypic plasticity in response to climate change: the importance of cue variation. *Phil. Trans. R. Soc. B*, 374:20180178, 2019.

600

601

[9] N. Champagnat and S. Méléard. Invasion and adaptive evolution for individual-based spatially structured populations. *Journal of Mathematical Biology*, 55(2):147-188, 2007.

602

603

[10] J. Chen, V. Nolte, and C. Schlötterer. Temperature stress mediates decanalization and dominance of gene expression in *Drosophila melanogaster*. *PLoS Genet.*, 11(2):e1004883, 2015.

604

605

[11] P.L. Chesson. Coexistence of competitors in spatially and temporally varying environments: a look at the combined effects of different sorts of variability. *Theor. Popul. Biol.*, 28:263-287, 1985.

606

607

[12] L.-M. Chevin, Z. Gompert, and P. Nosil. Frequency dependence and the predictability of evolution in a changing environment. *Evolution Letters*, 6(1):21-33, 2022.

608

609

610

[13] R. Cogni, C. Kuczynski, S. Koury, E. Lavington, E.L. Behrman, K.R. O'Brien, P.S. Schmidt, and W.F. Eanes. The intensity of selection acting on the *couch potato* gene – Spatial-temporal variation in a diapause cline. *Evolution*, 68(2):538-548, 2013.

611

612

[14] P. Collet, S. Méléard, and J.A.J. Metz. A rigorous model study of the adaptive dynamics of Mendelian diploids. *Journal of Mathematical Biology*, 67:569-607, 2013.

613

614

[15] T. Connallon and S.F. Chenoweth. Dominance reversals and the maintenance of genetic variation for fitness. *PLOS Biol.*, 17:e3000118, 2019.

615

616

[16] C. Coron. Slow-fast stochastic diffusion dynamics and quasi-stationarity for diploid populations with varying size. *Journal of Mathematical Biology*, 72:171-202, 2016.

617

618

[17] C. Coron, M. Costa, H. Leman, and C. Smadi. A stochastic model for speciation by mating preferences.

619

620

621

[18] J.R. David, P. Capy, and J.-P. Gauthier. Abdominal pigmentation and growth temperature in *Drosophila melanogaster*: Similarities and differences in the norms of reaction of successive segments. *Journal of Evolutionary Biology*, 3(5-6):429-445, 1990.

622

623

624

[19] S. De Castro, F. Peronnet, J.F. Gilles, E. Mouchel-Vielh, and J.-M. Gibert. *bric à brac (bab)*, a central player in the gene regulatory network that mediates thermal plasticity of pigmentation in *Drosophila melanogaster*. *PLoS Genet.*, 14(8):e1007573, 2018.

625

626

[20] L. Dekens. Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal? *Journal of Mathematical Biology*, 84(3):15, 2022.

627

628

629

[21] L.M. Dembeck, W. Huang, M.M. Magwire, F. Lawrence, R.F. Lyman, and T.F.C. Mackay. Genetic architecture of abdominal pigmentation in *Drosophila melanogaster*. *PLoS Genetics*, 11(5):1-22, 2015.
588

589

590

[4] H. Bastide, A. Yassin, E.J. Johanning, and J.E. Pool. Pigmentation in *Drosophila melanogaster* reaches its maximum in Ethiopia and correlates most strongly with ultra-violet radiation in sub-Saharan Africa. *BMC Evolutionary Biology*, 14:179, 2014.

591

592

593

[5] A.O. Bergland, E.L. Behrman, K.R. O'Brien, P.S. Schmidt, and D.A. Petrov. Genomic evidence of rapid and stable adaptive oscillations over seasonal time scales in *Drosophila*. *PLOS Genetics*, 10(11):1-19, 2014.

594

595

[6] J. Bertram and J. Masel. Different mechanisms drive the maintenance of polymorphism at loci subject to strong versus weak fluctuating selection. *Evolution*, 73(5):883-896, 2019.

596

597

[7] B. Bolker and S.W. Pacala. Using moment equations to understand stochastically driven spatial pattern formation in ecological systems. *Theoretical Population Biology*, 52(3):179-197, 1997.

598

599

[8] S. Bonamour, L.-M. Chevin, A. Charmantier, and C. Teplitsky. Phenotypic plasticity in response to climate change: the importance of cue variation. *Phil. Trans. R. Soc. B*, 374:20180178, 2019.

600

601

[9] N. Champagnat and S. Méléard. Invasion and adaptive evolution for individual-based spatially structured populations. *Journal of Mathematical Biology*, 55(2):147-188, 2007.

602

603

[10] J. Chen, V. Nolte, and C. Schlötterer. Temperature stress mediates decanalization and dominance of gene expression in *Drosophila melanogaster*. *PLoS Genet.*, 11(2):e1004883, 2015.

604

605

[11] P.L. Chesson. Coexistence of competitors in spatially and temporally varying environments: a look at the combined effects of different sorts of variability. *Theor. Popul. Biol.*, 28:263-287, 1985.

606

607

[12] L.-M. Chevin, Z. Gompert, and P. Nosil. Frequency dependence and the predictability of evolution in a changing environment. *Evolution Letters*, 6(1):21-33, 2022.

608

609

610

[13] R. Cogni, C. Kuczynski, S. Koury, E. Lavington, E.L. Behrman, K.R. O'Brien, P.S. Schmidt, and W.F. Eanes. The intensity of selection acting on the *couch potato* gene – Spatial-temporal variation in a diapause cline. *Evolution*, 68(2):538-548, 2013.

611

612

[14] P. Collet, S. Méléard, and J.A.J. Metz. A rigorous model study of the adaptive dynamics of Mendelian diploids. *Journal of Mathematical Biology*, 67:569-607, 2013.

613

614

[15] T. Connallon and S.F. Chenoweth. Dominance reversals and the maintenance of genetic variation for fitness. *PLOS Biol.*, 17:e3000118, 2019.

615

616

[16] C. Coron. Slow-fast stochastic diffusion dynamics and quasi-stationarity for diploid populations with varying size. *Journal of Mathematical Biology*, 72:171-202, 2016.

617

618

[17] C. Coron, M. Costa, H. Leman, and C. Smadi. A stochastic model for speciation by mating preferences.

619

620

621

[18] J.R. David, P. Capy, and J.-P. Gauthier. Abdominal pigmentation and growth temperature in *Drosophila melanogaster*: Similarities and differences in the norms of reaction of successive segments. *Journal of Evolutionary Biology*, 3(5-6):429-445, 1990.

622

623

624

[19] S. De Castro, F. Peronnet, J.F. Gilles, E. Mouchel-Vielh, and J.-M. Gibert. *bric à brac (bab)*, a central player in the gene regulatory network that mediates thermal plasticity of pigmentation in *Drosophila melanogaster*. *PLoS Genet.*, 14(8):e1007573, 2018.

625

626

[20] L. Dekens. Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal? *Journal of Mathematical Biology*, 84(3):15, 2022.

627

628

629

[21] L.M. Dembeck, W. Huang, M.M. Magwire, F. Lawrence, R.F. Lyman, and T.F.C. Mackay. Genetic architecture of abdominal pigmentation in *Drosophila melanogaster*. *PLoS Genetics*, 11(5):1-22, 2015.
630

631

632

633

[22] U. Dieckmann and R. Law. Relaxation projections and the method of moments. In U. Dieckmann, R. Law, and J.A.J. Metz, editors, *The Geometry of Ecological Interactions: Simplifying Spatial Complexity*, Cambridge Studies in Adaptive Dynamics, pages 412–455. Cambridge University Press, 2000.

634

635

[23] T. Dobzhansky. Genetics of natural populations IX. Temporal changes in the composition of populations of *Drosophila pseudoobscura*. *Genetics*, 28(2):162–186, 1943.

636

637

[24] I. Dombeck and J. Jaenike. Ecological genetics of abdominal pigmentation in *Drosophila falleni*: A pleiotropic link to nematode parasitism. *Evolution*, 58(3):587–596, 2004.

638

639

[25] J.A. Endler and A.E. Houde. Geographic variation in female preferences for male traits in *Poecilia reticulata*. *Evolution*, 49(3):456–468, 1995.

640

641

[26] R.A. Fisher. The correlation between relatives on the supposition of Mendelian inheritance. *Proc. Roy. Soc. Edinburgh*, 52:399–433, 1918.

642

643

[27] T. Flatt. Life-history evolution and the genetics of fitness components in *Drosophila melanogaster*. *Genetics*, 214:3–48, 2020.

644

645

646

[28] N. Fournier and S. Méléard. A microscopic probabilistic description of a locally regulated population and macroscopic approximations. *The Annals of Applied Probability*, 14(4):1880–1919, 2004.

647

648

649

[29] L. Freoa, L.-M. Chevin, P. Christol, S. Méléard, M. Rera, A. Véber, and J.-M. Gibert. Drosophilids with darker cuticle have higher body temperature under light. *Scientific Reports*, 13(1):3513, 2023.

650

651

652

[30] J.-M. Gibert, E. Mouchel-Vielh, S. De Castro, and F. Peronnet. Phenotypic plasticity through transcriptional regulation of the evolutionary hotspot gene *tan* in *Drosophila melanogaster*. *PLOS Genetics*, 12(8):1–22, 2016.

653

654

655

[31] P. Gibert, B. Moreteau, and J.R. David. Developmental constraints on an adaptive plasticity: reaction norms of pigmentation in adult segments of *Drosophila melanogaster*. *Evolution & Development*, 2(5):249–260, 2000.

656

657

658

[32] P. Gibert, B. Moreteau, J.C. Moreteau, and J.R. David. Growth temperature and adult pigmentation in two *Drosophila* sibling species: An adaptive convergence of reaction norms in sympatric populations? *Evolution*, 50(6):2346–2353, 1996.

659

660

661

[33] P. Gibert, B. Moreteau, A. Munjal, and J.R. David. Phenotypic plasticity of abdominal pigmentation in *Drosophila kikkawai*: Multiple interactions between a major gene, sex, abdomen segment and growth temperature. *Genetica*, 105:165–176, 1999.

662

663

[34] A.R. Gibson and J.B. Falls. Thermal biology of the common garter snake *Thamnophis sirtalis* (L.): II. The effects of melanism. *Oecologia*, 43(1):99–109, 1979.

664

665

[35] D. Giorgi, S. Kaakai, and V. Lemaire. Efficient simulation of individual-based population models: the R Package IBMPopSim. *arXiv preprint arXiv:2303.06183*, 2023.

666

667

668

[36] A. Glaser-Schmitt, T.J.S Ramnarine, and J. Parsch. Rapid evolutionary change, constraints and the maintenance of polymorphism in natural populations of *Drosophila melanogaster*. *Molecular Ecology*, 00:1–14, 2023.

669

670

[37] D. Goulson. Determination of larval melanization in the moth, *Mamestra brassicae*, and the role of melanin in thermoregulation. *Heredity*, 73:471–479, 1994.

671

672

[38] K. Grieshop, E.K.H. Ho, and K.R. Kasimatis. Dominance reversals: the resolution of genetic conflict and maintenance of genetic variation. *Proc. Biol. Sci.*, 291(2018):20232816, 2024.
630

631

632

633

[22] U. Dieckmann and R. Law. Relaxation projections and the method of moments. In U. Dieckmann, R. Law, and J.A.J. Metz, editors, *The Geometry of Ecological Interactions: Simplifying Spatial Complexity*, Cambridge Studies in Adaptive Dynamics, pages 412–455. Cambridge University Press, 2000.

634

635

[23] T. Dobzhansky. Genetics of natural populations IX. Temporal changes in the composition of populations of *Drosophila pseudoobscura*. *Genetics*, 28(2):162–186, 1943.

636

637

[24] I. Dombeck and J. Jaenike. Ecological genetics of abdominal pigmentation in *Drosophila falleni*: A pleiotropic link to nematode parasitism. *Evolution*, 58(3):587–596, 2004.

638

639

[25] J.A. Endler and A.E. Houde. Geographic variation in female preferences for male traits in *Poecilia reticulata*. *Evolution*, 49(3):456–468, 1995.

640

641

[26] R.A. Fisher. The correlation between relatives on the supposition of Mendelian inheritance. *Proc. Roy. Soc. Edinburgh*, 52:399–433, 1918.

642

643

[27] T. Flatt. Life-history evolution and the genetics of fitness components in *Drosophila melanogaster*. *Genetics*, 214:3–48, 2020.

644

645

646

[28] N. Fournier and S. Méléard. A microscopic probabilistic description of a locally regulated population and macroscopic approximations. *The Annals of Applied Probability*, 14(4):1880–1919, 2004.

647

648

649

[29] L. Freoa, L.-M. Chevin, P. Christol, S. Méléard, M. Rera, A. Véber, and J.-M. Gibert. Drosophilids with darker cuticle have higher body temperature under light. *Scientific Reports*, 13(1):3513, 2023.

650

651

652

[30] J.-M. Gibert, E. Mouchel-Vielh, S. De Castro, and F. Peronnet. Phenotypic plasticity through transcriptional regulation of the evolutionary hotspot gene *tan* in *Drosophila melanogaster*. *PLOS Genetics*, 12(8):1–22, 2016.

653

654

655

[31] P. Gibert, B. Moreteau, and J.R. David. Developmental constraints on an adaptive plasticity: reaction norms of pigmentation in adult segments of *Drosophila melanogaster*. *Evolution & Development*, 2(5):249–260, 2000.

656

657

658

[32] P. Gibert, B. Moreteau, J.C. Moreteau, and J.R. David. Growth temperature and adult pigmentation in two *Drosophila* sibling species: An adaptive convergence of reaction norms in sympatric populations? *Evolution*, 50(6):2346–2353, 1996.

659

660

661

[33] P. Gibert, B. Moreteau, A. Munjal, and J.R. David. Phenotypic plasticity of abdominal pigmentation in *Drosophila kikkawai*: Multiple interactions between a major gene, sex, abdomen segment and growth temperature. *Genetica*, 105:165–176, 1999.

662

663

[34] A.R. Gibson and J.B. Falls. Thermal biology of the common garter snake *Thamnophis sirtalis* (L.): II. The effects of melanism. *Oecologia*, 43(1):99–109, 1979.

664

665

[35] D. Giorgi, S. Kaakai, and V. Lemaire. Efficient simulation of individual-based population models: the R Package IBMPopSim. *arXiv preprint arXiv:2303.06183*, 2023.

666

667

668

[36] A. Glaser-Schmitt, T.J.S Ramnarine, and J. Parsch. Rapid evolutionary change, constraints and the maintenance of polymorphism in natural populations of *Drosophila melanogaster*. *Molecular Ecology*, 00:1–14, 2023.

669

670

[37] D. Goulson. Determination of larval melanization in the moth, *Mamestra brassicae*, and the role of melanin in thermoregulation. *Heredity*, 73:471–479, 1994.

671

672

[38] K. Grieshop, E.K.H. Ho, and K.R. Kasimatis. Dominance reversals: the resolution of genetic conflict and maintenance of genetic variation. *Proc. Biol. Sci.*, 291(2018):20232816, 2024.
673

674

[39] D. Gulisija and Y. Kim. Emergence of long-term balanced polymorphism under cyclic selection of spatially variable magnitude. *Evolution*, 69(4):979-992, 2015.

675

676

[40] D. Gulisija, Y. Kim, and J.B. Plotkin. Phenotypic plasticity promotes balanced polymorphism in periodic environments by a genomic storage effect. *Genetics*, 202(4):1437-1448, 2016.

677

678

[41] G.E. Hill. Male mate choice and the evolution of female plumage coloration in the house finch. *Evolution*, 47(5):1515-1525, 1993.

679

680

[42] R.E. Irwin. The evolution of plumage dichromatism in the New World blackbirds: Social selection on female brightness. *The American Naturalist*, 144(6):890-907, 1994.

681

682

683

684

[43] M. Kapun, J.C.B. Nunez, M. Bogaerts-Márquez, J. Murga-Moreno, M. Paris, J. Outten, M. Coronado-Zamora, C. Tern, O. Rota-Stabelli, M.P. García Guerreiro, et al. Drosophila evolution over space and time (DEST): a new population genomics resource. *Molecular Biology and Evolution*, 38(12):5782-5805, 2021.

685

686

687

[44] H. B. D. Kettlewell, C. J. Cadbury, and D. R. Lees. Recessive melanism in the moth *Lasiocampa quercus* L. in industrial and non-industrial areas. In R. Creed, editor, *Ecological Genetics and Evolution*, pages 175-201. Springer US, New York, NY, 1971.

688

689

[45] J.G. Kingsolver and D.C. Wiernasz. Seasonal polyphenism in wing-melanin pattern and thermoregulatory adaptation in *Pieris* butterflies. *The American Naturalist*, 137(6):816-830, 1991.

690

691

692

[46] I.C. Kutch, H. Sevgili, T. Wittman, and K.M. Fedorka. Thermoregulatory strategy may shape immune investment in *Drosophila melanogaster*. *The Journal of Experimental Biology*, 217(20):3664-3669, 2014.

693

694

[47] E. Lafuente, D. Duneau, and P. Beldade. Genetic basis of variation in thermal developmental plasticity for *Drosophila melanogaster* body pigmentation. *Molecular Ecology*, 00:e17294, 2024.

695

696

697

[48] R.C. Lederhouse and J.M. Scriber. Intrasexual selection constrains the evolution of the dorsal color pattern of male black swallowtail butterflies, *Papilio polyxenes*. *Evolution*, 50(2):717-722, 1996.

698

699

[49] S. Lion and S. Gandon. Evolution of class-structured populations in periodic environments. *Evolution*, 76(8):1674-1688, 2022.

700

701

702

[50] H.E. Machado, A.O. Bergland, R. Taylor, S. Tilk, E. Behrman, K. Dyer, D.K. Fabian, T. Flatt, J. González, T.L. Karasov, et al. Broad geographic sampling reveals the shared basis and environmental correlates of seasonal adaptation in *Drosophila*. *eLife*, 10:e67577, 2021.

703

704

[51] T.A. Markow and P. O'Grady. *Drosophila: A guide to species identification and use*. Amsterdam Heidelberg: Academic Press, 2005.

705

706

707

[52] J.H. Massey, N. Akiyama, T. Bien, K. Dreisewerd, P.J. Wittkopp, J.Y. Yew, and A. Takahashi. Pleiotropic effects of *ebony* and *tan* on pigmentation and cuticular hydrocarbon composition in *Drosophila melanogaster*. *Frontiers in Physiology*, 10:518, 2019.

708

709

[53] S. Méléard and V. Bansaye. *Stochastic Models for Structured Populations – Scaling Limits and Long Time Behavior*. Mathematical Biosciences Institute Lecture Series. Springer Cham, 2015.

710

711

[54] R. Neukirch and A. Bovier. Survival of a recessive allele in a Mendelian diploid model. *Journal of Mathematical Biology*, 75:145-198, 2017.

712

713

[55] R. Parkash, S. Rajpurohit, and S. Ramniwas. Pigmentation and fitness trade-offs through the lens of artificial selection. *Journal of Insect Science*, 9(1):49, 2009.
673

674

[39] D. Gulisija and Y. Kim. Emergence of long-term balanced polymorphism under cyclic selection of spatially variable magnitude. *Evolution*, 69(4):979-992, 2015.

675

676

[40] D. Gulisija, Y. Kim, and J.B. Plotkin. Phenotypic plasticity promotes balanced polymorphism in periodic environments by a genomic storage effect. *Genetics*, 202(4):1437-1448, 2016.

677

678

[41] G.E. Hill. Male mate choice and the evolution of female plumage coloration in the house finch. *Evolution*, 47(5):1515-1525, 1993.

679

680

[42] R.E. Irwin. The evolution of plumage dichromatism in the New World blackbirds: Social selection on female brightness. *The American Naturalist*, 144(6):890-907, 1994.

681

682

683

684

[43] M. Kapun, J.C.B. Nunez, M. Bogaerts-Márquez, J. Murga-Moreno, M. Paris, J. Outten, M. Coronado-Zamora, C. Tern, O. Rota-Stabelli, M.P. García Guerreiro, et al. Drosophila evolution over space and time (DEST): a new population genomics resource. *Molecular Biology and Evolution*, 38(12):5782-5805, 2021.

685

686

687

[44] H. B. D. Kettlewell, C. J. Cadbury, and D. R. Lees. Recessive melanism in the moth *Lasiocampa quercus* L. in industrial and non-industrial areas. In R. Creed, editor, *Ecological Genetics and Evolution*, pages 175-201. Springer US, New York, NY, 1971.

688

689

[45] J.G. Kingsolver and D.C. Wiernasz. Seasonal polyphenism in wing-melanin pattern and thermoregulatory adaptation in *Pieris* butterflies. *The American Naturalist*, 137(6):816-830, 1991.

690

691

692

[46] I.C. Kutch, H. Sevgili, T. Wittman, and K.M. Fedorka. Thermoregulatory strategy may shape immune investment in *Drosophila melanogaster*. *The Journal of Experimental Biology*, 217(20):3664-3669, 2014.

693

694

[47] E. Lafuente, D. Duneau, and P. Beldade. Genetic basis of variation in thermal developmental plasticity for *Drosophila melanogaster* body pigmentation. *Molecular Ecology*, 00:e17294, 2024.

695

696

697

[48] R.C. Lederhouse and J.M. Scriber. Intrasexual selection constrains the evolution of the dorsal color pattern of male black swallowtail butterflies, *Papilio polyxenes*. *Evolution*, 50(2):717-722, 1996.

698

699

[49] S. Lion and S. Gandon. Evolution of class-structured populations in periodic environments. *Evolution*, 76(8):1674-1688, 2022.

700

701

702

[50] H.E. Machado, A.O. Bergland, R. Taylor, S. Tilk, E. Behrman, K. Dyer, D.K. Fabian, T. Flatt, J. González, T.L. Karasov, et al. Broad geographic sampling reveals the shared basis and environmental correlates of seasonal adaptation in *Drosophila*. *eLife*, 10:e67577, 2021.

703

704

[51] T.A. Markow and P. O'Grady. *Drosophila: A guide to species identification and use*. Amsterdam Heidelberg: Academic Press, 2005.

705

706

707

[52] J.H. Massey, N. Akiyama, T. Bien, K. Dreisewerd, P.J. Wittkopp, J.Y. Yew, and A. Takahashi. Pleiotropic effects of *ebony* and *tan* on pigmentation and cuticular hydrocarbon composition in *Drosophila melanogaster*. *Frontiers in Physiology*, 10:518, 2019.

708

709

[53] S. Méléard and V. Bansaye. *Stochastic Models for Structured Populations – Scaling Limits and Long Time Behavior*. Mathematical Biosciences Institute Lecture Series. Springer Cham, 2015.

710

711

[54] R. Neukirch and A. Bovier. Survival of a recessive allele in a Mendelian diploid model. *Journal of Mathematical Biology*, 75:145-198, 2017.

712

713

[55] R. Parkash, S. Rajpurohit, and S. Ramniwas. Pigmentation and fitness trade-offs through the lens of artificial selection. *Journal of Insect Science*, 9(1):49, 2009.
714

715

716

[56] M. Pfenninger and Q. Foucault. Population genomic time series data of a natural population suggests adaptive tracking of fluctuating environmental changes. *Integrative and Comparative Biology*, 62(6):1812-1826, 2022.

717

718

[57] M. Pigliucci. *Phenotypic Plasticity: Beyond Nature and Nurture*. Johns Hopkins University Press, 2001.

719

720

[58] S. Rajpurohit, R. Richardson, J. Dean, R. Vazquez, G. Wong, and P.S. Schmidt. Pigmentation and fitness trade-offs through the lens of artificial selection. *Biology Letters*, 12:1220160625, 2016.

721

722

[59] M.F. Rodrigues, M.D. Vibranovski, and R. Cogni. Clinal and seasonal changes are correlated in *Drosophila melanogaster* natural populations. *Evolution*, 75(8):2042-2054, 2021.

723

724

725

[60] S.M. Rudman, S.I. Greenblum, S. Rajpurohit, N.J. Betancourt, J. Hanna, S. Tilk, T. Yokoyama, D.A. Petrov, and P. Schmidt. Direct observation of adaptive tracking on ecological time scales in *Drosophila*. *Science*, 375:1246, 2022.

726

727

[61] T. Slagsvold and J.T. Lifjeld. Plumage color is a condition-dependent sexual trait in male pied flycatchers. *Evolution*, 46(3):825-828, 1992.

728

729

[62] C. Smadi. An eco-evolutionary approach of adaptation and recombination in a large population of varying size. *Stochastic Processes and their Applications*, 125(5):2054-2095, 2015.

730

731

[63] N.W. Timofeeff-Ressovsky. Zur analyse des polymorphismus bei *Adalia bipunctata*. *Biol. Zbl.*, 60:130-137, 1940.

732

733

[64] J.R. True. Insect melanism: the molecules matter. *Trends in Ecology and Evolution*, 18(12):640-647, 2003.

734

735

[65] B. Walsh and M. Lynch. *Evolution and Selection of Quantitative Traits*. Oxford University Press, 2018.

736

737

[66] A. Walter and S. Lion. Epidemiological and evolutionary consequences of periodicity in treatment coverage. *Proc. R. Soc. B*, 288:20203007, 2021.

738

739

740

[67] W.B. Watt. Adaptive significance of pigment polymorphisms in *Colias* butterflies, II. Thermoregulation and photoperiodically controlled melanin variation in *Colias eurytheme*. *Proc Natl Acad Sci USA*, 63(3):767-774, 1969.

741

742

[68] P.G. Willmer and D.M. Unwin. Field analyses of insect heat budgets: Reflectance, size and heating rates. *Oecologia*, 50(2):250-255, 1981.

743

744

745

[69] P.J. Wittkopp and P. Beldade. Development and evolution of insect pigmentation: genetic mechanisms and the potential consequences of pleiotropy. *Seminars in Cell and Developmental Biology*, 20(1):65-71, 2009.

746

747

[70] P.J. Wittkopp, K. Vaccaro, and S.B. Carroll. Evolution of *yellow* gene regulation and pigmentation in *Drosophila*. *Current Biology*, 12(18):1547-1556, 2002.

748

749

750

[71] M.J. Wittmann, A.O. Bergland, M.W. Feldman, P.S. Schmidt, and D.A. Petrov. Seasonally fluctuating selection can maintain polymorphism at many loci via segregation lift. *Proceedings of the National Academy of Sciences*, 114(46):E9932-E9941, 2017.

751

752

753

[72] S. Wright and T. Dobzhansky. Genetics of natural populations XII. Experimental reproduction of some of the changes caused by natural selection in certain populations of *Drosophila pseudoobscura*. *Genetics*, 31(2):125-156, 1946.

754

755

[73] M. Yamamichi, A.D. Letten, and S.J. Schreiber. Eco-evolutionary maintenance of diversity in fluctuating environments. *Ecol Lett.*, Suppl. 1:S152-S167, 2023.
714

715

716

[56] M. Pfenninger and Q. Foucault. Population genomic time series data of a natural population suggests adaptive tracking of fluctuating environmental changes. *Integrative and Comparative Biology*, 62(6):1812-1826, 2022.

717

718

[57] M. Pigliucci. *Phenotypic Plasticity: Beyond Nature and Nurture*. Johns Hopkins University Press, 2001.

719

720

[58] S. Rajpurohit, R. Richardson, J. Dean, R. Vazquez, G. Wong, and P.S. Schmidt. Pigmentation and fitness trade-offs through the lens of artificial selection. *Biology Letters*, 12:1220160625, 2016.

721

722

[59] M.F. Rodrigues, M.D. Vibranovski, and R. Cogni. Clinal and seasonal changes are correlated in *Drosophila melanogaster* natural populations. *Evolution*, 75(8):2042-2054, 2021.

723

724

725

[60] S.M. Rudman, S.I. Greenblum, S. Rajpurohit, N.J. Betancourt, J. Hanna, S. Tilk, T. Yokoyama, D.A. Petrov, and P. Schmidt. Direct observation of adaptive tracking on ecological time scales in *Drosophila*. *Science*, 375:1246, 2022.

726

727

[61] T. Slagsvold and J.T. Lifjeld. Plumage color is a condition-dependent sexual trait in male pied flycatchers. *Evolution*, 46(3):825-828, 1992.

728

729

[62] C. Smadi. An eco-evolutionary approach of adaptation and recombination in a large population of varying size. *Stochastic Processes and their Applications*, 125(5):2054-2095, 2015.

730

731

[63] N.W. Timofeeff-Ressovsky. Zur analyse des polymorphismus bei *Adalia bipunctata*. *Biol. Zbl.*, 60:130-137, 1940.

732

733

[64] J.R. True. Insect melanism: the molecules matter. *Trends in Ecology and Evolution*, 18(12):640-647, 2003.

734

735

[65] B. Walsh and M. Lynch. *Evolution and Selection of Quantitative Traits*. Oxford University Press, 2018.

736

737

[66] A. Walter and S. Lion. Epidemiological and evolutionary consequences of periodicity in treatment coverage. *Proc. R. Soc. B*, 288:20203007, 2021.

738

739

740

[67] W.B. Watt. Adaptive significance of pigment polymorphisms in *Colias* butterflies, II. Thermoregulation and photoperiodically controlled melanin variation in *Colias eurytheme*. *Proc Natl Acad Sci USA*, 63(3):767-774, 1969.

741

742

[68] P.G. Willmer and D.M. Unwin. Field analyses of insect heat budgets: Reflectance, size and heating rates. *Oecologia*, 50(2):250-255, 1981.

743

744

745

[69] P.J. Wittkopp and P. Beldade. Development and evolution of insect pigmentation: genetic mechanisms and the potential consequences of pleiotropy. *Seminars in Cell and Developmental Biology*, 20(1):65-71, 2009.

746

747

[70] P.J. Wittkopp, K. Vaccaro, and S.B. Carroll. Evolution of *yellow* gene regulation and pigmentation in *Drosophila*. *Current Biology*, 12(18):1547-1556, 2002.

748

749

750

[71] M.J. Wittmann, A.O. Bergland, M.W. Feldman, P.S. Schmidt, and D.A. Petrov. Seasonally fluctuating selection can maintain polymorphism at many loci via segregation lift. *Proceedings of the National Academy of Sciences*, 114(46):E9932-E9941, 2017.

751

752

753

[72] S. Wright and T. Dobzhansky. Genetics of natural populations XII. Experimental reproduction of some of the changes caused by natural selection in certain populations of *Drosophila pseudoobscura*. *Genetics*, 31(2):125-156, 1946.

754

755

[73] M. Yamamichi, A.D. Letten, and S.J. Schreiber. Eco-evolutionary maintenance of diversity in fluctuating environments. *Ecol Lett.*, Suppl. 1:S152-S167, 2023.