# Fine-scale genetic differentiation in the bee-specialized Antirrhinum charidemi covaries more strongly with spatial isolation than with corolla colour

Myriam Heuertz, Marcial Escudero, José María Gómez, Pablo Vargas

► To cite this version:

Myriam Heuertz, Marcial Escudero, José María Gómez, Pablo Vargas. Fine-scale genetic differentiation in the bee-specialized Antirrhinum charidemi covaries more strongly with spatial isolation than with corolla colour. AoB Plants, 2025, 17 (4), plaf017. 10.1093/aobpla/plaf017. hal-05504967

HAL Id: hal-05504967

https://hal.inrae.fr/hal-05504967v1

Submitted on 11 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Studies

# Fine-scale genetic differentiation in the bee-specialized
*Antirrhinum charidemi* covaries more strongly with spatial
isolation than with corolla colour

Myriam Heuertz<fnref n="1" />,‡, ID, Marcial Escudero<fnref n="2" />,†, ID, José María Gómez<fnref n="3" />,<fnref n="4" />, ID, Pablo Vargas<fnref n="5" />,*, ID

**Abstract** The snapdragon *Antirrhinum majus* has been a model species for genetics, plant development, and evolution since the 19th century. Recent studies have expanded the focus to the entire *Antirrhinum* genus as a model system for rapid evolution (26 species in < 5 million years). However, in-depth studies to reconstruct microevolution in additional snapdragon species are lacking. This study aimed to explore to what extent potential pollinators, flower colour morphs, spatial and environmental factors contribute to differentiation in a small population of the Mediterranean *A. charidemi* (south-eastern Spain). We studied a population of approximately 200 *A. charidemi* individuals with either pink or white corollas, characterized by strong topographic heterogeneity (horizontal extension of 120 x 80 meters; 40-meter altitude difference) and diversity in environmental factors (substrate, vegetation). The study analysed pollinator preference for either white or pink corollas, genetic diversity using 13 nuclear SSR loci and three plastid haplotypes, and the spatial population structure. Flower visitors displayed some indication of preference for pink corollas (five of ten bee species) and flower colour morphs were genetically differentiated. However, the strongest pattern of genetic differentiation was associated with a fine-scale spatio-topographic isolation in the population, with five topo-genetic subpopulations and a pollen-to-seed dispersal distance ratio of 4.32. Our results agree with similar patterns of strong spatial genetic isolation found in *A. charidemi* at larger scales: phylogeographic differentiation of populations and phylogenetic relationships within a south-eastern Iberian *Antirrhinum* clade. Despite the extreme corolla specialization for bee pollination, spatial isolation appears to be the predominant factor driving short- and long-term differentiation in *A. charidemi*. We argue that a comprehensive understanding of early stages of rapid evolution requires detailed investigation of fine-scale evolutionary drivers, including both spatial isolation (topography) and ecological factors (e.g. pollination fauna).

**Keywords:** *Antirrhinum charidemi*, bee specialization; flower colours; genetic structure; nuclear SSRs; pink morph; plastid haplotypes; spatial isolation; white morph; topography.

## Introduction

In most angiosperms, flower colouration has historically been attributed to a decisive role in driving speciation through pollinator-mediated reproductive isolation (Narbona et al. 2021). Indeed, we can find a wide range of colour patterns among flowering plants. They range from large clades of species with a single colour morph (e.g. white, purple, or yellow in *Cistus*, Guzmán and Vargas 2005) to two or more colour morphs within the same population of a single species (Streisfeld et al. 2013; Del Valle et al. 2019; Jiménez-López et al. 2023). However, it has been argued that pollinators are not always responsible for evolutionary transitions in flower colour variation (Rausher 2008). A recent review showed that there is little evidence for significant pollinator-mediated selection among existing studies (Trunschke et al. 2021).

Colour transitions unrelated to pollinator-mediated selection can result from climatic variation, antagonist-mediated selection, genetic drift, pleiotropic effects, or genes influencing anthocyanin pathways in vegetative rather than reproductive traits (Rausher 2008; Koski and Ashman 2016; Smith 2016; Del Valle et al. 2019).

Anthocyanins play a crucial role in plant physiology, providing antioxidant functions and protecting plants from various biotic and abiotic stresses, but are also key pigments responsible for the diverse hues of red, purple, and blue in plant organs like fruits and flowers (Torres et al. 2002; Landi et al. 2015; Mannino et al. 2021). The occurrence of anthocyanins associated with purple flowers is commonplace in many floras. For instance, the flora of the Mediterranean floristic region (around 25,000 species) has a high percentage

<footnote n="1">¹Univ. Bordeaux, INRAE, BIOGECO, 69 route d'Arcachon, F-33610 Cestas, France</footnote>

<footnote n="2">²Department of Plant Biology and Ecology, University of Seville, Reina Mercedes 6, E-41012 Seville, Spain</footnote>

<footnote n="3">³Estación Experimental de Zonas Aridas (EEZA-CSIC), Ctra Sacramento s/n, La Cañada de San Urbano, E-04120 Almería, Spain</footnote>

<footnote n="4">⁴Research Unit Modeling Nature (MNat), Universidad de Granada, E-18071 Granada, Spain</footnote>

<footnote n="5">⁵Real Jardín Botánico, RJB-CSIC, Plaza de Murillo 2, E-28014 Madrid, Spain</footnote>

<footnote n="1">¹Equal contributions.</footnote>

<footnote n="2">²*Corresponding author. Pablo Vargas, Real Jardín Botánico, RJB-CSIC, Plaza de Murillo 2, E-28014 Madrid, Spain. E-mail: vargas@rjb.csic.es*</footnote>

<footnote>Associate Editor: Antonio Castilla</footnote>
(87.7%) of species with flower colour polymorphisms (Kantsa et al. 2017; Narbona et al. 2018). In the Mediterranean basin, snapdragons (species of *Antirrhinum*, Plantaginaceae) form a plant group that exhibits one of the highest levels of colour polymorphism of flowers. *Antirrhinum majus* L. has been historically used as a model plant in genetics and flower development (Schwarz-Sommer et al. 2003), including studies of flower colour (purple, yellow, white) transitions at the populational level (Whibley et al. 2006; Bradley et al. 2017; Tavares et al. 2018). Despite the well-studied role of colouration in *Antirrhinum majus*, the evolutionary significance of flower colour variation in other *Antirrhinum* species remains unclear.

In addition to colouration, flower shape has been proposed as a strong physical barrier to pollinators of all the species of *Antirrhinum* (Sutton 1988; Magallón and Vargas 2014). Indeed, the constant flower phenotype of snapdragons displays one of the most occluded corollas among angiosperms (Fig. 1), which is only accessible by bees (Torres et al. 2002; Vargas et al. 2010, 2017; Magallón and Vargas 2014). As a result, a tight mutualism between snapdragons and bumblebees has been historically proposed (Macior 1967; Odell et al. 1999; Dyer et al. 2007; Whitney et al. 2013). Once bee pollination was confirmed, researchers investigated in detail the relationship among occlusion, corolla colour variation, and pollination by bumblebees in *Antirrhinum majus* variants in order to elucidate the genetic pathways underlying the evolution of these floral traits (Whibley et al. 2006; Shang et al. 2011; Bradley et al. 2017; Tavares et al. 2018; Qiao et al.

2023). Numerous feeding experiments using *A. majus* have demonstrated that: (a) *Bombus* species discriminate flower colouration (Comba et al. 2000); (b) *Bombus terrestris* associates warmth with corolla floral colour (Dyer et al. 2007); (c) single mutations affecting flower signal (colour, reflectance) can have profound effects on pollinator behaviour (Dyer et al. 2007; Moss et al. 2024); (d) a corolla with conical epidermal cells can provide fitness benefits (Glover and Martin 1998; Whitney et al. 2011); e) the corolla venation patterning might have evolved as a response to bees learning of their use as nectar guides (Whitney et al. 2013; Tavares et al. 2018). To link the evolutionary patterns between *A. majus* and relatives sharing the distinctive occluded flower shape, numerous studies have contributed to a broader understanding of this unique floral morphology, which appears to have emerged between 30 and 45 million years ago (Vargas et al. 2009, 2014; Gorospe et al. 2020).

Taxonomic treatments of *Antirrhinum* considered morphological characters including flower colours to group the approximately 26 currently accepted species into three subsections: Antirrhinum (11 species with purple, magenta, and yellow corollas); Streptosepalum (two species with yellow corollas); and Kickxiella (13 species displaying mainly whitish flowers) (Rothmaler 1956; Otero et al. 2021). Phylogenetic results do not support these groupings and a more dynamic transition of corolla colours appears to have occurred in the course of evolution of *Antirrhinum* (Otero et al. 2021; Durán-Castillo et al. 2022). At a macroevolutionary scale, the white corolla is the predominant phenotype in the majority

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Studied population of *Antirrhinum charidemi* in the Barranco del Dragoncillo Blanco (Almería, Spain). View of the ravine where the c. 200 individuals occurred (2010) (A); one of the 10 individuals with white corollas (B); close-up of the two corolla types found on each individual showing the morphs pink (C) or white (D). Photo credit: Myriam Heuertz (A, B) and Pablo Vargas (C, D).</figcaption></figure>
of *Antirrhinum* species, which are intermingled in the three main clades (Otero *et al.* 2021; Durán-Castillo *et al.* 2022). In contrast, a predominant pattern of geographic speciation in *Antirrhinum* was inferred using both genetic (Vargas *et al.* 2009) and genomic (Otero *et al.* 2021) data. The question remains as to what extent white corollas have been favoured at early stages of population differentiation (microevolutionary scale).

Based on more than 30 studies on *Antirrhinum*, two non-exclusive hypotheses are proposed to account for the rapid differentiation (26 species) in the last five million years (Vargas *et al.* 2009). On the one hand, a tight snapdragon-bumblebee dependence in self-incompatible species (Vargas *et al.* 2017) indicates that bees may have triggered shifts in flower morphology and corolla colouration in the course of evolution. On the other hand, the narrow distribution of a considerable number of snapdragons (9 species) across the Iberian Peninsula (Güemes 2009; Vargas *et al.* 2009; Vargas 2024), together with a biogeographic pattern for the whole genus (Otero *et al.* 2021), suggests that geography has been key for speciation. The two hypotheses were addressed in an analysis of pollinators and geographic distribution of the species, primarily 11 bumblebee species across 18 species of *Antirrhinum* across the Iberian Peninsula (Vargas *et al.* 2017), suggesting that both hypotheses were likely at play. The extent to which early stages of evolutionary divergence (i.e. population differentiation) are primarily driven by pollinator preference for flower colouration or by spatial isolation is, however, unclear. In other words, is isolation by ecology (pollinator preference) or isolation by distance (spatial differentiation) predominantly responsible for triggering early genetic differentiation within *Antirrhinum*?

To address the question of the drivers of incipient evolutionary differentiation, we examined flower colouration, pollinator behaviour, and population genetic structure in a topographically complex population formed by plants with either pink or white corollas of a narrow-endemic snapdragon species, *Antirrhinum charidemi* Lange in southeastern Spain. In particular, we addressed the following hypotheses:

* Hypothesis 1: Under a scenario of pollinator preference for flower colour, pollinators are expected to discriminate by corolla colouration and will thus preferentially visit and pollinate plants of their favoured colour morph. Over time, assortative mating due to pollinator preference is expected to lead to genetic divergence patterns associated with individuals with a particular flower colour (Rymer *et al.* 2010).

* Hypothesis 2: Alternatively, under a spatial isolation scenario, early genetic differentiation is associated with distance and topography rather than flower colouration.

* Hypothesis 3: Since genetic structure is not only shaped by patterns of pollen dispersal (mediated by pollinators) but also by spatially limited seed dispersal (barochory, wind), complementary information on the relative roles of both dispersal vectors can be gained by examining genetic structure and gene dispersal estimates from maternally and biparentally inherited markers (plastid sequences vs. SSRs).

In addition, some other ecological characteristics were addressed to test the contribution of local conditions such as altitude, substrate, and co-occurring vegetation.

# Materials and methods

## Study species and focal population

*Antirrhinum charidemi* Lange is a perennial shrub endemic to an 8 × 4 km region in the volcanic mountains of Cabo de Gata (Almería, south-eastern Spain), where populations are formed by few to hundreds of individuals (Forrest *et al.* 2017). The flower colour is generally pink, with weak variation in colour intensity. Some flowers are observed all year round in certain individuals, although all individuals are flowering between March and June (Cueto *et al.* 2008; Vargas *et al.* 2010; M. Heuertz, personal observation). We focussed our study on the Barranco del Dragoncillo Blanco (BDB, 36°46'55'' N, 2°07'54'' W) population of ca. 200 individuals distributed in a ravine (barranco in Spanish), on rocky outcrops and on scree slopes bordering the ravine (Figs. 1A and 2). This is the population where the highest genetic diversity was found in *A. charidemi* (Forrest *et al.* 2017). We visited the area regularly between 2009 and 2010. The flowering peak in 2009 was reached in May, including 10 individuals displaying white corollas (Fig. 1B) and 172 individuals displaying pink corollas (Fig. 1C).

## Flower colour variation

The five petals of any *Antirrhinum* species display a personate flower, i.e. the corolla has a tube, an adaxial or upper lip (two fused upper petals), and a lower or abaxial lip (three fused lower petals) that typically develops a basal convexity (palate) strongly occluding the mouth of the tube (see Vargas 2024 at www.rjb.csic.es/snapdragons) (Fig. 1C and D). All the flowers of *A. charidemi* display the same occluded corolla shape and no size differences between pink and white corollas were observed (typically corolla tubes are about 19 mm long and 6 mm wide; see Vargas *et al.* 2010). Both colour morphs have yellow spots on the tube sides and the lower lip where it touches the upper lip (corolla access) (Fig. 1C and D). Every flower lasts 5-7 days (P. Vargas personal observation) and produces one of the highest amounts of sugar (4.5 mg per flower in 24 h) among the flowers of eastern Andalusia (Herrera 1988). Reflection of UV light was measured for the flowers of two individuals of each colour morph, using an USB4000 miniature fiber optic spectrometer with a USB-DT deuterium tungsten halogen source (Ocean Optics, Dunedin, Florida, USA).

## Preference pattern of pollinators for flower colour

To test whether pollinators exhibited any preference pattern for flower colour (Hypothesis 1), we conducted a pollinator observation experiment. During 2007, we established three subplots, in each of which we selected one *A. charidemi* individual with white and one of similar size with pink corollas growing side-by-side. During the main flowering period (March to June) we examined for 15-minute periods (trials hereafter) all insects visiting the flowers of the experimental plants, counting only those that successfully opened the corolla and contacted the reproductive organs. Although no direct tests were conducted to estimate the success of pollen deposition in visited flowers, we considered that flower visitors (typically bees) are likely actual pollinators of *A. charidemi*. This consideration was based on (i) low flower visit ratio by insects other than bees in *Antirrhinum* (Vargas *et al.* 2017), (ii) high efficiency of bees in specialized flowers (Fenster *et al.* 2004), and (iii) high levels of fruit and seed sets
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Spatial distribution of adult plants of *Antirrhinum charidemi* in the BDB (*Barranco del Dragoncillo Blanco*) population in 2009-2010. A. Satellite image from Google Earth: white-flowered individuals are represented by white symbols and pink-flowered individuals by pink symbols. ©2017 Google Data SIO, NOAA, U.S. Navy NGA, GEBCO, Image Landsat/Copernicus, Imagery date 6/22/2016, lat. 36.78°, lon -2.14°, elev. 316m, eye alt. 403 m. B. Distribution of plastid DNA haplotypes (H1, H2, H3) represented on a terrain model in QGis. Note that the direction of North is not the same in A and B.</figcaption></figure>

when bees touched flower reproductive parts with their large
body sizes and suitable structures (hairy scutum) (Vargas et
al. 2010). Visiting insects were identified using local faunas
and taxonomic treatments (Ornosa and Ortiz-Sánchez 2004;

Michener 2007; Vargas *et al.* 2010). To control for the dif-
ference in flower availability, the number of flowers on the
experimental plants was counted during each trial. We per-
formed all trials on sunny days, with a temperature of over
20°C and without wind. We also visited populations of the species when weather conditions were different (rain, wind, low temperature) and hardly observed insects visiting the plants (see Vargas *et al.* 2017). In total, we conducted 140 trials, with a total of 2,100 minutes (35 h) of observation.

The preference of each insect visiting the flowers of the ex-
perimental plants was estimated using the Jacobs' D index
(Jacobs 1974; Lechowicz 1982). This index compares the
relative availability of each flower colour and their rela-
tive visits by each pollinator, taking a value of zero under
random feeding and deviating symmetrically from zero be-
tween 1 and -1 for preferred and avoided items, respectively.
Jacob's D index was calculated using the function *ivlev* in the
R package *selectapref* (Richardson 2020). The statistical sig-
nificance of the preferences was obtained by testing for devi-
ation from random visitation rate using a test for goodness
of fit (Lechowicz 1982). Under random visitation, the null
hypothesis is that the proportion of insects visiting a given
flower colour equals the proportion of the flower colour in
the environment (Lechowicz 1982). In addition, to analyse if
flowers with different colour attract different potential pol-
linators, we tested whether the composition of the visiting in-
sect assemblage varied between plants with different flower
colours by means of a Permutational Analysis of Variance
(Anderson 2001). This analysis was performed using the
function *adonis2* in the R package *vegan* on a Bray-Curtis
dissimilarity matrix (Oksanen *et al.* 2001).

Sampling of plant material and ecological data
collection

Leaf samples were collected from 182 A. charidemi indi-
viduals at BDB between 2009 and 2010 and samples were
silica-dried. The sampling process involved meticulously
scouting the studied area where plants were located within
ca. 120 × 80 m, with 40 m difference in altitude, and its bor-
dering area. For each individual, we recorded geographic
position and altitude using a Topcon geospatial station (map-
ping ca. 80% of the individuals with a precision of 2 mm for
distances up to 500m) or a handheld GPS device (Garmin
Oregon, used in subsequent field visits with multiple position
measurements per individual to reach a GPS precision of ca.
2m). We also recorded plant height and diameter, as well as
vegetation at immediate contact with the target plant possibly
preventing herbivory and providing a favourable microenvir-
onment for germination. DNA was extracted from ca. 15 mg
of dried leaf tissue using the Invisorb DNA Plant HTS 96 Kit
(Invitek, Berlin, Germany).

Nuclear microsatellite analyses

We amplified 14 microsatellite (SSR) loci in two multiplexes using the QIAGEN multiplex PCR kit (Qiagen) following the manufacturer's protocol, as described in Forrest *et al.* (2017). Forward primers were fluorescently labelled and fragments were resolved on an ABI3730 DNA analyzer (Applied Biosystems). Allele sizes were determined using Genemapper 3.7 (Applied Biosystems) by comparison with the GeneScan™ -500 LIZ size standard.

Plastid DNA analyses

Two plastid DNA (pDNA) spacers (*trn*D-T and *trn*S-fM) were chosen based on *Forrest et al.* (2017), PCR amplified and sequenced in all individuals using the services of

Macrogen (dna.macrogen.com). Sequences were assembled and edited with CodonCode Aligner ver. 3.7.1 (CodonCode Corporation, MA).

Genetic data analyses

**Genetic diversity** To assess SSR marker quality we checked for null alleles using Micro-checker version 2.2.3 (van Oosterhout *et al.* 2004) and assessed linkage disequilibria between loci using Genepop version 4.7.5 (Rousset 2007), correcting for multiple testing with a sequential Bonferroni correction. We estimated genetic (pDNA and SSRs) diversity parameters for the 182 individuals of the BDB population using the software SPAGeDi ver. 1.5a (Hardy and Vekemans 2002), as well as for groupings according to (i) flower colour or (ii) geo-genetic subpopulations, defined as geographically delimited (spatio-topographic) subpopulations in which individuals generally shared genetic cluster membership (see STRUCTURE analysis described below). For pDNA, we assessed allelic (haplotype) richness using rarefaction, AR, and gene diversity, $H_e$, corrected for sample size. For SSRs, we estimated AR, $H_e$, and the individual inbreeding coefficient $F_{IS}$ and tested for deviation from Hardy-Weinberg genotypic proportions using permutation tests (10 000 permutations).

**Determinants of within-population genetic structure** To assess the effect of ecological predictors on genetic differentiation within BDB, we used a matrix-based modelling approach that allowed us to evaluate the effect of flower colour (through pollinator preference, Hypothesis 1), the effect of spatial isolation (Hypothesis 2), and effects of additional predictors, alone or in combination. Genetic similarity between pairs of individuals was assessed by means of the kinship coefficient $F_{ij}$. (Loiselle *et al.* 1995) between individuals i and j using SPAGeDi. The nuclear kinship matrix was set as the response variable (the plastid matrix was not used because of too low variation). We ran eleven models using five predictor variables and combinations thereof using the Bayesian generalized linear mixed modelling approach implemented in the *MCMCglmm* package (Hadfield 2010) in R: (i) NULL, a null model without any predictor, (ii) GEO, with a predictor matrix representing spatial distance between individuals, computed from easting and northing (latitude and longitude), (iii) ALT, with a distance matrix based on altitude only, (iv) SUB, with a matrix based on sharing topography (ravine vs. scree slopes/outcrops, pairs of individuals on the same substrate are coded as ‘1’; pairs on distinct substrates as ‘0’), (v) COL, with a matrix based on sharing of corolla colour, (vi) VEG, with a matrix predictor based on sharing vegetation at individual plant contact, i.e. conditions potentially facilitating seedling establishment in harsh environmental conditions. To construct this matrix, we assessed whether individuals were in immediate contact with other plants; 16% missing data were imputed using the *mice* package (van Buuren and Groothuis-Oudshoorn 2011) in R. We also applied the following models with combinations of these matrices. (vii). GEOLAlt, 8. GEO))^Alt))^SUB, 9. COL^VEG, 10. GEO))^Alt)^COL, 11. GEO))^COL matrices. MCMCglmm was initiated with standard priors and run with a burn-in of 5000 followed by 20 000 iterations with a thinning interval of 750. Model comparison was conducted by means of comparing the Deviation Information Criterion, DIC.
**Characterization of within-population genetic structure**
To test for spatial autocorrelation of genetic relatedness in each plastid and nuclear data, the matrices of inter-individual kinship coefficients (see above) were regressed on the logarithm of geographic distance to yield regression slopes, $b_p$ for plastid and $b_n$ for nuclear DNA. These slopes were compared to their expected distributions obtained from 10 000 permutations of rows and columns of the spatial distance matrix. To quantify the strength of SGS, we computed *Sp* statistics ($Sp_p$ for plastid and $Sp_n$ for nuclear DNA) from the regression slopes as $-b/(1 - F_1)$, where $F_1$ is the average kinship coefficient for 'neighbouring' individuals, i.e. for all pairs of individuals in the first distance class (Vekemans and Hardy 2004). Seven distance classes were defined for computation of $F_1$ and graphical representation, with maximal inter-individual distances of 3, 6, 12, 24, 48, 96, and 192 m.

To additionally explore non-linear SGS patterns, we used a spatial principal component analysis (sPCA) in the *adegenet* package (Jombart 2008) in R. Spatial genetic autocorrelation (e.g. patches of related individuals, allele frequency gradients) was tested using a G-test, whereas inter-individual genetic repulsion was tested with an L-test. The strength of SGS was estimated as the eigenvalue of the first sPCA axis, eig.sPCA.

We further used genetic clustering methods; specifically, the STRUCTURE (Pritchard et al. 2000) and GENELAND (Guillot et al. 2005) programmes, respectively ignoring or including the spatial position of individuals as prior information. STRUCTURE ver. 2.3.2.1 (Pritchard et al. 2000) was run using the admixture model with correlated allele frequencies and $K = 1$ to $K = 12$ genetic clusters. For each $K$, ten runs were performed with 100,000 iterations discarded as burnin, followed by 500,000 iterations. Analyses were parallelized on the Genotoul computer cluster, INRAE Toulouse, France, using the StrAuto software (Chhatre and Emerson 2017). The $K$ that best explained the data was inferred using maximization of lnP(data) and the deltaK criterion in StructureHarvester (Earl 2012). STRUCTURE results were combined across runs using Clumpak (Kopelman et al. 2015) and the ancestry proportions of each individual in each gene pool were mapped as pie charts and using spatial interpolation (Inverse Distance Weighting) with QGIS Desktop (QGIS 2017). In GENELAND (Guillot et al. 2005), five MCMC chains were run using a spatial model with correlated allele frequencies, simulating $K$ populations in Hardy-Weinberg equilibrium and unlinked loci (allowing $K$ from 1 to 20; Guillot et al. 2005; Guillot 2008). We used 100,000 MCMC iterations with a thinning of 100 and setting the burn-in to 20,000 iterations.

To comparatively address the support of our data to Hypotheses 1 and 2, we computed genetic differentiation as $F_{ST}$ between flower colour morphs, among STRUCTURE gene pools assigning individuals to gene pools based on $q > 0.5$ gene pool membership, and among geo-genetic subpopulations to which individuals were assigned based on their spatiotopographic position and the results from STRUCTURE.

**Inference of gene flow patterns**
To address Hypothesis 3 on the intra-population gene flow and its components due to seed and pollen dispersal, we first used coalescent simulations and Bayesian inference in the program MIGRATE (Beerli 2008), and second, the *Sp* estimates from SPAGeDi. For the MIGRATE analysis, we considered the geo-genetic subpopulations previously defined. We assumed a constant

mutation rate across loci for the mutation-scaled effective population size $\Theta$ and ran migration matrix models with a variable $\Theta$ and, average (model 'mean_m'), symmetric (model 'symmetric_m') or freely variable (model 'free') migration rates between all combinations of subpopulations. We also implemented models with a fix population parameter $\Theta$ and an average (model 'fixtheta_mean') or symmetric (model 'fixtheta_sym') migration rates. We used uniform-wide priors and ran two replicates of one long Markov chain with 10 000 steps discarded as burnin and 5000 steps retained for estimates, with a sampling increment of 100 steps. The heating scheme was static with four chains of distinct temperatures (1, 1.5, 3, 1,000,000), swapping at each step. Model choice was conducted using Bayes Factors based on the Bezier approximated score for Log-Probability of the data given the model. Estimates of migration rates were plotted on the population map using the Flow mapping plugin (Gulluoglu 2016) in QGIS Desktop.

We further estimated historical dispersal distances from *Sp* under the assumption of drift-dispersal equilibrium in the population. In hermaphroditic plants, the *Sp* statistics are expected to approach $1/(4\pi D_e \sigma_g^2)$ for nuclear (biparentally inherited) and $1/(2\pi D_e \sigma_s^2)$ for plastid (maternally inherited) markers at equilibrium, where $\sigma_g^2$ is half the mean squared dispersal distance of genes, $\sigma_s^2$ is half the mean squared dispersal distance of seeds, and $D_e$ is the effective population density (Vekemans and Hardy 2004). Thus, *Sp* computed separately from nuclear and plastid markers can provide information on the relative distances of gene vs. seed dispersal, and by inference, of the pollen to seed dispersal ratio. We used $D_e = 0.01$ individual per m², corresponding to half the census population density, thus assuming a relatively balanced contribution of individuals to reproduction. This value was justified by the long flowering season and observation of flowers in all adults. We used an iterative procedure for the estimation of $\sigma_g$ and $\sigma_s$ from the kinship-In distanceratio (Heuertz et al. 2003). The relative contribution of seed dispersal to overall gene dispersal is then given by the $\sigma_g/\sigma_s$ ratio and applying $\sigma_g^2 = \sigma_s^2 + 1/2 \sigma_p^2$ (Crawford 1984), we obtained the pollen-to-seed dispersal distance ratio as $\sigma_g/\sigma_s$.

## Results

### Corolla features and pollinator behaviour

The spectral profiles of the flowers were remarkably similar between pink and white corollas, with no ultraviolet signals detected in any of the samples [Supplementary Information—Fig. S1]. We recorded a total of 878 visits by ten bee species (*Apis mellifera*, *Anthophora plumipes*, *Anthophora dispar*, *Eucera alternans*, *Ceratina cucurbitina*, *Megachile deceptoria*, *Rhodanthidium sticticum*, *Chalicodoma lefebreyi*, *Chalidocoma sicula*, and *Osmia submicans*) and one beetle species (*Oxytyrhea funesta*) entering the flowers of the observed plants (Table 1). Altogether, the visiting insects assemblage significantly preferred to visit pink flowers than white flowers (Table 1), indicating support for our Hypothesis 1. This preference for pink flowers remained the same when splitting the visiting insects into two groups: territorial and non-territorial species (Table 1). In addition, three of the five species with sufficient sample size also showed a significant preference for pink flowers, whereas the other two
(Rhodantidium sticticum and Osmia submicans) did not prefer to visit any floral morph (Table 1). Flower colour preference by bee species, including results of the Jacob's D index and proportion of corolla colour preference, are shown in Supplementary Information—Figs. S2 and S3, respectively. The permutation-based multivariate analysis of variance (PERMANOVA) for the two corolla morphs showed similar relative frequencies of visits by bee species [Supplementary Information—Table S1].

## Genetic diversity

Plastid DNA spacer sequences displayed one point mutation each, A or C at the matrix position 137 in trnD-trnT and T or G at position 35 in trnS-trnfM, which resulted in three haplotypes (H1: AT, H2: AG, H3: CG). Haplotype data were obtained for 176 individuals (Fig. 2). Haplotype diversity statistics are reported in Table 2. All white-flowered individuals shared the same haplotype, while plants with pink flowers had the three haplotypes. SSR genotypes were obtained for 182 individuals and 13 loci; locus MSAT56 could not be interpreted as a biparentally inherited locus and was therefore excluded. Null alleles were observed at low frequencies

at three loci [Supplementary Information—Table S2] and one pair of loci (MSAT55-MAAC4) was in linkage disequilibrium (LD). In the full dataset of 13 loci, 62 alleles were observed, overall genetic diversity was $H_E = 0.462$, with little variation between categories based on flower colour or genetic (STRUCTURE-based) subpopulations (Table 3). The inbreeding coefficient $F_I$ was not significantly different from zero within categories or subpopulations, as expected for a self-incompatible plant species (Navascués et al. 2010). The full dataset of 13 loci was kept for all analyses because removing four loci to correct for null alleles and LD resulted in very similar genetic estimates [Supplementary Information—Tables S3, S4].

## Determinants of within-population genetic structure

Using MCMCglmm models with a single predictor matrix, we found that models with spatially defined matrices (GEO, ALT, or SUB) predicted inter-individual relatedness in the BDB population better, as shown by their lower DIC values than models including information on flower colour (COL) or vegetation at plant contact (VEG) (Table 4). This result

<table><caption>**Table 1.** Preference patterns of *A. charidemi* pollinators.</caption><thead><tr><th rowspan="2">Pollinators (T, territorial or floral constancy; N, non-territorial)</th><th colspan="2">Number of insect visits</th><th colspan="2">Visitation rate (Visits per flower per day)</th><th colspan="3">Jacob's D index of preference</th><th rowspan="2">P-value</th></tr><tr><th>White morph</th><th>Pink morph</th><th>White morph</th><th>Pink morph</th><th>White morph</th><th>Pink morph</th><th>Chi-square</th></tr></thead><tbody><tr><td colspan="9"><strong>HYMENOPTERA</strong></td></tr><tr><td><em>Apis mellifera</em> (T)</td><td>3</td><td>0</td><td>0.003</td><td>0.000</td><td>1.00</td><td>-1.00</td><td>4.95</td><td>0.026</td></tr><tr><td><em>Anthophora plumipes</em> (T)</td><td>3</td><td>6</td><td>0.003</td><td>0.003</td><td>-0.10</td><td>0.10</td><td>0.08</td><td>0.784</td></tr><tr><td><em>Anthophora dispar</em> (N)</td><td>9</td><td>42</td><td>0.008</td><td>0.022</td><td>-0.48</td><td>0.48</td><td>8.78</td><td>0.003</td></tr><tr><td><em>Eucera alternans</em> (N)</td><td>0</td><td>3</td><td>0.000</td><td>0.002</td><td>-1.00</td><td>1.00</td><td>1.82</td><td>0.177</td></tr><tr><td><em>Ceratina cucurbitina</em> (T)</td><td>9</td><td>59</td><td>0.008</td><td>0.032</td><td>-0.60</td><td>0.60</td><td>17.40</td><td>0.000</td></tr><tr><td><em>Megachile deceptoria</em> (N)</td><td>3</td><td>0</td><td>0.003</td><td>0.000</td><td>1.00</td><td>-1.00</td><td>4.95</td><td>0.026</td></tr><tr><td><em>Rhodanthidium sticticum</em> (T)</td><td>72</td><td>144</td><td>0.063</td><td>0.077</td><td>-0.10</td><td>0.10</td><td>1.80</td><td>0.180</td></tr><tr><td><em>Chalicodoma lefebreyi</em> (N)</td><td>133</td><td>308</td><td>0.117</td><td>0.165</td><td>-0.17</td><td>0.17</td><td>10.83</td><td>0.001</td></tr><tr><td><em>Chalicodoma sicula</em> (N)</td><td>3</td><td>0</td><td>0.003</td><td>0.000</td><td>1.00</td><td>-1.00</td><td>4.95</td><td>0.026</td></tr><tr><td><em>Osmia submicans</em> (N)</td><td>26</td><td>54</td><td>0.023</td><td>0.029</td><td>-0.12</td><td>0.12</td><td>0.94</td><td>0.332</td></tr><tr><td colspan="9"><strong>COLEOPTERA</strong></td></tr><tr><td><em>Oxythyrea funesta</em> (T)</td><td>1</td><td>0</td><td>0.001</td><td>0.000</td><td>1.00</td><td>-1.00</td><td>1.65</td><td>0.199</td></tr><tr><td>Territorial species</td><td>87</td><td>209</td><td>0.077</td><td>0.112</td><td>-0.19</td><td>0.19</td><td>8.81</td><td>0.003</td></tr><tr><td>Non-territorial species</td><td>171</td><td>407</td><td>0.151</td><td>0.217</td><td>-0.18</td><td>0.18</td><td>16.42</td><td>0.000</td></tr><tr><td>All species</td><td>262</td><td>616</td><td>0.231</td><td>0.329</td><td>-0.18</td><td>0.18</td><td>23.41</td><td>0.000</td></tr></tbody></table>

Total number of insects visiting the flowers of the observed plants, visitation rate calculated as number of insects visiting each flower per day, magnitude of Jacob's D index of preference and results of the goodness of fit testing the statistical significance of the preference pattern are shown. In bold those species and sample groups in which tests had a robust sample size. Pollinator's names are followed by pollination behaviour in brackets (T, territorial or floral constancy; N, nomadic) (see Vargas et al. 2017).

<table><caption>**Table 2.** Plastid DNA diversity. *n*, sample size, NA, number of alleles, $A_R$ ($k = 7$) allelic richness in sample of seven gene copies (individuals), $H_E$, gene diversity corrected for sample size.</caption><thead><tr><th>Category</th><th>n</th><th>% missing</th><th>NA</th><th>A<sub>R</sub> (k = 7)</th><th>H<sub>E</sub></th></tr></thead><tbody><tr><td>White</td><td>10</td><td>0</td><td>1</td><td>1</td><td>0.000</td></tr><tr><td>Pink</td><td>172</td><td>3.5</td><td>3</td><td>2.02</td><td>0.364</td></tr><tr><td>All</td><td>182</td><td>3.3</td><td>3</td><td>1.99</td><td>0.348</td></tr></tbody></table>
<table><caption>**Table 3.** Genetic diversity statistics of SSR markers for flower-colour-based categories, and topographic-genetic groups based on the STRUCTURE analysis.</caption><thead><tr><th>Category</th><th>n</th><th>missing %</th><th>NA</th><th>A<sub>R</sub>(k = 14)</th><th>H<sub>E</sub></th><th>F<sub>I</sub></th><th>P</th></tr></thead><tbody><tr><td>White</td><td>10</td><td>0.8</td><td>2.92</td><td>2.65</td><td>0.387</td><td>0.056</td><td>0.496</td></tr><tr><td>Pink</td><td>172</td><td>3.1</td><td>4.69</td><td>3.02</td><td>0.462</td><td>0.024</td><td>0.145</td></tr><tr><td>Subpop1</td><td>73</td><td>2.6</td><td>3.85</td><td>2.78</td><td>0.434</td><td>0.023</td><td>0.405</td></tr><tr><td>Subpop2</td><td>18</td><td>0.4</td><td>3.38</td><td>2.8</td><td>0.425</td><td>-0.004</td><td>0.918</td></tr><tr><td>Subpop3</td><td>16</td><td>5.3</td><td>3.46</td><td>2.93</td><td>0.461</td><td>-0.012</td><td>0.869</td></tr><tr><td>Subpop4</td><td>13</td><td>2.4</td><td>3.31</td><td>2.94</td><td>0.462</td><td>-0.005</td><td>0.972</td></tr><tr><td>Subpop5</td><td>62</td><td>3.7</td><td>3.92</td><td>2.94</td><td>0.462</td><td>-0.009</td><td>0.746</td></tr><tr><td>All</td><td>182</td><td>3</td><td>4.77</td><td>3.02</td><td>0.462</td><td>0.033</td><td>0.037</td></tr></tbody></table>

$n$, sample size; NA, number of alleles; $A_R(k = 14)$ allelic richness in sample of 14 gene copies (seven individuals); $H_I$, gene diversity corrected for sample size); $F_I$, individual inbreeding coefficient; $P$, probability that $F_I$ is different from zero based on 10 000 permutations.

<table><caption>**Table 4.** MCMCglmm modelling results for predicting inter-individual genetic relatedness.</caption><thead><tr><th>Predictor variables</th><th>DIC</th><th>deltaDIC</th><th>DICweight</th></tr></thead><tbody><tr><td>GEO绸</td><td>-25594.7</td><td>0</td><td>0.4044</td></tr><tr><td>GEO绸_SUB</td><td>-25594.6</td><td>0.11</td><td>0.3834</td></tr><tr><td>GEO绸_COL</td><td>-25593.5</td><td>1.29</td><td>0.2123</td></tr><tr><td>GEO</td><td>-25549.9</td><td>44.82</td><td>7.49E-11</td></tr><tr><td>GEO_绸</td><td>-25547.9</td><td>46.89</td><td>2.66E-11</td></tr><tr><td>ALT</td><td>-25423.1</td><td>171.66</td><td>2.14E-38</td></tr><tr><td>SUB</td><td>-25355.3</td><td>239.45</td><td>4.09E-53</td></tr><tr><td>NULL</td><td>-25344.3</td><td>250.48</td><td>1.65E-55</td></tr><tr><td>VEG</td><td>-25343.3</td><td>251.47</td><td>1E-55</td></tr><tr><td>COL</td><td>-25342.4</td><td>252.31</td><td>6.59E-56</td></tr><tr><td>COL_VEG</td><td>-25341.4</td><td>253.37</td><td>3.88E-56</td></tr></tbody></table>

Predictor variables listed include: NULL = null model, GEO = spatial distances, ALT = altitude differences, SUB = substrate differences (ravine vs. scree slopes/outcrops), COL = differences in corolla colour, VEG = differences with regard to vegetation at contact with the focal plant (with or without vegetation). For further details on model definitions, see Materials and Methods. The model with the lowest deltaDIC is the one that best explains genetic relatedness.

suggests a stronger support for our Hypothesis 2 (spatial determinants of genetic structure) than for our Hypothesis 1 (flower colour affects genetic structure). Specifically, the GEO model had a lower DIC than the COL model, showing that spatial distance explained genetic relatedness better than corolla colour. Of all models compared, including those with multiple predictor matrices, the GEO\_ALT model best predicted inter-individual relatedness (deltaDIC = 0). Predictive power remained similar when adding information on flower colour or substrate (GEO\_ALT\_COL or GEO\_ALT\_SUB).

## Characterization of within-population genetic structure

Permutation tests revealed a significant decrease of inter-individual kinship as a function of the logarithm of distance for both pDNA and SSR data ($P < 0.001$ in both tests), with a markedly stronger spatial genetic structure at pDNA, $Sp_n = 0.332$ (see haplotype distribution in Fig. 2B), than at SSRs, $Sp_n = 0.0198 \pm 0.0001$ SE (SSR kinship-distance plot in Fig. 3A). Spatial PCA confirmed the presence of a spatial genetic autocorrelation pattern through a significant G-test ($P < 0.001$), and revealed the absence of spatial genetic repulsion (non-significant L-test), with an eigenvalue of the first sPCA axis of eig.sPCA = 0.147. The geographic distribution

of scores on the first sPCA axis is illustrated in Fig. 3B, where patches of similar shading indicate genetic similarity. Spatial interpolation of scores on the first sPCA axis is shown in Fig. 3C, revealing the spatial distribution and discontinuities in genetic similarity.

Using the STRUCTURE software, $K = 4$ gene pools best explained the BDB data as this model had the highest lnP(P) and the highest deltaK [Supplementary Information—Fig. S4]. Spatial interpolation of cluster membership (Fig. 4) revealed the spatial clustering of related individuals, which led to the following definition of geo-genetic subpopulations (spatio-topographic groups of individuals, generally sharing genetic similarity): Subpop1, a large subpopulation dwelling in the narrow ravine gully in the eastern part of the population, where individuals were predominantly assigned to cluster K1 although other gene pools were also represented; Subpop2, a subpopulation on a rocky outcrop in the SE of the population with individuals mainly assigned to K2; Subpop3, a subpopulation on an ancient road cutting a scree slope in the north of the population, with individuals mainly assigned to K3; Subpop4, a subpopulation in the southern part of the population, representing mostly K4; and Subpop5, a large subpopulation corresponding to the uphill portion of the valley in the western part of the population, where all
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Isolation by distance pattern identified using nuclear SSRs. Kinship—distance regression; the 95% confidence interval for absence of isolation by distance is indicated with dotted lines (A). Individual scores on the first sPCA axis (B). Spatial interpolation of scores on the first sPCA axis (C).</figcaption></figure>

gene pools were represented (Figs. 2 and 4). GENELAND analyses [Supplementary Information—Fig. S5] were largely congruent with STRUCTURE analyses revealing a total of as much as *K* = 9 clusters, with larger patches of related individuals in the periphery of the population and a greater mixture in its centre, in the ravine.

The genetic differentiation between flower morphs was significantly different from zero with *F*<sub>ST</sub> = 0.067. The white morph occurred in STRUCTURE gene pools 1, 2, and 4, and in subpopulations 1, 4, and 5, and the pink morph occurred in all gene pools and subpopulations. Differentiation among the four STRUCTURE gene pools was *F*<sub>ST</sub> = 0.086, with pairwise *F*<sub>ST</sub> ranging from 0.069 to 0.120 (Table 5). Differentiation

among subpopulations was *F*<sub>ST</sub> = 0.045; pairwise *F*<sub>ST</sub> was the smallest, while significant, between the central Subpop1 and Subpop5, and the largest between peripheral Subpop2, Subpop3 and Subpop4 (maximum pairwise *F*<sub>ST</sub> = 0.106, Table 5).

## Inference of gene flow

In Migrate, the model with freely varying population parameters and migration rates was found to best explain the data, followed by the model with a single average migration rate and then the model with symmetric migration rates (Table 6). In the free model, theta estimates were similar for all subpopulations, with 1.45 ≤ Θ ≤ 1.79. The central
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Genetic clustering based on STRUCTURE, for K = 4 clusters, and definition of subpopulations for Migrate analysis (outlined with continuous lines) (A). Graphic representation of migration rate estimates between subpopulations (MIGRATE program) using subpopulation definitions based on STRUCTURE analysis (B). Subpopulation nodes that are net sources of migrants are indicated in green; net sinks are indicated in blue (B).</figcaption></figure>
<table><caption>**Table 5.** Genetic differentiation (pairwise $F_{ST}$) between colour morphs, between STRUCTURE gene pools, and between subpopulations defined based spatio-topographic position and STRUCTURE results. All values are significant with *P* < 0.0001.</caption><tbody><tr><td colspan="5">Pink-white</td><td>0.068</td></tr><tr><td>Gene pools</td><td>GP2</td><td>GP3</td><td>GP4</td><td></td><td></td></tr><tr><td>GP1</td><td>0.069</td><td>0.079</td><td>0.101</td><td></td><td></td></tr><tr><td>GP2</td><td></td><td>0.069</td><td>0.079</td><td></td><td></td></tr><tr><td>GP3</td><td></td><td></td><td>0.120</td><td></td><td></td></tr><tr><td>Subpopulations</td><td>Subpop2</td><td>Subpop3</td><td>Subpop4</td><td>Subpop5</td><td></td></tr><tr><td>Subpop1</td><td>0.045</td><td>0.079</td><td>0.065</td><td>0.020</td><td></td></tr><tr><td>Subpop2</td><td></td><td>0.106</td><td>0.071</td><td>0.063</td><td></td></tr><tr><td>Subpop3</td><td></td><td></td><td>0.095</td><td>0.047</td><td></td></tr><tr><td>Subpop4</td><td></td><td></td><td></td><td>0.052</td><td></td></tr></tbody></table>

<table><caption>**Table 6.** Natural log marginal likelihoods (lmL) and natural log Bayes factors (LBF) for the tested migration modes.</caption><thead><tr><th>Model</th><th>Thermodyn. lmL</th><th>Bezier lmL</th><th>LBF (Bezier)</th><th>Model prob.</th><th>Choice</th></tr></thead><tbody><tr><td>free</td><td>-1160950.4</td><td>-187250.3</td><td>0.0</td><td>1</td><td>1 (best)</td></tr><tr><td>mean_m</td><td>-1312621.9</td><td>-211995.0</td><td>-49489.3</td><td>0</td><td>2</td></tr><tr><td>symmetric_m</td><td>-1811656.0</td><td>-291591.5</td><td>-208682.3</td><td>0</td><td>3</td></tr><tr><td>fixtheta_sym_m</td><td>-2878131.5</td><td>-461468.2</td><td>-548435.7</td><td>0</td><td>4</td></tr><tr><td>fixtheta_mean_m</td><td>-2886602.7</td><td>-462812.0</td><td>-551123.3</td><td>0</td><td>5</td></tr></tbody></table>

The thermodynamic and Bezier log marginal likelihoods are reported as given by the migrate software. For model definitions, see materials and methods.

<table><caption>**Table 7.** Migration rate estimates from Migrate software for a model with freely variable population parameters and migration rates.</caption><thead><tr><th rowspan="2">Source population</th><th colspan="5">Receiving population</th><th rowspan="2">Θ</th></tr><tr><th>Subpop1</th><th>Subpop2</th><th>Subpop3</th><th>Subpop4</th><th>Subpop5</th></tr></thead><tbody><tr><td>Subpop1</td><td>0</td><td>21.6</td><td>19</td><td>19.3</td><td>14.6</td><td>1.45</td></tr><tr><td>Subpop2</td><td>1.6</td><td>0</td><td>9.6</td><td>9</td><td>2.2</td><td>1.68</td></tr><tr><td>Subpop3</td><td>1.6</td><td>4.1</td><td>0</td><td>8.2</td><td>2</td><td>1.54</td></tr><tr><td>Subpop4</td><td>1.6</td><td>6.9</td><td>5.5</td><td>0</td><td>2.1</td><td>1.79</td></tr><tr><td>Subpop5</td><td>6.4</td><td>15.8</td><td>19.7</td><td>15.9</td><td>0</td><td>1.52</td></tr></tbody></table>

Source populations for gene flow are represented on separate rows, receiving populations are represented in columns. The population parameter theta, Θ, is also represented.

subpopulations Subpop1 and Subpop5 were sources whereas the marginal subpopulations were mostly sinks for gene flow (Fig. 4B, Table 7). The number of migrants per generation received per subpopulation was always much larger than 1, a minimum rule-of-thumb threshold for conservation of genetic diversity (Mills and Allendorf 1996).

For inference of gene dispersal distances from spatial genetic structure, we need to assume that the population is in drift-dispersal equilibrium. Despite evidence of population substructure and thus a likely violation of this assumption, inter-individual kinship decreased near-linearly with the logarithm of distance, and the iterative estimation of dispersal distance converged for both SSRs and cpDNA, resulting in an estimate of gene dispersal distance of $σ_g = 21.707$ meters and an estimate of seed dispersal distance of $σ_s = 6.761$ meters. The inferred ratio of seed-to-gene dispersal distance was $σ_s/σ_g = 0.311$, inferred pollen dispersal distance was $σ_p = 29.172$ meters and the pollen-to-seed dispersal distance was $σ_p/σ_s = 4.315$. These estimates should be interpreted with

caution due to likely violation of assumptions, although the ratio estimates should be more robust.

## Discussion

In ecological and evolutionary contexts, the prevailing view is that pollinators have driven flower adaptations to novel pollinator niches since the origin of angiosperms (Rauscher 2008; Magallón and Vargas 2014; Peris and Condamine 2024). Within this framework, the *Antirrhinum majus-Bombus* system of flower-bee specialization has been proposed as an efficient driver of corolla evolution at the population level based on experimental (Dyer et al. 2007; Shang et al. 2011; Whitney et al. 2013) and field (Tastard et al. 2008) studies. However, over-generalizing the role of this pollination system in the evolution of any *Antirrhinum* species is problematic (see Romero et al. 2020). Indeed, most results focus primarily on the *Antirrhinum majus-Bombus terrestris* system, but they may not be representative of the interactions between other
*Antirrhinum* species and their pollinators (Vargas et al. 2017). Accumulating evidence suggests that not only pollinator behaviour patterns but also spatial isolation have played crucial roles in the diversification of snapdragons (Vargas 2009; Otero et al. 2021; Durán-Castillo et al. 2022). Fine-scale analyses that concurrently examine both spatial and environmental factors are lacking in the field. Our results from *A. charidemi* contribute new insights into the relative contributions of these elements to the early stages of differentiation. We showed that in the topographically heterogeneous Barranco del Drangoncillo Blanco population, the genetic structure was best explained by a model that included topographic position only (latitude, longitude, and altitude). The strongest patterns of differentiation identified in the population corresponded to partially isolated peripheral groups of individuals. Insect visitors displayed a preference for pink over white flowers and flower morphs were significantly differentiated. However, flower morphs were not associated with inferred gene pools or any geographically defined subpopulations. Overall, our results suggested primarily contributions of both spatial and pollinator-mediated processes to genetic differentiation in this *A. charidemi* population, with a stronger role of the former.

## Bee diversity provides multiple pollination opportunities

Our characterization of flower visitors on *A. charidemi*, including ten bee and one beetle species, is in agreement with previous field results. (i) Bees are the most important insect guild for pollen and nectar collection in the genus *Antirrhinum* (Vargas et al. 2017), (ii) large bumblebees are not necessary to prise open the occluded corolla (but see Davies et al. 2006), and (iii) medium to large bees appear to be responsible for effective pollination (Vargas et al. 2010) (Table 1). In our study, we recorded ten bee species visiting *A. charidemi*, which represents one of the richest faunas of pollinators associated with *Antirrhinum* species (Vargas et al. 2017). *Antirrhinum* species are clustered into two pollinator niches with respect to the diversity and functions of their visiting pollinators. The first niche contains *Antirrhinum majus* and seven other *Antirrhinum* species mostly visited by large nomadic long-tongued bees (essentially of the genus *Bombus*). *Antirrhinum charidemi*, together with six more snapdragon species, forms part of the second niche, which exhibits a higher diversity of pollinator species including long-tongued nomadic, long-tongued territorial, and short-tongued bees (Vargas et al. 2017). A rich pollination fauna entails a high number of bee species from different taxonomic families, body sizes and shapes, and behaviour patterns (Vargas et al. 2017), which can bring about conflicting selection pressures in relation to pollen transfer (Lankinen and Larsson 2009).

The fauna of flower visitors in *A. charidemi* showed a significant preference for pink flowers over white flowers (Table 1). Only five of the ten bee species displayed a significant preference for a flower morph (three preferred pink, two preferred white), but these results at the species level hide strong differences in visitation rates by the different species. The two main pollinators, *Chalicodoma lefevrei* (441 visits, significant preference for pink) and *Rhodanthidium sticticum* (216 visits, non-significant test), represented 75% of flower visits observed (878 in total). Preference for pink flowers was also pronounced in two additional bee species with a robust sample size, *Anthophora dispar* (Apidae, large, nomadic) and

*Ceratina cucurbitina* (Apidae, small, territorial), whereas significant preference for white colour relied on very small samples in two species, *Megachile deceptoria* and *Chalicodoma sicula* (Table 1). To assess whether such flower colour preferences contribute to genetic differentiation, it remains to be tested whether flower visits by these species lead to effective pollen deposition and seed set in the visited flowers. Further analyses at the species level, coupled with that at the individual level, over time might reveal even greater complexity in flower visitor behaviour between the two colour morphs.

The insect guild visiting *A. charidemi* flowers illustrates that bee movements could be significantly different among species and individuals. The two main flower visitors are Megachilidae species (*Rhodanthidium sticticum* and *Chalicodoma lefevrei*) which display foraging behaviours and pollination systems considerably different from those of *Bombus* species (typical of the *A. majus* pollinator niche). *Bombus* individuals deposit pollen by brooming into a specialised leg structure (corbicula); instead, Megachilid bees broom and deposit pollen onto the ventral surface of the abdomen (scopa) (Michener 1999). Pollination can also be different within the same bee species. For example, the territorial behaviour pattern observed for males of *Rhodanthidium sticticum* potentially results in more spatially restricted pollen flow than in drones or queens of nomadic *Bombus* species (see Fragoso and Brunet 2023). Conversely, females of *R. sticticum* showed nomadic movements among plants, which prevented from considering a single behaviour pattern for each bee species (Romero et al. 2020). Interestingly, short-distance gene flow by pollen mediated by *R. sticticum* males may be reinforced by the observed territorial, small bees (*Ceratina cucurbitina*) and the social honeybees (*Apis mellifera*), which have preference for floral constancy and pollination of nearby plants (Vargas et al. 2010; P. Vargas personal observation).

## Fine-scale spatial isolation by topography

We observed a strong effect of topography on genetic structure in the Barranco del Dragoncillo Blanco population, where the combination of spatial distance and terrain differences (altitude) were the strongest predictors of genetic structure (Table 4). This strong isolation signal is remarkable given the small population size (ca. 120 × 80m, e.g. about 1 ha), which is significantly divided into four-five distinct subpopulations (Fig. 4). Topography, particularly substrate type (classified as rock, wall, or soil), also played a crucial role in the differentiation of two additional *Antirrhinum* species endemic to south-eastern Spain (Barnbrook et al. 2023). In a population of the snapdragon *A. microphyllum* inhabiting rugged terrain in central Spain, Torres et al. (2003) also identified clusters of genetically similar individuals within small patches. To assess whether such differentiated subpopulations may represent early stages potentially leading to speciation, it would be interesting to examine their genomic patterns of divergence to search for any signatures of adaptive divergence.

The fine-scale spatial genetic structure observed in our population $Sp_n = 0.020$ ranks in the upper 30% of species with the strongest spatial structure and was typical for plants with gravity-dispersed seeds (mean $Sp = 0.028$) and animal-dispersed pollen (mean $Sp = 0.017$, Vekemans and Hardy 2004). Our estimate of the pollen-to-seed dispersal distance ratio of 4.31 m suggests a strong effect of gravity-mediated seed dispersal which contributed to isolating groups of tens
of individuals on a microspatial scale (1 ha). Our estimated gene dispersal distance was σ<sub>g</sub> = 21.7 meters, pollen dispersal distance σ<sub>p</sub> = 29.2 meters, and gene flow via seed dispersal distance σ<sub>s</sub> = 6.8 meters. These values should be interpreted with caution due to likely model violations. However, they appear further-ranging than estimates in other plants with similar or stronger *Sp*, e.g. *Silene acaulis*, *Medicago trunctula*, or *Hibiscus moscheutos*, which have gene dispersal estimates, σ<sub>g</sub>, that do not exceed 5 metres. The lower effective density of *A. charidemi* compared to those species can at least partially explain this difference (see Vekemans and Hardy 2004). In summary, topography appears to be key for microevolutionary patterns of differentiation in *Antirrhinum charidemi*. This result is congruent with that of Marrot and collaborators (2022), who found that natural selection acting on several vegetative traits varied both in magnitude and in direction in response to environmental variables at the fine spatial scale of even a meter.

Isolation by distance has played a main role in differenti-
ation also at larger scales. Phylogeography and genetic diver-
sity of *A. charidemi* across the narrow volcanic mountains of
Almería revealed a strong effect of geographic isolation on the
population genetic structure (Forrest et al. 2017). At a bio-
geographic scale, spatial isolation leading to speciation is also
interpreted because *A. charidemi* exhibits a well-supported
phylogenetic association with two species (*A. mollissimum*,
*A. rupestre*) also endemic to south-eastern Spain (Otero et al.
2021; Durán-Castillo et al. 2022).

Conclusions

Our population-scale study in the narrow endemic *Antirrhinum charidemi* identified a preference of most flower visitors for the pink flower morph, along with a significant genetic differentiation between the pink and white flower morphs. However, fine-scale spatial isolation in a topographically heterogeneous environment was the strongest driver of genetic differentiation in our study. These results on the very local scale support the hypothesis of spatial isolation as a prime force for evolutionary divergence in *Antirrhinum* previously verified on the biogeographic (macroevolutionary) scale and on the intra-specific (microevolutionary) scale in some species. We consider that spatial isolation followed by particular reproductive systems (self-incompatibility, exclusive reliance on bee pollination, pollen gene flow dependent on bee behaviour patterns, unassisted seed dispersal) and additional ecological conditions (topography, climate, substrate types, and altitude) have been particularly determinant for the microevolutionary patterns found in *A. charidemi*. Further research quantifying how spatial and ecological factors interact to drive differentiation of each *Antirrhinum* species is needed to determine their relative contributions in the early stages of differentiation.

Supplementary data

Supplementary data is available at AoB Plants online.

Acknowledgements

We acknowledge the Parque Natural de Cabo de Gata-Níjar
and Michelin for giving permission to work at the Barranco

del Dragoncillo Blanco (Barranco de Santa Bárbara) for 10
years. We thank Katharina Budde and Santiago González-
Martínez for help with field work, and Emilio Cano, Alan
Forrest and Mario Fernández-Mazuecos for help in the lab
and discussing results. A particular acknowledgement is
given to all the Ecoflor members who provided helpful re-
recommendations when one of the members (PV) presented the
*Antirrhinum* results of his 'second PhD' at the last 15 annual
Ecoflor conferences.

Author contributions

P.V. and M.H. designed the study, P.V. conducted flower vis-
itor observations, M.H. and P.V. conducted sampling and
mapping of the population, P.V., M.E., and J.M.G conducted
flower spectrophotometry analyses, M.H. produced the gen-
etic data. P.V., M.H., M.E., and J.M.G. analysed the data.
M.H. and P.V. wrote the first draft of the manuscript; all au-
thors interpreted results and critically commented and revised
successive versions of the manuscript.

Conflict of interest

None declared.

Funding

This study was supported by Banco Santander and Fundación General CSIC (FGCSIC) through a Proyecto Cero 2010 ('Do all endangered species hold the same value?: origin and conservation of living fossils of flowering plants endemic to Spain'). MH acknowledges a JAE-Doc 2009 postdoctoral fellowship of the Spanish Ministry for Science and Innovation.

Data availability

The data matrix with microsatellite and plastid DNA geno-
types as well as metadata is available on Zenodo (https://doi.
org/10.5281/zenodo.14198421). Flower visitor data is found
in Table 1.

References

Anderson MJ. A new method for non-parametric multivariate analysis of variance. *Austral Ecol* 2001;**26**:32–46.

Barnbrook M, Durán-Castillo M, Critchley J *et al*. Recent parallel speciation in *Antirrhinum* involved complex haplotypes and multiple adaptive characters. *Mol Ecol* 2023;**32**:5305–22. https://doi.org/10.1111/mec.17101

Beerli P. 2008. Migrate version 3.0: A maximum likelihood and bayesian estimator of gene flow using the coalescent. https://peterbeerli.com/ migrate-html5/download_version3/ (16 April 2025, date last accessed).

Bradley D, Xu P, Mohorianu II *et al*. Evolution of flower color pattern through selection on regulatory small RNAs. *Science* 2017;**358**:925–8.

Chhatre VE, Emerson KJ. StrAuto: automation and paralleliza-
tion of structure analysis. BMC Bioinf 2017;18:192. https://doi.
org/10.1186/s12859-017-1593-0

Comba L, Corbet SA, Hunt H *et al*. The role of genes influencing the corolla in pollination of *Antirrhinum majus*. *Plant, Cell & Environment* 2000;**23**:639–47. https://doi.org/10.1046/j.1365-3040.2000.00580.x
Crawford TJ. The estimation of neighbourhood parameters for plant populations. *Heredity* 1984;52:273-83. https://doi.org/10.1038/hdy.1984.29

Cueto M, Vargas P, Güemes J, Schwarzer H, Carrió E, Cabello J. 2008. *Antirrhinum charidemi*. Atlas y Libro Rojo de La Flora Vascular Amenazada de España. In: Bañares A, Blanca G, Güemes J, Moreno JC, Ortiz S. Adenda p.22-23. Dirección General de Medio Natural y Política Forestal (Ministerio de Medio Ambiente y Medio Rural y Marino)-Sociedad Española de Biología de la Conservación de Plantas Madrid.

Davies P, Compton SG, Macnair MR. The relative importance of floral color and scent to the pollination of *Antirrhinum majus*. *Am J Bot* 2006;93:399-408.

Del Valle JC, Alcalde-Eon C, Escribano-Bailón MT *et al*. Stability of petal color polymorphism: the significance of anthocyanin accumulation in photosynthetic tissues. *BMC Plant Biol* 2019;19:1-13.

Durán-Castillo M, Hudson A, Wilson Y *et al*. A phylogeny of *Antirrhinum* reveals parallel evolution of alpine morphology. *New Phytol* 2022;233:1426-39. https://doi.org/10.1111/nph.17581

Dyer AG, Whitney HM, Arnold SE *et al*. Mutations perturbing petal cell shape and anthocyanin synthesis influence bumblebee perception of *Antirrhinum majus* flower colour. *Arthropod Plant Interact* 2007;1:45-55.

Earl DA. 2012. Structure harvester: a Website and program for visualizing structure output and implementing the Evanno Method. *Conservation Genetics Resources* 4. Springer, 359-61.

Fenster CB, Armbruster WS, Wilson P *et al*. Pollination syndromes and floral specialization. *Annu Rev Ecol Evol Syst* 2004;35:375-403. https://doi.org/10.1146/annurev.ecolsys.34.011802.132347

Forrest A, Escudero M, Heuertz M *et al*. Testing the hypothesis of low genetic diversity and population structure in narrow endemic species: the endangered *Antirrhinum charidemi* (Plantaginaceae). *Bot J Linn Soc* 2017;183:260-70. https://doi.org/10.1093/botlinnean/bow002

Fragoso FP, Brunet J. Differential ability of three bee species to move genes via pollen. *PLoS One* 2023;18:e0271780. https://doi.org/10.1371/journal.pone.0271780

Glover BJ, Martin C. The role of petal cell shape and pigmentation in pollination success in *Antirrhinum majus*. *Heredity* 1998;80:778-84. https://doi.org/10.1038/sj.hdy.6883450

Gorospe JM, Monjas D, Fernández-Mazuecos M. Out of the Mediterranean region: worldwide biogeography of snapdragons and relatives (tribe Antirrhineae, Plantaginaceae). *J Biogeogr* 2020;47:2442-56. https://doi.org/10.1111/jbi.13939

Güemes J. *Antirrhinum* L. In: Benedí CC, Rico E, Güemes J, Herrero A. *Flora iberica*. Madrid: CSIC, 2009, 134-66.

Guillot G. Inference of structure in subdivided populations at low levels of genetic differentiation--the correlated allele frequencies model revisited. *Bioinformatics* (Oxford, England) 2008;24:2222-8. https://doi.org/10.1093/bioinformatics/btn419

Guillot G, Mortier F, Estoup A. Geneland: a computer package for landscape genetics. *Mol Ecol Notes* 2005;5:712-5. https://doi.org/10.1111/j.1471-8286.2005.01031.x

Gulluoglu C. 2016. *FlowMapper* (v0.2) *Plugin for Quantum GIS*. https://plugins.qgis.org/plugins/FlowMapper/version/0.2/ (16 April 2025, date last accessed).

Guzmán B, Vargas P. Systematics, character evolution, and biogeography of *Cistus* L. (*Cistaceae*) based on ITS, *trnL-trnF*, and *matK* sequences. *Mol Phylogenet Evol* 2005;37:644-60. https://doi.org/10.1016/j.ympev.2005.04.026

Hadfield JD. MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package. *Journal of Statistical Software* 2010;33:1-22.

Hardy OJ, Vekemans X. Spagedi: a versatile computer program to analyse spatial genetic structure at the individual or population levels. *Mol Ecol Notes* 2002;2:618-20. https://doi.org/10.1046/j.1471-8286.2002.00305.x

Herrera Maliani FJ. Datos sobre biología floral en la flora de Andalucía Morskaya. *Lagascalia* 1988;15:607-14.

Heuertz M, Vekemans X, Hausman J-F *et al*. Estimating seed vs. pollen dispersal from spatial genetic structure in the common ash. *Mol Ecol* 2003;12:2483-95. https://doi.org/10.1046/j.1365-294x.2003.01923.x

Jacobs J. Quantitative measurement of food selection. *Oecologia* 1974;14:413-7. https://doi.org/10.1007/BF00384581

Jiménez-López FJ, Arista M, Talavera M *et al*. Multiple pre-and postzygotic components of reproductive isolation between two co-occurring *Lysimachia* species. *New Phytol* 2023;238:874-87. https://doi.org/10.1111/nph.18767

Jombart T. Adegenet: A R package for the multivariate analysis of genetic markers. *Bioinformatics* 2008;24:1403-5. https://doi.org/10.1093/bioinformatics/btn129

Kantsa A, Raguso RA, Dyer AG *et al*. Community-wide integration of floral colour and scent in a Mediterranean scrubland. *Nat. Ecol. Evol* 2017;1:1502-10.

Kopelman N, Mayzel MJ, Jakobsson M *et al*. Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. *Mol Ecol Resour* 2015;15:1179-91.

Koski MH, Ashman TL. Macroevolutionary patterns of ultraviolet floral pigmentation explained by geography and associated climatic factors. *New Phytol* 2016;211:708-18. https://doi.org/10.1111/nph.13921

Landi M, Tattini M, Gould KS. Multiple functional roles of anthocyanins in plant-environment interactions. *Environ Exp Bot* 2015;119:4-17. https://doi.org/10.1016/j.envexpbot.2015.05.012

Lankinen A, Larsson MC. Conflicting selection pressures on reproductive functions and speciation in plants. *Evol Ecol* 2009;23:147-57. https://doi.org/10.1007/s10682-007-9227-z

Lechowicz MJ. The sampling characteristics of electricity indices. *Oecologia* 1982;52:22-30. https://doi.org/10.1007/BF00349007

Loiselle BA, Sork VL, Nason J *et al*. Spatial genetic structure of a typical understory shrub *Psychotria officinalis* (Rubiaceae). *Am J Bot* 1995;82:1420-5. https://doi.org/10.1002/j.1537-2197.1995.tb12679.x

Macior LW. Pollen-foraging behavior of *Bombus* in relation to pollination of nototribic flowers. *Am J Bot* 1967;54:359-64. https://doi.org/10.1002/j.1537-2197.1967.tb06930.x

Magallón S, Vargas P. Eudicotyledons: the greatest flower diversity in angiosperms. In: Vargas P. and Zardoya R. *The Tree of Life*. Sunderland, MA: Sinauer Associates, 2014, 157-66.

Mannino G, Gentile C, Ertani A *et al*. Anthocyanins: Biosynthesis, distribution, ecological role, and use of biostimulants to increase their content in plant foods—a review. *Agriculture* 2021;11:212. https://doi.org/10.3390/agriculture11030212

Michener CD. The corbiculae of bees. *Apidologie* 1999;30:67-74. https://doi.org/10.1051/apido:19990108

Michener CD. 2007. *The Bees of the World*, 2nd ed. Baltimore: Johns Hopkins University Press.

Mills LS, Allendorf F. The one-migrant-per-generation rule in conservation and management. *Conserv Biol* 1996;10:1509-18.

Moss SM, Zhou Y, Butelli E *et al*. Painted flowers: Eluta generates pigment patterning in *Antirrhinum*. *New Phytol* 2024;243:738-752.

Narbona E, Wang H, Ortiz PL *et al*. Flower colour polymorphism in the Mediterranean Basin: occurrence, maintenance and implications for speciation. *Plant Biol* 2018;20:8-20.

Narbona E, Arista M, Whittall JB *et al*. The role of flower color in angiosperm evolution. *Front Plant Sci* 2021;12:736998. https://doi.org/10.3389/fpls.2021.736998

Navascués M, Stoeckel S, Mariette S. Genetic diversity and fitness in small populations of partially asexual, self-incompatible plants. *Heredity* 2010;104:482-92. https://doi.org/10.1038/hdy.2009.159

Odell E, Raguso RA, Jones KN. Bumblebee foraging responses to variation in floral scent and color in snapdragons (*Antirrhinum: Scrophulariaceae*). *Amer. Mid. Naturalist* 1999;142:257-65. https://doi.org/10.1674/0003-0031(1999)142[0257:bfrtvii]2.0.co;2

Oksanen J, Simpson G, Blanchet F *et al*. 2025. vegan: Community Ecology Package. R package version 2.6-10. https://CRAN.R-project.org/package=vegan (16 April 2025, date last accessed).
Ornosa C, Ortiz-Sánchez FJ. 2004. Hymenoptera, Apoidea I. In: Ramos MA, editor. *Fauna Ibérica* vol. 23. Madrid: CSIC.

Otero A, Fernández-Mazuecos M, Vargas P. Evolution in the model genus *Antirrhinum* based on phylogenomics of topotypic material. *Front Plant Sci* 2021;12:631178. https://doi.org/10.3389/fpls.2021.631178

Peris D, Condamine FL. The angiosperm radiation played a dual role in the diversification of insects and insect pollinators. *Nat Commun* 2024;15:552. https://doi.org/10.1038/s41467-024-44784-4

Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. *Genetics* 2000;155:945-59. https://doi.org/10.1093/genetics/155.2.945

Qiao Z, Song X, Kong Y *et al*. Molecular mechanisms regulating ornamental traits and scent production in snapdragon (*Antirrhinum majus* L.). *Horticulture Advances* 2023;1:15.

QGIS.org. 2017. QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.org (16 April 2025, date last accessed).

Rausher MD. Evolutionary transitions in floral color. *Int J Plant Sci* 2008;169:7-21. https://doi.org/10.1086/523358

Richardson J. 2020. selectapref: Analysis of Field and Laboratory Foraging. R package version 0.1.2. https://CRAN.R-project.org/package=selectapref (16 April 2025, date last accessed).

Romero D, Ornosa C, Vargas P *et al*. Solitary bees (Hymenoptera, Apoidea) as connectors in pollination networks: the case of *Rhodanthidium*. *Apidologie* 2020;51:844-54. https://doi.org/10.1007/s13592-020-00765-2

Rothmaler W. Taxonomische monographie der gattung *Antirrhinum*. *Feddes Repert. Spec. Nov. Regni Veg* 1956;136:1-124.

Rousset F. *genepop'007*: a complete re-implementation of thegenepopsoftware for Windows and Linux. *Mol Ecol Resour* 2007;8:103-106. https://doi.org/10.1111/j.1471-8286.2007.01931.x

Rymer PD, Johnson SD, Savolainen V. Pollinator behaviour and plant speciation: can assortative mating and disruptive selection maintain distinct floral morphs in sympatry? *New Phytol* 2010;188:426-36. https://doi.org/10.1111/j.1469-8137.2010.03438.x

Schwarz-Sommer Z, Davies B, Hudson A. An everlasting pioneer: the story of *Antirrhinum* research. *Nat Rev Genet* 2003;4:655-64. https://doi.org/10.1038/nrg1127

Shang Y, Venail J, Mackay S *et al*. The molecular basis for venation patterning of pigmentation and its effect on pollinator attraction in flowers of *Antirrhinum*. *New Phytol* 2011;189:602-15. https://doi.org/10.1111/j.1469-8137.2010.03498.x

Smith SD. Pleiotropy and the evolution of floral integration. *New Phytol* 2016;209:80-5. https://doi.org/10.1111/nph.13583

Streisfeld MA, Young WN, Sobel JM. Divergent selection drives genetic differentiation in an R2R3-MYB transcription factor that contributes to incipient speciation in *Mimulus aurantiacus*. *PLoS Genet* 2013;9:e1003385. https://doi.org/10.1371/journal.pgen.1003385

Sutton DA. 1988. *A Revision of the Tribe Antirrhineae*. Oxford: Oxford University Press, pp. 575

Tastard E, Andalo C, Giurfa M *et al*. Flower colour variation across a hybrid zone in *Antirrhinum* as perceived by bumblebee pollinators.

Arthropod Plant Interact 2008;2:237-46. https://doi.org/10.1007/s11829-008-9046-3

Tavares H, Whibley A, Field DL *et al*. Selection and gene flow shape genomic islands that control floral guides. *Proc Natl Acad Sci USA* 2018;115:11006-11.

Torres E, Iriondo JM, Pérez C. Vulnerability and determinants of reproductive success in the narrow endemic *Antirrhinum microphyllum* (Scrophulariaceae). *Am J Bot* 2002;89:1171-9. https://doi.org/10.3732/ajb.89.7.1171

Torres E, Iriondo JM, Pérez C. Genetic structure of an endangered plant, *Antirrhinum microphyllum* (Scrophulariaceae): allozyme and RAPD analysis. *Am J Bot* 2003;90:85-92. https://doi.org/10.3732/ajb.90.1.85

Trunschke J, Lunau K, Pyke GH *et al*. Flower color evolution and the evidence of pollinator-mediated selection. *Front Plant Sci* 2021;12:17851.

van Buuren S, Groothuis-Oudshoorn K. Mice: multivariate imputation by chained equations in R. *Journal of Statistical Software* 45. *American Statistical Association* 2011;45:1. http://doc.utwente.nl/78938/

Van Oosterhout C, Hutchinson WF, Wills DPM *et al*. Micro-checker: Software for identifying and correcting genotyping errors in microsatellite data. *Mol Ecol Notes* 2004;4:535-8. https://doi.org/10.1111/j.1471-8286.2004.00684.x

Vargas P. 2024. Evolution in Snapdragons. https://www.rjb.csic.es/snapdragons/Paginas/. 21 July 2024.

Vargas P, Carrió E, Guzmán B *et al*. A geographical pattern of *Antirrhinum* (Scrophulariaceae) speciation since the Pliocene based on plastid and nuclear DNA polymorphisms. *J Biogeogr* 2009;36:1297-312. https://doi.org/10.1111/j.1365-2699.2008.02059.x

Vargas P, Ornosa C, Ortiz-Sánchez FJ *et al*. Is the occluded corolla of *Antirrhinum* bee-specialized? *J Nat Hist* 2010;44:1427-43. https://doi.org/10.1080/00222930903383552

Vargas P, Valente M, Blanco-Pastor JL *et al*. Testing the biogeographical congruence of palaeofloras using molecular phylogenetics: snapdragons and the Madrean–Tethyan Flora. *J Biogeogr* 2014;41:932-43.

Vargas P, Liberal I, Ornosa C *et al*. Flower specialisation: the occluded corolla of snapdragons (*Antirrhinum*) exhibits two pollinator niches of large long-tongued bees. *Plant biology* (Stuttgart, Germany) 2017;19:787-97. https://doi.org/10.1111/plb.12588

Vekemans X, Hardy OJ. New insights from fine-scale spatial genetic structure analyses in plant populations. *Mol Ecol* 2004;13:921-35. https://doi.org/10.1046/j.1365-294x.2004.02076.x

Whibley AC, Langlade NB, Andalo C *et al*. Evolutionary paths underlying flower color variation in *Antirrhinum*. *Science* 2006;313:963-6. https://doi.org/10.1126/science.1129161

Whitney HM, Bennett KMV, Dorling M *et al*. Why do so many petals have conical epidermal cells?. *Ann Bot* 2011;108:609-616. https://doi.org/10.1093/aob/mcr065

Whitney HM, Milne G, Rands SA *et al*. The influence of pigmenta-
tion patterning on bumblebee foraging from flowers of *Antir-
rhinum majus*. *Naturwissenschaften* 2013;100:249-56. https://doi.
org/10.1007/s00114-013-1020-y