The effect of environmental enrichment on behavioral variability depends on genotype, behavior, and type of enrichment

Stable behavioral differences among conspecifics are seen in a wide array of species. These differences, caused (by definition) by some combination of genetic and environmental factors, are commonly referred to as individuality (Dall et al., 2004; Sih et al., 2004; Wolf and Weissing, 2010). Yet, even after experimentally homogenizing genotype and environment, individuality still persists (Gärtner, 2012), often undiminished or even increased (Ayroles et al., 2015; Buchanan et al., 2015). Multiple studies across different organisms demonstrate non-genetic individuality, which we refer to as intragenotypic variability (fruit flies: Ayroles et al., 2015; Kain et al., 2012; pea aphid: Schuett et al., 2011; nematodes: Stern et al., 2017; fish: Bierbach et al., 2017; crayfish: Vogt et al., 2008; mice: Freund et al., 2013; Kurikawa et al., 2018).

This intragenotypic variability may originate in sensitivity to stochastic microenvironmental effects that persist even when large-scale differences in environment across individuals are removed (Debat and David, 2001; Honegger and de Bivort, 2018; Willmore et al., 2007). Along these lines, environmental causes of phenotypic differences can be subdivided into deterministic (macro) and stochastic (internal or micro) aspects (Clarke, 1998; Willmore et al., 2007). Examples of macroenvironmental effects are different levels of fertilizer or different temperatures across treatments. Examples of microenvironmental effects include whether an individual animal ate more food in the morning or evening, or (in the case of flies) whether they pupated on the plastic vial or food media surface. Generally, microenvironmental effects exist within a treatment regime (Debat and David, 2001) and are hard to measure (Honegger and de Bivort, 2018). For individuals of the same genotype raised in a homogenous experimental environment, trait differences would be primarily due to microenvironmental effects, and the propensity to this variation is known as microenvironmental plasticity (Morgante et al., 2015). For a given trait, it would seem that intragenotypic variability would be maladaptive since some individuals are far from the trait optimum. Yet, in unpredictable and/or fluctuating environments, having intragenotypic variability can be advantageous (Simons, 2011). In such so-called diversifying bet-hedging strategies, variability can protect against sudden environmental changes, by increasing the likelihood that at any time a subset of the population has high fitness (Hopper, 1999; Simons, 2009).

Many studies have focused on characterizing the intragenotypic variability in morphological, physiological and behavioral traits (Abley et al., 2016; Ayroles et al., 2015; Bierbach et al., 2017; Blasco et al., 2017; Dworkin, 2005; Freund et al., 2013; Kain et al., 2012, 2015; Mellert et al., 2016; Morgante et al., 2015; Sørensen et al., 2015; Sztepanacz et al., 2017; Tonsor et al., 2013). In Drosophila melanogaster, intragenotypic variability was found in chill coma recovery time, starvation resistance, sternopleural bristles, wing traits and neuronal morphology in the larval ventral nerve cord (Dworkin, 2005; Mellert et al., 2016; Morgante et al., 2015; Sørensen et al., 2015; Sztepanacz et al., 2017). Morphological variations present in the ventral nerve cord and optic lobes are of particular note as they correlate with the timing of flight initiation and visually-guided locomotor biases, respectively, providing a link between morphological and behavioral intragenotypic variability (Linneweber et al., 2019 preprint; Mellert et al., 2016). Our research has identified intragenotypic variability in isogenic lines of D. melanogaster for turning bias, phototaxis and thermotaxis (Ayroles et al., 2015; Kain et al., 2012, 2015). Outside flies, intragenotypic variability in behavior has been studied in inbred mice and clonal fish (Amazon molly), with mice showing variation in exploratory behavior and fish showing variation in activity (Bierbach et al., 2017; Freund et al., 2013). If these examples of intragenotypic variability have their basis in microenvironmental differences, it may be hard to attribute the behavioral outcomes of specific individuals to their micro-causal underpinnings. It is, however, possible to test whether changing the degree of microenvironmental variation predicts changes in the amount of intragenotypic variability.

As most lab organisms are already raised in heavily standardized environments where microenvironmental variation is minimized, it is feasible to increase microenvironmental variation and examine the effects on behavior. ‘Enrichment’ treatments include a variety of modifications to regular laboratory housing, such as opportunities for exercise, novel object interaction and socialization (van Praag et al., 2000). Enrichment may add microenvironmental variation to a particular treatment, potentially affecting both the mean and variance of phenotypic traits (Körholz et al., 2018). Typically, enrichment treatments are hypothesized to more closely match an animal's natural habitat, increasing mean well-being and cognition (while perhaps increasing intragenotypic variability). For mice and rats, enrichment has been shown to enhance mean gliogenesis, neurogenesis, and synapse formation in the cortex, hippocampus and cerebellum, leading to improved memory and cognition (Bruel-Jungerman et al., 2005; Garthe et al., 2016; Leger et al., 2015; Mohammed et al., 2002; van Praag et al., 2000). While there are some consistent enrichment effects in rats and mice, physiological and behavioral strain differences in response to enrichment have been observed (Konkle et al., 2010; Toth et al., 2011; van de Weerd et al., 1994). There is conflicting evidence for the effect of enrichment on brain size in fish, with enrichment having no effect on three-spined sticklebacks but decreasing brain size in eastern mosquitofish (Toli et al., 2017; Turschwell and White, 2016). Early studies in D. melanogaster have found that changing the social milieu affects the size of brain structures, with social isolation leading to decreased sizes and numbers of Kenyon cell fibers (Barth and Heisenberg, 1997; Heisenberg et al., 1995; Technau, 1984). In addition, social isolation in D. melanogaster leads to faster cancer progression, suggesting that a stimulating social environment buffers against stresses (Dawson et al., 2018). In crickets, mushroom body neurogenesis is higher in enriched environments with complex visual, olfactory and auditory stimuli as compared to impoverished environments (Scotto Lomassese et al., 2000). On the other hand, enriched olfactory environments did not change mushroom body calyx size or affect odor learning in D. melanogaster (Wang et al., 2018).

The effect of enrichment on trait variability has been primarily studied in mice and rats, with researchers asking whether enriched rearing and housing conditions would decrease the statistical power to detect treatment effects by increasing within-sample variance (Toth et al., 2011). The evidence presented from behavioral and physiological studies has been conflicting: studies have shown that enrichment can increase, decrease or have no effect on variability depending on the trait in question (Freund et al., 2015; Toth et al., 2011; van de Weerd et al., 1994, 1997, 2002; Wolfer et al., 2004). A recent study by Körholz et al. (2018) chose to focus more directly on whether intragenotypic variability in behavior and brain plasticity is influenced by the diversity of experiences that results from an enriched environment. They found that enrichment increases variation in specific domains: mice from enriched environments showed higher variation in exploratory behavior (object interaction times, habituation, but not locomotion), adult neurogenesis and motor cortex thickness (Körholz et al., 2018). They attributed this increase in variation directly to the diversity of experiences or diversity of microenvironments that individuals could explore in an enriched environment.

Given the conflicting evidence for the effects of enrichment on variability, we propose two hypotheses for how enrichment may influence variation in traits. The first hypothesis is that enrichment introduces microenvironmental differences, which, in turn, increase trait differences through microenvironmental plasticity. Our second hypothesis is that by more closely matching natural conditions to which organisms are adapted, enrichment increases the robustness of development and somatic maintenance, with a corresponding reduction in variation (due to the removal of aberrant phenotypes that may appear in impoverished laboratory conditions). Even though they predict opposite outcomes, both of these hypotheses are intuitive, and have some support in the literature. We chose to test them by measuring intragenotypic variability in D. melanogaster under control and enriched treatments. This species is suitable for this work because of the ease of rearing large experimental groups from isogenic lines, and its suitability for automated behavioral phenotyping. We were also interested in testing whether the observed effects of enrichment were dependent on genotype. We measured a variety of behavioral metrics associated with spontaneous locomotion and phototaxis (Ayroles et al., 2015; Buchanan et al., 2015; Kain et al., 2012) in one isogenic line across two enrichment treatments and five isogenic lines in a single enrichment treatment to examine the effects of enrichment and the interaction of genotype and enrichment on behavioral variability. We found that while enrichment often caused a small increase in intragenotypic variability, the predominant determinants of behavioral means and variabilities were genotype and its interaction with the particular enrichment treatment.

Fly stocks were cultured in vials on Caltech formula medium at 25°C in temperature-controlled incubators on a 12 h:12 h light:dark cycle. For our isogenic populations, we used Canton-S and four lines from the Drosophila Genetic Reference Panel (DGRP; Mackay et al., 2012). These lines (numbers 45, 105, 535 and 796) were derived from different wild-caught gravid females and then inbred for 20 generations. Thus, there is significant genetic variation between the lines, but not within them. We chose to work with these particular lines because we have previously observed that they vary in intragenotypic variability in Y-maze turn bias (Ayroles et al., 2015). Flies were subjected to two enrichment treatments: mild enrichment (the addition of a small ‘jungle gym’ to each culture vial) and intense enrichment (growth of the flies in a 1 m³ cage filled with many rotting fruit substrates, plants, rocks etc.) (Fig. 1).

For mild enrichment, in each vial of medium where experimental animals were to develop, 3 female and 2 male parental flies were housed for 3–5 days. The parents were removed and the jungle gym enrichment was inserted. The enrichment object consisted of plastic tubing, pipe cleaners and fuzzy pom-poms that were hot-glued to a balsa wood applicator stick that was inserted into the medium (Fig. 1C). These were identically constructed for ∼30 vials, with the exception of the pom-pom color, which in some vials was white and in the others pink at random. F₁ experimental progeny developed in this enriched environment for around 10 days in incubators. Once they started eclosing, they were allowed to accumulate for 1–2 days, after which they were removed and mixed under cold anesthetization with other flies from the same genotype. They were then sorted into cohorts of 40 males and 40 females and placed in mildly enriched vials for 3–5 days. At this point, their behavior was measured in the Y-mazes for 2 h according to the methods in Buchanan et al. (2015) (although here we loaded anesthetized experimental animals on ice to transfer them into the Y-mazes, rather than CO₂). After the Y-maze assay, flies were anesthetized on ice and returned to their mildly enriched vials for 5–10 days at which point their phototactic preferences were measured using FlyVac according to the methods in Kain et al. (2012). This mild enrichment procedure was used for flies of all genotypes. See Fig. 2 for representations of the six behavioral metrics we acquired in these two assays, and their distributions across Canton-S flies.

For intense enrichment, we prepared a population cage using 1 m wooden dowels to make a cubic frame, with sheer white polyester drapery material as walls (Fig. 1B). A tube of this material, normally held closed by binder clips, provided access to the inside of the cage. The items shown in Fig. 1D were introduced to the cage at the time of its construction: six kinds of fly food (a variety of decomposing fruits as well as bottles of standard cornmeal medium), houseplants, stones, varied plastic objects etc. The cage was placed outside on a deck where it experienced natural fluctuations in luminance, temperature, rainfall, wind, humidity, etc. during the course of our experiment (Fig. 1E). For the experiment, a parental Canton-S population of 200 males and 200 females was placed in the cage on 5 September 2013, and removed 9 days later. F₁ flies were collected and assayed in the Y-mazes on 28 and 30 September, and 5 October 2013. Flies were recovered from the Y-mazes using cold anesthetization and stored in unenriched standard medium tubes in groups of ∼30 individuals until testing with FlyVac on 6 October 2013 (Fig. 1F). For all assays, males and females were tested in equal proportions.

For Y-maze enrichment experiments with Canton-S, we assayed 151 control flies, 203 mildly enriched flies and 206 intensely enriched flies. For Canton-S FlyVac experiments, we assayed 175 control flies, 140 mildly enriched flies and 86 intensely enriched flies. For Y-maze enrichment experiments with the DGRP lines, we assayed: DGRP 45: 166 control, 133 mildly enriched; DGRP 105: 130 control, 148 mildly enriched; DGRP 535: 113 control, 111 mildly enriched; DGRP 796: 132 control, 128 mildly enriched. For FlyVac enrichment experiments with DGRP lines, we assayed: DGRP 45: 157 control, 144 mildly enriched; DGRP 535: 122 control, 140 mildly enriched.

Behavior measures from the Y-maze assay (turn bias, number of turns, turn-direction switchiness and turn-timing clumpiness) were calculated from the vectors of turn directions and times that each fly produced in experiment. Behavior measures from FlyVac (light-choice probability and inter-choice interval) were calculated from the FlyVac data output file (Kain et al., 2012). These measures were computed and/or collected into a common data structure in MATLAB 2013a (The MathWorks, Inc., Natick, MA). With Y-maze arrays (Buchanan et al., 2015), we measured the left-vs.-right free locomotion turning bias of individual flies. With FlyVac, we measured the locomotory response to light cues of agitated flies (‘fast phototaxis’; Scott, 1943). We have previously used both of these assays to detect genetic and neural circuit regulators of intragenotypic variability (Ayroles et al., 2015; Buchanan et al., 2015; Kain et al., 2012), and between them, we examined both spontaneous and stimulus-evoked behaviors. In addition to turn bias in the Y-maze assay, we assessed: (1) the number of turns completed by individual flies within the 2 h trials, (2) flies' tendencies to alternate between left and right turns successively (‘switchiness’), and (3) the extent to which their turning events were clustered in time (‘clumpiness’). For the FlyVac dataset, we measured the average interval between phototactic choices in addition to the light-choice probability. While the Y-maze arrays and FlyVac were primarily designed to measure locomotor turning and phototaxis, respectively, they also produce precise estimates of these other individual behavioral measures, so for the purpose of this study, we do not emphasize any of these measures over the others.

Null hypothesis distributions (Fig. 2) were generated in MATLAB 2013a by resampling (with replacement) a million values for each distribution as follows. (1) For turn bias and light-choice probability: (i) all observed choice values (i.e. left versus right and light versus dark) were pooled across individuals, (ii) an individual was chosen at random from all tested, and a vector of length equal to the number of behavioral choices performed by that individual during the experiment was populated randomly with values from the pool, and (iii) the turn bias or light-choice probability for that vector was recorded. (2) For number of turns: (i) the observed inter-turn intervals (ITIs) were pooled across individuals, (ii) ITIs were chosen randomly one at a time until their cumulative sum exceeded 7,200,000 ms, the length of an experiment, and (iii) the number of turns in this sequence was recorded. The moderate discrepancy in mean between the null hypothesis distribution and experimental distribution in this analysis arises from the disproportionate number of short intervals contributed to the total pool by more active animals. (3) For switchiness, which arises from slight dependence between consecutive turns in the L–R (left-right) turn sequence (Ayroles et al., 2015): (i) we implemented a Markov chain in which the L–L (=R–R) transition probabilities yielded L–R sequences with mutual information between successive turns equal to the observed mutual information [0.018 bits, P(L–L)=0.592]; (ii) an individual was chosen at random from all tested, and a choice sequence of length equal to the number of choices performed by that individual during the experiment was generated using the Markov chain, and (iii) the switchiness of this sequence [=(#L–R+#R–L)/2×turn bias×(1–turn bias)×no. turns] was recorded. (4) For clumpiness: (i) the observed ITIs were pooled across individuals, (ii) an individual was chosen at random from all tested, and a vector of length equal to the number of choices performed by that individual during the experiment was populated randomly with values from the ITI pool, and (iii) clumpiness [=MAD(ITIs)/mean(ITIs)] was recorded, where MAD is the median absolute deviation from the median. Thus, this approach reflects sampling variation in number of turns and mean ITI as well as clumpiness. (5) For inter-choice interval (ICI): (i) the observed ICIs were pooled across individuals, (ii) an individual was chosen at random from all tested, and a vector of length equal to the number of choices performed by that individual during the experiment was populated randomly with values from the ICI pool, (iii) the mean ICI across this vector was recorded.

where y_n is the behavioral outcome of an individual n that comes from a normal distribution with parameters μ_n (mean) and σ_n (s.d.). μ and σ² are vectors specified via linear models, where a and v₀ are intercepts, X is a logical predictor matrix specifying the genotype and/or enrichment treatment for each individual, and b and v are vectors of coefficients of the linear model.

Here, y_n is the vector of number of turns made by individual n modeled by a negative binomial distribution with parameters μ_n (mean) and φ_n (dispersion). Both parameter vectors are related to the coefficients of a generalized linear model (vectors β and γ) via a log-link function. X is the experimental design matrix, as above; σ² is the variance vector calculated from the mean and dispersion parameter vectors.

where y_n is the inter-choice interval of individual n, and a_n and b_n are the shape and rate parameters of the gamma distribution, respectively, and the rest of the parameters are as above.

For the estimation of all posterior distributions, we set our priors on the coefficients to broad Cauchy distributions centered at 0 to allow them to be weakly informative (Gelman et al., 2008). Our qualitative findings were robust to the choice of prior, as uninformative uniform told the same story (Fig. S1). To sample posteriors, we used four chains and 50,000–100,000 iterations per chain, with the target average proposal acceptance probability of 0.8–0.9 and maximum tree depth of 10–15, to generate a posterior distribution of 100,000–200,000 samples (50% of chain iterations were used for tuning the Hamiltonian Monte Carlo sampler parameters and were discarded as the burn-in period). The ratio of chain effective sample size to sample size was in the range of 0.7–1.2, indicating that posterior estimate error due to autocorrelation was minimal. To get posterior distributions for the coefficients of variation, we took the square root of the variance and divided it by the mean at each step in the chain. To check our model fits, we carried out graphical posterior predictive checks (Fig. S2), and found that our models fit the data well.

We adapted the methodology used by Kruschke's BEST method (Kruschke, 2013) to our own posterior distributions in order to determine which effects were inferred to differ from zero. To estimate the posterior distribution of a treatment effect, we subtracted the parameter values of one treatment condition from the control (or the other treatment condition) at each step in the chain and took the distribution of that difference. We calculated the 99% highest density interval (the credible interval) of the posterior of treatment effects to evaluate whether the treatments had an effect: if the 99% highest density interval excluded 0, we inferred an effect between the treatments. Since this approach is subject to multiple comparisons concerns, we chose the 99% credible interval (rather than e.g. 95%) as a more stringent indicator of effects. While we think that the 99% credible interval is a useful guide to pulling out the strongest effects we observe, we believe that the strength of Bayesian inference lies in being able to examine the posterior distributions as a whole and observing their trends (rather than applying a threshold to identify effects).

where CV_ij and μ_ij are the variability and mean, respectively, of genotype i in treatment group j. CV₀ and μ₀ are the grand variability and mean, averaged over all genotypes and treatments. G_i is the deviation of the variability or mean of genotype i from the grand parameter in question, calculated over all treatment groups. E_j is the deviation of treatment group j from the grand parameter, calculated over all genotypes. The treatment groups in this experiment were mild enrichment or control vials. G_i×E_j is the specific deviation of genotype i in treatment group j after accounting for the main effects of genotype i and treatment group j. All deviations were standardized by dividing them by the grand parameter value in order to interpret them as effect sizes.

To measure intragenotypic variability, we employed two automated assays (which measure spontaneous locomotion and locomotor responses to light) to rapidly collect many behavioral observations from many individual flies. We first confirmed that intragenotypic variability was present in a standard lab wild-type strain, Canton-S, in left–right turn bias and light-choice probability (Fig. 2). Indeed, the observed distributions of these measures were significantly broader (P<0.001 by bootstrapping, χ² and Kolmogorov–Smirnov tests) than expected under null models in which all individuals behaved identically, i.e. sequences of behavior drawn from identical distributions (see Materials and Methods).

We next asked if there was evidence of intragenotypic variability in other measurements taken while measuring turn bias and phototactic preference in these assays. As with turn bias, the observed distributions of number of turns, turn switchiness and turn clumpiness were significantly broader than expected under null models in which all flies behaved identically (Fig. 2A). Using FlyVac data, we observed that the distribution across flies of the average interval between phototactic choices also was broader than expected if all flies were behaving identically (Fig. 2B). Thus, intragenotypic variability was evident in all six behavioral traits examined.

We developed enrichment protocols that were either ‘mild’ or ‘intense’ (Fig. 1). Our mild enrichment treatment, the fly jungle gym, was designed to provide a variety of textures, colors and light conditions, whereas the intense cage enrichment was designed for flies to experience natural weather conditions, in addition to several different foods and a greater variety of biotic and abiotic substrates.

To confirm that the enrichment treatments had an effect on our flies, we examined the mean values of our six behavioral phenotypes under each treatment. We used a Bayesian framework with a weakly informative prior to estimate the posterior distributions of the means and variances of each behavioral metric under the mild and intense enrichment treatments. We used the 99% highest density interval, also termed the credible interval (Kruschke, 2013; see Materials and Methods) to assess whether the posterior distributions of the means of each behavioral metric were different from each other. Intense enrichment caused strong decreases in the mean of number of turns and inter-choice interval, and a strong increase in the turn switchiness when compared with mild enrichment and the control (Fig. 3A). For these behaviors, intense enrichment had a larger effect on the mean than mild enrichment. Mild enrichment had a less pronounced effect on the mean number of turns and turn clumpiness. There was no apparent effect of enrichment on turn bias and light-choice probability, although the FlyVac assay has lower power than the Y-maze assay. We viewed these observed mean changes as a positive control that the flies were sensitive to our enrichment treatments. Our results were supported by a non-parametric test of mean differences (Fig. S3A).

We estimated the effects of genotype, mild enrichment and genotype×mild enrichment on behavioral means using four DGRP lines (45, 105, 535, 796) and Canton-S (Fig. 4). Genotype had an effect (i.e. zero was not in the 99% credible interval of the posterior distribution) on all behaviors except turn bias. Mild enrichment caused a genotype-independent increase in number of turns and switchiness, but had no effect on any of the other behaviors. There were genotype×mild enrichment effects on number of turns, switchiness and clumpiness, and the direction of those effects was variable.

We examined the effect of mild and intense enrichment on intragenotypic variability in our behavioral measures (Fig. 3B). We chose to look at the coefficient of variation as our measure of intragenotypic variability in order to standardize it across multiple types of measures and control for mean effects (estimates of the posterior distributions of variance effects, not normalized by the treatment means, are included in Figs S4, S5). For nearly all behaviors, intense enrichment had a larger effect on variability than mild enrichment, but these effects were not all in the same direction (results of a non-parametric test of differences in variance are shown in Fig. S3B). Intense enrichment decreased the variability of turn bias and turn direction switchiness but increased the variability of number of turns and inter-choice interval, when compared with the control and mild enrichment treatments. Intense enrichment increased variability in clumpiness, although the effect was more pronounced upon comparison to mild enrichment as opposed to the control.

Mild enrichment caused small or no differences (zero effect was within the 99% credible interval) when compared with the control treatment for all the behavioral measures from both assays, with turn-direction switchiness and turn-timing clumpiness the most likely behavioral measures to be affected by mild enrichment. Variability in turn bias and number of turns probably increased slightly under mild enrichment, whereas clumpiness and switchiness probably decreased slightly. In two of these cases, the direction of the effect matched the direction of the intense enrichment effect; in the other two cases, it did not. To summarize, intense enrichment had stronger effects on variability than mild enrichment and the direction of these effects was behavior dependent.

In our analysis of five isogenic lines (four DGRP and Canton-S) under unenriched and mildly enriched conditions, we found that genotype, mild enrichment and genotype×mild enrichment all had effects on intragenotypic variability (Fig. 5). We found that the variability of all behavioral measurements, except the number of turns, were affected by genotype. The variability of number of turns, clumpiness, and light choice increased in a genotype-independent manner under mild enrichment. Variability of switchiness and inter-choice interval were probably also increased in a genotype-independent manner by mild enrichment (a large majority of their respective posterior distributions was above zero). We observed genotype×mild enrichment effects for number of turns, switchiness, clumpiness and light-choice probability. Of all the behavioral measures, switchiness showed the most variable and the strongest genotype×enrichment effects. To summarize, mild enrichment often increased variability in a genotype-independent fashion, but there were also frequently genotype×mild enrichment effects.

Interestingly, we found that the average magnitudes of mean effects were smaller than the average magnitudes of variability effects (Fig. 6A,B). For both mean and variability, genotype effects tended to be larger than the mild enrichment or genotype×mild enrichment effects. This pattern was especially prominent for mean effects. We also found that the size of the effects on behavioral means was uncorrelated with the size of the effects on variability (Fig. 6C).

The goal of our study was to test two opposing hypotheses about the effects of enrichment on intragenotypic behavioral variability. We hypothesized that enrichment could increase variability due to the increase in microenvironmental diversity or decrease variability due to enriched environments more closely mimicking natural conditions (resulting in more robust development of behaviors and elimination of extremes that can occur in impoverished conditions) (Körholz et al., 2018). To this end, we examined six behaviors in D. melanogaster across several genotypes and employed two levels of enrichment. We found that for five of the six behaviors, when examined across several genotypes, mild enrichment via a fly jungle gym probably led to an increase in intragenotypic variability, supporting the hypothesis that enrichment causes an increase in behavioral variability as a result of the increase in microenvironmental diversity. However, these genotype-independent effects were generally smaller than the effects of genotype or genotype×enrichment interactions. Therefore, the effect of enrichment on intragenotypic variability appears to depend on the particular genotype and behavior being assayed. When we examined the effects of enrichment on variability within a genotype, we found support for both hypotheses depending on the behavioral measure and enrichment treatment examined (mild or intense). Therefore, while it was broadly true that mild enrichment caused a small increase in the intragenotypic variability, a more granular look at the effects of enrichment revealed both strong increases and decreases in variability. We also found that the effects of enrichment on behavioral means and variabilities were largely independent of each other, with variability effects having larger magnitudes than mean effects (Fig. 6). This finding confirms that our enrichment paradigm was able to affect both mean and variability, and that these effects are potentially independent. From our experiments, it remains uncertain which aspects of enrichment influence mean and which influence variability, or indeed, if these aspects are one and the same. Behavioral variability being more strongly affected by enrichment may underscore a biological flexibility that is not present in determining mean behavior.

With respect to both mean and variability, we found that genotype usually had a larger effect than the mild enrichment. This was especially obvious when looking at the genotype effects on behavioral means, where all behaviors except turn bias showed large genotype effects (Figs 4 and 6). The lack of effect of genotype on turn bias is consistent with previous work that found no differences in the mean turn biases of 159 DGRP lines (Ayroles et al., 2015). Genotype also had strong effects on intragenotypic variability (Fig. 5), as expected (Ayroles et al., 2015; Buchanan et al., 2015; Kain et al., 2012).

We found evidence of interactions between genotype and enrichment for practically all the behavioral measures examined, although the magnitude of these interactions was behavior dependent (Fig. 6). For example, turn bias and inter-choice interval showed very little genotype-by-enrichment effect for variability, but large effects were seen for turn switchiness. Dependence of variability on the particular parameter measured was previously noted in mouse enrichment studies (Körholz et al., 2018; van de Weerd et al., 2002). Behavioral parameters may fall into different categories with respect to their response to enrichment. For example, switchiness is a measure of intra-individual variability (Fig. 2), hinting at a link between the biological mechanism controlling variability from trial to trial and individual to individual (Stamps et al., 2013). Our results also make it clear that in assessing the effects of enrichment on a particular measure of behavior, genotype cannot be ignored. These interactions are consistent with previous findings in rats and mice (Konkle et al., 2010; Toth et al., 2011; van de Weerd et al., 1994), where the effects of enrichment differed between strains. Our results also imply that there is genotype-dependent plasticity in variability. In essence, phenotype variability is not a static feature of a genotype, but depending on the trait measured, the environment can have a large effect. Evolution of plasticity has usually been examined in trait means. Our findings suggest such inquiries should extend to variability. In the future, it may be necessary to consider variability as a flexible, evolvable trait to understand how phenotype distributions arise.

Our fly jungle gym enrichment (mild enrichment) featured an array of perching sites and materials. The goal was that flies on the jungle gyms would experience a diversity of perching sites and textures, as well as be forced to navigate a more complex environment. Flies would therefore be subject to a diversity of experiences closer to what they might have in the wild while still under the constraints of laboratory conditions. We expected that the mild enrichment would mostly affect locomotion and activity behaviors, such as turn bias and inter-choice interval. We found that mild enrichment caused an increase in the mean of number of turns and switchiness, as well as an increase in the variability of all behaviors except turn bias (Figs 4 and 5). Surprisingly, we saw that variability in light-choice probability also increased under the mild enrichment treatment, which led us to believe that our jungle gym construction may have also created differential light conditions in the vial or stimulated phototaxis variability via more indirect means.

The six behavior measures examined were chosen largely because they could be measured at scale across many individuals, which is a requirement for measuring effects on variability. However, we can speculate on the ecological relevance of variability in several of the phenotypes measured. We suspect that variability in turn bias could be potentially advantageous for exploration, dispersal and/or foraging via a bet-hedging mechanism. Individuals with stronger turning biases move through the environment with lower effective diffusion constants. If the spatial scale of resources in the environment fluctuates unpredictably, variation in turn bias could reflect a matched strategy of diversifying diffusion constants. Light-choice bias in our assay may be reflective of an escape response since the fly is startled prior to the light choice (Kain et al., 2012). Variability in moving toward light upon being startled could reflect a bet-hedging mechanism as well, if the threats faced by flies are variable, i.e. if it is alternatively advantageous to seek light or dark after a startle. Predator escape behavior has been considered to be a possible bet-hedging trait; for example, clonal pea aphids show variability in predator escape behavior among individuals (Schuett et al., 2011). We have also found that variability in light preference under non-startled conditions can influence the thermal experience of a fruit fly in nature and therefore could also be part of a bet-hedging strategy (Kain et al., 2015).

We examined whether the effects on variability would change with a different type of enrichment. We raised one cohort of flies in a naturalistic setting, subject to the environmental fluctuations of the outdoors and with access to numerous organic and inorganic substrates (Fig. 1). Compared with the jungle gyms of the mild enrichment, this intense enrichment treatment had more structural complexity and diversity of biotic (fruits, plants, spider predators) and abiotic (sunlight, temperature) factors. By increasing the microenvironmental diversity along several different axes, we expected to observe stronger effects on variability on flies reared in this treatment. In general, this intense enrichment did have stronger effects on our behavioral measures than the mild enrichment. Even though the effects of intense enrichment were more pronounced, the direction of these effects was behavior dependent (Fig. 3). For example, we saw a decrease in intragenotypic variability for turn bias and turn switchiness under intense enrichment, but an increase in the variability for number of turns and turn clumpiness. The directions of these effects varied by behavior, even relative to the direction of the mild enrichment effect.

The intense enrichment treatment was created to have a higher level of microenvironmental diversity than the mild enrichment treatment, yet we observed that the behavioral variability did not always change in the same direction between these treatments. For example, we saw a small increase in variability under mild enrichment but a large decrease under intense enrichment, such as for turn bias (Fig. 3). This leads us to believe that the relationship between the mild and intense enrichment is not just a simple increase in ‘enrichment intensity’. However, we recorded some direct evidence that flies in the intense enrichment indeed experienced at least one dimension of increased microenvironmental variation compared with the mild enrichment: flies recovered from the intense enrichment cage exhibited variation in the color of their gut contents, consistent with their having recently fed on different food sources (Fig. S6). This variation was absent in flies subject to mild enrichment. One of our predictions was that a more naturalistic enrichment treatment could lead to a decrease in variability because of an increase in robustness, but it could also be that naturalistic enrichments cause fly populations to exhibit more natural behaviors in general, whether or not that corresponds to a decrease in behavioral variability. Future studies could address what constitutes natural fly behavior in more detail, whether by making field-deployable assays or bringing wild flies directly to the lab for testing, although any comparisons with our current enrichment paradigm would need to carefully consider population genotypic variance.

Overall, our results support the hypothesis that enrichment increases intragenotypic variability, although this effect is highly dependent on the particular genotype, enrichment and behavior in question. We also conclude that genotype is likely the main determinant of intragenotypic variability. Our findings make it apparent that the genotype used and behavior measured will affect the inferred relationship between environmental variability and behavioral variability. Moreover, the type of enrichment (e.g. mild vs. intense enrichment) can qualitatively and quantitatively alter this relationship. This, and the effects of genotype and behavioral measure, could be why effects observed in one enrichment study may not be replicated in another (Toth et al., 2011). While the multifactorial nature of enrichment provides challenges, its specific effects continue to be of great interest in behavioral research, and high-throughput data-driven approaches have the potential to illuminate the complex relationships between environmental variability and behavioral variability.

We are grateful to Kyle Honegger for discussions of Bayesian inference and Nico O. Wagner for feedback on the manuscript. We also thank Carolyn Elya for moral support when the sampler misbehaved.