The opposition between gene flow and selection across heterogeneous landscapes governs both local adaptation and the emergence of new species. Extrinsic barriers to gene flow (e.g., distance, dispersal barriers, or gaps in suitable habitat) allow allopatric populations to drift or adapt independently toward speciation (reviewed in Gavrilets 2003;Coyne and Orr 2004). However, a small number of migrants between populations can be enough to homogenize genetic variation genome-wide E-mail: lila.fishman@umontana.edu (Slatkin 1973;Lenormand 2002), even if environmental heterogeneity causes divergent selection. At one extreme, this interplay between selection and gene flow may result in a handful of adaptive outlier loci standing out against a background of low genome-wide differentiation (Jones et al. 2018). At the other, local adaptation itself may generate strong barriers to gene flow among habitats, extending genetic differentiation beyond the loci under divergent selection and potentially initiating speciation. Theory shows that speciation in the absence of geographic barriers is particularly possible when adaptation involves so-called "magic traits," which simultaneously confer habitat-specific fitness and enforce strongly assortative mating (Servedio et al. 2011;Smadja and Butlin 2011;Kopp et al. 2018). However, while adaptation and speciation often progress in contexts including gene flow, it remains unclear whether and how local adaptation can generate reproductive isolation and genetic differentiation in the absence of geographic barriers.
In flowering plants, soil mosaics provide an ideal context in which to jointly consider adaptation and speciation without geographic barriers (reviewed in Rajakaruna 2018). Differences in soil moisture, nutrient availability, and mineral concentrations between geologically distinct soils can impose strong divergent selection over distances within the range of seed and pollen movement (Brady et al. 2005;Richardson et al. 2014). Selection across such edaphic mosaics may result in differentiation of adaptive traits controlled by major loci, but little neutral genomic divergence, as in the classic Linanthus flower color polymorphism (Schemske and Bierzychudek 2007) or copper tolerance in monkeyflowers (Lee and Coop 2017). Nonetheless, the abundance of edaphically specialized endemic plant species suggest that soil adaptation may often generate divergence in complex traits that cause reproductive isolation (reviewed in Baack et al. 2015;Rajakaruna 2018). Occupation of different environments may directly influence reproductive traits (Levin 2009;2010), contributing to isolation by environment (IBE; Wang and Bradburd 2014;Sexton et al. 2014) even without adaptation. However, stressful edaphic patches may often favor local shifts in flowering time (e.g., Dittmar and Schemske 2017;Sianta and Kay 2021) and/or mating system (Ferris et al. 2014). In particular, the repeated evolution of the floral "selfing syndrome" (Sicard and Lenhard 2011;Opedal 2018), often in marginal habits (Baker 1955), profoundly affects patterns of genomic variation in plant populations (Schoen and Brown 1991;Charlesworth and Wright 2001;Barrett et al. 2014). Shifts in reproductive traits are thought to be a major drivers of plant speciation (e.g., Barrett and Harder 1996), and phenological and mating system evolution associated with adaptation to soil mosaics have even been implicated in putative cases of sympatric speciation (Savolainen et al. 2006) and reinforcement (Antonovics 2006). Thus, investigating how traits, habitats, and geography interact to shape genome-wide patterns of plant differentiation across edaphic mosaics is central to understanding adaptation and speciation.
We investigate how life-history and mating system divergence influence population structure and individual genomic variation in yellow monkeyflowers across geothermal soil mosaic in Yellowstone National Park (YNP). The yellow monkeyflowers of the Mimulus guttatus species complex are a model system for understanding the evolutionary genetics of edaphic adaptation (Lee and Coop 2017;Selby and Willis 2018) and life history (Lowry and Willis 2010;Coughlan et al. 2021), phenology (Friedman and Willis 2013;Fishman et al. 2014;Kooyers et al. 2019), and mating system (Fishman et al. 2002) divergence, cre-ating a robust comparative framework. Across a latitudinal range from Alaska to Baja California, M. guttatus populations exhibit alternative annual and perennial ecotypes, with annuals found in ephemeral wet soils and perennials wherever soils remain saturated through the summer (Lowry and Willis 2010;Oneal et al. 2014;Twyford and Friedman 2015;Peterson et al. 2020). Numerous traits associated with life history (including flowering time, flower and plant size, and stolon production) map to a diagnostic chromosomal inversion on Chromosome 8 of M. guttatus (DIV1; Hall et al. 2006;Lowry and Willis 2010), and widespread annuals and perennials exhibits strongly elevated differentiation across this large region of suppressed recombination (Oneal et al. 2014;Twyford and Friedman 2015;Gould, et al. 2017). Within annual M. guttatus, several distinct species have arisen via the evolution of obligate selfing (Fenster and Ritland 1994;Brandvain et al. 2014). The most widespread, M. nasutus, occurs in sympatry with M. guttatus across a broad range (Martin and Willis 2007;Kenney and Sweigart 2016), while others are endemic to distinct edaphic substrates (MacNair and Cumbes 1989;Ferris et al. 2014). Edaphic adaptation (e.g., to serpentine and copper mine soils) has also occurred within M. guttatus without mating system divergence and/or speciation (MacNair 1983;Lee and Coop 2017;Selby and Willis 2018). Thus, edaphic, life-history, and mating system adaptation have distinct outcomes even within this one species complex.
In YNP (home to >50% of Earth's geysers), geothermal inputs recapitulate on a microgeographic scale the soil conditions that shape adaptation across the range of M. guttatus. Yellow monkeyflowers are one of a few vascular plants that live in extreme geothermal crusts, which create a ground-level oasis when cooled by winter/spring snowmelt but become deadly hot (45°C-65°C) and dry by early summer (Lekberg et al. 2012). Throughout YNP, summer-flowering M. guttatus also occur in nonthermal wetland habitats (bogs, riversides) that are under snow much of the year, and in diverse intermediate habitats near geothermally influenced water sources (Fig. 1A). Previous work focused on two sites, Agrostis Headquarters (AHQ) and Rabbit Creek (RBC), which share similar thermal and nonthermal extremes but strongly differ in habitat-associated population genetic structure (Lekberg et al. 2012;Hendrick et al. 2016). In common gardens, plants from focal thermal sites exhibited shifts in floral morphology, flowering cues, and growth form consistent with adaptation to the narrow spatial/temporal window for growth and reproduction in this extreme habitat (but life-history traits were not measured; Lekberg et al. 2012). Complex microgeographic adaptation in YNP creates an ideal opportunity to investigate the buildup of genetic and phenotypic differentiation in the absence of external barriers to gene flow.
In this study, we characterize heritable phenotypic variation, genome-wide population differentiation, and genetic diversity of plants sampled across the complex matrix of M. guttatus habitat in YNP. We defined each collection site as either an annual or perennial habitat (AH and PH, respectively) based on summer dry-down and grew progeny from wild plants in a greenhouse common garden to extract heritable variation for life-history, phenological, and floral traits. We used a reduced representation sequencing protocol (ddRADSeq) to genotype the wild individuals at thousands of loci genome-wide, then examined genomic variation in the context of habitat, geography, and phenotype. First, we assessed whether plants from AH sites have evolved annuality and characterize the full spectrum of phenotypes from AH and PH habitats. Second, we characterize patterns of genetic variation across the landscape, assign individuals to population clusters, and ask how variation is structured by environment (AH vs. PH) and geography. Third, we ask how habitat affects isolationby-distance and F IS (inbreeding coefficient) and whether heritable phenotypes predict F IS within each habitat. Finally, we examine patterns of AH-PH differentiation genome-wide to test whether life-history adaptation in Yellowstone thermal-soil habitats reuses widespread annual-perennial chromosomal polymorphisms or has an independent genetic basis.
Mimulus guttatus Fisch. Ex. DC (Phrymaceae) is an herbaceous plant common in diverse wet soil habitats throughout western North America. Here, we focus on M. guttatus variation within YNP, on the far eastern edge of the species range. In YNP, M. guttatus plants are abundant across a complex edaphic mosaic generated by geothermal activity warming soils in winter and making them inhospitable in summer (Fig. 1A). Previous work characterized genetic and phenotypic variation at just two focal sites with adjacent thermal and nonthermal extremes, AHQ and RC (Lekberg et al. 2012;Hendrick et al. 2016). Here, we sampled M. guttatus in diverse nonthermal and thermal habitats in all major geyser basins except Mammoth and Mud Volcano, where high and low pH, respectively, may inhibit M. guttatus occurrence.
Mimulus guttatus individuals and habitat data were collected from >170 geo-referenced sites across YNP's geyser basins (n = 2-3 plants ∼1 m apart per site, Fig. 1B, Table S1).
We visited each site twice, timed as logistically feasible due to Park restrictions on research activity in thermal areas. The first visit (mostly in 2017) was as early in the site-specific growing season as possible to collect adult plants for transplanting into a greenhouse (see below). We strove to sample individuals from the full range of thermal and nonthermal habitats within each geyser basin. In August 2018, we revisited and categorized each collection location (at same 1 m scale) as annual habitat (AH) or perennial habitat (PH) based on soil moisture and the presence of live M. guttatus: AH = dry with dead plants; PH = wet with live adult plants. Although this binary categorization does not capture the complexity and diversity of environments generated by geothermal activity in YNP, it encapsulates the key factor selecting on M. guttatus life history there and range-wide.
To extract heritable phenotypes from our field plants under common garden conditions, we hand-pollinated one wild plant from each site to generate inbred full sibships (N = 174 sites; 154 from 2017 transplants, plus 20 grown from previous AHQ and RC wild seed collections). Seeds were stratified in 24-well microplates on wet sand for two weeks (4 °C, dark), then germinated at 24°C/12°C (16/8 h day/night cycle). Seedlings were transplanted to 2 2 pots (two pots per maternal family, randomized positions) and maintained under ∼30°C/15°C 16/8 h cycle in a greenhouse with supplemental lighting, daily bottom-watering to saturation, and monthly fertilization. When the first flower on each plant opened, we recorded days to flower, the first flower height (mm) and node number (plus internode length calculated as the first flower height/node number), and its corolla width, corolla length, and stigma-anther distance (as defined in Fishman et al. 2002). We collected all four anthers of the first flower into 50 μL of aniline blue lactophenol dye and used a hemocytometer to estimate pollen viability and pollen number (Fishman and Saunders 2008). At 85 days posttransplant (∼3 weeks after the latest first flower), we counted the number of stolons and reproductive branches and measured total plant height. Stolons were defined as vegetative branches >1 cm long arising below the second leaf pair and remaining horizontal at least to the pot edge, whereas reproductive branches were defined as branches (>1 cm long) arising from above the second leaf pair and bearing floral buds.
To test whether plants from different habitats (AH vs. PH) differed, on average, YNP-wide and to estimate broadsense heritability (H 2 ; proportion of the total within-habitat phenotypic variance explained by family) and we used mixed model REML with habitat as a fixed effect and selfed family as a random effect nested within habitat, implemented in JMP 14 (SAS Institute 2018).
Genomic DNA was extracted from the greenhouse-grown wild plants (and the using a CTAB-chloroform protocol modified for 96-well plates (dx.doi.org/10.17504/protocols.io.bgv6jw9e). We used a double-digest restriction-site associated DNA sequencing (ddRADSeq) protocol to generate genome-wide sequence clusters (tags), following the BestRAD library preparation protocol (Ali et al. 2016), using restriction enzymes PstI and BfaI (New England Biolabs, Ipswich, MA). An in silico digestion of the M. guttatus reference genome with this pair generated highly genic and evenly dispersed tags. After enzymatic digestion, sets of 48 individual DNAs were labeled with unique in-line barcoded oligos, and each set pooled into a single tube. Pools were barcoded using NEBNext indexing oligos with a degenerate barcode, PCR amplified, and size-selected (200-500 bp fragments) using agarose gels. Libraries were prepared using NEBNext Ultra II library preparation kits for Illumina (New England BioLabs, Ipswich, MA). The libraries were paired-end sequenced (150 bp) on two 1/2 lanes of an Illumina HiSeq4000 sequencer.
Illumina reads were demultiplexed and true PCR duplicates (same restriction site, plus same degenerate barcode) removed using a custom Python script (dx.doi.org/10.17504/protocols.io.bjnbkman). Adapters were trimmed and low-quality reads removed using Trimmomatic (Bolger et al. 2014). Reads were mapped to the M. guttatus v2 reference genome (https://phytozome.jgi.doe.gov) using BWA MEM (Li and Durbin 2009) and indexed with SAMtools (Li et al. 2009). We retained 58 million reads after alignment and quality filtering. After filtering of individuals for mean coverage (>7×) and other quality metrics, we retained genomic data for 231 individuals (N = 119 AH and 112 PH) representing 165 collection sites. SNPs were called using the STACKS pipeline (Catchen et al. 2011(Catchen et al. , 2013) ) for aligned sequence data. The "Gstacks" module was used to identify loci and SNPs and to genotype individuals, and we used "Populations" to filter SNPs (maximum missing data = 0.2, minor allele frequency = 0.05) and export genotypes to vcf format. For analyses of population genetic structure and individual genetic diversity, we pruned the dataset to one random SNP per ddRAD tag across the 14 main V2 chromosomes, for a total of 14,919 variable sites. For scans of YNP-wide AH-PH differentiation, we included all variant SNPs (n = 78,479) across those same scaffolds.
Individuals were geo-referenced and named by geyser basin (Fig. 1B), site number within geyser basin, and plant number within site (Table S1), but not assigned to populations a priori. To explore broad patterns of genetic structure and establish populations for downstream analysis, we used discriminant analysis
of principal components (DAPC) implemented in the R package adegenet (Jombart 2008). We used the find.clusters function on the first 20 PCs to identify the number of clusters (k) that minimized the Bayesian information content and then used the dapc function to cluster individuals (retained PCs = 7, retained axes = 4). To visualize the correlations underlying the DAPC, we performed genetic principal components analysis (PCA) using plink2.0 (Purcell et al. 2007). We also recalculated the PCA removing all but one individual from the extreme AHQT population. Finally, we used the Bayesian clustering method conStruct, which includes continuous geographic distance in model generation, to infer population structure and individual admixture proportions (Bradburd et al. 2018). The threefold cross-validation method in conStruct was used to evaluate k = 2-8 layers, with 3000 Markov-chain Monte Carlo generations assessed for independent convergence for each value of k.
To assess the hierarchical partitioning of genetic variation, we stratified individuals by regions and habitat (excluding AHQT) and conducted analysis of molecular variance (AMOVA) using the R package poppR (Kamvar et al. 2014). Significance of region and habitat effects were tested with the function randtest in package ade4 (Dray and Dufour 2007) using 999 permutations of the data matrix. We also calculated mean pairwise Weir and Cockerham's F ST (plus bootstrapped 95% confidence intervals) among the nine region × habitat subpopulations using the R program Hierfstat (Goudet 2005).
To investigate the landscape scale of genetic differentiation among individuals, we calculated pairwise Euclidean genetic and geographic distance matrices among all individuals (excluding AHQT) using the "vegan" package in R (Oksanen et al. 2017). We did not implement Mantel tests (Matesanz et al. 2011) or use Bayesian IBE approaches built on habitat dissimilarity matrices (Bradburd et al. 2013) to estimate the effects of environmental distance versus geographical distance on genetic distance. AMOVA indicated that relatively little genetic variation was explained by habitat category and replication of habitats in multiple regions (vs. a single joint gradient) is better captured by a nested model (Sexton et al. 2014). Furthermore, we expect habitat (and associated traits) to affect patterns of genetic distance within habitats via mating system, as well as to influence gene flow across habitat boundaries. To evaluate this interaction directly, we regressed genetic distance versus geographical distance for the three categories of habitat pair (AH-AH, AH-PH, PH-PH), then tested for effects of habitat on genetic distance in four geographic distance categories (0-1 km, 1-3 km, 3-10 km, >10 km) using ANOVA with habitat-pair, distance category, and their two-way interaction using JMP14 (SAS Institute 2018). Unequal distance categories were chosen because we expect the strongest effects of reproductive ecology at within-patch scales (<1 km), collections were not continuously distributed, and the relationship between geographic and genetic distance was nonlinear.
To examine genome-wide patterns of diversity, we calculated pairwise nucleotide diversity (π) in vcftools (Danecek et al. 2011), for all genotyped individuals, regional populations, and the nine region × habitat subpopulations using variant and invariant sites at 23,545 ddRAD tags (2,409,908 sites) that passed coverage and missing data filters. 95% confidence intervals for mean π values were obtained by bootstrapping.
To infer individual inbreeding history, we calculated inbreeding coefficient (F IS ) using the -het command in vcftools (Danecek et al. 2011), with expected homozygosity calculated within the five regional populations. We compared mean F IS of AH and PH YNP-wide with ANOVA, as well as among regions and all AH-PH partitions within regions. Three individuals with highly negative (<-0.20) values of F IS were excluded from this and the following analyses as potentially contaminated or recent hybrids between distinct populations; this exclusion did not materially affect the results or interpretation. To test for associations between heritable phenotypes (progeny) and inbreeding history (wild plants) within habitat, we excluded all but one AHQT individual, split the dataset by habitat (AH vs. PH) and conducted forward-backward model selection (AICc difference > 2 criterion) on all 14 phenotypes as predictors of F IS (traits and F IS standardized to mean and variance = 1). We then examined effects in the best-fit models for each habitat using a linear regression framework. F IS statistical analyses were performed in JMP14 (SAS Institute 2018).
Although dense at >50 tags/Mb or ∼10/cM (∼100 cM and 15-25 Mb per chromosome; Flagel et al. 2019), our ddRAD dataset does not provide resolution to characterize the genomic basis of geothermal adaptation (Lowry et al. 2016). However, genome scans can rule in or out the involvement of widespread inversions diagnostic of life-history divergence across the species range. Two genomic rearrangements, inv8 (aka DIV1) on Chromosome 8 and inv5 on Chromosome 5 exhibit elevated F ST in broad-scale annual-perennial comparisons (Oneal et al. 2014;Twyford and Friedman 2015;Gould et al. 2017). inv8/DIV1 (but not inv5) is also a major QTL for numerous traits defining the alternative lifehistory strategies (Hall et al. 2006;Lowry and Willis 2010). To ask whether either inversion was differentiated between AH and PH plants in YNP, we calculated mean F ST in nonoverlapping 50-SNP windows across the all-SNP dataset. We then compared mean F ST in windows from genomic scaffolds in inv8 and inv5 as defined in Gould et al. (2017) to the remainder of the genome using ANOVA in JMP14 (SAS Institute 2018). Consistent with the selection for autogamous selfpollination and reduced pollinator attraction in annual habitats, AH plants had ∼50% lower stigma-anther separation and slightly narrower corollas than PH plants (Fig. 1, Table 1). Parallel reductions in corolla length and width-length ratio were only marginally significant (P = 0.05 and 0.09, respectively). The pollen number did not differ among habitats or families, but the pollen viability was moderately heritable (H 2 = 0.38; Table 1) and AH plants had, on average, higher pollen viability than PH plants (0.89 vs. 0.83, P < 0.005; Table 1).
Both DAPC and conStruct primarily clustered individuals by geography (geyser basin) rather than habitat (Fig. 2). In the DAPC, values of k (clusters) from five to nine minimized variance among clusters (Fig. 2B, inset), so we conservatively assigned individuals to five populations (Fig. 2). One consisted entirely of plants from an extreme thermal population AHQT, which was also isolated in previous analyses (Lekberg et al. 2012;Hendrick et al. 2016). The other four clusters correspond to broad geographic regions (designated Central, Northern, Southern, and Southeastern; Fig. 2A-C). AHQT's distinctness does not reflect novel variation; a genomic PCA including only a single AHQT sample places it with Central plants (Fig. S1) and AHQT and the large Central population had similarly low numbers of private alleles per sample (1.3 and 1.1, respectively). conStruct also resolved primarily geographical structure (Fig. 2D), and models including a spatial component within layers were always better supported (Fig. S2). At the k = 5 supported by cross-validation (Fig. S2), only AHQT individuals were >85% assigned to a single layer and individuals within each DAPC-defined region showed substantial admixture (Fig. 2D). Furthermore, at k = 2-8 layers, habitat never defined a discrete component (Fig. S3).
The four DAPC-defined regional populations were moderately differentiated from each other (F ST = 0.10 -0.18) and AHQT was highly differentiated from all other populations 2). Geographical region (AHQT excluded) explained a larger portion of the variance (populations; ∼12%) than habitat within region (subpopulations; ∼2%) in AMOVA, but within-subpopulation variation was greatest (∼28%) (all P < 0.001; Table 3). The small but significant habitat effect may reflect sampling of a subset of adaptive loci or neutral effects of the environment on gene flow (Wang and Bradburd 2014).
As expected, genetic and geographic distance were highly (and nonlinearly) related regardless of habitat (Fig. S4). However, the relationship between genetic and geographic distance varied with the habitat at local spatial scales, with distance class, habitat-pair category (AH-AH, AH-PH, PH-PH) and their interaction all significant by ANOVA (all P < 0.0001; Fig. 3). as did AH (π = 0.008-0.010) and PH subsets (π = 0.009-0.011) within each region. AHQT was about half as diverse as other partitions (π = 0.005) and relatively invariant at sites segregating YNP-wide (55% of sites polymorphic vs. 87-98% for the regional mixed populations). Despite its lower expected heterozygosity, the mean inbreeding coefficient (F IS ) at AHQT was significantly higher than all other regional populations (Table 4). AH plants overall (excluding AHQT) also had significantly higher F IS (mean = 0.41 ± 0.02 SE, N = 103) relative to PH individuals (mean = 0.29 ± 0.02 SE, N = 110; habitat r 2 = 0.065, P < 0.0001; Fig. 4; Table 4). This pattern was driven by the larger majority-annual Southeastern and Central populations; in the majority-perennial Northern and Southern regions, AH partitions did not have greater F IS (Table 4). AH plants (even excluding AHQT) exhibited bimodality of F IS values (Fig. 4) and parental inbreeding coefficients were associates with offspring traits. The AH F IS model with minimum AICc (r 2 = 0.48) included four traits (Table 5; N = 40 families with all traits measured), with higher F IS values associated with flowering closer to the ground, reduced herkogamy, relatively wide corollas, and lower pollen viability (in progeny). Among PH plants (n = 44), no phenotypic variables improved prediction of F IS over the null model.
The DIV1/inv8 inversion, which is phenotypically and genetically diagnostic of life history range-wide, did not show elevated F ST between AH and PH plants in Yellowstone (Fig. 5A). Indeed, DIV1 had slightly reduced differentiation relative to the rest of the genome (mean F ST = 0.02 vs. 0.03, Fig. 5B). However, the inv5 region exhibited substantially elevated F ST (mean = 0.09; Fig. 5B) (Gould et al. 2017), and contained 25% (10/40) of the top 2.5% outlier windows. The evolution of annuality from perenniality is common in wild and crop plants, but life-history transitions are poorly understood from both fundamental and applied perspectives (Lundgren and Des Marais 2020). Life histories are complex trait syndromes with multiple components that often show consistent correlations (reviewed in Friedman and Rubin 2015; Friedman 2020; Lundgren and Des Marais 2020). In droughtintolerant Mimulus, perennial stolons (or rhizomes) are the key for overwintering after survival through summer, and loss of stolons/rhizomes is diagnostic of annuality in multiple species (Hall and Willis 2006;Lowry and Willis 2010;Nelson et al. 2021). Yellowstone thermal AH plants were dramatically more annual than PH plants, producing fewer stolons, more reproductive branches, and a twofold higher ratio of sexual to asexual reproduction (Table 1, Fig. 1C). A strong negative correlation between stolon and flowering shoot number (Fig. 1C) may reflect a genetic tradeoff in meristem allocation (Friedman and Rubin 2015;Friedman 2020) . However, unlike many systems, greater allocation to sexual reproduction in annuals did not reflect a consistent shift to earlier flowering (i.e., a live-fast, die-young strategy).
Flowering time under summer-like conditions showed high family-level heritability and was significantly bimodal within annuals but exhibited only weak AH-PH differentiation (Table 1). This contrasts with widespread M. guttatus ecotypes, where annuals generally exhibit a "fast" strategy characterized by rapid flowering (Kooyers et al. 2015;Troth et al. 2018) and stolon number and flowering time are positively genetically correlated in hybrids (Lowry and Willis 2010;Friedman et al. 2015). High variation among AH plants may reflect experience relaxed or reversed selection on intrinsic flowering speed, as the most extreme thermal sites (e.g., AHQT) experience an extended winter-spring growing season (November-June; Lekberg et al. 2012). In addition, genetic variation in critical photoperiod cues (Lekberg et al. 2012), likely uncouples developmental speed from realized flowering time in the wild. As with other studies of edaphic adaptation where phenology in common gardens do not necessarily reflect patterns of selection (e.g., Sianta and Kay 2021), abundant intrinsic flowering time variation within AH M. guttatus underlines the complexity of phenological traits. Nonetheless, the disconnect between annuality and fast flowering in YNP suggests that common correlations among life-history components may reflect finely tuned natural selection rather than intrinsic tradeoffs (see also Lundgren and Des Marais 2020;Friedman 2020;Nelson et al. 2021).
Strong reductions in stigma-anther separation (aka herkogamy) in AH flowers (Fig. 1, Table 1) likely promote autogamous self-pollination (Barrett and Harder 1996;reviewed in Sicard and Lenhard 2011;Opedal 2018), providing reproductive assurance during spring flowering prior to pollinator emergence. Low stigma-anther separation as the first step toward routine selfing is predicted by natural selection (Fishman and Willis 2008) and experimental evolution (Bodbyl Roels and Kelly 2011) in the absence of pollinators, as well as by the high evolvability of herkogamy across plant populations (reviewed in Opedal et al. 2017). Across the M. guttatus complex, derived selfer species also exhibit early flowering in ephemeral habitats, even relative to sympatric M. guttatus annuals (Martin and Willis 2007;Ferris et al. 2014;Kenney and Sweigart 2016). Other close selfer-outcrosser pairs show similar habitat/phenology associations (e.g., Pettengill and Moeller 2012;Briscoe Runquist et al. 2014). However, unlike these derived selfer taxa, plants likely retain the option of pollinator attraction and outcrossing-even the smallest AH corollas are absolutely larger than those of many highly outcrossing annual M. guttatus (e.g., Bodbyl Roels and Kelly 2011). Two nonexclusive mechanisms, adaptive shifts in energy investment after selfing becomes routine (Charlesworth and Charlesworth 1981) or fixation of formerly deleterious recessive "small" alleles by drift once pollination attraction is unnecessary (Fishman et al. 2002), likely contribute to the slight reductions (on average) in corolla size in AH plants.
In some extreme early-drying AH sites (e.g., AHQT) plants may rarely encounter pollinators, promoting corolla reduction via both mechanisms. At others (e.g., Shoshone Geyser Basin, where AH and PH plants were co-flowering in late July), AH plants likely interact with pollinators and PH mates, and thus maintain larger corollas through selection and/or gene flow.
Correlated floral and environmental variation in YNP provides a rare opportunity to investigate the causes and consequences of mating system variation in its very early stages (Whitehead et al. 2018). Our estimates of F IS suggest abundant and structured variation-individual inbreeding coefficients were shaped by habitat (Fig. 4), presumably via diverse mechanisms, and (in AH plants) also associated with heritable phenotypes for four traits (Table 5). The negative correlation of stigma-anther separation with parental inbreeding coefficient in YNP annuals supports herkogamy's importance to autogamous self-pollination (reviewed in Opedal 2018). Similarly, the negative association of F IS with first flower height is likely a proxy for extreme phenology, as short stature adapts plants to the narrow boundary layer of winter-warm air in thermal crusts (Lekberg et al. 2012). The positive association of corolla width/length ratio with F IS is unexpected, as selection for efficient autogamy favors relatively narrow corollas (Fishman and Willis 2008), and perhaps reflects inbreeding depression of this trait in a subset of PH plants. The negative correlation between pollen viability and F IS within AH plants (Table 5) also suggests variation in genetic load, as did the higher pollen viability of experimentally inbred AH versus PH plants (Table 1). However, because AH-PH differentiation largely matches previous common gardens with outbred plants (Lekberg et al. 2012), inbreeding depression is unlikely to account for the overall phenotypic patterns.
The phenotypic diversity of thermal annual M. guttatus (Fig. 1) may in part arise from the distinct genetic basis of their life-history shift, which is independent of the widespread DIV/inv8 inversion polymorphism (Fig 4). DIV1 controls range-wide variation in numerous traits beyond stolon production (Hall et al. 2010;Lowry and Willis 2010;Lowry et al. 2019). By packaging co-favored variation into alternative haplotypes spanning 100s genes, such inversions can accelerate local adaptation in the face of gene flow (Kirkpatrick and Barton 2006). However, suppressed recombination in heterozygotes for established alternative haplotypes may also constrain adaptation to novel habitats. Although we cannot yet say why the widespread DIV1 inversion polymorphism was not used to adapt to novel annual habitats in YNP or what underlies the convergent evolution of annuality there, key life-history traits clearly evolved independently from both genetic and ecological perspectives. Unlike DIV1, inv5 has not been associated with life-history traits in QTL mapping (Hall et al. 2006;Lowry and Willis 2010;Friedman 2014;Friedman et al. 2015), but elevated F ST in YNP mirrors inland annualcoastal perennial differentiation in this region (Gould et al. 2017). The only known phenotypic association of inv5 region is with vernalization requirement in hybrids between inland and coastal perennials (Friedman and Willis 2013), suggesting it may contribute to elevational or climatic adaptation within perennials. Selection on similar phenological cues, and/or unmeasured below-ground or physiological traits, may explain YNP differentiation across inv5. Finally, because our ddRAD scan did not reveal extended elevation of F ST elsewhere in the genome (Fig. 5), correlated trait shifts in YNP are likely not coordinated by a novel large (>5 Mb, like inv5 and DIV1) inversion.
Patterns of genome-wide genetic variation in YNP M. guttatus illustrate three distinct outcomes of strong natural selection on life-history and reproductive traits across a complex habitat matrix. First, genetic variation was largely structured by geography rather than ecology, with AH-PH adaptation occurring on a background of abundant genetic exchange within each of multiple regional populations. Second, despite the general lack of habitatassociated population structure, local adaptation and direct environmental effects generate habitat-specific differences in individual inbreeding coefficient and isolation by distance. Third, one large thermal annual population (AHQT) exhibits the signatures of incipient (but still reversible) speciation. These population genetic outcomes of life-history and mating system divergence are a microcosm of those seen across the entire M. guttatus species complex, as well as broadly across plant species. Manifested over short spatial (and likely temporal) scales in this unique system, these results provide a window into evolutionary processes often confounded by geography and multiple layers of reproductive isolation.
Patterns of genomic variation across all major thermal basins demonstrate that shared AH-PH phenotypic divergence generally represents adaptation in the face of gene flow. Multiple clustering approaches (Fig. 2) and estimation of the effects of habitat within the context of geographical (Tables 2 and3), found little genome-wide differentiation between AH-PH plants. This parallels the predominantly geographical structure of annual and perennial ecotypes M. guttatus across the species range (Oneal et al. 2014;Twyford and Friedman 2015;Gould et al. 2017). Low habitat-specific differentiation also characterizes the emergence of dune-specific ecotypes in Helianthus and Senecio (Andrew et al. 2012;James et al. 2021), abundant hybridization between microhabitat-adapted serpentine Monardella species (Kay et al. 2018) and weak habitat-associated population structure in succulents from a South African edaphic mosaic (Musker et al. 2021). However, this finding contrasts with the preponderance of IBE found in a review of many systems (Sexton et al. 2014), though IBD was predominant in plants. Of course, environments (and the traits they shape) matter for gene flow in this system (as evidenced by the significant effects in Table 3 and Fig. 4); however, at the spatial scales that we sampled, geography is clearly the largest barrier to (neutral) genetic exchange.
In the absence of geographic barriers, theory predicts that major (and antagonistically pleiotropic) loci are necessary to generate and maintain phenotypic differentiation resistant to gene flow (Yeaman and Whitlock 2011). For Mendelian traits, such as pigmentation (Schemske and Bierzychudek 2001;Jones et al. 2018) or copper tolerance (Lee and Coop 2017), phenotypic sorting by strong selection in adjacent habitats may generate only low migration load. However, maintaining polygenic multitrait syndromes (such as life-history strategies) in the face of gene flow can generate high load or local maladaptation. Thus, as theory predicts (Kirkpatrick and Barton 2006), complex alternative strategies are often associated with chromosomal inversions, which suppress recombination (in heterozygotes) across large genomic regions containing many functionally diverse genes (Tigano and Friesen 2016;reviewed in Mérot et al. 2020). For example, in the dune Helianthus petiolaris system (Andrew et al. 2012), multiple chromosomal inversions appear to be key to the maintenance of locally adapted ecotypes (Huang et al. 2020). Here, we can rule out re-use of the widespread annual-perennial DIV1/inv8 inversion and do not identify similarly large outlier regions, but the genetic basis of correlated trait divergence (Fig. 1c) awaits further work.
Despite the predominance of geographic structure, reproductive adaptation to the unique temporal niche created by thermal AH habitats nonetheless leaves genomic signatures (Figs. 3 and4, Tables 4 and5). Shifts in mating system profoundly should affect patterns of individual genetic variation, as well as population genetic diversity and connectivity (reviewed in Wright et al. 2008Wright et al. , 2013;;Barrett et al. 2014;Hu 2015). Over the short term, a shift to exclusive self-fertilization distributes existing variation among increasingly homozygous individuals. However, in established plant populations or species that have evolved the floral selfing syndrome (Sicard and Lenhard 2011), genetic diversity is often further minimized by low effective population size, an ecological propensity for bottlenecks, and accentuated background selection (St Onge et al. 2011;Barrett et al. 2014;Brandvain et al. 2014;Barrett and Harder 2017;López-Villalobos and Eckert 2018). In most AH M. guttatus, we see something closer to the expectation for a recent shift in mating system. In aggregate, each AH "subpopulation" carries nearly the same genetic variation (both in identity and quantity) as nearly PH subpopulations (Table 4), but that variation is (on average) more homozygous in AH individuals (Fig. 4). Furthermore, AH plants were relatively isolated from each other and PH plants at short spatial scales (<1 km), while AH-PH and PH-PH pairs were relatively isolated at intermediate distances (Fig. 3). Habitat effects at short spatial scales may also reflect clonality in PH plants but persisted when the 10 most close genetic pairs (possible clones; 9 PH, 1 AH, all within 10 m) were removed (data not shown). Few studies have examined within-population variation in outcrossing rate (but see Takebayashi et al. 2006) and direct measures of individual outcrossing rates (e.g., Colicchio et al. 2020) will be necessary to tease apart the contributions of the temporal/spatial environment and reproductive phenotypes on patterns of mating and genetic relatedness in YNP M. guttatus. However, our finding that individual inbreeding history is phenotypically and environmentally predictable (Fig. 3, Table 5) suggests novel opportunities for investigating the early dynamics of mating system evolution.
Overall, our results indicate that a multi-trait adaptive syndrome, including floral traits promoting 100% autogamous selfpollination, has evolved and is maintained within the context of regional genetic exchange. Even the most extreme AH individuals appear to be built from essentially random subsets of local PH variation, plus alleles at a limited number of loci under selection. Given the spatial and temporal patchiness of the YNP geothermal matrix (Figs. 1 and2), AH-adapted alleles likely have an extensive history of transiting through relatively contiguous PH populations, and then reassembling (with help from selection as well as inbreeding) in distant or newly erupted thermal annual sites. This might resemble, on a vastly smaller spatial and temporal scale, the global dynamics of sticklebacks, where (to a large extent) local adaptation to newly appearing freshwater habitats occurs from standing variation persisting as rare maladaptive variants in the marine environment (Colosimo et al. 2005;Nelson and Cresko 2018). This scenario strongly contrasts with the idea that mating system and flowering phenology act as "magic
traits" whose divergence can generate sufficient reproductive isolation to initiate plant speciation in sympatry (e.g., Savolainen et al. 2006). Here, despite shifts in mating system and phenological traits sufficient to generate predictable habitat differences in individual outcrossing rate and connectivity, we do not see the build-up of reproductive isolation. Although definitions of species vary (Harrison 2012), this is the exact opposite of speciation (i.e., the maintenance of genome-wide genetic continuity despite strong natural selection on reproductively isolating traits).
The one exception to the overall pattern, AHQT, provides a rare window into an alternative evolutionary outcome of strong microgeographic selection in plants-speciation by mating system shift. Although plants from this population share a common thermal annual phenotype with many other AH plants, they are highly differentiated genome-wide from all other populations (F ST ∼ 0.2) and showed reduced π, polymorphism, and heterozygosity despite being sampled over a similar spatial scale. Indeed, the high differentiation of AHQT in this and previous studies (Lekberg et al. 2012;Hendrick et al. 2016) is comparable to values for genome-wide F ST for cryptic species of Eucalyalptus and Silene recently identified at much larger spatial scales (Binks et al. 2021;Kahl et al. 2021). High differentiation and low diversity of AHQT also parallels that seen in a (likely selfing) solftara specialist (Nagasawa et al. 2020); however, unlike that sedge and a highly selfing Agrostis grass that also colonizes our YNP thermal areas (Tercek and Whitbeck 2004), AHQT exhibits de novo mating system evolution and reproductive isolation over a scale of meters.
What makes the AHQT population so special? AHQT's genetic distinctness undoubtedly depends on the extreme ecology of this large and uniformly geothermal site, where soils reach 65°C in summer and there is no source of soil moisture other than melting snow. However, the novel trigger of AHQT's differentiation seems to be its discrete geographic setting-even the mildest thermal areas at AHQT are separated by ∼200 m gap of no M. guttatus habitat from the nearest nonthermal population, the late-flowering bog perennial AHQN (Hendrick et al. 2016). This sharp habitat boundary plus geographical gap prevents gene flow via intermediate populations, facilitating both unimpeded thermal-soil specialization and the neutral differentiation by drift typical of established selfers. A perfect storm of habitat-specific selection for reproductive assurance, drift via high selfing rates, and very low migration likely recapitulates the early stages of speciation in numerous derived selfing taxa (Barrett et al. 2014). Alone, each factor would be insufficient to fully isolate AHQT over such a short spatial scale, but all together they appear to allow the first steps toward speciation. Thus, while divergence in "magic traits" like mating system or flowering time undoubtedly contributes to the evolution of reproductive isolation in diverse geographical contexts, a boost from space (plus time) may be necessary for barriers to develop to the point that isolation can persist in full sympatry. In obligate selfer M. nasutus, which often self-pollinates in the bud and doesn't open flowers to pollinators at all, that point was likely reached before its spread into extensive range overlap with progenitor M. guttatus (Martin and Willis 2007;Kenney and Sweigart 2016). Similarly, allopatric isolation appears to have been key to the evolution of a selfing subspecies of Clarkia (Pettengill and Moeller 2012). At AHQ, speciation may continue down that path to the point of no return but would still be easily reversed by the emergence of a warm stream allowing intermediate-flowering M. guttatus to bridge the temporal/ spatial gap. Nonetheless, the current isolation of AHQT (in the context of rampant gene flow elsewhere) suggests that both geographic discontinuity and reproductive trait divergence may be necessary to effectively shut down gene flow genome-wide and lead to speciation.
The yellow monkeyflowers of YNP exhibit dramatic life history, growth form, and reproductive trait adaptation to the novel edaphic environments created by geothermal activity. This microgeographic divergence parallels differentiation between widespread annual and perennial ecotypes range-wide but has a novel genetic basis and reflects the selective pressures of the unique YNP habitat(s). Except for the extreme and isolated AHQT population, thermal annual plants maintain multi-trait adaptation in the face of homogenizing gene flow across the habitat matrix, with geography the primary factor structuring genetic variation on the landscape. Nonetheless, thermal annual individuals exhibit elevated inbreeding coefficients and isolation from neighbors, providing rare insight into the early stages of mating system evolution. Together, these findings have important implications for understanding the (parallel) evolution of trait syndromes, the structuring of genetic variation by intrinsic and extrinsic factors, and the limits to sympatric speciation.
KMK conducted most of the field, greenhouse, and lab work, as well as the bioinformatic analyses, as part of his M.S. thesis under the supervision of LF. CB performed some of the population genetic analyses and assisted with data curation. TCN developed genomics protocols and assisted with lab work and bioinformatics. LF conceived of the research, assisted with field, greenhouse, and lab work, conducted phenotypic and downstream genetic analyses, and constructed most figures. The manuscript was primarily written by LF and KMK, with contributions and editing from all authors.
EVOLUTION APRIL 2022