Abstract
Adaptive evolution requires both raw genetic material and an accessible path of high fitness from one fitness peak to another. In this study, we used an introgression line (IL) population to map quantitative trait loci (QTL) for leaf traits thought to be associated with adaptation to precipitation in wild tomatoes (Solanum sect. Lycopersicon; Solanaceae). A QTL sign test showed that several traits likely evolved under directional natural selection. Leaf traits correlated across species do not share a common genetic basis, consistent with a scenario in which selection maintains trait covariation unconstrained by pleiotropy or linkage disequilibrium. Two large effect QTL for stomatal distribution colocalized with key genes in the stomatal development pathway, suggesting promising candidates for the molecular bases of adaptation in these species. Furthermore, macroevolutionary transitions between vastly different stomatal distributions may not be constrained when such large-effect mutations are available. Finally, genetic correlations between stomatal traits measured in this study and data on carbon isotope discrimination from the same ILs support a functional hypothesis that the distribution of stomata affects the resistance to CO2 diffusion inside the leaf, a trait implicated in climatic adaptation in wild tomatoes. Along with evidence from previous comparative and experimental studies, this analysis indicates that leaf traits are an important component of climatic niche adaptation in wild tomatoes and demonstrates that some trait transitions between species could have involved few, large-effect genetic changes, allowing rapid responses to new environmental conditions.
THE course of phenotypic evolution depends upon both the genetic basis of functionally important traits and the mechanistic relationship between traits. All else being equal, traits governed by mutations of large phenotypic effect with low pleiotropy can respond to selection more readily than traits that are complex, highly pleiotropic, and/or strongly correlated with other traits. Hence, the genetic basis of trait differences and covariation influences which phenotypic axes are most likely to respond to selection in new environmental conditions; understanding this underlying architecture therefore provides insight into the genetic mechanisms constraining or facilitating phenotypic evolution (Lee et al. 2014). It can also indicate which evolutionary forces are responsible for trait differences within and between species. When adaptation to a complex ecological niche involves many traits, for example, genetic analysis can determine whether correlations between traits are caused by a shared genetic basis (pleiotropy) or whether natural selection favors trait combinations because they function well together (coselection). The magnitude and direction of allelic effects between phenotypically divergent genotypes can also be used to infer whether this divergence is consistent with directional natural selection (Orr 1998b; Rieseberg et al. 2002), and thus can provide evidence for adaptation, independent of inference from other approaches, such as comparative methods. The genetic architecture of putatively adaptive traits can therefore reveal why some traits evolve rapidly and others more gradually over micro- and macroevolutionary timescales.
With these goals, the genetics of adaptive traits have been successfully examined in a number of cases, most notably of qualitative transitions in pigmentation in response to variation in habitat coloration (e.g., Nachman et al. 2003; Hoekstra et al. 2006) or pollinators (e.g., Schemske and Bradshaw 1999; Streisfeld et al. 2013). In comparison, the genetics of adaptive structural (morphological) traits are less well understood, but have been described in the context of, for example, body armature (skeletal plates) in sticklebacks (Colosimo et al. 2005) and leaf shape evolution in Solanum (Kimura et al. 2008). Despite these important advances, we still know little about the genetic basis and evolutionary forces underlying physiological and morphological evolution to complex environmental variation (Des Marais et al. 2014).
Here we examine the genetic architecture, evolutionary forces, and physiological mechanisms that contribute to divergence in leaf traits connected to adaptive differentiation. Leaf traits resemble other functionally important quantitative traits in that they face numerous selective pressures and functional constraints. As the primary photosynthetic organs in terrestrial plants, leaves mediate several fundamental fitness tradeoffs, including water lost through transpiration and carbon (CO2) gained through photosynthesis. Quantitative traits such as leaf size, shape, and anatomy affect these important physiological functions (Table 1). However, although these traits represent important axes of phenotypic variation across all land plants, respond to natural selection, and appear to vary adaptively in nature (Smith and Hughes 2009; Nicotra et al. 2011; Vogel 2012), relatively little is known about their genetic architecture. Our first goal in this study was to understand the genetic basis of variation in a set of these leaf traits. Specifically, we map QTL underlying interspecific divergence between two closely related tomatoes, the domesticated tomato (Solanum lycopersicum var. esculentum) and a wild relative, S. pennellii. Our focal traits were chosen based on a priori functional considerations, their observed associations with climatic variation, and/or their known ecological significance in wild tomatoes or more generally (Table 1).
Our second goal was to use these data to evaluate evidence for the role of natural selection in shaping leaf trait variation between Solanum species. By revealing the direction of allelic effects underlying functionally important traits, QTL analysis can provide evidence that traits evolved under directional natural selection rather than genetic drift. In particular, Orr (1998b) showed that it is unlikely that many QTL affecting a trait in the same direction (e.g., smaller leaf size) will fix in the absence of directional selection. His proposed statistical test of phenotypic selection—the QTL sign test (QTLST)—has since been used to show that the distribution of QTL effects for many trait differences between taxa are consistent with the action of selection (Rieseberg et al. 2002). We applied this test to assess the evidence that natural selection was involved in fixing leaf trait differences between tomato species. We also used our genetic analysis to assess whether observed trait correlations are more likely due to genetic associations or to the action of natural selection. Several of our leaf traits covary among wild tomato species: species from drier habitats generally have smaller leaf size, lower leaf mass per area, and greater stomatal ratio (more even distribution of stomata on upper and lower leaf surfaces) (C. D. Muir, unpublished data). This covariation could be the product of coselection, whereby each trait improves thermoregulation and photosynthetic efficiency in dry environments. Alternatively, trait covariation across species could be due to a shared genetic basis; traits are constrained to evolve in concert because of underlying molecular/developmental pleiotropy. Thus, we also aimed to determine whether genetic correlations caused by pleiotropy or linkage recapitulate, and therefore could explain, phenotypic correlations across species. One caveat is that increased linkage itself can be the product of selection to keep locally adapted alleles together in disequilibrium (Noor et al. 2001; Rieseberg 2001; Kirkpatrick and Barton 2006), as suggested in some species (Lowry and Willis 2010; Jones et al. 2012). However, if traits are not genetically correlated, this supports an inference that they covary across species because of natural selection.
Of the leaf traits examined in this study, stomatal density and stomatal ratio (Table 1) show the strongest evidence for the action of natural selection (Results). Our third goal was to more directly dissect the mechanistic significance of variation in these stomatal traits by evaluating their association with variation in carbon isotope discrimination (δ13C), a widely used indicator of integrated water-use efficiency (WUE) in quantitative genetic studies (McKay et al. 2003, 2008; Juenger et al. 2005). Both traits potentially affect water-use efficiency and carbon isotope discrimination; however, because each trait influences different steps in the CO2 diffusion pathway (atmosphere → inside of leaf vs. inside of leaf → chloroplasts; see Materials and Methods and Discussion), the correlation between these different stomatal traits and δ13C can be used to infer which portion of this diffusion pathway is more functionally significant in these species. We aimed to use associations between these structural traits and integrated water-use efficiency (previously measured on the same lines) (Xu et al. 2008) to more precisely pinpoint the physiological significance of putatively adaptive leaf trait variation in these species.
Our final goal was to identify candidate genes involved in leaf and stomatal development that colocalize with QTL for our stomatal traits. Identifying specific candidates underlying quantitative trait loci generates targets for future functional analysis (Dalziel et al. 2009) but can also potentially indicate specific developmental and physiological mechanisms underlying observed trait variation. This is especially fruitful for traits associated with a relatively discrete and well-characterized developmental network, such as for stomatal development (Pillitteri and Torii 2012). Patterns of molecular evolution in candidate genes, especially those consistent with positive selection, can also provide an additional method for narrowing specific functional candidates for specific QTL.
To address these four goals, we measured leaf phenotypes (Table 1) on an extensively examined set of introgression lines (ILs) consisting of homozygous introgressions from donor parent S. pennellii (Spen) in the genetic background of domestic tomato, S. lycopersicum var. esculentum (Slyc). Previous studies using this IL population have identified numerous QTL (Lippman et al. 2007) for traits ranging from fruit metabolite concentrations (Schauer et al. 2006), yield (Semel et al. 2006), reproductive isolation (Moyle and Nakazato 2010), and leaf traits (Holtan and Hake 2003; Chitwood et al. 2013). With these data, we confirm evidence in support of adaptive evolution of leaf traits, examine the physiological and ecological significance of trait variation in several clearly adaptive traits, and identify potential molecular developmental mechanisms responsible for this trait variation. In doing so, we demonstrate the power of genetic analyses of adaptation, especially in combination with detailed physiological, ecological, and genomic data, to address the likely evolutionary dynamics of adaptive trait evolution.
Materials and Methods
Plant material and cultivation
For our QTL analysis, we used the previously developed introgression line library (Eshed and Zamir 1994), consisting of short known segments of the wild tomato species, S. pennellii (donor) serially introgressed into the genetic background of the domesticated tomato, S. lycopersicum var. esculentum (recipient). Donor introgressions are homozygous in a homozygous recipient background. Development of these ILs has been described elsewhere (Eshed and Zamir 1994). Briefly, all lines are the advanced-generation backcross progeny of a single F1 plant produced by crossing Spen as the male parent to Slyc. Original backcross lines were genotyped using 350 RFLP markers (Eshed and Zamir 1994), several subsequent analyses increased marker density and refined the location of introgression breakpoints (Lippman et al. 2007), and these were recently confirmed using whole-genome resequencing (Chitwood et al. 2013). We grew 70 of the 76 introgression lines (ILs) as well as the recipient and donor species from seeds obtained from the Tomato Genetics Resource Center (TGRC) at University of California Davis (http://tgrc.ucdavis.edu). Between January 9 and January 13, 2010, seeds were soaked in 50% household bleach for 30 min, rinsed thoroughly, and placed on moist paper in plastic boxes to germinate in a growth chamber. Half of the seeds were cold treated as part of a separate experiment on germination response. We did not attempt to characterize the effect of seed cold treatment on adult plant traits, which were minor; however, we describe the treatment for completeness. For the cold treatment, seeds were placed on wet filter paper and placed in a cold room (−4°) in the dark to mimic winter for 3 weeks. Control seeds were placed on dry filter paper, stored in the dark, and hydrated 3 weeks later. We included seed treatment (cold vs. control) in statistical analyses to account for any plasticity caused by this manipulation. Following germination, all seedlings were transplanted to 7.56-liter pots, and grown on benches in a fully randomized design at the Indiana University greenhouse with supplemental lighting to maintain a constant 16:8 h light:dark cycle. Plants were irrigated to field capacity daily to prevent drought stress and fertilized weekly.
Trait measurements
We measured three leaf morphological traits: leaf area (LA), leaf mass per area (LMA), and leaflet dry matter content (LDMC) and four leaf epidermal traits: stomatal density (SD), stomatal ratio (SR), trichome density (TD), and trichome ratio (TR). LA was measured on all available plants (Slyc sample size nSL = 51, Spen sample size nSP = 18, and mean IL sample size 
), whereas other traits were measured on a subset of these. Sample sizes for LMA and LDMC were: nSL = 50, nSP = 18, and 
. Sample sizes for stomatal and trichome traits were: nSL = 48, nSP = 7, and 
. For all traits, two leaves per plant were collected at two developmental time points: “early” (63–65 days after germination) and “late” (77–82 days after germination). Leaves were marked early in their development (<2 cm in length) and measured only after they had fully expanded. The mean leaf age at time of measurement was 33.1 days (range 24–46 days) and 40.3 days (range 21–69 days) for early and late leaves, respectively. LMA increased slightly with leaf age; therefore, we accounted for leaf age statistically while estimating genotypic values. Leaf morphological traits LA, LMA, and LDMC were measured as described in Muir and Moyle (2009). Briefly, LDMC is the dry mass divided by the fresh mass; LMA is the dry mass divided by the projected area of fresh leaf tissue. Individual leaves were detached and scanned using an Epson GT-20000 flat bed scanner. From scanned leaves, we measured the area of the entire compound leaf using custom macros in ImageJ (Abràmoff et al. 2004). After scanning, we measured the fresh mass of only the terminal leaflet using a fine balance. Whole leaves were dried in a desiccating oven at 60° and weighed to determine dry mass. Thus, we measured LDMC on the terminal leaflet while LMA is for the entire leaf. Some leaves developed mold in the drying oven because of overcrowding. We statistically removed its effect before estimating the phenotype of a given genotype. LDMC and LMA were used to generate a proxy for leaf thickness (T), as defined by Vile et al. (2005):
(1)T was not measured directly because a thickness gauge is unreliable for wild tomatoes (C. D. Muir, personal observation) and sectioning samples to accurately measure T was not feasible for a project of this scale. However, in wild tomatoes, Equation 1 is highly correlated with T (r2 = 0.67) measured directly from sections (Muir et al. 2014a). For simplicity, we hereafter refer to leaf thickness (T) even though we are actually reporting a proxy of T.
Stomatal and trichome traits were measured on a subset of the same leaves as morphological traits, on all ILs. Leaf surface impressions were made from fresh tissue by applying a thin layer of nail polish to a medial portion of the terminal leaflet away from major veins. Nail polish layers were removed using adhesive tape and mounted on a microscope slide. Stomata were counted from three microscope fields each of the abaxial (“lower”) and adaxial (“upper”) leaf surfaces. The SR was calculated as ratio of the SDs on the adaxial and abaxial side:
(2)Trichome traits were measured by direct counts on the same three microscope fields and similarly used to calculate TD and TR. These traits were included to evaluate the strength of general associations between leaf epidermal traits and to assess an hypothesis that changes in SR could be a byproduct of selection on other developmentally integrated epidermal traits, such as trichomes. Tomato leaves are bifacial, with distinct internal (mesophyll) layers (Kebede et al. 1994; Galmés et al. 2013) and epidermal anatomies on each surface. However, if general developmental pleiotropy determines the anatomy of abaxial (lower) and adaxial (upper) epidermises, then species with high SR (i.e., more even distribution of stomata between upper and lower leaf surfaces) should show similar patterns in other epidermal traits such as trichomes.
QTL analysis
We determined QTL using (generalized) linear models. Where appropriate, we transformed variables or used models with non-Gaussian residual distributions. In particular, we log-transformed leaf morphological traits (LA, LMA, and T) to remove heteroscedasticity and normalize residuals. Therefore, we used linear models assuming normally distributed residuals on log-transformed traits. We modeled counts (SD and TD) and ratios (SR and TR) residuals as Poisson and binomial, respectively. In our data, SD, SR, and TD residuals were overdispersed, so we used quasipoisson and quasibinomial distributions that include a dispersion factor. Using a quasilikelihood method reduced the number of statistically significant QTL compared to using standard likelihood, but did not change estimated effect sizes. We fit the data on a log-link and logit-link scale for Poisson and binomial families, respectively. In addition to genotype, we estimated the effects of seedling treatment (control vs. cold), developmental stage (early vs. late), and leaf age (continuous). For traits based on dry mass (LMA and T), we also accounted for effects of mold on dry mass by including a factor with three levels (“no mold,” “slight mold,” and “mold”). We calculated IL effect size in terms of percentage of phenotypic variance explained using type II sum of squares, percentage of interspecific difference, and standard deviations of the Slyc parent. For non-Gaussian traits, we calculated percentage of deviance accounted for, the analog of percentage of variance explained. Analyses were conducted using the lm or glm function in R (R Development Core Team 2012).
Coefficients of the generalized linear mixed [(G)LM] model provide estimates of the effect of an introgression on the trait compared to Slyc. We determined statistical significance using the t-statistic. Because of multiple comparisons (70 IL v. Slyc comparisons per trait) involved in this analysis, we used strong control of false discovery rate (Benjamini and Hochberg 1995) as implemented in the R package multtest (Pollard et al. 2013) to calculate adjusted P-values. Only ILs that were statistically significant (adjusted P < 0.05) were used to identify QTL. QTL were identified as the minimum number of chromosomal regions needed to account for the empirical pattern of IL phenotypes, taking into account that some introgressed regions are shared between different ILs. When two ILs with overlapping introgressions both have a significant trait effect, then a parsimonious explanation is that a single QTL located in the overlapping region is responsible (also called “bin mapping”) (Pan et al. 2000). Following the criteria of Holtan and Hake (2003), who used this logic to define 107 nonoverlapping bins in the Slyc genome, we used an analogous procedure to define bins for our measured ILs, using IL boundaries obtained from the Sol Genomics Network (Bombarely et al. 2011). This procedure is conservative, since it identifies the minimum number of QTL needed to account for ILs with significant phenotypes.
We also calculated broad sense heritability (H2) as the proportion of phenotypic variance attributable among ILs after accounting for seedling treatment, developmental stage, leaf age, and mold. Variances were estimated using (G)LM models treating IL as a random effect. We used the same distributions to model residuals as in the QTL analysis (Gaussian for LA, LMA, and T; Poisson for SD and TD; and binomial for SR and TR). Models were fit using a Bayesian MCMC method implemented in the package MCMCglmm (Hadfield 2010). We used default priors for all traits except TR, where a more diffuse prior yielded reasonable estimates and better mixing. Posterior distributions were estimated from 5000 samples taken every 10 steps in a 50,000 step chain following a 5000-step burn in.
QTL sign test
The QTLST calculates the probability that the observed distribution of effect sizes and directions would have arisen in the absence of directional selection (Orr 1998b). We implemented QTLST rather than the more commonly used QTLEE (QTL equal effects) method (e.g., Rieseberg et al. 2002) because we had quantitative estimates of the effect size distribution (see below) and QTLST is less susceptible to ascertainment bias than QTLEE (Anderson and Slatkin 2003). For comparison, we also report probabilities from QTLEE, since most studies use this method. We first fit an exponential effect size distribution using maximum likelihood to all significant QTL for each trait separately. An exponential distribution is predicted by theory and often approximates real data well (Orr 1998a, but see Rockman 2012). For some traits, the exponential distribution did not closely fit the observed distribution of effect sizes, likely because we had insufficient power to detect small-effect QTL. However, QTLST assuming a log-normal rather than exponential distribution of effect sizes yielded very similar probabilities (data not shown), indicating that our results are not sensitive to a specific distribution. Calculating the critical probability for the QTLST with exponentially distributed effect sizes requires integrating the cumulative distribution of the gamma distribution (Orr 1998b), for which there is no closed form solution. We therefore numerically integrated probability densities using the distr package (Ruckdeschel et al. 2006) and integrate function in R (R Development Core Team 2012). This code is available upon request from the authors.
Genetic correlations
We estimated genetic covariation between traits by regressing IL means. Nonindependence between ILs due to overlapping introgression segments should increase spurious covariation between traits. However, this bias toward increased genetic covariation goes against our hypothesis that trait correlations are caused by natural selection rather than pleiotropy or linkage. Genetic correlations might also be evolutionarily relevant if one or a few loci affect multiple traits. Therefore, we also tested for significant colocalization of QTL for different traits using a one-tailed Fisher’s exact test (Fisher 1922). For both tests, we adjusted P-values to correct for multiple comparisons using the same Benjamini and Hochberg (1995) procedure described above. For carbon isotope discrimination (δ13C), we used previously published data from 49 of the same ILs discussed here (Xu et al. 2008). We assessed the correlation between SD, SR, and δ13C using linear regression of IL means, as with other traits. The correlation between these different stomatal traits and δ13C can be used to infer which portion of the CO2 diffusion pathway is more functionally significant in these species. Although it is generally assumed that carbon isotope discrimination varies with stomatal conductance, which is proportional to the density and aperture of stomata (Franks and Farquhar 2001; Sack et al. 2003), many other traits can also influence δ13C (Farquhar et al. 1982, 1989; Evans et al. 1986). Stomatal density (which we use as a proxy for stomatal conductance) affects the change in CO2 concentration from the atmosphere to the inside of the leaf. Stomatal ratio affects the CO2 path length from stomata to chloroplasts (greater SR reduces path length). Therefore, both SR and SD potentially affect water-use efficiency and carbon isotope discrimination, but lead to different resulting correlations between WUE and δ13C. Lower SD increases WUE and reduces discrimination, whereas greater SR increases WUE but also increases discrimination (see Discussion). Hence, the correlation between these different stomatal traits and δ13C can be used to infer whether diffusion from atmosphere → inside of leaf (SD) vs. from inside of leaf → chloroplasts (SR) is more functionally significant in these species.
Stomatal development candidate gene analysis: colocalization with QTL and patterns of molecular evolution
Focusing on traits whose QTL show the strongest evidence for the action of natural selection, we looked for colocalization between stomatal trait QTL and stomatal development candidate genes that were identified without prior knowledge of QTL location. Stomatal development candidate genes were identified from the well-characterized pathway in Arabidopsis (Pillitteri and Torii 2012) supplemented with additional select candidate genes from sterol and abscisic acid (ABA) biosynthesis pathways. In Arabidopsis, sterol biosynthesis affects stomatal development (Qian et al. 2013) through BIN2, which responds to brassinosteroid and phosphorylates multiple stomatal development genes (Gudesblat et al. 2012; Kim et al. 2012); therefore we included BIN2 orthologs in our analysis. ABA is a stress response hormone that regulates stomatal aperture, but might also regulate stomatal development (Tanaka et al. 2013); therefore, we also investigated ABA biosynthetic and catabolic enzymes (Nambara and Marion-Poll 2005). We considered 25 candidate genes in total (Table 4).
We determined orthologs in tomato for the Arabidopsis stomatal genes by a BLASTn search of the tomato reference transcripts (ITAG v2.4, www.solgenomics.net). The top match for each gene was selected as the ortholog with a minimum e-value of 1 × 10−20. Genes where several top matches were equally likely (often the case in large gene families) were removed from consideration. S. lycopersicum reference orthologs determined for each candidate stomatal gene were located on the tomato physical map using the SGN browser (http://solgenomics.net/gb2/gbrowse/ITAG2.4_genomic/). Based on known breakpoints for each of our introgression lines, we determined whether the physical position of each candidate fell anywhere within the introgressed Spen region of each IL that had significant effects on our stomatal traits, regardless of physical overlap with other significant ILs. We used this more liberal threshold for calling QTL boundaries than used for counting QTL (see above) to ensure that we identified all possible candidate loci and to account for the possibility that multiple separate loci fall within phenotypically significant ILs.
S. lycopersicum candidates for each stomatal gene were also used to search for their orthologs in other Solanum species using a BLASTn search of a combined Solanum library we constructed using previously published sequence data and/or transcripts from close relatives of S. lycopersicum. S. pennellii (LA0716) and S. habrochaites (LA1777) RNA-Seq reads were obtained from Koenig et al. (2013). These paired-end Illumina reads were trimmed using both Scythe (https://github.com/vsbuffalo/scythe) to remove adapters and a custom Python quality/length FASTQ trimmer. Trimmed reads were then assembled de novo into transcripts using Trinity (Grabherr et al. 2011) with standard parameters. The resulting set of assembled transcripts was combined with coding sequences downloaded from ftp.solgenomics.net for S. corneliomulleri (labeled S. peruvianum), S. lycopersicum (SL2.40 reference), S. pimpinellifolium (reference v.1.0), S. dulcamara, and S. tuberosum (v3.4) (Potato Genome Sequencing Consortium 2011) to form a working coding sequence library.
Our S. lycopersicum reference orthologs for each candidate stomatal gene were used as the query for a BLASTn search of this combined Solanum library. Only unambiguous orthologs with e-value < 1 × 10−50 were retained for further analysis. In two cases (NCED1 and AO3), S. tuberosum sequences were not annotated and were derived from a strong BLAST match to unannotated chromosomal sequences in the potato reference. Coding sequences for each candidate gene were aligned separately using PRANK (codon mode, default settings) with minor manual corrections (Löytynoja and Goldman 2008). The low divergence between these Solanum sequences ensured confident and high-quality coding sequence alignments. We then analyzed these alignments using the codeml package in PAML (v4.7; Yang 2007) with the independent branches model (all branches dN/dS independent), frame-specific codon frequencies, fixed tree topology (Supporting Information, Figure S1) with initial branch lengths equal to 1. Since the sequences have extremely low divergence, we also used codeml ancestral sequence reconstructions to count nonsynonymous and synonymous changes on each branch of the tree. Collation of sequences and PAML results were carried out with custom Python scripts. Alignments are available from Dryad Digital Repository (Muir et al. 2014b).
Results
Numerous QTL underlie species differences in heritable leaf traits
The focal species (Slyc and Spen) differed significantly (nonoverlapping 95% confidence intervals) in most traits (Figure 1A). Spen had smaller, thicker leaves, and greater leaf mass per area. Spen also had slightly higher stomatal density, but a much greater stomatal ratio, meaning that stomata were distributed evenly across both ab- and adaxial surfaces compared to Slyc, where most stomata are on the abaxial (lower) surface (Figure 1). Differences in trichome traits were less pronounced, although Spen had a slightly greater trichome ratio than Slyc. Observed differences between Spen and Slyc in these traits are consistent with other studies from our lab and others (Kebede et al. 1994). In contrast, some of our detected trait differences between Slyc and Spen might reflect trait sensitivity (plasticity) to the specific growth conditions in our experiment. In particular, we observed that Spen had greater LMA than Slyc, which is inconsistent with some other studies from our lab. The discrepancy is probably a result of different light intensities during experimental cultivation; LMA typically responds very strongly to irradiance (Poorter et al. 2010), and we consistently find that Slyc has lower LMA compared to Spen in our greenhouse (where irradiance is 200–500 photosynthetic photon flux density, PPFD), but not in experimental fields (where irradiance is 1500–2000 PPFD) (C. D. Muir, personal observation). Because more trait plasticity is observed in the Slyc genotype—LMA in Spen is less variable across cultivation conditions, but Slyc has low LMA under low light, and high LMA under high light—we infer that LMA is more plastic with respect to irradiance in Slyc than Spen.
Species’ differences (A) and the distribution of IL phenotypes (B–H) for seven leaf traits: leaf area (LA), leaf mass per area (LMA), leaf thickness (T), stomatal density (SD), stomatal ratio (SR), trichome density (TD), and trichome ratio (TR). (A) The 95% confidence intervals around the mean trait value for the recipient S. lycopersicum var. esculentum (Slyc; open bars) and donor S. pennellii (Spen; solid bars) in units of Slyc standard deviations (σSL) with the Slyc mean set to zero. (B–H) Histograms of 70 IL phenotypes for seven leaf traits. Shaded bars represent ILs that were not significantly different from Slyc after adjusting for multiple comparisons. Solid bars are lines that were significantly different. Slyc and Spen point to species’ means, but these are not included in the histogram bars.
For our introgression lines, H2, the number of QTL, the distribution of effect sizes, and the amount of transgressive segregation varied considerably across traits. Leaf traits had low to moderate broad sense heritability (LA, 0.08; LMA, 0.17; T, 0.12; SD, 0.24; SR, 0.27; TD, 0.39; and TR, 0.53) but all were significantly different from zero based on 95% highest posterior density intervals. Between 2 and 30 ILs per trait were significantly different from Slyc after correcting for multiple comparisons. A detailed list of statistically significant ILs, their genomic location, additive genetic effect, and percentage of interspecific divergence explained are given in Table S1. From IL phenotypes, we identified between 2 and 23 QTL per trait (Table 2). A complete description of QTL, names, and locations are provided in File S1 and shown graphically in Figure 2 (SR QTL only) and Figure S1, Figure S2, Figure S3, Figure S4, Figure S5, and Figure S6 (remaining traits).
Stomatal ratio QTL positions and ILs observed. To the left of each chromosome is the position of S. pennellii introgressions in this IL mapping population. Only ILs that were significantly different than S. lycopersicum var. esculentum after adjusting P-values for multiple comparisons are labeled (e.g., “IL1-1”). Narrow regions (dark blue boxes) were determined using bin mapping; broad regions (light blue boxes) show the extent of all contiguous ILs that significantly affected the trait (see Materials and Methods). QTL names are given on the right side of chromosomes. Similar maps for other traits can be found in Figure S1, Figure S2, Figure S3, Figure S4, Figure S5, and Figure S6.
Only two QTL were detected for leaf area (Figure S3), and both were transgressive: they conferred increased leaf size, which is opposite to the direction of the Spen donor phenotype (Figure 1B). In contrast, there were a substantial number of QTL for leaf mass per area (Figure S4) and leaf thickness (Figure S5), the majority of which did not transgress parental species’ means (Figure 1, C–D). We detected a large number of QTL for stomatal density (Figure S6) and stomatal ratio (Figure 2) compared to trichome density (Figure S7) and trichome ratio (Figure S8), respectively. With one exception (sd3.1), stomatal density QTL were of small to moderate effect, always in the direction of Spen, and not highly transgressive (Figure 1E). Stomatal ratio was characterized by two large-effect QTL (sr3.1 and sr9.1, explaining 15.9 and 5.9 percent of the phenotypic variance, respectively; Table S1), in addition to several small-effect QTL (Figure 1F).
Some observed differences in the genetic complexity of leaf traits depended on genetic and experimental factors. In particular, the number of ILs that differed from Slyc (and therefore the number of QTL per trait) was influenced by the heritability, the phenotypic difference between species, and the amount of variation within Slyc. For example, although stomatal density did not differ greatly between species, there was very little variation within Slyc for this trait, allowing us to detect many small effect QTL (Figure 1E and Figure S6). Compared to stomatal density, there was greater variation in stomatal ratio within Slyc; however, because species differed by >3 standard deviations, several moderate and large effect QTL were still detected for SR (Figure 1F and Figure 2). In contrast to stomatal traits, we detected few QTL for leaf area, despite the fact that LA differs significantly between Spen and Slyc. This was either because LA was determined by many small-effect QTL below our detection threshold and/or because of low heritability (high variance among replicate genotypes) possibly due to environmentally induced plasticity (as discussed above, for LMA). In addition, because we used introgression lines, we were unable to measure epistasis among different chromosomal regions from the same genetic background; therefore, trait differences that are determined by strong epistatic effects among unlinked loci will likely have undetected QTL under this experimental design.
QTL sign test indicates directional selection on stomatal traits
For five of seven traits, all or almost all QTL effects were in the direction of the donor parent Spen (Table 2). A QTL sign test showed that the distribution of QTL effect sizes for stomatal traits SD and SR were unlikely to have arisen in the absence of directional selection (Table 2). We could not reject a neutral model for LA, LMA, T, TD, and TR, although using less stringent criteria for counting QTL would have made LMA and T significant (data not shown); thus, a more powerful study might have rejected the neutral hypothesis for these traits. The more commonly used QTLEE gave qualitatively similar results to QTLST, but P-values were approximately two times lower (Table 2), indicating that QTLST is more conservative than QTLEE (unlike stated in Orr 1998b).
Most leaf traits are not genetically correlated
There was little evidence for genetic correlations between traits (Table 3). Among our focal traits (LA, LMA, T, SD, and SR), only two correlations were significant. Leaf mass per area and leaf thickness were strongly positively correlated, which is expected because LMA is by definition the product of leaf thickness and bulk tissue density (Witkowski and Lamont 1991). Stomatal density and stomatal ratio were also positively correlated, but this finding was strongly influenced by a single pleiotropic QTL, sr3.1, that had a normal abaxial (lower) stomatal density, but high adaxial (upper) stomatal density, and thus high SR. The correlation between SD and SR is nonsignificant if this QTL is removed from the analysis (r = 0.04, P = 0.38). Interestingly, with the exception of sr3.1, abaxial and adaxial stomatal density appear to be genetically decoupled from one another. Across ILs, there was no correlation between ab- and adaxial SD (r = 0.18, P = 0.13). Because most stomata were on the abaxial surface (similar to Slyc) genetic variation in SD was determined primarily by variation in stomatal density on this lower leaf surface (Figure 3A). In contrast, variation in SR was determined almost entirely by variation in stomatal density on the upper (adaxial) leaf surface (Figure 3B).
(A) Abaxial (lower) stomatal density is the most important determinant of stomatal density (SD) of the entire leaf, but (B) adaxial (upper) stomatal density determines stomatal ratio (SR). Each point is the phenotype of one IL. In both A and B solid points are ILs that had significantly greater SR than Slyc after adjusting for multiple comparisons. A major QTL for both SR (sr3.1) and SD (sd3.1) found on IL3-2 is labeled to show that the same locus has a large effect on both traits. In A, the line was fit and r2 calculated using linear regression on untransformed stomatal counts. In B, the line was fit and r2 calculated using linear regression between untransformed stomatal counts and SR on a logit scale, but data are plotted on an untransformed scale.
Several observations indicate that stomatal patterning is genetically independent of trichome patterning. First, there was no correlation between stomatal and trichome traits (Table 3). Second, the relationship between trichome density and trichome ratio was very different from that between stomatal density and stomatal ratio: The strong correlation between TD and TR was not driven by a few influential QTL (data not shown) and ab- and adaxial TD were highly positively correlated (r = 0.83, P < 0.001). Although we had no a priori prediction, trichome ratio was negatively correlated with both leaf mass per area and thickness (Table 3), which could reflect an unknown connection between trichome patterning and other aspects of leaf development or could be a false positive. The results from Fisher’s exact tests to examine evidence of QTL colocalization closely mirrored genetic correlations, except that there were no cases of significant colocalization after adjustment for multiple comparisons (Table S2). We also did not find evidence that any particular IL harbored a significant excess or deficit of QTL; given the number of QTL and ILs, the distribution of QTL across ILs did not depart from chance 
.
Stomatal ratio, not stomatal density is correlated with carbon isotope discrimination
We found that stomatal ratio, but not stomatal density, was significantly correlated with carbon isotope discrimination (δ13C) in these ILs (Figure 4). Our data indicate that for these species, when grown under well-watered greenhouse conditions, the CO2 diffusion path length (mediated by SR) has a significant effect on δ13C, but variation in stomatal conductance (using SD as a proxy) does not. In particular, leaves with greater SR exhibited greater discrimination (i.e., more negative δ13C; Figure 4B), indicating that having stomata evenly distributed on both leaf surfaces significantly affects δ13C by shortening the average distance between stomata and chloroplasts. As the diffusion pathway shortens, more of the fractionation occurs during carboxylation by Rubisco rather than during gas diffusion. Since Rubisco discriminates against 13C much more than does gaseous diffusion (27–30‰ vs. 4‰) (Farquhar et al. 1989), leaves with shorter path lengths should have a more negative δ13C, as we observe. One exception to the overall pattern between stomatal traits and δ13C is IL5-4, which had much lower discrimination than expected based on the stomatal ratio observed in this genotype (Figure 2 and Figure 4B); the δ13C estimates in this particular genotype might be due instead to lower stomatal conductance, as concluded by Xu et al. (2008).
Carbon isotope data suggest that amphistomy (greater stomatal ratio) decreases the distance between stomata and chloroplast. (A) Carbon isotope discrimination is uncorrelated with stomatal density (SD), but (B) positively correlated with stomatal ratio (SR) across ILs (more discrimination equals more negative δ13C). Each point is an IL. Slyc and Spen means are indicated for comparison, but these were not used in the regression. Carbon isotope data were obtained from Xu et al. (2008); stomatal data are from this study.
We are uncertain why stomatal density did not generally affect δ13C, but two possible explanations are that SD was not a good proxy for stomatal conductance, which is also affected by stomatal aperture, or that the range of variation in stomatal conductance among lines was too small to have an appreciable affect on δ13C. When there are large differences in average stomatal conductance between plants (e.g., when drought induces stomatal closure), variation in δ13C correlates strongly with WUE. However, when plants are grown under nonstressed conditions, stomatal conductance is uniformly high; therefore, differences in δ13C of a few permil between genotypes can be affected by several traits other than stomatal conductance. As with our analysis here, in their study Xu et al. (2008) measured δ13C under nonstressed conditions.
Multiple candidate genes colocalize with stomatal QTL
Of the candidate genes drawn from the stomatal developmental biology literature (Pillitteri and Torii 2012), several colocalized with stomatal QTL (Table 4), including QTL with the largest effect sizes (sr3.1, sr9.2, and sd3.1). In addition, we observed elevated rates of nonsynonymous compared to synonymous substitution (dN/dS > 1) in two of our candidate genes (ER and STOMAGEN) specifically on the branches leading to our focal species, Slyc and Spen, perhaps indicating a history of positive selection. Interestingly, neither of these genes colocalized with the largest effect QTL. We found no colocalization between QTL and ABA biosynthesis protein genes, suggesting that constitutively altered ABA levels are not responsible for stomatal trait evolution.
Discussion
It is tempting to ascribe the phenotypic differences between species to the action of natural selection, but mutation and drift can also cause phenotypic divergence without adaptation. In this study, we looked at several leaf traits that are believed to be involved in climatic adaption and ecological divergence between species. In particular, we considered two adaptive hypotheses: first, that phenotypes diverged because of directional natural selection, and second, that traits are correlated across species because of selection rather than genetic correlations caused by pleiotropy or linkage. We used QTL mapping to identify the genetic architecture of phenotypic divergence and found evidence supporting our adaptive hypotheses, especially for stomatal traits. In addition, by combining our trait variation data with information on physiological variation and genomic position, we were able to directly assess the physiological significance of stomatal trait variation and uncover candidate genes that point toward particular developmental mechanisms underlying adaptation.
Genetic architecture of species differences suggests different evolutionary dynamics and trait lability
We detected between 2 and 23 QTL for each of our putatively adaptive leaf traits, each of which altered the phenotype between 0.12 and 2.10 Slyc trait standard deviations (Table S1). This large range of genetic complexity could result from several genetic factors, including the magnitude of phenotypic difference between species, the amount of variation within Slyc, and trait differences in heritability (Results). In some cases, low heritability appeared to be affected by environmentally induced plasticity, especially for traits related to leaf area (LA and LMA; Results), which are known to vary plastically with growth environment. As such, QTL identified for these traits in our study might be associated with genetic differences in plasticity rather than constitutive differences. A better understanding of G × E interactions for these traits will clarify whether these QTL also underlie species differences in nature.
The different genetic architectures we observe suggest that certain traits might evolve more readily in response to new environmental conditions. Many traits in tomato and other species, including some traits in this study, are highly polygenic (Lippman et al. 2007), which could be consistent with either weak phenotypic selection, significant pleiotropy, and/or lack of large-effect mutations. In contrast, the two large-effect loci for one of our traits (stomatal ratio) more closely resemble the genetic architecture of traits that are known to have evolved under strong selection with limited pleiotropy such as flower color in Mimulus (YUP and MYB transcription factors) (Schemske and Bradshaw 1999; Streisfeld et al. 2013; Yuan et al. 2013), armor plating in stickleback (EDA) (Colosimo et al. 2005), and coat/skin color in mammals (MC1R) (Nachman et al. 2003; Hoekstra et al. 2006). Like these traits, few genetic changes appear necessary to achieve large phenotypic transitions in SR; together, our two large effect SR loci account for 98% of the phenotypic difference between our two parental species (assuming they act additively). The similarity in genetic architecture suggests that SR is likewise an ecologically important trait under divergent selection in tomatoes, although we have not yet identified the specific selective agent.
Another inference from the striking genetic architecture of SR is that the availability of large-effect loci might make this trait more labile in response to changes in selection. There is some qualitative support for this hypothesis. Variation in SR among wild tomato species spans the entire range of variation observed across the Solanaceae, which contains numerous hypo- and amphistomatous species (C. D. Muir, unpublished data). In comparison, variation in LA among wild tomatoes, for which we found no large-effect loci, comprises a small fraction of the total variation in Solanum and has high phylogenetic signal (C. D. Muir, unpublished data). In the future, this kind of hypothesis can be tested quantitatively by examining macroevolutionary patterns of trait lability across broader taxonomic groups, such as all Solanum (∼1500 species), in response to environmental niche transitions. In wild tomatoes, the evolution of amphistomy is closely correlated with transitions to drier environments, an important component of niche evolution between sister species (Nakazato et al. 2010). Although the ancestral SR of all wild tomatoes is ambiguous (outgroups are both hypo- and amphistomatous; C. D. Muir, unpublished data), amphistomy likely evolved in Spen as it colonized desert habitats; Spen has the highest SR of any tomato species and its sister species (S. habrochaites) is hypostomatous (Muir et al. 2014a and unpublished data), indicating that a trait shift to more equal distribution of stomata on both leaf surfaces occurred specifically within the Spen lineage.
Genetic architecture reveals evidence for natural selection on individual traits and trait combinations
A clear preponderance of QTL had allelic effects in the direction of the Spen donor phenotype (Table 2), and a model assuming random, nondirectional fixation was strongly rejected for stomatal traits (SD and SR), consistent with trait differentiation due to divergent natural selection. Because the number and size of detected QTL limits the power of the QTLST test, failure to reject a model of random fixation for several of our other traits does not rule out the action of natural selection. For example, smaller leaf area correlates with precipitation specifically in wild tomatoes, across the entire genus Solanum (C. D. Muir, unpublished data) and in many other clades (McDonald et al. 2003). However, because we did not detect QTL that reduced LA in the direction of Spen, QTLST could not be performed for this trait. In such cases, especially where there is substantial independent evidence that a trait is adaptive, alternative methods such as selection gradient analysis will be needed (Dudley 1996). Regardless, our analysis clearly demonstrated selection on stomatal traits, either because selection on stomatal traits has been stronger than that on leaf morphology or because the larger effect sizes of stomatal trait QTL made them easier to detect and hence increased our power to reject the random fixation model.
Comparison of genetic architecture across traits also provided evidence that trait correlations across species are due to selection rather than pleitropy or linkage. In wild tomatoes, leaf area, leaf mass per area, and stomatal ratio are correlated with each other and with climate, especially precipitation (C. D. Muir, unpublished data). While one explanation for these correlations is shared underlying mechanisms, we found no evidence that leaf traits are correlated across species because of genetic correlations. The only evidence for pleiotropy/linkage came from a single major QTL on chromosome 3 that affected both stomatal density and stomatal ratio (Figure 3). In the absence of this QTL, however, the remaining SD and SR QTL are uncorrelated in our analysis. Similarly, across wild tomatoes in general, SD and SR are uncorrelated (C. D. Muir, unpublished data). In fact, the wild tomato species with the highest stomatal density (S. habrochaites) has very low stomatal ratio (most of the stomata are on the abaxial surface). Overall, from these data we infer that leaf traits phenotypically covary across species because they function well together in certain ecological contexts; that is, these suites of covarying traits are the product of selection. In particular, small leaves, low LMA, and high SR are likely more common in wild tomato species from drier habitats because they confer greater photosynthetic efficiency and reduce water lost relative to carbon gain.
Finally, we also examined whether species differences in stomatal traits could be a pleiotropic effect of selection on other features of the plant epidermis—specifically trichomes—rather than the direct target of selection. We found that the density and distribution of stomata is independent of the density and distribution of trichomes. Although we did not measure additional epidermal traits, such as cell size or shape, these data suggest that stomatal traits did not evolve as a pleiotropic byproduct of changes in other traits, providing a second piece of evidence that stomatal traits per se are the target of divergent selection.
Structure–function relationship reveals the mechanistic significance of adaptive stomatal trait variation
By examining the association between variation in our stomatal traits and a previously measured metric of whole plant water use efficiency (Xu et al. 2008), we were able to address the precise physiological significance of leaf trait variation under our experimental conditions. The positive correlation we observed between stomatal ratio and carbon isotope discrimination suggests both physiological and ecological inferences: first, that having stomata on both sides of the leaf shortens the CO2 diffusion path length inside the leaf, thereby potentially influencing water-use efficiency, and second, that interspecific variation in SR is ecologically significant in wild tomatoes.
The physiological inference, that SR decreases the diffusion path length, differs from that typically made using δ13C in quantitative genetic studies of water-use efficiency. Generally, lower stomatal conductance increases water-use efficiency and reduces discrimination against the heavier isotope 13C, which is why carbon isotope discrimination is often used as a proxy for WUE. In contrast to this general expectation, here we found that variation in stomatal ratio, but not stomatal density (a proxy for conductance), predicted carbon isotope discrimination among ILs. This is surprising because it suggests that the CO2 diffusion path length from stomata to the chloroplast, rather than stomatal conductance, determines δ13C in these ILs.
The ecological inference is that selection has acted specifically on photosynthetic efficiency (rather than solely water conservation) during adaptation to low precipitation across wild tomato species. In wild tomatoes, stomatal ratio is typically greater in species native to drier climates, such as Spen (C. D. Muir, unpublished data). By reducing the distance between stomata and chloroplasts, greater SR can increase WUE by increasing the concentration of CO2 around the chloroplasts (and therefore the rate of photosynthesis) without incurring greater evaporative water loss. However, under this mechanism, greater discrimination would be indicative of greater WUE, which is exactly the opposite of how δ13C is interpreted when stomatal conductance is assumed to be the dominant mechanism affecting both δ13C and WUE. Our finding stresses the importance of identifying the mechanism underlying variation in δ13C in quantitative genetic studies (Juenger et al. 2005), when interpreting the ecological significance of this variation.
Candidate genes colocalized with stomatal trait QTL suggest specific molecular mechanisms of trait variation
Overall, our stomatal traits were most clearly and consistently associated with both ecological and physiological roles in climatic adaptation of wild tomatoes. To begin assessing the molecular mechanisms responsible for stomatal trait evolution, we identified candidate genes that colocalize with quantitative trait loci underlying these traits, especially for QTL of large effect for stomatal ratio (sr3.1 and sr9.1) and stomatal density (sd3.1) (Table 4). Our phenotypic data indicate that sr3.1 and sr9.1 both changed SR specifically by greatly increasing adaxial stomatal density (2 vs. 53 and 22 per field, respectively) without a concomitant reduction in abaxial density. These phenotypes imply a particular developmental mechanism underlies trait variation at both these loci, wherein Spen alleles confer adaxial expression of the entire stomatal developmental pathway but this expression is nearly absent in Slyc. This adaxial-specific change could arise from structural changes in proteins regulating stomatal development, cis-regulatory changes altering binding by ab-adaxial polarity factors, and/or altered sensitivity to hormonal regulators of stomatal development in the adaxial epidermis. Based on the identity and molecular evolution of our colocalized candidates, we find evidence in support of the latter two mechanisms only. In the first case, none of the candidate genes that colocalized with our large-effect QTL had a significant signature of adaptive evolution in either of the branches leading to Slyc or Spen, arguing against structural protein changes, although the analysis of selection used here (dN/dS) will not detect adaptive differences due to single or few nonsynonymous substitutions. In contrast, the specific molecular identity of some of our candidates suggest the second alternative—cis-regulatory changes—could underlie our observed phenotypic differences. In particular, we find that the tomato homolog of SPCH resides within sr3.1, our largest effect SR locus (sr3.1). In Arabidopsis, SPCH activity initiates stomatal development; therefore loss (gain) of adaxial stomata could be as simple as SPCH losing (gaining) regulation by a factor responsible for ab-adaxial polarity. Moreover, several described mutations in the cis-regulatory region of SPCH in Arabidopsis can completely eliminate stomatal development specifically on one side of the leaf, but not the other, indicating that SPCH interacts with ab-adaxial polarity factors (Dow et al. 2014). We also found that SDD1 and FACKEL colocalized with sr9.1. The products of these genes modulate stomatal density, and hence altered expression in adaxial tissue could alter stomatal ratio.
Therefore, based on current functional information on colocalized stomatal development genes, we find that mutation(s) in cis-regulatory elements are the most likely molecular mechanism underlying the trait change associated with our largest effect SR QTL. These hypotheses will require future functional analyses, which are well developed in tomato (Gupta et al. 2009). Regardless, despite intense interest in stomata because of their role in adaptation to climate change, ours is among the first studies on the genetics of natural variation in stomata density and distribution (Hall et al. 2005; Rae et al. 2006; Gailing et al. 2008), and the first study of stomatal ratio genetics.
Conclusions
Based on fundamental biophysical principles and empirical trait–environment correlations, we hypothesized that interspecific variation in several ecologically important leaf traits reflects adaptation to climate in wild tomatoes. In all leaf traits, there was significant heritability and we identified QTL underlying phenotypic differences between two tomato species. Most QTL affected traits in the direction of the donor parent, the wild species S. pennellii, suggesting that loci were fixed by natural selection rather than genetic drift, and this pattern was strongly supported for stomatal traits. We found no evidence that focal leaf traits that correlate across species are genetically correlated, or that stomatal traits were shaped by selection on other epidermal features. We infer from this that association between traits in nature likely reflects natural selection. We also found evidence that variation in the distribution of stomata is functionally significant and could affect ecologically important traits like water-use efficiency. Collectively, these data along with ecophysiological and comparative data on variation in leaf traits, make a compelling argument for adaptation between recently diverged species in this ecologically and physiological diverse clade. For key traits such as SD and SR, our analysis also suggests that these ecologically important transitions can involve relatively few changes of large phenotypic effect, that allow rapid responses to new environmental conditions. Because this leaf trait variation also captures important axes of phenotypic variation across all land plants, these results also suggest mechanisms that might underlie the evolution of these traits across much broader macroevolutionary scales.
Acknowledgments
M. Yakub and S. Vosters assisted with germinating, transplanting, and maintaining plants in the greenhouse. E. Josephs, S. Josway, and E. Levin assisted with phenotyping. This research was supported by a National Science Foundation (NSF) graduate research fellowship to C.D.M. and NSF grants DEB-0841957 and DEB-1135707 to L.C.M.
Footnotes
Supporting information is available online at http://www.genetics.org/lookup/suppl/doi:10.1534/genetics.114.169276/-/DC1.
Communicating editor: A. H. Paterson
- Received August 1, 2014.
- Accepted October 4, 2014.
- Copyright © 2014 by the Genetics Society of America
