The genetic basis of traits associated with the evolution of serpentine endemism in monkeyflowers

Abstract

The floras on chemically and physically challenging soils, such as gypsum, shale, and serpentine, are characterized by narrowly endemic species. The evolution of edaphic endemics may be facilitated or constrained by genetic correlations among traits contributing to adaptation and reproductive isolation across soil boundaries. The yellow monkeyflowers in the Mimulus guttatus species complex are an ideal system in which to examine these evolutionary patterns. To determine the genetic basis of adaptive and prezygotic isolating traits, we performed genetic mapping experiments with F2 hybrids derived from a cross between a serpentine endemic, M. nudatus, and its close relative M. guttatus. Few large effect and many small effect QTL contribute to interspecific divergence in life history, floral, and leaf traits, and a history of directional selection contributed to trait divergence. Loci contributing to adaptive traits and prezygotic reproductive isolation overlap, and their allelic effects are largely in the direction of species divergence. These loci contain promising candidate genes regulating flowering time and plant organ size. Together, our results suggest that genetic correlations among traits can facilitate the evolution of adaptation and speciation and may be a common feature of the genetic architecture of divergence between edaphic endemics and their widespread relatives.

Introduction

Geology influences the composition, abundance, and traits of species, exemplified by dramatic floristic transitions across soils. Soils that are chemically and physically challenging (e.g., shale barrens, granite outcrops, serpentine) typically harbor unique floras characterized by narrowly endemic species (McVaugh, 1943; Wherry, 1930; Whittaker, 1954). These harsh soils can facilitate the persistence of older species as refugia and the evolution of new species (Stebbins, 1942; Stebbins & Major, 1965). Many edaphic (soil-related) endemics have close relatives that occur nearby and are not restricted to specific soil types, suggesting that new species often evolve from widespread progenitors across soil boundaries (Anacker & Strauss, 2014; Grossenbacher et al., 2014). Yet edaphic endemics are typically outnumbered by widespread species on their specialized soils, which may reflect decreased diversification rates in endemics relative to widespread species (Anacker et al., 2011) or the difficulty of evolving reproductive isolation with gene flow (Caisse & Antonovics, 1978; Felsenstein, 1981). The genetic architecture of divergence may influence the probability of speciation across soil boundaries and coexistence upon secondary contact. If adaptive and reproductive isolating traits have a shared genetic basis, genetic correlations, either through pleiotropy or linkage, can facilitate the formation and maintenance of new species across soil boundaries (Hawthorne & Via, 2001). However, few studies have investigated the number, location, and effect size of loci contributing to reproductive isolation and adaptive divergence between edaphic endemics and their close relatives (Ferris et al., 2017; Macnair & Cumbes, 1989).

Whether endemics arise through the divergence of individuals that establish in a new edaphic environment from a widespread ancestral species or depletion of a previously widespread species composed of edaphic ecotypes, both routes to endemism begin with the evolution of locally adapted populations. Local adaptation is common in plant populations inhabiting different edaphic environments (Hereford, 2009; Sambatti & Rice, 2006; Wright et al., 2006) and can evolve rapidly (Davies & Snaydon, 1976; Wu et al., 1975). Classic evolutionary theory predicts that most adaptive substitutions have small or intermediate effects on fitness due to pleiotropy (Fisher, 1930) and drift (Kimura, 1983). More recent theory predicts an exponentially decreasing distribution of effect sizes when populations are far from their phenotypic optimum (Orr, 1998a), and that large effect loci play a major role in adaptation in the presence of gene flow (Griswold, 2006; Yeaman & Whitlock, 2011; reviewed in Dittmar et al., 2016), conditions that are likely in the evolution of local adaptation across soil boundaries. In addition, soil boundaries are often discrete, transitioning abruptly between soil types with no intermediates, a scenario where loci with small effects may be selected against because intermediate phenotypes have low fitness (Macnair et al., 1989; Selby & Willis, 2018). While edaphic adaptation may initially favor the fixation of alleles with large effects (e.g., loci contributing to soil tolerance), adaptive alleles subsequently fixed within edaphic endemics may have smaller effects. As locally adapted populations evolve reproductive isolation, gene flow becomes less likely to swamp divergence at smaller effect alleles. Furthermore, environmental differences within soil types may no longer be as extreme and abrupt as between soil types, and thus smaller effect alleles may be favored.

Edaphic ecotypes often differ in traits contributing to prezygotic reproductive isolation, like flowering time and mating system (Dittmar & Schemske, 2018; McNeilly & Antonovics, 1968; Sianta & Kay, 2021). Genetic correlations among adaptive and reproductive isolating traits may facilitate divergence between soil ecotypes and may be common between taxa inhabiting chemically and physically different soils. Loci contributing to serpentine soil adaptation (nickel tolerance, leaf area, and succulence) and premating isolation (flowering time) overlap between Silene vulgaris ecotypes (Bratteler et al., 2006). Similarly, loci contributing to divergence in adaptive (leaf shape) and prezygotic isolating traits (flowering time and flower size) overlap between a granite endemic and its widespread close relative (Ferris et al., 2017).

The Mimulus guttatus (syn. Erthranthe guttata (Fisch. DC.) G.L. Nesom) species complex, a group of closely related yet morphologically and ecologically diverse wildflowers, is an excellent system for investigating the genetic basis of adaptation and reproductive isolation across soil boundaries. Mimulus guttatus is widespread throughout western North America and is closely related to multiple recently derived geographically restricted species nested within its range. These close relatives are often edaphic endemics, including on granite, serpentine, and copper mines, suggesting that speciation occurs repeatedly across soil boundaries in this group (Ferris et al., 2014; MacNair et al., 1989; Selby et al., 2014). One example is M. nudatus, a diploid annual mixed-mating serpentine soil endemic restricted to the northern coast range of California. Serpentine soils, derived from metamorphic ultramafic rock, have low levels of Ca relative to Mg, are deficient in essential plant nutrients, often have an excess of toxic heavy metals, are rocky and shallow, and have a low water holding capacity (Brady et al., 2005). Populations of M. guttatus have repeatedly evolved serpentine tolerance (Selby & Willis, 2018), and these diploid annual mixed-mating serpentine-tolerant populations co-occur within meters of populations of M. nudatus (Gardner & MacNair, 2000; Oneal et al., 2016; Toll & Willis, 2018). Throughout its range, M. guttatus lives in streams and meadows, whereas M. nudatus grows in gravel washes and rocky outcrops. Consistent with their habitat affinities, in lab experiments, M. nudatus is more drought tolerant than M. guttatus (Hughes et al., 2001; Wu et al., 2010). The drier sites inhabited by M. nudatus tend to have fewer competitors and more bare ground than the more densely vegetated and wetter sites inhabited by M. guttatus, a common pattern observed in serpentine soil endemics when compared to species that live on and off serpentine soils (Sianta & Kay, 2019).

Plants adapted to serpentine soils usually have smaller leaves and shorter statures compared to non-serpentine populations or closely related species (Pichi-Sermolli, 1948; Rune, 1953; reviewed in Krukeberg, 1954). While serpentine adapted populations of M. guttatus are slightly shorter and have smaller rosettes than non-serpentine populations off serpentine soils, they are phenotypically very similar to off-serpentine populations (Selby & Willis, 2018). Consistent with broader floristic patterns, the serpentine-endemic M. nudatus is smaller than both non-serpentine and serpentine-adapted populations of M. guttatus and differs in several ecologically important traits: M. nudatus flowers under shorter critical daylengths and often flowers earlier in the greenhouse, has smaller flowers, and smaller and narrower leaves than M. guttatus (Figure 1, Ferris et al., 2015; Friedman & Willis, 2013; Selby et al., 2014; Sianta & Kay, 2021). Mimulus nudatus transitions to reproductive growth at an earlier developmental stage, reflected by earlier bolting (elongation between earlier forming leaf pairs) and flowering at earlier nodes (leaf pairs), resulting in reduced allocation to vegetative growth compared to M. guttatus. Small leaves and small flowers transpire at a slower rate than large leaves and consequently are associated with arid environments (Galen et al., 1999; Nicotra et al., 2011), early flowering is a common drought escape strategy in annual plants, and plants that avoid or escape drought often invest less in vegetative growth (Kooyers, 2015). Finally, while hybrid seed inviability is an extremely important and strong postzygotic barrier between these species, other reproductive barriers maintain species boundaries and facilitate coexistence in this species pair (Gardner & McNair, 2000; Oneal et al., 2016; Toll & Willis, 2018). Specifically, flowering time and flower size differences contribute to temporal and mechanical prezygotic isolation between them in sympatry (Gardner & MacNair, 2000; Sianta & Kay, 2021).

To determine the genetic architecture of traits contributing to adaptive divergence and prezygotic reproductive isolation between M. nudatus and M. guttatus, we performed a genetic mapping experiment with F2 hybrids derived from a cross between them. We crossed lines from a sympatric serpentine-tolerant population of M. guttatus with M. nudatus. Since serpentine and non-serpentine populations of M. guttatus are phenotypically very similar, this allowed us to focus on traits associated with occupation of drier and sparser microhabitats typical of serpentine endemics. Furthermore, prezygotic isolating traits between these species are more important in sympatry, where soil boundaries do not restrict hybridization. We considered traits as putatively adaptive if they have repeatedly evolved in serpentine endemics (xeromorphic leaves and short stature) and/or could potentially contribute to drought avoidance in the sparse dry outcrops inhabited by edaphic endemics (reduced allocation to vegetative growth via differences in developmental timing), but do not directly contribute to reproductive isolation in sympatry. Although flowering time and flower size could contribute to drought escape and avoidance, we considered these as traits contributing to reproductive isolation based on prior work demonstrating their role in temporal and mechanical prezygotic isolation (Gardner & MacNair, 2000; Sianta & Kay, 2021). With this experiment, we asked the following questions: What is the genetic architecture of traits associated with prezygotic reproductive isolation and putatively adaptive divergence associated with serpentine endemism? Do the loci contributing to traits associated with reproductive isolation and putatively adaptive divergence overlap? Are allelic effects at overlapping loci all in the direction of species divergence (facilitating divergence) or do they oppose each other (constraining divergence)? Did a history of consistent directional selection contribute to species divergence?

Materials and Methods

Plant material and experimental design

We created inbred lines from maternal families collected at a sympatric site at the UC Donald and Sylvia McLaughlin Natural Reserve in Lower Lake, California (W 122° 24.614, N 38° 51.528). Although we only used a single parental line per species in this experiment, these parental lines closely resembled other inbred lines of these species (Ferris et al., 2015; Friedman & Willis, 2013), and were very different from each other (resembling Figure 1). We created F1-hybrids by reciprocally crossing parental lines that were self-fertilized for two generations in the greenhouse. However, only F1-hybrid seeds derived from crosses where M. guttatus was the maternal parent successfully germinated. We created F2-hybrids for genetic mapping by self-fertilizing a single M. guttatus ♀x M. nudatus ♂ F1-hybrid. We planted seeds of the M. guttatus line REM122 (n = 21), the M. nudatus parental line REMn129 (n = 24), F1-hybrids (n = 10), and F2-hybrids (n = 576) in 2.5” Kord pots with Fafard 4P soil, cold stratified them for 5 days, and placed in an 18 hr light/6 hr dark greenhouse at Duke University. We fully randomized pots and flats were rotated weekly. Eighteen M. guttatus, 16 M. nudatus, 3 F1-hybrids, and 398 F2-hybrids germinated and were included in our phenotype analysis. Of these F2s, we extracted sufficient high-quality DNA from 379 F2-hybrids to generate reduced representation libraries for QTL mapping.

Phenotypic measurements and analyses

We recorded germination and flowering dates and calculated the days to flowering from germination. On the date of the first flower, we recorded the node of the first flower and measured the length of the first internode, plant height, the first open corolla width and length, and the first leaf length, width, and area. The node of first flower is the leaf pair at which a given plant produced its first flower, counting from the first true leaf pair after the cotyledons (0 = cotyledon node, 1 = first true leaf node, and so on). The length of the first internode is the length between the cotyledons and the first leaf pair and describes whether a plant bolted (transitioned to reproductive growth) early in development or not. We measured leaf traits (length, width, and area) with images produced using a flatbed scanner in ImageJ (Schneider et al., 2012). We calculated leaf shape by dividing leaf width from leaf length.

We tested whether the parental lines differed in the traits using Welch’s two-sample t-tests and assessed significance using Bonferroni-adjusted p-values. To estimate broad-sense heritabilities, we used parental and F2s variances to calculate environmental, phenotypic, and genetic variances. First, we calculated the environmental variance as the weighted average of the parental line variances, excluding F1s due to low sample size (F1 n = 3) (V_E = (var(REM122) + var(REMn129))/2). Then, we calculated genetic variance by subtracting V_E from V_P, the F2 phenotypic variance (V_G = V_P−V_E; Falconer & MacKay, 1996). Finally, we estimated broad-sense heritability by dividing the genetic variance by the F2 variance (H² = V_G/V_P). We estimated trait correlations in the F2 hybrids by calculating Pearson correlation coefficients.

DNA extractions and library generation

We extracted DNA from F2-hybrids and parents using a modified CTAB protocol (Lin & Ritland, 1995). We prepared reduced representation libraries for sequencing using a modified multiplexed shotgun genotyping protocol (Andolfatto et al., 2011; Supplementary Methods S1). All F2-hybrid samples were pooled in equal amounts (50 ng) and sequenced on a single lane of Illumina HiSeq 2000/2500 with 50-bp single-end reads. Parental samples were pooled in equal amounts (200 ng) and sequenced with four other samples on a single lane of Illumina Hi-Seq 4000 with 50-bp single-end reads.

Marker generation and linkage map

We created markers for genetic mapping from our sequencing data using the TASSEL 5.0 GBS v2 pipeline (Glaubitz et al., 2014; Swarts et al., 2014; Supplementary Methods S2). We constructed a linkage map using onemap in R (Margarido et al., 2007), ordering 1,011 markers along 14 linkage groups using a two-point-based algorithm, Rapid Chain Delineation (Doerge, 1996; Supplementary Methods S2). We tested for segregation distortion using a chi-square test and assessed significance using a Bonferroni correction.

QTL mapping

We initially performed standard interval mapping using the scanone function in R/qtl, then used the scantwo function to perform a two-dimensional two-QTL scan using Haley–Knott regression (Broman et al., 2003). We calculated the main effect and interaction penalties with 1,000 permutations of the scantwo function to conduct stepwise fitting of multiple QTL models. We used the stepwiseqtl function with a maximum of 10 QTL using Haley–Knott regression to conduct a forward/backward search of models that optimized the penalized LOD score criterion. We determined the 1.5 LOD interval for each QTL using the lodint function in R/qtl.

We used the fitqtl function in R/qtl to determine the percent variance explained by each QTL, peak LOD scores, and additive (a) and dominance (d) effects for each significant QTL from final stepwise models for each trait. For traits where the parental lines significantly differed, we also estimated QTL effect sizes by calculating the proportion of the mean parental divergence explained by each locus. We divided the difference between the mean phenotypes of the alternate homozygotes by the difference between the parental phenotype means and multiplied by 100.

We calculated the degree of dominance for each QTL by dividing the absolute value of the dominance deviation (d) by the additive effect (a) for each QTL. We binned QTL in categories based on the degree of dominance estimate: additive (0–0.20), partially dominant (0.21–0.80), dominant (0.81–1.20), and under- or overdominant (>1.20) (Stuber et al., 1987). For partially dominant and dominant loci, we examined which parental species was dominant.

Tests of QTL overlap

We tested whether loci underlying traits associated with reproductive isolation and putatively adaptive divergence overlapped using a permutation test (n = 1,000 permutations) in the R package regioneR (Gel et al., 2016). We performed two permutation tests: First, using all measured traits, and second using all traits associated with reproductive isolation, but only testing their overlap with putatively adaptive traits where we rejected a null hypothesis of neutral divergence in a trait-based test of selection (Fraser, 2020).

In genomic regions where QTL overlap, we examined QTL effects to determine whether genetic correlations constrain or facilitate species divergence. If overlapping QTL all have effects in the direction of the parental divergence, these genetic correlations may have facilitated divergence. On the other hand, if some overlapping QTL effects are in the direction of parental divergence, while others are in the opposite direction, these correlations may have constrained divergence between this species pair. To examine whether QTL tended to co-localize in regions with low recombination, we compared genetic distance to physical positions in the M. guttatus v2.0 reference genome (Hellsten et al., 2013).

Tests of directional natural selection between species

We tested whether directional natural selection could explain species divergence using two tests: a QTL sign test of equal effects (QTLST-EE; Orr, 1998b) and the v-test (Fraser, 2020). The QTLST-EE assumes that if a trait has a continuous history of directional selection, allelic effects at QTL will largely be in the same direction as species divergence for a given trait (Orr, 1998b). Since at least six QTL must be detected to reject the null hypothesis of neutral divergence, we were only able to apply this test for plant height and leaf shape. The v-test assumes that if a trait has a continuous history of directional selection, the variance between the parental lines is expected to be larger than the segregating variance in the F2-hybrids. To calculate our test statistic, v, we estimated the between-parent variance using a one-way ANOVA, calculated the variance within each parental line and the F2s, and used the corrected Equation 2 in Fraser (2020). We used a conservative value of c = 2, indicating full additivity of QTL.

Identifying candidate genes

We used the base pair positions of markers at the start and end of 1.5 LOD intervals for each QTL to identify annotated genes in the Mimulus guttatus v2.0 reference genome (Hellsten et al., 2013). We then used The Arabidopsis Information Resource (TAIR) to examine GO annotations for genes in our QTL regions, scanning for terms related to flowering time, growth hormones, organ growth, size, and shape (Berardini et al., 2015).

Results

Parental species differed in traits associated with serpentine endemism and reproductive isolation

Compared with M. guttatus, M. nudatus flowered at significantly earlier nodes and had longer first internodes (putatively adaptive traits), narrower and shorter corollas (traits associated with reproductive isolation), and smaller and narrower leaves (putatively adaptive traits) (Table 1, Figure 2). However, they did not significantly differ in days to flowering (a trait associated with reproductive isolation) or height (a putatively adaptive trait). F1 hybrids were typically more similar phenotypically to the M. guttatus parent (Table 1, Figure 2). Trait heritabilities ranged from 0.23 for leaf area to 0.78 for the first internode length (Table 1).

Traits associated with serpentine endemism and reproductive isolation are correlated in the F2 hybrids

Many traits measured in the F2 hybrids were significantly correlated, but correlation coefficients were not very high (Figure 3, Pearson correlation coefficients |r| < .7, per Dormann et al., 2013). The highest correlations were between node of first flower and days to flowering (r = .6), node of first flower and first internode length (r = −.53), and corolla width and length (r = .53). Except for leaf shape, which was only correlated with height (r = −.12), most traits were correlated with multiple (5–6) traits, suggesting a shared genetic basis.

Individual QTL explain a low percentage of the segregating F2 variance but a large percentage of parental divergence

Our genetic map consisted of 14 chromosomes with a total map length of 1,017.5 cM (range 43.4–125.8 cM per chromosome) and average spacing of 1 cM. Roughly 22% (222/1,011) of markers significantly deviated at an α of 0.05 from Mendelian ratios (1:2:1), particularly on chromosome 6 (Supplementary Figure S1). Almost all distorted markers showed a deficiency of M. nudatus alleles relative to M. guttatus alleles.

We detected a total of 26 QTL across all traits, ranging from one QTL for days to flowering to eight QTL for plant height (Table 2). The average percent variance explained by an individual locus was 7.6% (range: 3%–36%), while the average total variance explained by all significant loci was 24.6% (range: 10%–40%) for each trait. We did not detect significant epistatic interactions for any QTL. Five QTL were additive, 14 QTL were partially dominant, 3 were dominant, and 4 QTL were overdominant. Of the 17 QTL that were partially dominant or dominant, 14 had an M. guttatus dominant allele, while only 3 had an M. nudatus dominant allele.

For most traits, the sum of the variance explained by all loci was less than half of the trait heritabilites (average: 41%, range: 24%–57%), suggesting that additional loci with small effects contribute to each of our measured traits. For example, days to flowering had a broad sense heritability of 0.41 (Table 1) and a single QTL explained ~10% of the segregating F2 variance (Table 3), leaving 76% of the heritability unexplained. On the other hand, many loci individually explained a large percentage of the parental divergence, suggesting that few substitutions between species are necessary to explain dramatic phenotypic differences between them (Table 3, Figure 4). For traits where the parental lines differed, QTL effects were typically in the direction of the parental divergence (i.e., the percentage of the parental divergence explained was positive for 13 of 17 QTL in Table 3). The percentage of the parental divergence explained by loci with effects in the direction of the parental divergence averaged 59% (range: 9%–267%). The sum of these loci exceeded the parental divergence in node of first flower and first internode length (381% and 111%, respectively) and explained most in corolla length (96%), corolla width (66%), leaf shape (88%), and a modest percentage of leaf area (24%).

QTL contributing to traits associated with serpentine endemism and reproductive isolation overlapped

As predicted from the phenotypic correlations in the F2-hybrids (Figure 3), loci contributing to multiple traits had overlapping 1.5-LOD intervals in five genomic regions (Table 2, Figure 4). Loci contributing to putatively adaptive traits (traits that repeatedly evolve in serpentine endemics and/or potentially contribute to drought avoidance but not reproductive isolation in sympatry: xeromorphic leaves, short stature, and decreased allocation to vegetative growth) and traits associated with reproductive isolation (flowering time and flower size) overlapped more than expected by chance (overlaps: 11, Z-score = 4.09, p-value = .001). When we constrained our overlap test to putatively adaptive traits where we rejected the null hypothesis of neutral divergence (via the v-test) and traits associated with reproductive isolation, we again find that these loci overlapped more than expected by chance (overlaps: 5, Z-score = 3.00, p-value = .003).

Four genomic regions had overlapping loci contributing to two to three traits, while one region had loci contributing to five traits. QTL overlap was often, but not always, between putatively adaptive traits and traits associated with reproductive isolation. For example, QTL contributing to corolla width (a trait associated with reproductive isolation), and leaf area (a putatively adaptive trait) overlapped on chromosome 5 (the corolla width QTL also overlaps with a height QTL, but the height QTL does not overlap with the leaf area QTL on chromosome 5; Table 2, Figure 4). QTL contributing to height (a putatively adaptive trait) and corolla width overlapped on one arm of chromosome 8, QTL contributing to height, leaf shape (a putatively adaptive trait), and corolla width overlapped on a separate arm of chromosome 8, and QTL contributing to first internode length and node of first flower (putatively adaptive traits) overlapped on chromosome 10 (Table 2, Figure 4). One region on chromosome 11 contained 5 overlapping loci contributing to putatively adaptive traits (height, node of first flower, and first internode length) and traits associated with reproductive isolation (flowering time, corolla length, and width, Figure 4). Most overlapping QTL were in regions with high recombination (Supplementary Figure S3), except for the overlapping corolla width and length QTL on chromosome 11, which spanned a region with relatively low recombination (large physical distance/low genetic distance, Supplementary Figure S4) and partially inverted marker orders (Supplementary Figures S2 and S4).

Allelic effects at overlapping QTL were mostly in the direction of species divergence

In the five genomic regions where QTL overlapped, QTL effects were largely in the direction of species divergence, suggesting that genetic correlations facilitated divergence in this species pair (Table 3). In regions where multiple QTL overlap, individual loci (n = 15) on average explained 9% of the segregating F2 variance (65% of the parental difference), while loci outside of these overlapping regions (n = 11) on average explained 6% of the segregating F2 variance (|27|% of the parental difference). Loci with effects in the opposite direction of the parental divergence (n = 4) did not overlap with any other QTL. In the overlapping region on chromosome 10, F2s homozygous for M. guttatus alleles flowered at later nodes and had shorter first internodes. In the overlapping region on chromosome 5, F2s homozygous for M. guttatus alleles were taller, had wider flowers, and larger leaves. Chromosome 8 had two overlapping regions on opposite ends of the chromosome. In the first overlapping region on chromosome 8, allelic effects differed relative to the direction as species divergence: F2s homozygous for M. guttatus alleles were shorter and had wider flowers. In the second overlapping region on chromosome 8, allelic effects were all in the direction of species divergence: F2s homozygous for M. guttatus alleles were taller, had wider corollas, and had rounder leaves. Finally, in the overlapping region on chromosome 11, F2s homozygous for M. guttatus alleles flowered later, flowered at later nodes, had shorter first internodes, wider and longer corollas.

Floral and leaf trait divergence is likely due to directional selection

We were only able to apply the QTL sign test with equal effects on traits where we identified more than six QTL: plant height and leaf shape. For plant height, allelic effects at four of eight QTL were in the direction of species divergence. For leaf shape, allelic effects at three of six QTL were in the direction of species divergence. We could not reject the null hypothesis of neutral divergence between species for these two traits (p > .05).

We were not constrained by the number of QTL for the v-test and found support for a history of directional selection on the length of the first internode, corolla width, corolla length, leaf area, and leaf shape (Bonferroni adjusted p < .05; Supplementary Table S1). We could not reject our null hypothesis of neutral divergence for days to flowering, node of first flower, and plant height (Bonferroni adjusted p > .05; Supplementary Table S1).

QTL contain homologs of genes regulating flowering time, flower size, and leaf traits

The 1.5 LOD regions for QTL we identified contain hundreds of annotated genes in the M. guttatus v2 reference genome (average: 403, range: 102–1,065, Supplementary Table S2). These regions contained homologs of Arabidopsis thaliana genes with known functions in regulating flowering time, growth, flower size, and leaf size and shape (Supplementary Table S2).

Discussion

Edaphic variation generates strong divergent selection that contributes to local adaptation and the evolution of new species. Edaphic variation may be more likely to contribute to the formation of new species if traits associated with edaphic adaptation simultaneously contribute to reproductive isolation or if these traits are genetically correlated. Here we found that differences in putatively adaptive traits and traits associated with reproductive isolation between a serpentine soil endemic, M. nudatus, and its closely related widespread congener, M. guttatus, are genetically correlated through linkage or pleiotropy, and overlap more than expected by chance (Figure 4). Our two effect size estimates, the percent of F2 trait variance explained and the percent of the parental trait divergence explained yield complementary information about the genetic architecture of adaptive divergence (Table 3). In terms of the percent variance explained, few large effect loci contributed to divergence between species in a subset of traits, but many small effect loci contributed to divergence in all traits. In terms of the percent of the parental divergence explained, few loci had intermediate to large effects on all traits, but parental divergence was not fully explained for most traits, suggesting that additional small effect loci contribute. The distribution of effect sizes we observed may reflect the initial discrete, followed by gradual environmental variation experienced by populations after adaptation to serpentine soils (MacNair, 1983; Selby & Willis, 2018), or the initially large and subsequently decreasing influence of gene flow on the genetic architecture of divergence as populations adapt to serpentine and evolve reproductive isolation (Griswold, 2006; Yeaman & Whitlock, 2011). Our analyses of selection supported the hypothesis that a history of directional selection likely contributed to species divergence in traits associated with reproductive isolation (flower size) and traits contributing to drought avoidance (early bolting, leaf size, and leaf shape) (Supplementary Table S1). Alleles from serpentine-endemic M. nudatus were often recessive, which could reflect adaptive loss of function alleles in M. nudatus (Monroe et al., 2021), or may be related to M. nudatus’ demographic history, since speciation likely began with the initial colonization of small population on serpentine soils (Slatkin, 1996).

Major and minor QTL contribute to traits associated with serpentine endemism

Many soil endemics, including M. nudatus, grow readily on surrounding nontoxic soils, suggesting that specialization is not due to unique abiotic requirements for growth (Baskin & Baskin, 1988; Krukeberg, 1954). Thus, biotic factors, like competition and herbivory, are hypothesized to restrict their distributions. The repeated evolution of similar traits, like xeromorphic (small, narrow, or succulent) leaves and short statures, in the floras of unique soil types suggest that these traits are adaptive (Escudero et al., 2015; Krukeberg, 1954), and may partially explain endemism, as smaller sizes and higher light requirements may reduce competitive abilities off serpentine soils (Baskin & Baskin, 1988; Krukeberg, 1954). Our analyses suggest that a history of consistent directional natural selection is likely responsible for species divergence in leaf area, although we could not reject neutrality for plant height (Supplementary Table S1). Our tests of selection on leaf shape provided conflicting results—while the v-test suggested that divergence is due to a history of directional selection (Supplementary Table S1), the QTLST-EE did not reject the null hypothesis of neutral divergence. However, the sign test is conservative; all six leaf shape QTL would need to be in the direction of species divergence to reject the null hypothesis. Furthermore, given the remaining unexplained heritability of leaf shape (H² = 0.72, percent variance explained by QTL = 32%), we know that additional loci with small effects contribute that could not be included in the sign test.

Mimulus nudatus’ small size might also contribute to its endemism, especially if its small size reduces its competitive ability off serpentine soils. Furthermore, the evolution of serpentine tolerance could directly contribute to endemism in M. nudatus if alleles that increase serpentine tolerance decrease competitive ability off or on less harsh serpentine soils. We identified eight small to large effect QTL contributing to plant height, the largest of which overlapped with a genomic region associated with serpentine tolerance between locally adapted populations of M. guttatus (Selby & Willis, 2018). However, allelic effects at this locus are in the opposite direction of species divergence; the M. nudatus homozygote at the chromosome 13 QTL is taller at flowering than the M. guttatus homozygote. Given that our hybrids were grown in individual pots without above or below-ground competition, and that our parental lines did not significantly differ in height, it is likely that allelic effects at loci influencing plant size differ in natural environments or that different loci may contribute altogether (Mutic & Wolf, 2007; Weinig et al., 2002, 2003). In addition, our parental lines were both from serpentine-tolerant populations, and height differences are greater between M. nudatus and off-serpentine populations of M. guttatus and may be more pronounced later in the life cycle (Selby & Willis, 2018).

Serpentine soils are not only chemically challenging but also physically challenging due to their shallowness, rockiness, and low water-holding capacity. It is perhaps not surprising that plants that live on serpentine have repeatedly evolved xeromorphic foliage. Consistent with this pattern, M. nudatus’ has very small narrow leaves, a unique trait in the M. guttatus species complex. Three of the leaf area and shape loci we identified in our current study, on chromosomes 4, 5, and 8, overlap with loci contributing to leaf shape differences between a lobed leaf granite endemic, M. laciniatus, and round-leaved M. guttatus (Ferris et al., 2015, 2017). Allelic effects at the chromosome 5 and 8 QTL were in the same direction as parental divergence (M. nudatus homozygotes had smaller and narrower leaves), while allelic effects at the chromosome 4 QTL were in the opposite direction of parental divergence. Although these regions are associated with different leaf phenotypes (lobed vs. small vs. narrow leaves), QTL overlap suggests that these regions contain genes involved in the repeated evolution of leaf size and shape in edaphic endemics. For example, our chromosome 5 leaf area QTL contains a homolog of CYCLIN D1;1, a cyclin-dependent kinase that regulates leaf cell division and size (Cho et al., 2004). Our chromosome 8 leaf shape QTL contains a homolog of TCP3 (Migut.H01974), a transcription factor that regulates heterochronic leaf differentiation (Efroni et al., 2008), and a homolog of ROT3, a cytochrome P-450 enzyme that regulates polar elongation of leaf cells (Kim et al., 1998).

Major and minor QTL contribute to prezygotic reproductive isolation and putatively adaptive divergence in allocation to vegetative growth

Although flowering time differences between M. nudatus and M. guttatus have been observed in the greenhouse (Sianta & Kay, 2021) and in the field (Selby et al., 2014), our parental lines did not significantly differ in our greenhouse common garden. However, our parental line averages were within the range of previously observed phenotypic variation in each species. In a previously published study using nearby populations of M. nudatus (CSH) and M. guttatus (CSS), M. nudatus did not flower earlier than other members of the M. guttatus species complex under long (16-hr) days (Friedman & Willis, 2013). Future genetic mapping experiments under short day lengths may identify different genetic regions underlying species differences in phenology, especially since M. nudatus has a shorter critical photoperiod requirement for flowering compared to sympatric populations of M. guttatus (8 vs. 12 hr of light per day required for flowering, Friedman & Willis, 2013). Flowering time under short day lengths may be more reflective of natural variation in flowering time in the field, since these populations initiate flowering in spring when day lengths are shorter (12–13 hr) than the 18-hr day lengths experienced by our experimental plants.

Despite no significant difference in flowering time, we observed a significant difference in the developmental timing of allocation to reproductive growth between the M. nudatus and M. guttatus parental lines. Mimulus nudatus transitioned to reproductive growth at an earlier developmental stage, reflected by elongation of the first internode (bolting) and by the production of flowers at earlier nodes (Table 1). This discrepancy suggests that M. nudatus grows more slowly (produces fewer leaf pairs) than M. guttatus prior to transitioning to reproductive growth, and/or takes longer between the transition to reproductive growth and anthesis (flower opening). In our selection analysis, we found support for consistent directional natural selection on earlier bolting but could not reject the null hypothesis of neutrality for flowering time or the node of first flower. All three traits shared a major QTL on chromosome 11, suggesting that while interspecific flowering time differences are context dependent, shifts in developmental timing are genetically correlated with flowering time variation.

Flower size differences between M. nudatus and M. guttatus contribute to floral mechanical isolation between these species (Gardner & MacNair, 2000). Mimulus nudatus is primarily pollinated by small-bodied bees, whereas M. guttatus is primarily pollinated by larger-bodied bees; however, pollinators frequently transition between species (Gardner & MacNair, 2000; Koski et al., 2015). Gardner and MacNair (2000) observed that small-bodied bees transitioning from M. nudatus to M. guttatus flowers were able to avoid touching M. guttatus’ stigma while collecting pollen due to M. guttatus’ large flower size. Our selection analysis supported the hypothesis that a history of consistent divergent natural selection contributed to species divergence in corolla width and length, although the source of this selection is not known. Divergence in flower size may have been driven by selection to reduce costly interspecific hybridization (Servedio & Noor, 2003), or as a drought adaptation to serpentine soils. Given that corolla width and length had one of the highest observed correlations in the F2s (Figure 3), and we only detected a single QTL contributing to corolla length with a similar LOD peak and identical LOD interval to a corolla width QTL, it is unlikely that directional selection contributed to divergence in these traits independently. However, we cannot distinguish which trait was the direct target of directional selection with this experiment, whether correlations with a different unmeasured trait resulted in indirect selection on flower size, or whether these correlated changes are due to developmental constraints.

We found that major and minor QTL contribute to divergence in flower size between M. nudatus and M. guttatus, consistent with genetic mapping studies across the M. guttatus species complex (Ferris et al., 2017; Fishman et al., 2002; Hall et al., 2006). The corolla width QTL we identified contains promising candidate genes (Supplementary Table S2), including a homolog of BIGPETALp, a basic helix-loop-helix gene that controls petal size (Szécsi et al., 2006), a homolog of AINTEGUMENTA a transcription factor that regulates floral organ cell number and cell size (Mizukami & Fischer, 2000), and KLU a cytochrome P-450 enzyme that promotes petal growth (Anastasiou et al., 2007). All the flower size QTL we detected co-localize with previously identified loci and all allelic effects at these overlapping loci are in the direction of parental species or ecotype divergence, suggesting that the same genes might be involved in intra- and interspecific divergence in the M. guttatus species complex. The chromosome 5 QTL we identified co-localizes with an interspecific flower size QTL between Mimulus laciniatus and M. guttatus (Ferris et al., 2017) and Mimulus nasutus and M. guttatus (Fishman et al., 2002). Both chromosome 8 QTL we identified overlap with flower size QTL between annual and perennial ecotypes of M. guttatus (Hall et al., 2006) and between M. laciniatus and M. guttatus (Ferris et al., 2017). Finally, the chromosome 11 QTL we identified co-localizes with flower size QTL between M. nasutus and M. guttatus (Fishman et al., 2002) and between ecotypes of M. guttatus (Hall et al., 2006).

Genetic correlations may have facilitated adaptation and speciation

Soils that are chemically and physically harsh cover a low percentage of the earth’s surface but harbor a high percentage of endemic species (Brady et al., 2005; Escudero et al., 2015). Despite the high diversity and endemism on serpentine soils, tolerator species (species with populations on and off serpentine soils) outnumber endemics (Anacker, 2014). One reason why serpentine endemics may be less common than tolerator species is that strong divergent selection across soil boundaries can readily generate genetic divergence between populations but is unlikely to result in the evolution of reproductive isolation on its own (Bolnick & Fitzpatrick, 2007). A period of allopatry, where small populations are genetically isolated on edaphic islands, is likely necessary to complete the process of speciation. Genetic correlations between traits that increase assortative mating and fitness under divergent selection across soil boundaries may facilitate the evolution of adaptation and reproductive isolation between ecotypes (Felsenstein, 1981; Hawthorne & Via, 2001). Following a period of allopatry, genetic correlations between adaptive and reproductive isolating traits also increase the probability of coexistence by maintaining species boundaries upon secondary contact (Rundle & Nosil, 2005).

Mimulus nudatus and M. guttatus are reproductively isolated through several prezygotic isolating barriers, including microhabitat segregation (Toll & Lowry, 2022; Toll & Willis, 2018; Toll et al., 2021), flowering time differences (Sianta & Kay, 2021), flower size and pollinator differences (Gardner & MacNair, 2000), and hybrid seed inviability (Oneal et al., 2016). We found that traits associated with the evolution of serpentine endemism and reproductive isolation are genetically correlated, overlap more than expected by chance, and may have facilitated divergence and/or coexistence in this species pair. Despite similar peak positions and LOD profiles of overlapping QTL, our experiment cannot distinguish between pleiotropy and linkage (Table 2, Figure 4). However, most of the candidate genes we identified in the overlapping 1.5-LOD regions on chromosomes 5, 8, and 11 have no published pleiotropic effects on flowering time, leaf, and floral traits. Thus, QTL overlap may be due to linkage among genes influencing these traits individually.

The two genomic regions with the greatest number of overlapping QTL co-localize with loci contributing to intra- and interspecific divergence across the M. guttatus species complex (Coughlan et al., 2021; Fishman et al., 2002; Hall et al., 2006). The chromosome 8b QTL overlaps with an adaptive chromosomal inversion that contributes to floral, life history, and vegetative trait divergence between coastal perennial and annual ecotypes of M. guttatus (DIV2, Hall et al., 2010; Lowry & Willis, 2010). Mimulus nudatus and co-occurring annual populations of M. guttatus are collinear in this region, and this region is known to contribute to floral and/or life history divergence between other collinear species (Coughlan & Willis, 2019; Coughlan et al., 2021; Ferris et al., 2017). In at least one species pair, this region can be separated into two distinct QTL, suggesting that multiple adaptive alleles within this region contribute to divergence and predate the inversion (Coughlan & Willis, 2019). Additionally, this region contributes to variation in allocation to vegetative and reproductive growth within a single natural population of M. guttatus where the inversion is at very low frequency (Kelly, 2022). Like the chromosome 8b QTL, the overlapping chromosome 11 QTL identified in this study appears to be largely collinear between M. nudatus and M. guttatus (Supplementary Figure S4). However, two flower size QTL spanned a region with low recombination and partially reversed marker orders relative to M. guttatus v2.0 reference genome (Supplementary Figures S2 and S4). Since not all QTL span this region, structural variation is unlikely to explain the overlap of loci underlying traits associated with reproductive isolation and putatively adaptive divergence.

Transmission ratio distortion

We found that F2 genotype frequencies significantly deviated from Mendelian expectations and M. guttatus alleles were typically overrepresented. Transmission ratio distortion is common in genetic mapping experiments and has been observed in multiple intra- and interspecific Mimulus hybrids (Bradshaw et al., 1998; Hall & Willis, 2005; Sweigart et al., 2006), including an interspecific cross where allele frequencies were biased toward M. guttatus (Fishman et al., 2001). Transmission ratio distortion may be caused by inbreeding depression, conspecific pollen precedence, female meiotic drive, or differential zygote mortality caused by epistatic interactions among loci in hybrid genomes (Fishman et al., 2001, 2008). Although we cannot distinguish among these mechanisms in our current study, inbreeding depression may have contributed. Since our parental lines were only inbred in the greenhouse for two generations, initial heterozygosity was only decreased by ¾, so they may have still harbored recessive lethal alleles.

Mimulus nudatus alleles were more likely to be recessive

For most partially dominant and dominant QTL, the M. nudatus allele was recessive (14 of 17, Table 3). This pattern contrasts with genetic mapping studies between predominantly self-fertilizing and outcrossing species, where no directional dominance was observed (Fishman et al., 2002; Ferris et al., 2017; Kay & Surget-Groba, 2022), and is unexpected because our species have similar mixed-mating systems. Loss of function alleles are likely to be recessive and are consistent with the direction of phenotypic differences between species, as many studies describe gene knockouts with smaller organ sizes and dwarf phenotypes (Supplementary Table S2). Although we do not know what the causal genes are in our QTL, much less the causal mutations, a promising direction for future studies to explain the preponderance of recessives in M. nudatus is investigating whether candidate genes have obvious knockouts, including premature stops, radical AA substitutions, and indels in coding regions.

A non-mutually exclusive hypothesis is that founder events in M. nudatus’ demographic history resulted in the increased fixation of beneficial recessive alleles. In outcrossing populations, adaptation from new mutations is biased against the fixation of recessive alleles (Haldane, 1927), while adaptation from standing variation is unbiased with respect to dominance (except for fully dominant loci, Orr & Betancourt, 2001). The probability of fixation of beneficial recessive alleles increases with self-fertilization (Charlesworth, 1992), but M. nudatus and M. guttatus have similar mixed-mating systems (Toll et al., 2021). For a given initial allele frequency, the probability of fixation of a beneficial recessive allele decreases with decreasing population size (Kimura, 1962). However, fixation probabilities increase with initial allele frequencies and allele frequencies are initially higher in small populations provided those alleles are not eliminated by drift. When followed by a period of rapid growth, those recessive alleles are at higher frequency than they would have been if population sizes remained constant, increasing their probability of fixation (Slatkin, 1996). This demographic scenario is likely during the evolution of edaphic endemics, where small populations initially colonize new edaphic environments, followed by a period of rapid population growth (Kruckeberg & Rabinowitz, 1985). Future population genetic studies in M. nudatus may clarify whether this scenario can explain the observed patterns of dominance.

Conclusion

Physically and chemically harsh soils harbor floras characterized by narrowly endemic species. The evolution of these edaphic endemics is influenced by the genetic architecture of divergence across soil boundaries. Correlations between traits that increase edaphic adaptation and assortative mating can increase the probability of speciation and coexistence upon secondary contact. Consistent with theory about the effect size distribution of adaptive substitutions, we found that few large effect and many small effect loci contribute to divergence between the serpentine endemic, M. nudatus, and co-occurring serpentine adapted populations of its widespread close relative M. guttatus. Selection analyses support the hypothesis that consistent directional natural selection contributed to species divergence in the developmental timing of bolting, flower size, and leaf size and shape. For partially dominant and dominant QTL, alleles from the serpentine endemic were mostly recessive, potentially due to the types of substitutions fixed during divergence or founder events in its demographic history. We also found that QTL contributing to traits associated with serpentine endemism and reproductive isolation overlap more than expected by chance, either due to pleiotropy or linkage. These overlapping QTL contain many genes, and thus, future fine-mapping studies are necessary to distinguish between pleiotropy and linkage. Overall, our study suggests that genetic correlations facilitated the evolution of adaptation and speciation in this pair of species.

Data availability

Raw sequence reads are deposited in the SRA (BioProject PRJNA1032852). Individual genotype data, phenotype data, and code is available on Dryad (https://doi.org/10.5061/dryad.t1g1jwt64).

Author contributions

K.T. and J.H.W. conceived of the study. K.T. performed the QTL mapping experiment, analyzed the data, and wrote the manuscript.

Funding

This work was supported by National Science Foundation (NSF) Grants Nos. 1354688 and 1558113 awarded to J.H. Willis and an NSF Doctoral Dissertation Improvement Grant No. 1406952 was awarded to K. Toll and J.H. Willis.

Conflict of interest: The authors declare no conflict of interest.

Acknowledgments

The authors wish to thank Laryssa Baldridge, Jennifer Coughlan, Kyle Christie, Meng Chen, Kathleen Donohue, Katie Ferris, Annie Jeong, Eric LoPresti, David Lowry, Paul Magwene, Tom Mitchell-Olds, Elen Oneal, Jessica Selby, Mark Rausher, Ashley Troth, and Carrie Wessinger for advice and feedback throughout this project. We thank Natalie Knox and Brian Weill for assisting with plant phenotyping, and the Duke University greenhouse staff, particularly John Mays, for plant care. We also thank the directors of the UC McLaughlin reserve, Cathy Koehler and Paul Aigner, for facilitating field seed collection (DOI: 10.21973/N3W08D).

References

Author notes