Genetic and environmental integration of the hawkmoth pollination syndrome in Ruellia humilis (Acanthaceae)

Abstract

Background and Aims The serial homology of floral structures has made it difficult to assess the relative contributions of selection and constraint to floral integration. The interpretation of floral integration may also be clouded by the tacit, but largely untested, assumption that genetic and environmental perturbations affect trait correlations in similar ways. In this study, estimates of both the genetic and environmental correlations between components of the hawkmoth pollination syndrome are presented for chasmogamous flowers of Ruellia humilis, including two levels of control for serial homology.

Methods A greenhouse population for quantitative genetic analysis was generated by a partial diallel cross between field-collected plants. An average of 634 chasmogamous flowers were measured for each of eight floral traits that contribute to the hawkmoth syndrome. Genetic correlations (across parents) and environmental correlations (across replicate flowers) were estimated by restricted maximum likelihood.

Key Results Stigma height, anther height and floral tube length were very tightly integrated in their responses to both genetic and environmental perturbations. The inclusion of floral disc width as a control for serial homology suggests this integration is an adaptive response to correlational selection imposed by pollinators. In contrast, integration of non-homologous traits was low. Furthermore, when comparisons between the dimensions of serially homologous structures were excluded, the genetic and environmental correlation matrices showed little congruence.

Conclusions The results suggest that hawkmoths have imposed strong correlational selection on floral traits involved in the deposition and removal of pollen, and that this is a consequence of stabilizing selection on the relative positions of stigmas and anthers in the face of substantial flower size variation. Low integration of other floral traits, and conflicting patterns of genetic and environmental correlations among these traits, suggest weak or no correlational selection within the range of variability expressed within a population.

Correlational selection, floral integration, flower size variation, genetic correlation, hawkmoth pollination, herkogamy, phenotypic integration, Ruellia humilis

INTRODUCTION

Selection imposed by pollinating animals is thought to have played a major role in the diversification of flowers among angiosperm species, resulting in different suites of floral traits (pollination syndromes) that are consistently associated with different guilds of pollinators (Darwin, 1862; Stebbins, 1950; Baker, 1963; Grant and Grant, 1965; Fægri and Van der Pijl, 1966; Fenster et al., 2004). As a consequence, when comparisons are made across species, correlations are apparent between floral traits, a pattern that has been referred to as selective correlation (Stebbins, 1950), selective covariance (Felsenstein, 1988) or evolutionary integration (Cheverud, 1996a). Couched in genetic terms, diversifying multivariate selection generates linkage disequilibrium across species (evolutionary integration), and reproductive isolation prevents the disequilibrium from decaying.

Evolutionary integration does not depend on physiological or developmental interactions among traits and therefore provides little information about such interactions (Felsenstein, 1988; Armbruster and Schwaegerle, 1996). In contrast, substantial correlations between traits within individual populations of partially or fully outcrossing species necessarily reflect either tight linkage between genes affecting the traits (preventing the decay of linkage disequilibrium) or physiological/developmental interactions among traits (phenotypic integration). Thus, in an effort to understand floral integration, much recent work has focused on correlations between floral traits within individual populations (see review by Armbruster et al., 2014, and accompanying articles).

Most work on floral integration has been motivated by two competing hypotheses for the significance of phenotypic integration. One hypothesis proposes that integration is an evolved mechanism for maintaining optimal function in the face of genetic and environmental perturbations, and is therefore adaptive (e.g. Olson and Miller, 1958; Cheverud, 1982; Conner and Via, 1993; Wagner, 1996). The alternative hypothesis proposes that correlations among traits are non-adaptive, reflecting homologies, developmental side-effects or the multivariate distribution of mutational effects (e.g. Lande, 1984; Lynch 2007; Armbruster et al., 2014). According to this second hypothesis, covariances among traits are of evolutionary significance primarily as constraints on the response to multivariate selection (Arnold, 1992; Cheverud, 1984). Because conserved developmental pathways are expected to generate correlations between the dimensions of serially homologous flower parts, significant correlations between such traits do not imply adaptive integration (Conner and Via, 1993; Waitt and Levin, 1993; Conner and Sterling, 1996; Herrera et al., 2002). Unfortunately, studies of floral integration have focused primarily on such traits (see reviews of Conner et al., 2014 and Fornoni et al., 2016), making it difficult to distinguish between selection and homology as explanations for floral integration.

One method of controlling for serial homology is to compare patterns of floral trait correlation between populations or species that ostensibly differ in the pattern or strength of selection imposed by pollinators. Comparative studies suggest that plant species specializing on relatively few pollinator guilds, and therefore presumably experiencing relatively strong stabilizing selection on flower structure, tend to display higher levels of floral integration than species that are pollinator generalists or wind-pollinated (Ashman and Majetic, 2006; Rosas-Guerrero et al., 2011; Gómez et al., 2014; Pérez-Barrales et al., 2014), suggesting some degree of adaptive floral integration. However, the observed differences in the level of floral integration between pollinator specialists and pollinator generalists are often modest (e.g. Rosas-Guerrero et al., 2011; Gómez et al., 2014) and sometimes lacking (Armbruster et al., 1999).

Several refinements to this comparative approach have failed to add much clarity. By comparing patterns of integration with patterns of selection imposed by current pollinators, Lázaro and Santamaría (2016) found significant congruence, whereas Nattero et al. (2010) did not. Pérez et al. (2007) found significant changes in the pattern of floral integration across species of Schizanthus that are associated with transitions to autogamy and changes in pollinator behaviour, but comparisons among populations within a species have either found no evidence for divergence in the variance–covariance matrix (Waitt and Levin, 1993; Armbruster et al., 2004) or significant divergence that is either not associated (Herrera et al., 2002) or only weakly associated (Pérez-Barrales et al., 2007) with pollinator divergence.

A more direct method of controlling for homology is to examine a set of serially homologous traits within a single population, some of which have clear pollination functions and others of which lack pollination functions. Using this approach, Hansen et al. (2007), Pélabon et al. (2011) and Wanderley et al. (2016) found evidence for higher integration of floral traits with pollination functions. Another direct method for removing the effects of homology is to quantify integration among non-homologous floral traits with clear pollination functions, but the heavy focus of prior work on the dimensions of floral traits has largely precluded such comparisons.

In addition to the problem of homology, a failure to separate phenotypic correlations into their genetic and environmental components may further obscure the role of selection in shaping floral integration. Although phenotypic and genetic correlations are frequently concordant (Cheverud, 1988; Roff, 1996; Waitt and Levin, 1998; Dochtermann, 2011), this is a mathematical necessity when heritabilities are high (Willis et al., 1991), as is often the case in laboratory and greenhouse experiments. Direct comparisons of genetic and environmental correlations in animal species have revealed moderate to weak concordance in some cases (Hegman and DeFries, 1970; Cheverud, 1996b; Klingenberg and McIntyre, 1998; Willmore et al., 2005; Breuker et al., 2006; Klingenberg et al., 2010) and pronounced discordance in others (Cheverud, 1982; Hadfield et al., 2007; Breno et al., 2011). Similar studies in plants have revealed weak concordance in one case (Venable and Búrquez, 1990) and discordance in another (Pélabon et al., 2013). Thus, the extent to which genetic and environmental perturbations affect the same developmental or physiological pathways, and thereby generate similar patterns of covariation among traits, may vary across traits and species in ways that are as yet poorly understood (Willis et al., 1991; Lynch and Walsh, 1998), especially for plants. In animal studies, estimates of environmental correlations are usually either extracted from data on fluctuating asymmetry (e.g. Klingenberg and McIntyre, 1998) or calculated indirectly from estimates of phenotypic and genetic correlations (Cheverud, 1982). Although the relation between genetic and environmental correlations is virtually unstudied in plants, the relative ease with which environmental correlations can be estimated across clones (Venable and Búrquez, 1990) or metamers (Williams and Conner, 2001) makes plants particularly amenable to such studies.

Here we present data on variances and correlations within a suite of eight floral traits that are predicted to influence the ability of hawkmoths to effect cross-pollination of chasmogamous flowers in Ruellia humilis, a long-lived herbaceous perennial with mixed mating. Our purpose is to test the hypothesis that this guild of pollinators has selected for canalization and integration of the hawkmoth pollination syndrome. Our study of floral integration is unique in that (1) we include important components of the hawkmoth syndrome that are not serially homologous and are therefore less likely to display genetic correlations due to conserved developmental pathways, and (2) we estimate both genetic and environmental correlations in order to compare patterns of floral integration in the face of these two sources of developmental perturbation.

MATERIALS AND METHODS

Study species

Ruellia humilis Nutt. (Acanthaceae) is a highly variable, herbaceous perennial that is native to the eastern and midwestern USA from Florida to Pennsylvania and from Texas to Minnesota, where it may be encountered in a variety of habitats, including dry prairies, dry woodlands, glades, bluffs, sandy soils and roadways (Long, 1961; Steyermark, 1963). Like many of its congeners, R. humilis is self-compatible and has a mixed mating system, producing both large, open chasmogamous (CH) flowers that presumably evolved to promote outcrossing by attracting animal pollinators, and small, closed cleistogamous (CL) flowers that automatically self-pollinate. The corollas of CH flowers are weakly zygomorphic and funnelform, with five distinct lobes, and are lavender in appearance to the human eye with a dark purple nectar guide on the lower corolla lobe (Fig. 1). Nectar extracted from CH flowers has a sweet, fruity aroma when placed directly under the nose, but the fragrance is very subtle in the vicinity of intact flowers.

Fig. 1.

Morphology of CH flowers of Ruellia humilis in Missouri. (A) The corolla is weakly zygomorphic with anthers appressed to the top of the throat and a nectar guide on the bottom of the throat. (B) The lower three corolla lobes have been removed to reveal the four didynamous stamens with their filaments inserted at the top of the floral tube, where they encircle the style. (C) A corolla of average size (bottom) compared with one at the low end of the size range (top).

The shape and pigmentation of CH flowers are highly variable in expression, as are vegetative traits, motivating several authors to recognize subspecific taxa (Fernald, 1945; Tharp and Barkley, 1949; Steyermark, 1963; Turner, 1991). Geographical variation in floral traits might reflect local adaptation to different pollinators, but little is known about the pollination ecology of this species. Tripp and Manos (2008) reported that R. humilis in Pennsylvania is visited by bees, butterflies and diurnal moths, but the effectiveness of these visitors as pollinators was not investigated. In south-western Missouri, CH corollas are characterized by pale coloration, long floral tubes, erect presentation, nocturnal anthesis and significant nectar volumes, with ∼20 % sugar content, a suite of traits associated with pollination by hawkmoths (Lepidoptera: Sphingidae) (Percival, 1969; Cruden et al., 1983; Grant, 1993; Willmer, 2011). Visitation by the sphingid Hyles lineata has been reported on Konza Prairie in Kansas (D. K. Townsend, Manhattan, KS, USA, pers. comm.), but over 20 years of field observations suggest a permanent loss of hawkmoth visitation in the highly fragmented remnants of tallgrass prairie in Missouri (J. S. Heywood, unpubl. data). Although molecular markers suggest that extant individuals of this long-lived perennial were generated by mixed mating, current mating system estimates indicate that CH flowers receive primarily self-pollen as a consequence of the activities of pollen-harvesting halictid bees (Smith, 2002).

In our study population the CH corollas persist for less than a day, typically dehiscing before noon in the field. The four filaments in a CH flower enclose the style at their point of insertion at the top of the floral tube, ensuring that the stigma is held directly above the anthers and is drawn across the anthers when the corolla dehisces, thus promoting delayed selfing (J. S. Heywood, unpubl. data). Chasmogamous and CL fruits are identical in appearance and contain a maximum of four seeds. Greenhouse plants produce both CH and CL flowers in their first year and display the same flowering phenology as plants in the field, facilitating short-term studies in the greenhouse despite a potential life span of 100 years or more (J. S. Heywood, unpubl. data).

Source population and partial diallel crossing design

Seeds of R. humilis were collected from Mount Vernon Prairie Natural Area (MVP) in Lawrence County, MO, USA (37°8′21″ N, 93°47′21″ W). The soils at this site are derived from cherty limestone and support an upland dry prairie (Davit, 1999). Mature CH fruits were collected in September 2010 from 16 plants widely dispersed within the south-western quarter (∼4 ha) of MVP. After winter stratification at 4 °C, seeds were germinated in the spring of 2011, transplanted individually into 10·6-L plastic tree pots filled with Fafard Growing Mix No. 2 (Sun Gro Horticulture, Agawam MA, USA), and subsequently maintained in a greenhouse on the campus of Missouri State University in Springfield, MO. In the summer of 2012 we selected one healthy plant from each family to be used as parents in a partial diallel crossing design. Because the families were derived from widely spaced natural plants, it is assumed that these 16 individuals are unrelated for the purposes of data analysis.

Each of the 16 parents was crossed with four others in reciprocal, yielding 64 families of full sibs. Prior to pollen transfer, pollen recipients were emasculated and examined with a hand lens to ensure the absence of pollen on the stigma. Based on the success of trial crosses performed in 2011, all crosses were completed between 0500 and 0700 h. We repeated each cross as many times as feasible during the flowering season. Cross-compatibility was highly variable, with the number of outcross seeds produced ranging from 0 to 20 (mean 8·6) per family.

Seed germination was initiated on 24 May 2013. In early June, five seedlings (or all seedlings if fewer than five were available) were selected from each family and transplanted to 2·5-L plastic pots filled with Fafard No. 2. If fewer than five seedlings were available for a family, additional seedlings (as available) were selected from the reciprocal family up to a maximum total of ten for both families, with the goal of keeping the design as balanced as possible when reciprocals were combined. The final offspring population consisted of 264 plants representing 55 of the 64 crosses performed and 30 of the 32 outcross families (reciprocals combined). These 264 seedlings were randomized in the greenhouse.

Chasmogamous flowers produced by these plants in the summer of 2013 were used to quantify floral traits. Flowers intermediate between CH and CL are occasionally produced, becoming more common as CL flowers become frequent in September. These were avoided by including only CH flowers with fully expanded corollas of typical shape. Throughout the 2013 flowering season, the temperature and humidity inside the greenhouse were kept close to external ambient values by exchanging the air at a high rate.

Floral traits

Eight floral traits that are components of the hawkmoth pollination syndrome were measured on an average of 2·4 flowers per plant (Table 1). These included four dimensional measurements taken on serially homologous floral whorls [corolla disc width (DW), floral tube length (TL), stigma height (SH), anther height (AH)], one relational measurement [stigma–anther separation (SA)] and three traits not expected to show homologies with the dimensional traits [nectar volume (NV), nectar concentration (NC), bud break (BB)]. Flower dimensions were recorded from digital photographs using ImageJ (Rasband, 1997–2014). Each photograph included a scale placed in the plane of the desired measurement and was taken with the camera perpendicular to the plane of the scale. Corolla disk width was measured as the greatest distance between tips of the lateral corolla lobes (Fig. 1A). To measure SA, the lower and two lateral lobes of the corolla disk and throat were cut away to the top of the corolla tube without detaching the flower from the top of the ovary. We recorded SA as the distance between the top of the stigma and the top of the anther on the longest stamen, parallel to the floral tube (Fig. 1B). Floral tube length was measured from the bottom of the detached corolla to the point of insertion of the filaments, which is the point where the throat flares (Fig. 1C). Anther height was measured from the bottom of the floral tube to the top of the anther of the longest stamen. Stigma height was calculated as AH+SA. After photographing a flower, nectar was extracted by inserting 562-μm internal diameter acrylic capillary tubing (CT562-750-5, Paradigm Optics, Vancouver, WA, USA) into the base of the floral tube and measuring the length of the column to the nearest 0·1 mm with digital callipers. The nectar was then expressed onto the prism of a temperature-corrected Brix refractometer (model RF15, Extech Instruments, Nashua NH, USA) and the sucrose equivalent was recorded to the nearest 0·1%.

Table 1.

Sample mean, sample standard deviation, sample CV and range for eight CH floral traits measured on first-year shoots of 264 individuals of R. humilis grown in a greenhouse. Time of bud break (BB) is reported as minutes past midnight, Central Daylight Savings Time. CV was not calculated for traits with negative values (SA) or an arbitrary zero point (BB)

Floral trait	n	Mean	s.d.	CV	Minimum	Maximum
DW (corolla disk width), mm	644	37·8	6·7	0·177	18·2	56·9
TL (floral tube length), mm	650	36·2	5·9	0·163	19·6	53·7
SH (stigma height), mm	634	52·9	7·4	0·141	32·8	75·9
AH (anther height), mm	634	48·5	7·0	0·145	28·0	69·6
SA (stigma–anther separation), mm	636	4·4	2·3	–	−4·3	14·1
NV (nectar volume), μL	651	2·3	1·0	0·447	0·0	6·0
NC (nectar concentration), %	633	19·3	1·8	0·095	12·4	27·6
BB (time of floral bud break), min	590	119	60	–	−90	300

The time (Central Daylight Savings Time) at which CH flowers begin to open (BB) was determined by placing plants with mature buds under a time-lapse infrared camera (PC800 Hyperfire Professional, Reconyx, Holmen WI, USA) set to record an image every 5 min. Up to 16 plants could be accommodated under the camera on the same evening. Bud break occurs rapidly and can be assigned to a 5-min interval without ambiguity.

To assess whether NV or NC varies with time of day, the diurnal dynamics of these variables were examined by sampling multiple CH flowers from each of 11 large plants between 0200 h (mean BB) and 0900 h on two separate days in July 2014. To examine the effect of plant age on NV and NC, two to four CH flowers (mean 3·0) were sampled from each of 18 plants in late July 2015 (third-year shoots) and compared with data from the same plants in 2013 (first-year shoots).

Statistical analysis

We estimated phenotypic variance and covariance components for pairs of traits using restricted maximum likelihood (REML) as implemented in PROC MIXED of SAS 9·2 (www.sas.com). Following Möhring et al. (2011), we used dummy variables to ensure that the same random effects were estimated for a parent as both female and male. We used the multivariate approach described by Holland (2006) to estimate phenotypic variance and covariance components for pairs of traits. The variance–covariance structure for a pair of traits was estimated simultaneously across parents, across plants within full-sib families, and across flowers within plants, thus providing an estimate of the correlation between the traits at each of these levels. Although the inbreeding coefficient f of the parent population is unknown, the additive genetic variances and covariances are both inflated by a factor of (1+f) (Kempthorne, 1957) so that the genetic correlations are not affected. Thus, the phenotypic correlation across parents [the correlation in general combining ability (GCA)] approximates the additive genetic correlation, assuming non-additive, maternal and paternal effects are negligible. The correlation across plants within a family includes both genetic and environmental effects, and it is not possible to separate these components without knowledge of parental inbreeding coefficients. The phenotypic correlation across flowers within a plant is entirely environmental. Patterns of trait variation and covariation across flowers often reflect regular patterns of ontogenetic variability associated with flower age, plant age or location within an inflorescence (Diggle, 2014), in which case the environmental variation that is responsible for floral trait variation is orchestrated by the plant and thus non-random. However, these architectural effects are unlikely to contribute to floral trait variability in R. humilis since flowers are borne individually and last <1 d. Our failure to detect any relation between flower size and shoot age (see Results section) further supports this conclusion. Thus, among-flower variances and covariances probably reflect random, local variations in internal and/or external environmental variables, i.e. a lack of canalization.

Each of these correlations (genetic, across siblings and environmental) was tested for statistical significance by fitting a reduced model in which the covariance of interest was set to zero and then performing a likelihood ratio test. Prior to analysis, all traits were standardized to the same mean and variance to facilitate convergence of the Newton–Raphson algorithm (Holland, 2006). In the few instances where the algorithm failed to converge on a solution for a reduced model, we report the less powerful Wald z-test that is provided by PROC MIXED. Complete SAS code and explanations can be found in Supplementary Data Table S1.

Because of the unbalanced design, an estimate of the total phenotypic correlation between two traits was obtained by first summing the across-parents, across-siblings and within-plants REML estimates for the first trait variance, the second trait variance, and the covariance, and then using these summed variances and covariance to calculate an estimated correlation. The squared standard error of the estimated total phenotypic covariance was calculated as the sum of the squared asymptotic standard errors provided by PROC MIXED. A z-test was then used to test for significance of the total phenotypic covariance (and correlation).

Two measures of phenotypic integration have commonly been reported in the literature: the sample mean of the absolute values of all pairwise correlation coefficients, $\overline{\left|r\right|}$⁠, and the sample variance of the eigenvalues of the sample correlation matrix, V(λ). We calculated both statistics for the genetic correlation matrix, the within-plant environmental correlation matrix and the total phenotypic correlation matrix. Reported values for V(λ) are corrected for bias and expressed as a fraction of their theoretical maximum (Wagner, 1984; Cheverud et al., 1989). The value for V(λ) obtained from the total phenotypic correlation matrix derived from the REML estimates was identical (to three significant digits) to the value obtained by calculating sample correlations across all flowers and jackknifing V(λ) across the 16 original parents, so we also report the jackknifed estimate of the standard error of V(λ) for the total phenotypic correlation matrix based on the latter approach.

We tested the null hypothesis of independence between the genetic and environmental correlation matrices by calculating the correlation between their off-diagonal elements (matrix correlation) and obtaining a P-value by permutation of the rows and columns of one of the matrices (Mantel and Valand, 1970; Dietz, 1983). Although many procedures have been proposed for comparing matrices (Roff et al., 2012), matrix correlation followed by permutation is the most commonly employed test for equality in studies of phenotypic integration, and it performs well (Marroig and Cheverud, 2001; Roff et al., 2012). To avoid biasing the matrix correlation we omitted SA from this analysis due to its redundancy with AH and SH. We performed all of the 7! = 5040 possible permutations rather than sampling at random from this set.

For purposes of comparison with other species, we report estimated coefficients of variation (CV) for individual traits. For each trait an estimate of the total phenotypic variance was obtained using a general linear model and the method of moments to partition the variance into between-cross, between-plant (within cross) and between-flower (within plant) components, and then summing these three components.

All analyses other than the REML variance–covariance estimates were performed with Minitab 17 (www.minitab.com). All conclusions are based on a type-I error rate of 0·05.

RESULTS

Nectar dynamics

The sample mean of NV increased significantly between midnight and 0900 h, primarily due to lower than average values prior to 0400 h (Fig. 2). An analysis of covariance utilizing data collected after 0400 h suggests a weak effect of time, with most of the variance in NV attributed to differences among plants (46 %), among flowers within plants (43 %) and a plant × time interaction (8 %) (Supplementary Data Table S2). These data suggest that nectar secretion is essentially complete by 0400 h and therefore that NV values were not affected much by the time of measurement. The sample mean NV was significantly greater in the third year of life than in the first year (3·87 versus 2·41 µL) for 18 unrelated individuals (repeated-measures ANOVA, P<10⁻⁵ for main effect of age), suggesting that data from first-year plants underestimate typical nectar production.

Fig. 2.

Nectar volume as a function of time of day for 112 CH flowers on 11 greenhouse plants, with a different symbol for each plant. The line is the least-squares linear regression model for all flowers (R²=17·5 %, P <10⁻⁵).

Nectar concentration increased slightly between midnight and 0900 h (linear regression, P = 0·023, R² = 4·9 %), but no trend was apparent after 0400 h (P = 0·40, R² = 1·2 %). The sample mean NC was slightly lower in the third year of life than in the first year (18·3 versus 19·4 %) for 18 unrelated individuals (repeated-measures ANOVA, P = 0·039 for main effect of age).

Variability of floral traits

Chasmogamous flower size was highly variable in the greenhouse, with the largest flowers exceeding the smallest by a factor of nearly 3 in DW, TL, SH and AH (Table 1, Fig. 1C). Mean flower size was not significantly correlated with time of year for the sample of flowers that were measured (e.g. r  = −0·043 for DW, n = 644, P = 0·72).

Stigma–anther separation ranged from −4·3 to 14·1 mm, with a mean of 4·4 mm (Table 1). Because of the way SA was measured, negative and zero values within the observed range represent stigmas in direct contact with some portion of one or more anthers, whereas positive values represent stigmas elevated above the most distal anther (‘approach herkogamy’; Webb and Lloyd, 1986). Of the 636 flowers for which SA was measured, 3·9 % lacked herkogamy and 96·1 % displayed approach herkogamy. Sample mean sizes of herkogamous flowers (⁠$\overline{DW}$ = 38·0 mm, s.e. = 0·28) and non-herkogamous flowers (⁠$\overline{DW}$ = 34·9 mm, s.e. = 1·03) were significantly different (two-sample t-test with Satterthwaite’s correction for unequal variances, P = 0·007), suggesting that smaller flowers are slightly more likely to lack herkogamy (Fig. 3).

Fig. 3.

Stigma–anther separation (SA) as a function of flower size for 628 CH flowers. The solid line is the least-squares linear regression model (R² = 5·8 %, P < 10⁻⁸). When SA is negative the stigma is in direct contact with one or more anthers.

Nectar volume was highly variable (CV = 0·45) while NC was less so (CV = 0·10), although the latter still displayed a large range of values (Table 1).

Flower buds began to open (BB) as early as 2230 h and as late as 0500 h, but nearly 96 % of all flowers opened between midnight and 0400 h, with a circular mean time of ∼0159 h (Table 1). Bud break was negatively correlated with temperature at 0200 h (r = −0·37, n = 590, P<0·001), the latter ranging from 23 to 27 °C. On average, BB was delayed 24 min for each 1 °C reduction in temperature.

Genetic and environmental correlations

The total phenotypic correlations between the four serially homologous dimensional traits (AH, SH, TL, DW) were all strong, positive and statistically significant (Table 2, group 1). Both the genetic and environmental correlations between AH, SH and TL were very high (0·922–0·972) and statistically significant. The environmental correlations between these three traits and DW were also significant but lower in magnitude (0·684–0·718), and the genetic correlations were lower still (0·342–0·488) and not statistically significant. The only other significant genetic correlations were between SA and DW (r = 0·809) and between SA and NC (r = 0·617) (Table 2, group 2). However, the genetic correlations were estimated with relatively low precision due to the modest number of parents.

Table 2.

Estimated phenotypic correlations generated by different sources of variation. Sample sizes for correlations involving BB range from 462 to 476; sample sizes for all other correlations range from 615 to 649. Statistically significant correlations (P < 0·05) are shown in bold font. Standard errors are not reported for the estimated correlations because the computer time required to jackknife the REML variance–covariance estimates was prohibitive. The five groups of trait pairs are referenced in the text. Floral traits are identified as in Table 1

Trait pair	Genetic (GCA)			Among full sibs			Within plant			Total phenotypic
	r	P		r	P		r	P		r	P
Group 1
AH–SH	0·968	<10−6		0·942	0·0001		0·954	<10−6		0·953	<10−6
SH–TL	0·954	<10−6		0·935	0·0001		0·922	<10−6		0·926	<10−6
AH–TL	0·972	<10−6		0·983	4 × 10−5		0·970	<10−6		0·972	<10−6
DW–SH	0·488	0·098		0·710	0·0018		0·718	<10−6		0·700	<10−6
DW–AH	0·349	0·26		0·640	0·0026		0·711	<10−6		0·674	<10−6
DW–TL	0·342	0·26		0·602	0·0033		0·684	<10−6		0·644	<10−6
Group 2
SA–AH	0·293	0·38		−0·219	0·24		0·042	0·41		0·010	0·86
SA–SH	0·524	0·09		0·102	0·62		0·340	<10−6		0·311	<10−6
SA–TL	0·304	0·34		−0·169	0·36		0·025	0·63		0·009	0·88
SA–DW	0·809	0·010		0·121	0·53		0·170	0·0008		0·198	0·0006
SA–NV	0·288	0·40		−0·278	0·14		0·162	0·0011		0·088	0·11
SA–NC	0·617	0·037		0·148	0·57		0·212	3 × 10−5		0·231	0·0001
SA–BB	–0·056	0·86		0·455	0·042		–0·074	0·21		0·016	0·81
Group 3
NV–NC	0·243	0·45		−0·171	0·53		0·263	<10−6		0·207	0·0004
NV–AH	–0·110	0·74		0·500	0·014		0·521	<10−6		0·477	<10−6
NV–SH	−0·042	0·90		0·423	0·048		0·535	<10−6		0·478	<10−6
NV–TL	0·003	0·99		0·567	0·0043		0·472	<10−6		0·455	<10−6
NV–DW	−0·173	0·60		0·517	0·012		0·557	<10−6		0·507	<10−6
Group 4
NC–AH	−0·118	0·72		−0·165	0·53		−0·036	0·48		−0·059	0·32
NC–SH	0·051	0·87		−0·109	0·70		0·038	0·46		0·021	0·72
NC–TL	−0·020	0·95		−0·251	0·33		−0·054	0·29		−0·077	0·19
NC–DW	0·545	0·087		−0·196	0·46		0·084	0·098		0·080	0·17
Group 5
BB–NV	−0·265	0·40		−0·715	0·0004		−0·322	<10−6		−0·381	<10−6
BB–NC	0·133	0·65		0·283	0·35		−0·375	<10−6		−0·249	0·0001
BB–AH	−0·320	0·31		−0·645	0·0021		−0·082	0·17		−0·194	0·0024
BB–SH	−0·292	0·35		−0·555	0·014		−0·097	0·10		−0·183	0·0044
BB–TL	−0·484	0·11		−0·618	0·0023		−0·071	0·23		−0·202	0·0019
BB–DW	−0·293	0·33		−0·634	0·0020		−0·205	0·0003		−0·280	7 × 10−6

Nectar volume had a significant and positive environmental correlation with NC and with each of the four dimensional traits (Table 2, group 3). In contrast, the genetic correlations for these five pairs suggest no such pattern.

Significant negative phenotypic correlations were found between BB and all of the other traits except SA (Table 2, group 5). Although the negative genetic correlations between BB and the four dimensional traits were not statistically significant, the large and significant correlations across full sibs, combined with weak environmental correlations within plants, suggest significant genetic contributions to these correlations. The negative correlations indicate that larger flowers tended to open earlier.

To eliminate bias when comparing genetic and environmental correlations, we removed SA from the data set for this analysis because it is redundant with AH and SH. Based on 21 pairwise comparisons between the remaining seven traits, there was a significant matrix correlation between estimated genetic correlations and estimated within-plant environmental correlations (r = 0·736, P = 0·008; Fig. 4). However, if the correlations between serially-homologous traits (shown as open circles in Fig. 4) are eliminated, the relationship is much weaker (r = 0·195). Within this latter group of 15 trait pairs, the sample mean magnitude of the estimated correlation was similar for the genetic (0·21) and environmental (0·25) components.

Fig. 4.

Relation between estimates of the genetic correlation and the within-plant environmental correlation for 21 pairs of traits. Stigma–anther separation (SA) was excluded because of its redundancy with AH and SH. Correlations between the dimensions of serially homologous floral structures (AH, SH, TL, DW) are shown as open circles.

On average, the magnitude of the within-plant covariance was 2·9 times the magnitude of the combined across-parents and across-sibs covariances (Fig. 5), the latter quantity being equal to (1−f)/2 times the genetic covariance plus the between-plant environmental covariance. Thus, on average, the within-plant environmental covariance was >3 times the between-plant environmental covariance, and the total environmental covariance was 3–4 times the genetic covariance. For example, 74·7–75·6 % of the total phenotypic covariance between AH, SH and TL was due to environmental variation within plants (last three trait pairs in Fig. 5), even though the estimated genetic correlations between these traits were very high (Table 2). As a consequence, the total phenotypic correlation reflects primarily the within-plant environmental correlation for most pairs of traits, as can be seen by comparing these two columns in Table 2. The covariances between BB and the three highly-correlated dimensional traits (AH, SH, TL) are apparent exceptions to this pattern, with the genetic and environmental contributions being comparable (the three trait pairs immediately to the left of the dashed box in Fig. 5).

Fig. 5.

Relative magnitudes of the three covariance components for 28 pairs of traits. All traits were standardized to a sample variance of 1·0 before estimating covariance components. The total phenotypic covariance was significantly different from zero for all trait pairs outside of the dashed box.

Floral integration

Measures of genetic, environmental and total phenotypic integration of the hawkmoth floral syndrome were comparable and not significantly different (Table 3). Because of the large heterogeneity in correlation coefficients, these measures of integration are very sensitive to the choice of floral traits to be measured. For example, if only the four dimensional traits had been included then the measures of total phenotypic integration would be $\overline{\left|r\right|}$ = 0·812 and V(λ) = 0·678. The high integration of AH, SH and TL generates low variances in anther position, stigma position and the distance between them, despite high variation in flower sizes (Fig. 3).

Table 3.

Measures of phenotypic integration for the hawkmoth pollination syndrome in R. humilis. Standard errors are not reported for the genetic and environmental components of V(λ) because the computer time required to jackknife the REML variance–covariance estimates was prohibitive

Component	$\overline{\left|r\right|}$ ± s.e.	V(λ) ± s.e.
Genetic	0·359 ± 0·055	0·179
Environmental (within plants)	0·346 ± 0·058	0·209
Total phenotypic	0·342 ± 0·056	0·200 ± 0·019

DISCUSSION

Low canalization of traits in the hawkmoth pollination syndrome

The variability of dimensional floral traits in R. humilis, as measured by the CV, is very close to the average of values reported for a large sample of other species, including both pollinator specialists and generalists (compare Tables 1 and 4). The most variable trait measured in R. humilis was NV, a trait that also tends to be highly variable in other species (Cresswell, 1998). Thus, our data do not suggest a higher than average degree of canalization of floral traits as a result of specialization on a single pollinator guild. However, NC is a possible exception (CV = 0·10). The rate of nectar transport along the hawkmoth proboscis decreases with increasing sugar concentration due to increased viscosity, so there may be an optimal NC that maximizes the rate of caloric uptake (Cruden et al., 1983).

Table 4.

Summary of floral trait variability within populations of animal-pollinated flowering plants, quantified by the CV. Standard deviations were estimated from information provided in the source, if available

Floral traits measured	Species included	Total number of CV estimates	$\overline{CV}$ ± s.d.	Source
Sex organ dimensions	25	54	0·18	Cresswell, 1998
Vector matching dimensions	99	106	0·14	Cresswell, 1998
Advertisement size or shape	41	54	0·22	Cresswell, 1998
Nectar traits	26	58	0·54 ± 0·26	Cresswell, 1998
Petal length	14 zygomorphic	14	0·09 ± 0·02	Wolfe and Krstolic, 1999
Petal length	17 actinomorphic	17	0·12 ± 0·02	Wolfe and Krstolic, 1999
Flower or pseudanthia width	43	43	0·17 ± 0·06	Lázaro and Totland, 2014
Petal length	15 zygomorphic	15	0·11 ± 0·06	Nikkeshi et al., 2015
Petal length	21 actinomorphic	21	0·17 ± 0·08	Nikkeshi et al., 2015

Integration of the hawkmoth pollination syndrome

When the entire set of measured floral traits is taken into consideration, the degree of integration found in CH flowers of R. humilis is comparable to (and slightly less than) the mean of values reported for a large sample of other species, including both pollinator specialists and generalists (Table 5). However, the degree of integration between the four serially homologous dimensional traits that we measured is very high and is well above average when compared with measures of flower size reported in other species (Table 5). These high correlations contrast sharply with the generally low correlations involving other floral traits. Thus, the hawkmoth pollination syndrome as a whole is not highly integrated in R. humilis.

Table 5.

Summary of four reviews of floral integration. Reported statistics are as follows: |$\overline{r}$|, absolute value of the correlation between a pair of traits within a population; $\overline{\left|r\right|}$⁠, average of the absolute values of all pairwise correlations for a set of traits within a population; r_g, genetic correlation between a pair of traits within a population; V(λ), variance of the eigenvalues of the phenotypic correlation matrix for a set of traits within a population, standardized by its theoretical maximum. Standard deviations were estimated from information presented in the sources or obtained from the authors

Floral traits included	Species included	Floral integration				Source
		Statistic	n	Mean ± s.d.
Linear dimensions	28	|$\overline{r}$|	299	0·39 ± 0·23		Conner et al., 2014
Linear dimensions	41 self-compatible	$\overline{\left|r\right|}$	41	0·41 ± 0·16		Fornoni et al., 2016
Linear dimensions	23 self-incompatible	$\overline{\left|r\right|}$	23	0·35 ± 0·12		Fornoni et al., 2016
Sexual and attractive	8 zygomorphic	rg	8	0·51 ± 0·15		Ashman and Majetic, 2006
Sexual and attractive	13 actinomorphic	rg	13	0·23 ± 0·11		Ashman and Majetic, 2006
Not reported	36	V(λ)	80	0·22 ± 0·15		Ordano et al., 2008

If the correlations between floral dimensions reflect nothing more than conserved developmental constraints, then we would predict a higher correlation between the two measures of corolla size (TL and DW) than between either of these and AH or SH, which is contrary to what was observed. Therefore, the very high correlations between TL, AH and SH, compared with their correlations with DW, suggest adaptive integration in response to correlational selection. The likely source of this correlational selection is strong stabilizing selection imposed by pollinators on SA, and on the elevations of these two structures above the floral tube, each of which affects the processes of pollen removal and deposition. Thus, these three traits appear to constitute an adaptive intrafloral module that regulates the deposition and removal of pollen, similar to what has been reported in other species (reviewed by Diggle, 2014). The elevations of the anthers and stigmas above the floral tube are highly canalized against both genetic and environmental variation. The canalization against environmental variation is especially strong since most of the flower size variance is found within plants.

When pollinators exert stabilizing selection on the difference δ between two dimensional traits x and y, selection is imposed on two variance–covariance parameters, which can be seen as follows. Assuming for simplicity that the two dimensional traits have the same variance (⁠${\mathit{\sigma }}_{\mathit{x}}^2$ = ${\mathit{\sigma }}_{\mathit{y}}^2$ = σ²), the variance in δ is given by ${\mathit{\sigma }}_{\mathit{\delta }}^2$ = 2σ²(1−ρ_xy), where ρ_xy is the correlation between the dimensional traits. Thus, selection to reduce ${\mathit{\sigma }}_{\mathit{\delta }}^2$ generates both stabilizing selection on flower size (selection to reduce σ²) and correlational selection on x and y (selection to increase ρ_xy). Because flower size has remained highly variable in R. humilis (σ² is large), low ${\mathit{\sigma }}_{\mathit{\delta }}^2$ has evolved primarily in response to correlational selection, resulting in very high correlations between TL, AH and SH. Furthermore, because ${\mathit{\sigma }}_{\mathit{\delta }}^2$ is directly proportional to σ², the strength of stabilizing selection on δ increases with increasing σ². In other words, tight integration of TL, AH and SH is strongly selected precisely because flower size is so variable.

It is possible that the moderately high correlations between flower size and NV in R. humilis (r = 0·46–0·51) reflect a second intrafloral module that coordinates allocation to pollinator attraction (Bissell and Diggle, 2010; Ellis et al., 2014). For example, there would be correlational selection on these traits if smaller flowers tend to rely on alternative pollinators that harvest pollen instead of nectar. Alternatively, these correlations may simply reflect conserved developmental constraints that cause larger flowers to produce larger nectaries.

Because of the high correlations between the four measures of flower size, the observed negative correlations between flower size and time of BB (r = −0·18 to − 0·28) essentially represent a single phenotypic correlation (between BB and flower size) of low magnitude (r² ≤ 0·08) and thus should be interpreted with caution. It is possible that this is a purely developmental phenomenon, with larger flower buds producing greater quantities of signalling molecules that contribute to the ‘decision’ to initiate BB. Alternatively, an adaptive hypothesis might emerge from data on the foraging schedules and proboscis lengths of alternative native hawkmoth pollinators.

The general lack of integration among floral traits not directly involved in pollen removal and deposition suggests that pollinators have not exerted strong correlational selection on these trait combinations. Although divergent pollinator guilds are expected to promote evolutionary integration across populations or species, there is little reason to expect a single pollinator guild to promote the same pattern of integration within a population or species. For example, within the range of variation expressed in our study population, it is not clear why NV and time of BB would need to be coordinated to ensure successful cross-pollination by hawkmoths. Thus, hypotheses about adaptive integration within populations should be based on models or measurements of selection imposed by specific pollinators, not on unjustified extrapolations from patterns of evolutionary integration.

Genetic versus environmental integration

The estimated genetic and environmental correlations between the four dimensional traits (TL, AH, SH and DW) were positive, high and concordant, suggesting that genetic and environmental perturbations affect these traits through common physiological and developmental pathways, as would be expected for serially homologous traits. Our results suggest that both the genetic and environmental correlations among TL, AH and SH have been further amplified in response to correlational selection.

In contrast, we found no concordance between estimated genetic and environmental correlations for all remaining pairs of traits, i.e. those not comparing serially homologous structures. Estimated environmental variances and covariance were generally much larger than their genetic counterparts, so that phenotypic correlations reflect predominantly the environmental correlations. Thus, when comparisons between serially homologous traits are excluded, there was also no concordance between estimated genetic and phenotypic correlations. Cheverud (1988) suggested that lack of concordance between genetic and phenotypic correlation matrices may reflect the low precision with which the former is estimated, and he argued that the significantly higher variance often observed among estimated genetic correlations, as opposed to phenotypic correlations, supports this hypothesis. However, contrary to Cheverud’s prediction, our estimated genetic and environmental correlations have similar mean magnitudes and variances. Furthermore, Willis et al. (1991) pointed out that Cheverud’s analysis was based almost entirely on morphometric traits that are expected to share developmental pathways, and are therefore expected to show significant correlations for sufficient sample sizes.

Thus, we conclude that the low concordance between estimated genetic and environmental correlations for non-homologous floral traits in R. humilis is likely to reflect low concordance between the actual correlations among these traits. This implies that, for non-homologous traits, the patterns of floral integration depend very much on whether perturbations are genetic or environmental in origin, suggesting that different physiological and/or developmental pathways are being affected. For example, NV and flower size are significantly integrated in the face of environmental variation (r_e ≈ 0·5) but not in the face of genetic variation (r_g ≈ 0).

The relatively low correlations between non-homologous traits, combined with the conflicting patterns of genetic and environmental correlations, suggest an absence of strong correlational selection on these traits, or an inability to respond to such selection. The fact that only serially homologous traits appear to be functionally integrated in this species may be a coincidence, but it is also possible that homologous traits are more amenable to tight integration because their shared developmental architecture results in coordinated responses to genetic and environmental perturbations.

Why is CH flower size so variable in R. humilis?

Hawkmoths that transfer pollen on their bodies have been shown to exert strong stabilizing selection on the length of the floral tube or nectar spur (Nilsson, 1988, 1998; Maad, 2000; Whittall and Hodges, 2007), so the high variability of tube length in R. humilis is surprising. It is possible that inbreeding disrupts developmental homeostasis (Levin, 1970), although this seems unlikely for a species with an evolutionary history of mixed mating and is inconsistent with the observed homeostasis of SA. An alternative possibility is that the physiological switch between CH and CL flower production is not strictly binary, resulting in intermediate flowers. Intermediate flowers are common late in the summer as CL flower production is ramped up, and it is possible that this occurs to a lesser extent throughout the flowering season.

Another possible explanation for variability in flower size and other floral traits in R. humilis is a history of disruptive selection imposed by multiple pollinators. The mixed mating system of R. humilis is consistent with an evolutionary history of unreliable visitation by hawkmoths (Kalisz et al., 2004). Hawkmoth abundance can be highly variable in space and time (Miller, 1981; Martinez del Rio and Burquez, 1986; Pettersson, 1991; Hodges 1995), so secondary visitors may impose substantial selection even if they are less efficient at promoting outcrossing. Although sphingophily is a classic example of pollinator specialization (Darwin, 1862; Fægri and Van der Pijl, 1966; Grant and Grant, 1983), most sphingophilous species have secondary floral visitors that may effect some pollination (Baker, 1961; Gregory, 1963–64; Grant and Grant, 1983; Grant, 1983, 1985, 1993; Pettersson, 1991; Barthell and Knops, 1997; Boyd, 2004; Perez-Barrales et al., 2007; Artz et al., 2010), and competing floral visitors have been shown to generate disruptive selection in several sphingophilous species (Inoue, 1986; Anderson et al., 2010; Kulbaba and Worley, 2012, 2013). In R. humilis, self-pollination promoted by halictid bees may select for reduced SA in the absence of hawkmoths, providing a possible explanation for the incomplete herkogamy. Interestingly, within the Brassicaceae, the highest variation and lowest integration in corolla shape are found in species pollinated by nocturnal moths (Gómez et al., 2016).

Finally, patterns of trait integration and variability in Missouri populations of R. humilis may have a historical component. Palaeobotanical evidence and the absence of endemism suggest the current tallgrass prairie ecosystem is post-glacial (Axelrod, 1985), and it is possible that these fire-climax communities are anthropogenic and thus extremely young (Ladd, 1991; Steuter, 1991). For the long-lived perennials that dominate this ecosystem, the communities are <1000 generations old and are likely to be far from evolutionary equilibrium. Thus, it would not be surprising if many plant–pollinator relations are in evolutionary flux. Rather than promoting visual detection by hawkmoth pollinators, blue floral pigments might be vestigial in R. humilis, reflecting a recent evolutionary transition to hawkmoth pollination. Traits under positive selection, such as lengthened floral tubes, earlier anthesis and reduced NC, are expected to evolve more quickly than vestigial traits under relaxed selection, such as corolla pigmentation and nectar guides. A recent transition to sphingophily might also explain the modest NV in R. humilis, either because the trait has low heritability or because increased NV would allow potential pollinators to access the nectar without contacting the anthers. Furthermore, the geographical distribution of R. humilis may have been split into two glacial refugia during the Pleistocene, with south-western populations adapting to hawkmoth pollinators and south-eastern populations adapting to a different pollinator community. Consistent with this hypothesis, CH flowers from populations east of the Mississippi River have much shorter floral tubes (E. Tripp, University of Colorado, USA, pers. comm.). The post-glacial admixture of populations from different refugia could explain high variability within extant populations (Soltis et al., 1997; Lagercrantz and Ryman, 1990; Cheddadi et al., 2006; Barnard-Kubow et al., 2015). Future work will use molecular markers to address these historical hypotheses.

SUPPLEMENTARY DATA

Supplementary data are available online at https://academic.oup.com/aob and consist of the following. Table S1: SAS code for the REML estimation of variance and covariance components from a diallel crossing design. Table S2: results of GLM models assessing the contributions of individuals and time of day to variance in NV.

ACKNOWLEDGEMENTS

We are grateful to Sarah Dewey, Ben Smith and LeAnn Hubbert for assistance. This work was supported by Missouri State University (Faculty Research Grant F07208-162001-022 awarded to J.S.H.).