{'fraction_non_alphanumeric': 0.02935070043927705, 'fraction_numerical': 0.019013484275629957, 'mean_word_length': 4.8619011017166285, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-pollinator selection for a floral homeotic mutant conferring loss of nectar reward in Aquilegia coerulea

Here, we describe a polymorphic population of Aquilegia coerulea with a naturally occurring floral homeotic mutant, A. coerulea var. daileyae, where the characteristic petals with nectar spurs are replaced with a second set of sepals. Although it would be expected that this loss of pollinator reward would be disadvantageous to the mutant, we find that it has reached relatively high frequency (∼25%) and is under strong, positive selection across multiple seasons (s = 0.17-0.3) primarily due to reduced floral herbivory. We identify the underlying locus (APETALA3-3) and multiple causal loss-of-function mutations indicating an ongoing soft sweep. Elevated linkage disequilibrium around the two most common causal alleles indicates that positive selection has been occurring for many generations. Lastly, genotypic frequencies at AqAP3-3 indicate a degree of positive assortative mating by morphology. Together, these data provide both a compelling example that large-scale discontinuous morphological changes differentiating taxa can occur due to single mutations and a particularly clear example of linking genotype, phenotype, and fitness.


In brief
Cabin et al. present a case of natural selection in action. They show that a floral homeotic mutant of Aquilegia is under strong positive selection primarily due to reduced floral herbivory. Loss-offunction alleles at APETALA3-3 underlie the phenotype and molecular signatures at the locus indicate a soft sweep and assortative mating by morphology. INTRODUCTION Darwin famously argued that evolution proceeds gradually through many changes of small effect, slowly improving one form over another. 1 Early theory supported this view, 2 but later it was recognized that intermediate-effect size mutations are also likely to contribute to adaptation. 3 Most recently it has been suggested that larger sized mutations may be selected particularly when they occur early during the process of adaptation, and the identification of large-effect QTL supports this prediction. [4][5][6] However, we usually do not know whether or not these large-effect QTL are caused by single mutations or multiple, smaller-effect mutations. The most extreme mutations known are homeotic mutations where whole organs are misplaced during development and as such they are usually manifestly unfit (e.g., the replacement of antennae with legs in the Antennapedia mutant). 7 However, comparative analyses show that changes in the number or expression patterns of genes involved in organ development (e.g., homeotic patterning genes [8][9][10][11][12] or organ symmetry genes 13,14 ) are correlated with major morphological differences between taxa. Although whether these differences evolved through single macromutations or the accumulation of many smaller changes is unknown. In addition, it is unclear whether these changes occurred at the time of speciation (resulting in cladogenesis) or subsequently. In fact, we have very little evidence that single homeotic mutations actually can and do survive in nature.
In plant populations, a few examples of naturally occurring floral homeotic mutants have been documented, but they either represent very rare individuals in a population where a selective advantage could not be determined 15 or occur in predominantly selfing or clonal species where they are unlikely to promote reproductive isolation. 16,17 Perhaps the best example to date is a floral homeotic mutant in Capsella bursa-pastoris where petals are converted into stamens. 16 The mutant form is found intermixed with wild-type plants in vineyards and ruderal sites in Europe, and although the mutation maps to a single QTL, the gene(s) and mutation(s) have not been determined. 18 Although some aspects of reproduction, such as flowering time, seed set, and germination rates differ between the morphs, 19 their fitness appears to be similar and QTL analysis suggests that these fitness associations are not due to pleiotropy with the floral mutation. Instead, the homeotic shift likely arose in a lineage that had diverged in these other traits, and the associations are maintained by the high degree of selfing in the species. 18 Given that no clear example of a known single macromutation being favored in a natural population has been observed, whether homeotic mutations could account for some patterns of evolution has remained controversial. [20][21][22] Here, we present strong evidence of the success of a naturally occurring homeotic mutant in a natural population.
The columbine genus, Aquilegia, is noted for having petals with long, tubular nectar spurs. These nectar spurs have been considered as a ''key innovation'' 23 driving diversification in the genus as the length, shape, and color of these nectar spurs are correlated with the major pollinators of each species. 24 A. coerulea (''WT''; Figure 1A) follows the basic floral bauplan of the genus having five flat, petaloid sepals in the outermost whorl and five petals in the second whorl, each with a tubular outgrowth (nectar spur). Pollination has been extensively studied in this species, and there are two major pollinators: hawkmoths, which probe the long spurs for nectar, and bumblebees, which collect pollen. [25][26][27][28][29] The abundance of hawkmoths (but not bumblebees) in a population is significantly associated with higher outcrossing rates. 26 As A. coerulea suffers from severe inbreeding depression, 30 traits promoting outcrossing (e.g., hawkmoth visitation) are expected to be favored.
In 1897, a homeotic mutant where the petals are converted into a second set of sepals (10 sepals, 0 petals) was described as A. coerulea var. daileyae (''d''; Figure 1B). 31 As noted in the original description, apart from the replacement of nectar bearing petals with non-nectar bearing sepals, there appear to be no differences in floral structure between the two morphs ( Figures 1A and 1B). Despite this radical change in 2 nd whorl organ identity, it was noted at the time that the mutant form appeared to be increasing in frequency, 31 and 100 years later, the high abundance of the d morph was noted $24 km south of the original locality in central Colorado (Reynolds Park Open Space). 32 The relatively high frequency of the d mutant suggests that it is at a minimum selectively neutral or nearly so and may in fact confer a fitness advantage. This would be surprising as pollinators are usually thought to be primary agents of selection on floral form, but the loss of petals and nectar reward would likely discourage visitation by a major pollinator, hawkmoths. Furthermore, even if hawkmoths attempted to visit d plants, their long tongues would prevent their bodies from coming into close contact with the reproductive organs of the flower ( Figure 2A) and effect pollination. However, non-pollinator agents of selection, both biotic and abiotic, can also be important factors shaping floral morphology. [33][34][35] The Reynolds Park population ( Figure 1C) provided us with a unique opportunity to study the apparent success of this homeotic mutant in a naturally occurring common garden. We use a combination of field observations and molecular techniques to show relatively strong, consistent, and ongoing positive selection for a radical floral homeotic mutant and its likely causes.

RESULTS
Habitat, pollination, and resource allocation In each of three flowering seasons (2014-2016), we mapped the position of every flowering plant in one portion of Reynolds Park ( Figure 1C) and recorded their morphology, the number of flowers produced, any signs of floral herbivory and whether fruits matured. The vast majority of plants had flowers of either WT or d morphology, although rare individuals ($1%) with intermediate morphology (or variable morphology within and/or between flowers) were observed ( Figures S2A and S2B), which we included as mutants in analyses. The two morphs did not differ in their spatial distribution in the population and grew as close as 5 cm apart (Figures 1C and S1), indicating that environmental factors, such as soil chemistry, sunlight, and water availability, are unlikely to differ between the morphs. The two morphs also did not differ in flower number in any year nor fruit size, which is highly correlated with seed number (Figures 2B-2D). However, we found in each flowering season that d flowers were significantly more likely to produce fruits than WT flowers ( Figure 2E; s WT = 0.17-0.3; Table S1). We measured fruit set as the proportion of flowers that successfully set seed. The large differences in fruit set among seasons is likely the result of interannual variation in environmental factors, such as water availability and temperature (Data S5), and although the annual variation is notable, both morphs are affected in a similar manner. We observed both major pollinators (hawkmoths and bumblebees) visiting A. coerulea in Reynolds Park in all three field seasons. In 2014 we used video cameras to quantify visitation to flowers of both morphs growing close together and observed hawkmoths and bumblebees visiting both morphs and that individual pollinators often transitioned between the morphs (Video S1; Table S2). When hawkmoth visited the d morph, their body was held away from the reproductive organs in contrast to their visits to WT flowers (Figure 2A; Video S1). Thus, WT flowers are likely predominantly pollinated by both  hawkmoths and bumblebees, whereas d flowers are likely predominantly pollinated by bumblebees. Lastly, we tested for differences in pollen dispersal using four unlinked, neutral (4-fold degenerate) SNPs and found that d plants had significantly higher inbreeding coefficients (F) than WT plants (paired t test, t = 3.07, p = 0.027; Figure 2F).
For open-pollinated flowers, the first flower to open on an inflorescence is significantly more likely to set a fruit and to be larger than subsequent fruits. 27 For plants without signs of floral herbivory, the morphs did not differ in fruit size for the first flower produced, although first-flower fruit size was larger for plants producing more total flowers (ANOVA, p = 2.36 3 10 À5 ; Table S4). We further controlled for these architectural and plant size effects by comparing fruit set between morphs for plants producing the same number of flowers (1, 2, or 3) and lacking any signs of floral herbivory during each of the three flowering seasons. We found only one out of nine comparisons to be marginally significantly different between the morphs (Table S3) indicating overall similar levels of fruit set for the morphs. However, in 8 of the 9 comparisons d flowers set proportionally more fruit than WT flowers, suggesting a role for resource allocation favoring d plants.

Floral herbivory favors homeotic mutant
Across all three seasons we found significantly greater floral herbivory on WT flowers ( Figure 3A). We observed three floral herbivores-aphids, caterpillars, and deer. Aphids feed on the pedicels of flowers (and sometimes the flowers themselves), often causing wilted and stunted floral development ( Figures S2C and S2D). Caterpillars (primarily Platypolia anceps; Figures S2E and S3) preferentially grazed on reproductive organs (as well as sepals and petals; Figure S2E), usually resulting in complete loss of fruit production. Mule deer (Odocoileus hemionus) were observed to eat entire flowers ( Figure S2F; Video S2). In 2015 and 2016 we recorded damage due to each specific herbivore and found that although caterpillars showed no significant morph-preference in either year, aphids were significantly associated with WT flowers in both years and Mule deer showed a preference for WT flowers, although this was significant only in 2016 ( Figure 3B).
Multiple independent loss-of-function mutations cause the floral mutant phenotype An obvious candidate gene underlying the d morph is the B-class floral identity gene APETALA3-3 (AqAP3-3), as previous work has shown knockdown of expression of this transcription factor in A. coerulea results in a homeotic shift from petals to sepals without affecting other organs. 36 Thus, we sought to determine if likely loss-of-function AqAP3-3 mutations are associated with the d morph. We assessed sequence variation using Illumina sequencing of an approximately 3.6 kb PCR product, which included the entire coding region, nearly 1.1 kb upstream of the transcriptional start site and thus likely the entire cis-regulatory region 37 and the entire downstream intergenic region. We  ) and identified 55 variable positions, four of which were deletions likely to create loss-of-function mutations: two single bp deletions that would each cause frameshifts and a major truncation of the protein, a 9 bp deletion that causes the loss of highly conserved amino acids in the DNA-binding domain, 37 and a 682-bp deletion immediately upstream of the gene, eliminating a known cis-regulatory element ( Figure 4A). 38 Genotypes at these loci were highly correlated with flower morphology. All 32 d plants were either homozygous for one of these mutations or heterozygous for two of them. In contrast, 39 of the 49 WT plants either had none or were heterozygous for only one of these mutations. The remaining 10 WT plants appeared homozygous for one of the loss-of-function mutations but were also homozygous at every other site, suggesting that one allele might not have amplified ( Figure S4). We tested this possibility using PCR-based genotyping probes 39 for the three most common mutations and found that all 10 plants were actually heterozygous for a single loss-of-function mutation. The d morph therefore appears to be caused by fully recessive loss-of-function alleles at AqAP3-3, and we collectively refer to them as d alleles (specifically as d 1 , d 9 , d 1b , and d 682 ) and dominant functional alleles as WT ( Figure 4A). There was one other mutation of interest (a 10 bp deletion) in the last exon that would cause a frameshift, changing the last 25 amino acids and adding an additional 17 amino acids ( Figure 4B; Data S4, marked as ''frameshift''). This mutation occurred as a single heterozygote with WT morphology; thus, we could not assess its effect when homozygous or heterozygous with d alleles and was excluded from further analyses.
To further test that these loss-of-function alleles underlie the d phenotype, we genotyped all 850 flowering plants (d = 198, WT = 652) from 2016 using the genotyping probes for the d 1 , d 9 , and d 1b alleles. No more than two d alleles were found in any individual confirming complete linkage disequilibrium among these loss-of-function alleles and that individual haplotypes carry only one d mutation. We therefore assigned an overall genotype to each individual as d/d (homozygous or heterozygous), d/WT, or WT/WT. These genotypes were exceptionally strongly associated with floral morphology (p = 1.49 3 10 À153 ; Tables 1 and S5; Data S3). We note that nine out of the 10 rare plants with intermediate morphologies were genotyped as d/d and the other as d/WT (Table S5) but were not associated with any specific d allele or allelic combination. Fifteen d plants did not carry two of the assayed d alleles and therefore may carry one or two unassayed d alleles (e.g., d 682 ). The 13 WT plants genotyped as d/d could be the result of rare recombinants between d haplotypes, compensatory mutations, allelic dropout due to primer site variation, or mistakes in sample labeling or processing. In addition, there is a second copy of AqAP3-3 (AqAP3-3b) that could account for d/d plants having WT or intermediate morphs. In a horticultural line of A. coerulea, AqAP3-3b has very low expression in floral tissue, 36 although it is possible that variation producing higher expression at this locus in Reynolds Park occasionally rescues loss of function at AqAP3-3. It is also possible that epistatic interactions at other loci account for the d/d plants with WT and intermediate morphologies.
To determine whether different d mutations may occur on the same haplotype and whether any other mutations are associated with non-functional AqAP3-3 alleles, we phased the sequence data and constructed a haplotype network ( Figure 4B). 40 We identified functional haplotypes (WT) as those that when heterozygous with a haplotype carrying a d mutation produce WT morphology. No haplotype has more than a single d mutation and haplotypes with different d mutations are all derived from functional haplotypes. Haplotypes with the d 9 and d 1b mutations differ from functional WT haplotypes by just these mutations indicating that they likely cause loss of function at AqAP3-3. Haplotype 12 carries the d 1 mutation but also differs from the functional haplotype 10 at another site. However, both states of the second site occur in functional haplotypes (7 and 8) indicating that neither state affects function and that the d 1 mutation is causal for loss of function. Similarly, haplotypes 21 and 22 are also separated by the d 1 mutation and an additional variable site (triple arrowhead). This second site, which occurs upstream of the AqAP3-3 transcriptional start position but not at a known regulatory element is unlikely to cause a loss of function. Since the d 1 allele alone is causal on haplotype 12, it is likely that the d 1 mutation is also causal for phenotype on haplotypes 2, 3, and 21. We also note that the high degree of interconnectedness of the haplotype network indicates a great deal of recombination has given rise to different haplotypes. For example, the site separating haplotypes 23 and 24 is the same site separating 19 and 20; thus, this quartet of haplotypes forms a classic 4-gamete signal of recombination. Similarly, the d 1 haplotypes (2, 3, 12, and 21) are most likely the result of a single origin for the d 1 mutation followed by recombination events with other haplotypes (Data S4).
We were also able to assess whether different genotypes causing the same phenotype differed in fruit set or herbivory.

Molecular signatures of selection
Variable selection in the past or at other life stages acting directly on AqAP3-3 or closely linked genes could enhance, eliminate, or reverse the effect of positive selection we observed for the d morph. 41 If, however, directional selection has consistently favored d alleles for many generations, then we would expect to see heightened linkage disequilibrium surrounding highfrequency d alleles compared with ancestral WT alleles. 42 We assessed linkage disequilibrium using extended haplotype homozygosity (EHH) for the haplotypes identified carrying the two high-frequency d mutations, d 1 (28%) and d 9 (13%), compared with the WT haplotypes ( Figure 5). 42,43 EHH for WT haplotypes rapidly declines within 1 kb, reflecting previous findings for a rapid decline in linkage disequilibrium for AqAP3-3 and other genes in Aquilegia. 44 In contrast, both d alleles have much higher levels of EHH extending beyond 1.1 kb and 2.5 kb up-and downstream of the mutations, respectively ( Figure 5). These results indicate that both of these d alleles likely rose to high frequencies due to positive selection over many generations.
Morph-based assortative mating The pollinator observations described earlier revealed that both hawkmoths and bumblebees visit both morphs, but differential pollen placement could still be driving some amount of assortative mating ( Figure 2E). Assortative mating between morphs would specifically cause a deficit of d/WT heterozygotes above that observed for d/d heterozygotes (e.g., d 1 /d 9 ). The observed genotypes strongly departed from Hardy-Weinberg expectations (c 2 = 23.97, p = 2.5 3 10 À5 , d.f. = 3; Table 2), and standardized residuals indicated that this was largely the result of a lack of d/WT heterozygotes (but not d/d heterozygotes) along with an excess of d/d homozygotes ( Table 2). The excess of d/d homozygotes corresponds with the elevated F found for d plants ( Figure 2F).

DISCUSSION
Here, we have identified a homeotic mutant with a strong fitness advantage in a natural population, as well as the likely ecological drivers of selection, the underlying locus and multiple causative mutations. Our measurement of selection against the WT (s = 0.17-0.3) is similar to or greater than measures of strong selection in notable, classic examples of evolution, such as the pepper moth (Biston betularia) and pocket mice (Chaetodipus intermedius). 45,46 Furthermore, knowledge of the genetic underpinnings of the mutant allowed us to identify molecular signatures, suggesting that selection has been positively acting on multiple independent alleles over many generations and the presence of assortative mating by morphology. Although our direct observations were limited, we found no evidence that hawkmoths or bumblebees preferentially visit either morph. We did find a significant increase for inbreeding in the d morph, and as the morphs do not differ in their spatial structure or seed dispersal, reduced pollen dispersal for d plants likely explains this finding. This would match what is known for both hawkmoth and bumblebee pollination patterns. Multiple studies have shown that when tongues are longer than spurs pollen transfer is reduced. 24  See Table S5 and Data S3. Genotype is highly associated with phenotype at AqAp3-3.
ll OPEN ACCESS pollen movement ( Figure 2A). Alternatively, bumblebees actively collect pollen and groom themselves resulting in pollen moving shorter distances. 49,50 As bumblebees act as the main pollinator for the d morph, this behavior would result in higher inbreeding, which would most likely incur a fitness cost through inbreeding depression. 30 The fact that WT flowers have effective visitation by both hawkmoths and bumblebees makes it difficult to imagine a scenario where pollinators are likely to cause selection favoring the d morph. Given the likely differences between the morphs in pollen dispersal dynamics further studies of pollen movement in this system would be particularly interesting. Another possible selective advantage for the d morph could be through saved resources. For example, nectar production has been shown to be a significant energetic cost in some species 51,52 and can reduce seed set. 53 However, in other species, nectar costs appear to be negligible with no reproductive or vegetative costs. 54-56 Petals themselves could also be more energetically costly to produce or maintain than sepals. We found that neither flower production overall, fruit set, nor fruit size significantly differed between the morphs on plants that lacked herbivory. However, in eight of the nine comparisons of fruit set, d plants set a higher proportion of fruit than WT plants indicating a role for resource differences. We also note that these comparisons were only made within a season, and it is possible that in these perennial plants reallocation of resources occurs between seasons and affects life-time reproductive output. However, we doubt that resource allocation plays a major role as then we would expect many instances of petal or spur loss across the entire genus, which does not occur. Nonetheless, further studies of how resource allocation influences this system are warranted.
In contrast to pollination and resource allocation, floral herbivory (specifically aphids and deer) was consistently more abundant on WT flowers and thus favored the d morph. There are a number of ways that the loss of petals and/or the gain of an additional whorl of sepals could directly lead to herbivores either learning to discriminate against the d morph or having an innate preference for WT. For example, deer may prefer WT flowers that have sweet nectar in their spurs and learn to avoid d flowers. WT petals offer clear visual differences from the d morph having both large protruding spurs compared with flat sepals and contrasting white petal blades set against purple sepals compared with two whorls of identical purple sepals. Although we found a significant preference by deer for WT plants in only one year, it has previously been shown that large mammal herbivory can be inconsistent across generations and still elicit a strong evolutionary response, even cancelling out and reversing selection driven by pollinators. 33,34 Petals are also often a primary source of volatile compounds 57 ; thus, the two morphs may very well differ in olfactory cues. Aphids identify their host plant through a cocktail of chemical cues and will leave an incorrect host for an alternate one and, in some cases, refuse to feed. 58,59 Additionally, as aphids feed on phloem, lack of nectar production in the d morph may result in lower phloem flow, which could decrease aphid population growth rate and result in less damage to d flowers. Regardless, increased reproductive output due to reduced floral herbivory appears to be a major factor favoring the d morph in Reynolds Park. Future studies testing specific hypotheses on how the homeotic shift influences herbivores will be especially intriguing.
An alternative hypothesis to the homeotic shift directly affecting herbivory and fitness is that some other locus is responsible and is tightly linked to d alleles causing them to hitchhike to high frequency. This seems unlikely for a number of reasons. First, within 20 kb downstream and 30 kb upstream of AqAP3-3 there are just three identified genes ( Figure 5). One of these genes, Aqcoe5G180700, has been functionally annotated as a xylulose kinase-related protein, which has no obvious involvement with herbivory, whereas other two genes are of unknown function. Although possible, it seems unlikely that any of these genes actually cause differences in herbivory or fitness. Additionally, for the highest frequency d mutation, d 1 , linkage disequilibrium is breaking down before reaching any of these genes; therefore, we might expect a lower correlation between d 1 and herbivory or fitness caused at a linked locus as compared with other d alleles with higher linkage disequilibrium. However, we found no evidence for such differences (Data S5).  This linkage hypothesis would also require that the presently high-frequency d 1 and d 9 mutations both independently arose on haplotypes carrying favorable mutations or that favorable mutations arose separately on haplotypes carrying the d 1 and d 9 mutations. For the first scenario, selection would presumably favor the alternative haplotypes before the d mutations occurred; therefore, we would expect to find at least one WT haplotype both closely related to d 1 and d 9 haplotypes and at high frequency, but this is not the case ( Figure 4B). The second scenario requires favorable alleles arising only on separate d haplotypes when they were rare and then causing them to rise in frequency. Again, this seems unlikely. Thus, the simplest explanation of our data is that the loss of AqAP3-3 functionality, and the corresponding homeotic conversion of petals to sepals also results in increased fitness through reduced herbivory and possibly saved resources.
This study is also a clear example of a multiple origin soft sweep, which occurs when multiple beneficial alleles arise before any single allele is able to reach fixation. 60 These types of sweeps are far more likely when there is a large mutational target or the selected allele is recessive as it takes time for the initial allele to rise to high enough frequency to occur as a homozygote and selection can act, giving functionally equivalent alleles a chance to appear. 60,61 For the d morph, both of these conditions are met as many mutations can cause lossof-function alleles at AqAP3-3, and the phenotypic effect appears to be completely recessive (i.e., we find no evidence of a fitness difference between homozygous and heterozygous WT plants). Thus, the time between the original mutation and selection favoring d plants may be relatively long and allow for additional d alleles to accumulate before selection can bring them to high frequency. This would be especially true if selection does not begin to act until there are some number of d plants present in a population as would be the case for a learned response by an herbivore. Such dynamics may also limit or slow the spread of d alleles to nearby populations where the process would have to repeat. We also note that the relative frequencies of the most common d alleles also support the inference of a soft sweep. Theory predicts that a hard sweep can also produce multiple haplotypes, although in this situation, the most common haplotype carrying the adaptive mutation will be far more common than the next most common haplotype (e.g., 503) but that in a soft sweep the frequencies will be far more similar. 61 In Reynolds Park the most common d haplotype is haplotype 3 at 26% (carrying the d 1 mutation), and the second most common is haplotype 11 at 19% (carrying the d 9 mutation), further supporting the pattern of a soft sweep ( Figure 4B).
Although the establishment of d individuals is likely to be slow initially, the subsequent increase would be enhanced by assortative mating. Our data at AqAP3-3 indicate a degree of assortative mating by morph, which is likely caused by the mechanics of plant-pollinator interactions. In particular, as noted earlier, hawkmoths likely primarily cause mating among WT plants. In contrast, bumblebees have access to anthers of both morphs, but d plants may be more attractive because of greater pollen availability via reduced removal by hawkmoths. 26,28,29 Therefore, bumblebees may spend more time collecting pollen at d plants. Again, additional studies documenting pollen transfer in this population would be especially useful in clarifying the degree of assortative mating.
The changes in mating patterns suggest that this mutation could lead to a cladogenetic split from A. coerulea. This would likely require some combination of a reduction in gene flow between morphs and divergent selection. Although this is unlikely to occur within Reynolds Park where we find no evidence of divergent selection, speciation may occur if, for instance, floral herbivory varies among populations. In populations without floral herbivores WT plants likely have an advantage over d plants because WT plants likely have higher outcrossing rates (Figure 2F; Table 2). If populations become fixed for alternate morphs, pollinators (particularly hawkmoths) may strongly discriminate between them. We note that Reynolds Park is at the very edge of the species' range where selection pressures are often unusual. 62 Reynolds Park is at a lower elevation and is forested unlike the high altitude and open rocky habitats of typical A. coerulea populations. 63 Thus, this population may be experiencing unusual selective pressures, such as floral herbivory, compared with other populations.
Although we cannot know for sure whether the d mutant will ultimately result in speciation, comparative studies show multiple plant lineages are associated with the loss of nectariferous petals. Within the Ranunculaceae there are nine independently evolved apetalous lineages, 64 suggesting that selection has favored these losses despite the presumed mutualistic advantage of producing a food reward for pollinators. Similarly, the loss of the nectar spur, but not the entire petal, has occurred in a single species of Aquilegia, A. ecalcarata. 65, 66 Thus, there are multiple examples of petal loss being associated with cladogenesis. Furthermore, seven of these nine petal losses are associated with loss of function at AP3-3. 67 However, the majority of these petal losses are not due to a replacement with sepals and thus may evolve through a different pathway. 68 Regardless, these examples show that loss of nectariferous petals can be an evolutionary stable and successful transition.
In summary, we have described a clear example of a homeotic mutant that confers a major morphological shift (petal-to-sepal transition) and that has been under positive selection largely due to non-pollinator selection (herbivory avoidance and resource allocation). In addition, the radical change in floral morphology alters pollination dynamics causing a degree of assortative mating that could ultimately promote speciation. As such, these findings fit well into the concept of a ''hopeful monster'' coined by Goldschmidt where major morphological shifts that delineate taxa arise in a single step. 69 Regardless, our results do show that homeotic changes in morphology due to single mutations can be selected for in nature.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
The authors declare no competing interests.

INCLUSION AND DIVERSITY
One or more of the authors of this paper self-identifies as an underrepresented ethnic minority in science. One or more of the authors of this paper self-identifies as a member of the LGBTQ+ community.

RESOURCE AVAILABILITY
Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Zachary Cabin (zcabin@ucsb.edu)

Materials availability
This study did not generate new, unique reagents.

Data and code availability
Amplicon sequencing data is deposited on the SRA under Bioproject: PRJNA660745. CO1 sequences are available in the GenBank Database: MW302492-MW032499. All statistics were calculated in R v4.0.3. All code is available upon request. Any additional information required to reanalyze the data reported in this study is available from the lead contact upon request. . Fruit production was assessed in early August and 2-3 seeds were collected from fruits. During this time, evidence for herbivore damage was again noted. Aphid evidence was simply presence/absence. Caterpillars were either identified on the flower, or through evidence such as partially eaten spurs, holes in floral tissue, or partial removal of anthers, staminodia, or carpels ( Figure S2E). During seed collection, caterpillar evidence could be identified through partially eaten carpels or, when no carpels remained, the floral receptacle was intact atop the inflorescence. Mule deer were observed to eat entire flowers leaving characteristic stumps on the flower pedicles or main inflorescence ( Figure S2F; Video S2). Intermediate morphology flowers (''short spur'', ''SS'', Figures S2A and S2B) were seen occasionally throughout the population (<1% frequency). Generally, flowers of the SS morphology had typical outer-whorl sepals, but 2nd whorl organs tended to be entirely blue (similar to d morph) but had some presence of nectar spurs. These spurs were inconsistent in length and were not uniformly present (not every 2nd whorl organ had a spur) both between and within plants. Nectary size, which is highly correlated with nectar volume, 78 was highly reduced. Since a high majority of these plants were genotyped as d/d (homozygous and heterozygous, see Table S5), these plants were included as mutant morphology in all analyses. In 2014, a video camera (GoPro) was set to simultaneously film both WT and d flowers in an attempt to capture pollinator visitation (Video S1; Table S2). Patches of flowers were chosen that included both morphs in a similar developmental phase (e.g. male (some but not all anthers dehisced)).

Population mapping
To create the spatial map, two reflective targets (T1, T2) were set a fixed distance apart (D T ). Then, using a laser distance measurer (Leica Disto D2 Laser Distance Measurer), the distance from each plant in a ''target group'' to T1 (D 1 ) and T2 (D 2 ) was measured.
It was impossible to measure the distance for all plants to a single pair of targets due to obstructions such as trees and thus sub-maps of smaller portions of the population were created and then linked together to create a map of the entire population. At least three plants were included in both adjacent sub-maps in order to link them together. For example, there were 16 sub-maps used to create the map of the entire population in 2016 ( Figure 1C). Simple algebra (Equation 1) was used to convert measured distances (D 1 , D 2 , D T ) into unique x-y coordinates for each sub-map (Data S2). In order to link these sub-maps, two neighboring sub-maps (e.g. A, B) would be plotted at the same scale. One sub-map (B) would then be rotated q degrees until the plants that spanned both sub-maps overlapped. Equation 2 was then used to transform the x-y coordinates of sub-map B to match the coordinates of sub-map A, creating a unified map. This process was repeated until all sub-maps were unified. This was done for 2014, 2015 and 2016, but only the map from 2016 is presented here. Coordinates for all years are available in Data S2. More detailed methods and equations available upon request. Equation 1: converting distances to x-y coordinates : rotating x-y coordinates x 0 = ½x Ã cosðqÞ À ½y Ã sinðqÞ Caterpillar identification Caterpillars ( Figure S2E) were collected in 100% EtOH. DNA was extracted following Gloor and Engles. 79 The extraction buffer consisted of 10 mM TRIS (pH 8.2), 1 mM EDTA (pH 8.0), 25 mM NaCl, and 200 mg/mL Proteinase K. A small portion of each caterpillar (5 mm) was submerged in 50 mL of the extraction buffer in 0.5 mL tubes and smashed with a pipette tip for 10 seconds. Each 0.5 mL tube would then incubate at 37 C for 25 minutes followed by a 2 minutes at 95 C to deactivate Proteinase K. DNA concentration ranged from 3-40 ng/mL Prior to amplification, DNA was diluted 1:100. Approximately 650 bp of the mitochondrial gene cytochrome oxidase 1 was PCR amplified using the LepF and LepR primers from Hajibabaei 80 and an initial denaturation at 95 C for 2 min followed by 35 cycles at 95 C for 30 sec, then 57.6 C for 30 sec, then up to 72 C for 1 min. Sanger sequencing was conducted at the UC Berkeley Sequencing Facilities. Sequences were trimmed using Geneious version 9.1.6 and consensus sequences for each sample were generated (http://www.geneious.com) 77 and then blasted against the NCBI database ( Figure S3). Sequences are available on the GenBank Database: MW302492-MW032499.

DNA extraction
Field collected leaf tissue was collected every year (2014-2016) in Reynolds Park for every flowering plant. Leaf tissue was stored in silica gel desiccant until processing at UCSB. Samples grown in the greenhouse were snap-frozen in liquid nitrogen prior to DNA extraction. Approximately 20 mg of dried or frozen tissue was placed in a 2 mL tube with 3-4 ball bearings. Samples were placed in a bead-beater to grind leaf tissue into a fine powder. Subsequently, all DNA was extracted using the MagAttract 96 DNA Plant Core Kit (Qiagen) on the BioSprint 96 (Qiagen) with accompanying reagents and consumables. Samples were quantified using a Qubit 2.0 Fluorometer (Thermo Fisher Scientific).

PCR amplification and library preparation
Primers were designed (F: 5 0 GAGAGACCTTGGTGGGGAGA 3 0 , R: 5 0 AGCCAGCTTTACCGTACACC 3 0 ) around AqAP3-3 (locus identifier: Aqcoe5G180800.1) to amplify the entire gene and $1 kb up and downstream (total amplicon size = $3.6 kb) with Qiagen Multiplex PCR Kit (Qiagen). The PCR amplification protocol consisted of initial denaturation at 95 C for 15 min followed by 36 cycles at 94 C for 30 sec, then 61.9 C for 1:30, then 72 C for 4 min. PCR products were quantified using a Qubit 2.0 fluorometer (ThermoFisher Scientific). Dual-indexed libraries were prepared using ½ reactions of NEBNext Ultra II (New England Biolabs) for Illumina and 10ng of input DNA. Libraries were sequenced on the Illumina NextSeq 500 (UCSB CNSI BNL).

SNP genotyping
In order to get accurate allele frequencies for d alleles, real-time PCR allele specific SNP Genotyping was used. 39

Spatial analyses
Spatial distribution of all plants in 2016 was modeled in R 70 following the methods of Kellenberger. 81 Initially, the map was subset into 16 windows of the same size (17m x 30m, Figure S1A). A Kolmogorov-Smirnov test was used to determine overall randomness of the spatial distribution of plants in the population. Differences in spatial distribution patterns between D and WT morphs were computed with a Studentised Permutation Test with 999 random permutations ( Figure S1B). All spatial analyses were implemented in R with the package spatstat v.1.64-1. 82

Floral measurements
A Wilcoxon rank sum test was used to assess whether either morph produced more flowers per plant for all three flowering seasons ( Figure 2B). In 2016 entire, undehisced fruits lacking evidence of herbivores but having a range of sizes within 1 km of the study population were collected and the height (h) and width in two directions (w 1 , w 2 ) of each fruit was measured (Data S1) and the total number of developed seeds were counted. The volume of the fruit was modeled as an ellipsoid with the Equation 3: Equation 3: calculating fruit volume Fruit volume and seed number was correlated using a simple linear regression ( Figure 2C). In 2016 fruit size was measured for all fruits in the study population (Data S1). The fruit volume had to be square root transformed to approximate normality before it could be tested whether floral morphs differed in fruit volume ( Figure 2D).

Variant identification
Sequences were aligned to the A. coerulea 'Goldsmith' v3.1 reference genome 83 (https://phytozome.jgi.doe.gov) using the Burrows-Wheeler aligner 73 and then sorted and indexed with samtools 01. 1.19. 72 Sequence data is available in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/sra) under BioProject: PRJNA660745. Sequences were initially scanned for large structural mutations using the IGV browser and subsequently excluded those regions from the analysis. Small indels and SNPs were identified using GATK (HaplotypeCaller (-dontusesoftclippedbases), CombineGVCFs, GenotypeGVCFs). 71 For variants in the coding region, it was determined if they altered amino acid positions and if any appeared only in d plants. See Data S4 for a full list of samples, variable positions, and excluded regions.
Sample and variant filtration for haplotype phasing Average sequence coverage was 1278x (GATK DepthOfCoverage). Many samples showed clear amplification bias for one allele ( Figure S4). Thus, samples with %75x (n = 22) were excluded from the haplotype phasing as their genotype calls were often incorrect. Variant sites that fell within large structural mutations were also excluded (Data S4) and the remaining variants were filtered to include only biallelic sites with a MAF R 5% and high mapping quality (MBQ > 30). This MAF cutoff resulted in the two rarest d mutations to be excluded (d 1b and d 682 ) though samples carrying these alleles were flagged for downstream analyses. One variable position was added to represent a 280bp deletion (Data S4).

Haplotype phasing and analysis
Filtered VCFs were formatted with GATK and PLINK 74 and haplotypes were constructed using SHAPEIT2 (v2.r904, default settings, UCSB KNOT). 75 Since some greenhouse plants were offspring of sequenced field-collected plants, for parent-offspring pairs that were both sequenced and phased haplotypes were removed that were identical by descent. To be conservative, if both the parent and offspring were heterozygous for the same haplotypes, both of the offspring haplotypes were removed from the analysis. Lastly, 11 haplotypes were removed as they were the clear result of allelic drop out.
Allelic drop out occurs when one chromosome (allele) fails to amplify during the initial PCR reaction and therefore does not get sequenced. Allelic drop out could be identified by finding individuals that, according to the sequence data, were homozygous for a single d allele and at every position across the entire amplicon, but WT morph. To be conservative, all of these individuals were kept in the analysis but only contributed a single haplotype to the EHH calculations. The authors note that allelic drop out could only be identified in d individuals, raising the possibility that it occurred in WT plants and went unnoticed, possibly inflating homozygosity measures for wt haplotypes. This implies that the estimation of WT haplotype decay might be rather conservative. If a sample suffered from allelic drop out and the remaining haplotype was identical-by-descent, then that haplotype was removed as well (effectively removing that individual from the analysis). Lastly, a minimum spanning haplotype network was created (http://popart.otago. ac.nz) 40 and then specific positions that separated d haplotypes and WT haplotypes were identified ( Figure 4B). Each of these, with the exception of one position, also separated functional WT haplotypes. The position in question (10,481,304, marked as the ''triple arrowhead'' in Figure 4B) only appears on haplotypes 2, 3, and 21 (all d 1 haplotypes). This mutation is a C/T transition and occurs upstream of the gene. While it only occurs on those three haplotypes, it is not present on haplotype 12, the other d 1 haplotype suggesting it does not have a substantial effect on function.

Genotypic linkage disequilibrium
To further confirm the apparent complete linkage disequilibrium between d alleles, the SNP genotyping data was to calculate pairwise genotypic linkage disequilibrium between the three most common d alleles (d 1 , d 9 , d 1b ). The -r2 tool in PLINK (v1.07) calculates r 2 between multiple sites with unphased genotypic data. 74 ll OPEN ACCESS e5 Current Biology 32, 1332-1341.e1-e5, March 28, 2022