Genetic architecture of multiple mutualisms and mating system in Turnera ulmifolia

Abstract

Plants often associate with multiple arthropod mutualists. These partners provide important services to their hosts, but multiple interactions can constrain a plant's ability to respond to complex, multivariate selection. Here, we quantified patterns of genetic variance and covariance among rewards for pollination, biotic defence and seed dispersal mutualisms in multiple populations of Turnera ulmifolia to better understand how the genetic architecture of multiple mutualisms might influence their evolution. We phenotyped plants cultivated from 17 Jamaican populations for several mutualism and mating system‐related traits. We then fit genetic variance–covariance (G) matrices for the island metapopulation and the five largest individual populations. At the metapopulation level, we observed significant positive genetic correlations among stigma–anther separation, floral nectar production and extrafloral nectar production. These correlations have the potential to significantly constrain or facilitate the evolution of multiple mutualisms in T. ulmifolia and suggest that pollination, seed dispersal and defence mutualisms do not evolve independently. In particular, we found that positive genetic correlations between floral and extrafloral nectar production may help explain their stable coexistence in the face of physiological trade‐offs and negative interactions between pollinators and ant bodyguards. Locally, we found only small differences in G among our T. ulmifolia populations, suggesting that geographic variation in G may not shape the evolution of multiple mutualisms.

Abstract

We found genetic correlations between stigma–anther separation and investment in pollination and biotic defence mutualisms in Turnera ulmifolia. These correlations have the capacity to affect the evolution of multiple mutualisms by linking the evolutionary fates of traits associated with mating system, pollination and biotic defence.

biotic defence, G matrix, multiple mutualisms, pollination, seed dispersal

INTRODUCTION

Organisms typically interact with multiple mutualists and antagonists, either simultaneously or sequentially (Strauss, 2014; TerHorst et al., 2018). The outcomes of these multiple interactions are often difficult to predict (Strauss, 2014), given the potential for different community members to exert opposing fitness effects (Sletvold et al., 2015) or to interact nonadditively (Afkhami et al., 2014; Afkhami & Stinchcombe, 2016). Incorporating the multiple interactions that capture the complexity of real communities has improved our understanding of the ecology of species interactions (Strauss, 2014; TerHorst et al., 2018), but relatively few studies have explored how populations adapt to the multiple species with which they interact.

The translation of ecological complexity into the evolution of species interactions is complicated by the fact that organisms do not respond to selection imposed by individual partners in isolation. Rather, selection is a multivariate process acting on many traits, and its effects depend on genetic features such as pleiotropy and linkage disequilibrium among loci targeted by selection (Kopp & Matuszewski, 2014; Lande, 1979; TerHorst et al., 2018). These genetic features may result in genetic covariances among traits that affect interactions with different community members and can either facilitate or constrain phenotypic evolution, depending on whether the principal axes of covariation are aligned with those of selection (Agrawal & Stinchcombe, 2009; Arnold, 1992; Assis et al., 2020; Teplitsky et al., 2014; TerHorst et al., 2018; Wise & Rausher, 2013). More formally, Lande's equation, R = Gβ, predicts that the vector of responses to selection on multiple traits, R, equals the product of the genetic variance–covariance matrix, G, and the vector of selection gradients acting on those traits, β (Hansen & Houle, 2008; Lande & Arnold, 1983). Traits that mediate species interactions will therefore have a correlated response to selection whenever they are genetically correlated, such that multiple species interactions may evolve in tandem. The extent to which genetic correlations affect the evolution of multiple species interactions, however, remains poorly understood.

While genetic covariances may in general have only limited effects on rates of adaptation (Agrawal & Stinchcombe, 2009), the evolution of traits involved in multiple species interactions may be especially susceptible to this form of constraint. Multiple species interactions are often mediated by several traits in a focal species in ways that can increase the prevalence of correlational selection (TerHorst et al., 2018; Wise & Rausher, 2013). Correlational selection occurs when fitness depends on interactive effects of two or more traits and is known to influence G matrix evolution (Svensson et al., 2021). When multiple partner species interact and affect the fitness landscape of multiple interactions (e.g. multiple mutualist effects, Afkhami et al., 2014), they fundamentally change the nature of multivariate selection, creating new patterns of correlational selection dependent on the presence of specific partner species (Afkhami et al., 2021; Friman & Buckling, 2013; Guimarães et al., 2017; Keith & Mitchell‐Olds, 2019; Wise, 2010; Wise & Rausher, 2013, 2016). Moreover, phenotypes that affect more than one ecological interaction (either directly or through diffuse interactions among partner species) also create the possibility of antagonistic or synergistic effects of multiple species interactions on the responses of even single traits to selection (Afkhami et al., 2021; Cuautle et al., 2005; Sletvold et al., 2015; Villamil et al., 2019; Wise & Rausher, 2016). Patterns of correlational selection created by these context‐dependent ecological interactions can lead to the evolution of strong genetic covariance among traits (Guimarães et al., 2017; Jones et al., 2004, 2012; Penna et al., 2017; Svensson et al., 2021), which in turn may constrain the direction and magnitude of future evolution (Hangartner et al., 2020; Keith & Mitchell‐Olds, 2019; Silva et al., 2020).

To date, precious few studies have investigated how the interplay between ecological interactions and genetic architecture affects the evolution of multiple mutualisms (Guimarães et al., 2017). Multiple mutualisms are ubiquitous; organisms interact with several partners, often from different guilds (Dáttilo et al., 2016; Guimarães et al., 2017; Stanton, 2003). Plants associate with many beneficial organisms, including animal pollinators (Ollerton et al., 2011), seed dispersers, and biotic defenders such as ant bodyguards (Marazzi et al., 2013; Weber & Keeler, 2013) in addition to mycorrhizal fungi and diverse microbiomes (which are themselves a type of multiple mutualism). Here, we focus on the aboveground mutualistic interactions of Turnera ulmifolia (Passifloraceae), a weedy neotropical plant that associates with a diversity of arthropod mutualists, including pollinators and both defensive and seed‐dispersing ants (Barrett & Shore, 1987; Dutton, Luo, et al., 2016; Dutton, Shore, et al., 2016). In T. ulmifolia, these interactions are mediated by multiple plant traits, principally floral nectar that is consumed by insect pollinators, and elaiosomes (fleshy, lipid‐rich appendages attached to seeds) and extrafloral nectar that are consumed by ants that disperse seeds and defend the plant against herbivores (Barrett & Shore, 1987; Cuautle et al., 2005; Dutton, Luo, et al., 2016; Dutton, Shore, et al., 2016; Villamil et al., 2019).

Conflict between ants and pollinators, which consume similar nectar rewards often in close proximity on plants (Dutton, Luo, et al., 2016; Keith & Mitchell‐Olds, 2019; Martíınez‐Bauer et al., 2015; Villamil et al., 2019), can give rise to nonadditive effects of pollinators and ants on plants (e.g. multiple mutualist effects, Afkhami et al., 2014). Direct interactions between ants and pollinators are common, as aggressive ant mutualists can deter pollinators from visiting in addition to reducing herbivory (Cembrowski et al., 2014; Ness, 2006; Villamil et al., 2018). Biotic defence and pollination are further linked through herbivores, which can be attracted to showy floral displays. This can create an ecological trade‐off between attracting mutualists and repelling antagonists (Knauer & Schiestl, 2017), such that the evolution of floral traits affecting pollination may be driven not only by pollinators, but also by herbivores and the effectiveness of plant defences (Agren et al., 2013; Cabin et al., 2022; Gervasi & Schiestl, 2017; Kessler et al., 2015; Ramos & Schiestl, 2019; Sletvold et al., 2015; Soper Gorden & Adler, 2016; Strauss et al., 1996; Thompson & Johnson, 2016). We might therefore expect to see correlational selection favouring some trait combinations (e.g. greater floral rewards together with high investment in antiherbivore defence), potentially leading to genetic covariances among plant traits. Furthermore, the existence of genetic covariances, such as that suggested by the genetic and physiological links between the production of biochemically similar rewards like extrafloral and floral nectar (Dutton, Luo, et al., 2016; Heil, 2011), may feed back to affect the evolution of mutualistic plant traits.

In addition to potentially linking the fates of pollination and biotic defence through patterns of correlational selection and genetic covariance, multiple mutualisms may shape the coordinated evolution of phenotypes associated with dispersal, defence and reproduction. The evolution of selfing, often described as a syndrome that involves a set of linked phenotypes including reductions in floral size, nectar volume and stigma–anther separation (Belaoussoff & Shore, 1995; Sicard & Lenhard, 2011), can be expected to covary with mutualistic seed dispersal and antiherbivore traits. Selfing is often a mechanism of reproductive assurance in good colonizers (Baker, 1955; Pannell & Barrett, 1998), and if natural selection favours selfing in highly dispersive genotypes, these traits may covary genetically. Similarly, in the small, isolated populations established by good dispersers at range edges or in ephemeral habitats, plants have the potential to escape their predators and parasites and may experience reduced selection to maintain specialized, constitutive defences, sometimes in favour of phenotypically plastic, inducible defences (Brudvig et al., 2015; Campbell & Kessler, 2013; Chevin & Lande, 2011). Over time, this may lead genetic covariances between dispersal and defence traits to build up. In our study species, we have further reason to think that dispersal and antiherbivore defence traits may be linked; T. ulmifolia's seeds are dispersed by ants that are attracted not only to elaiosomes, but to extrafloral nectar, which also has a defensive function (Cuautle et al., 2005; Dutton, Shore, et al., 2016).

Ecological and evolutionary links among reproduction, defence and dispersal in plants raise the possibility that traits governing the multiple mutualisms mediating these processes do not evolve independently. Indeed, similarity in the composition of plant rewards like floral and extrafloral nectar suggests a shared genetic architecture (Afkhami & Stinchcombe, 2016; Dutton, Luo, et al., 2016), and previous studies of the evolution of biotic defence and pollination at broad phylogenetic scales have hinted at correlated trait evolution (Chamberlain & Rudgers, 2012). Whether plant traits associated with multiple mutualisms are genetically correlated to one another is however not adequately understood, and traits that intuitively seem linked may not always be (Ossler & Heath, 2018). Furthermore, while phenotypic covariances between ecologically important traits are common, these can be environmentally induced, rather than genetic, and parsing genetic from environmental covariances requires an appropriate study design (Stinchcombe et al., 2002). More generally, evolutionary ecology requires a better understanding of how the genetic architecture of these traits interacts with their complex ecology to determine the outcome and evolutionary trajectories of multiple mutualisms (Assis et al., 2020; Guimarães et al., 2017; Lau & Terhorst, 2015; TerHorst et al., 2018).

Here, we estimate genetic variance–covariance (G) matrices for phenotypic traits associated with multiple mutualisms in T. ulmifolia. We grew plants in a common greenhouse environment and measured plant investment in rewards to pollinators, bodyguards and seed dispersers, as well as variation in mating system, and fit G matrices for the 5 largest populations we sampled as well as across the island of Jamaica. We fit both local and island‐wide G matrices to assess trait covariances at different spatial scales; T. ulmifolia is a good colonizer that readily spreads to new regions, suggesting few barriers to gene flow at small spatial scales. Because we grew plants from seeds of known parentage in a common environment, we could estimate the total genetic variance for each trait as well as genetic covariances between each pair of traits without the confounding effects of environmental variation. Specifically, we addressed the following questions: (1) Does the production of floral nectar (FN), extrafloral nectar (EFN) and elaiosomes covary genetically in T. ulmifolia, and do these mutualistic traits covary with stigma–anther separation? (2) How variable are the G matrices of different T. ulmifolia populations across a small geographic range? and (3) How might G facilitate or constrain the evolution of multiple mutualisms under simulated, but realistic selection scenarios?

MATERIALS AND METHODS

Study system

Turnera ulmifolia is a continuously flowering, ruderal, herbaceous perennial that readily colonizes roadsides and disturbed habitats in its native range of the Caribbean (Barrett, 1978; Barrett & Shore, 1987). Formerly identified as T. ulmifolia var. angustifolia within the Turnera ulmifolia L. complex characterized by extensive chromosomal variation (Shore & Barrett, 1985; Urban, 1883), all known populations of T. ulmifolia are hexaploid (López et al., 2013). Turnera ulmifolia is self‐compatible and its populations are generally small, ephemeral and moderately inbred (Barrett & Shore, 1987; Belaoussoff & Shore, 1995; Dutton, Luo, et al., 2016). Though self‐compatible, its flowers are conspicuous (Figure 1) and are visited by pollinating insects and hummingbirds that increase seed set (Barrett, 1978; Barrett & Shore, 1987; Cuautle & Rico‐Gray, 2003). These flowers are ephemeral, opening in the morning around 7 AM and then closing and wilting in the early afternoon (Barrett, 1978). Turnera ulmifolia's out‐crossing rate is correlated with stigma–anther separation, a variable trait in Jamaica (our study location) that suggests directional selection for increased out‐crossing in some populations (Barrett & Shore, 1987; Belaoussoff & Shore, 1995; Shore & Barrett, 1990).

FIGURE 1

Photographs of Turnera ulmifolia. Top row, left: elaiosomes (indicated by white arrow) are visible as fleshy, beige structures affixed to brown seeds within a dehisced capsule. Top row, centre: extrafloral nectar (EFN) accumulating at nectaries on the petioles of leaves. Top row, right: flowers, showing variation in stigma–anther separation. Bottom: plants growing in the greenhouse. Photographs by C. Reid.

In addition to producing large, nectar‐rich flowers, T. ulmifolia produces extrafloral nectar (EFN) from glands located on the petioles of leaves (Dutton, Luo, et al., 2016, Figure 1). This EFN attracts ant and wasp bodyguards that defend plants against herbivores and increase reproductive success (Cuautle & Rico‐Gray, 2003; Dutton, Luo, et al., 2016). Turnera ulmifolia's seeds are dispersed by mutualistic ants which collect seeds from dehisced capsules and feed the attached lipid‐rich appendages (elaiosomes, Figure 1) to their larvae. The recruitment of seed‐dispersing ants to T. ulmifolia's seed capsules is also affected by the production of EFN, which they collect and consume (Dutton, Shore, et al., 2016; Salazar‐Rojas et al., 2012). Prior to this study, little was known about variation in mutualistic traits other than stigma–anther separation in T. ulmifolia, and very little is currently known about the extent of gene flow among Jamaican T. ulmifolia populations or the relevant scale at which selection acts on this species.

Plant collection and cultivation

In 2016, we collected seeds from 218 wild plants distributed across 17 populations of T. ulmifolia in Jamaica (Figure 2). At each site, we collected seeds from up to 20 maternal plants, though only six of the populations (Brown's Town, Cave, Haining, Mosquito Cove, Murdock and Saint Ann's Bay) comprised 20 or more mature plants (Table S1). We germinated seeds in a 1:1:1:1 mix of topsoil, peat loam, manure and concrete sand and transplanted seedlings into 3.5″ pots in the greenhouse at the University of Toronto (15L:9D cycle). We watered plants daily and fertilized them bi‐weekly. The greenhouse temperatures were set to 20–24°C and 18–19°C during the day and night, respectively, although the temperatures inside varied somewhat with season.

FIGURE 2

Map of Turnera ulmifolia collection sites in Jamaica. Population size is represented by circle diameter, and the five well‐sampled populations where we collected seeds from more than 20 maternal plants are labelled.

The majority of plants reached reproductive maturity after approximately 8 months. To reduce maternal environmental effects, we generated selfed seed by hand‐pollinating flowers from each plant with pollen from the same flower and collected mature seeds after approximately 20–25 days (Dutton, Luo, et al., 2016). In March 2017, we germinated seeds from 197 maternal families (at most one per plant from which we collected seeds in the field) and planted 10 seedlings per maternal line, for a total of 1970 plants. Each seedling was planted in a 5″ pot. We randomly assigned plants to two greenhouses and continued to care for them as described above.

Trait measurements

Once plants began flowering, we measured stigma–anther separation using callipers. To quantify plant investment in pollination, defence and seed dispersal mutualisms, we measured floral nectar (FN) production, EFN production and elaiosome mass, respectively (Figure 1). As the ephemeral flowers of T. ulmifolia are a highly dynamic and variable resource for pollinators, we quantified floral resources at the pollinator‐relevant level of individual flowers. We used microcapillaries to collect FN from a single open flower per plant. The number of extrafloral nectaries is substantially more stable over time, with new nectaries being produced only with the growth of new leaves. Furthermore, given that ants forage and collect EFN from the entirety of individual plants in the field, we measured the abundance of nectar rewards available for biotic defenders at the plant level. Twenty‐four hours prior to collecting EFN, we thoroughly washed plants with water to dislodge any nectar that had previously accumulated or crystallized on extrafloral nectaries. The next day, we collected all the EFN that had accumulated on the plants' active nectaries in 24 h using microcapillaries; we later divided them by the number of nectaries to get a measure of EFN production per nectary. To determine the total sugar content of FN and EFN in each sample, we diluted samples with known volumes of distilled water and measured the total volume with callipers. We measured the sugar concentration of FN and EFN in degrees Brix using a hand‐held refractometer. We then calculated the total sugar mass of the original samples (see Appendix S2 for full details). Finally, we collected autogamous seed from mature seed capsules and dried them for several months at room temperature. We separated intact elaiosomes from three seeds per plant and weighed them on a microbalance. Some plants died or failed to mature; in total, we were able to collect phenotypic data from 1892 plants across 197 maternal families.

Data analysis

We performed all statistical analyses in R version 3.5.1 (R Core Team, 2022). Briefly, we fit a common genetic variance–covariance matrix (G_W, Chenoweth et al., 2010) for the Jamaican metapopulation of T. ulmifolia using Bayesian methods. Because the extent of gene flow among populations and the scale at which natural selection acts on T. ulmifolia is unknown, we also estimated G for our five largest individual populations and evaluated the extent to which they differed from one another using three different approaches. We also fit a D matrix summarizing multivariate divergence in trait means among our multiple populations and assessed whether it is aligned with G_W to determine whether G may have constrained the direction of phenotypic divergence in T. ulmifolia. Finally, we used our local estimates of G to assess evolutionary constraint and facilitation in the multiple mutualisms of T. ulmifolia and calculated heritabilities for elaiosome mass, EFN production, FN production and stigma–anther separation.

Estimation and significance testing of G matrices

We estimated the common genetic variance–covariance (G_W) matrix for the island metapopulation using all plants for which we had trait data (1892 plants from 197 maternal inbred families). G matrices for this metapopulation, and those later fit for individual populations, were estimated through the implementation of multivariate ‘animal models’ that included elaiosome mass, EFN production, FN production and stigma–anther separation as response variables. In these models, a vector of individual phenotypes (Y) is fitted as a function of the mean phenotypes of a population (μ, fixed effects), the genetic breeding values of individuals (g, random effects) and residual error (e). A fuller description of the multivariate animal model we implemented is:

$$Y=\mu +Z_{g}g+Z_{b}b+e$$

where Y is the vector of mean‐scaled phenotypes for all plants and μ is the vector of mean phenotypes. The genetic effects are represented by the vector g and the associated design matrix Z_g, while the random effects of block are represented by the vector of block effects b and its corresponding design matrix, Z_b (Kruuk, 2004; Teplitsky et al., 2014). Variance and covariances among the estimated breeding values of individuals within a population (g) then yield the genetic variance–covariance matrix.

We fit our models with the ‘MCMCglmm’ package in R (Hadfield, 2010; R Core Team, 2022). For all individuals, we scaled phenotypic data against the grand mean of the Jamaican metapopulation. We scaled trait data against the island‐wide mean (rather than within each population separately) in order to estimate G_W as if the whole island of Jamaica were a single T. ulmifolia population not subject to limited dispersal (which it may be). To explore the potential influence of population‐level differences on patterns of covariation, we also estimated the matrix of divergence among population trait means (D). We fit D in a second model where we also estimated the average G matrix within populations—that is, an estimate of the metapopulation G that accounts for population structure. This model included covariance structures for both the random effects of site and genetic lineage (D and G matrices, respectively). For our estimation of G_W, we applied a set of weakly informative inverse Wishart priors that differed in their specification of predicted genetic variance (V = 0.25, 0.333, and 0.5) and degree of belief (nu = n, n‐1, and n‐2, where n is the number of traits). All priors produced similar G matrices; we present results for G matrices estimated with the prior selected by comparison of the Deviance Information Criterion (DIC) scores (V = 0.5, nu = 2, see Table S2). We iterated models 10 100 000 times and recorded values every 1000 generations with a burn‐in period of 100 000 generations, for a total of 10 000 posterior estimates. We confirmed model convergence by examining trace and posterior distribution plots and calculated autocorrelations between samples (all of which were below the recommended level of 0.1) using the R package ‘coda’ (Plummer et al., 2006).

We calculated the mean genetic variances and genetic correlations from the posterior distribution of the 10 000 MCMC samples. Estimates of genetic correlations among traits can range from −1 to 1, and we regarded genetic correlations as significant if their 95% highest posterior density intervals (HDPI) did not overlap 0. Estimates of genetic variance cannot be negative, so we determined the significance of these terms by comparing our estimates to those generated from null distributions where G matrices were fit to randomized data. First, we generated 1000 randomized data sets by randomly assigning phenotypes to individuals. We then fit G matrices to each of these randomized data sets and calculated the 95% confidence intervals of mean genetic variance for each of our four traits from these 1000 randomized G matrices. Estimates of genetic variance were considered significant if they exceeded the upper confidence limit of the null distribution. In addition to estimating genetic correlations among mutualistic traits, we calculated phenotypic correlation matrices as Pearson correlations among paired traits using the ‘psych’ package (Revelle, 2020), correcting for multiple tests in the determination of significance.

We calculated broad‐sense trait heritabilities in the Jamaican metapopulation of T. ulmifolia by fitting linear mixed‐effect models using restricted maximum likelihood (REML) in the ‘lme4’ package (Bates et al., 2015), and by implementing univariate Bayesian methods in MCMCglmm (Hadfield, 2010), using univariate versions of the animal models described above with modified priors (nu = 1). For the linear mixed models, we log‐transformed total FN sugar mass but did not re‐scale any trait data. We modelled maternal line and block as random effects, and we estimated broad‐sense heritability as the proportion of total variance explained by maternal family. We used likelihood ratio tests implemented with the ‘RLRsim’ package (Scheipl et al., 2008) to test the statistical significance of our heritability estimates by evaluating the significance of the random effect of maternal line.

Finally, to investigate the potential consequences of having selfed plants in the greenhouse for one generation on the expression of mutualistic traits, we measured out‐crossed, selfed and autogamous seed set in a subset of plants used to fit our G matrices. We fit a linear mixed‐effects model regressing seed set against stigma–anther separation that included pollination treatment as a fixed effect and population and maternal line as random effects to investigate the possibility that our breeding design had a disproportionate effect on genotypes from out‐crossing populations (see Appendix S2 for detailed methods). We followed up on this analysis by fitting a model in MCMCglmm to estimate G_w for the subset of plants used in the seed set experiments that included selfed seed set as a fifth phenotype. As with other G matrices, we mean‐scaled all phenotypes.

Comparison of G matrices

To investigate differences in the patterns of genetic variance and covariance among individual populations of T. ulmifolia, we estimated G matrices for our five largest field populations (Brown's Town, Cave, Mosquito Cove, Murdock and Saint Ann's Bay). We continued to scale phenotype data to the grand mean rather than within each population separately, because this method maintains differences in trait means across populations (Figure 3) that may be biologically meaningful (Hine et al., 2009; Puentes et al., 2016). We assessed potential differences in the overall size, or amount of total genetic variation among G matrices, by calculating the trace (the sum of diagonal elements) for each matrix and compared their 95% HPD intervals. We further explored variation in the orientation of G by calculating the eigenvectors and eigenvalues of our G matrices. Eigenvectors are linear vectors that capture the axes of greatest variation in multivariate space, and their corresponding eigenvalues reflect the relative amount of total variation they capture. By examining the eigenstructure of G, one can determine the combination of phenotypic traits that contributes most to multivariate variation in a dataset. To quantitatively test for differences among our G matrices, we used three methods: random skewers, Krzanowski's common subspace analysis and the fourth‐order genetic covariance tensor (Aguirre et al., 2014; Chakrabarty & Schielzeth, 2020; Hine et al., 2009). Here, we discuss only the specifics of the fourth‐order genetic covariance tensor, but detailed descriptions of the random skewers and Krzanowski's common subspace analyses are in Appendix S2.

FIGURE 3

Variation in trait means (with standard error) across 17 Turnera ulmifolia populations.

The genetic covariance tensor method of Hine et al. (2009) allows for the comparison of G matrices directly by decomposing axes of variation among multiple matrices through eigenanalysis. Briefly, a tensor of the 4th order is used to compare G matrices, which are themselves 2nd‐order objects. This tensor is represented by Σ and describes the variance and covariances among multiple covariance matrices. From this 4th‐order covariance tensor, the second‐order eigentensors, E that describe variation among matrices summarized by Σ, can be explored and used to obtain the eigenvectors (e) and eigenvalues of E. The eigenvalues of an eigentensor can then be examined to explore which trait combinations contribute to detected changes in the covariance structures of individual G matrices. We employed the genetic covariance tensor method to test for and examine patterns of population divergence in G using code modified from Aguirre et al. (2014). We used Bayesian posterior estimates of G to project the eigentensors (E) onto Σ to determine 95% HPD intervals for α, the variance of Σ accounted for by each eigentensor. We then tested for significant differences among our underlying G matrices by comparing these estimates of α with those estimated from randomized G matrices. We generated randomized null G matrices for the purpose of significance testing by sampling and assigning breeding values (Aguirre et al., 2014) before reconstructing randomized individual phenotypes by incorporating empirical estimates of environmental variance. We then fit null G matrices using these phenotypes (Morrissey et al., 2019) and used the resultant randomized estimates of G to calculate null summary statistics for the genetic covariance tensor analysis and other G matrix comparison methods. Finally, we projected the leading eigenvectors (e) of the eigentensors (E) onto our observed G matrices in order to estimate the amount of genetic variance in the direction of the principal axes of multivariate genetic variation in each population, using all posterior estimates of G to fit 95% HPD intervals.

Evolutionary constraint and facilitation through G

To investigate the effects of genetic covariances on responses to selection in the island metapopulation and five largest populations of T. ulmifolia, we multiplied our estimates of G by three vectors of directional selection gradients representing simplified but plausible selection regimes. Measurements of selection on mutualistic traits in the field are rare, but our skewers reflect conservative estimates of selection based on published accounts. There are no estimates of selection on elaiosome mass, but their presence increases seed removal by ants by up to 20%, a large increase in fitness given that dispersal to ant nests can double germination rates in T. ulmifolia (Cuautle et al., 2005; Salazar‐Rojas et al., 2012). Weak positive selection on EFN is also likely common, given that its presence boosts seed dispersal by a factor of up to 2.5 (Cuautle et al., 2005; Dutton, Shore, et al., 2016), and selection through plant defence can range from weakly negative (β = −0.096) to strongly positive (β = 0.695) depending on the presence of ants and herbivores (Rutter & Rausher, 2004). Evidence of selection on floral traits is more mixed, as papers report little or no selection on FN volume and sugar content in flowering plants (Benitez‐Vieyra et al., 2010; Gijbels et al., 2014; Kulbaba & Worley, 2011), and positive selection on stigma–anther separation can be sex‐specific (Kulbaba & Worley, 2011). In all cases, we specified conservative skewers that featured either weak positive selection (β = 0.1), neutral selection (β = 0) or weak negative selection (β = −0.1) consistent with the range of published values.

In the first scenario, we specified universal, weak positive selection for all traits, while in the second selection regime, we simulated weak selection for defence (EFN) and dispersal (elaiosome mass), and against pollinator attraction (FN) and out‐crossing (stigma–anther separation). In the third scenario, we simulated positive selection for defence and dispersal and against stigma–anther separation, but made FN production selectively neutral. We chose these three selection regimes because they represented ecologically plausible patterns of selection; namely, universal selection for investment in multiple mutualisms (1), selection for traits associated with dispersal, including concomitant selection for inducible defences (EFN) and against floral traits associated with pollination and out‐crossing (2), and selection for traits associated with selfing, with neutral selection on FN production given some evidence that it may be a relatively inexpensive resource (Harder & Barrett, 1992), and may not affect seed production in T. ulmifolia (Dutton, Luo, et al., 2016) (3).

We projected selection vectors through all 10 000 posterior estimates to produce response vectors and calculated the 95% HPD intervals around the selection responses of individual traits for the purpose of significance testing. Additionally, we followed Agrawal and Stinchcombe (2009) to quantify the effect of genetic covariances on the rate of adaptation in our five best sampled populations. Briefly, we compared the rate of adaptation with the observed G matrices to the rate of adaptation with the genetic covariances (off‐diagonals) set to zero. In addition to calculating the responses to the selection of our observed G matrices under our simulated selection scenarios, we calculated response vectors for their counterparts without genetic covariances. We then calculated R values by dividing the response vectors obtained with the full G matrix by those obtained with genetic covariances set to 0. R values > 1 indicate that genetic correlations facilitate adaptation, while values <1 indicate that adaptation is constrained by genetic correlations among traits (Agrawal & Stinchcombe, 2009).

In addition to assessing the impact of genetic correlations on the evolvability of mutualistic traits by measuring their effects on the response to selection, we also estimated the matrix of divergence among population trait means (D, Blows & Higgie, 2003). We estimated D as the covariance structure explained by the random effect of population in our model and compared it to the estimate of G fit in the same model for the island‐wide metapopulation. This estimate of G represents the average G matrix within populations, accounting for differences in trait means and sample size across populations. The eigenvectors or principal components of D represent the major multivariate axes of phenotypic variation across populations and can be compared with the eigenvectors of G to determine whether or not the multivariate axes of genetic variation (e.g. genetic lines of least resistance, Schluter, 1996) may have constrained phenotypic divergence between populations. To compare G and D, we employed both one and two‐dimensional comparison methods. We began by calculating the angle between the leading eigenvectors of each matrix (g_max and d_max), to investigate how closely aligned the vectors describing maximal genetic and phenotypic differentiation are (Blows & Higgie, 2003; Schluter, 1996). Additionally, we explored the extent to which G has constrained phenotypic divergence among populations by comparing variation in the size and shape of G and D. For both summary matrices, we calculated the size of the matrices by summing their eigenvectors and assessed their eccentricity by calculating the amount of variation captured by their leading eigenvectors. For each of these measures, we calculated the 95% HPD intervals from the posterior estimates of our MCMCglmm objects. Finally, we employed the Krzanowski subspace analysis (described in Appendix S2) to determine whether G and D share a common multivariate substructure (Silva et al., 2020).

RESULTS

Genetic variances and covariances among traits

We found significant genetic variation for all four traits in the Jamaican metapopulation of T. ulmifolia; all genetic variance estimates in the metapopulation G matrix G_W (Table 1), and all trait heritabilities calculated from the full data set using both REML and univariate Bayesian approaches (Table S3) were statistically significant. The genetic variance in most of these traits was considerable, ranging from 0.036 to 0.485 depending on the trait (Table 1). Both Bayesian and REML approaches found that EFN production was the least heritable trait (0.17–0.22), albeit still with significant genetic variation among maternal lines, and stigma–anther separation was the most heritable trait (0.61–0.77), with FN production and elaiosome mass intermediate (0.48–0.61 and 0.46–0.55, respectively) (Table 1 and Tables S3 and S4).

TABLE 1

Genetic variance–covariance (G) matrices of the Jamaican metapopulation and five largest individual populations of Turnera ulmifolia

	Elaiosome m	EFN	FN	S‐A separation
(a) Jamaica
Elaiosome m	0.036 (0.030, 0.043)
EFN	0.007 (−0.001, 0.016)	0.080 (0.059, 0.103)
FN	0.005 (−0.010, 0.020)	0.052 (0.022, 0.083)	0.334 (0.261, 0.410)
S‐A sep.	−0.008 (−0.024, 0.008)	0.053 (0.021, 0.087)	0.085 (0.031, 0.139)	0.485 (0.411, 0.565)
(b) Brown's Town
Elaiosome m	0.036 (0.025, 0.048)
EFN	0.002 (−0.015, 0.021)	0.112 (0.053, 0.182)
FN	−0.001 (−0.071, 0.064)	0.04 (−0.154, 0.247)	1.511 (0.707, 2.391)
S‐A sep.	0.003 (−0.021, 0.028)	0.01 (−0.057, 0.081)	0.30 (0.069, 0.555)	0.215 (0.106, 0.341)
(c) Cave
Elaiosome m	0.049 (0.034, 0.065)
EFN	0.014 (−0.008, 0.038)	0.131 (0.067, 0.211)
FN	−0.012 (−0.038, 0.014)	−0.01 (−0.061, 0.042)	0.212 (0.14, 0.289)
S‐A sep.	0.019 (−0.024, 0.059)	0.0 (−0.056, 0.122)	−0.11 (−0.186, −0.027)	0.418 (0.249, 0.613)
(d) Mosquito Cove
Elaiosome m	0.056 (0.036, 0.077)
EFN	−0.005 (−0.042, 0.028)	0.199 (0.086, 0.344)
FN	−0.035 (−0.12, 0.045)	0.096 (−0.162, 0.38)	1.046 (0.22, 2.023)
S‐A sep.	−0.03 (−0.095, 0.027)	0.046 (−0.122, 0.234)	0.399 (−0.047, 0.891)	0.662 (0.336, 1.018)
(e) Murdock
Elaiosome m	0.052 (0.036, 0.07)
EFN	0.009 (−0.018, 0.038)	0.177 (0.072, 0.303)
FN	−0.009 (−0.038, 0.019)	0.013 (−0.074, 0.111)	0.215 (0.103, 0.339)
S‐A sep.	0.007 (−0.018, 0.032)	0.025 (−0.051, 0.103)	0.029 (−0.041, 0.102)	0.189 (0.108, 0.271)
(f) Saint Ann's Bay
Elaiosomem	0.06 (0.041, 0.083)
EFN	0.001 (−0.019, 0.02)	0.088 (0.05, 0.135)
FN	−0.021 (−0.041, −0.003)	0.006 (−0.028, 0.04)	0.115 (0.07, 0.163)
S‐A sep.	−0.054 (−0.088,  −0.021)	0.017 (−0.041, 0.071)	0.079 (0.019, 0.137)	0.31 (0.198, 0.437)

Note: Reported values are the mean and 95% highest posterior density intervals of estimates. Bolded values are significantly different from 0.

For the Jamaican metapopulation, all significant genetic correlations among traits were positive. We found that EFN and FN production was significantly genetically correlated, and both covaried significantly with stigma–anther separation (Table 1). In other words, maternal lines that secrete more EFN also make more herkogamous flowers that secrete more FN. In contrast, the genetic correlations between elaiosome mass and stigma–anther separation, EFN, or FN production were small in magnitude (−0.008 to 0.007) and not significantly different from zero.

We found no evidence that selfing plants differentially affected genotypes from populations with higher stigma–anther separation. Although stigma–anther separation and seed set were negatively correlated across all treatments (Tables S5 and S6), we did not observe a significant interaction between pollen source and stigma–anther separation (Figure S1; Table S5). This analysis revealed a (positive) genetic correlation between selfed seed set and elaiosome mass, but no significant genetic correlations between seed set and the production of either EFN or FN, albeit in analysis with low power for estimating G (Table S6).

Comparison of G matrices

Our fourth‐order genetic covariance tensor analysis found no significant divergence among our five largest populations of T. ulmifolia. For each of the four eigentensors of the genetic covariance tensor, the 95% HPD intervals of α (the amount of genetic variance in the direction of each eigentensor) of the observed and randomized G matrices overlapped (Figure 4a). Thus, these eigentensors do not describe significant variation among the observed G matrices of our five largest populations. These results are consistent with the results of both the application of random skewers and Krzanowski's common subspace analysis, which both failed to detect significant differences among the G matrices of our largest populations (see Appendix S1, Figures S2 and S3). Taken together, these results suggest little divergence among our populations in the size and orientation of G.

FIGURE 4

Results of the 4th‐order genetic covariance tensor analysis comparing the G matrices of our five largest Turnera ulmifolia populations. (a) The amount of variance (α) in G explained by each eigentensor E of the 4th‐order covariance tensor Σ for observed (blue) and randomized (red) G matrices. Genetic variance in the direction of (b) the first eigenvector of E1 and (c) the first eigenvector of E2. All error bars are 95% HPD intervals.

These comparisons, however, are frequently limited by issues of statistical power (Puentes et al., 2016; Teplitsky et al., 2014), particularly when G matrix estimates are fit with wide confidence intervals. Greater sampling may have shown differences in G among populations, given that projecting the leading eigenvectors of the first and second eigentensors of our covariance tensor (e₁₁ and e₂₁) back onto population‐specific estimates of G revealed greater genetic variance in the Mosquito Cove and Brown's Town populations in the directions of both e₁₁ and e₂₁ (Figure 4b,c) (see also Appendix S2 for more details). Likewise, our examination of the total amount of genetic variance in each of our population estimates of G hinted at differences among populations. G matrices fit for Brown's Town (mean and 95% HPD intervals of matrix traces: 1.87, 0.89–2.96) and Mosquito Cove (1.96, 0.68–3.46) trended towards greater total variance than other populations, but wide confidence intervals rendered differences nonsignificant with the sole exception of Brown's Town having greater genetic variance than Saint Ann's Bay (0.57, 0.36–0.82, see Table 1). However, greater sampling in most T. ulmifolia populations was not possible, because of very small population sizes (Table S1).

Our estimates of the phenotypic correlation (P) matrices for the individual populations and metapopulation of T. ulmifolia were broadly aligned with G. In 32 of 36 paired trait combinations, the confidence intervals of our estimates of phenotypic and genetic correlations overlapped with one another (Table 1; Table S9). Three of the differences in G and P we observed were correlations involving elaiosome mass, which displayed weaker genetic than phenotypic correlations across the board. Because we grew plants and measured phenotypes in a common greenhouse environment, we expected that P and G matrices would be similar, and they were; we observed only one case where the sign of trait correlations was significant and different from one another.

Evolutionary constraint and facilitation through G

Applying selection gradients onto the G matrices of the Jamaican metapopulation and our five largest populations revealed that the patterns of genetic correlation we observed can affect the adaptive potential of T. ulmifolia. When we simulated weak positive selection on all traits (scenario 1) in the Jamaican metapopulation, we found that positive correlations among extrafloral and floral nectar production and stigma–anther separation facilitated their response to simulated selection (mean R value with 95% HPD intervals: EFN—2.41 (1.84–2.97), FN—1.43 (1.22–1.64), stigma–anther separation—1.27 (1.12–1.42), Figure 5a; Figure S4a). Unsurprisingly, these same correlations have the potential to dramatically reduce the visibility of mutualistic traits to multivariate selection when selection acts in opposing directions on floral traits and extrafloral nectar production. When we simulated a selection regime associated with selection for dispersal, inducible biotic defences and increased selfing (2), we found that the evolution of extrafloral nectar was significantly constrained (mean R value = −0.230, upper 95% HPD = 0.348, Figure S4b) to the point that it showed no response to positive selection (Figure 5b).

FIGURE 5

Responses of mutualistic traits to multivariate selection gradients in the Jamaican metapopulation and the five largest Turnera ulmifolia populations. We simulated selection scenarios in which (a) all traits are under positive directional selection (β = 0.1), (b) elaiosome mass (Ela m) and EFN are under positive directional selection while FN production and stigma–anther separation (S‐A sep) are selected against (β = −0.1), and (c) elaiosome mass and EFN are under positive directional selection, S‐A sep is selected against, and FN is selectively neutral (β = 0).

Among the five largest populations of T. ulmifolia, we found that slight but nonsignificant differences in G had the capacity to alter their responses to multivariate selection. In some cases, our analyses showed that population‐specific estimates of G nearly completely inhibited their capacity to respond to positive directional selection, as was the case with floral nectar production in the Cave population and elaiosome mass in Saint Ann's Bay (Figure 5a). In both cases, this result reflects a combination of low genetic variance for these traits coupled with significant negative genetic correlations with other traits under positive selection (Table 1). Other populations, however, had more capacity to respond to selection on stigma–anther separation and floral nectar production due to the existence of positive genetic correlations that accelerated their response to selection (Figure 5a). These results are consistent with estimated R values (see Figure S4a), which provided evidence of genetic constraint in the evolution of FN production in Cave and elaiosome mass in the Saint Ann's Bay population. We observed similar variation in the capacity of our five populations to respond to multivariate selection that exerted opposing pressures on mutualistic traits. Here again, we found that misalignment between axes of genetic correlation and selection vectors reduced the capacity for trait evolution in certain populations (e.g. floral nectar in Cave, Murdock and Saint Ann's Bay), while the responses to selection for this same trait in other populations (Brown's Town and Mosquito Cove) were facilitated by the existence of either negative correlations among traits experiencing conflicting selection pressures or positive correlations among traits under the same selection pressures (Figure 5b,c, Table 1). We observed similar variation in the effects of genetic correlations on the evolution of stigma–anther separation, where alignment between axes of selection and covariation facilitated its response to selection in the Brown's Town, Mosquito Cove and Saint Ann's Bay populations and misalignment constrained its evolution in the Cave population (Figure 5b,c; Figure S4).

When we compared the size, structure and orientation of G, the average G matrix within populations that accounts for variation in sample sizes and trait means across sites, and D, the matrix summarizing multivariate divergence among the mutualistic traits of all 17 T. ulmifolia populations, we found evidence that the structure of G may have constrained the direction, but not magnitude, of phenotypic divergence at the metapopulation level. While we observed large, significant differences in the size of the D (mean with 95% HPD intervals: 1.09 [0.44–1.99]) and average G (0.24 [0.17–0.32]) matrices, the measure of their eccentricities was nearly identical (0.556 (0.454–0.634) and 0.533 (0.453–0.634) for D and G, respectively). Similarly, the multivariate directions of maximal genetic (g_max) and phenotypic (d_max) variation were nearly indistinguishable from one another, with an angle of only 0.327°(0.164°–2.95°) between them. When we assessed the extent to which our estimates of G and D shared a common substructure using Krzanowski's shared subspace analysis, we found further evidence of similarity. The eigenvalues of the eigenvectors of H, the summary matrix encapsulating the shared multivariate subspace between matrices, overlapped with the maximal value of 2 (p, the number of matrices being compared) for both interpretable eigenvectors (h₁: 1.94 (1.86–2.00), h₂: 1.70 (1.07–2.00)). This result strongly suggests a shared subspace among our G and D matrices (Aguirre et al., 2014; Krzanowski, 1979).

DISCUSSION

Genetic correlations among mutualist reward traits in T. ulmifolia have the capacity to affect the evolution of mating system and rewards for both pollinators and biotic defenders. We found significant positive genetic correlations among the production of FN, EFN and stigma–anther separation in the common genetic variance–covariance matrix G_W for the Jamaican metapopulation of T. ulmifolia. In other words, floral nectar production was coupled to extrafloral nectar production, and both were linked to variation in plant mating system, while elaiosome mass was not genetically correlated with other mutualistic traits. At the island‐wide level, these positive genetic correlations are predicted to accelerate T. ulmifolia's response to selection when selection favours increased investment in pollination, out‐crossing and biotic defence. However, opposing selection pressures on floral traits and biotic defence can constrain the evolution of mutualistic reward traits in T. ulmifolia. When we compared the multivariate axes of phenotypic divergence among our T. ulmifolia populations and those of multivariate genetic variation in the island metapopulation, we found that they were closely aligned, suggesting that the structure of G may have constrained the direction of phenotypic evolution in T. ulmifolia on Jamaica (but see below for alternative explanations for the alignment of G and D). Likewise, when we compared our estimates of G fit on individual populations from across the island, we found little evidence of differences in the orientation of G among populations.

Genetic variances and covariances among mutualistic traits

Selection on traits associated with multiple mutualisms likely varies across space and time, given that ecological outcomes are affected by community composition and nonadditive interactions among partners affect fitness outcomes of their hosts (Afkhami et al., 2014). The outcome of multiple mutualisms involving pollinators, biotic defenders and seed‐dispersing ants may be particularly prone to nonadditive interactions that alter the pattern of multivariate selection. Floral and extrafloral nectar are biochemically or physiologically linked (Dutton, Luo, et al., 2016; Heil, 2011), and selection on floral traits associated with out‐crossing and pollination can depend on both the presence of ant bodyguards that deter pollinators (Cembrowski et al., 2014; Villamil et al., 2019) and herbivores (Knauer & Schiestl, 2017; Ramos & Schiestl, 2019) that are in turn affected by the presence of ants. As such, our results may help explain the ubiquity and stable coexistence of both pollination by insects (Ollerton et al., 2011) and ant defence (Marazzi et al., 2013; Weber & Keeler, 2013), despite physiological trade‐offs and negative ecological interactions among partner species. A positive genetic correlation between floral and extrafloral nectar production may serve to maintain both mutualisms by reducing their visibility to negative selection and indirectly leading to the evolution of increased reward production through selection on correlated traits, even when the benefits of a particular mutualism may be context‐dependent or transient.

Our findings also hint at potential mechanisms responsible for generating genetic correlations between floral and extrafloral nectar production. Given the biochemical similarity of floral and extrafloral nectar (Heil, 2011), a positive correlation suggests shared genes are involved in nectar synthesis. Alternatively, a positive genetic correlation between floral and extrafloral nectar production might itself be adaptive. While extrafloral nectaries are thought to serve several adaptive functions, namely the attraction of bodyguards that protect against herbivory and, in our system, the attraction of seed‐dispersing ants (Cuautle et al., 2005; Cuautle & Rico‐Gray, 2003; Dutton, Shore, et al., 2016; Heil, 2011; Janzen, 1966), the distraction hypothesis has long been proposed as an alternate explanation for the evolution of extrafloral nectaries (Kerner, 1878). The distraction hypothesis posits that extrafloral nectar production is adaptive because it incentivizes ant recruitment away from flowers, thereby reducing ant consumption of floral nectar and ant‐pollinator conflict (Cembrowski et al., 2014; Villamil et al., 2018, 2019). Hence, a positive genetic correlation between floral and extrafloral nectar production in T. ulmifolia could be the result of pleiotropy, or of strong selection on plants with showy and nectar‐rich flowers to prevent ants from interfering with the pollination process.

Similarly, there may be an adaptive explanation for the positive genetic correlations we observed between stigma–anther separation and both floral and extrafloral nectar production in the Jamaican metapopulation of T. ulmifolia. Plants face a trade‐off between attracting pollinators and floral antagonists with showy floral displays, and multiple studies have shown that variation in floral morphology and mating system is more visible to selection when herbivores are absent (Knauer & Schiestl, 2017; Ramos & Schiestl, 2019; Thompson & Johnson, 2016). Both the deterrence of herbivores and the isolation of ant bodyguards away from floral tissues may thus be of greater selective value to genotypes with high stigma–anther separation, a quantitative trait correlated with out‐crossing rate in T. ulmifolia (Belaoussoff & Shore, 1995). Correlated selection on suites of traits such as this can lead to the evolution of stronger genetic correlations among them (Jones et al., 2004, 2012; Penna et al., 2017). Positive genetic correlations between stigma–anther separation, and pollinator and biotic defender rewards in T. ulmifolia are also consistent with the positive correlations between floral traits and extrafloral nectar production across species of Gossypium (Chamberlain & Rudgers, 2012). Such positive correlations may reflect a broader pattern of correlated selection on mating system and biotic defence in flowering plants.

Comparison of G matrices

As with other genetic entities, variation in G can arise through familiar evolutionary forces including heterogeneous selection, drift, migration and mutation (Björklund & Gustafsson, 2015; Guillaume & Whitlock, 2007; Hangartner et al., 2020; Jones et al., 2004; Roff, 2000). That we did not observe divergence among the G matrices of our five largest populations of T. ulmifolia suggests that the structure of G is relatively stable to these perturbations across Jamaica. These results are consistent with those of other studies that have reported stability in the shape and size of G within a single species, even across large geographic and climatic ranges (e.g. Hangartner et al., 2020; Puentes et al., 2016; Sniegula et al., 2018; Teplitsky et al., 2011).

We must, however, interpret evidence of similarity in G with some caution given the relatively low number of individuals (153–167) and maternal families (19–20) used to fit the G matrices of our largest populations and the conservative nature of G matrix comparison methods (Puentes et al., 2016; Sniegula et al., 2018). Local populations of T. ulmifolia, a highly opportunistic weed, tend to be small and ephemeral (Table S1, Barrett, 1978; Barrett & Shore, 1987); in many cases, we collected seeds from all or nearly all individuals in a natural population of T. ulmifolia, meaning we could not have phenotyped plants from additional maternal families for these populations. Although all methods found no significant differences in G among populations, our populations nonetheless differed in notable ways. G matrices fit for the Brown's Town and Mosquito Cove populations were larger than those for Cave, Murdock and Saint Ann's Bay and displayed more genetic variance along the principal axes of multivariate genetic variation among populations captured by the fourth‐order covariance tensor (Figure 4b,c). We also detected meaningful differences in evolutionary responses to selection among some populations (e.g. FN production in Cave and Mosquito Cove, Figure 5a; Figure S4). It is possible that such patterns may be explained by metapopulation dynamics in T. ulmifolia. Its seeds are often prone to extensive movement across the landscape given their secondary dispersal by ants and its use as an ornamental by humans (Bottcher et al., 2016; Gómez & Espadaler, 1998). Bottleneck events such as those associated with stochastic colonization events can affect the structure of G (Whitlock et al., 2002) and affect the relative genetic contribution of local populations to a connected metapopulation (e.g. hard selection; Wade, 1985; De Lisle & Svensson, 2017). However, local variation in G may not persist in the face of extensive gene flow in a connected metapopulation.

Evolutionary constraint and facilitation through G

When ecologically relevant traits are genetically correlated, as in T. ulmifolia, the evolutionary responses of populations are determined not only by genetic variance in and selection on individual phenotypes, but also by the genetic covariances among traits, which can facilitate or constrain adaptive evolution depending on whether axes of genetic covariance are aligned with those of selection (Agrawal & Stinchcombe, 2009; Teplitsky et al., 2014; TerHorst et al., 2018; Wise & Rausher, 2013). We found that the evolution of FN, EFN and stigma–anther separation will be facilitated when selection favours increased investment in pollination, out‐crossing and biotic defence (Figures 5a; Figure S4a). Conversely, under selection regimes consistent with selection for selfing, favouring increased inducible defences and dispersal and reduced FN production and stigma–anther separation (Baker, 1955; Belaoussoff & Shore, 1995; Campbell et al., 2014; Campbell & Kessler, 2013; Pannell & Barrett, 1998), we found that both EFN and FN showed little or no response to simulated selection in the Jamaican metapopulation of T. ulmifolia (Figure 5b,c; Figure S4b,c). In other words, plant mutualisms with pollinators and ants are unlikely to evolve independently in T. ulmifolia, but whether adaptive evolution is accelerated or hindered by G will depend on the nature of multivariate selection. The response to directional selection acting on one trait would, for example, differ depending on whether stabilizing selection is resisting change in other correlated traits (i.e. ‘conditional evolvability’, Hansen & Houle, 2008).

Although G did not differ significantly among populations, there was extensive variation in the means of multiple mutualism‐associated traits among populations (Figure 3), a pattern that could be indicative of genetic constraint on the direction of phenotypic evolution. The matrix of multivariate phenotypic divergence among population trait means, D, and G was nearly identical in eccentricity and orientation, with only a 0.33° angle between the directions of maximal genetic variation (g_max) and phenotypic divergence (d_max). These results may be evidence of genetic constraint, with phenotypic divergence among populations being restricted to genetic lines of least resistance, that is, the prominent multivariate axes of genetic variation (Schluter, 1996). Multiple previous studies have demonstrated that G can constrain phenotypic divergence among populations, including in populations separated by greater evolutionary (e.g. McGlothlin et al., 2018) and geographic (e.g. Silva et al., 2020) distances than ours. Conversely, gene flow across the island of Jamaica may also explain the alignment we observed between patterns of multivariate genetic variance and the direction of phenotypic divergence; simulations have shown that high gene flow can realign G by increasing its genetic variance in the direction of phenotypic divergence (Guillaume & Whitlock, 2007). It is, however, also possible that the close alignment of D and G we observed is not a product of one's influence on the other, but rather of both being shaped by an adaptive landscape that is capable over time of shaping both G and D along selective lines of least resistance (Arnold et al., 2001). Distinguishing between these alternative hypotheses remains difficult in our case, especially given that the traits we studied can mediate multiple ecological interactions (e.g. EFN affects defence and seed dispersal), and multiple mutualisms are often context‐dependent, which could result in precisely the kind of correlated selection regime capable of bringing both G and D into line with the adaptive landscape (Arnold et al., 2001; McGlothlin et al., 2018; Silva et al., 2020).

CONCLUSION

The evolution of traits mediating multiple species interactions will be influenced by both ecological complexity and their underlying genetic architectures. Positive genetic correlations among FN, EFN and stigma–anther separation in T. ulmifolia have the potential to link the evolutionary fates of these traits in a metapopulation. These findings highlight the importance of estimating the genetic correlations among traits characterized by ecological (Cembrowski et al., 2014; Ramos & Schiestl, 2019; Villamil et al., 2018, 2019) and physiological trade‐offs (Dutton, Luo, et al., 2016). Despite the potential for conflict between pollination and biotic defence, we found that the production of nectar rewards for these two mutualisms was positively correlated at the metapopulation level and uncorrelated within individual smaller populations. Genetic correlations among mutualistic traits may be a result of pleiotropy, but may also plausibly reflect a process of coordinated positive selection on a suite of phenotypes. These patterns appear to be relatively stable in Jamaican T. ulmifolia, as we found weak evidence of local variation in G among populations across the island. The stability of G suggests that local variation is unlikely to impact the evolution of multiple mutualisms at the metapopulation level and that phenotypic divergence in mutualistic traits may be constrained by the prominent axes of multivariate genetic variation in a metapopulation. Our results provide rare insight into the genetic architecture of multiple mutualisms, highlighting the potential for genetic correlations among reward traits to affect the evolution and ecology of these commonplace and essential interaction networks.

AUTHOR CONTRIBUTIONS

JRL and MEF designed the study. JRL collected seeds and grew plants; JRL, CGR, CB, TW, and CK collected trait data. CGR, JRL, and MEF designed the hand‐pollination experiment; CGR collected data. JRL analysed data and drafted the manuscript; JRL and MEF finalized the text.

ACKNOWLEDGEMENTS

We would like to thank many undergraduate students for their assistance in plant care and data collection, including A. Copeland, A. Wichert, B. Garden‐Smith, C. Chen, C. Qiao, E. Zhang, F. Samada‐zada, J. Han, K. Ong, M. Giblon, M. Gorchkova, S. Fernandes, and S. Ravoth; S. Barrett for guidance in planning field collections; J. Shore and E. Dutton for advice on working with Turnera ulmifolia; J. Stinchcombe, J. Sztepanacz, and C. Wood for assistance in fitting and comparing G matrices; J. Stinchcombe for providing comments on our manuscript; and A. O'Brien and B. McGoey for helping with coding. M.E.F. received funding from an NSERC Discovery Grant and the University of Toronto; J.R.L. was supported by an NSERC CGS‐D Alexander Graham Bell scholarship; C.G.R. was supported by an NSERC Undergraduate Student Research Award.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

DATA AVAILABILITY STATEMENT

The data and sample code that support the findings of this study are openly available on Zenodo, at https://zenodo.org/badge/latestdoi/471031478.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.14098.