# Contemporary Ecological Interactions Improve Models of Past Trait Evolution

Contemporary Ecological Interactions Improve Models of Past Trait Evolution Abstract Despite the fact that natural selection underlies both traits and interactions, evolutionary models often neglect that ecological interactions may, and in many cases do, influence the evolution of traits. Herein, we explore the interdependence of ecological interactions and functional traits in the pollination associations of hawkmoths and flowering plants. Specifically, we develop an adaptation of the Ornstein–Uhlenbeck model of trait evolution that allows us to study the influence of plant corolla depth and observed hawkmoth–plant interactions on the evolution of hawkmoth proboscis length. Across diverse modelling scenarios, we find that the inclusion of contemporary interactions can provide a better description of trait evolution than the null expectation. Moreover, we show that the pollination interactions provide more-likely models of hawkmoth trait evolution when interactions are considered at increasingly fine-scale groups of hawkmoths. Finally, we demonstrate how the results of best-fit modeling approaches can implicitly support the association between interactions and trait evolution that our method explicitly examines. In showing that contemporary interactions can provide insight into the historical evolution of hawkmoth proboscis length, we demonstrate the clear utility of incorporating additional ecological information to models designed to study past trait evolution. Macroevolution, mutualism, Ornstein-Uhlenbeck, pollination, Sphingidae For interactions to occur between species, both of the species involved must possess traits conducive to their occurrence. In the case of some ecological interactions, such as reproductive mutualisms, the interdependence of species across the interaction is long thought to lead to adaptation “in the most perfect manner” (Darwin 1859) such that the morphological traits of mutualists become tightly coupled (Thompson 1994). Indeed, morphological trait matching between mutualists is realized across a variety of qualitatively different ecological systems and taxa: from hummingbird-Heliconia pollination and bird-palm seed dispersal to nitrogen fixation mutualisms (Stiles 1975; Brouat et al. 2001; Miklashevichs et al. 2001; Stang et al. 2009; Galetti et al. 2013). Widespread morphological matching implies that in the most closely-coupled associations—whether at the scale of entire guilds (à la interaction syndromes; Janzen and Martin 1982; Fenster et al. 2004) or that of specific species pairs (Darwin 1862)—the morphology of one species may in fact be informative enough to infer the relevant traits and identities of their interaction partners. The evolutionary trajectories of traits on either side of an interaction are in fact influenced by multiple selection pressures, even in purported tightly coupled interactions (Schluter et al. 1991; Kessler and Halitschke 2009). In particular, this may come about since one-to-one interaction specialization is exceedingly rare; species long thought to be hyper-specialized often have more than one interaction partner in reality (Feinsinger et al. 1986; Pellmyr 2003; Machado et al. 2005), and therefore, experience more than one biotic selection pressure on the trait, or traits, relevant to that interaction. Additionally, trait evolution can be influenced by environmental factors, energy budget, and ancestry, amongst others (Hansen 1997; Slabbekoorn and Smith 2002; Niven and Laughlin 2008). Putting this all together, it remains open to debate whether and when an extant species’ interactions and interaction partners’ traits are informative with regard to past trait evolution. The pollination associations between hawkmoths and flowering plants are an ideal system in which to examine the utility of ecological interactions for understanding trait evolution given that close coupling of hawkmoth proboscis length and flower corolla depth is often observed in these interactions (Nilsson 1988; Haverkamp et al. 2016; Sazatornil et al. 2016). That said, morphological matching across interactions in this system is most typically examined from an ecological perspective—and with function as the primary focus—whereas the role of evolutionary processes in building congruence is potentially overlooked (Stang et al. 2009; Dehling et al. 2014; Sazatornil et al. 2016). Consequently, one approach to better explore this question would be to investigate whether a hawkmoth’s interaction partners and those partners’ corolla depths prove informative when modeling the evolution of hawkmoth proboscis length. Contemporary models of trait evolution tend to be adaptations of the random-walk Brownian motion and adaptive Ornstein–Uhlenbeck models (Felsenstein 1985; Hansen 1997; Butler and King 2004; Hansen et al. 2008; Ingram and Mahler 2013; Uyeda and Harmon 2014). Many of these models allow researchers to test hypotheses about how trait evolution relates to species’ ecology, from environmental factors (Hansen et al. 2008; Labra et al. 2009; Muschick et al. 2012; Mahler et al. 2013; Moen et al. 2015) to behavior (Scales et al. 2009; Monteiro and Nogueira 2011). Additionally, more-recent approaches have broadened the scope of these methods by allowing some consideration of species interactions (Nuismer et al. 2013; Nuismer and Harmon 2015; Drury et al. 2016; Manceau et al. 2016; Bartoszek et al. 2017), introducing Bayesian approaches that allow additional data to inform the model (Uyeda and Harmon 2014; Kostikova et al. 2016), and by allowing variation of model parameters across the tree (Beaulieu et al. 2012; Ingram and Mahler 2013; Uyeda and Harmon 2014). Herein, we explicitly incorporate empirical data on these interactions—the biological milieu in which proboscides evolve—into models of hawkmoth trait evolution. This allows us to examine both the utility of contemporary interactions in informing the evolution of this trait as well as delineating the ecological groupings, or selective regimes, where they are the most informative. In taking a hypothesis-driven approach to both the groupings of species and the values of parameters within the model, we explore the degree to which flowering plants and their contemporary floral traits can describe the proboscis-length evolution of their hawkmoth pollinators. If contemporary interactions between plants and hawkmoths do in fact capture some degree of proboscis-length evolution, our minimum expectation is that models with explicit incorporation of hawkmoth interactions will be better supported than an interaction-free model of trait evolution. Moreover, different hypothesis-based models will almost always differ in their degree of support—based on their underlying assumptions—and thus should provide additional information about the evolution of these traits. Finally, we examine the results of an uninformed, best-fit modeling approach to assess the importance of our initial hypothesis—that ecological interactions have a discernible impact on trait evolution—in these pollination systems. Consequently, we demonstrate that there is unambiguous evidence that contemporary interactions can shed light on hawkmoth trait evolution. We show that this evidence can be found in the comparisons of specific species-interaction-based hypotheses and when comparing these scenarios to unconstrained approaches to trait evolution. Methods Hawkmoth Dataset For the purposes of this study, we put together a dataset that encompasses $$75$$ hawkmoth species from $$25$$ genera and $$92$$ plant species from $$60$$ genera based on the findings of two recent studies (Kawahara and Barber 2015; Sazatornil et al. 2016). To use these data to explore the utility of ecological interactions for modeling trait evolution, we first compiled hawkmoth–plant interactions and relevant traits from four of the quantitative pollination networks presented by (Sazatornil et al., 2016). We combined the data from these four networks so that we were left with one set of realized interactions per hawkmoth species and one value for the length of a hawkmoth’s proboscis—when a hawkmoth occurred in multiple networks these were the set of all interactions observed across networks and the mean of recorded trait values, respectively. Second, the phylogeny we used to model trait evolution was built from the time-calibrated, molecular phylogeny of hawkmoths from Kawahara and Barber (2015). Kawahara and Barber (2015) used five nuclear loci and a single mitochondrial locus that totaled 7449 bp to examine the evolutionary relationships between hawkmoth species. This phylogeny included all but one of the genera present in our dataset (a single species in our data) and thus provided, at a minimum, genus-level divergence dates for 74 out of 75 species. Below the genus level, we collapsed intra-genus dates into polytomies as not all congeners in our dataset were present in the phylogeny which prevented us from reliably inferring intra-genus relationships between taxa. The presence of unresolved polytomies in this phylogeny is a limitation of the data in this study. Descriptive analyses of the data. Given that our focus in this study was to examine how well the evolution of hawkmoth proboscis length may be explained by the relevant traits (i.e., corolla length) of the plants they pollinate, we also described each hawkmoth species in one additional way. We calculated the effective interaction trait (EIT) of the plants that each hawkmoth pollinates. The EIT of each hawkmoth is the weighted mean corolla length of the set of plant species, it was observed to pollinate where weights are given by the frequency with which the hawkmoth visited each plant. Ecologically, this quantity for each hawkmoth represents the corolla length that best describes their interaction patterns and interaction preferences. In the remainder of the text, we will focus on analyses and results for hawkmoth EIT so defined; however, our results do not differ qualitatively if we use an unweighted EIT based solely on partner identity (see Supplementary Material S1 available on Dryad at http://dx.doi.org/10.5061/dryad.3vb73). To get an initial indication of the relationship between hawkmoth traits and phylogenetic relatedness, we estimated the level of phylogenetic signal of both proboscis length and EIT using Blomberg’s K (Blomberg et al. 2003). In each case, we tested whether the observed phylogenetic signal was significant with a null model comprising 10,000 shuffled trait assignments with the phytools::phylosig function in R Core Team (2013) and Revell (2012). We also examined the relationship between hawkmoth proboscis length and EIT with a phylogenetic least-squares approach (Orme 2013) to explore the extent to which the raw data captures a relationship between interactions and traits. Finally, we tested whether any morphological matching between hawkmoth proboscis length and EIT—measured as the absolute difference between the two—can be explained by the potentially confounding factors of hawkmoth lineage age (i.e., the phylogenetic branch length pertaining to a distinct species) and hawkmoth generalism (i.e., the number of interaction partners a hawkmoth has) using the same phylogenetic least-squares approach as above. Interaction-Informed Evolutionary Models of Proboscis Evolution To explore the extent to which pollination interactions can be informative descriptors of hawkmoth trait evolution, we started with a baseline evolutionary model and then constructed five additional models that explicitly incorporate both ecological interactions and specific hypotheses regarding the biological scale at which those interactions are informative. In the baseline model, traits evolve according to Brownian motion; i.e., the evolution of a trait $$X$$ over time is modeled as $$dX(t) = \sigma d\beta(t)$$, where $$d\beta (t)$$ is a random deviate and $$\sigma$$ is the volatility of that deviation. Herein, a pure-drift process acts and all species evolve independently from each other. The resulting differences in traits between species are a product of their phylogenetic relatedness (Felsenstein 1985; Butler and King 2004). We fit this baseline model with the ouch::brown function in (R Core Team 2013) and (King and Butler, 2009). In contrast to this baseline model, all of our interaction-informed models are versions of an Ornstein–Uhlenbeck process where trait evolution proceeds in the presence of adaptive peaks and potential selection towards such peaks (Hansen 1997; Butler and King 2004). The most simple Ornstein–Uhlenbeck model of trait evolution takes the form $$dX(t) = \alpha[\theta - X(t)]dt + \sigma d \beta(t),$$ (1) where $$X$$ is the trait of interest, $$\theta$$ is the evolutionarily optimal trait value, $$\alpha$$ is the rate of attraction to the optimum, and $$d \beta(t)$$ and $$\sigma$$ are as before (i.e., the components of the Brownian model). The key difference between Equation 1 and our interaction-informed models relates to the optimal trait value $$\theta$$. To address the question of how a hawkmoth’s current interactions may reflect the evolution of its proboscis length, we wanted to test hypotheses about specific values of $$\theta$$ that can be calculated from the traits of the plants they pollinate. In contrast, previous approaches maintain $$\theta$$ as a free parameter that is also estimated during model fitting with a Bayesian or Maximum Likelihood approach (Butler and King 2004; Mahler et al. 2013; Uyeda and Harmon 2014; Khabbazian et al. 2016). Similar to these previous approaches, we also group species on the phylogeny into selective regimes (i.e., groups of species hypothesized to be evolving towards the same evolutionary optima). In current approaches, each selective regime is associated with an estimated optimum trait value (Butler and King 2004; Mahler et al. 2013; Uyeda and Harmon 2014). Rather than examining how latent optima can be best-fit given a phylogenetic tree, trait distribution, and hypothesized selective regimes, our approach pre-specifies $$\theta$$ values for each selective regime directly from the empirical data (in this case, data on species interactions). To implement this model, we used a variation on the classic Ornstein–Uhlenbeck model—known as the Hansen (1997) model —and drew on the implementation of Butler and King (2004) and King and Butler (2009). All R code and data used in these analyses are freely available on Data Dryad (doi:10.5061/dryad.3vb73). Building on this approach, we calculated each group $$g$$’s optimal value for proboscis length ($$\theta_{g}$$) from its constituent hawkmoths’ pollination interactions and the traits of their flowering-plant interaction partners (i.e., an EIT value for each group). As with the simplest Ornstein–Uhlenbeck models, our model also included an ancestrally optimal trait value ($$\theta_N$$; estimated during model fitting) that allowed for the fact that contemporary pollination interactions, and the optima based on them, may not have exerted persistent selection across the full depth of the phylogeny. Finally, we included a dummy variable $$\omega_{i}(t)$$ for each species $$i$$ in group $$g$$ that governed its evolution towards either $$\theta_{g}$$ or $$\theta_{N}$$ at any point in time. Species $$i$$’s trait evolution is attracted towards the null optimum ($$\theta_N$$) when $$\omega_{i}(t) = 0$$ and towards $$\theta_{g}$$ when $$\omega_{i}(t) = 1$$. Our implementation of the Hansen (1997) model thus takes the form: $$dX_i (t) = \alpha \left[ \left(1 - \omega_i(t)\right)\theta_N + \omega_i(t)\, \theta_{g} -X_i(t) \right] + \sigma d \beta(t).$$ (2) As mentioned, we used this model to test ecological hypotheses regarding the specific values of evolutionary optima for proboscis length based on hawkmoth interactions. While our overall objective was to assess how well evolutionary optima derived from contemporary interactions can describe past trait evolution, we also leveraged this approach to address a related question regarding the biological scale at which these interactions are most informative. To do this, we compared and contrasted several sets of hypothesized selective regimes and optima. These sets fall broadly into two categories—taxonomic and functional—both of which are important ecological groupings of pollinators (Olesen et al. 2007; Rezende et al. 2009; Gómez et al. 2010) and have allowed us to assess whether proboscis-length evolution is better described by hawkmoth interactions at one or both scales. For example, to test the hypothesis that the evolution of hawkmoth proboscis length is best explained by pollination interactions at the genus level, we took the following four steps. First, we estimated the interaction-based optimum for each genus as the EIT of all hawkmoth species in that genus (where weights are given by the set of hawkmoth–plant interactions of members of the genus). In the model, the EIT of a genus represents the evolutionarily optimal proboscis length for that genus (i.e., the trait value that members of said genus will evolve towards when $$\omega_{i}(t) = 1$$); the genus itself constitutes the selective regime—a group of taxa that evolve towards the same $$\theta_{g}$$ (Fig. 1a–d). Second, we fixed the point in time at which selection toward these genus-specific optima can begin; that is, the time $$t_i^*$$ at which $$\omega_{i}(t > t_i^*) = 1$$ (Fig. 1e). Note that a hawkmoth family composed of three genera would have three different evolutionary optima in this model (i.e., one specific to each genus). Therefore, a hawkmoth species in genus A could not experience selection toward the optimal trait value of its genus until the point in time at which it diverges from genera B and C. Until that point, all species evolve according to the ancestral optimum $$\theta_N$$ (i.e. $$\omega_{i}(t < t_i^*) = 0$$). We then fit the model—with its set of hypothesized optima $$\{\theta_{g}\}$$ and appropriate $$\omega_i(t)$$ and $$t_i^*$$ values—with a maximum likelihood approach to find the best-fit values of the three remaining parameters ($$\alpha$$, $$\sigma$$, $$\theta_N$$). We implemented this procedure with the sbplx function (based on Rowan 1990) in the nloptr package in R Core Team (2013) and Ypma (2014). Finally, we calculated the model-selection statistic corrected Akaike Information Criterion (AICc) to facilitate comparison across hypotheses (Akaike 1998). Note that the number of $$\theta_{g}$$ does not have an impact on AICc as these parameters are not estimated in model-fitting (Burnham and Anderson 2001). Figure 1. View largeDownload slide Conceptual depiction of the process of how we determine hypothetical selective regimes and their corresponding optima for our interaction-informed Ornstein–Uhlenbeck model. First, we combine (a) data on pollination interactions between hawkmoths and plants with (b) observed plant corolla lengths to generate (c) an effective interaction trait (EIT) for each hawkmoth. In (a), interactions are represented by lines joining a hawkmoth and plant where line thickness represents interaction intensity. In (c), we show the EIT for each hawkmoth and they are grouped into selective regimes based on a hypothesis of how species evolve relative to each other. Herein, hawkmoth species are assigned to three functional selective regimes (orange, purple, and green) based on their EIT values. d) For each selective regime, an optimum trait value is calculated as the mean EIT of the group. e) Finally, the selective regimes and their optima are painted onto the phylogeny and the model parameters can be fit to determine the degree to which those selective regimes and optima explain observed variation in hawkmoth proboscis length. Figure 1. View largeDownload slide Conceptual depiction of the process of how we determine hypothetical selective regimes and their corresponding optima for our interaction-informed Ornstein–Uhlenbeck model. First, we combine (a) data on pollination interactions between hawkmoths and plants with (b) observed plant corolla lengths to generate (c) an effective interaction trait (EIT) for each hawkmoth. In (a), interactions are represented by lines joining a hawkmoth and plant where line thickness represents interaction intensity. In (c), we show the EIT for each hawkmoth and they are grouped into selective regimes based on a hypothesis of how species evolve relative to each other. Herein, hawkmoth species are assigned to three functional selective regimes (orange, purple, and green) based on their EIT values. d) For each selective regime, an optimum trait value is calculated as the mean EIT of the group. e) Finally, the selective regimes and their optima are painted onto the phylogeny and the model parameters can be fit to determine the degree to which those selective regimes and optima explain observed variation in hawkmoth proboscis length. We followed the above procedure for five models of proboscis-length evolution: a single optimum for all species (i.e., all hawkmoths in the same selective regime; referred to as OU1), two models based on taxonomic groupings, and two models based on functional groupings. The taxonomic models assessed how well hawkmoth interactions at the genus- and species-specific scales (i.e., separate $$\theta_{g}$$ for each genus or species, respectively) explain observed hawkmoth proboscis lengths (OU2 and OU3). The functional models aimed to assess how hawkmoth interactions can explain observed traits when hawkmoths are grouped based on the similarity of their interactions at both a coarse scale and finer scale. In the former, hawkmoths fall into two groups based their EIT value ($$x < 38\,mm$$ and $$x > 38\,mm$$; OU4). In the latter, species were assigned to one of six groups ($$x < 15\,mm$$, $$15\,mm < x < 30\,mm$$, $$30\,mm < x < 38\,mm$$, $$38\,mm < x < 50\,mm$$, $$50\,mm < x < 70\,mm$$, and $$x > 70\,mm$$; OU5). In both cases, we chose trait bins based on natural groupings of EIT values. In the functional-group models, all values $$t_{i}^*$$ were given by the time that a species diverged from other taxa (i.e., the species’ age) since functional groupings were scattered across the phylogeny. For each model, we assessed how their hypothesized $$\theta_g$$ values related to the observed hawkmoth traits. To do so, we calculated the mean-squared error (MSE) for each model where observed hawkmoth proboscis lengths were the “observed” values and the optimum proboscis length associated to each species was the corresponding “predicted” value. The MSE for each model therefore represents how close the matching between trait and optimum value was across the tree. Closer matching between the two may explain differences in model-fit and provide insights to the selective regimes where contemporary interactions best describe traits. Comparison to Best-Fit Evolutionary Model The primary purpose of our hypothesis-based approach was to examine the degree to which incorporating specific ecological data (species interactions) into models of trait evolution can make their results more informative. Our hypothesis-based approach provides insights into the ability of contemporary interactions to capture past trait evolution and the biological scale at which that ability is maximized. However, the flip side of this approach is that the selective regimes of relatively few candidate models can realistically be biologically-justified and tested. Consequently, this makes it unlikely that any one of the few scenarios tested here exactly represents the most-likely model of proboscis-length evolution. Though this does not compromise our ability to test different hypotheses (Burnham and Anderson 2001), it does imply a lack of an upper bound to which the interaction-informed models can be compared. Identifying an upper bound can thus provide added insight into how good the fit of our best interaction-informed model is and how well its selective regimes capture patterns of trait evolution. Therefore, to further examine the results, we also aimed to identify the statistically-best model of proboscis-length evolution with the SURFACE algorithm (Ingram and Mahler 2013). In SURFACE, both species’ groups (selective regimes) and the optimal trait values that describe them are estimated by the model (Ingram and Mahler 2013). Using an iterative AICc procedure, the method assesses different sets of optima and shifts (analogous to the $$\theta_{g}$$, $$\omega_i$$, and $$t_i^*$$ in our model), calculates the AICc for the model with those parameter sets, and converges on the best-fit model—as defined by the smallest AICc value (Ingram and Mahler 2013). SURFACE has previously been used to successfully identify best-fit OU models and convergent evolution in other taxa (Mahler et al. 2013; Grundler and Rabosky 2014; Almécija et al. 2015). SURFACE is also just one of multiple, related methods that we could have used to make this comparison; for brevity, we focus on its comparison here in the main text but also compare the results of our approach to SLOUCH (Hansen et al. 2008) in the Supplementary Material available on Dryad. We also assessed whether or not a fingerprint of the interactions that were made explicit in our informed models could also be identified in the SURFACE results. Here, we used a linear model to explore the relationship between the best-fit optima of the SURFACE model and hawkmoth EIT. If hawkmoth interactions are reflected in the SURFACE optima, we would expect to see a positive relationship between the two. To assess the likelihood that the strength of the observed relationship between EIT and optima may be observed by chance, we also compared the observed relationship to a null distribution of 9999 random predictor test statistics. Each value in the null model was found by randomizing a hawkmoth’s pollination interactions, recalculating EIT, and then re-examining the relationship between EIT and optima with the linear model. To generate the randomized EIT values, we shuffled the hawkmoth-plant interaction network such that a hawkmoth’s interaction partners may change but overall interaction frequency and degree (i.e., number of interaction partners) are maintained (algorithm swsh_both_c in the vegan::nullmodel function; Oksanen et al. 2016). As our hypothesis suggests that the observed relationship between EIT and optima should be stronger than expected by chance, we performed a one-tailed test ($$\alpha = 0.05$$). Next, we examined the results of the SURFACE model with respect to the groupings of hawkmoths used in our interaction-informed models. The degree to which our hypothesized groupings are captured by the best-fit model can indicate how well a basic assumption about the evolution of a trait can approximate the best statistical groupings. We first quantified how much the selective regimes of the SURFACE model overlap with each of our interaction-informed models by measuring their mutual information (Danon et al. 2005). As before, we sought to compare the observed mutual information to what might be expected by chance, so we shuffled the SURFACE regime assignments 1000 times across the hawkmoth species and reassessed the mutual information. We again performed a one-tailed test and hence considered the observed overlap significantly different than expected by chance if the observed mutual information was greater than the null distribution ($$\alpha = 0.05$$). We also compared the SURFACE selective regimes to our six-functional-group model to see how these two scenarios captured biological groupings of hawkmoths. First, we compared the phylogenetic signal of each model’s regimes to assess the degree to which hawkmoth relatedness predicted selective regime co-membership. To do so, we converted group participation into a binary design matrix upon which we were able to estimate phylogenetic signal as the multivariate K statistic of (Adams, 2014) with the function geomorph::physignal and evaluate it based on the Brownian expectation $$K = 1$$ (Adams and Otárola-Castillo 2013). Second, we estimated the degree to which the regimes represent groupings of morphologically similar hawkmoths. To do so, we examined the variation of hawkmoth proboscis length within selective regimes. We used an ANOVA to compare hawkmoth proboscis length across selective regimes, where a larger F-statistic—a relatively larger amount of variation in proboscis length explained by selective regime—would suggest tighter functional groupings of hawkmoths. Results In this group of hawkmoths, we observed distinct and significant phylogenetic signal in their proboscis lengths ($$K = 1.169$$, $$P < 0.001$$, $$n = 10000$$), suggesting that more closely related hawkmoths tend to have more similar feeding traits (Fig. 2a). We saw a similar, but weaker, pattern when we expanded our focus to include the plants those hawkmoths visit. Phylogenetic signal of hawkmoth EIT was lower than what we observed for proboscis length but still significant ($$K = 0.412$$, $$P = 0.008$$, $$n = 10000$$; Fig. 2a). On the other hand, there was a significant positive relationship between hawkmoth proboscis length and EIT (PGLS: $$\beta = 0.501$$, $$t = 3.386$$, $$P = 0.001$$, $$R^2 = 0.136$$; Fig. 2b) suggesting, as expected, that hawkmoths with longer proboscides tend to visit longer flowers. Interestingly, the degree of morphological matching we observed between an hawkmoth’s proboscis length and EIT is not explained by either the age of the hawkmoth lineage on our phylogeny or its number of observed interaction partners (PGLS: $$\beta = -0.095$$, $$t = -0.187$$, $$P = 0.852$$, $$R^2 < 0.001$$ & $$\beta = 0.287$$, $$t = 1.479$$, $$P = 0.143$$, $$R^2 = 0.030$$, respectively). Figure 2. View largeDownload slide Relationship between hawkmoth phylogeny, proboscis length (PL; mm), and effective interaction trait (EIT; mm). a) There is distinct phylogenetic signal of proboscis length compared with EIT ($$K = 1.169$$, $$P < 0.001$$, $$n = 10000$$ and $$K = 0.412$$, $$P = 0.008$$, $$n = 10000$$, respectively; where $$K > 1$$ indicates phylogenetic clustering). b) There is a significant, positive relationship between proboscis length and EIT when accounting for phylogenetic relatedness ($$t = 3.386$$, $$P = 0.001$$, $$R^2 = 0.136$$). The dashed grey line in (b) represents the fitted PGLS relationship between proboscis length and EIT where $$\beta = 0.501$$. Figure 2. View largeDownload slide Relationship between hawkmoth phylogeny, proboscis length (PL; mm), and effective interaction trait (EIT; mm). a) There is distinct phylogenetic signal of proboscis length compared with EIT ($$K = 1.169$$, $$P < 0.001$$, $$n = 10000$$ and $$K = 0.412$$, $$P = 0.008$$, $$n = 10000$$, respectively; where $$K > 1$$ indicates phylogenetic clustering). b) There is a significant, positive relationship between proboscis length and EIT when accounting for phylogenetic relatedness ($$t = 3.386$$, $$P = 0.001$$, $$R^2 = 0.136$$). The dashed grey line in (b) represents the fitted PGLS relationship between proboscis length and EIT where $$\beta = 0.501$$. Given this baseline picture of the interplay between hawkmoth proboscis length, interactions, and evolutionary history, it was still unclear how informative ecological interactions might ultimately be for describing proboscis-length evolution or the biological scale at which this would be maximal. Our models, however, suggested that trait evolution was best captured by contemporary ecological interactions when hawkmoths are grouped at a fine functional scale or on their own as individual species. Indeed, only models with optima that were fitted at the species (OU3) and six-functional-group scale (OU5) performed better than Brownian motion at explaining proboscis-length evolution (Table 1). In contrast, the models with a single global optimum (the collective EIT of all hawkmoths in the phylogeny), with genus-specific optima, and with two functional groups all had higher AICc values—and were therefore poorer candidate models—than the Brownian-motion model (Table 1). In terms of AICc, the six-functional-group model and the species-specific model were indistinguishable (OU3 and OU5; Table 1). Table 1 Summary of model fit for each scenario for hawkmoth proboscis evolution Model Type $$\alpha$$ $$\sigma$$ $$\theta_{N}$$ nSR lnL $$n$$ AICc MSE BM — — 25.240 40.604 — $$-$$306.515 2 617.196 — OU1 Global 0.468 27.391 43.245$$^{\mathrm{\dagger}}$$ 1 $$-$$307.591 2 621.521 590.562 OU2 Tax. 0.030 24.770 43.329 25 $$-$$306.935 3 620.208 349.916 OU3 Tax. 0.025 23.367 28.410 75 $$-$$302.702 3 611.742 404.395 OU4 Func. 0.012 25.177 17.840 2 $$-$$308.856 3 624.050 477.653 OU5 Func. 0.026 23.427 28.993 6 $$-$$302.855 3 612.048 395.813 SF — 0.069 2.856 27.599* 7 $$-$$253.744 9 528.258 119.258 Model Type $$\alpha$$ $$\sigma$$ $$\theta_{N}$$ nSR lnL $$n$$ AICc MSE BM — — 25.240 40.604 — $$-$$306.515 2 617.196 — OU1 Global 0.468 27.391 43.245$$^{\mathrm{\dagger}}$$ 1 $$-$$307.591 2 621.521 590.562 OU2 Tax. 0.030 24.770 43.329 25 $$-$$306.935 3 620.208 349.916 OU3 Tax. 0.025 23.367 28.410 75 $$-$$302.702 3 611.742 404.395 OU4 Func. 0.012 25.177 17.840 2 $$-$$308.856 3 624.050 477.653 OU5 Func. 0.026 23.427 28.993 6 $$-$$302.855 3 612.048 395.813 SF — 0.069 2.856 27.599* 7 $$-$$253.744 9 528.258 119.258 Notes: $$\alpha$$ = strength of attraction; $$\sigma$$ = magnitude of variation; $$\theta_{N}$$ = intrinsic optimum; nSR = number of selective regimes in model (in interaction-informed models this is also the number of $$\theta_{g}$$ values fixed); lnL = log-likelihood; n = number of free parameters; AICc = small-sample-size-corrected Akaike Information Criterion; MSE = mean-squared error between $$\theta_{g}$$ values and observed traits; BM = Brownian motion; OU1 = global-optimum; OU2 = genus-specific optima; OU3 = species-specific optima; OU4 = two-functional-group-based optima; OU5 = six-functional-group-based optima; SF = SURFACE implementation. “Tax.” refers to taxonomic groupings of species, “Func.” refers to groupings based on similarity of interactions, “Global” refers to no groupings. A “—” indicates that the value is not applicable to the model, a “†” indicates that $$\theta_{N}$$ was equal to the fitted optima, and the “*” after the $$\theta_{N}$$ value for the SURFACE model refers to the estimated $$\theta$$ value at the root. Table 1 Summary of model fit for each scenario for hawkmoth proboscis evolution Model Type $$\alpha$$ $$\sigma$$ $$\theta_{N}$$ nSR lnL $$n$$ AICc MSE BM — — 25.240 40.604 — $$-$$306.515 2 617.196 — OU1 Global 0.468 27.391 43.245$$^{\mathrm{\dagger}}$$ 1 $$-$$307.591 2 621.521 590.562 OU2 Tax. 0.030 24.770 43.329 25 $$-$$306.935 3 620.208 349.916 OU3 Tax. 0.025 23.367 28.410 75 $$-$$302.702 3 611.742 404.395 OU4 Func. 0.012 25.177 17.840 2 $$-$$308.856 3 624.050 477.653 OU5 Func. 0.026 23.427 28.993 6 $$-$$302.855 3 612.048 395.813 SF — 0.069 2.856 27.599* 7 $$-$$253.744 9 528.258 119.258 Model Type $$\alpha$$ $$\sigma$$ $$\theta_{N}$$ nSR lnL $$n$$ AICc MSE BM — — 25.240 40.604 — $$-$$306.515 2 617.196 — OU1 Global 0.468 27.391 43.245$$^{\mathrm{\dagger}}$$ 1 $$-$$307.591 2 621.521 590.562 OU2 Tax. 0.030 24.770 43.329 25 $$-$$306.935 3 620.208 349.916 OU3 Tax. 0.025 23.367 28.410 75 $$-$$302.702 3 611.742 404.395 OU4 Func. 0.012 25.177 17.840 2 $$-$$308.856 3 624.050 477.653 OU5 Func. 0.026 23.427 28.993 6 $$-$$302.855 3 612.048 395.813 SF — 0.069 2.856 27.599* 7 $$-$$253.744 9 528.258 119.258 Notes: $$\alpha$$ = strength of attraction; $$\sigma$$ = magnitude of variation; $$\theta_{N}$$ = intrinsic optimum; nSR = number of selective regimes in model (in interaction-informed models this is also the number of $$\theta_{g}$$ values fixed); lnL = log-likelihood; n = number of free parameters; AICc = small-sample-size-corrected Akaike Information Criterion; MSE = mean-squared error between $$\theta_{g}$$ values and observed traits; BM = Brownian motion; OU1 = global-optimum; OU2 = genus-specific optima; OU3 = species-specific optima; OU4 = two-functional-group-based optima; OU5 = six-functional-group-based optima; SF = SURFACE implementation. “Tax.” refers to taxonomic groupings of species, “Func.” refers to groupings based on similarity of interactions, “Global” refers to no groupings. A “—” indicates that the value is not applicable to the model, a “†” indicates that $$\theta_{N}$$ was equal to the fitted optima, and the “*” after the $$\theta_{N}$$ value for the SURFACE model refers to the estimated $$\theta$$ value at the root. We saw that MSE of proboscis length and $$\theta_g$$ somewhat corresponded to the pattern of model-fit. In the species-specific (OU3) and six-functional group model (OU5), matching between the two was close while it was most discordant at the global scale (OU1; Table 1). Most interestingly, the model with genus-specific optima (OU2) showed better matching than expected from its model fit as its MSE was the lowest of our interaction-informed models despite it showing a fit worse than Brownian Motion based on both $$lnL$$ and AICc (OU2; Table 1). While some of our interaction-informed models provided an improved description of hawkmoth proboscis-length evolution relative to Brownian motion, the unconstrained model of proboscis-length evolution (identified by SURFACE) fit substantially better than our best interaction-informed model, as expected ($$lnL = -253.744$$, $$\alpha = 0.069$$, $$\sigma = 2.856$$, $$AICc = 528.258$$, $$\#\theta = 7$$; Table 1). The MSE of optima and traits for this model was also the lowest of all the models that we fit (Table 1). There was a clear positive relationship between hawkmoth EIT and optima for the SURFACE model ($$\beta = 0.436$$, $$t = 7.018$$, $$P < 0.001$$, $$R^2 = 0.403$$; Fig. 3a). The observed hawkmoth EIT values showed a stronger relationship with the SURFACE optima than any of the models of the null distribution (Fig. 3b). The positive relationship between EIT and optima suggested that the regimes and optima upon which SURFACE converges are at least correlated with hawkmoth interactions. Perhaps unsurprisingly, this relationship is weaker than the same for the six-functional-group model where the optima were drawn from the EIT values of each functional group. Figure 3. View largeDownload slide The relationship between a hawkmoth species’ optimal proboscis length (PL) and its observed effective interaction trait (EIT). a) There is a significant, positive correlation between a species’ EIT and its optimal proboscis length as determined by SURFACE. The solid line shows the slope from a linear regression of the values while the dashed line shows the same for our best-fitting interaction-informed model for comparative purposes. The observed slope between EIT and SURFACE optima is smaller ($$\beta= 0.44$$) than that from EIT and the optima of the interaction-informed model, where optima are based on EIT values ($$\beta = 1.01$$; dashed line). b) We show the degree to which the observed relationship between EIT and SURFACE optima is stronger than expected by chance. The vertical line indicates the observed t-value and the distribution represents 9999 null models in which hawkmoth EIT is shuffled between species and the relationship between SURFACE optima and EIT reassessed. Figure 3. View largeDownload slide The relationship between a hawkmoth species’ optimal proboscis length (PL) and its observed effective interaction trait (EIT). a) There is a significant, positive correlation between a species’ EIT and its optimal proboscis length as determined by SURFACE. The solid line shows the slope from a linear regression of the values while the dashed line shows the same for our best-fitting interaction-informed model for comparative purposes. The observed slope between EIT and SURFACE optima is smaller ($$\beta= 0.44$$) than that from EIT and the optima of the interaction-informed model, where optima are based on EIT values ($$\beta = 1.01$$; dashed line). b) We show the degree to which the observed relationship between EIT and SURFACE optima is stronger than expected by chance. The vertical line indicates the observed t-value and the distribution represents 9999 null models in which hawkmoth EIT is shuffled between species and the relationship between SURFACE optima and EIT reassessed. We quantified the similarity of selective regimes between SURFACE and each of our models as the mutual information shared by the two groupings of hawkmoths into selective regimes. The regimes of SURFACE and the six-functional-group model have mutual information $$MI = 0.232$$ (Monte Carlo test; $$P < 0.004$$), where $$1$$ indicates perfect overlap and $$0$$ is no overlap. Interestingly, when comparing the SURFACE regimes to our other candidate interaction-informed models, the genus-wise model showed a higher mutual information than our best-fitting models ($$MI = 0.666$$, $$P < 0.001$$) while the global ($$MI = 0.000$$, $$P = 1.00$$) and two-functional-group ($$MI = 0.144$$, $$P = 0.324$$) showed lower values. The species-specific model showed a higher value than the six-functional-group model ($$MI = 0.544$$, $$P < 0.001$$). Although both our species-specific and six-functional-group scenarios performed as well as each other, we focused on the SURFACE and six-functional-group models when comparing the similarities of selective regimes as comparisons to the species-specific scenario would be uninformative. We saw that the SURFACE selective regimes show significantly greater phylogenetic signal than expected under Brownian motion ($$K = 1.518$$, $$P < 0.001$$, $$n = 10000$$; Supplementary Material S2 available on Dryad) as opposed to the regimes of the six-functional-group model which showed significantly less ($$K = 0.289$$, $$P = 0.040$$, $$n = 10000$$; Fig. 4a). Moreover, the SURFACE regimes also represented tighter morphological groupings of hawkmoths (in terms of the degree to which selective regime explains proboscis-length variation) than the six-functional-group model ($$F = 91.76$$, $$P < 0.001$$ & $$F = 10.4$$, $$P < 0.001$$, respectively). These results suggest that the SURFACE regimes captured taxonomic and morphological similarity between hawkmoths which is not captured by the functional groupings that we define based on hawkmoth EIT values and might explain the discrepancy in fit between the two models (Table 1). Figure 4. View largeDownload slide Summary of the more informative of our best-fitting interaction-informed models: the six-functional-group Ornstein–Uhlenbeck model (OU5). a) The phylogeny of hawkmoths in our dataset. Black lines indicate the periods in which selection acted toward the ancestral optimum $$\theta_N$$ and the different colors represent distinct functional groups and the distinct selective regimes. b) The hawkmoth trait space encompassed by each selective regime. Each distribution is made up of the effective interaction trait (EIT) values of hawkmoths in each of the six selective regimes. c) The optimal hawkmoth proboscis length for each selective regime that we used in the model; these values are the mean EIT of each functional group. Figure 4. View largeDownload slide Summary of the more informative of our best-fitting interaction-informed models: the six-functional-group Ornstein–Uhlenbeck model (OU5). a) The phylogeny of hawkmoths in our dataset. Black lines indicate the periods in which selection acted toward the ancestral optimum $$\theta_N$$ and the different colors represent distinct functional groups and the distinct selective regimes. b) The hawkmoth trait space encompassed by each selective regime. Each distribution is made up of the effective interaction trait (EIT) values of hawkmoths in each of the six selective regimes. c) The optimal hawkmoth proboscis length for each selective regime that we used in the model; these values are the mean EIT of each functional group. Discussion Altogether, our results provide a strong indication that contemporary ecological interactions can improve our understanding of past evolution of hawkmoth proboscis length. In demonstrating this link, we contribute an alternative approach to modeling trait evolution where the consideration of ecological interactions, or other aspects of natural history, can be built directly into the model parameters (see also Hansen et al. 2008; Uyeda and Harmon 2014; Kostikova et al. 2016) and potentially provide a more complete story of trait evolution. Furthermore, the biological realism incorporated into our model—in the form of ecological interactions—suggests that the evolution of proboscis length in this group of hawkmoths is explained by the contemporary pollination interactions of hawkmoths that interact with similar plants just as well as it is when each species and their interactions are considered separately. Lastly, our finding that contemporary interactions can be informative facets of trait evolution is supported by the fingerprint of these interactions that can be identified in an uninformed, unconstrained, and statistically best model. That interactions between species are predicated on their possession of traits selected to facilitate or avoid that interaction is a basic tenet of evolutionary ecology. Testing the assumption—that contemporary traits arose from the feedback of interaction-based selection—in the absence of historical data is difficult to impossible (however, see Gervasi and Schiestl 2017). Modeling evolutionary hypotheses is one approach to test questions and assumptions, such as these, that cannot be manipulated empirically (Butler and King 2004; Ingram et al. 2012; Mahler et al. 2013; Uyeda and Harmon 2014; Moen et al. 2015; Ingram et al. 2016). Our results show that, for a group of New World hawkmoths, contemporary pollination interactions can provide an informative lens through which the evolution of relevant pollination traits can be examined. In the dataset that we examine here, we see that hawkmoth traits correlate with their interactions even if one-to-one matching between proboscis and EIT is not observed. Intriguingly, our results suggest that matching between interactions and proboscis length is poor on a species-by-species scale since the interactions of functionally similar groups provide just as close of a match between trait and interactions as the species themselves do. These discrepancies may arise from the right-skewed trait distributions for pollinators in this data (Sazatornil et al. 2016) and others (Stang et al. 2009), but also deserves further attention. The biological scale at which traits evolve concordantly (i.e., converge) is an interesting question in its own right. For the hawkmoths studied here, we see that it is just as good to assume that interactions explain the evolution of functionally similar hawkmoth species as it is to assume that each hawkmoth is unique. Additionally, when these functional groups are too coarse-grained they can prove no more informative than a Brownian motion model. A result such as this might seem obvious based on the idea that species with similar interaction patterns should possess similar traits (Jordano 1995; Fenster et al. 2004; Olesen et al. 2007). However, phylogenetic conservation of interactions is a widely investigated phenomenon in evolutionary ecology (Rezende et al. 2007, 2009; Gómez et al. 2010) which would indicate that ecological interactions may well have captured evolution at a taxonomic scale such as the genus-scale that we examine. An hypothesis of phylogenetic conservation would also be supported by the phylogenetic signal of hawkmoths proboscis length and EIT that we see in this study. Our results show identical support between a finer-scale functional approach and species-specific approach which may question the role of taxonomic patterns in a wider evolutionary context and beggar further investigation. On the other hand, the best model from SURFACE, which fit better than any of our models, seemed to capture both taxonomic and morphological similarity in its selective regimes. Altogether, this suggests that future hypotheses of how species evolve in relation to each other require more nuance to disentangle the interplay between morphology and ancestry when it comes to the role of species interactions in directing adaptation. Models of trait evolution are rapidly moving on from Brownian motion and the basic Ornstein-Uhlenbeck approaches (Butler and King 2004; Hansen et al. 2008; Ingram and Mahler 2013; Uyeda and Harmon 2014; Moen et al. 2015; Khabbazian et al. 2016). However, as methods improve and computational power increases, it is important to consider the biological insight that can be drawn out of optimization-heavy approaches for modeling adaptation by natural selection. While these approaches can undeniably identify the statistically-best result—as demonstrated here—the biological implications of such results can be obscured by the variability of latent optima (Ho and Ané 2014). Spurious parameter estimates can arise due to flat likelihood surfaces that are common in evolutionary models (Ho and Ané 2014). At the same time, there is a movement to do a better job of incorporating species interactions into models of trait evolution (Nuismer and Harmon 2015; Drury et al. 2016; Manceau et al. 2016; Bartoszek et al. 2017). Our results use empirical data to demonstrate that even small alterations to established models can allow the incorporation of ecological interactions, change the way trait evolution modeling is approached, and provide insights into hypothetical parameter values that are rooted in biological reality. The obvious comparison within the scope of this study is the difference between the two modeling approaches. In the approach that we introduce—where empirical data is used to inform model parameters in a rigid manner—the potential biological relevance of the model is explicit. Hence this is an alternative to the recently proposed Bayesian approaches which can also consider additional data in the form of informative priors (Uyeda and Harmon 2014; Kostikova et al. 2016). On the other hand, the link between trait evolution and hawkmoth pollination interactions in the SURFACE results, despite a substantially-better fit, is only detectable a posteriori. SURFACE, the best-fitting model, and our approach, the biologically-explicit model, may well represent two ends of a spectrum, somewhere along which the best approach lies. In developing the approach that we present here, we add some ecological reality to a model of trait evolution. However, there are several areas in which this approach could be expanded beyond its current limitations. The major limitation of our approach currently is that it is rather static. First, the optima that we employ are fixed based on contemporary species interactions. We are therefore limited to the strong hypothesis that contemporary interactions are representative of—at least recent—evolutionary history even though species interactions can change on the order of weeks and months (Olesen et al. 2008; CaraDonna et al. 2017; Ponisio et al. 2017) and traits in long-lived species can evolve at quickest on the order of tens to hundreds of years (Galetti et al. 2013). Second, we only focus on one side of an interaction type that is regarded up as a textbook example of coevolution (Darwin 1862; Herre 1989; Anderson and Johnson 2008; Johnson and Anderson 2010). Presently, our approach does not consider how reciprocal selection has altered proboscis length, corolla length, and the interactions themselves through evolutionary time. To bring additional realism into an approach such as ours, a future adaptation could simultaneously model the evolution of traits on both sides of an interaction where optima on each side are updated through time based on how the traits of their interaction partners evolve. This would also require a realistic model of how interactions rewire through time but would provide a new approach to studying coevolution. Undeniably, species interactions are a foundational feature of ecological communities and are a source of selection pressure on trait evolution (Paine 1966; Thompson 2005). The approach we take here is amenable to the incorporation of data from a wide range of these interactions not just highly-specific mutualisms. Predator-prey relationships are ubiquitous and there are several examples of arms-races where trait change is thought to be driven by these interactions (Koskela and Ylönen 1995; Geffeney et al. 2002; Bro-Jørgensen 2013). For example, our method could be easily used to assess how hunting traits in coursing predators, such as top speed and endurance, track those same traits in their prey species. In this case, predators and their traits would be the focal group while predator EIT would be comprised of interaction data and the fleeing traits of their prey. Aside from species interactions, our approach could also be used to assess traits relating to the abiotic environment. For instance, it is common to assess optimal values or tolerance levels for temperature and nutrients in sessile organisms (Coles et al. 1976; Coles and Jokiel 1978; Ingestad 1979; Chapin III et al. 1983; Quesnel et al. 2006). These species are thought to be at risk from rapidly changing environmental conditions and their anthropogenic sources (Davis and Shaw 2001; McKenney et al. 2007; Ainsworth et al. 2016). Our method could be used to assess how groups such as corals, micro-algae, and trees may respond to these pressures. For example, species’ preferred nutrient levels or temperatures may used as the focal trait while current or previous values of the abiotic factor under consideration may be used for calculating model parameters. In all, the possibilities inherent in an interaction-informed modeling approach—where the question is explicit in the model—may represent a way in which evolutionary models can be used to test hypotheses about the links between traits and various aspects of a species’ ecology from interactions to abiotic conditions. A complete understanding of trait evolution in a group of species can require a wealth of information. Despite adding an additional layer of required data to comparative studies, we have demonstrated that knowledge of a group’s contemporary interaction partners can be useful for examining specific hypotheses about their trait evolution. Our results contribute to a growing emphasis on the need to take ecological interactions into account when examining species’ traits and their evolution (Nuismer et al. 2013; Nuismer and Harmon 2015; Manceau et al. 2016; Bartoszek et al. 2017). In particular for studies of systems that purportedly show interaction-selected morphology, our implementation of the Hansen model (Hansen 1997) may provide a new angle from which to explore the evolutionary pressure exerted on populations through their interactions with other species. Within the current biodiversity crisis (Barnosky et al. 2011; Ceballos et al. 2015), the loss of species interactions—or the appearance of novel ones—is guaranteed. Indeed, dramatic and rapid trait evolution has already been documented when interaction partners are lost (Galetti et al. 2013). If we are to grasp the consequences and understand the outcomes of the imminent—if not current—changes that ecological systems face, developing a better conceptual understanding of the interplay between ecological interaction and evolutionary adaptation will be key. Supplementary Material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.3vb73. Funding This work was supported by the University of Canterbury Summer Scholarship program and a Princeton University First Year Natural Sciences and Engineering Fellowship to M.C.H; São Paulo Research Foundation [FAPESP; 2013/13319-5 and 2014/20572-1 to M.P.G]; and Rutherford Discovery Fellowship, administered by the Royal Society of New Zealand to D.B.S. Acknowledgments The authors acknowledge the Associate Editor and three anonymous Reviewers for their valuable comments on this manuscript. Additionally, the authors acknowledge both sets of scientists whose work built the hawkmoth datasets that we incorporate as the example in this study. References Adams D.C. 2014 . A generalized K statistic for estimating phylogenetic signal from shape and other high-dimensional multivariate data. Syst. Biol. 63 : 685 – 697 . Google Scholar CrossRef Search ADS PubMed Adams D.C., Otárola-Castillo E. 2013 . geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods Ecol. Evol. 4 : 393 – 399 . Google Scholar CrossRef Search ADS Ainsworth, T.D., Heron S.F., Ortiz J.C., Mumby P.J., Grech A., Ogawa D., Eakin C.M., Leggat W. 2016 . Climate change disables coral bleaching protection on the Great Barrier Reef. Science 352 : 338 – 342 . Google Scholar CrossRef Search ADS PubMed Akaike H. 1998 . Information theory and an extension of the maximum likelihood principle. In: Selected Papers of Hirotugu Akaike. New York : Springer . p. 199 – 213 . Almécija, S., Smaers J.B., Jungers W.L. 2015 . The evolution of human and ape hand proportions. Nat. Commun. 6 : 7717 . Google Scholar CrossRef Search ADS PubMed Anderson, B., Johnson S.D. 2008 . The geographical mosaic of coevolution in a plant–pollinator mutualism. Evolution 62 : 220 – 225 . Google Scholar CrossRef Search ADS PubMed Barnosky, A.D., Matzke N., Tomiya S., Wogan G.O., Swartz B., Quental T.B., Marshall C., McGuire J.L., Lindsey E.L., Maguire K.C., Mersey B. 2011 . Has the earth’s sixth mass extinction already arrived? Nature 471 : 51 – 57 . Google Scholar CrossRef Search ADS PubMed Bartoszek, K., Glémin S., Kaj I., Lascoux M. 2017 . Using the Ornstein–Uhlenbeck process to model the evolution of interacting populations. J. Theor. Biol. 429 : 35 – 45 . Google Scholar CrossRef Search ADS PubMed Beaulieu, J.M., Jhwueng D.-C., Boettiger C., O’Meara B.C. 2012 . Modeling stabilizing selection: expanding the Ornstein–Uhlenbeck model of adaptive evolution. Evolution 66 : 2369 – 2383 . Google Scholar CrossRef Search ADS PubMed Blomberg, S.P., Garland T., Ives A.R. 2003 . Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution 57 : 717 – 745 . Google Scholar CrossRef Search ADS PubMed Bro-Jørgensen, J. 2013 . Evolution of sprint speed in african savannah herbivores in relation to predation. Evolution 67 : 3371 – 3376 . Google Scholar CrossRef Search ADS PubMed Brouat, C., Garcia N., Andary C., McKey D. 2001 . Plant lock and ant key: pairwise coevolution of an exclusion filter in an ant–plant mutualism. Proc. R. Soc. Lond. B Biol. Sci. 268 : 2131 – 2141 . Google Scholar CrossRef Search ADS Burnham, K.P., Anderson D.R. 2001 . Kullback-Leibler information as a basis for strong inference in ecological studies. Wildlife Res. 28 : 111 – 119 . Google Scholar CrossRef Search ADS Butler, M.A., King A.A. 2004 . Phylogenetic comparative analysis: a modeling approach for adaptive evolution. Am. Nat. 164 : 683 – 695 . Google Scholar CrossRef Search ADS PubMed CaraDonna, P.J., Petry W.K., Brennan R.M., Cunningham J.L., Bronstein J.L., Waser N.M., Sanders N.J. 2017 . Interaction rewiring and the rapid turnover of plant–pollinator networks. Ecol. Lett. 20 : 385 – 394 . Google Scholar CrossRef Search ADS PubMed Ceballos, G., Ehrlich P.R., Barnosky A.D., García A., Pringle R.M., Palmer T.M. 2015 . Accelerated modern human–induced species losses: Entering the sixth mass extinction. Sci. Adv. 1 : e1400253 . Google Scholar CrossRef Search ADS PubMed Chapin F.S. Tryon III, P.R., Cleve K.V. 1983 . Influence of phosphorus on growth and biomass distribution of Alaskan taiga tree seedlings. Can. J. For. Res. 13 : 1092 – 1098 . Google Scholar CrossRef Search ADS Coles, S., Jokiel P.L. 1978 . Synergistic effects of temperature, salinity and light on the hermatypic coral Montipora verrucosa. Mar. Biol. 49 : 187 – 195 . Google Scholar CrossRef Search ADS Coles, S.L., Jokiel P.L., Lewis C.R. 1976 . Thermal tolerance in tropical versus subtropical Pacific reef corals 30 : 159 – 166 . Danon, L., Diaz-Guilera A., Duch J., Arenas A. 2005 . Comparing community structure identification. J. Stat. Mech. Theory Exp. 2005 : P09008 . Darwin C. 1859 . On the Origin of the Species. London : Murray . Darwin C. 1862 . On the various contrivances by which British and foreign orchids are fertilised by insect, and the good effects of intercrossing. London, UK : John Murray . Davis, M.B., Shaw R.G. 2001 . Range shifts and adaptive responses to Quaternary climate change. Science 292 : 673 – 679 . Google Scholar CrossRef Search ADS PubMed Dehling, D.M., Töpfer T., Schaefer H.M., Jordano P., Böhning-Gaese K., Schleuning M. 2014 . Functional relationships beyond species richness patterns: trait matching in plant–bird mutualisms across scales. Glob. Ecol. Biogeogr. 23 : 1085 – 1093 . Google Scholar CrossRef Search ADS Drury, J., Clavel J., Manceau M., Morlon H. 2016 . Estimating the effect of competition on trait evolution using maximum likelihood inference. Syst. Biol. 65 : 700 – 710 . Google Scholar CrossRef Search ADS PubMed Feinsinger, P., Murray K.G., Kinsman S., Busby W.H. 1986 . Floral neighborhood and pollination success in four hummingbird-pollinated cloud forest plant species. Ecology 67 : 449 – 464 . Google Scholar CrossRef Search ADS Felsenstein J. 1985 . Phylogenies and the comparative method. Am. Nat. 125 : 1 – 15 . Google Scholar CrossRef Search ADS Fenster, C.B., Armbruster W.S., Wilson P., Dudash M.R., Thomson J.D. 2004 . Pollination syndromes and floral specialization. Annu. Rev. Ecol. Evol. Syst. 35 : 375 – 403 . Google Scholar CrossRef Search ADS Galetti, M., Guevara R., Côrtes M.C., Fadini R., Von Matter S., Leite A.B., Labecca F., Ribeiro T., Carvalho C.S., Collevatti R.G., Pires M.M., Guimarães P.R., Jr Brancalion P.H., Ribeiro M.C., Jordano P. 2013 . Functional extinction of birds drives rapid evolutionary changes in seed size. Science 340 : 1086 – 1090 . Google Scholar CrossRef Search ADS PubMed Geffeney, S., Brodie E.D., Ruben P.C. 2002 . Mechanisms of adaptation in a predator-prey arms race: Ttx-resistant sodium channels. Science 297 : 1336 – 1339 . Google Scholar CrossRef Search ADS PubMed Gervasi, D.D., Schiestl F.P. 2017 . Real-time divergent evolution in plants driven by pollinators. Nat. Commun. 8 : 14691 . Google Scholar CrossRef Search ADS PubMed Gómez, J.M., Verdú M., Perfectti F. 2010 . Ecological interactions are evolutionarily conserved across the entire tree of life. Nature 465 : 918 – 921 . Google Scholar CrossRef Search ADS PubMed Grundler, M.C., Rabosky D.L. 2014 . Trophic divergence despite morphological convergence in a continental radiation of snakes. Proc. R. Soc. Lond. B Biol. Sci. 281 : 20140413 . Google Scholar CrossRef Search ADS Hansen, T.F. 1997 . Stabilizing selection and the comparative analysis of adaptation. Evolution 51 : 1341 – 1351 . Google Scholar CrossRef Search ADS PubMed Hansen, T.F., Pienaar J., Orzack S.H. 2008 . A comparative method for studying adaptation to a randomly evolving environment. Evolution 62 : 1965 – 1977 . Google Scholar PubMed Haverkamp, A., Bing J., Badeke E., Hansson B.S., Knaden M. 2016 . Innate olfactory preferences for flowers matching proboscis length ensure optimal energy gain in a hawkmoth. Nat. Commun. 7 : 11644 . Google Scholar CrossRef Search ADS PubMed Herre, E.A. 1989 . Coevolution of reproductive characteristics in 12 species of new world figs and their pollinator wasps. Cell. Mol. Life Sci. 45 : 637 – 647 . Google Scholar CrossRef Search ADS Ho, L.S.T., Ané C. 2014 . Intrinsic inference difficulties for trait evolution with Ornstein-Uhlenbeck models. Methods Ecol. Evol. 5 : 1133 – 1146 . Google Scholar CrossRef Search ADS Ingestad, T. 1979 . Nitrogen stress in birch seedlings. Physiol. Plant. 45 : 149 – 157 . Google Scholar CrossRef Search ADS Ingram, T., Harmon L., Shurin J. 2012 . When should we expect early bursts of trait evolution in comparative data? Predictions from an evolutionary food web model. J. Evol. Biol. 25 : 1902 – 1910 . Google Scholar CrossRef Search ADS PubMed Ingram, T., Harrison A., Mahler D.L., del Rosario Castañeda, M. Glor R.E., Herrel A., Stuart Y.E., Losos J.B. 2016 . Comparative tests of the role of dewlap size in Anolis lizard speciation. Proc. Biol. Sci. 283 . pii: 20162199 . Google Scholar CrossRef Search ADS PubMed Ingram, T., Mahler D.L. 2013 . SURFACE: detecting convergent evolution from comparative data by fitting Ornstein-Uhlenbeck models with stepwise Akaike Information Criterion. Methods Ecol. Evol. 4 : 416 – 425 . Google Scholar CrossRef Search ADS Janzen, D.H., Martin P.S. 1982 . Neotropical anachronisms: the fruits the gomphotheres ate. Science 215 : 19 – 27 . Google Scholar CrossRef Search ADS PubMed Johnson, S.D., Anderson B. 2010 . Coevolution between food-rewarding flowers and their pollinators. Evol. Educ. Outreach 3 : 32 – 39 . Google Scholar CrossRef Search ADS Jordano, P. 1995 . Angiosperm fleshy fruits and seed dispersers: a comparative analysis of adaptation and constraints in plant-animal interactions. Am. Nat. 145 : 163 – 191 . Google Scholar CrossRef Search ADS Kawahara, A.Y., Barber J.R. 2015 . Tempo and mode of antibat ultrasound production and sonar jamming in the diverse hawkmoth radiation. Proc. Natl. Acad. Sci. USA 112 : 6407 – 6412 . Google Scholar CrossRef Search ADS Kessler, A., Halitschke R. 2009 . Testing the potential for conflicting selection on floral chemical traits by pollinators and herbivores: predictions and case study. Funct. Ecol. 23 : 901 – 912 . Google Scholar CrossRef Search ADS Khabbazian, M., Kriebel R., Rohe K., Ané C. 2016 . Fast and accurate detection of evolutionary shifts in Ornstein–Uhlenbeck models. Methods Ecol. Evol. 7 : 811 – 824 . Google Scholar CrossRef Search ADS King, A., Butler M. 2009 . ouch: Ornstein-Uhlenbeck models for phylogenetic comparative hypotheses. R Package Version 2. Available from: https://cran.r-project.org/web/packages/ouch/index.html. Koskela, E., Ylönen H. 1995 . Suppressed breeding in the field vole (Microtus agrestis): an adaptation to cyclically fluctuating predation risk. Behav. Ecol. 6 : 311 – 315 . Google Scholar CrossRef Search ADS Kostikova, A., Silvestro D., Pearman P.B., Salamin N. 2016 . Bridging inter-and intraspecific trait evolution with a hierarchical Bayesian approach. Syst. Biol. 65 : 417 – 431 . Google Scholar CrossRef Search ADS PubMed Labra, A., Pienaar J., Hansen T.F. 2009 . Evolution of thermal physiology in liolaemus lizards: adaptation, phylogenetic inertia, and niche tracking. Am. Nat. 174 : 204 – 220 . Google Scholar CrossRef Search ADS PubMed Machado, C.A., Robbins N., Gilbert M.T.P., Herre E.A. 2005 . Critical review of host specificity and its coevolutionary implications in the fig/fig-wasp mutualism. Proc. Natl. Acad. Sci. USA 102 : 6558 – 6565 . Google Scholar CrossRef Search ADS Mahler, D.L., Ingram T., Revell L.J., Losos J.B. 2013 . Exceptional convergence on the macroevolutionary landscape in island lizard radiations. Science 341 : 292 – 295 . Google Scholar CrossRef Search ADS PubMed Manceau, M., Lambert A., Morlon H. 2016 . A unifying comparative phylogenetic framework including traits coevolving across interacting lineages. Syst. Biol. 66 : 551 – 568 . McKenney, D.W., Pedlar J.H., Lawrence K., Campbell K., Hutchinson M.F. 2007 . Potential impacts of climate change on the distribution of North American trees. AIBS Bull. 57 : 939 – 948 . Miklashevichs, E., Röhrig H., Schell J., Schmidt J. 2001 . Perception and signal transduction of rhizobial NOD factors. Crit. Rev. Plant Sci. 20 : 373 – 394 . Google Scholar CrossRef Search ADS Moen, D.S., Morlon H., Wiens J.J. 2015 . Testing convergence versus history: convergence dominates phenotypic evolution for over 150 million years in frogs. Syst. Biol. 65 : 146 – 160 . Google Scholar CrossRef Search ADS PubMed Monteiro, L.R., Nogueira M.R. 2011 . Evolutionary patterns and processes in the radiation of phyllostomid bats. BMC Evol. Biol. 11 : 137 . Google Scholar CrossRef Search ADS PubMed Muschick, M., Indermaur A., Salzburger W. 2012 . Convergent evolution within an adaptive radiation of cichlid fishes. Curr. Biol. 22 : 2362 – 2368 . Google Scholar CrossRef Search ADS PubMed Nilsson L.A. 1988 . The evolution of flowers with deep corolla tubes. Nature 334 : 147 – 149 . Google Scholar CrossRef Search ADS Niven, J.E., Laughlin S.B. 2008 . Energy limitation as a selective pressure on the evolution of sensory systems. J. Exp. Biol. 211 : 1792 – 1804 . Google Scholar CrossRef Search ADS PubMed Nuismer, S.L., Harmon L.J. 2015 . Predicting rates of interspecific interaction from phylogenetic trees. Ecol. Lett. 18 : 17 – 27 . Google Scholar CrossRef Search ADS PubMed Nuismer, S.L., Jordano P., Bascompte J. 2013 . Coevolution and the architecture of mutualistic networks. Evolution 67 : 338 – 354 . Google Scholar CrossRef Search ADS PubMed Oksanen, J., Blanchet F.G., Kindt R., Legendre P., Minchin P.R., O’Hara R.B., Simpson G.L., Solymos P., Stevens M.H.H., Wagner H. 2016 . vegan: Community Ecology Package. R Package Version 2.3-3. Available from: https://cran.r-project.org/web/packages/vegan/index.html. Olesen, J.M., Bascompte J., Dupont Y.L., Jordano P. 2007 . The modularity of pollination networks. Proc. Natl. Acad. Sci. USA 104 : 19891 – 19896 . Google Scholar CrossRef Search ADS Olesen, J.M., Bascompte J., Elberling H., Jordano P. 2008 . Temporal dynamics in a pollination network. Ecology 89 : 1573 – 1582 . Google Scholar CrossRef Search ADS PubMed Orme, D. 2013 . The caper package: comparative analysis of phylogenetics and evolution in R. R package version 5. Available from: https://cran.r-project.org/web/packages/caper/index.html. Paine, R.T. 1966 . Food web complexity and species diversity. Am. Nat. 100 : 65 – 75 . Google Scholar CrossRef Search ADS Pellmyr, O. 2003 . Yuccas, yucca moths, and coevolution: a review. Ann. Mo. Bot. Gard. 90 : 35 – 55 . Google Scholar CrossRef Search ADS Ponisio, L.C., Gaiarsa M.P., Kremen C. 2017 . Opportunistic attachment assembles plant–pollinator networks. Ecol. Lett. 20 : 1261 – 1272 . Google Scholar CrossRef Search ADS PubMed Quesnel, P.-O., Côté B., Fyles J.W., Munson A.D. 2006 . Optimum nutrient concentrations and CND scores of mature white spruce determined using a boundary-line approach and spatial variation of tree growth and nutrition. J. Plant Nutr. 29 : 1999 – 2018 . Google Scholar CrossRef Search ADS R Core Team. 2013 . R: A language and environment for statistical computing. R Core Team . Available from: http://www.R-project.org/. Revell, L.J. 2012 . phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3 : 217 – 223 . Google Scholar CrossRef Search ADS Rezende, E.L., Albert E.M., Fortuna M.A., Bascompte J. 2009 . Compartments in a marine food web associated with phylogeny, body mass, and habitat structure. Ecol. Lett. 12 : 779 – 788 . Google Scholar CrossRef Search ADS PubMed Rezende, E.L., Lavabre J.E., Guimaraes P.R., Jordano P., Jr., Bascompte J. 2007 . Non-random coextinctions in phylogenetically structured mutualistic networks. Nature 448 : 925 – U6 . Google Scholar CrossRef Search ADS PubMed Rowan, T.H. 1990 . Functional stability analysis of numerical algorithms. Ph.D. thesis. University of Texas (Austin) . Sazatornil, F.D., Moré M., Benitez-Vieyra S., Cocucci A.A., Kitching I.J., Schlumpberger B.O., Oliveira P.E., Sazima M., Amorim F.W. 2016 . Beyond neutral and forbidden links: morphological matches and the assembly of mutualistic hawkmoth-plant networks. J. Anim. Ecol. 85 : 1586 – 1594 . Google Scholar CrossRef Search ADS PubMed Scales, J.A., King A.A., Butler M.A. 2009 . Running for your life or running for your dinner: what drives fiber-type evolution in lizard locomotor muscles? Am. Nat. 173 : 543 – 553 . Google Scholar CrossRef Search ADS PubMed Schluter, D., Price T.D., Rowe L. 1991 . Conflicting selection pressures and life history trade-offs. Proc. R. Soc. Lond. B Biol. Sci. 246 : 11 – 17 . Google Scholar CrossRef Search ADS Slabbekoorn, H., Smith T.B. 2002 . Habitat-dependent song divergence in the Little Greenbul: An analysis of environmental selection pressures on acoustic signals. Evolution 56 : 1849 – 1858 . Google Scholar CrossRef Search ADS PubMed Stang, M., Klinkhamer P.G., Waser N.M., Stang I., van der Meijden. E. 2009 . Size-specific interaction patterns and size matching in a plant–pollinator interaction web. Ann. Bot. 103 : 1459 – 1469 . Google Scholar CrossRef Search ADS Stiles, F.G. 1975 . Ecology, flowering phenology, and hummingbird pollination of some Costa Rican Heliconia species. Ecology 56 : 285 – 301 . Google Scholar CrossRef Search ADS Thompson, J.N. 1994 . The coevolutionary process. Chicago, IL : University of Chicago Press . Thompson, J.N. 2005 . The geographic mosaic of coevolution. Chicago, IL : University of Chicago Press . Uyeda, J.C., Harmon L.J. 2014 . A novel Bayesian method for inferring and interpreting the dynamics of adaptive landscapes from phylogenetic comparative data. Syst. Biol. 63 : 902 – 918 . Google Scholar CrossRef Search ADS PubMed Ypma, J. 2014 . nloptr: R interface to NLopt. R package version 1 . Available from: https://cran.r-project.org/web/packages/index.html. © The Author(s) 2018. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. All rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Systematic Biology Oxford University Press

# Contemporary Ecological Interactions Improve Models of Past Trait Evolution

, Volume 67 (5) – Sep 1, 2018
12 pages

/lp/ou_press/contemporary-ecological-interactions-improve-models-of-past-trait-21c3d5Dvw5
Publisher
Oxford University Press
ISSN
1063-5157
eISSN
1076-836X
D.O.I.
10.1093/sysbio/syy012
Publisher site
See Article on Publisher Site

### Abstract

Abstract Despite the fact that natural selection underlies both traits and interactions, evolutionary models often neglect that ecological interactions may, and in many cases do, influence the evolution of traits. Herein, we explore the interdependence of ecological interactions and functional traits in the pollination associations of hawkmoths and flowering plants. Specifically, we develop an adaptation of the Ornstein–Uhlenbeck model of trait evolution that allows us to study the influence of plant corolla depth and observed hawkmoth–plant interactions on the evolution of hawkmoth proboscis length. Across diverse modelling scenarios, we find that the inclusion of contemporary interactions can provide a better description of trait evolution than the null expectation. Moreover, we show that the pollination interactions provide more-likely models of hawkmoth trait evolution when interactions are considered at increasingly fine-scale groups of hawkmoths. Finally, we demonstrate how the results of best-fit modeling approaches can implicitly support the association between interactions and trait evolution that our method explicitly examines. In showing that contemporary interactions can provide insight into the historical evolution of hawkmoth proboscis length, we demonstrate the clear utility of incorporating additional ecological information to models designed to study past trait evolution. Macroevolution, mutualism, Ornstein-Uhlenbeck, pollination, Sphingidae For interactions to occur between species, both of the species involved must possess traits conducive to their occurrence. In the case of some ecological interactions, such as reproductive mutualisms, the interdependence of species across the interaction is long thought to lead to adaptation “in the most perfect manner” (Darwin 1859) such that the morphological traits of mutualists become tightly coupled (Thompson 1994). Indeed, morphological trait matching between mutualists is realized across a variety of qualitatively different ecological systems and taxa: from hummingbird-Heliconia pollination and bird-palm seed dispersal to nitrogen fixation mutualisms (Stiles 1975; Brouat et al. 2001; Miklashevichs et al. 2001; Stang et al. 2009; Galetti et al. 2013). Widespread morphological matching implies that in the most closely-coupled associations—whether at the scale of entire guilds (à la interaction syndromes; Janzen and Martin 1982; Fenster et al. 2004) or that of specific species pairs (Darwin 1862)—the morphology of one species may in fact be informative enough to infer the relevant traits and identities of their interaction partners. The evolutionary trajectories of traits on either side of an interaction are in fact influenced by multiple selection pressures, even in purported tightly coupled interactions (Schluter et al. 1991; Kessler and Halitschke 2009). In particular, this may come about since one-to-one interaction specialization is exceedingly rare; species long thought to be hyper-specialized often have more than one interaction partner in reality (Feinsinger et al. 1986; Pellmyr 2003; Machado et al. 2005), and therefore, experience more than one biotic selection pressure on the trait, or traits, relevant to that interaction. Additionally, trait evolution can be influenced by environmental factors, energy budget, and ancestry, amongst others (Hansen 1997; Slabbekoorn and Smith 2002; Niven and Laughlin 2008). Putting this all together, it remains open to debate whether and when an extant species’ interactions and interaction partners’ traits are informative with regard to past trait evolution. The pollination associations between hawkmoths and flowering plants are an ideal system in which to examine the utility of ecological interactions for understanding trait evolution given that close coupling of hawkmoth proboscis length and flower corolla depth is often observed in these interactions (Nilsson 1988; Haverkamp et al. 2016; Sazatornil et al. 2016). That said, morphological matching across interactions in this system is most typically examined from an ecological perspective—and with function as the primary focus—whereas the role of evolutionary processes in building congruence is potentially overlooked (Stang et al. 2009; Dehling et al. 2014; Sazatornil et al. 2016). Consequently, one approach to better explore this question would be to investigate whether a hawkmoth’s interaction partners and those partners’ corolla depths prove informative when modeling the evolution of hawkmoth proboscis length. Contemporary models of trait evolution tend to be adaptations of the random-walk Brownian motion and adaptive Ornstein–Uhlenbeck models (Felsenstein 1985; Hansen 1997; Butler and King 2004; Hansen et al. 2008; Ingram and Mahler 2013; Uyeda and Harmon 2014). Many of these models allow researchers to test hypotheses about how trait evolution relates to species’ ecology, from environmental factors (Hansen et al. 2008; Labra et al. 2009; Muschick et al. 2012; Mahler et al. 2013; Moen et al. 2015) to behavior (Scales et al. 2009; Monteiro and Nogueira 2011). Additionally, more-recent approaches have broadened the scope of these methods by allowing some consideration of species interactions (Nuismer et al. 2013; Nuismer and Harmon 2015; Drury et al. 2016; Manceau et al. 2016; Bartoszek et al. 2017), introducing Bayesian approaches that allow additional data to inform the model (Uyeda and Harmon 2014; Kostikova et al. 2016), and by allowing variation of model parameters across the tree (Beaulieu et al. 2012; Ingram and Mahler 2013; Uyeda and Harmon 2014). Herein, we explicitly incorporate empirical data on these interactions—the biological milieu in which proboscides evolve—into models of hawkmoth trait evolution. This allows us to examine both the utility of contemporary interactions in informing the evolution of this trait as well as delineating the ecological groupings, or selective regimes, where they are the most informative. In taking a hypothesis-driven approach to both the groupings of species and the values of parameters within the model, we explore the degree to which flowering plants and their contemporary floral traits can describe the proboscis-length evolution of their hawkmoth pollinators. If contemporary interactions between plants and hawkmoths do in fact capture some degree of proboscis-length evolution, our minimum expectation is that models with explicit incorporation of hawkmoth interactions will be better supported than an interaction-free model of trait evolution. Moreover, different hypothesis-based models will almost always differ in their degree of support—based on their underlying assumptions—and thus should provide additional information about the evolution of these traits. Finally, we examine the results of an uninformed, best-fit modeling approach to assess the importance of our initial hypothesis—that ecological interactions have a discernible impact on trait evolution—in these pollination systems. Consequently, we demonstrate that there is unambiguous evidence that contemporary interactions can shed light on hawkmoth trait evolution. We show that this evidence can be found in the comparisons of specific species-interaction-based hypotheses and when comparing these scenarios to unconstrained approaches to trait evolution. Methods Hawkmoth Dataset For the purposes of this study, we put together a dataset that encompasses $$75$$ hawkmoth species from $$25$$ genera and $$92$$ plant species from $$60$$ genera based on the findings of two recent studies (Kawahara and Barber 2015; Sazatornil et al. 2016). To use these data to explore the utility of ecological interactions for modeling trait evolution, we first compiled hawkmoth–plant interactions and relevant traits from four of the quantitative pollination networks presented by (Sazatornil et al., 2016). We combined the data from these four networks so that we were left with one set of realized interactions per hawkmoth species and one value for the length of a hawkmoth’s proboscis—when a hawkmoth occurred in multiple networks these were the set of all interactions observed across networks and the mean of recorded trait values, respectively. Second, the phylogeny we used to model trait evolution was built from the time-calibrated, molecular phylogeny of hawkmoths from Kawahara and Barber (2015). Kawahara and Barber (2015) used five nuclear loci and a single mitochondrial locus that totaled 7449 bp to examine the evolutionary relationships between hawkmoth species. This phylogeny included all but one of the genera present in our dataset (a single species in our data) and thus provided, at a minimum, genus-level divergence dates for 74 out of 75 species. Below the genus level, we collapsed intra-genus dates into polytomies as not all congeners in our dataset were present in the phylogeny which prevented us from reliably inferring intra-genus relationships between taxa. The presence of unresolved polytomies in this phylogeny is a limitation of the data in this study. Descriptive analyses of the data. Given that our focus in this study was to examine how well the evolution of hawkmoth proboscis length may be explained by the relevant traits (i.e., corolla length) of the plants they pollinate, we also described each hawkmoth species in one additional way. We calculated the effective interaction trait (EIT) of the plants that each hawkmoth pollinates. The EIT of each hawkmoth is the weighted mean corolla length of the set of plant species, it was observed to pollinate where weights are given by the frequency with which the hawkmoth visited each plant. Ecologically, this quantity for each hawkmoth represents the corolla length that best describes their interaction patterns and interaction preferences. In the remainder of the text, we will focus on analyses and results for hawkmoth EIT so defined; however, our results do not differ qualitatively if we use an unweighted EIT based solely on partner identity (see Supplementary Material S1 available on Dryad at http://dx.doi.org/10.5061/dryad.3vb73). To get an initial indication of the relationship between hawkmoth traits and phylogenetic relatedness, we estimated the level of phylogenetic signal of both proboscis length and EIT using Blomberg’s K (Blomberg et al. 2003). In each case, we tested whether the observed phylogenetic signal was significant with a null model comprising 10,000 shuffled trait assignments with the phytools::phylosig function in R Core Team (2013) and Revell (2012). We also examined the relationship between hawkmoth proboscis length and EIT with a phylogenetic least-squares approach (Orme 2013) to explore the extent to which the raw data captures a relationship between interactions and traits. Finally, we tested whether any morphological matching between hawkmoth proboscis length and EIT—measured as the absolute difference between the two—can be explained by the potentially confounding factors of hawkmoth lineage age (i.e., the phylogenetic branch length pertaining to a distinct species) and hawkmoth generalism (i.e., the number of interaction partners a hawkmoth has) using the same phylogenetic least-squares approach as above. Interaction-Informed Evolutionary Models of Proboscis Evolution To explore the extent to which pollination interactions can be informative descriptors of hawkmoth trait evolution, we started with a baseline evolutionary model and then constructed five additional models that explicitly incorporate both ecological interactions and specific hypotheses regarding the biological scale at which those interactions are informative. In the baseline model, traits evolve according to Brownian motion; i.e., the evolution of a trait $$X$$ over time is modeled as $$dX(t) = \sigma d\beta(t)$$, where $$d\beta (t)$$ is a random deviate and $$\sigma$$ is the volatility of that deviation. Herein, a pure-drift process acts and all species evolve independently from each other. The resulting differences in traits between species are a product of their phylogenetic relatedness (Felsenstein 1985; Butler and King 2004). We fit this baseline model with the ouch::brown function in (R Core Team 2013) and (King and Butler, 2009). In contrast to this baseline model, all of our interaction-informed models are versions of an Ornstein–Uhlenbeck process where trait evolution proceeds in the presence of adaptive peaks and potential selection towards such peaks (Hansen 1997; Butler and King 2004). The most simple Ornstein–Uhlenbeck model of trait evolution takes the form $$dX(t) = \alpha[\theta - X(t)]dt + \sigma d \beta(t),$$ (1) where $$X$$ is the trait of interest, $$\theta$$ is the evolutionarily optimal trait value, $$\alpha$$ is the rate of attraction to the optimum, and $$d \beta(t)$$ and $$\sigma$$ are as before (i.e., the components of the Brownian model). The key difference between Equation 1 and our interaction-informed models relates to the optimal trait value $$\theta$$. To address the question of how a hawkmoth’s current interactions may reflect the evolution of its proboscis length, we wanted to test hypotheses about specific values of $$\theta$$ that can be calculated from the traits of the plants they pollinate. In contrast, previous approaches maintain $$\theta$$ as a free parameter that is also estimated during model fitting with a Bayesian or Maximum Likelihood approach (Butler and King 2004; Mahler et al. 2013; Uyeda and Harmon 2014; Khabbazian et al. 2016). Similar to these previous approaches, we also group species on the phylogeny into selective regimes (i.e., groups of species hypothesized to be evolving towards the same evolutionary optima). In current approaches, each selective regime is associated with an estimated optimum trait value (Butler and King 2004; Mahler et al. 2013; Uyeda and Harmon 2014). Rather than examining how latent optima can be best-fit given a phylogenetic tree, trait distribution, and hypothesized selective regimes, our approach pre-specifies $$\theta$$ values for each selective regime directly from the empirical data (in this case, data on species interactions). To implement this model, we used a variation on the classic Ornstein–Uhlenbeck model—known as the Hansen (1997) model —and drew on the implementation of Butler and King (2004) and King and Butler (2009). All R code and data used in these analyses are freely available on Data Dryad (doi:10.5061/dryad.3vb73). Building on this approach, we calculated each group $$g$$’s optimal value for proboscis length ($$\theta_{g}$$) from its constituent hawkmoths’ pollination interactions and the traits of their flowering-plant interaction partners (i.e., an EIT value for each group). As with the simplest Ornstein–Uhlenbeck models, our model also included an ancestrally optimal trait value ($$\theta_N$$; estimated during model fitting) that allowed for the fact that contemporary pollination interactions, and the optima based on them, may not have exerted persistent selection across the full depth of the phylogeny. Finally, we included a dummy variable $$\omega_{i}(t)$$ for each species $$i$$ in group $$g$$ that governed its evolution towards either $$\theta_{g}$$ or $$\theta_{N}$$ at any point in time. Species $$i$$’s trait evolution is attracted towards the null optimum ($$\theta_N$$) when $$\omega_{i}(t) = 0$$ and towards $$\theta_{g}$$ when $$\omega_{i}(t) = 1$$. Our implementation of the Hansen (1997) model thus takes the form: $$dX_i (t) = \alpha \left[ \left(1 - \omega_i(t)\right)\theta_N + \omega_i(t)\, \theta_{g} -X_i(t) \right] + \sigma d \beta(t).$$ (2) As mentioned, we used this model to test ecological hypotheses regarding the specific values of evolutionary optima for proboscis length based on hawkmoth interactions. While our overall objective was to assess how well evolutionary optima derived from contemporary interactions can describe past trait evolution, we also leveraged this approach to address a related question regarding the biological scale at which these interactions are most informative. To do this, we compared and contrasted several sets of hypothesized selective regimes and optima. These sets fall broadly into two categories—taxonomic and functional—both of which are important ecological groupings of pollinators (Olesen et al. 2007; Rezende et al. 2009; Gómez et al. 2010) and have allowed us to assess whether proboscis-length evolution is better described by hawkmoth interactions at one or both scales. For example, to test the hypothesis that the evolution of hawkmoth proboscis length is best explained by pollination interactions at the genus level, we took the following four steps. First, we estimated the interaction-based optimum for each genus as the EIT of all hawkmoth species in that genus (where weights are given by the set of hawkmoth–plant interactions of members of the genus). In the model, the EIT of a genus represents the evolutionarily optimal proboscis length for that genus (i.e., the trait value that members of said genus will evolve towards when $$\omega_{i}(t) = 1$$); the genus itself constitutes the selective regime—a group of taxa that evolve towards the same $$\theta_{g}$$ (Fig. 1a–d). Second, we fixed the point in time at which selection toward these genus-specific optima can begin; that is, the time $$t_i^*$$ at which $$\omega_{i}(t > t_i^*) = 1$$ (Fig. 1e). Note that a hawkmoth family composed of three genera would have three different evolutionary optima in this model (i.e., one specific to each genus). Therefore, a hawkmoth species in genus A could not experience selection toward the optimal trait value of its genus until the point in time at which it diverges from genera B and C. Until that point, all species evolve according to the ancestral optimum $$\theta_N$$ (i.e. $$\omega_{i}(t < t_i^*) = 0$$). We then fit the model—with its set of hypothesized optima $$\{\theta_{g}\}$$ and appropriate $$\omega_i(t)$$ and $$t_i^*$$ values—with a maximum likelihood approach to find the best-fit values of the three remaining parameters ($$\alpha$$, $$\sigma$$, $$\theta_N$$). We implemented this procedure with the sbplx function (based on Rowan 1990) in the nloptr package in R Core Team (2013) and Ypma (2014). Finally, we calculated the model-selection statistic corrected Akaike Information Criterion (AICc) to facilitate comparison across hypotheses (Akaike 1998). Note that the number of $$\theta_{g}$$ does not have an impact on AICc as these parameters are not estimated in model-fitting (Burnham and Anderson 2001). Figure 1. View largeDownload slide Conceptual depiction of the process of how we determine hypothetical selective regimes and their corresponding optima for our interaction-informed Ornstein–Uhlenbeck model. First, we combine (a) data on pollination interactions between hawkmoths and plants with (b) observed plant corolla lengths to generate (c) an effective interaction trait (EIT) for each hawkmoth. In (a), interactions are represented by lines joining a hawkmoth and plant where line thickness represents interaction intensity. In (c), we show the EIT for each hawkmoth and they are grouped into selective regimes based on a hypothesis of how species evolve relative to each other. Herein, hawkmoth species are assigned to three functional selective regimes (orange, purple, and green) based on their EIT values. d) For each selective regime, an optimum trait value is calculated as the mean EIT of the group. e) Finally, the selective regimes and their optima are painted onto the phylogeny and the model parameters can be fit to determine the degree to which those selective regimes and optima explain observed variation in hawkmoth proboscis length. Figure 1. View largeDownload slide Conceptual depiction of the process of how we determine hypothetical selective regimes and their corresponding optima for our interaction-informed Ornstein–Uhlenbeck model. First, we combine (a) data on pollination interactions between hawkmoths and plants with (b) observed plant corolla lengths to generate (c) an effective interaction trait (EIT) for each hawkmoth. In (a), interactions are represented by lines joining a hawkmoth and plant where line thickness represents interaction intensity. In (c), we show the EIT for each hawkmoth and they are grouped into selective regimes based on a hypothesis of how species evolve relative to each other. Herein, hawkmoth species are assigned to three functional selective regimes (orange, purple, and green) based on their EIT values. d) For each selective regime, an optimum trait value is calculated as the mean EIT of the group. e) Finally, the selective regimes and their optima are painted onto the phylogeny and the model parameters can be fit to determine the degree to which those selective regimes and optima explain observed variation in hawkmoth proboscis length. We followed the above procedure for five models of proboscis-length evolution: a single optimum for all species (i.e., all hawkmoths in the same selective regime; referred to as OU1), two models based on taxonomic groupings, and two models based on functional groupings. The taxonomic models assessed how well hawkmoth interactions at the genus- and species-specific scales (i.e., separate $$\theta_{g}$$ for each genus or species, respectively) explain observed hawkmoth proboscis lengths (OU2 and OU3). The functional models aimed to assess how hawkmoth interactions can explain observed traits when hawkmoths are grouped based on the similarity of their interactions at both a coarse scale and finer scale. In the former, hawkmoths fall into two groups based their EIT value ($$x < 38\,mm$$ and $$x > 38\,mm$$; OU4). In the latter, species were assigned to one of six groups ($$x < 15\,mm$$, $$15\,mm < x < 30\,mm$$, $$30\,mm < x < 38\,mm$$, $$38\,mm < x < 50\,mm$$, $$50\,mm < x < 70\,mm$$, and $$x > 70\,mm$$; OU5). In both cases, we chose trait bins based on natural groupings of EIT values. In the functional-group models, all values $$t_{i}^*$$ were given by the time that a species diverged from other taxa (i.e., the species’ age) since functional groupings were scattered across the phylogeny. For each model, we assessed how their hypothesized $$\theta_g$$ values related to the observed hawkmoth traits. To do so, we calculated the mean-squared error (MSE) for each model where observed hawkmoth proboscis lengths were the “observed” values and the optimum proboscis length associated to each species was the corresponding “predicted” value. The MSE for each model therefore represents how close the matching between trait and optimum value was across the tree. Closer matching between the two may explain differences in model-fit and provide insights to the selective regimes where contemporary interactions best describe traits. Comparison to Best-Fit Evolutionary Model The primary purpose of our hypothesis-based approach was to examine the degree to which incorporating specific ecological data (species interactions) into models of trait evolution can make their results more informative. Our hypothesis-based approach provides insights into the ability of contemporary interactions to capture past trait evolution and the biological scale at which that ability is maximized. However, the flip side of this approach is that the selective regimes of relatively few candidate models can realistically be biologically-justified and tested. Consequently, this makes it unlikely that any one of the few scenarios tested here exactly represents the most-likely model of proboscis-length evolution. Though this does not compromise our ability to test different hypotheses (Burnham and Anderson 2001), it does imply a lack of an upper bound to which the interaction-informed models can be compared. Identifying an upper bound can thus provide added insight into how good the fit of our best interaction-informed model is and how well its selective regimes capture patterns of trait evolution. Therefore, to further examine the results, we also aimed to identify the statistically-best model of proboscis-length evolution with the SURFACE algorithm (Ingram and Mahler 2013). In SURFACE, both species’ groups (selective regimes) and the optimal trait values that describe them are estimated by the model (Ingram and Mahler 2013). Using an iterative AICc procedure, the method assesses different sets of optima and shifts (analogous to the $$\theta_{g}$$, $$\omega_i$$, and $$t_i^*$$ in our model), calculates the AICc for the model with those parameter sets, and converges on the best-fit model—as defined by the smallest AICc value (Ingram and Mahler 2013). SURFACE has previously been used to successfully identify best-fit OU models and convergent evolution in other taxa (Mahler et al. 2013; Grundler and Rabosky 2014; Almécija et al. 2015). SURFACE is also just one of multiple, related methods that we could have used to make this comparison; for brevity, we focus on its comparison here in the main text but also compare the results of our approach to SLOUCH (Hansen et al. 2008) in the Supplementary Material available on Dryad. We also assessed whether or not a fingerprint of the interactions that were made explicit in our informed models could also be identified in the SURFACE results. Here, we used a linear model to explore the relationship between the best-fit optima of the SURFACE model and hawkmoth EIT. If hawkmoth interactions are reflected in the SURFACE optima, we would expect to see a positive relationship between the two. To assess the likelihood that the strength of the observed relationship between EIT and optima may be observed by chance, we also compared the observed relationship to a null distribution of 9999 random predictor test statistics. Each value in the null model was found by randomizing a hawkmoth’s pollination interactions, recalculating EIT, and then re-examining the relationship between EIT and optima with the linear model. To generate the randomized EIT values, we shuffled the hawkmoth-plant interaction network such that a hawkmoth’s interaction partners may change but overall interaction frequency and degree (i.e., number of interaction partners) are maintained (algorithm swsh_both_c in the vegan::nullmodel function; Oksanen et al. 2016). As our hypothesis suggests that the observed relationship between EIT and optima should be stronger than expected by chance, we performed a one-tailed test ($$\alpha = 0.05$$). Next, we examined the results of the SURFACE model with respect to the groupings of hawkmoths used in our interaction-informed models. The degree to which our hypothesized groupings are captured by the best-fit model can indicate how well a basic assumption about the evolution of a trait can approximate the best statistical groupings. We first quantified how much the selective regimes of the SURFACE model overlap with each of our interaction-informed models by measuring their mutual information (Danon et al. 2005). As before, we sought to compare the observed mutual information to what might be expected by chance, so we shuffled the SURFACE regime assignments 1000 times across the hawkmoth species and reassessed the mutual information. We again performed a one-tailed test and hence considered the observed overlap significantly different than expected by chance if the observed mutual information was greater than the null distribution ($$\alpha = 0.05$$). We also compared the SURFACE selective regimes to our six-functional-group model to see how these two scenarios captured biological groupings of hawkmoths. First, we compared the phylogenetic signal of each model’s regimes to assess the degree to which hawkmoth relatedness predicted selective regime co-membership. To do so, we converted group participation into a binary design matrix upon which we were able to estimate phylogenetic signal as the multivariate K statistic of (Adams, 2014) with the function geomorph::physignal and evaluate it based on the Brownian expectation $$K = 1$$ (Adams and Otárola-Castillo 2013). Second, we estimated the degree to which the regimes represent groupings of morphologically similar hawkmoths. To do so, we examined the variation of hawkmoth proboscis length within selective regimes. We used an ANOVA to compare hawkmoth proboscis length across selective regimes, where a larger F-statistic—a relatively larger amount of variation in proboscis length explained by selective regime—would suggest tighter functional groupings of hawkmoths. Results In this group of hawkmoths, we observed distinct and significant phylogenetic signal in their proboscis lengths ($$K = 1.169$$, $$P < 0.001$$, $$n = 10000$$), suggesting that more closely related hawkmoths tend to have more similar feeding traits (Fig. 2a). We saw a similar, but weaker, pattern when we expanded our focus to include the plants those hawkmoths visit. Phylogenetic signal of hawkmoth EIT was lower than what we observed for proboscis length but still significant ($$K = 0.412$$, $$P = 0.008$$, $$n = 10000$$; Fig. 2a). On the other hand, there was a significant positive relationship between hawkmoth proboscis length and EIT (PGLS: $$\beta = 0.501$$, $$t = 3.386$$, $$P = 0.001$$, $$R^2 = 0.136$$; Fig. 2b) suggesting, as expected, that hawkmoths with longer proboscides tend to visit longer flowers. Interestingly, the degree of morphological matching we observed between an hawkmoth’s proboscis length and EIT is not explained by either the age of the hawkmoth lineage on our phylogeny or its number of observed interaction partners (PGLS: $$\beta = -0.095$$, $$t = -0.187$$, $$P = 0.852$$, $$R^2 < 0.001$$ & $$\beta = 0.287$$, $$t = 1.479$$, $$P = 0.143$$, $$R^2 = 0.030$$, respectively). Figure 2. View largeDownload slide Relationship between hawkmoth phylogeny, proboscis length (PL; mm), and effective interaction trait (EIT; mm). a) There is distinct phylogenetic signal of proboscis length compared with EIT ($$K = 1.169$$, $$P < 0.001$$, $$n = 10000$$ and $$K = 0.412$$, $$P = 0.008$$, $$n = 10000$$, respectively; where $$K > 1$$ indicates phylogenetic clustering). b) There is a significant, positive relationship between proboscis length and EIT when accounting for phylogenetic relatedness ($$t = 3.386$$, $$P = 0.001$$, $$R^2 = 0.136$$). The dashed grey line in (b) represents the fitted PGLS relationship between proboscis length and EIT where $$\beta = 0.501$$. Figure 2. View largeDownload slide Relationship between hawkmoth phylogeny, proboscis length (PL; mm), and effective interaction trait (EIT; mm). a) There is distinct phylogenetic signal of proboscis length compared with EIT ($$K = 1.169$$, $$P < 0.001$$, $$n = 10000$$ and $$K = 0.412$$, $$P = 0.008$$, $$n = 10000$$, respectively; where $$K > 1$$ indicates phylogenetic clustering). b) There is a significant, positive relationship between proboscis length and EIT when accounting for phylogenetic relatedness ($$t = 3.386$$, $$P = 0.001$$, $$R^2 = 0.136$$). The dashed grey line in (b) represents the fitted PGLS relationship between proboscis length and EIT where $$\beta = 0.501$$. Given this baseline picture of the interplay between hawkmoth proboscis length, interactions, and evolutionary history, it was still unclear how informative ecological interactions might ultimately be for describing proboscis-length evolution or the biological scale at which this would be maximal. Our models, however, suggested that trait evolution was best captured by contemporary ecological interactions when hawkmoths are grouped at a fine functional scale or on their own as individual species. Indeed, only models with optima that were fitted at the species (OU3) and six-functional-group scale (OU5) performed better than Brownian motion at explaining proboscis-length evolution (Table 1). In contrast, the models with a single global optimum (the collective EIT of all hawkmoths in the phylogeny), with genus-specific optima, and with two functional groups all had higher AICc values—and were therefore poorer candidate models—than the Brownian-motion model (Table 1). In terms of AICc, the six-functional-group model and the species-specific model were indistinguishable (OU3 and OU5; Table 1). Table 1 Summary of model fit for each scenario for hawkmoth proboscis evolution Model Type $$\alpha$$ $$\sigma$$ $$\theta_{N}$$ nSR lnL $$n$$ AICc MSE BM — — 25.240 40.604 — $$-$$306.515 2 617.196 — OU1 Global 0.468 27.391 43.245$$^{\mathrm{\dagger}}$$ 1 $$-$$307.591 2 621.521 590.562 OU2 Tax. 0.030 24.770 43.329 25 $$-$$306.935 3 620.208 349.916 OU3 Tax. 0.025 23.367 28.410 75 $$-$$302.702 3 611.742 404.395 OU4 Func. 0.012 25.177 17.840 2 $$-$$308.856 3 624.050 477.653 OU5 Func. 0.026 23.427 28.993 6 $$-$$302.855 3 612.048 395.813 SF — 0.069 2.856 27.599* 7 $$-$$253.744 9 528.258 119.258 Model Type $$\alpha$$ $$\sigma$$ $$\theta_{N}$$ nSR lnL $$n$$ AICc MSE BM — — 25.240 40.604 — $$-$$306.515 2 617.196 — OU1 Global 0.468 27.391 43.245$$^{\mathrm{\dagger}}$$ 1 $$-$$307.591 2 621.521 590.562 OU2 Tax. 0.030 24.770 43.329 25 $$-$$306.935 3 620.208 349.916 OU3 Tax. 0.025 23.367 28.410 75 $$-$$302.702 3 611.742 404.395 OU4 Func. 0.012 25.177 17.840 2 $$-$$308.856 3 624.050 477.653 OU5 Func. 0.026 23.427 28.993 6 $$-$$302.855 3 612.048 395.813 SF — 0.069 2.856 27.599* 7 $$-$$253.744 9 528.258 119.258 Notes: $$\alpha$$ = strength of attraction; $$\sigma$$ = magnitude of variation; $$\theta_{N}$$ = intrinsic optimum; nSR = number of selective regimes in model (in interaction-informed models this is also the number of $$\theta_{g}$$ values fixed); lnL = log-likelihood; n = number of free parameters; AICc = small-sample-size-corrected Akaike Information Criterion; MSE = mean-squared error between $$\theta_{g}$$ values and observed traits; BM = Brownian motion; OU1 = global-optimum; OU2 = genus-specific optima; OU3 = species-specific optima; OU4 = two-functional-group-based optima; OU5 = six-functional-group-based optima; SF = SURFACE implementation. “Tax.” refers to taxonomic groupings of species, “Func.” refers to groupings based on similarity of interactions, “Global” refers to no groupings. A “—” indicates that the value is not applicable to the model, a “†” indicates that $$\theta_{N}$$ was equal to the fitted optima, and the “*” after the $$\theta_{N}$$ value for the SURFACE model refers to the estimated $$\theta$$ value at the root. Table 1 Summary of model fit for each scenario for hawkmoth proboscis evolution Model Type $$\alpha$$ $$\sigma$$ $$\theta_{N}$$ nSR lnL $$n$$ AICc MSE BM — — 25.240 40.604 — $$-$$306.515 2 617.196 — OU1 Global 0.468 27.391 43.245$$^{\mathrm{\dagger}}$$ 1 $$-$$307.591 2 621.521 590.562 OU2 Tax. 0.030 24.770 43.329 25 $$-$$306.935 3 620.208 349.916 OU3 Tax. 0.025 23.367 28.410 75 $$-$$302.702 3 611.742 404.395 OU4 Func. 0.012 25.177 17.840 2 $$-$$308.856 3 624.050 477.653 OU5 Func. 0.026 23.427 28.993 6 $$-$$302.855 3 612.048 395.813 SF — 0.069 2.856 27.599* 7 $$-$$253.744 9 528.258 119.258 Model Type $$\alpha$$ $$\sigma$$ $$\theta_{N}$$ nSR lnL $$n$$ AICc MSE BM — — 25.240 40.604 — $$-$$306.515 2 617.196 — OU1 Global 0.468 27.391 43.245$$^{\mathrm{\dagger}}$$ 1 $$-$$307.591 2 621.521 590.562 OU2 Tax. 0.030 24.770 43.329 25 $$-$$306.935 3 620.208 349.916 OU3 Tax. 0.025 23.367 28.410 75 $$-$$302.702 3 611.742 404.395 OU4 Func. 0.012 25.177 17.840 2 $$-$$308.856 3 624.050 477.653 OU5 Func. 0.026 23.427 28.993 6 $$-$$302.855 3 612.048 395.813 SF — 0.069 2.856 27.599* 7 $$-$$253.744 9 528.258 119.258 Notes: $$\alpha$$ = strength of attraction; $$\sigma$$ = magnitude of variation; $$\theta_{N}$$ = intrinsic optimum; nSR = number of selective regimes in model (in interaction-informed models this is also the number of $$\theta_{g}$$ values fixed); lnL = log-likelihood; n = number of free parameters; AICc = small-sample-size-corrected Akaike Information Criterion; MSE = mean-squared error between $$\theta_{g}$$ values and observed traits; BM = Brownian motion; OU1 = global-optimum; OU2 = genus-specific optima; OU3 = species-specific optima; OU4 = two-functional-group-based optima; OU5 = six-functional-group-based optima; SF = SURFACE implementation. “Tax.” refers to taxonomic groupings of species, “Func.” refers to groupings based on similarity of interactions, “Global” refers to no groupings. A “—” indicates that the value is not applicable to the model, a “†” indicates that $$\theta_{N}$$ was equal to the fitted optima, and the “*” after the $$\theta_{N}$$ value for the SURFACE model refers to the estimated $$\theta$$ value at the root. We saw that MSE of proboscis length and $$\theta_g$$ somewhat corresponded to the pattern of model-fit. In the species-specific (OU3) and six-functional group model (OU5), matching between the two was close while it was most discordant at the global scale (OU1; Table 1). Most interestingly, the model with genus-specific optima (OU2) showed better matching than expected from its model fit as its MSE was the lowest of our interaction-informed models despite it showing a fit worse than Brownian Motion based on both $$lnL$$ and AICc (OU2; Table 1). While some of our interaction-informed models provided an improved description of hawkmoth proboscis-length evolution relative to Brownian motion, the unconstrained model of proboscis-length evolution (identified by SURFACE) fit substantially better than our best interaction-informed model, as expected ($$lnL = -253.744$$, $$\alpha = 0.069$$, $$\sigma = 2.856$$, $$AICc = 528.258$$, $$\#\theta = 7$$; Table 1). The MSE of optima and traits for this model was also the lowest of all the models that we fit (Table 1). There was a clear positive relationship between hawkmoth EIT and optima for the SURFACE model ($$\beta = 0.436$$, $$t = 7.018$$, $$P < 0.001$$, $$R^2 = 0.403$$; Fig. 3a). The observed hawkmoth EIT values showed a stronger relationship with the SURFACE optima than any of the models of the null distribution (Fig. 3b). The positive relationship between EIT and optima suggested that the regimes and optima upon which SURFACE converges are at least correlated with hawkmoth interactions. Perhaps unsurprisingly, this relationship is weaker than the same for the six-functional-group model where the optima were drawn from the EIT values of each functional group. Figure 3. View largeDownload slide The relationship between a hawkmoth species’ optimal proboscis length (PL) and its observed effective interaction trait (EIT). a) There is a significant, positive correlation between a species’ EIT and its optimal proboscis length as determined by SURFACE. The solid line shows the slope from a linear regression of the values while the dashed line shows the same for our best-fitting interaction-informed model for comparative purposes. The observed slope between EIT and SURFACE optima is smaller ($$\beta= 0.44$$) than that from EIT and the optima of the interaction-informed model, where optima are based on EIT values ($$\beta = 1.01$$; dashed line). b) We show the degree to which the observed relationship between EIT and SURFACE optima is stronger than expected by chance. The vertical line indicates the observed t-value and the distribution represents 9999 null models in which hawkmoth EIT is shuffled between species and the relationship between SURFACE optima and EIT reassessed. Figure 3. View largeDownload slide The relationship between a hawkmoth species’ optimal proboscis length (PL) and its observed effective interaction trait (EIT). a) There is a significant, positive correlation between a species’ EIT and its optimal proboscis length as determined by SURFACE. The solid line shows the slope from a linear regression of the values while the dashed line shows the same for our best-fitting interaction-informed model for comparative purposes. The observed slope between EIT and SURFACE optima is smaller ($$\beta= 0.44$$) than that from EIT and the optima of the interaction-informed model, where optima are based on EIT values ($$\beta = 1.01$$; dashed line). b) We show the degree to which the observed relationship between EIT and SURFACE optima is stronger than expected by chance. The vertical line indicates the observed t-value and the distribution represents 9999 null models in which hawkmoth EIT is shuffled between species and the relationship between SURFACE optima and EIT reassessed. We quantified the similarity of selective regimes between SURFACE and each of our models as the mutual information shared by the two groupings of hawkmoths into selective regimes. The regimes of SURFACE and the six-functional-group model have mutual information $$MI = 0.232$$ (Monte Carlo test; $$P < 0.004$$), where $$1$$ indicates perfect overlap and $$0$$ is no overlap. Interestingly, when comparing the SURFACE regimes to our other candidate interaction-informed models, the genus-wise model showed a higher mutual information than our best-fitting models ($$MI = 0.666$$, $$P < 0.001$$) while the global ($$MI = 0.000$$, $$P = 1.00$$) and two-functional-group ($$MI = 0.144$$, $$P = 0.324$$) showed lower values. The species-specific model showed a higher value than the six-functional-group model ($$MI = 0.544$$, $$P < 0.001$$). Although both our species-specific and six-functional-group scenarios performed as well as each other, we focused on the SURFACE and six-functional-group models when comparing the similarities of selective regimes as comparisons to the species-specific scenario would be uninformative. We saw that the SURFACE selective regimes show significantly greater phylogenetic signal than expected under Brownian motion ($$K = 1.518$$, $$P < 0.001$$, $$n = 10000$$; Supplementary Material S2 available on Dryad) as opposed to the regimes of the six-functional-group model which showed significantly less ($$K = 0.289$$, $$P = 0.040$$, $$n = 10000$$; Fig. 4a). Moreover, the SURFACE regimes also represented tighter morphological groupings of hawkmoths (in terms of the degree to which selective regime explains proboscis-length variation) than the six-functional-group model ($$F = 91.76$$, $$P < 0.001$$ & $$F = 10.4$$, $$P < 0.001$$, respectively). These results suggest that the SURFACE regimes captured taxonomic and morphological similarity between hawkmoths which is not captured by the functional groupings that we define based on hawkmoth EIT values and might explain the discrepancy in fit between the two models (Table 1). Figure 4. View largeDownload slide Summary of the more informative of our best-fitting interaction-informed models: the six-functional-group Ornstein–Uhlenbeck model (OU5). a) The phylogeny of hawkmoths in our dataset. Black lines indicate the periods in which selection acted toward the ancestral optimum $$\theta_N$$ and the different colors represent distinct functional groups and the distinct selective regimes. b) The hawkmoth trait space encompassed by each selective regime. Each distribution is made up of the effective interaction trait (EIT) values of hawkmoths in each of the six selective regimes. c) The optimal hawkmoth proboscis length for each selective regime that we used in the model; these values are the mean EIT of each functional group. Figure 4. View largeDownload slide Summary of the more informative of our best-fitting interaction-informed models: the six-functional-group Ornstein–Uhlenbeck model (OU5). a) The phylogeny of hawkmoths in our dataset. Black lines indicate the periods in which selection acted toward the ancestral optimum $$\theta_N$$ and the different colors represent distinct functional groups and the distinct selective regimes. b) The hawkmoth trait space encompassed by each selective regime. Each distribution is made up of the effective interaction trait (EIT) values of hawkmoths in each of the six selective regimes. c) The optimal hawkmoth proboscis length for each selective regime that we used in the model; these values are the mean EIT of each functional group. Discussion Altogether, our results provide a strong indication that contemporary ecological interactions can improve our understanding of past evolution of hawkmoth proboscis length. In demonstrating this link, we contribute an alternative approach to modeling trait evolution where the consideration of ecological interactions, or other aspects of natural history, can be built directly into the model parameters (see also Hansen et al. 2008; Uyeda and Harmon 2014; Kostikova et al. 2016) and potentially provide a more complete story of trait evolution. Furthermore, the biological realism incorporated into our model—in the form of ecological interactions—suggests that the evolution of proboscis length in this group of hawkmoths is explained by the contemporary pollination interactions of hawkmoths that interact with similar plants just as well as it is when each species and their interactions are considered separately. Lastly, our finding that contemporary interactions can be informative facets of trait evolution is supported by the fingerprint of these interactions that can be identified in an uninformed, unconstrained, and statistically best model. That interactions between species are predicated on their possession of traits selected to facilitate or avoid that interaction is a basic tenet of evolutionary ecology. Testing the assumption—that contemporary traits arose from the feedback of interaction-based selection—in the absence of historical data is difficult to impossible (however, see Gervasi and Schiestl 2017). Modeling evolutionary hypotheses is one approach to test questions and assumptions, such as these, that cannot be manipulated empirically (Butler and King 2004; Ingram et al. 2012; Mahler et al. 2013; Uyeda and Harmon 2014; Moen et al. 2015; Ingram et al. 2016). Our results show that, for a group of New World hawkmoths, contemporary pollination interactions can provide an informative lens through which the evolution of relevant pollination traits can be examined. In the dataset that we examine here, we see that hawkmoth traits correlate with their interactions even if one-to-one matching between proboscis and EIT is not observed. Intriguingly, our results suggest that matching between interactions and proboscis length is poor on a species-by-species scale since the interactions of functionally similar groups provide just as close of a match between trait and interactions as the species themselves do. These discrepancies may arise from the right-skewed trait distributions for pollinators in this data (Sazatornil et al. 2016) and others (Stang et al. 2009), but also deserves further attention. The biological scale at which traits evolve concordantly (i.e., converge) is an interesting question in its own right. For the hawkmoths studied here, we see that it is just as good to assume that interactions explain the evolution of functionally similar hawkmoth species as it is to assume that each hawkmoth is unique. Additionally, when these functional groups are too coarse-grained they can prove no more informative than a Brownian motion model. A result such as this might seem obvious based on the idea that species with similar interaction patterns should possess similar traits (Jordano 1995; Fenster et al. 2004; Olesen et al. 2007). However, phylogenetic conservation of interactions is a widely investigated phenomenon in evolutionary ecology (Rezende et al. 2007, 2009; Gómez et al. 2010) which would indicate that ecological interactions may well have captured evolution at a taxonomic scale such as the genus-scale that we examine. An hypothesis of phylogenetic conservation would also be supported by the phylogenetic signal of hawkmoths proboscis length and EIT that we see in this study. Our results show identical support between a finer-scale functional approach and species-specific approach which may question the role of taxonomic patterns in a wider evolutionary context and beggar further investigation. On the other hand, the best model from SURFACE, which fit better than any of our models, seemed to capture both taxonomic and morphological similarity in its selective regimes. Altogether, this suggests that future hypotheses of how species evolve in relation to each other require more nuance to disentangle the interplay between morphology and ancestry when it comes to the role of species interactions in directing adaptation. Models of trait evolution are rapidly moving on from Brownian motion and the basic Ornstein-Uhlenbeck approaches (Butler and King 2004; Hansen et al. 2008; Ingram and Mahler 2013; Uyeda and Harmon 2014; Moen et al. 2015; Khabbazian et al. 2016). However, as methods improve and computational power increases, it is important to consider the biological insight that can be drawn out of optimization-heavy approaches for modeling adaptation by natural selection. While these approaches can undeniably identify the statistically-best result—as demonstrated here—the biological implications of such results can be obscured by the variability of latent optima (Ho and Ané 2014). Spurious parameter estimates can arise due to flat likelihood surfaces that are common in evolutionary models (Ho and Ané 2014). At the same time, there is a movement to do a better job of incorporating species interactions into models of trait evolution (Nuismer and Harmon 2015; Drury et al. 2016; Manceau et al. 2016; Bartoszek et al. 2017). Our results use empirical data to demonstrate that even small alterations to established models can allow the incorporation of ecological interactions, change the way trait evolution modeling is approached, and provide insights into hypothetical parameter values that are rooted in biological reality. The obvious comparison within the scope of this study is the difference between the two modeling approaches. In the approach that we introduce—where empirical data is used to inform model parameters in a rigid manner—the potential biological relevance of the model is explicit. Hence this is an alternative to the recently proposed Bayesian approaches which can also consider additional data in the form of informative priors (Uyeda and Harmon 2014; Kostikova et al. 2016). On the other hand, the link between trait evolution and hawkmoth pollination interactions in the SURFACE results, despite a substantially-better fit, is only detectable a posteriori. SURFACE, the best-fitting model, and our approach, the biologically-explicit model, may well represent two ends of a spectrum, somewhere along which the best approach lies. In developing the approach that we present here, we add some ecological reality to a model of trait evolution. However, there are several areas in which this approach could be expanded beyond its current limitations. The major limitation of our approach currently is that it is rather static. First, the optima that we employ are fixed based on contemporary species interactions. We are therefore limited to the strong hypothesis that contemporary interactions are representative of—at least recent—evolutionary history even though species interactions can change on the order of weeks and months (Olesen et al. 2008; CaraDonna et al. 2017; Ponisio et al. 2017) and traits in long-lived species can evolve at quickest on the order of tens to hundreds of years (Galetti et al. 2013). Second, we only focus on one side of an interaction type that is regarded up as a textbook example of coevolution (Darwin 1862; Herre 1989; Anderson and Johnson 2008; Johnson and Anderson 2010). Presently, our approach does not consider how reciprocal selection has altered proboscis length, corolla length, and the interactions themselves through evolutionary time. To bring additional realism into an approach such as ours, a future adaptation could simultaneously model the evolution of traits on both sides of an interaction where optima on each side are updated through time based on how the traits of their interaction partners evolve. This would also require a realistic model of how interactions rewire through time but would provide a new approach to studying coevolution. Undeniably, species interactions are a foundational feature of ecological communities and are a source of selection pressure on trait evolution (Paine 1966; Thompson 2005). The approach we take here is amenable to the incorporation of data from a wide range of these interactions not just highly-specific mutualisms. Predator-prey relationships are ubiquitous and there are several examples of arms-races where trait change is thought to be driven by these interactions (Koskela and Ylönen 1995; Geffeney et al. 2002; Bro-Jørgensen 2013). For example, our method could be easily used to assess how hunting traits in coursing predators, such as top speed and endurance, track those same traits in their prey species. In this case, predators and their traits would be the focal group while predator EIT would be comprised of interaction data and the fleeing traits of their prey. Aside from species interactions, our approach could also be used to assess traits relating to the abiotic environment. For instance, it is common to assess optimal values or tolerance levels for temperature and nutrients in sessile organisms (Coles et al. 1976; Coles and Jokiel 1978; Ingestad 1979; Chapin III et al. 1983; Quesnel et al. 2006). These species are thought to be at risk from rapidly changing environmental conditions and their anthropogenic sources (Davis and Shaw 2001; McKenney et al. 2007; Ainsworth et al. 2016). Our method could be used to assess how groups such as corals, micro-algae, and trees may respond to these pressures. For example, species’ preferred nutrient levels or temperatures may used as the focal trait while current or previous values of the abiotic factor under consideration may be used for calculating model parameters. In all, the possibilities inherent in an interaction-informed modeling approach—where the question is explicit in the model—may represent a way in which evolutionary models can be used to test hypotheses about the links between traits and various aspects of a species’ ecology from interactions to abiotic conditions. A complete understanding of trait evolution in a group of species can require a wealth of information. Despite adding an additional layer of required data to comparative studies, we have demonstrated that knowledge of a group’s contemporary interaction partners can be useful for examining specific hypotheses about their trait evolution. Our results contribute to a growing emphasis on the need to take ecological interactions into account when examining species’ traits and their evolution (Nuismer et al. 2013; Nuismer and Harmon 2015; Manceau et al. 2016; Bartoszek et al. 2017). In particular for studies of systems that purportedly show interaction-selected morphology, our implementation of the Hansen model (Hansen 1997) may provide a new angle from which to explore the evolutionary pressure exerted on populations through their interactions with other species. Within the current biodiversity crisis (Barnosky et al. 2011; Ceballos et al. 2015), the loss of species interactions—or the appearance of novel ones—is guaranteed. Indeed, dramatic and rapid trait evolution has already been documented when interaction partners are lost (Galetti et al. 2013). If we are to grasp the consequences and understand the outcomes of the imminent—if not current—changes that ecological systems face, developing a better conceptual understanding of the interplay between ecological interaction and evolutionary adaptation will be key. Supplementary Material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.3vb73. Funding This work was supported by the University of Canterbury Summer Scholarship program and a Princeton University First Year Natural Sciences and Engineering Fellowship to M.C.H; São Paulo Research Foundation [FAPESP; 2013/13319-5 and 2014/20572-1 to M.P.G]; and Rutherford Discovery Fellowship, administered by the Royal Society of New Zealand to D.B.S. Acknowledgments The authors acknowledge the Associate Editor and three anonymous Reviewers for their valuable comments on this manuscript. Additionally, the authors acknowledge both sets of scientists whose work built the hawkmoth datasets that we incorporate as the example in this study. References Adams D.C. 2014 . A generalized K statistic for estimating phylogenetic signal from shape and other high-dimensional multivariate data. Syst. Biol. 63 : 685 – 697 . Google Scholar CrossRef Search ADS PubMed Adams D.C., Otárola-Castillo E. 2013 . geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods Ecol. Evol. 4 : 393 – 399 . Google Scholar CrossRef Search ADS Ainsworth, T.D., Heron S.F., Ortiz J.C., Mumby P.J., Grech A., Ogawa D., Eakin C.M., Leggat W. 2016 . Climate change disables coral bleaching protection on the Great Barrier Reef. Science 352 : 338 – 342 . Google Scholar CrossRef Search ADS PubMed Akaike H. 1998 . Information theory and an extension of the maximum likelihood principle. In: Selected Papers of Hirotugu Akaike. New York : Springer . p. 199 – 213 . Almécija, S., Smaers J.B., Jungers W.L. 2015 . The evolution of human and ape hand proportions. Nat. Commun. 6 : 7717 . Google Scholar CrossRef Search ADS PubMed Anderson, B., Johnson S.D. 2008 . The geographical mosaic of coevolution in a plant–pollinator mutualism. Evolution 62 : 220 – 225 . Google Scholar CrossRef Search ADS PubMed Barnosky, A.D., Matzke N., Tomiya S., Wogan G.O., Swartz B., Quental T.B., Marshall C., McGuire J.L., Lindsey E.L., Maguire K.C., Mersey B. 2011 . Has the earth’s sixth mass extinction already arrived? Nature 471 : 51 – 57 . Google Scholar CrossRef Search ADS PubMed Bartoszek, K., Glémin S., Kaj I., Lascoux M. 2017 . Using the Ornstein–Uhlenbeck process to model the evolution of interacting populations. J. Theor. Biol. 429 : 35 – 45 . Google Scholar CrossRef Search ADS PubMed Beaulieu, J.M., Jhwueng D.-C., Boettiger C., O’Meara B.C. 2012 . Modeling stabilizing selection: expanding the Ornstein–Uhlenbeck model of adaptive evolution. Evolution 66 : 2369 – 2383 . Google Scholar CrossRef Search ADS PubMed Blomberg, S.P., Garland T., Ives A.R. 2003 . Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution 57 : 717 – 745 . Google Scholar CrossRef Search ADS PubMed Bro-Jørgensen, J. 2013 . Evolution of sprint speed in african savannah herbivores in relation to predation. Evolution 67 : 3371 – 3376 . Google Scholar CrossRef Search ADS PubMed Brouat, C., Garcia N., Andary C., McKey D. 2001 . Plant lock and ant key: pairwise coevolution of an exclusion filter in an ant–plant mutualism. Proc. R. Soc. Lond. B Biol. Sci. 268 : 2131 – 2141 . Google Scholar CrossRef Search ADS Burnham, K.P., Anderson D.R. 2001 . Kullback-Leibler information as a basis for strong inference in ecological studies. Wildlife Res. 28 : 111 – 119 . Google Scholar CrossRef Search ADS Butler, M.A., King A.A. 2004 . Phylogenetic comparative analysis: a modeling approach for adaptive evolution. Am. Nat. 164 : 683 – 695 . Google Scholar CrossRef Search ADS PubMed CaraDonna, P.J., Petry W.K., Brennan R.M., Cunningham J.L., Bronstein J.L., Waser N.M., Sanders N.J. 2017 . Interaction rewiring and the rapid turnover of plant–pollinator networks. Ecol. Lett. 20 : 385 – 394 . Google Scholar CrossRef Search ADS PubMed Ceballos, G., Ehrlich P.R., Barnosky A.D., García A., Pringle R.M., Palmer T.M. 2015 . Accelerated modern human–induced species losses: Entering the sixth mass extinction. Sci. Adv. 1 : e1400253 . Google Scholar CrossRef Search ADS PubMed Chapin F.S. Tryon III, P.R., Cleve K.V. 1983 . Influence of phosphorus on growth and biomass distribution of Alaskan taiga tree seedlings. Can. J. For. Res. 13 : 1092 – 1098 . Google Scholar CrossRef Search ADS Coles, S., Jokiel P.L. 1978 . Synergistic effects of temperature, salinity and light on the hermatypic coral Montipora verrucosa. Mar. Biol. 49 : 187 – 195 . Google Scholar CrossRef Search ADS Coles, S.L., Jokiel P.L., Lewis C.R. 1976 . Thermal tolerance in tropical versus subtropical Pacific reef corals 30 : 159 – 166 . Danon, L., Diaz-Guilera A., Duch J., Arenas A. 2005 . Comparing community structure identification. J. Stat. Mech. Theory Exp. 2005 : P09008 . Darwin C. 1859 . On the Origin of the Species. London : Murray . Darwin C. 1862 . On the various contrivances by which British and foreign orchids are fertilised by insect, and the good effects of intercrossing. London, UK : John Murray . Davis, M.B., Shaw R.G. 2001 . Range shifts and adaptive responses to Quaternary climate change. Science 292 : 673 – 679 . Google Scholar CrossRef Search ADS PubMed Dehling, D.M., Töpfer T., Schaefer H.M., Jordano P., Böhning-Gaese K., Schleuning M. 2014 . Functional relationships beyond species richness patterns: trait matching in plant–bird mutualisms across scales. Glob. Ecol. Biogeogr. 23 : 1085 – 1093 . Google Scholar CrossRef Search ADS Drury, J., Clavel J., Manceau M., Morlon H. 2016 . Estimating the effect of competition on trait evolution using maximum likelihood inference. Syst. Biol. 65 : 700 – 710 . Google Scholar CrossRef Search ADS PubMed Feinsinger, P., Murray K.G., Kinsman S., Busby W.H. 1986 . Floral neighborhood and pollination success in four hummingbird-pollinated cloud forest plant species. Ecology 67 : 449 – 464 . Google Scholar CrossRef Search ADS Felsenstein J. 1985 . Phylogenies and the comparative method. Am. Nat. 125 : 1 – 15 . Google Scholar CrossRef Search ADS Fenster, C.B., Armbruster W.S., Wilson P., Dudash M.R., Thomson J.D. 2004 . Pollination syndromes and floral specialization. Annu. Rev. Ecol. Evol. Syst. 35 : 375 – 403 . Google Scholar CrossRef Search ADS Galetti, M., Guevara R., Côrtes M.C., Fadini R., Von Matter S., Leite A.B., Labecca F., Ribeiro T., Carvalho C.S., Collevatti R.G., Pires M.M., Guimarães P.R., Jr Brancalion P.H., Ribeiro M.C., Jordano P. 2013 . Functional extinction of birds drives rapid evolutionary changes in seed size. Science 340 : 1086 – 1090 . Google Scholar CrossRef Search ADS PubMed Geffeney, S., Brodie E.D., Ruben P.C. 2002 . Mechanisms of adaptation in a predator-prey arms race: Ttx-resistant sodium channels. Science 297 : 1336 – 1339 . Google Scholar CrossRef Search ADS PubMed Gervasi, D.D., Schiestl F.P. 2017 . Real-time divergent evolution in plants driven by pollinators. Nat. Commun. 8 : 14691 . Google Scholar CrossRef Search ADS PubMed Gómez, J.M., Verdú M., Perfectti F. 2010 . Ecological interactions are evolutionarily conserved across the entire tree of life. Nature 465 : 918 – 921 . Google Scholar CrossRef Search ADS PubMed Grundler, M.C., Rabosky D.L. 2014 . Trophic divergence despite morphological convergence in a continental radiation of snakes. Proc. R. Soc. Lond. B Biol. Sci. 281 : 20140413 . Google Scholar CrossRef Search ADS Hansen, T.F. 1997 . Stabilizing selection and the comparative analysis of adaptation. Evolution 51 : 1341 – 1351 . Google Scholar CrossRef Search ADS PubMed Hansen, T.F., Pienaar J., Orzack S.H. 2008 . A comparative method for studying adaptation to a randomly evolving environment. Evolution 62 : 1965 – 1977 . Google Scholar PubMed Haverkamp, A., Bing J., Badeke E., Hansson B.S., Knaden M. 2016 . Innate olfactory preferences for flowers matching proboscis length ensure optimal energy gain in a hawkmoth. Nat. Commun. 7 : 11644 . Google Scholar CrossRef Search ADS PubMed Herre, E.A. 1989 . Coevolution of reproductive characteristics in 12 species of new world figs and their pollinator wasps. Cell. Mol. Life Sci. 45 : 637 – 647 . Google Scholar CrossRef Search ADS Ho, L.S.T., Ané C. 2014 . Intrinsic inference difficulties for trait evolution with Ornstein-Uhlenbeck models. Methods Ecol. Evol. 5 : 1133 – 1146 . Google Scholar CrossRef Search ADS Ingestad, T. 1979 . Nitrogen stress in birch seedlings. Physiol. Plant. 45 : 149 – 157 . Google Scholar CrossRef Search ADS Ingram, T., Harmon L., Shurin J. 2012 . When should we expect early bursts of trait evolution in comparative data? Predictions from an evolutionary food web model. J. Evol. Biol. 25 : 1902 – 1910 . Google Scholar CrossRef Search ADS PubMed Ingram, T., Harrison A., Mahler D.L., del Rosario Castañeda, M. Glor R.E., Herrel A., Stuart Y.E., Losos J.B. 2016 . Comparative tests of the role of dewlap size in Anolis lizard speciation. Proc. Biol. Sci. 283 . pii: 20162199 . Google Scholar CrossRef Search ADS PubMed Ingram, T., Mahler D.L. 2013 . SURFACE: detecting convergent evolution from comparative data by fitting Ornstein-Uhlenbeck models with stepwise Akaike Information Criterion. Methods Ecol. Evol. 4 : 416 – 425 . Google Scholar CrossRef Search ADS Janzen, D.H., Martin P.S. 1982 . Neotropical anachronisms: the fruits the gomphotheres ate. Science 215 : 19 – 27 . Google Scholar CrossRef Search ADS PubMed Johnson, S.D., Anderson B. 2010 . Coevolution between food-rewarding flowers and their pollinators. Evol. Educ. Outreach 3 : 32 – 39 . Google Scholar CrossRef Search ADS Jordano, P. 1995 . Angiosperm fleshy fruits and seed dispersers: a comparative analysis of adaptation and constraints in plant-animal interactions. Am. Nat. 145 : 163 – 191 . Google Scholar CrossRef Search ADS Kawahara, A.Y., Barber J.R. 2015 . Tempo and mode of antibat ultrasound production and sonar jamming in the diverse hawkmoth radiation. Proc. Natl. Acad. Sci. USA 112 : 6407 – 6412 . Google Scholar CrossRef Search ADS Kessler, A., Halitschke R. 2009 . Testing the potential for conflicting selection on floral chemical traits by pollinators and herbivores: predictions and case study. Funct. Ecol. 23 : 901 – 912 . Google Scholar CrossRef Search ADS Khabbazian, M., Kriebel R., Rohe K., Ané C. 2016 . Fast and accurate detection of evolutionary shifts in Ornstein–Uhlenbeck models. Methods Ecol. Evol. 7 : 811 – 824 . Google Scholar CrossRef Search ADS King, A., Butler M. 2009 . ouch: Ornstein-Uhlenbeck models for phylogenetic comparative hypotheses. R Package Version 2. Available from: https://cran.r-project.org/web/packages/ouch/index.html. Koskela, E., Ylönen H. 1995 . Suppressed breeding in the field vole (Microtus agrestis): an adaptation to cyclically fluctuating predation risk. Behav. Ecol. 6 : 311 – 315 . Google Scholar CrossRef Search ADS Kostikova, A., Silvestro D., Pearman P.B., Salamin N. 2016 . Bridging inter-and intraspecific trait evolution with a hierarchical Bayesian approach. Syst. Biol. 65 : 417 – 431 . Google Scholar CrossRef Search ADS PubMed Labra, A., Pienaar J., Hansen T.F. 2009 . Evolution of thermal physiology in liolaemus lizards: adaptation, phylogenetic inertia, and niche tracking. Am. Nat. 174 : 204 – 220 . Google Scholar CrossRef Search ADS PubMed Machado, C.A., Robbins N., Gilbert M.T.P., Herre E.A. 2005 . Critical review of host specificity and its coevolutionary implications in the fig/fig-wasp mutualism. Proc. Natl. Acad. Sci. USA 102 : 6558 – 6565 . Google Scholar CrossRef Search ADS Mahler, D.L., Ingram T., Revell L.J., Losos J.B. 2013 . Exceptional convergence on the macroevolutionary landscape in island lizard radiations. Science 341 : 292 – 295 . Google Scholar CrossRef Search ADS PubMed Manceau, M., Lambert A., Morlon H. 2016 . A unifying comparative phylogenetic framework including traits coevolving across interacting lineages. Syst. Biol. 66 : 551 – 568 . McKenney, D.W., Pedlar J.H., Lawrence K., Campbell K., Hutchinson M.F. 2007 . Potential impacts of climate change on the distribution of North American trees. AIBS Bull. 57 : 939 – 948 . Miklashevichs, E., Röhrig H., Schell J., Schmidt J. 2001 . Perception and signal transduction of rhizobial NOD factors. Crit. Rev. Plant Sci. 20 : 373 – 394 . Google Scholar CrossRef Search ADS Moen, D.S., Morlon H., Wiens J.J. 2015 . Testing convergence versus history: convergence dominates phenotypic evolution for over 150 million years in frogs. Syst. Biol. 65 : 146 – 160 . Google Scholar CrossRef Search ADS PubMed Monteiro, L.R., Nogueira M.R. 2011 . Evolutionary patterns and processes in the radiation of phyllostomid bats. BMC Evol. Biol. 11 : 137 . Google Scholar CrossRef Search ADS PubMed Muschick, M., Indermaur A., Salzburger W. 2012 . Convergent evolution within an adaptive radiation of cichlid fishes. Curr. Biol. 22 : 2362 – 2368 . Google Scholar CrossRef Search ADS PubMed Nilsson L.A. 1988 . The evolution of flowers with deep corolla tubes. Nature 334 : 147 – 149 . Google Scholar CrossRef Search ADS Niven, J.E., Laughlin S.B. 2008 . Energy limitation as a selective pressure on the evolution of sensory systems. J. Exp. Biol. 211 : 1792 – 1804 . Google Scholar CrossRef Search ADS PubMed Nuismer, S.L., Harmon L.J. 2015 . Predicting rates of interspecific interaction from phylogenetic trees. Ecol. Lett. 18 : 17 – 27 . Google Scholar CrossRef Search ADS PubMed Nuismer, S.L., Jordano P., Bascompte J. 2013 . Coevolution and the architecture of mutualistic networks. Evolution 67 : 338 – 354 . Google Scholar CrossRef Search ADS PubMed Oksanen, J., Blanchet F.G., Kindt R., Legendre P., Minchin P.R., O’Hara R.B., Simpson G.L., Solymos P., Stevens M.H.H., Wagner H. 2016 . vegan: Community Ecology Package. R Package Version 2.3-3. Available from: https://cran.r-project.org/web/packages/vegan/index.html. Olesen, J.M., Bascompte J., Dupont Y.L., Jordano P. 2007 . The modularity of pollination networks. Proc. Natl. Acad. Sci. USA 104 : 19891 – 19896 . Google Scholar CrossRef Search ADS Olesen, J.M., Bascompte J., Elberling H., Jordano P. 2008 . Temporal dynamics in a pollination network. Ecology 89 : 1573 – 1582 . Google Scholar CrossRef Search ADS PubMed Orme, D. 2013 . The caper package: comparative analysis of phylogenetics and evolution in R. R package version 5. Available from: https://cran.r-project.org/web/packages/caper/index.html. Paine, R.T. 1966 . Food web complexity and species diversity. Am. Nat. 100 : 65 – 75 . Google Scholar CrossRef Search ADS Pellmyr, O. 2003 . Yuccas, yucca moths, and coevolution: a review. Ann. Mo. Bot. Gard. 90 : 35 – 55 . Google Scholar CrossRef Search ADS Ponisio, L.C., Gaiarsa M.P., Kremen C. 2017 . Opportunistic attachment assembles plant–pollinator networks. Ecol. Lett. 20 : 1261 – 1272 . Google Scholar CrossRef Search ADS PubMed Quesnel, P.-O., Côté B., Fyles J.W., Munson A.D. 2006 . Optimum nutrient concentrations and CND scores of mature white spruce determined using a boundary-line approach and spatial variation of tree growth and nutrition. J. Plant Nutr. 29 : 1999 – 2018 . Google Scholar CrossRef Search ADS R Core Team. 2013 . R: A language and environment for statistical computing. R Core Team . Available from: http://www.R-project.org/. Revell, L.J. 2012 . phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3 : 217 – 223 . Google Scholar CrossRef Search ADS Rezende, E.L., Albert E.M., Fortuna M.A., Bascompte J. 2009 . Compartments in a marine food web associated with phylogeny, body mass, and habitat structure. Ecol. Lett. 12 : 779 – 788 . Google Scholar CrossRef Search ADS PubMed Rezende, E.L., Lavabre J.E., Guimaraes P.R., Jordano P., Jr., Bascompte J. 2007 . Non-random coextinctions in phylogenetically structured mutualistic networks. Nature 448 : 925 – U6 . Google Scholar CrossRef Search ADS PubMed Rowan, T.H. 1990 . Functional stability analysis of numerical algorithms. Ph.D. thesis. University of Texas (Austin) . Sazatornil, F.D., Moré M., Benitez-Vieyra S., Cocucci A.A., Kitching I.J., Schlumpberger B.O., Oliveira P.E., Sazima M., Amorim F.W. 2016 . Beyond neutral and forbidden links: morphological matches and the assembly of mutualistic hawkmoth-plant networks. J. Anim. Ecol. 85 : 1586 – 1594 . Google Scholar CrossRef Search ADS PubMed Scales, J.A., King A.A., Butler M.A. 2009 . Running for your life or running for your dinner: what drives fiber-type evolution in lizard locomotor muscles? Am. Nat. 173 : 543 – 553 . Google Scholar CrossRef Search ADS PubMed Schluter, D., Price T.D., Rowe L. 1991 . Conflicting selection pressures and life history trade-offs. Proc. R. Soc. Lond. B Biol. Sci. 246 : 11 – 17 . Google Scholar CrossRef Search ADS Slabbekoorn, H., Smith T.B. 2002 . Habitat-dependent song divergence in the Little Greenbul: An analysis of environmental selection pressures on acoustic signals. Evolution 56 : 1849 – 1858 . Google Scholar CrossRef Search ADS PubMed Stang, M., Klinkhamer P.G., Waser N.M., Stang I., van der Meijden. E. 2009 . Size-specific interaction patterns and size matching in a plant–pollinator interaction web. Ann. Bot. 103 : 1459 – 1469 . Google Scholar CrossRef Search ADS Stiles, F.G. 1975 . Ecology, flowering phenology, and hummingbird pollination of some Costa Rican Heliconia species. Ecology 56 : 285 – 301 . Google Scholar CrossRef Search ADS Thompson, J.N. 1994 . The coevolutionary process. Chicago, IL : University of Chicago Press . Thompson, J.N. 2005 . The geographic mosaic of coevolution. Chicago, IL : University of Chicago Press . Uyeda, J.C., Harmon L.J. 2014 . A novel Bayesian method for inferring and interpreting the dynamics of adaptive landscapes from phylogenetic comparative data. Syst. Biol. 63 : 902 – 918 . Google Scholar CrossRef Search ADS PubMed Ypma, J. 2014 . nloptr: R interface to NLopt. R package version 1 . Available from: https://cran.r-project.org/web/packages/index.html. © The Author(s) 2018. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. All rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

### Journal

Systematic BiologyOxford University Press

Published: Sep 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations