Larval performance assay

Caterpillars used for the performance assay were obtained from eggs laid in July 2013 by females taken from Silver Lake, NV (SLA, Fig 1). These females utilize Astragalus canadensis (Fabaceae). Alfalfa does not occur at Silver Lake, and because L. melissa tend to be very localized in their movements [45], the females collected were unlikely to have previously encountered alfalfa. Females were placed in oviposition arenas (described below) containing only A. canadensis, and eggs were collected after two days and pooled. Within hours of hatching, caterpillars were placed singly in 20 x 90 mm petri dishes containing alfalfa from one of the seven source populations (Fig 1). Caterpillars were kept under lamps at room temperature (approximately 20°–23°C) and allowed to eat ad libitum, with leaves replaced at the first sign of wilting. Approximately sixty larvae were reared on each alfalfa population (n = 61 to 63, depending on the population).

Given the geographic distances separating alfalfa locations, it was not logistically feasible to keep all caterpillars supplied with sufficient fresh foliage to complete development. Therefore, most caterpillars were sacrificed at 14 days (after recording mass and survival), approximately halfway between egg hatch and adult eclosion. Caterpillars being fed alfalfa from three locations (VUH, AWFS, and APPL) were reared to adults, then sexed and weighed. These three locations were chosen because two of them (VUH and AWFS) are close to each other (~8.5 km apart) and close to the research laboratory at the University of Nevada, Reno, and initial results (after 14 days) from the third (APPL; which was extralimital to the range of L. melissa) suggested that it was potentially an informative, extreme contrast to the other two populations. All mass measurements were taken with a Mettler-Toledo XP26 microbalance to the nearest microgram.

Oviposition preference assay

Oviposition preference of adult female butterflies was assayed as per [35], with females taken from either Silver Lake, NV (SLA) or Verdi, NV (VUH). As mentioned above, females from SLA are very unlikely to have encountered M. sativa prior to capture. Female L. melissa from VUH use alfalfa, and no known native host plants occur at this location. Females were placed in oviposition arenas containing alfalfa foliage from two out of the three sources that were used in the full-development rearing described above; specifically, either AWFS and APLL, or AWFS and VUH were paired in each oviposition arena. As a negative control, leaves from Lotus nevadensis were also included in each arena. L. nevadensis is a Fabaceous plant that L. melissa does not consume in the wild, but can subsist upon [34] L. nevadensis was collected near Yuba Gap, CA. These choice tests using wild-caught females are efficient and effective: they provide results that are similar to no-choice tests, and similar to results from lab-reared females [31].

Tests were conducted using plastic cups (12 cm x 9.5 cm) as oviposition arenas, with each cup containing three different branches (one branch each from the two alfalfa populations being compared, and the negative control) that the female could choose between. Branch stems protruded from small holes in the bottom of each cup. A second cup was filled with water and placed underneath the branch-containing cup, with stems submerged in water. Female butterflies were sealed inside these cups with fine mesh, and mesh was misted every two to four hours with tap water from 10 a.m. to 5 p.m. and swabbed with fruit punch Gatorade^® twice daily as a nectar substitute [35]. Oviposition arenas were kept outside in dappled shade on the University of Nevada, Reno campus at ambient temperature (approximately 20°–30°C). Assays were conducted from August 19, 2013 through September 5, 2013. After 48 hours within the arena, females were removed and the number of eggs laid on each plant was counted.

Analyses of preference and performance data

Preference data were analyzed in a hierarchical Bayesian framework using the bayesPref package [46] in the R computing environment [47]. In contrast to a frequentist analysis of preference data that provides information regarding the rejection of a null of no difference in preference, this approach estimates preference for each of the different hosts (along with 95% equal-tail probability intervals [ETPIs]). Models employ a Markov chain Monte Carlo (MCMC) approach to characterize posterior probability distributions describing female preference for a plant at two hierarchies, that of the individual and the population (from which the female was drawn). For a full description of the form of the likelihood function and conditional priors for individual preference see [46]. Models were run for 20,000 MCMC iterations, with a 5,000 iteration burn-in, and output was examined to ensure adequate mixing. Larval survival was also analyzed with this methodology using counts of surviving individuals on different alfalfa sources as data (5,000 iterations with 1,000 iteration burn-in). This allowed us to estimate the probability of survival on a particular alfalfa source (in the same way preference for a particular alfalfa source is modeled with counts of eggs from preference experiments). Means of samples characterizing posterior distributions were used as point estimates for population level preference or survival. ETPIs of posterior probability distributions were examined for overlap in a pairwise fashion to determine differences in preference and survival among alfalfa populations. If ETPIs did not overlap for a given pairwise comparison, then preference or survival was inferred to differ between alfalfa populations. Finally, a Bayesian analysis of variance was used to determine differences in larval weight between populations [48]. Deflections from the mean for each population were sampled from normal distributions with a mean of zero and precisions were independently modeled for each population from folded t distributions (μ = 0, τ = 0.001). The grand mean was modeled from a normal distribution centered at zero with precision 0.001. Within-group variation was modeled independently to account for non-homogeneous variances among groups by sampling from a gamma distribution with shape and rate parameters of one. The model was run for 1,000,000 MCMC iterations divided between three chains with a burn-in of 5,000 iterations and a thinning rate of 1/100. ETPIs of posterior probability distributions of mean larval mass estimates were examined for overlap among populations as described above.

Phytochemical analysis

Foliar tissue from 19–20 individuals from each alfalfa population was collected in late August 2013 (while larval performance experiments were ongoing), stored in paper bags, and frozen at -20°C until extraction. Plants selected were mature and collected haphazardly from throughout each population; the newest growth was avoided in phytochemical analyses. Approximately 150 mg of foliar tissue was extracted twice using a 70% meOH solution followed by 20 minutes of sonication. Samples were dried under reduced pressure, resuspended in 1 ml meOH, and passed through a 0.45 μm filter into an autosampler vial for high-performance liquid chromatography (HPLC) analysis. HPLC analyses were performed using a Waters Alliance HPLC system with a 2996 diode array detector and Empower Pro Software. Each injection was 10 μl eluted on a gradient (90:10 (water:acetonitrile), reaching 60:40 at 40 minutes and 5:95 at 60 minutes) at the rate of 1 ml/min on a Symmetry C-18 reverse phase column (3.5 μm, 4.6 x 75 mm) (Waters Corp.). Phytochemical variation was characterized by retention time and UV absorbance between 230 and 400 nm. This restricted our characterization of compounds to those recoverable by our extraction and HPLC protocol, and those with chromophores for UV absorption (usually those with double bonds), thereby providing a fingerprint of plant phytochemistry based on anonymous compounds. Many saponins (which are compounds known to occur within M. sativa) do not possess chromophores and consequently our method is not sensitive to these compounds. However, there are some saponins detectable via UV spectroscopy; additionally, our method should recover members of many other compound classes reported from M. sativa including flavonoids [49] and phytoestrogens [50]. Data were standardized by dry weight and Hellinger-transformed. This transformation is recommended prior to ordination of datasets containing many zeros, as was the case for our data [51].

Spectral data for all samples were ordinated using non-metric multidimensional scaling (NMDS) on a Manhattan distance matrix created using the vegan package in R (version 2.2 [52]). Examination of stress scree plots showed that an ordination across five dimensions provided the best compromise between stress reduction and complexity (with stress at 10.5). Linear discriminant analysis (LDA) was also employed using the MASS package (version 7.3 [53]) in R to explore phytochemical differences among alfalfa populations. LDA seeks to find the linear combination of predictor variables that best separate data by group, based on a pre-defined set of groups. We employed this approach to determine if phytochemical differences among alfalfa populations could predict L. melissa colonization of those populations. Also, we tested if observed phytochemical differences among populations could predict concomitant variation in larval performance.

To test if phytochemistry could predict L. melissa occupancy of an alfalfa population we first constructed a model trained using a randomly selected subset of the phytochemistry data (n = 90, for each iteration), and examined how well this model could correctly assign colonization status to validation data (n = 45, for each iteration). This process was repeated 10,000 times and the proportion of correct assignments was extracted at each iteration and tabulated. The mean of this vector of results was used as an estimate of model success. For this analysis we only used data for the 28 compounds that occurred across all alfalfa populations, as linear separation often occurred when using the whole dataset that included many rare compounds. We used Monte Carlo simulations to confirm that the results of the LDA would not be expected given a random distribution of chemotypes among populations, or a random distribution of compound concentrations among individuals (see S1 Appendix for details [54]).

Exploring phytochemical diversity

Phytochemical variation within and between alfalfa populations was also examined using the complexity-as-diversity approach as described by [55]. This approach extends the use of numbers equivalents [56, 57] to phenotypic complexity using (Eq. 1 in [57]) where the diversity (^qD) represents the effective number of distinct phenotypes in a group, the diversity order q influences how sensitive D is to rare compounds, and p is the proportion of the data composed of the i-th compound (out of s total compounds). When q = 0, all compounds are weighted equally (i.e., chemical richness). When q = 1, compounds are weighted by their relative abundance. As q increases, the relative importance of rare compounds is reduced. Diversity profile plots across orders of q were used to assess the importance of rare and common compounds in determining the number of distinct chemical phenotypes (see [55]). This approach was taken at three nested hierarchical levels: plant, population, and global. At the within-population level, β-diversity is the effective number of distinct chemical phenotypes among plants (chemotypes), whereas at the among-population level, β-diversity is the effective number of distinct chemical phenotypes among populations. Diversity equivalencies were calculated using the vegetarian and hierDiversity packages in R [58, 59].

Finally, phytochemical covariance matrices were constructed using data specific to each population. These matrices were correlated with site-specific genetic covariance matrices (described below) and the significance of correlations was tested with Mantel tests (1,000 permutations; vegan package in R) to determine the extent to which differences between plant genotypes could explain phytochemical differences between those same individuals.

Foliar tissue total protein determination

A Bradford assay [60] was used to ascertain the protein content of foliar tissue taken from the plants used in the phytochemistry assay (the same 19–20 individuals per population). See supplemental methods for a full description of protocol used (S1 Appendix). Prior to statistical analyses, absorbency data were standardized by sample mass. Site-specific distance matrices (Euclidean) of absorbencies were correlated with genetic and phytochemical distance matrices to test for relationships between these three axes of host variation at each alfalfa population (significance of correlations was determined using Mantel tests).

Alfalfa population genetics

DNA was isolated and purified from desiccated (oven dried) leaf tissue sampled from 132 alfalfa plants using Qiagen's DNAeasy 96 Plant Kit (Qiagen Inc.). DNA was taken from the plants used in the phytochemistry and protein assays described above, however, DNA was only taken from individuals belonging to five of our seven focal populations. DNA was not successfully extracted from APLL and BWP samples; insufficient yield in these cases likely resulted from compromised DNA quality associated with drying. We generated DNA fragment libraries for genotyping-by-sequencing using our established protocol [61, 62, 63]. For details of our library preparation, sequencing, and bioinformatics protocol see supplemental methods. Briefly, genomic DNA was first enzymatically digested with the restriction enzymes EcoRI and MseI and double-stranded adaptor oligonucleotides ligated onto the digested DNA fragments. These adaptors included 8–10 base pair (bp) barcode sequences that were used to match sequences to individual plants. Fragment libraries were amplified using PCR and size-selected to fragments 250 and 350 bps in length using a BluePippin (Sage Science). Libraries were sequenced on an Illumina HiSeq 2500 (one lane, 1 x 100 bp reads) at the University of Texas Genome Sequence and Analysis Facility.

Sequences were aligned (236 million DNA sequences total) to a draft genome that we generated from the diploid progenitor of alfalfa (total scaffold length = 673 Mbp, N50 scaffold size = 37 kbp, number of scaffolds = 41319; we will more fully describe this genome sequence in a future publication, [64]). We then used the Unified Genotyper in GATK [65] to identify variable nucleotide positions and calculate genotype likelihoods for each individual and variable position (this is a Bayesian genotype and variant caller). We assumed a ploidy of four as alfalfa is a tetraploid. We set the minimum base quality to 20 and set the prior expectation for heterozygosity to 0.001. A custom Perl script was then used to filter the initial set of variants to those that met our quality criteria (see supplemental methods). Seventy-one plants (five populations) and 16,920 single nucleotide variants (SNVs) were retained for population genetic analysis.

We estimated the posterior probabilities of each genotype for each individual at each locus using a Bayesian approach. To do this we took the genotype likelihoods from GATK and multiplied them by the prior probabilities of each genotype assuming Hardy-Weinberg genotype frequencies and a maximum likelihood estimate of global non-reference allele frequency, which was also obtained using GATK. We then took the mean of the posterior distribution for each locus and individual as the genotype estimate for downstream analysis (this value is between zero and four, is not constrained to be an integer, and is an estimate of the number of non-reference allele copies at a locus). This general approach, which has been used previously [33, 63] allowed us to make better use of the information in low to moderate coverage population genomic data than if we had simply called genotypes directly from the sequence data for each individual.

We used several methods to quantify and summarize patterns of genetic variation within and among populations. We generated a genetic covariance matrix where each element in the matrix measured the genetic similarity (covariance in genotypes) for a pair of individuals. We then used principal components analysis (PCA) to visualize patterns of genetic similarity based on this matrix. Next, to quantify genetic variation within populations, we calculated the average genetic variance (1 − [p² + (1 − p)²]) and the variance in PC 1 and 2 scores. Genome-average pairwise and global Fst values were then estimated from the sample allele frequencies to assess the extent of population genetic structure. Finally, we tested for a positive correlation between genome-average Fst and the geographic distance between pairs of populations to determine whether alfalfa showed signs of isolation-by-distance (a Mantel test with 1,000 permutations was used to test whether the observed correlation was significantly different from zero at α = 0.05).

Larval performance

A total of 434 caterpillars were fed alfalfa from seven different alfalfa source locations. Survival in the first half of development was quite variable across populations, ranging from less than 50% to greater than 80% (Fig 2a). Interestingly, the alfalfa source (APLL) conferring the lowest larval survival in the first 14 days supported the highest survival (among a subset of three alfalfa sources) across the entire course of development (Fig 2b). Mass was also variable among treatment groups (Fig 2c), with individuals consuming alfalfa rom two of the sources (AFAL and APLL) being two or more times greater in mass than the individuals reared on the other alfalfa sources. The consumption of APLL alfalfa was associated with the largest butterflies at the end of development (Fig 2d), and fastest time to eclosion as compared with VUH and AWFS alfalfa (Wilcoxon rank sum test, p < 0.01) (S1 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Results from caterpillar performance trials: survival (a & b), and mass (c & d).

Performance in (a) and (c) are across all alfalfa source populations through 14 days of development (>60 larvae reared on each alfalfa population); performance in (b) and (d) involve caterpillars reared to adults on a subset of alfalfa sources (number of surviving, weighed larvae in parentheses). Bars in (a) and (b) are 95% equal-tail probability intervals (ETPIs) of survival rate estimates; bars in (c) and (d) are 95% ETPIs of estimates of mean mass. AWFS in panel (d) was not included in analyses because we only have mass from a single individual (no other larvae survived on alfalfa from this population).

https://doi.org/10.1371/journal.pone.0147971.g002

Oviposition preference

We challenged females presumed naïve to alfalfa (females were from SLA, where alfalfa is not known to occur) and females associated in the wild with alfalfa (from VUH) with plants from either VUH and AWFS, or from AWFS and APLL, the three alfalfa sources used in extended larval rearing (see S2 Table for the number of females assayed and other details). Females from SLA showed greater preference for AWFS alfalfa in both trials, compared to APLL and VUH alfalfa and the negative control (Fig 3a and 3b). The behavior of alfalfa-associated females (from VUH) was more complex: they discriminated between AWFS and APLL (Fig 3c), but then preferred alfalfa from their home location when given a choice between VUH and AWFS (Fig 3d). In one of the experiments (VUH females choosing between AWFS and APLL; Fig 3c), the negative control received more eggs than one of the alfalfa sources (APLL). This willingness to lay eggs on the negative control is consistent with behavior previously observed in females from this area [35].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Preference of Silver Lake (a & b) and Verdi females (c & d) for different combinations of alfalfa sources (gray symbols) and a negative control (open symbol).

Silver Lake (SLA) females were from a native-host population (no previous exposure to alfalfa), and Verdi females were from an alfalfa-feeding population (VUH). Bars denote 95% equal-tail probability intervals for population preference and symbols indicate the mean of posterior probability distribution. Different letters denote non-overlapping probability intervals.

https://doi.org/10.1371/journal.pone.0147971.g003

Protein content

Variation in total protein content was observed among surveyed alfalfa populations (Kruskal-Wallis test, χ² = 17.7, df = 6, p < 0.01; a rank-based test was used because of non-normality of data) (Fig 4a). Plants from VUH and AFAL had significantly higher protein content compared with those from GVL (post-hoc multiple comparison test after Kruskal-Wallis, p ≤ 0.05). We failed to detect a significant correlation between protein content at the population level and larval weight gain or survival: Spearman’s rank correlation, rho = 0.61, p = 0.16, and rho = -0.75, p = 0.07, for weight gain and survival respectively (Fig 4b and 4c). Neither did we detect significant correlations between distance matrices of protein content, phytochemistry, or genetic covariance (S3 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Results from Bradford assay of foliar protein in seven populations of alfalfa.

(a) Absorbance (standardized by mass) by population, with significant differences denoted via superscripts (Kruskal-Wallis). Absorbance is directly proportional to protein content. (b & c) Relationships between protein concentration and larval survival and mass gain at 14 days were not significant. Lines denote standard errors of mean estimates.

https://doi.org/10.1371/journal.pone.0147971.g004

Phytochemistry

Our HPLC protocol resulted in data for 49 compounds. Of these, 28 compounds occurred in all alfalfa populations, five occurred only in populations colonized by L. melissa, and two occurred only in uncolonized populations. Consequently, phytochemical α-diversity was very similar across populations (S2 Fig). Calculation of phytochemical β-diversity across all populations suggested few distinct chemical phenotypes were present (~1.2 chemotypes; S2 Fig). However, when grouping populations in terms of L. melissa colonization, we found that colonized populations exhibited less among-population chemical heterogeneity (Fig 5a), but had similar levels of within-population phytochemical heterogeneity, which also tended to be higher than in uncolonized populations (Fig 5b). AFAL was an exception to this pattern, as it exhibited high within-population phytochemical variation, but was uncolonized by L. melissa (Fig 5b).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Phytochemistry and genetic structure of surveyed alfalfa populations.

(a) Number of distinct chemical phenotypes for populations colonized by L. melissa (circles), and uncolonized populations (triangles). Higher values of q reflect higher order diversity equivalents. Abundant compounds are more heavily weighted at values of q > 1. Underlying phytochemical data was collected via HPLC. (b) Uncolonized populations were phytochemically dissimilar, but typically had low within-population phytochemical variation. On the other hand, colonized populations were phytochemically similar, despite the fact that these populations tended to exhibit more within-population phytochemical variation than uncolonized populations (also see S2 Fig). (c) NMDS of phytochemical data showed overall similarity in phytochemistry across surveyed alfalfa populations. (d) However, LDA showed that important phytochemical variation did exist, and could well separate populations differing in colonization by L. melissa. Each point represents an individual plant. Values plotted were calculated by substituting each datum into the discriminant function. The x-axis is not labeled because it serves only to spread samples for visualization. (e) PCA of genetic covariance matrix constructed using ~17,000 SNVs.

https://doi.org/10.1371/journal.pone.0147971.g005

NMDS ordination suggested considerable overlap in phytochemistry among populations (Fig 5c, S3 Fig). However, LDA demonstrated that important phytochemical variation among populations did exist. Indeed, phytochemical differences between colonized and uncolonized populations allowed for near complete separation by LDA (Fig 5d). This result held both when analyzing data for all 49 compounds, and when analysis was limited to those compounds occurring at all alfalfa populations (28 compounds) (S4 Fig). When limiting the analysis to compounds occurring at all populations, the mean estimate of the ability of LDA to successfully predict colonization status was 67% (95% confidence intervals: 53%–80%). The predictive ability of the LDA was confirmed by comparison to output of simulated null models (see S1 Appendix). LDA also demonstrated differences in phytochemistry among all seven populations (S5 Fig). For this analysis, the top four functions output by LDA provided 28.4%, 21.2%, 18.1%, and 15.7% separation respectively (S5 Fig). However, while providing good separation of alfalfa populations in terms of phytochemistry, these discriminant functions did not predict larval performance (S5 Fig). Finally, we found little correlation between our measures of genetic and phytochemical variation among individual plants. Genetic covariance and phytochemical covariance were not correlated at any population, save for a weakly-positive signal at GVL (Mantel tests, S4 Table).

Population structure

A total of 16,920 single nucleotide variants were identified and used to construct a genetic covariance matrix, which was analyzed via PCA. PCs 1 and 2 explained 16.2% and 5.3% of the variation in pairwise genetic similarities among individuals (Fig 5e). VUH and GVL differed most with respect to PC1, whereas PC2 separated these populations from AWFS. The AWFS population was notably more variable with respect to both PC1 (var = 0.011) and PC2 (var = 0.0052) scores (var PC1: AFAL = 0.003, GVL = 0.006, SCC = 0.005, VUH = 0.008; var PC2: AFAL < 0.0001, GVL = 0.0008, SCC = 0.0004, VUH = 0.0011). All five populations exhibited similar levels of genetic variance (0.27–0.29). Overall, little genetic variation was partitioned among populations (global Fst = 0.020, 95% bootstrap CIs 0.020–0.021). However, population pairs varied in the degree to which they differed from each other (mean pairwise Fst = 0.013, minimum = 0.008, maximum = 0.019) (S5 Table). Finally, we found no evidence that genetic differences were greater between more geographically distant populations (Mantel r = -0.463, p = 0.95).