This study used the herbaceous plant Common Evening Primrose (Oenothera biennis L., Onagraceae), which is native to open habitats of eastern North America (Dietrich et al., 1997). Plants form a rosette before bolting and flowering, and O. biennis populations frequently genetically vary in life-history strategy, reproducing in their first (annual) or second (biennial) year of growth (Johnson, 2007). Reproduction is fatal, so that plants flower once and then die. Oenothera biennis is dominated by a diverse class of phenolic secondary metabolites (Zinsmeister & Bartl, 1971; Howard & Mabry, 1972) and plants are frequently damaged by herbivores that reduce plant fitness (Johnson & Agrawal, 2005). Oenothera biennis is also functionally asexual, a consequence of its permanent translocation heterozygote (PTH) genetic system, which occurs in ca. 57 species from three plant families (Cleland, 1972; Holsinger & Ellstrand, 1984). The life-history and genetic system of O. biennis make it possible to estimate total male and female lifetime fitness of individual plants and genotypes.
As we discuss below, the functional asexuality of O. biennis is predicted to influence the evolution of plant traits. Because all genes in the genome are effectively linked, the evolutionary consequences of selection on one trait are dependent on the strength and direction of selection on other traits (Otto & Lenormand, 2002). As such, data on the genetic (co)variance and selection on multiple traits can help to identify: 1) traits that are likely to reach an optimum, 2) traits on which interference selection constrains their evolution, and 3) the nature of correlated evolution.
Plant collections and growth
We obtained 39 genotypes of O. biennis from old-fields and roadside populations by collecting ripe fruits from single plants separated by 0.5 km or more (Tompkins County, NY, USA, average distance between populations = 12 km; maximum = 30 km). Delimiting populations of functionally asexual organisms is inherently difficult because individuals do not interbreed, but they may nonetheless interact. In this study we follow previous convention of viewing the genotypes studied as a random sample from a single large population with multiple subpopulations (Johnson 2007). Because of O. biennis’ PTH genetic system, seeds from each collection were genetically identical, and six of the nine microsatellite markers developed for this study were sufficient to distinguish 36 of the 39 genotypes (Larson et al., 2008).
In February 2006, we germinated seeds on moistened filter paper in sealed petri dishes exposed to direct sunlight. Seedlings were transplanted to 500 ml plastic pots filled with potting soil (Metro-Mix, Sun Gro Horticulture, Bellevue, WA, USA) and grown in a glasshouse with supplemental light at Cornell University. Plants were provided ad libitum water and fertilized weekly with a dilute solution of 21:5:20 N:P:K. After 10 weeks of growth, we transferred plants to shaded cages outside for 10 days and then transplanted them into the ground of an abandoned agricultural field. The site was chosen because it lacked naturally occurring O. biennis but had rocky soil characteristic of O. biennis habitat.
We created ten complete experimental blocks, each one having one to two randomized representatives from each of the 39 genotypes (total N = 400, 9-11 replicate plants per genotype). Plants were separated by 1m among rows and columns and experimental blocks were arranged into columns and rows separated by 5m. Plants were irrigated upon planting, but not for the rest of the season. We sprayed five of the blocks (half the plants) with esfenvalerate insecticide (trade name: Asana XL, DuPont Agricultural Products, Wilmington, Delaware) every 2-3 weeks to reduce insect herbivory. To understand whether the insecticide could confound our results by directly influencing the performance of plants, we conducted a six week growth assay prior to our experiment in which plants were grown in the absence of herbivores and treated three times with esfenvalerate or a water control; we did not detect any effects of the treatment on shoot (F1,40 = 0.18, P = 0.28) or root (F1,40 = 0.26, P = 0.61) biomass. In the field, our insecticide treatment reduced the abundance of herbivores (F1,352 = 54.31, P < 0.001). This reduction had no significant effect on lifetime fruit production (F1,38 = 0.24, P = 0.63) or any life-history trait (effect of insecticide: P > 0.30 for all traits; see “Traits measured” below), and there were no significant genotype-by-insecticide interactions (P>0.10). Therefore, we combined the life-history data collected from insecticide and control plots when calculating the genotype breeding values of fitness and life-history traits. To estimate constitutive trait values, we only measured physiology, chemistry and morphology from sprayed plots (N = 200 plants), while herbivory was only measured from unsprayed plots (N = 200 plants).
We measured four types of traits: ecophysiological, life-history, herbivore resistance, and morphological. Ecophysiological traits were defined based on their role in the primary function and metabolism of plants (e.g., C:N ratio, leaf water content). Life-history traits were those most closely related to reproduction, mortality and growth (e.g., flowering time, lifetime biomass). Herbivore resistance traits included concentrations of secondary compounds and direct measures of herbivory. Morphology refers to physical plant traits, of which we only measured trichome density. In late July and early August, we measured traits from all plants in the experiment. Foliar herbivory was almost exclusively due to the exotic Japanese Beetle (Popillia japonica Scarabidae) and was estimated following population outbreak as the proportion of leaves with leaf damage (most plants had well over 30 leaves). This introduced beetle is also a frequent pest of fruit crops in the rose family. In Ithaca, P. japonica infests plants in June and July, and it has been the dominant leaf-chewing herbivore on most local O. biennis populations in recent years (MJ and AA, Pers. Obs.), and common since at least the late-1970s (Kinsman, 1982). Plants were also point censused for beetles on the day that leaf herbivory was measured; despite low counts there was a genetic correlation between our measure of plant level herbivory and beetle abundance (N = 39, r = 0.659, P<0.001).
Two young, fully expanded leaves were taken from every plant for assessment of phenolics (see below) and the ratio of leaf carbon to nitrogen (C:N). Leaves were kept in a cooler on ice, freeze-dried in the laboratory, and ground to a fine powder. C:N ratio was assessed in an elemental combustion system (Cornell University Stable Isotope Laboratory).
To measure leaf water content, specific leaf area (SLA), and trichome density, we took a single leaf disc (28 mm2) from the tip of the youngest fully-expanded leaf on each plant and placed each disk in a sealed plastic tube in a cooler with ice. We weighed each disc to the nearest microgram on the day of collection (fresh mass) and after 24 hours of drying at 40°C (dry mass). We calculated water content as the ratio of water mass (fresh - dry mass) to the fresh mass of the leaf disc. We counted trichomes on both sides of each leaf disc under a dissection microscope and divided this number by disc area to estimate trichome density per cm2. We assessed specific leaf area (SLA) of the disc as leaf area per unit dry mass.
Life-history measures were taken throughout the growing seasons. Plants exhibited one of two life-history strategies: they either remained as a rosette (score = 0) or bolted into a flowering stalk (score = 1) during summer 2006. We measured the number of days from germination (February 1) until bolting (bolting date) and also until the first flowers (flowering date) opened. Surveys in spring 2007 revealed that all plants that remained as rosettes in 2006 died during winter due to mammalian herbivory, pathogens, and other unknown causes. After fruits had finished ripening in November 2006, we harvested all plants, dried the tissue in forced air-drying ovens and weighed the total dry mass. We also counted the number of fruits on each plant, which provided a measure of total male and female fitness.
Leaf sample extraction
We characterized and quantified phenolic composition of O. biennis leaves using 20 mg of leaf powder from each sample, extracted in 500 μl acetone/water (7/3, v/v, containing 0.1% ascorbic acid, m/v) for 1h with a Vortex-Genie 2T mixer (Scientific Industries, N.Y., USA). After centrifugation (14000 rpm, 10 mins; Eppendorf centrifuge 5402, Eppendorf AG, Hamburg, Germany) we stored the supernatant in a separate vial and repeated the extraction on the original tissue three times with fresh solvent (total of four extractions). Acetone was removed from the pooled extracts using an Eppendorf concentrator 5301 (Eppendorf AG, Hamburg, Germany) and the sample was freeze-dried. Prior to analysis, we dissolved the sample in 1400 μl of water and filtered it through a 0.45 μm PTFE filter.
We analyzed phenolics from each sample using high-performance liquid chromatography with a diode array detector (HPLC-DAD). The LaChrom HPLC system (Merck-Hitachi, Tokyo, Japan) consisted of a pump L-7100, a diode array detector L-7455, a programmable autosampler L-7200, and an interface D-7000. Individual phenolics were separated using the Merck Chromolith Performance RP-18e (100 4.6 mm i.d.) column with CH3CN (A) and 0.05M H3PO4 (B) as eluents. After an extensive comparison of solvent and flow rate gradients, we found the following to be most effective. Solvent gradient: 0-2 min, 2% A (isocratic); 2-16 min, 2-30% A in B (linear gradient); 16-18 min, 30-80% A in B (linear gradient); 18-21 min, 80% A (isocratic); 21-23 min, 80-2% A in B (linear gradient); 23-30 min, 2% A (isocratic). Flow-rate gradient: 0-16 min, 1.5 ml/min (isocratic); 16-18 min, 1.5-2.0 ml/min (linear gradient); 18-28 min, 2.0 ml/min (isocratic); 28-28.5 min, 2.0-1.5 ml/min (linear gradient); 28.5-30 min, 1.5 ml/min (isocratic). We used four acquisition wavelengths (200 nm, 280 nm, 315 nm, 349 nm) and UV spectra were recorded for each peak between 195 – 600 nm. For quantitative analyses, 10 µl of the O. biennis extract was injected into the Merck Chromolith Performance column where hydrolyzable tannins (ellagitannins) were quantified in pentagalloyl glucose equivalents (280 nm), flavonoid glycosides in quercetin equivalents (349 nm), and caffeoyl tartaric acid in caffeoyl quinic acid equivalents (315 nm).
Characterization of O. biennis phenolics
O. biennis phenolics separated by HPLC-DAD were first characterized on the basis of their UV spectra. More specific characterization was achieved by HPLC–ESI-MS in the negative ion mode using a Perkin-Elmer Sciex API 365 triple quadrupole mass spectrometer (Sciex, Toronto, Canada). The HPLC system consisted of two Perkin-Elmer Series 200 micro pumps (Perkin-Elmer, Norwalk, CT, USA) connected to a Series 200 autosampler (Perkin-Elmer, Norwalk, CT, USA) and a 785A UV/VIS detector (Perkin-Elmer, Norwalk, CT, USA). ESI-MS conditions were the same as those described in a previous paper (Salminen et al., 1999). Chromatographic conditions were as described above for HPLC-DAD, but the 0.05M H3PO4 was replaced with 0.4% HCOOH. After UV detection at 280 nm, 20% of the eluate was split off and introduced into the ESI-MS system.
Genetic variance, heritability, and coefficients of variation
We used restricted maximum likelihood (REML) in Proc Mixed of SAS (SAS Institute, Cary NC) to estimate the variance explained by plant genotype for each trait. The statistical model included plant genotype and spatial block as random effects, where the significance of genotype was tested using a log-likelihood ratio test (Littell et al., 1996). Because of O. biennis’ functional asexuality, we calculated broad-sense heritability for each trait as H2 = Vg/VT (Lynch & Walsh, 1998), where Vg is the total genetic variance (additive and non-additive) and VT is the total phenotypic variance (genetic and environmental) in the trait. The coefficient of genetic variation was calculated as Vg0.5/μi, where μ is the mean for trait i. All analyses were performed on untransformed data as recommended by Houle (1992).
Genetic associations among traits
Genetic covariation among traits was characterized using Pearson correlation coefficients and genetic covariances. Pearson correlations were determined for all pairwise combinations of traits using the best linear unbiased predictors (BLUPs) (similar to mean values) of genotype breeding values (N = 39 genotypes per correlation). BLUPs are more accurate than family means because they are less biased by environmental effects and more robust to unbalanced replication than are family means (Shaw et al., 1995; Littell et al., 1996). The genetic covariance among traits was calculated according to the equation: covg = rg(G11*G22)0.5, where rg is the genetic Pearson correlation coefficient between two traits, and Gii is the genetic variance of each trait from REML. In the case of life-history strategy (rosette vs. flowering), we used generalized linear mixed models in Proc Glimix to estimate the genetic variance and the binomial equation’s estimate of variance for total phenotypic variance. The statistical significance of genetic covariances was assessed as the P-value from the t-statistic of rg (Lynch & Walsh, 1998; p. 641). To assess the role of multiple tests in producing spurious significant effects, we assessed whether the frequency of significant correlations deviated from the random expectation using the binomial expansion test (Zar, 1996).
To determine how covarying traits were related to one another, we performed hierarchical cluster analysis. This was done by first calculating a square matrix of genetic Pearson correlations coefficients among all traits, followed by implementing Ward’s minimum variance method (Ward, 1963) in Systat (Vers. 9, Systat Software, Chicago IL, USA) to define linkages among traits and groups of traits (Wilkinson, 1999). We objectively defined a “cluster” when all traits were linked by < 2 sums-of-squares variance (i.e. the measure of distance) (Wilkinson, 1999).
Genetic variation in plant traits that predict herbivory
We used two methods to identify how genetic variation in plant traits affected herbivory by P. japonica. First, we performed forward stepwise multiple regression in which we regressed the BLUPs for herbivory against the BLUPs of all plant traits, using Type II (partial) error in Proc Reg of SAS. Stepwise regression was used instead of a fully parameterized model because of the limited statistical power associated with the latter. The model was built by allowing a variable to enter the equation if the linear regression coefficient had a partial P-value <0.10, while we excluded variables that had a partial P-value >0.10. Once the best linear coefficients model was found, we used forward stepwise regression again to explore whether any of the traits already in the model exhibited significant quadratic effects on herbivory. A significant quadratic effect of a plant trait on herbivory indicates that there are non-linear effects of a trait on resistance to herbivores, which could be caused by intermediate levels of a trait conferring the highest or lowest resistance to a herbivore.
Because of the extensive covariation among traits, we also used principal components analysis (PCA) to reduce the dimensionality of trait variation into a smaller number of variables. We performed PCA in Systat using the Varimax rotation method and Pearson correlations as the distance measure among genotypes. We retained all principal components (PC) with an eigenvalue >1 (Legendre & Legendre, 1998), which resulted in nine PCs. We then used the principal components in a forward stepwise regression to determine which components best predicted herbivory.
Natural selection on plant traits
We estimated selection on plant traits using conventional covariance measures of selection differentials (Price, 1970), as well as multivariate genotypic selection analyses (Lande & Arnold, 1983; Rausher, 1992). The selection differential on each trait was measured as the covariance between relative fitness of genotypes and normally standardized trait variation among genotypes.
We then used forward stepwise regression to regress the BLUPs of relative fitness against the BLUPs of all traits including herbivory. This process was started by only including linear coefficients in the model, as these coefficients estimate the strength of directional selection (β). We then included all quadratic coefficients for the traits in the best linear coefficients model to estimate the strength of curvilinear selection, and multiplied these quadratic coefficients by two to estimate γ, the strength of quadratic selection (Lande & Arnold, 1983; Stinchcombe et al., 2008). We visually inspected the partial regression plots of relative fitness versus trait variation, and used the equation for computing maximum/minimum values (i.e. -β/γ) to determine whether there was stabilizing or disruptive selection on traits (Mitchell-Olds & Shaw, 1987). We inferred that stabilizing selection acted on a trait when the quadratic coefficient was negative and exhibited a maximum value within the range of data, and disruptive selection when the quadratic coefficient was positive and showed a minimum value within the range of data. Finally, principal components analysis was used again, here including herbivory, to reduce the dimensionality of the data; nine PCs had eigenvalues >1. As before, we used stepwise multiple regression to determine how relative fitness related to the PCs.
Characterization of phenolics
Eleven HPLC peaks were uniformly distributed across all plants, which included six ellagitannins, four flavonoids and caffeoyl tartaric acid (Fig. 1 a,b). Other peaks were present in only some plants and/or were poorly chromatographically separated (e.g., peaks at retention times 8.5 – 11.0 min). These latter peaks were characterised as either ellagitannins or flavonoid glycosides based on their UV spectra and included into the subgroups “other ellagitannins” and “other flavonoid glycosides”. Detailed descriptions of the chemical profiling by HPLC-DAD and HPLC-ESI-MS will be presented in a forthcoming article (Salminen et al., unpublished results).