Behavior is one of the most variable and rapidly evolving phenotypes, with notable differences even between closely related species (Lorenz, 1958; Martins, 1996). Variable behaviors and rapid behavioral evolution likely facilitates adaptation to new or varying environments and speciation (Baier and Hoekstra, 1914; West-Eberhard, 2003). Despite the importance of animal behavior, progress in revealing the genetic basis of behavioral evolution has been slow (Gleason and Ritchie, 2004; Yamamoto and Ishikawa, 2013; Ellison et al., 2011; Shaw and Lesnick, 2009). In contrast, recent decades have seen significant progress in understanding the genetic causes of morphological evolution (Williams and Carroll, 2009; Shubin et al., 2009; Levine and Davidson, 2005; Stern and Frankel, 2013).

While there are many potential reasons for the discrepancy between studies of behavioral and morphological evolution, including the lack of a fossil record for behavior, a key difficulty has been identifying which aspects of an animal’s development and physiology are the proximate causes of behavior evolution. Evolutionary changes in behavior could emerge from alterations in the developmental patterning of neural circuits (e.g., brain networks, descending commands, central pattern generators), changes in hormonal regulation that influence neural activity, or even from changes in non-neuronal morphology (Baker et al., 2001; Massey et al., 2019). Each of these possibilities could result in behavioral effects at different, yet overlapping, timescales – from muscle twitches to stereotyped suites of behaviors to longer-lived states like foraging or courtship or aging that may control the relative frequency of a given behavior. This complexity may make it difficult to identify the precise aspects of behavior that have evolved.

To address these difficulties, the standard approach in the genetic study of behavioral evolution has been to identify focal behaviors that exhibit robust differences between species, such as courtship behavior in fruit flies (Cande et al., 2012; Cande et al., 2014; Ding et al., 2019) or burrow formation in deermice (Weber et al., 2013; Hu and Hoekstra, 2017). It has been possible to identify genomic regions that correlate with quantitative changes in focal behaviors. However, usually multiple genomic regions are identified, each containing many genes. Given the large number of putative genes involved, combined with the possibility of epistatic interactions between loci, identification of the contributions of individual genes to behavioral evolution has progressed slowly.

An alternative approach to focusing on single behaviors is to examine the full repertoire of movements that an animal performs. By identifying sets of behaviors that evolve together, as was recently performed for hand-tuned traits in a study of birds-of-paradise evolution (Ligon et al., 2018), it may be possible to identify regulators of these suites of behaviors. This approach has been thwarted by the challenge of robustly measuring multiple behavioral phenotypes simultaneously. Recent progress in the unsupervised identification of animal behaviors across length and time scales, however, has made this approach possible (Berman, 2018; Brown and de Bivort, 2018). In this study, we introduce a quantitative framework for studying the evolutionary dynamics of large suites of behavior. We have focused initially on fruit flies, which provide a convenient model for this problem – both because they exhibit a wide range of complex behaviors and because unsupervised approaches can be used to map all of the animal movements captured in video recordings (Berman et al., 2014; Cande et al., 2018; Berman et al., 2016).

We recorded movies of isolated male flies from six species in a nearly stimulus-deprived environment. Because we did not record flies experiencing social and other environmental cues, we did not observe many charismatic natural behaviors, such as courtship and aggression. Nevertheless, we found that the behaviors they performed, including walking and grooming, contain species-specific information. We thus hypothesized that our quantitative representations of behaviors could be studied in an evolutionary context. To infer the evolutionary trajectories of behavioral evolution, we estimated ancestral behavioral repertoires with a Generalized Linear Mixed Model (GLMM) approach (Hadfield, 2010), which builds upon Felsenstein’s approach to reconstructing ancestral states (Felsenstein, 1985; Felsenstein, 2005; Hadfield and Nakagawa, 2010; O’Meara, 2012). Using these results, we develop a framework that allows us to model the behavioral traits that co-vary both within a species and along the phylogeny. We found that within-species variance has a similar structure to long-lasting internal states of the animal that we characterized previously, and that inter-species variance can capture how disparate behaviors may have evolved together. This latter finding points toward the presence of higher order behavioral traits that would not have been detected by studying individual behaviors in isolation and that may be amenable to further evolutionary and genetic analysis.

Experiments and behavioral quantification

We captured video recordings of all behaviors performed by single flies isolated in a largely featureless environment for multiple individuals from six species of the Drosophila melanogaster species subgroup: D. mauritiana, D. melanogaster, D. santomea, D. sechellia, D. simulans, and D. yakuba (Cande et al., 2018). Although the animals could not jump or fly in these chambers and were not expected to exhibit social or feeding behaviors, the flies displayed a variety of complex behaviors, including locomotion and grooming. Each of these behaviors involves multiple body parts that move at varying time scales. The species studied here were chosen because their phylogenetic relationships are well understood (Clark et al., 2007; Obbard et al., 2012; Chyb and Gompel, 2013; Seetharam and Stuart, 2013) (summarized in the tree seen in Figure 3), and genetic tools are available for most of these species (Stern et al., 2017). Since a single strain represents a genomic ‘snapshot’ of each species, we assayed multiple individuals from each of multiple strains of each species to attempt to capture species-specific differences, and not variation specific to particular strains (see Materials and methods). In total, we collected data from 561 flies, each measured for an hour at a sampling rate of 100 Hz.

While previous studies have identified differences in specific behaviors, such as courtship behavior, between these species (Cande et al., 2012; Ding et al., 2019; Yamamoto and Ishikawa, 2013; Auer and Benton, 2016), here we assayed the full repertoire of behaviors the flies performed in the arena, with the aim of identifying combinations of behaviors that may be evolving together. To measure this repertoire, we used a previously described behavior mapping method (Berman et al., 2014; Cande et al., 2018) that starts from raw video images and finds each animal’s stereotyped movements in an unsupervised manner. The output of this method is a two-dimensional probability density function (PDF) that contains many peaks and valleys (Figure 1A), where each peak corresponds to a different stereotyped behavior (e.g., right wing grooming, proboscis extension, running, etc).

Figure 1

Download asset Open asset

Briefly, to create the density plots, raw video images were rotationally and translationally aligned to create an egocentric frame for the fly. The transformed images were decomposed using Principal Components Analysis into a low-dimensional set of time series. For each of these postural mode time series, a Morlet wavelet transform was applied, obtaining a local spectrogram between 1 Hz and 50 Hz (the Nyquist frequency). After normalization, each point in time was mapped using t-SNE (van der Maaten and Hinton, 2008) into a two-dimensional plane. Finally, convolving these points with a two-dimensional gaussian and applying the watershed transform (Meyer, 1994), produced 134 different regions, each of these containing a single local maximum of probability density that corresponds to a particular stereotypical behavior. We integrate over this local region of the probability density to calculate the probability that a fly is performing this behavior at a random point in time. Thus, we can associate each fly with a 134-dimensional real-valued vector that represents the probability of the fly performing a certain stereotyped behavior at a given time during the hour-long experimental session. We will refer to this quantity as the animal’s behavioral vector, $\overrightarrow{P}$.

The behavioral map averaged across all six species is shown in Figure 1A and displays a pattern of movements similar to those we found in previous work, where locomotion, idle/slow, anterior/posterior movements, etc. are segregated into different regions (Berman et al., 2014; Cande et al., 2018). Averaging across all individuals of each species, we found the mean behavioral vector for each species (Figure 1B) and observed that each species performs certain behaviors with different probabilities. For example, D. mauritiana individuals spend more time performing fast locomotion than all other species on average, and D. yakuba individuals spend much of their time performing an almost species-unique type of slow locomotion, but little time running quickly.

These average probability maps provide some insight into potential species differences, but to identify species-specific behaviors, we also need to account for variation in the probability that individuals of each species perform each behavior. One way to address this problem is to ask whether an individual’s species identity can be predicted solely from its multi-dimensional behavioral vector. To explore this question, we first used t-SNE to project all 561 individuals into a two-dimensional plane (Figure 2A), using the Jensen-Shannon divergence as the distance metric between individual behavioral vectors. In this plot, different colors represent different species, and different symbols with the same color represent different strains within the same species. Although species do not segment cleanly into separate regions of this plane, the distribution of species is far from random, with individuals from the same species tending to group near to one other. Given this structure, there is likely species-specific information in the behavioral vectors.

Figure 2

Download asset Open asset

To quantify this observation, we applied a multinomial logistic regression classifier to the data, performing a six-way classification based solely on the high-dimensional behavioral vectors. After training, the classifier correctly classified $89\pm .2\%$ of vectors in our test set (a randomly selected 30% of the entire data set that was not used during training). Moreover, the confusion matrix (Figure 2B) revealed no systematic misclassification bias amongst the species. Note that we have used a relatively simple classifier compared to modern deep learning methods (Goodfellow et al., 2016), so these results likely represent a lower bound on the distinguishability of the behavioral vectors. Thus, the behavioral vectors contain considerable species-specific information. We therefore proceeded to explore how these behavioral vectors may have evolved along the phylogeny.

Reconstructing ancestral behavioral repertoires

Multiple methods have been proposed for reconstructing ancestral states from data collected from extant species (Felsenstein, 1985; Felsenstein, 2005; Yang, 2006; O’Meara, 2012; Royer-Carenzi and Didier, 2016). These methods generally fall into two camps: parsimony reconstruction, which attempts to reconstruct evolutionary history with the fewest number of evolutionary changes (Cunningham et al., 1998), and diffusion-processes, which model evolution as a random walk on a multi-dimensional landscape (Hadfield and Nakagawa, 2010). Given the high-dimensional behavioral vectors that we are attempting to model, a diffusion process is more likely to capture the inter-trait correlations that we would like to understand. Thus, we focus on a diffusion-based model here.

Given a phylogeny for a collection of species, we modeled how species-specific complexes of behaviors might have emerged. We assumed that each animal’s behavior is a quantitative trait with an additive random effect, that is, each animal’s behavior is a trait that results from the additive effects of many genetic loci, each of small effect, that is combined with a non-genetic effect that represents inter-specific variation. We do not, however, assume that all behaviors evolve independently of each other. Thus, we are interested in predicting (1) whether intra- and inter-species variation can be separated to identify independently evolving sets or linear combinations of behaviors and (2) how behaviors co-vary along the phylogeny, potentially revealing co-evolving suites of behaviors.

We assumed that the observed flies’ behaviors evolved via a diffusion process with Gaussian noise from a common ancestor along the known phylogenetic tree. Note that this is a less restrictive assumption than neutrality, as multiple traits under selection may evolve in a correlated manner. Specifically, we fit a multi-response Generalized Linear Mixed Model (GLMM) to the data, using the approach described in Hadfield, 2010, modeling the evolutionary process such that the logarithm of the behavioral vector, $\overrightarrow{P}$, for each individual ($\overrightarrow{l}=(l_1,\mathrm{\dots },l_{K=134})$), is given by

(1) $\overrightarrow{l}=\overrightarrow{\mu }+\overrightarrow{\rho }+\overrightarrow{e},$

where $\overrightarrow{\mu }$ is the mean behavior of the common ancestor (treated as the fixed effects of this model), and $\overrightarrow{\rho }$ and $\overrightarrow{e}$ are the random effects corresponding to the phylogenetic and individual variability, respectively. We assume that these random effects are generated from the multi-dimensional normal distributions ${\displaystyle \mathcal{𝒩}(\overrightarrow{0},A\otimes V^{(a)})}$ (phylogenetic) and ${\displaystyle \mathcal{𝒩}(\overrightarrow{0},I\otimes V^{(e)})}$ (individual). Here, the matrix $A$ represents the information contained in the phylogenetic tree, with $A_{ij}$ being proportional to the length of the path from the most recent common ancestor of species $i$ and $j$ to the common ancestor. $A$ is normalized so that the diagonal elements are all equal to 1. Therefore, $A_{ij}$ represents the phylogenetic similarity between a pair of species. $I$ is the identity matrix, and $V^{(a)}$ and $V^{(e)}$ are the phylogenetic and within-species covariance matrices, respectively.

We fit μ, $V^{(a)}$, and $V^{(e)}$ using Markov Chain Monte Carlo (MCMC) simulations, confirming that the MCMC converged using the Gelman-Rubin diagnostic (see Materials and methods, Figure 3—figure supplement 1). In addition, our model is able to infer the mean endpoint behavioral repertoires (Figure 3—figure supplement 2), providing confidence that our model is consistent with our input data. In addition to the inferred behavioral states corresponding to the common ancestor, ${\overline{P}}^{Anc}$, we also reconstructed the mean behavioral representations for the intermediate ancestors (Figure 3).

Figure 3 with 3 supplements see all

Download asset Open asset

We also found that the model that allows behavioral co-evolution out-performs a model where each behavioral trait evolves independently. Specifically, we fit a model where behavioral correlations between individuals of different species were removed by enforcing that $V^{(a)}$ and $V^{(e)}$ must be diagonal matrices, a reduction of more than 17,500 parameters compared to the full model (see Materials and methods for details). The phylogenetic ancestral reconstruction was then made for each behavioral trait separately. To compare the relative performance of these models, we used the Deviance Information Criterion (DIC) (Spiegelhalter et al., 2002), a commonly used assessment tool for MCMC-fit hierarchical models that lack a good estimate for the number of effective parameters. Like other information-theoretic model selection criteria (e.g., the Akaike Information Criterion or the Bayesian Information Criterion), smaller values of the DIC imply a larger posterior probability of the model given the available data. Despite the large reduction in the number of parameters for the independent-trait model, the DIC for the independent-trait model was substantially higher ($\text{DIC}=(242\pm 2)\times {10}^3$) than the DIC for the full model ($\text{DIC}=(114\pm 2)\times {10}^3$). Moreover, the full model was able to predict the inter- and intra-species covariances between dissimilar pairs of behaviors (Figure 3—figure supplement 3). Hence, modeling the evolution of the full behavioral repertoire captures the structure of the observed data better than a single trait approach.

Individual variability and long timescale correlations

While it is not possible to directly test the accuracy of our ancestral state reconstructions (Figure 3), the inferred covariance matrices generate predictions about behavioral and genetic correlations that are, in principle, testable. We therefore focus on the fitted covariance matrices, $V^{(e)}$ and $V^{(a)}$ (each in $\mathrm{I}{\mathrm{R}}^{134\times 134}$), which account for within-species and phylogenetic random effects, respectively.

We will focus first on the intra-species covariance matrix, $V^{(e)}$. We note first that the matrix exhibits a modular structure (Figure 4A). After rearranging the behavior order via an information-based clustering procedure (Slonim et al., 2005), we see that a block diagonal pattern emerges, with positive correlations lying within the blocks and negative correlations lying off the diagonal. The details of this particular clustering approach are described in Materials and methods, but we find that the results are nearly identical for several different clustering methodologies (Figure 4—figure supplement 1). Quantifying the matrix’s modularity via the average within-cluster dissimilarity,

(2) $⟨d⟩=\sum _{k}p(C_k)\sum _{i,j\in C_k}\frac{1}{2}\left[1-\frac{V_{ij}^{(e)}}{\sqrt{V_{ii}^{(e)}{V}_{jj}^{(e)}}}\right],$

Figure 4 with 3 supplements see all

Download asset Open asset

where $C_k$ is the set of all behaviors belonging to the $k$ th cluster, we find that $⟨d⟩\approx 0.30$ and 0.22 for the 3- and 6-cluster solutions, respectively. These values are significantly smaller than the average distances obtained using random cluster assignments ($⟨d⟩=0.46\pm 0.03$ and $0.45\pm 0.04$ for 3 and 6 clusters respectively, see Figure 4—figure supplement 2 for other numbers of clusters and clustering methods). Thus, we can conclude that the intra-specific covariance matrix has a far-from-random modular structure, implying that between individuals of the same species, groups of behaviors tend to vary together in a stereotyped manner.

Moreover, these groups of behaviors that co-vary together within a species are not random collections of behaviors. Instead, we found that co-varying clusters are spatially contiguous in the behavioral map, implying that covariances of groups of similar behaviors (behaviors involving moving similar parts of the animals’ bodies at similar speeds) compose much of the observed intra-species variance. The clustering method does not take the spatial structure of the behavioral map into account at all (just the values in $V^{(e)}$), so the clusters of local behaviors in the behavior map reflect underlying similarity in the covariance of nearby behaviors, rather than an artifact of the algorithm. Moreover, co-varying clusters are hierarchically organized, where coarse-grained co-varying behaviors can be sub-divided into smaller co-varying clusters (Figure 4B), a feature that is not guaranteed by the information-based clustering algorithm.

This hierarchical structure of the behavioral map is reminiscent of the hierarchical temporal structure of behavior that was hypothesized originally by ethologists (Tinbergen, 1951; Deutsch et al., 2020) and was observed to optimally explain the history-dependent long timescale non-stationary structure of Drosophila melanogaster behavioral transitions (Berman et al., 2016). Thus, we hypothesized that the structure of the intra-species covariance matrix might be linked to deviations from statistical stationarity in the behavioral data that were not explicitly measured in the unsupervised clustering or modeled in the GLMM.

To explore this connection further, we performed an analysis that is analogous to the single-species study in Berman et al., 2016, finding coarse-grainings of the behavioral space (i.e., a description of the behavioral space using fewer behaviors) that are optimally predictive of the future behaviors that the flies perform. Specifically, if $b(t)$ is the behavior that a fly performs at time $t$, we would like to create a clustered version of our behavioral map, $Z$, such that we maximize the information that $z(t)\in Z$, the cluster that the fly is in at time $t$ contains as much information about the future behavior of the fly, $b(t+\tau )$ as possible. To keep $Z$ from separating each behavior into its own cluster, we also need to make sure that $Z$ is as simple a clustering as possible (i.e., a smaller number of clusters and a more even distribution of time spent within each clusters).

To be more precise, we calculated $Z$ using the the Deterministic Information Bottleneck (DIB) method (Strouse and Schwab, 2017). This approach minimizes the functional

(3) ${𝒥}_{\tau }=-I(Z(t);Z(b+\tau ))+\gamma \mathscr{H}(Z),$

where $b(t+\tau )$ is a fly’s behavior at time $t+\tau$, $Z(t)$ is the coarse-grained behavior visited at time $t$, $I(b(t);Z(t+\tau ))$ is the mutual information between these quantities, $\gamma$ is a positive constant, and $\mathscr{H}(Z)$ is the entropy of the coarse-grained representation (see Materials and methods). As $\gamma$ is increased, progressively simpler, but less predictive, representations are found.

Applying this method to the data, pooled across all six species and using $\tau =50$ (Figure 4C, Figure 4—figure supplement 3), we found the same hierarchical division of the behavioral map that was observed for freely moving D. melanogaster (Berman et al., 2016). Moreover, we found that the structure of the space using this approach closely mirrors the structure found via directly clustering the intra-species covariance matrix, $V^{(e)}$ (Figure 4C). Quantifying the similarity between both clustering partitions by calculating the Weighted Similarity Index (WSI), a modification of the Rand Index (Rand, 1971) (Materials and methods), the WSI between the information-based clustering method and the predictive information bottleneck for three clusters is $WSI=0.73$, and $WSI=0.87$ for six clusters. For random clusterings, we would expect to observe $0.51\pm 0.02$ and $0.70\pm 0.01$ for 3 and 6 clusters, respectively, indicating a non-random overlap between these two partitions. Figure 4—figure supplement 1, shows that this result is independent of the clustering method and the number of clusters.

The overlap between these two coarse-grainings indicates that most individual variability in the behaviors we observe results from non-stationarity in behavioral measurements, rather than from individual-specific variation. That is, much of the intraspecific variation appears to reflect flies recorded when they were experiencing different hidden behavioral states (e.g., circadian state, hunger, etc.), rather than reflecting fixed (environmental or genetic) differences between flies. This variation may have arisen because, although we controlled many variables (e.g., fly age, circadian cycle, temperature, and humidity), it is not possible to control for all internal factors (e.g., hunger, arousal, etc.) that affect an animal’s behavioral patterns (Anderson, 2016). The temporal coarse-graining of the behavioral space that we found via the DIB provides insight into these non-stationarities, as they are optimally predictive of the fly’s future behaviors. Given the contiguous nature of these regions, this result means that flies tended to stay within specific regions of the behavioral space much longer than one would assume from a Markov model, hinting that there is an important connection between variability across animals and variability between animals.

More precisely, these results imply that variation in behavior observed among individuals, especially in non-manipulated settings, may often reflect a large component of hidden behavioral states (Figure 5A). Thus, it may be possible to improve upon behavioral measurements in many settings by controlling for the variability associated with these hidden states. For example, just because one fly performs less anterior grooming than another may reflect that the animal is in a different long timescale behavioral state, rather than that the animal has a genetically encoded preference for reduced grooming.

Figure 5

Download asset Open asset

A potential method for accounting for these artifacts is to normalize each individual’s behavioral density such that the amount of time that the animal spends in each of the coarse-grained regions is equalized. In other words, the amount of time spent anterior grooming, locomoting, etc. are set to be the same for all animals, thus accounting for the variability associated with the inferred hidden states. Mathematically, if $P_i$ is the probability of observing behavior $i$, and $C_i$ is the clustering assignment of this behavior, we can define a normalized probability, ${\widehat{P}}_i$, via

(4) ${\widehat{P}}_i=\frac{{\overline{P}}^{(C_i)}}{P_i^{(C_i)}}{P}_i,$

where $P_i^{(C)}={\sum }_{k\in C}{P}_k$ is the total density in cluster $C$ for an individual fly and ${\overline{P}}^{(C)}$ is the average across all animals.

We found that applying this normalization to our data often results in substantial changes in the inferred distributions of behavioral densities. For example, Figure 5B displays how the difference in behavioral density between D. santomea and D. yakuba (as measured by the Mahalanobis distance between the distributions) alters as a result of normalization. For some behaviors, the signal increases (red points), and in some cases, it decreases (blue points). Thus, it is important to take these non-stationary effects into account when estimating how often single behaviors are performed in studies of behavioral evolution. To measure these non-stationary effects, many behaviors must be measured, not just a focal behavior, thus partially explaining the relative success of our multi-trait model compared to a model where each trait is analyzed independently.

Identifying phylogenetically linked behaviors

One of the advantages of our approach is that we separate variations in behavior corresponding to evolutionary patterns, the phylogenetic variability, from variations among individuals of the same species. By studying the properties of the phylogenetic covariance matrix ($V^{(a)}$), we can identify multiple behaviors that may have evolved together.

We first characterized the coarse-grained structure within $V^{(a)}$ through the information-based clustering used in the previous section (Slonim et al., 2005) (see Materials and methods). As seen in Figure 6A, the phylogenetically co-varying clusters are not spatially contiguous in the behavioral map. This finding is in contrast to the spatial contiguity we observed for the intra-species covariance matrix (Figure 4B). For example, the two-cluster solution (Figure 6A, left) groups the behavioral space into side legs movements (middle of the behavioral map) and certain locomotion gaits (far left of the behavioral map) versus the rest of behaviors. Similarly, when the matrix is clustered into a larger number of clusters, correlated groups are not contiguously arranged within the behavior map. Thus, our model predicts that many non-similar behaviors are evolving in a correlated manner.

Figure 6

Download asset Open asset

To quantify these patterns as traits, we decomposed $V^{(a)}$ via an eigendecomposition. As seen in Figure 6B, almost all of the variance within the matrix can be explained with only the first six eigenmodes. These eigenvectors (Figure 6C) share similar non-local structure to the clusterings described above. By projecting individual behavioral vectors onto these eigenvectors, the resulting dot products represent a meta-trait that is a linear combination of phylogenetically linked behaviors.

These evolving meta-traits may be suitable targets for further neurobiological or genetic studies. Three examples of these distributions are shown in Figure 6D–F for several pairs of closely related species. These three examples were not chosen at random, but instead because they showed significant differentiation between species. The aim of this analysis is not to show that all meta-traits would differ between all pairs of species, which is unlikely, but rather that it is possible to identify synthetic meta-traits that could be further interrogated with experimental methods.

We have developed a quantitative framework to study the evolution of behavioral repertoires, using fruit flies (Drosophila) as a model system. We started with observations of 561 individuals from six extant species behaving in an unremarkable environment. This assay did not include social behaviors, such as courtship and aggression, nor many foraging behaviors. Thus, at first glance, it might seem like we had excluded most species-specific behaviors from the analysis. Nonetheless, we found that other complex behaviors, like walking, running, and grooming, exhibit species-specific features that can be used to reliably assign individuals to the correct species. Thus, the motor patterns of behaviors that are not normally investigated for their species-specific features are likely evolving between even closely related species. It is not clear, however, if these differences reflect natural selection or genetic drift. All of these behaviors are critical to individual survival, however, so it is possible that these behaviors have evolved, at least in part, in response to natural selection. It is clear, however, that the underlying mechanisms, and perhaps the neural circuitry, controlling these behaviors must have evolved.

Inspired by these observations, we estimated patterns of behavioral evolution in the context of a well-understood phylogeny. We fit a Generalized Mixed Linear Model to our behavioral measurements and the given phylogeny to reconstruct ancestral behavioral repertoires and the intra- and inter-species covariance matrices. We found that the patterns of intra-species variability are similar to long timescale behavioral dynamics that violate statistical stationarity - a result we reported previously in a study of a single species (Berman et al., 2016). This result suggests that much of the intra-specific variability that emerged by sampling flies under well-controlled conditions reflects variability in the hidden behavioral states of individual flies. While it may be challenging to conceptualize that seemingly simple behaviors, like the pace of walking and running, are reflective of an underlying long time scale behavior state, many short time scale behaviors, such as the individual movements involved in grooming (Seeds et al., 2014), courtship (Calhoun et al., 2019; Deutsch et al., 2020) and aggression (Hoopfer, 2016; Duistermars et al., 2018) reflect behaviors performed only, or mainly, in the context of a longer lasting behavioral state. These types of long timescale variability may be a statistical confound for evolutionary and experimental studies of behavior. We therefore propose a method to control for these internal states by normalizing the frequency of behaviors relative to an estimate of an animal’s non-stationary states. This method improved the accuracy of behavioral phenotyping and dramatically altered estimates of some species-specific behaviors. For more focused studies, it may not be necessary to measure the full suite of behaviors to effectively normalize for behavioral state, since state can sometimes be estimated from a smaller number of behaviors. In fact, targeted studies of charismatic behaviors, including behaviors associated with aggression or courtship, often implicitly normalize by behavioral state.

Given our estimates for how suites of behaviors evolved, we examined whether the inter-species covariance matrix could be used to identify behavioral meta-traits that might be subjected to further evolutionary and experimental analysis. We identified multiple suites of behaviors that differed between closely related species, providing a starting point for further analysis of how the mechanisms underlying these suites of behaviors have evolved.

There are multiple possible interpretations of these phylogenetically correlated behaviors. For instance, at the neural level, each of these groups of movements may reflect a motor response to shared upstream commands (Cande et al., 2018). Here, for example, different types of locomotion might be controlled through the same descending neural circuitry, but due to evolutionary changes, the same commands could lead to different behavioral outputs, as has been observed in fly courtship patterns (Ding et al., 2019). Alternatively, at the genetic level, multiple behaviors may be linked by pleiotropic effects of individual genetic changes. Finally, groups of co-evolving traits may not be linked mechanistically, but co-evolution may instead reflect selection on suites of behaviors. For example, the male neurons that drive fly courtship song production in the ventral nerve cord are unlikely to be related to the female neurons in the central brain that perceive and interpret the courtship song. Nonetheless, these traits co-evolve such that females tend to prefer songs produced by males of their own species (Bennet-Clark and Ewing, 1969; Ding et al., 2019).

The analysis framework introduced here represents the first attempt to analyze full behavioral repertoires to gain insight into evolution. In principle, this approach could be applied to any data set where a large number of behaviors have been sampled in many species. However, there are several areas where one could add future improvements to this approach. First, we recorded behavior from only six species of flies. Adding additional species would place more constraints on the evolutionary dynamics, likely resulting in less variance in the ancestral state estimations and potentially adding more structure to the relatively low rank (i.e., highly modular) covariance matrices. Additionally, further work is required to determine the balance between sampling within and between strains and species that optimizes estimates of evolutionary dynamics.

Second, our framework assumes that all evolutionary changes in behavior resemble a diffusion process. Although this assumption is a reasonable initial hypothesis (Felsenstein, 1985), it may be possible to test this assumption. For example, deeper sampling of additional species may allow identification of specific behaviors on particular lineages where neutrality can be rejected (Tajima, 1993). If evidence emerges that the analyzed behaviors do not evolve under a diffusion process but under stabilizing selection, for example, the model for ancestral reconstruction can be changed from a Brownian motion to an Ornstein-Uhlenbeck process (Martins and Hansen, 1997; Hansen and Martins, 1996; Royer-Carenzi and Didier, 2016). Such a change can be implemented by altering the structure of the phylogenetic matrix, $A$(Martins and Hansen, 1997; Caetano and Beaulieu, 2020), but without other alterations to the overall methodology presented here.

Another potential limitation of our analysis is that some of the observed inter-specific differences may reflect species-specific responses to environmental factors like room temperature or humidity, rather than underlying genetic or developmental factors. Two observations mitigate against this possibility, however. First, there were no significant differences in overall activity level of the different species, which would be a key indicator of environment-induced covariance. Second, the intra-species covariance matrix (derived from data from all species) agrees well with previous findings within a single species (Berman et al., 2016), implying that many of the potential environmental co-varying factors are shared across all six species.

In addition, all of our current analyses ignored the temporal structure of behavior and sequences of movements. While we found that the intra-specific variance has a similar structure to temporal structure that we reported previously (Berman et al., 2016), the order in which behaviors occur may also provide important biological information, especially during events like courtship or aggression. It should be possible to incorporate temporal structure directly into the regression (Caetano and Beaulieu, 2020). Deciding exactly which quantities to measure and how they should be incorporated, however, are complex questions that are outside the scope of this initial study. In addition, the number of fit parameters, already over 18,000 here, would need to grow even larger to accommodate modeling transition rates between the behaviors as traits themselves. Thus, fitting such models would necessitate even larger data sets than the one collected here. Moreover, because the Perron-Frobenius Theorem mathematically couples the transition probabilities between behaviors and the probabilities of the fly performing a given behavior, additional care (and data) is required to ensure that observed differences in behavior are due to changes in temporal structure rather than changes in the frequencies of performing a given behavior.

Lastly, capturing the full range of animal behaviors for a large number of animals presents a number of technological challenges, which is why we focused on measuring behavior in a highly simplified environment. However, a more complete understanding of the structure of behavior will require more sophisticated ways to capture behavioral dynamics in more naturalistic settings and during complex social arrangements. While modern deep learning methods have made tracking animals in more realistic settings increasingly plausible (Pereira et al., 2019; Mathis and Mathis, 2020), there are still considerable hurdles to translating this information into a form that can be subjected to the kind of analysis we propose here.

Despite these limitations, this work represents a new way to quantitatively characterize the evolution of complex behaviors, which may provide new phenotypes that can be subjected to experimental analysis. In the absence of a behavioral fossil record, reconstructing ancestral behaviors requires an inferential approach like the one we present here. In addition, more complex models could be built to test assumptions of the diffusion-based model we employed. Finally, a strength of our approach is that it makes falsifiable predictions about how behaviors are linked mechanistically, providing predictions that can be tested experimentally to provide further insight into the genetic and neurobiological structure of behavior.

All behavioral region information is submitted with the article and is posted on GitHub (https://github.com/bermanlabemory/behavioral-evolution, copy archived at https://archive.softwareheritage.org/swh:1:rev:b01a6e3a2c7da193f38631dfe925c65229494d74). The original video data are too large to post (tens of TB), but will be made available upon request.

1. Conference

   1. Akaike H

   (1973)

   Information theory and an extension of the maximum likelihood principle

   Proceedings of the 2nd International Symposium on Information Theory. pp. 267–281.

   • Google Scholar
2. 1. Anderson DJ

   (2016) Circuit modules linking internal states and social behaviour in flies and mice

   Nature Reviews Neuroscience 17:692–704.

   https://doi.org/10.1038/nrn.2016.125

   • PubMed
   • Google Scholar
3. 1. Auer TO
   2. Benton R

   (2016) Sexual circuitry in Drosophila

   Current Opinion in Neurobiology 38:18–26.

   https://doi.org/10.1016/j.conb.2016.01.004

   • PubMed
   • Google Scholar
4. 1. Baier F
   2. Hoekstra HE

   (1914) The genetics of morphological and behavioural island traits in deer mice

   Proceedings of the Royal Society B 286:1697.

   https://doi.org/10.1098/rspb.2019.1697

   • Google Scholar
5. 1. Baker BS
   2. Taylor BJ
   3. Hall JC

   (2001) Are complex behaviors specified by dedicated regulatory genes? reasoning from Drosophila

   Cell 105:13–24.

   https://doi.org/10.1016/S0092-8674(01)00293-8

   • PubMed
   • Google Scholar
6. 1. Bennet-Clark HC
   2. Ewing AW

   (1969) Pulse interval as a critical parameter in the courtship song of Drosophila melanogaster

   Animal Behaviour 17:755–759.

   https://doi.org/10.1016/S0003-3472(69)80023-0

   • Google Scholar
7. 1. Berman GJ
   2. Choi DM
   3. Bialek W
   4. Shaevitz JW

   (2014) Mapping the stereotyped behaviour of freely moving fruit flies

   Journal of the Royal Society Interface 11:20140672.

   https://doi.org/10.1098/rsif.2014.0672

   • Google Scholar
8. 1. Berman GJ
   2. Bialek W
   3. Shaevitz JW

   (2016) Predictability and hierarchy in Drosophila behavior

   PNAS 113:11943–11948.

   https://doi.org/10.1073/pnas.1607601113

   • PubMed
   • Google Scholar
9. 1. Berman GJ

   (2018) Measuring behavior across scales

   BMC Biology 16:23.

   https://doi.org/10.1186/s12915-018-0494-7

   • PubMed
   • Google Scholar
10. 1. Brooks SP
    2. Gelman A

    (1998) General methods for monitoring convergence of iterative simulations

    Journal of Computational and Graphical Statistics 7:434–455.

    https://doi.org/10.1080/10618600.1998.10474787

    • Google Scholar
11. 1. Brown AEX
    2. de Bivort B

    (2018) Ethology as a physical science

    Nature Physics 14:653–657.

    https://doi.org/10.1038/s41567-018-0093-0

    • Google Scholar
12. 1. Caetano DS
    2. Beaulieu JM

    (2020) Comparative analyses of phenotypic sequences using phylogenetic trees

    The American Naturalist 195:E38–E50.

    https://doi.org/10.1086/706912

    • PubMed
    • Google Scholar
13. 1. Calhoun AJ
    2. Pillow JW
    3. Murthy M

    (2019) Unsupervised identification of the internal states that shape natural behavior

    Nature Neuroscience 22:2040–2049.

    https://doi.org/10.1038/s41593-019-0533-x

    • Google Scholar
14. 1. Cande J
    2. Andolfatto P
    3. Prud'homme B
    4. Stern DL
    5. Gompel N

    (2012) Evolution of multiple additive loci caused divergence between Drosophila yakuba and D. santomea in wing rowing during male courtship

    PLOS ONE 7:e43888.

    https://doi.org/10.1371/journal.pone.0043888

    • PubMed
    • Google Scholar
15. 1. Cande J
    2. Stern DL
    3. Morita T
    4. Prud’homme B
    5. Gompel N

    (2014) Looking Under the Lamp Post: Neither fruitless nor doublesex Has Evolved to Generate Divergent Male Courtship in Drosophila

    Cell Reports 8:363–370.

    https://doi.org/10.1016/j.celrep.2014.06.023

    • Google Scholar
16. 1. Cande J
    2. Namiki S
    3. Qiu J
    4. Korff W
    5. Card GM
    6. Shaevitz JW
    7. Stern DL
    8. Berman GJ

    (2018) Optogenetic dissection of descending behavioral control in Drosophila

    eLife 7:e34275.

    https://doi.org/10.7554/eLife.34275

    • Google Scholar
17. Book

    1. Chyb S
    2. Gompel N

    (2013)

    Atlas of Drosophila Morphology: Wild-Type and Classical Mutants

    London: Academic Press.

    • Google Scholar
18. 1. Clark AG
    2. Eisen MB
    3. Smith DR
    4. Bergman CM
    5. Oliver B
    6. Markow TA
    7. Kaufman TC
    8. Kellis M
    9. Gelbart W
    10. Iyer VN
    11. Pollard DA
    12. Sackton TB
    13. Larracuente AM
    14. Singh ND
    15. Abad JP
    16. Abt DN
    17. Adryan B
    18. Aguade M
    19. Akashi H
    20. Anderson WW
    21. Aquadro CF
    22. Ardell DH
    23. Arguello R
    24. Artieri CG
    25. Barbash DA
    26. Barker D
    27. Barsanti P
    28. Batterham P
    29. Batzoglou S
    30. Begun D
    31. Bhutkar A
    32. Blanco E
    33. Bosak SA
    34. Bradley RK
    35. Brand AD
    36. Brent MR
    37. Brooks AN
    38. Brown RH
    39. Butlin RK
    40. Caggese C
    41. Calvi BR
    42. Bernardo de Carvalho A
    43. Caspi A
    44. Castrezana S
    45. Celniker SE
    46. Chang JL
    47. Chapple C
    48. Chatterji S
    49. Chinwalla A
    50. Civetta A
    51. Clifton SW
    52. Comeron JM
    53. Costello JC
    54. Coyne JA
    55. Daub J
    56. David RG
    57. Delcher AL
    58. Delehaunty K
    59. Do CB
    60. Ebling H
    61. Edwards K
    62. Eickbush T
    63. Evans JD
    64. Filipski A
    65. Findeiss S
    66. Freyhult E
    67. Fulton L
    68. Fulton R
    69. Garcia AC
    70. Gardiner A
    71. Garfield DA
    72. Garvin BE
    73. Gibson G
    74. Gilbert D
    75. Gnerre S
    76. Godfrey J
    77. Good R
    78. Gotea V
    79. Gravely B
    80. Greenberg AJ
    81. Griffiths-Jones S
    82. Gross S
    83. Guigo R
    84. Gustafson EA
    85. Haerty W
    86. Hahn MW
    87. Halligan DL
    88. Halpern AL
    89. Halter GM
    90. Han MV
    91. Heger A
    92. Hillier L
    93. Hinrichs AS
    94. Holmes I
    95. Hoskins RA
    96. Hubisz MJ
    97. Hultmark D
    98. Huntley MA
    99. Jaffe DB
    100. Jagadeeshan S
    101. Jeck WR
    102. Johnson J
    103. Jones CD
    104. Jordan WC
    105. Karpen GH
    106. Kataoka E
    107. Keightley PD
    108. Kheradpour P
    109. Kirkness EF
    110. Koerich LB
    111. Kristiansen K
    112. Kudrna D
    113. Kulathinal RJ
    114. Kumar S
    115. Kwok R
    116. Lander E
    117. Langley CH
    118. Lapoint R
    119. Lazzaro BP
    120. Lee SJ
    121. Levesque L
    122. Li R
    123. Lin CF
    124. Lin MF
    125. Lindblad-Toh K
    126. Llopart A
    127. Long M
    128. Low L
    129. Lozovsky E
    130. Lu J
    131. Luo M
    132. Machado CA
    133. Makalowski W
    134. Marzo M
    135. Matsuda M
    136. Matzkin L
    137. McAllister B
    138. McBride CS
    139. McKernan B
    140. McKernan K
    141. Mendez-Lago M
    142. Minx P
    143. Mollenhauer MU
    144. Montooth K
    145. Mount SM
    146. Mu X
    147. Myers E
    148. Negre B
    149. Newfeld S
    150. Nielsen R
    151. Noor MA
    152. O'Grady P
    153. Pachter L
    154. Papaceit M
    155. Parisi MJ
    156. Parisi M
    157. Parts L
    158. Pedersen JS
    159. Pesole G
    160. Phillippy AM
    161. Ponting CP
    162. Pop M
    163. Porcelli D
    164. Powell JR
    165. Prohaska S
    166. Pruitt K
    167. Puig M
    168. Quesneville H
    169. Ram KR
    170. Rand D
    171. Rasmussen MD
    172. Reed LK
    173. Reenan R
    174. Reily A
    175. Remington KA
    176. Rieger TT
    177. Ritchie MG
    178. Robin C
    179. Rogers YH
    180. Rohde C
    181. Rozas J
    182. Rubenfield MJ
    183. Ruiz A
    184. Russo S
    185. Salzberg SL
    186. Sanchez-Gracia A
    187. Saranga DJ
    188. Sato H
    189. Schaeffer SW
    190. Schatz MC
    191. Schlenke T
    192. Schwartz R
    193. Segarra C
    194. Singh RS
    195. Sirot L
    196. Sirota M
    197. Sisneros NB
    198. Smith CD
    199. Smith TF
    200. Spieth J
    201. Stage DE
    202. Stark A
    203. Stephan W
    204. Strausberg RL
    205. Strempel S
    206. Sturgill D
    207. Sutton G
    208. Sutton GG
    209. Tao W
    210. Teichmann S
    211. Tobari YN
    212. Tomimura Y
    213. Tsolas JM
    214. Valente VL
    215. Venter E
    216. Venter JC
    217. Vicario S
    218. Vieira FG
    219. Vilella AJ
    220. Villasante A
    221. Walenz B
    222. Wang J
    223. Wasserman M
    224. Watts T
    225. Wilson D
    226. Wilson RK
    227. Wing RA
    228. Wolfner MF
    229. Wong A
    230. Wong GK
    231. Wu CI
    232. Wu G
    233. Yamamoto D
    234. Yang HP
    235. Yang SP
    236. Yorke JA
    237. Yoshida K
    238. Zdobnov E
    239. Zhang P
    240. Zhang Y
    241. Zimin AV
    242. Baldwin J
    243. Abdouelleil A
    244. Abdulkadir J
    245. Abebe A
    246. Abera B
    247. Abreu J
    248. Acer SC
    249. Aftuck L
    250. Alexander A
    251. An P
    252. Anderson E
    253. Anderson S
    254. Arachi H
    255. Azer M
    256. Bachantsang P
    257. Barry A
    258. Bayul T
    259. Berlin A
    260. Bessette D
    261. Bloom T
    262. Blye J
    263. Boguslavskiy L
    264. Bonnet C
    265. Boukhgalter B
    266. Bourzgui I
    267. Brown A
    268. Cahill P
    269. Channer S
    270. Cheshatsang Y
    271. Chuda L
    272. Citroen M
    273. Collymore A
    274. Cooke P
    275. Costello M
    276. D'Aco K
    277. Daza R
    278. De Haan G
    279. DeGray S
    280. DeMaso C
    281. Dhargay N
    282. Dooley K
    283. Dooley E
    284. Doricent M
    285. Dorje P
    286. Dorjee K
    287. Dupes A
    288. Elong R
    289. Falk J
    290. Farina A
    291. Faro S
    292. Ferguson D
    293. Fisher S
    294. Foley CD
    295. Franke A
    296. Friedrich D
    297. Gadbois L
    298. Gearin G
    299. Gearin CR
    300. Giannoukos G
    301. Goode T
    302. Graham J
    303. Grandbois E
    304. Grewal S
    305. Gyaltsen K
    306. Hafez N
    307. Hagos B
    308. Hall J
    309. Henson C
    310. Hollinger A
    311. Honan T
    312. Huard MD
    313. Hughes L
    314. Hurhula B
    315. Husby ME
    316. Kamat A
    317. Kanga B
    318. Kashin S
    319. Khazanovich D
    320. Kisner P
    321. Lance K
    322. Lara M
    323. Lee W
    324. Lennon N
    325. Letendre F
    326. LeVine R
    327. Lipovsky A
    328. Liu X
    329. Liu J
    330. Liu S
    331. Lokyitsang T
    332. Lokyitsang Y
    333. Lubonja R
    334. Lui A
    335. MacDonald P
    336. Magnisalis V
    337. Maru K
    338. Matthews C
    339. McCusker W
    340. McDonough S
    341. Mehta T
    342. Meldrim J
    343. Meneus L
    344. Mihai O
    345. Mihalev A
    346. Mihova T
    347. Mittelman R
    348. Mlenga V
    349. Montmayeur A
    350. Mulrain L
    351. Navidi A
    352. Naylor J
    353. Negash T
    354. Nguyen T
    355. Nguyen N
    356. Nicol R
    357. Norbu C
    358. Norbu N
    359. Novod N
    360. O'Neill B
    361. Osman S
    362. Markiewicz E
    363. Oyono OL
    364. Patti C
    365. Phunkhang P
    366. Pierre F
    367. Priest M
    368. Raghuraman S
    369. Rege F
    370. Reyes R
    371. Rise C
    372. Rogov P
    373. Ross K
    374. Ryan E
    375. Settipalli S
    376. Shea T
    377. Sherpa N
    378. Shi L
    379. Shih D
    380. Sparrow T
    381. Spaulding J
    382. Stalker J
    383. Stange-Thomann N
    384. Stavropoulos S
    385. Stone C
    386. Strader C
    387. Tesfaye S
    388. Thomson T
    389. Thoulutsang Y
    390. Thoulutsang D
    391. Topham K
    392. Topping I
    393. Tsamla T
    394. Vassiliev H
    395. Vo A
    396. Wangchuk T
    397. Wangdi T
    398. Weiand M
    399. Wilkinson J
    400. Wilson A
    401. Yadav S
    402. Young G
    403. Yu Q
    404. Zembek L
    405. Zhong D
    406. Zimmer A
    407. Zwirko Z
    408. Jaffe DB
    409. Alvarez P
    410. Brockman W
    411. Butler J
    412. Chin C
    413. Gnerre S
    414. Grabherr M
    415. Kleber M
    416. Mauceli E
    417. MacCallum I

    (2007) Evolution of genes and genomes on the Drosophila phylogeny

    Nature 450:203–218.

    https://doi.org/10.1038/nature06341

    • PubMed
    • Google Scholar
19. 1. Cunningham CW
    2. Omland KE
    3. Oakley TH

    (1998) Reconstructing ancestral character states: a critical reappraisal

    Trends in Ecology & Evolution 13:361–366.

    https://doi.org/10.1016/S0169-5347(98)01382-2

    • PubMed
    • Google Scholar
20. 1. Deutsch D
    2. Pacheco D
    3. Encarnacion-Rivera L
    4. Pereira T
    5. Fathy R
    6. Clemens J
    7. Girardin C
    8. Calhoun A
    9. Ireland E
    10. Burke A
    11. Dorkenwald S
    12. McKellar C
    13. Macrina T
    14. Lu R
    15. Lee K
    16. Kemnitz N
    17. Ih D
    18. Castro M
    19. Halageri A
    20. Jordan C
    21. Silversmith W
    22. Wu J
    23. Seung HS
    24. Murthy M

    (2020) The neural basis for a persistent internal state in Drosophila females

    eLife 9:e59502.

    https://doi.org/10.7554/eLife.59502

    • PubMed
    • Google Scholar
21. 1. Ding Y
    2. Lillvis JL
    3. Cande J
    4. Berman GJ
    5. Arthur BJ
    6. Long X
    7. Xu M
    8. Dickson BJ
    9. Stern DL

    (2019) Neural evolution of Context-Dependent fly song

    Current Biology 29:1089–1099.

    https://doi.org/10.1016/j.cub.2019.02.019

    • PubMed
    • Google Scholar
22. 1. Duistermars BJ
    2. Pfeiffer BD
    3. Hoopfer ED
    4. Anderson DJ

    (2018) A brain module for scalable control of complex, Multi-motor threat displays

    Neuron 100:1474–1490.

    https://doi.org/10.1016/j.neuron.2018.10.027

    • PubMed
    • Google Scholar
23. 1. Ellison CK
    2. Wiley C
    3. Shaw KL

    (2011) The genetics of speciation: genes of small effect underlie sexual isolation in the hawaiian cricket Laupala

    Journal of Evolutionary Biology 24:1110–1119.

    https://doi.org/10.1111/j.1420-9101.2011.02244.x

    • PubMed
    • Google Scholar
24. 1. Felsenstein J

    (1985) Phylogenies and the comparative method

    The American Naturalist 125:1–15.

    https://doi.org/10.1086/284325

    • Google Scholar
25. 1. Felsenstein J

    (2005) Using the quantitative genetic threshold model for inferences between and within species

    Philosophical Transactions of the Royal Society B: Biological Sciences 360:1427–1434.

    https://doi.org/10.1098/rstb.2005.1669

    • Google Scholar
26. 1. Gelman A
    2. Rubin DB

    (1992) Inference from iterative simulation using multiple sequences

    Statistical Science 7:457–472.

    https://doi.org/10.1214/ss/1177011136

    • Google Scholar
27. 1. Gleason JM
    2. Ritchie MG

    (2004) Do quantitative trait loci (QTL) for a courtship song difference between Drosophila simulans and D. sechellia coincide with candidate genes and intraspecific QTL?

    Genetics 166:1303–1311.

    https://doi.org/10.1534/genetics.166.3.1303

    • PubMed
    • Google Scholar
28. Book

    1. Goodfellow I
    2. Bengio Y
    3. Courville A

    (2016)

    Deep Learning

    Cambridge, USA: MIT press.

    • Google Scholar
29. 1. Hadfield JD

    (2010) MCMC Methods for Multi-Response Generalized Linear Mixed Models: The MCMCglmm R Package

    Journal of Statistical Software 33:1–22.

    https://doi.org/10.18637/jss.v033.i02

    • Google Scholar
30. 1. Hadfield JD
    2. Nakagawa S

    (2010) General quantitative genetic methods for comparative biology: phylogenies, taxonomies and multi-trait models for continuous and categorical characters

    Journal of Evolutionary Biology 23:494–508.

    https://doi.org/10.1111/j.1420-9101.2009.01915.x

    • PubMed
    • Google Scholar
31. 1. Hansen TF
    2. Martins EP

    (1996) Translating between microevolutionary process and macroevolutionary patterns: the correlation structure of interspecific data

    Evolution 50:1404–1417.

    https://doi.org/10.1111/j.1558-5646.1996.tb03914.x

    • Google Scholar
32. 1. Hoopfer ED

    (2016) Neural control of aggression in Drosophila

    Current Opinion in Neurobiology 38:109–118.

    https://doi.org/10.1016/j.conb.2016.04.007

    • PubMed
    • Google Scholar
33. 1. Hu CK
    2. Hoekstra HE

    (2017) Peromyscus burrowing: a model system for behavioral evolution

    Seminars in Cell & Developmental Biology 61:107–114.

    https://doi.org/10.1016/j.semcdb.2016.08.001

    • PubMed
    • Google Scholar
34. 1. Levine M
    2. Davidson EH

    (2005) Gene regulatory networks for development

    PNAS 102:4936–4942.

    https://doi.org/10.1073/pnas.0408031102

    • PubMed
    • Google Scholar
35. 1. Ligon RA
    2. Diaz CD
    3. Morano JL
    4. Troscianko J
    5. Stevens M
    6. Moskeland A
    7. Laman TG
    8. Scholes E

    (2018) Evolution of correlated complexity in the radically different courtship signals of birds-of-paradise

    PLOS Biology 16:e2006962.

    https://doi.org/10.1371/journal.pbio.2006962

    • PubMed
    • Google Scholar
36. 1. Lorenz KZ

    (1958) The evolution of behavior

    Scientific American 199:67–78.

    https://doi.org/10.1038/scientificamerican1258-67

    • PubMed
    • Google Scholar
37. Book

    1. Martins EP

    (1996)

    Phylogenies and the Comparative Method in Animal Behavior

    Oxford: Oxford University Press.

    • Google Scholar
38. 1. Martins EP
    2. Hansen TF

    (1997) Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data

    The American Naturalist 149:646–667.

    https://doi.org/10.1086/286013

    • Google Scholar
39. 1. Massey JH
    2. Chung D
    3. Siwanowicz I
    4. Stern DL
    5. Wittkopp PJ

    (2019) The yellow gene influences Drosophila male mating success through sex comb melanization

    eLife 8:e49388.

    https://doi.org/10.7554/eLife.49388

    • PubMed
    • Google Scholar
40. 1. Mathis MW
    2. Mathis A

    (2020) Deep learning tools for the measurement of animal behavior in neuroscience

    Current Opinion in Neurobiology 60:1–11.

    https://doi.org/10.1016/j.conb.2019.10.008

    • PubMed
    • Google Scholar
41. 1. Meyer F

    (1994) Topographic distance and watershed lines

    Signal Processing 38:113–125.

    https://doi.org/10.1016/0165-1684(94)90060-4

    • Google Scholar
42. 1. Obbard DJ
    2. Maclennan J
    3. Kim KW
    4. Rambaut A
    5. O'Grady PM
    6. Jiggins FM

    (2012) Estimating divergence dates and substitution rates in the Drosophila phylogeny

    Molecular Biology and Evolution 29:3459–3473.

    https://doi.org/10.1093/molbev/mss150

    • PubMed
    • Google Scholar
43. 1. O’Meara BC

    (2012) Evolutionary inferences from phylogenies: a review of methods annual review of ecology

    Evolution 43:267–285.

    https://doi.org/10.1146/annurev-ecolsys-110411-160331

    • Google Scholar
44. 1. Pereira TD
    2. Aldarondo DE
    3. Willmore L
    4. Kislin M
    5. Wang SS
    6. Murthy M
    7. Shaevitz JW

    (2019) Fast animal pose estimation using deep neural networks

    Nature Methods 16:117–125.

    https://doi.org/10.1038/s41592-018-0234-5

    • PubMed
    • Google Scholar
45. 1. Rand WM

    (1971) Objective criteria for the evaluation of clustering methods

    Journal of the American Statistical Association 66:846–850.

    https://doi.org/10.1080/01621459.1971.10482356

    • Google Scholar
46. 1. Royer-Carenzi M
    2. Didier G

    (2016) A comparison of ancestral state reconstruction methods for quantitative characters

    Journal of Theoretical Biology 404:126–142.

    https://doi.org/10.1016/j.jtbi.2016.05.029

    • PubMed
    • Google Scholar
47. 1. Seeds AM
    2. Ravbar P
    3. Chung P
    4. Hampel S
    5. Midgley FM
    6. Mensh BD
    7. Simpson JH

    (2014) A suppression hierarchy among competing motor programs drives sequential grooming in Drosophila

    eLife 3:e02951.

    https://doi.org/10.7554/eLife.02951

    • PubMed
    • Google Scholar
48. 1. Seetharam AS
    2. Stuart GW

    (2013) Whole genome phylogeny for 21 Drosophila species using predicted 2b-RAD fragments

    PeerJ 1:e226.

    https://doi.org/10.7717/peerj.226

    • PubMed
    • Google Scholar
49. 1. Shaw KL
    2. Lesnick SC

    (2009) Genomic linkage of male song and female acoustic preference QTL underlying a rapid species radiation

    PNAS 106:9737–9742.

    https://doi.org/10.1073/pnas.0900229106

    • PubMed
    • Google Scholar
50. 1. Shubin N
    2. Tabin C
    3. Carroll S

    (2009) Deep homology and the origins of evolutionary novelty

    Nature 457:818–823.

    https://doi.org/10.1038/nature07891

    • PubMed
    • Google Scholar
51. 1. Slonim N
    2. Atwal GS
    3. Tkacik G
    4. Bialek W

    (2005) Information-based clustering

    PNAS 102:18297–18302.

    https://doi.org/10.1073/pnas.0507432102

    • PubMed
    • Google Scholar
52. 1. Spiegelhalter DJ
    2. Best NG
    3. Carlin BP
    4. van der Linde A

    (2002) Bayesian measures of model complexity and fit

    Journal of the Royal Statistical Society: Series B 64:583–639.

    https://doi.org/10.1111/1467-9868.00353

    • Google Scholar
53. 1. Stern DL
    2. Crocker J
    3. Ding Y
    4. Frankel N
    5. Kappes G
    6. Kim E
    7. Kuzmickas R
    8. Lemire A
    9. Mast JD
    10. Picard S

    (2017) Genetic and transgenic reagents for Drosophila simulans, D. mauritiana, D. yakuba, D. santomea, and D. virilis

    G3: Genes, Genomes, Genetics 7:1339–1347.

    https://doi.org/10.1534/g3.116.038885

    • PubMed
    • Google Scholar
54. 1. Stern DL
    2. Frankel N

    (2013) The structure and evolution of Cis -regulatory regions: the shavenbaby story

    Philosophical Transactions of the Royal Society B: Biological Sciences 368:0028.

    https://doi.org/10.1098/rstb.2013.0028

    • Google Scholar
55. 1. Strouse DJ
    2. Schwab DJ

    (2017) The deterministic information bottleneck

    Neural Computation 29:1611–1630.

    https://doi.org/10.1162/NECO_a_00961

    • PubMed
    • Google Scholar
56. 1. Tajima F

    (1993) Simple methods for testing the molecular evolutionary clock hypothesis

    Genetics 135:599–607.

    https://doi.org/10.1093/genetics/135.2.599

    • PubMed
    • Google Scholar
57. Book

    1. Tinbergen N

    (1951)

    The Study of Instinct

    Oxford: Oxford University Press.

    • Google Scholar
58. Conference

    1. Tishby N
    2. Pereira FC
    3. Bialek W

    (1999)

    The information bottleneck method

    Proceedings of the 37th Annual Allerton Conference on Communication, Control and Computing Urbana-Champaign. pp. 368–377.

    • Google Scholar
59. 1. van der Maaten L
    2. Hinton G

    (2008)

    Visualizing data using t-SNE

    Journal of Machine Learning Research 9:2579–2605.

    • Google Scholar
60. 1. Weber JN
    2. Peterson BK
    3. Hoekstra HE

    (2013) Discrete genetic modules are responsible for complex burrow evolution in Peromyscus mice

    Nature 493:402–405.

    https://doi.org/10.1038/nature11816

    • PubMed
    • Google Scholar
61. Book

    1. West-Eberhard MJ

    (2003)

    Developmental Plasticity and Evolution

    Oxford: Oxford University Press.

    • Google Scholar
62. 1. Williams TM
    2. Carroll SB

    (2009) Genetic and molecular insights into the development and evolution of sexual dimorphism

    Nature Reviews Genetics 10:797–804.

    https://doi.org/10.1038/nrg2687

    • PubMed
    • Google Scholar
63. 1. Yamamoto D
    2. Ishikawa Y

    (2013) Genetic and neural bases for species-specific behavior in Drosophila species

    Journal of Neurogenetics 27:130–142.

    https://doi.org/10.3109/01677063.2013.800060

    • PubMed
    • Google Scholar
64. Book

    1. Yang Z

    (2006)

    Computational Molecular Evolution

    Oxford: Oxford University Press.

    • Google Scholar

Acceptance summary:

Different animal species exhibit distinct behavioral repertoires, even in cases where common ancestors are relatively recent. How do repertoires acquire species-specificity during evolution? This manuscript provides new methods to examine behavioral patterns of different Drosophila species, and even to reconstruct potentially ancestral behavioral modes. Consistent with past work, the authors find that behaviors are inherited in clusters. Overall, this stimulating paper presents an exciting method for the unbiased study of behavioral evolution based on large datasets that are increasingly common in animal behavior, which differs from previous work focusing on specific traits.

Decision letter after peer review:

Thank you for submitting your article "A framework for studying behavioral evolution by reconstructing ancestral repertoires" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Christian Rutz as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Iain D Couzin (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is “in revision at eLife”. Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

Summary:

In this paper, the authors build upon their previous (ground-breaking) work in which automated methodology was employed to quantify the structure of behaviour and its regulation in fruit flies. Specifically, here, they address whether such methodology can be insightful in the study of how behaviour evolves, again using fruit flies as convenient and powerful model species (their phylogenetic relationships are well-understood, they do well in the lab, and it has been previously shown how well the behaviour of these species can be decomposed using automated methodology). There are a host of reasons why such an analysis is of value, not least because revealing the hidden structure of behavioural regulation has the potential to be informative regarding ancestral behavioural repertoires, and thus how it evolves. Overall, there is consensus that this is a thought-provoking paper that provides a valuable starting point in the analysis of behavioural evolution via detailed quantitative ethology. With added consideration of the factors discussed below this work makes a very helpful and novel contribution.

Essential revisions:

1. Consider behavioural transitions

All three reviewers were confused as to why transitions between states were not considered. Given that the authors have previously shown their capacity to reconstruct the hierarchical nested structure of behavioural states, typically represented as a network with transition probabilities between states, which would seemingly be extremely informative about the evolution of behavioural modules, all reviewers were left wondering why they did not consider comparisons of these networks -- especially given the array of network analysis tools available. Organisms, as they are well aware, don't just jump randomly from one behaviour to another, but there is a hierarchical organisation by which one behaviour influences the probability that another is elicited, and so on. Why did the authors not do this? There may be very good reasons, but it will be important for the reader to be informed. This would seem to be the most direct way to address how behaviour is structured. It is not clear to the reviewers exactly how much extra work this would be -- or how the scope of the paper would be required to change if transitions were considered. Thus, while formal inclusion of transitions into the framework is strongly preferred for revision, it is not absolutely required. At the minimum, the authors need to clarify exactly why transitions were not considered as well as what may be gained from such an approach in future work.

2. Confirm utility of the approach, its discoveries, and how it compares to more focused inquiries into the evolution of behaviour

The work suffers from a lack of clarity over why/how the method is superior to the focused trait approach. The fact that behaviours co-evolve in suites is by itself not especially novel. The authors should bring to the forefront, if possible, exactly how and why this approach is superior to others. Ultimately, the paper is an argument for the utility of studying repertoires rather than specific behaviours. The main thrust the paper tries to make is that we learn something new by looking at the evolution of repertoires of behaviour rather than focused analyses of specific behaviours (really they are looking at time budgets and perhaps that's a bit of a disconnect from their message as well). But they put forth in Figure 3 a reconstruction that seems to generate a nonsense behavioural repertoire for the ancestral species. So the authors should address the validity of using the specific phylogenetic methods to reconstruct behaviours based on the repertoires. They can address this point in two ways:

2.1. Generate simulations of repertoire evolution on the known Drosophila phylogeny and then take the endpoint repertoires they have evolved and assess how well the reconstruction of the known ancestral state works. People do these sorts of studies regularly to assess phylogenetic reconstruction approaches and methods (e.g., Royer-Carenzi and Didier 2016, J. Theoretical Biol.). This might be done with a real or arbitrary repertoire as the starting point for the simulation. This is crucial since it appears the rest of their findings hinge on the ancestral reconstruction.

2.2. As a complement to 2.1, the authors might assess the reconstruction of the repertoires in comparison to ancestral trait reconstructions for specific behavioural traits that might be more conventionally measured, assuming this can be extracted from the existing data. Part of the argument the authors make is that by considering the repertoire we achieve a more complete understanding of behavioral evolution. The approach seems very promising, though it is difficult to assess how it compares to other approaches for measuring behaviour. However good or bad the reconstructions of ancestral behaviours are for the repertoires, how do they compare to a focused trait approach? Perhaps better? Perhaps worse? Unless a comparison is done it's hard to say. More broadly, the paper argues for the utility of this approach to studying behaviour but it is not clear how it compares to analyses that focus on specific elements of behaviour rather than time budget vectors.

3. Make the paper more accessible for a general audience

A major issue with this paper was its inaccessibility, which could, unfortunately, ultimately reduce its impact. Most readers are by now familiar with past work using clustering methods (e.g., Cande et al., 2018; Berman et al., 2014) -- but many of the analyses in this paper are novel, difficult to penetrate, and poorly introduced for a general audience. The authors should take a fine-grained comb to this paper to make each analysis as accessible as possible. Simple fixes could include more plain language around what each analysis is testing/achieving, as well as consistency in terminology. Specifically, the final sentence of each paragraph could state in plain terms what the preceding analysis has just demonstrated, ruled-in, or ruled out.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "A framework for studying behavioral evolution by reconstructing ancestral repertoires" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Jesse H Goldberg as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Christian Rutz as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Iain D Couzin (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this decision letter to help you prepare a revised submission.

Please address continued concerns from Reviewer #3, and please attempt to clarify key sections of the text as recommended by Reviewer #1.

1) Reviewer #3: The paper has improved in many respects but the revision failed to deal with what I see as the most glaring issue (and which I had raised in my previous review).

All the extant species show biases in behaviors as shown by the hotter/browner colors in their behavioral maps. This suggests that we should expect that any species of Drosophila, extant or extinct, would be highly likely to have biases behaviors as well. However, the ancestrally reconstructed behavioral repertoire is totally washed out – there is no clear biases in behavior in the ancestral graph that is estimated in figure 3. This suggests to me that the method potentially has flaws.

I had suggested that the authors start with a PDF from one of their species and then simulate its evolution along the known Drosophila phylogeny and then try to reconstruct that repertoire using ancestral state reconstruction. If that process ends up showing a similarly 'flat' PDF where everything is similarly low probability (i.e. a totally blue PDF as shown in figure 3) then there is some cause for concern. If it were to retrieve a PDF similar in structure to the one that initiated the simulation, then that would be very reassuring.

The author responses seem to say that is what they have done in using the MCMC method but where the ancestral reconstruction analysis of a simulated evolutionary trajectory is to be found in the paper is unclear to me. I searched 'simulat' in the text of the manuscript and cannot see where trait evolution was simulated and then used as the starting values to try to evaluate how closely the ancestral reconstruction gets to the known ancestral starting point.

There is a lot in this paper, so though I feel this issue is important to address, does not necessarily preclude it from being published in my view. Though I would like to see the authors address the 'flatness' and thus seemingly unrealistic (it would seem to me at least) ancestral state that they present in figure 3.

2) Reviewer #1: The writing is slightly improved, though some sections could still be improved by concluding paragraphs not just with a mathematical result but also what it means in plain terms. E.g. Lines 264-7 conclude with the finding that the covariance matrix has a 'far from random modular structure;' this does not explain how this result relates to the title of the section (Individual Variability and long timescale correlations).

3) Please note that eLife has recently adopted the STRANGE guidelines for animal behaviour research:

https://doi.org/10.1038/d41586-020-01751-5

https://reviewer.elifesciences.org/author-guide/journal-policies

In your revisions, please consider scope for sampling biases in your study and how these may limit the generalisability of your findings, and make declarations as necessary. A few sentences in the section "Data collection" (lines 462-470) may suffice.

Essential revisions:

1. Consider behavioural transitions

All three reviewers were confused as to why transitions between states were not considered. Given that the authors have previously shown their capacity to reconstruct the hierarchical nested structure of behavioural states, typically represented as a network with transition probabilities between states, which would seemingly be extremely informative about the evolution of behavioural modules, all reviewers were left wondering why they did not consider comparisons of these networks -- especially given the array of network analysis tools available. Organisms, as they are well aware, don't just jump randomly from one behaviour to another, but there is a hierarchical organisation by which one behaviour influences the probability that another is elicited, and so on. Why did the authors not do this? There may be very good reasons, but it will be important for the reader to be informed. This would seem to be the most direct way to address how behaviour is structured. It is not clear to the reviewers exactly how much extra work this would be -- or how the scope of the paper would be required to change if transitions were considered. Thus, while formal inclusion of transitions into the framework is strongly preferred for revision, it is not absolutely required. At the minimum, the authors need to clarify exactly why transitions were not considered as well as what may be gained from such an approach in future work.

We agree that studying the evolution of the transition structure between behaviors is an important subject, and we are currently developing methods to study it. The primary difficulty here is the number of parameters involved. Because we are fitting two covariance matrices and an ancestral mean, we need to fit at least (N + 1) ∗N + N parameters, where N is the number of traits being reconstructed. Currently, we study a system that has N = 134 traits, leading to 18,224 parameters that we must fit with the GLMM. If we were to add the transition probabilities (or transition rates), this would require another ≈N² parameters. Thus, instead of fitting a system with ≈ 18,000 parameters, our model would now have 327,284,280 parameters. While we would likely be able to perform the computation, we would be unlikely to believe the results at this point because our data set size is insufficient (the amount of time needed to accurately sample a transition matrix is a power of two longer than the amount of time needed to accurately sample the behavioral frequencies).

In addition, these traits (the transition matrix values) have different units than the behaviors, and the transition probabilities mathematically depend on the behavioral frequencies, since the first eigenvalue of the transition matrix is (by definition) proportional to the average behavioral frequency usage. Thus, it becomes necessary to add additional regularizing or normalizing factors, further complicating the analysis in a manner that is beyond the scope of what we aimed to demonstrate here.

Despite these difficulties, although transitions are not included in our actual model, one of the main points that we have highlighted is the similarity between the hierarchical structure that explains long timescale transitions, which we demonstrated previously, and the variability we observe here between individuals of the same species.

Again, though, this is a problem that we find extremely interesting, and we are actively working on exploring how to solve these technical and representational questions. We had some text describing this as a future direction in our Discussion section, but we have now expanded this section to further point out the importance of this question and some of the current technical challenges surrounding the study of the evolution of behavioral transitions (lines 428-441).

2. Confirm utility of the approach, its discoveries, and how it compares to more focused inquiries into the evolution of behaviour

The work suffers from a lack of clarity over why/how the method is superior to the focused trait approach. The fact that behaviours co-evolve in suites is by itself not especially novel. The authors should bring to the forefront, if possible, exactly how and why this approach is superior to others. Ultimately, the paper is an argument for the utility of studying repertoires rather than specific behaviours. The main thrust the paper tries to make is that we learn something new by looking at the evolution of repertoires of behaviour rather than focused analyses of specific behaviours (really they are looking at time budgets and perhaps that's a bit of a disconnect from their message as well). But they put forth in Figure 3 a reconstruction that seems to generate a nonsense behavioural repertoire for the ancestral species. So the authors should address the validity of using the specific phylogenetic methods to reconstruct behaviours based on the repertoires. They can address this point in two ways:

2.1. Generate simulations of repertoire evolution on the known Drosophila phylogeny and then take the endpoint repertoires they have evolved and assess how well the reconstruction of the known ancestral state works. People do these sorts of studies regularly to assess phylogenetic reconstruction approaches and methods (e.g., Royer-Carenzi and Didier 2016, J. Theoretical Biol.). This might be done with a real or arbitrary repertoire as the starting point for the simulation. This is crucial since it appears the rest of their findings hinge on the ancestral reconstruction.

This process is precisely what we are doing with our MCMC method (which was originally developed by Hadfield, as cited in the text). Specifically, given a known tree structure, we select a common ancestral mean behavior and two behavioral covariance matrix such that the likelihood of observing the measured behavioral distribution (including the covariances between the behaviors) is maximized. In addition to showing that our model has converged (Figure 3, Supplement 1), we also show that the predicted behavioral means at the endpoint species match the data (Figure 3, Supplement 2).

2.2. As a complement to 2.1, the authors might assess the reconstruction of the repertoires in comparison to ancestral trait reconstructions for specific behavioural traits that might be more conventionally measured, assuming this can be extracted from the existing data. Part of the argument the authors make is that by considering the repertoire we achieve a more complete understanding of behavioral evolution. The approach seems very promising, though it is difficult to assess how it compares to other approaches for measuring behaviour. However good or bad the reconstructions of ancestral behaviours are for the repertoires, how do they compare to a focused trait approach? Perhaps better? Perhaps worse? Unless a comparison is done it's hard to say. More broadly, the paper argues for the utility of this approach to studying behaviour but it is not clear how it compares to analyses that focus on specific elements of behaviour rather than time budget vectors.

To address these points, we have added an extra section in the main text, as well as in Materials and methods, that shows a comparison between our model and a simpler model where behaviors evolve independently (lines 214-229 and 491-524, as well as Figure 3, Supplement 3). In order to compare the performance of models with different complexity, we need a measure that not only takes into account how well the data is fit (i.e., maximizes likelihood), but that also penalizes the addition of extra parameters in more complex models (in this case, the off-diagonal terms of the two behavioral covariance matrices). For hierarchical models fit using MCMC methods (such as ours), the deviance information criterion (DIC) is the standard approach for model selection (Spiegelhalter et al., 2002). The DIC is a generalization of the Akaike criterion (see Materials and methods for additional details), and it is proportional to the negative log-likelihood evaluated at the mean of the posterior parameters, plus a penalty term proportional to the effective number of parameters in the model. The values of DIC for our model are much smaller than the models with a focused analysis, indicating that modeling covariance between behaviors explains the observed data much better than treating all behaviors independently. The implications of this result have been further expanded and discussed in the text, and we thank the reviewers for suggesting this validations approach.

3. Make the paper more accessible for a general audience

A major issue with this paper was its inaccessibility, which could, unfortunately, ultimately reduce its impact. Most readers are by now familiar with past work using clustering methods (e.g., Cande et al., 2018; Berman et al., 2014) -- but many of the analyses in this paper are novel, difficult to penetrate, and poorly introduced for a general audience. The authors should take a fine-grained comb to this paper to make each analysis as accessible as possible. Simple fixes could include more plain language around what each analysis is testing/achieving, as well as consistency in terminology. Specifically, the final sentence of each paragraph could state in plain terms what the preceding analysis has just demonstrated, ruled-in, or ruled out.

We agree with the reviewers that many aspects of the original manuscript text could have benefited from additional prose explaining the rationale and conclusions of the analysis. Following these comments, we have now included many additional clarifying comments to make the text more accessible.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Please address continued concerns from Reviewer #3, and please attempt to clarify key sections of the text as recommended by Reviewer #1.

(1) Reviewer #3: The paper has improved in many respects but the revision failed to deal with what I see as the most glaring issue (and which I had raised in my previous review).

All the extant species show biases in behaviors as shown by the hotter/browner colors in their behavioral maps. This suggests that we should expect that any species of Drosophila, extant or extinct, would be highly likely to have biases behaviors as well. However, the ancestrally reconstructed behavioral repertoire is totally washed out – there is no clear biases in behavior in the ancestral graph that is estimated in figure 3. This suggests to me that the method potentially has flaws.

I had suggested that the authors start with a PDF from one of their species and then simulate its evolution along the known Drosophila phylogeny and then try to reconstruct that repertoire using ancestral state reconstruction. If that process ends up showing a similarly 'flat' PDF where everything is similarly low probability (i.e. a totally blue PDF as shown in figure 3) then there is some cause for concern. If it were to retrieve a PDF similar in structure to the one that initiated the simulation, then that would be very reassuring.

The author responses seem to say that is what they have done in using the MCMC method but where the ancestral reconstruction analysis of a simulated evolutionary trajectory is to be found in the paper is unclear to me. I searched 'simulat' in the text of the manuscript and cannot see where trait evolution was simulated and then used as the starting values to try to evaluate how closely the ancestral reconstruction gets to the known ancestral starting point.

There is a lot in this paper, so though I feel this issue is important to address, does not necessarily preclude it from being published in my view. Though I would like to see the authors address the 'flatness' and thus seemingly unrealistic (it would seem to me at least) ancestral state that they present in figure 3.

We agree that the way in which we presented Figure 3 provides the impression that the ancestral repertoire is “washed out.” This is an effect of the way in which we plotted our data, and we thank the reviewer for pointing out this less-than-ideal visualization. The new version of the figure shows the ancestral distribution in non-logarithmic values, showing that there is indeed structure in this repertoire. We agree that it would be worrying if the distribution was indeed flat, but as can be seen, this isn’t the case here. All of the brown/hotter colors in the subsequent repertoires to the right are differences of the logs, which better shows the behavioral alterations along the phylogeny, as they are often subtle.

As to the proposed simulation, due to the reviewer’s excellent suggestion in the previous round of reviews, we added Figure 3—figure supplement 2, which shows that starting from the inferred ancestral state, the model’s mean predictions are in agreement with the data, providing a good double-check of our model’s self-consistency. In other words, when starting from the endpoints, we reconstruct the ancestral state (this is how the model is fit in the first place), and starting from the fit ancestral state, the model predicts, on average, the correct endpoint states. In addition, Figure 3 Figure supplement 3 (also emerging from the reviewers’ suggestions), shows that the model predicts the behavioral covariances well. We did not need to do simulations, per se, to show these results, since these values can be directly numerically calculated from the model (hence the lack of the word ’simulation’ connected to these efforts in the manuscript).

(2) Reviewer #1: The writing is slightly improved, though some sections could still be improved by concluding paragraphs not just with a mathematical result but also what it means in plain terms. E.g. Lines 264-7 conclude with the finding that the covariance matrix has a 'far from random modular structure;' this does not explain how this result relates to the title of the section (Individual Variability and long timescale correlations).

We thank the reviewer for their comments, and we have adapted several passages in the text to increase the clarity of our methodology and to make connections between the mathematical results and the according biological significance (see the latex-diff document for all changes).

(3) Please note that eLife has recently adopted the STRANGE guidelines for animal behaviour research:

https://doi.org/10.1038/d41586-020-01751-5

https://reviewer.elifesciences.org/author-guide/journal-policies

In your revisions, please consider scope for sampling biases in your study and how these may limit the generalisability of your findings, and make declarations as necessary. A few sentences in the section "Data collection" (lines 462-470) may suffice.

We thank the reviewers for pointing this out (even as an eLife reviewing editor, GJB was embarrassingly unaware of these guidelines). Although detailed information about how the flies were handled/housed was previously published in Cande et al., eLife, 2018 (and was pointed to accordingly via a citation), we added relevant details here to the Materials and methods section to make the manuscript more self-contained and to make our experimental details more clear.

As to sampling biases, we have provided a detailed list of the specific strains used in our experimental – all of which are readily available through the UCSD or Bloomington fly stocks. While we have no evidence that our precise choice of species/strains within the D. melanogaster species subgroup should generate bias, we do acknowledge in the Discussion section that our phylogentic reconstruction is likely under-constrained, and that adding more species/strains will “place more constraints on the evolutionary dynamics, likely resulting in less variance in the ancestral state estimations and potentially adding more structure to the relatively low rank (i.e., highly modular) covariance matrices. Additionally, further work is required to determine the balance between sampling within and between strains and species that optimizes estimates of evolutionary dynamics.”

Author details

1. Damián G Hernández

   1. Department of Physics, Emory University, Atlanta, United States
   2. Department of Medical Physics, Centro Atómico Bariloche and Instituto Balseiro, Bariloche, Argentina

   Contribution

   Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

   Contributed equally with

   Catalina Rivera

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-8995-7495

2. Catalina Rivera

   Department of Physics, Emory University, Atlanta, United States

   Contribution

   Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

   Contributed equally with

   Damián G Hernández

   Competing interests

   No competing interests declared

3. Jessica Cande

   Janelia Research Campus, Howard Hughes Medical Institute, Ashburn, United States

   Contribution

   Conceptualization, Data curation, Investigation, Methodology, Writing - review and editing

   Competing interests

   No competing interests declared

4. Baohua Zhou

   1. Department of Physics, Emory University, Atlanta, United States
   2. Department of Molecular, Cellular and Developmental Biology, Yale University, New Haven, United States

   Contribution

   Conceptualization, Software, Investigation, Methodology, Writing - review and editing

   Competing interests

   No competing interests declared

5. David L Stern

   Janelia Research Campus, Howard Hughes Medical Institute, Ashburn, United States

   Contribution

   Conceptualization, Resources, Funding acquisition, Investigation, Visualization, Methodology, Writing - review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-1847-6483

6. Gordon J Berman

   1. Department of Physics, Emory University, Atlanta, United States
   2. Department of Biology, Emory University, Atlanta, United States

   Contribution

   Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

   For correspondence

   gordon.berman@emory.edu

   Competing interests

   Reviewing editor, eLife

   "This ORCID iD identifies the author of this article:" 0000-0003-3588-7820

• 2,678

  Page views

• 333

  Downloads

• 6

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.