1. Introduction

Biotic pollination is a complex phenomenon whose results depend on the modulation of intrinsic characteristics of plants and their pollinators, as well as of multiple environmental factors (i.e., climate, accompanying floral species) that impact this animal–plant interaction. Given their multidimensional nature, the pollination niche arises as a concept that simultaneously describes several factors that contribute to pollen transfer among plants and determine plant reproduction (Phillips et al., 2020). The importance of the pollination niche concept lies in the fact that it enhances the understanding of ecological dynamics (e.g., plant coexistence; Pauw, 2013) and evolutionary trajectories (e.g., diversification; Johnson, 2010; Gómez et al., 2015).

Fine-tuning changes in pollination niches can drive reproductive isolation in plants, and in the long-term, promote pollinator-mediated speciation (Stebbins, 1970; van der Niet et al., 2014). To date, the most striking example of pollination niche shift with evolutionary consequences has been reported in a North American Erythranthe (= Mimulus) group of species. In these, a change in the YELLOW UPPER gene locus, that controls floral carotenoid pigmentation (Liang et al., 2023), trawls a change in pollinators generating pollinator-mediated reproductive isolation between a species pair (Bradshaw and Schemske, 2003). In addition, pollination niche shifts can also be triggered by ecological conditions. For example, a poorly assembled or unpredictable pollinator environment may be replaced by a more stable and effective group of pollinators, leading to increased pollen transfer (Johnson and Steiner, 2000). This shift may confer adaptive advantages for plant species under the new pollinator environment (Kalisz and Vogler, 2003; van der Niet et al., 2014).

Plant species in pollinator-poor environments may experience intense interspecific competition due to limited pollinator access. Under this scenario, it is expected that the pollination niche, including floral traits, will shift towards strategies that reduce pollen limitation, generating in the long term more diverse plant communities (Pauw, 2013). Changes in pollinator assemblages can influence flower phenology, promoting increased access to diverse pollinators, extended flowering periods, and enhanced opportunities for effective pollen transfer among individuals within a population (Griffin and Barrett, 2003; Thomann et al., 2013). Shared pollinators may exert similar selective pressures on a group of species (Gómez et al., 2015) leading to patterns of convergent evolution towards certain model plant species present in communities (Johnson, 2010; Gómez et al., 2022). In this way, when a group of plants within a population undergoes pollination niche convergence towards a second (model) species, the population may experience assortative mating among its members (Schliewen et al., 1994; Savolainen et al., 2006), a key mechanism for sympatric speciation (Stebbins, 1970; van der Niet et al., 2014). The evolution of the pollination niche intensifies assortative mating potentially leading to rapid divergence and reproductive isolation (Yuan et al., 2013).

Evolutionary convergence and divergence can simultaneously exist within a population, generating a trade-off between these processes that governs shifts in pollination niches (Anderson and Johnson, 2006; Gómez et al., 2022). Within a population, agonistic displacements in the pollination niche toward a second (model) species likely represent a case of "common advertising display" (Roy and Widmer, 1999; Jersáková et al., 2012). In this, a group of co-occurring rewarding flower species increases collectively their reproductive success as a result of convergence. We explore this idea by studying a group of sympatric globular cacti (Eriosyce) with contrasting floral morphology in some species but high similarity in others (Fig. 1), suggesting potential convergence and divergence in pollination niches among sympatric taxa.

We studied four coexistent taxa of Eriosyce belonging to sections Neoporteria and Horridocactus, which inhabit a narrow coastal strip (ca. 15 km²) in the northern portion of the Mediterranean-type biome of Central Chile (Fig. 1). These cacti show variation in flower morphology, tepal color, floral phenology and pollinator assemblages (Guerrero et al., 2011, 2019a). We hypothesized that shifts in pollination niches drive the speciation of Eriosyce chilensis (section Neoporteria) from the hummingbird-pollinated E. litoralis (section Neoporteria) with a convergence toward the bee-pollination mode observed in E. mutabilis (section Horridocactus). This shift likely provided adaptive value by reducing pollen limitation constraints associated with hummingbird pollination, thereby enhancing reproductive success. To evaluate whether their pollination niches were divergent or convergent, we integrated multiple lines of evidence. First, we analyzed floral morphology, phenology, pollinator assemblages, and the visual perception of flower colors by their pollinators. These analyses allowed us to assess whether niche traits differed among taxa, supporting divergence, or overlapped substantially, indicating convergence. Additionally, molecular phylogenetics and population genomic analyses were employed to estimate the evolutionary trajectory of pollination niches, divergence times, genetic clustering, and detect gene flow as well as potential genomic adaptations to different pollinators. The presence of genetic structure consistent with limited gene exchange between sister species would support the hypothesis that shifts in pollination niches reinforce reproductive isolation. Species distribution modeling further revealed the extent of geographic and environmental overlap among the taxa. If species sharing habitats exhibit divergent pollination niches and limited gene flow, this would provide evidence of pollination-mediated isolation in a sympatric setting. Finally, measurements of pollen limitation and seed production offered direct insights into reproductive success and the potential fitness consequences of specific pollination strategies. Reduced pollen limitation and increased seed output in species with certain pollination modes would suggest adaptive advantages, offering insights into shifts in pollinator preferences. By integrating these approaches, our study tested the central hypothesis that shifts in pollination niche drive both divergence from the sister and convergence to the model species in sympatric Eriosyce, influencing reproductive isolation and potentially promoting speciation.

2. Material and methods 2.1. Studied species and site

We studied all Eriosyce species members that inhabit the narrow and exclusive habitat of E. chilensis (Hildm. ex K. Schum.) Katt and its variety E. chilensis var. albidiflora (F. Ritter) Katt. (hereafter, E. albidiflora), which occupies a shore platform hosting rocky outcrops, located between Pichidangui (32°08′S) and Los Molles (32°14′S) sandy bays (Fig. 1). Both E. chilensis/E. albidiflora are members of the section Neoporteria, which have fuchsia and pale-yellow flowers, respectively, and display open cup-like shaped flowers that are mainly visited by bees (Fig. 1; Guerrero et al., 2019b). A third Neoporteria member in the area is E. litoralis (F. Ritter) P.C. Guerrero and Helmut Walter that exhibits fuchsia tepals and tubular flowers exclusively visited by the giant hummingbird Patagona gigas (Guerrero et al., 2019b). E. litoralis has an extensive distribution, stretching along a discontinuous coastal rocky habitat from 29°56′S to 32°30′S. A Horridocactus member that inhabits this area, E. mutabilis (F. Ritter) P.C. Guerrero and Helmut Walter, display reddish to pale-yellow tepals and cup-like flowers visited by bees (Guerrero et al., 2019b); this cactus is restricted to a coastal strip at 32°S (Guerrero et al., 2011; Walter et al., 2024). These four Eriosyce taxa overlap in their distribution only between Pichidangui and Los Molles (Fig. 1).

The sandy areas, devoid of Eriosyce cacti due to the soil's lack of necessary compaction and strength, act as a geographical barrier, thereby contributing to the isolation of the taxa under study. The habitat for cacti within the area has an estimated maximum age of approximately 0.5–0.8 million years, dating back to when the marine terrace emerged due to seismic movements (Rodríguez et al., 2013). To the west, a mountain range extends from the coast to the Andes, acting as a physical barrier for coastal species. This creates an 'island-like' landscape, which boasts a high level of plant diversity, with 57% of its flora being endemic to Chile (Muñoz-Schick and de Trenqualye, 2011).

2.2. Geographical distribution

To assess the extent of distributional overlap and habitat isolation among the taxa, we constructed species distribution models (SDMs) for the four taxa included in this study. This analysis was designed to evaluate whether the species occupy allopatric or sympatric distributions, providing insights into their spatial coexistence and the degree of overlap in their climatic niches. By modeling the environmental variables that define each species' potential distribution, we aimed to infer potential geographic barriers to gene flow. Species occurrences were obtained from field observations conducted from 2014 to 2017 (Fig. 1) for Eriosyce chilensis (N = 113 occurrences), E. albidiflora (N = 77 occurrences), E. litoralis (N = 123 occurrences), and E. mutabilis (N = 131 occurrences). A subsampling approach where only one occurrence per 1 km² was included in the final dataset. To perform SDM, we used the 'BIOMOD2' package in R, including bioclimatic variables as predictors obtained from WorldClim (Fick and Hijmans, 2017). We implemented variance inflation factor (VIF) analyses to eliminate collinearity among variables. All models available in the package were tested and replicated ten times. The true skill statistic (TSS) and receiver operating characteristic (ROC) verification algorithms were applied, which identified the random forest (RF) model as the best fit. The projections for each RF model were reclassified into a binary raster format (0: absence; 1: presence) using a presence threshold of 0.75 (Guerrero et al., 2011). This reclassification and its overlap analysis were performed using the 'raster' package in R (Hijmans, 2022).

2.3. Taxon sampling for evolutionary inferences

Sanger and restriction-site–associated DNA sequencing (RAD-seq) were conducted to infer phylogenetic relationships among sympatric taxa. Between 2014 and 2019, we collected tissue using a nonlethal protocol (Guerrero et al., 2019a). Within the distribution range of Eriosyce chilensis, we only collected flowering individuals to avoid any taxonomic confusion with E. litoralis, which looks similar when not flowering. Sequences obtained were added to a comprehensive dataset covering the entire geographic distribution of Eriosyce (Guerrero et al., 2019a). For RAD-seq sequencing, we sampled E. chilensis (N = 17 individuals), E. albidiflora (N = 17), E. mutabilis (N = 18), and E. litoralis (N = 26) from the coastal zone where the studied species coexist. In this region, the species grow in proximity, often within just a few meters of each other, creating a highly sympatric setting. For E. litoralis, we included two samples from an allopatric population, located 80 km north of their overlapping range with E. chilensis. We also included samples (two per species) from the other Neoporteria members (E. castanea, E. senilis, and E. subgibbosa) and samples (one per species) of four other Eriosyce members (E. aspillagae, E. armata, E. aurata, and E. marksiana). In addition, to assess effects of hybridization and introgression, we included hand-made hybrids into RAD-seq sequencing. Hybrids were produced in 2017 transferring E. litoralis pollen to E. mutabilis, which were included in the germination assay (see below). Hybrid seedlings were maintained under nursery conditions for 18 months to ensure that sufficient tissue was available for DNA extraction.

2.4. Phylogenetic inferences based on Sanger sequences

To test our hypothesis regarding pollination niche shifts and their role in the speciation of Eriosyce chilensis, we estimated divergence times among the taxa under study. The calibrated ages of key nodes offer critical insights into the temporal context of pollination niche shifts and their association with speciation events. To achieve this, cacti samples were ground to a fine powder, and DNA was extracted using the DNeasy Plant Kit (Qiagen, Valencia, California, U.S.A.). These molecular data provided the foundation for constructing a time-calibrated phylogeny, enabling us to trace evolutionary relationships and assess divergence patterns among the studied taxa. We amplified three noncoding chloroplast genes (rpl32–trnL, trnH–psbA, trnL–trnF), one plastid gene (ycf1), and a nuclear gene (PHYC). DNA extraction and PCR followed previously described protocols (Guerrero et al., 2019a). Sequence data were assembled and edited using Geneious Prime 2020 (https://www.geneious.com/) and aligned using MUSCLE implemented in MEGA v.7.0.18 (Kumar et al., 2016). A microsatellite region in the ycf1 dataset was excluded (450 bp) because the alignment was ambiguous for that region. The best models of molecular evolution were determined using PartitionFinder v.2.1.1 (Lanfear et al., 2017) with the '–raxml' function (Stamatakis, 2014). The best models were GTRG (rpl32–trnL and trnL–trnF), GTRINVGAMMA (trnH–psbA and ycf1), and HKYGAMMA (PHYC). Molecular dating and phylogenetic tree inferences were performed simultaneously in BEAST v.2.6.2 (Bouckaert et al., 2014). Preliminary analysis yielded a polytomy at the nodes of the target species. We forced monophyly to obtain the calibrated age of these nodes. Bayesian searches were run for 20 million generations sampling every 1000 generations. For each analysis, we evaluated stationarity using Tracer 1.5 (http://beast.bio.ed.ac.uk/Tracer). For the age of the root and Notocacteae nodes, we relied on a previously published phylogeny of Cactaceae (Hernández-Hernández et al., 2014), extracting secondary calibration nodes and their 95% highest posterior density intervals (Fig. S1): 17.15 Ma (12.67–24.46 Ma) for the separation of Blossfeldia liliputana, 12.44 Ma (8.59–17.95 Ma) for Echinopsis chiloensis subsp. litoralis, and 8.78 Ma (5.54–13.03 Ma) for the core Notocacteae. The utilized fasta matrix is available at https://github.com/pabloguerrero-cmd/Eriosyce-pollination. Phylogenetic inference using Sanger sequences is based on five concatenated markers (4800 aligned sites).

To evaluate whether pollination niche shifts followed a convergent or divergent evolutionary trajectory within Neoporteria, we reconstructed ancestral states of pollination modes using a maximum likelihood (ML) framework combined with stochastic mapping on the consensus Bayesian phylogenetic tree of Eriosyce. Morphological differences between funnel-shaped and tubular flowers were used to inform the reconstruction of pollination modes in Eriosyce. Flower morphology in Eriosyce is a reliable predictor of pollination strategies. Funnel-shaped flowers, characterized by widespread tepals and the absence of nectary chambers are consistent with bee pollination (Cuartas-Domínguez et al., 2002; Guerrero et al., 2019b; Meriño et al., 2024). In contrast, the tubular flowers observed in E. litoralis possess large nectary chambers and exhibit an inward curvature of the inner tepals, a morphology that restricts access to many bee species and aligns with ornithophilous (hummingbird-pollinated) floral traits (Walter, 2008; Cádiz-Véliz et al., 2021; Meriño et al., 2024). We generated a discrete character matrix coding each species' pollination mode, according to flower morphology. We employed stochastic character mapping with the make.simmap function to generate a posterior distribution of character histories (Pennel et al., 2014; Revell 2024). These simulated histories allowed us to quantify the uncertainty associated with ancestral reconstructions by summarizing the frequency of inferred states across simulated trees. Ancestral states were visualized as pie charts at nodes and tips to illustrate the probabilities of pollination modes through evolutionary time.

2.5. Phylogenetic relationships and gene flow inference based on genomic sequence data

To study in greater depth the evolutionary relationships between coexisting taxa, and to evaluate reproductive isolation associated with the origin of Eriosyce chilensis and assess the potential role of reticulated speciation due to hybridization between E. litoralis and E. mutabilis, we conducted a genome-wide analysis using next-RAD genotyping-by-sequencing. Genomic DNA was prepared as next-RAD libraries by SNPsaurus LLC, following the protocol described by Russello et al. (2015). Genomic libraries were assembled in ipyrad v.0.9.62 (Eaton and Overcast, 2020) using the 'reference' approach by mapping the Carnegiea gigantea genome (assembly SGP5_Cgig_v1.3; Copetti et al., 2017). All loci shared by at least four samples after filtering were retained. Downstream analyses were performed using ipyrad analysis in the Python API, excluding samples from outgroups, undetermined taxa, or nonfocal taxa. Phylogenetic analyses were conducted on 70 samples from the four aforementioned taxa, as well as on samples of E. armata, E. aspillagae, E. aurata, E. marksiana and two samples of E. aconcaguensis. The code used for assembly and analysis is available in Jupyter notebooks (https://github.com/eaton-lab/eriosyce-rad).

The phylogeny based on RAD loci assembly used the window_extracter tool (ipyrad-analysis), allowing at most 40% missing data at any site. This yielded a matrix of 1.28M sites (76 samples; 63, 817 SNPs; 21% missing data). All handmade hybrids were excluded from these analyses. The alignment was inputted to RAxML v.8.2.12 (Stamatakis, 2014) using the -f algorithm to infer a maximum likelihood tree using the rapid hill climbing algorithm and 100 nonparametric bootstrap replicates under the GTRCAT substitution model.

Population structure and admixture were determined using PCA (ipyrad-analysis). At most, 10% missing data was allowed at any site. The largest dataset (79 samples; 30, 161 SNPs after filtering) included all samples from the Neoporteria and Horridocactus clades, as well as suspected or known hybrids between them. To reduce the effects of genetic linkage, only a single SNP was subsampled from each RAD locus (2626 unlinked SNPs). To represent variation from subsamples, PCAs were replicated over 25 subsamples of unlinked SNPs from the total data. Two additional datasets were compiled by either excluding samples from the Horridocactus clade, except for E. mutabilis (70 samples; 29, 172 SNPs; 3042 unlinked SNPs), or by excluding all samples from the Horridocactus clade (54 samples; 24, 051 SNPs; 3357 unlinked SNPs).

To test for genetic admixture, we used the Durand's D-statistic (ABBA-BABA tests; Durand et al., 2011) on pooled population genotypes based on the genome-wide relative frequencies of discordant SNP patterns. In each test, SNPs were filtered such that at least 50% of the samples in each population had data. A single SNP met the filter requirement, which was randomly selected from each RAD locus in each bootstrap replicate; loci were resampled with replacement to estimate variance in the test statistic. All tests included 30, 000–50, 000 SNPs from 4000 to 7000 loci, with approximately 500–1200 discordant SNPs. The results showed a Z-score representing the number of standard deviations for the observed D-statistic, which deviated from the mean of 1000 bootstrap replicates. These tests were implemented in ipyrad-analysis.

2.6. Pollination niche characterization

We characterized the pollination niche for each taxon based on three axes: floral morphospace, pollinator assemblage, and flower phenology. We determined the pollination niche overlap (NO_{i, j}) between species pairs i and j by integrating the three aforementioned niche axes (Geange et al., 2011) and examined each niche component independently. The NO_{i, j} is based on the average niche overlap between species pairs over different niche axis t, and it was estimated as:

In this, T is the total number of niche axes t; NO_{i, j, t} is the overlap between species pairs for the t axis. NO_{i, j} ranges between 0 and 1 for species pairs that completely differ and coincide in their niches, respectively. We reported NO_{i, j} as a percentage (0–100). Null distributions based on species label permutations (N = 100, 000 replicates) were generated to determine whether niche spaces differed between pairs of taxa (Geange et al., 2011). Significant results (P < 0.05) indicate that the species pairs do not overlap in terms of their niches.

The first niche axis is floral morphospace, which was reconstructed using 11 traits (Fig. S2) measured in the field during 2016 from one flower per individual using a digital caliper. These traits included corolla length and width, as well as the dimensions of inner flower tube organs (Fig. S2). These traits were compiled in a 'species × trait' matrix and used to perform a permutational analysis of variance (PERMANOVA) using a non-metric multidimensional scaling (NMDS) with the Euclidean distance method ordered in two dimensions (Clarke, 1993). PERMANOVA was used to estimate the source of variation in distance matrices, using F tests (Oksanen, 2010). The first NMDS dimension scores were used as the flower trait descriptor in the multidimensional niche characterization described in the previous paragraph.

The second niche axis is floral phenology. This was determined by monitoring individual plants (N = 20 per taxa) between July and December in 2016 and 2017. Annually, during field campaigns, we recorded the number of individuals with open flowers each day starting from the winter solstice (Southern hemisphere). Data were pooled for the studied years. We fitted a non-linear distribution model per species to describe the phenological process, using the Levenberg–Marquardt algorithm (Steer et al., 2019). The model includes three parameters: the maximum proportional change rate of the phenological phase (r), the data concentration (c), and the time–lag parameter (t). For each taxon, we estimated the 95% confidence intervals (CI_95%) of distribution. Statistical differences among taxa were determined using the criterion of no overlap among CI_95%. Pseudo-R2 values, which measure the variance proportion explained by models, were estimated (Steer et al., 2019).

The last niche axis is pollinator assemblages. This was characterized based on flower visitor observations (all animals that contacted the stigmas) during 30-min periods from dawn to dusk, within 25 m² plots over three consecutive days per species, between July and December of 2016 and 2017. These plots contained 12.7 ± 1.5 adult individuals each (range: 2–17 flowering individuals; N = 17 plots). We only considered plots where a single taxon bloomed. Sampling effort per taxon consisted of four observers, except for Eriosyce mutabilis (two observers). We accumulated a total of 33, 960 min of observations. Samples of visitors were collected, and specimens were identified by a bee expert (Dr. Luisa Ruz, Pontifical Catholic University of Valparaíso, Chile). We estimated visitation rates for each cactus by pollinator taxa (visits·hour⁻¹), and differences between cacti were assessed using PERMANOVA based on NMDS (K = 2), employing the Bray–Curtis distance method. Pairwise multilevel comparisons among groups after PERMANOVA were performed using the function pairwise.adonis() in R (Martinez, 2020). This analysis complements the previously described niche overlap analysis for pollinators. Part of these results (pollinator assemblage identity) was previously published (Guerrero et al., 2019b); in this work, we use the updated name of E. litoralis instead of E. subgibbosa, as used in an aforementioned study, following the new phylogenetic inferences reported here. In addition, we classified visitors into three functional groups (hummingbird, small bees, and large bees) and compared the number of visits for each group across cacti. Bees were classified based on their body length, using the median flower aperture in E. litoralis (9.46 mm; N = 41 flowers from independent individuals analyzed) as a threshold, which might act as a physical entrance barrier for bees. For analyses, we used a generalized linear model with a Poisson distribution error (Zuur et al., 2009).

Pollination niche analyses were performed using R 4.0.0 (R Core Team, 2017), with the packages vegan (Oksanen, 2010) and nlstimedist (Eastwood, 2020) used for morphospace and phenological characterizations, respectively. All data and codes for analysis are available at https://github.com/pabloguerrero-cmd/Eriosyce-pollination.

2.7. Modeling flower color appearance for pollinators

We modeled flower color appearance from the perspective of pollinators to investigate the underlying sensory mechanisms and the phenotypic selection exerted by pollinators, providing insights into the link between pollinator preferences and reproductive output. To achieve this, flowers from the four sympatric taxa were collected in 2020 and 2021 and kept fresh until spectral reflectance measurements were performed. For this, we illuminated tepals with a deuterium–halogen lamp (Avantes AvaLight-DH-S) through an optical fiber, while a second optic fiber connected to a spectrometer (Ocean Optics USD2000, Dunedin, FL, U.S.A.) collected the light reflected by the tepals. The spectrometer was calibrated using a reference standard (Avantes WS-2).

Flower color was plotted in the perceptual space of pollinators using the receptor noise limited (RNL) model, a method that estimates similarities between colors and has been experimentally validated for the honey bee (Vorobyev and Osorio, 1998; Vorobyev et al., 2001) and several bird species (Goldsmith and Butler, 2005; Lind et al., 2013). In the RNL model, colors are depicted as points in a chromaticity diagram, and discriminability is given by the Euclidean distance (dS) between points. Greater distances (high dS) represent higher chromatic contrast, predicting better discrimination between colors, while overlap or shorter distances (low dS) indicate higher similarity in color appearance, predicting lower or no discrimination (Vorobyev and Osorio, 1998; Vorobyev et al., 2001). To assess the similarity in flower color appearance between Eriosyce species for bees and birds, we estimated the overlap in the chromaticity diagram among the flower color loci occupied by each species (Stoddard and Stevens, 2011). A bootstrap analysis was used to calculate the estimated color distance and their confidence intervals between species pairs (Maia and White, 2018). Analyses were performed with the PAVO package in R (Maia et al., 2019).

2.8. Pollen limitation and reproductive output

For each taxon, we obtained seed production per fruit for unmanipulated flowers (control) and pollen-supplemented flowers. We studied 24–34 individuals per taxon. Pollen-supplemented flowers were hand-pollinated over two consecutive days with pollen from donors of the same species, located at least 5 m away from each other to avoid geitonogamy. Donors were two-day-old flowers, which was considered a suitable time to keep pollen viable (Guerrero et al., 2019b). Eight weeks after hand pollination, we collected fruits to count seeds. We estimated the pollen limitation (PL) index for each taxon as [1 – (Po/Ps)], where Po is the seed number obtained from control flowers, and Ps is the seed number for pollen-supplemented flowers (Larson and Barrett, 2000). When PL > 0, pollen limitation is assumed (Larson and Barrett, 2000). Statistical significance of PL was estimated by comparing its null distribution with μ = 0 (no pollen limitation) using a two-sided distribution nonpaired t-test (N = 100, 000 iterations).

The adaptive value for each taxon was estimated based on fruit set, number of seeds, and seed viability. The fruit set was registered as 0 (no seeds produced) or 1 (at least one seed within the fruit); seed production was quantified for each fruit. These variables were estimated from the control group described in the previous paragraph. Viability was estimated for seeds collected in the field between October and December 2017 (N = 20 individuals per taxon). Because Eriosyce litoralis produced a low number of fruits at the studied site, we used seeds collected from Pullally, a population located 20 km south of the study site. We also assessed viability in hybrid seeds obtained from hand-pollinated E. mutabilis (pollen receptor) and E. litoralis (pollen donor). We only studied these hybrids because the reverse crossbreeding did not produce seeds. Seeds from the same taxon were pooled, and germination was assessed in an incubator (Pi–Tec model: BIOREF–19L) under a cycle of 12-h daylight (10 ℃) and 12-h dark (5 ℃). Prior to the assays, the seeds were immersed in fungicide (1 g/L of Captan) to prevent mildew and then washed with distilled water. The assays were conducted in closed 100 mm Petri dishes with absorbent paper as the substrate, sowing 20 seeds per dish, with 20 replicates per taxon except for E. chilensis (N = 40 replicates). Germination was monitored for three weeks. For each variable, we tested differences among taxa using generalized linear models in the lme4 for R 4.0.0 (Zuur et al., 2009) with quasibinomial (fruit set), Poisson (seed number), and binomial (germination) error distributions.

2.9. Adaptive inferences based on genomic sequence data

The purpose of this analysis was to investigate whether the genomic divergence observed between Eriosyce chilensis (including E. albidiflora) and E. litoralis represents adaptive responses driven by differences in pollination strategies. By addressing this, we aimed to provide insights into the potential convergence of pollination niches and their adaptive value. Specifically, we sought to identify non-neutral (adaptive) genomic regions associated with shifts in pollination modes, testing whether these shifts are evolutionary responses to selection pressures imposed by pollinator preferences. To achieve this, we examined the influence of flower visitors on the adaptive genomic divergence between E. chilensis and E. litoralis using RAD-seq loci assembly. Merging E. chilensis/E. albidiflora datasets enabled us to explore species-level divergence and control for potential effect of intraspecific variation. In the first stage, we performed a genotype-environment analysis to detect candidate loci, i.e., SNPs potentially influenced by natural selection (Forester et al., 2018; Capblancq and Forester, 2021). To do so, we conducted a redundancy analysis (RDA) (Capblancq et al., 2018; Shryock et al., 2021). The biallelic assembled genotypic dataset was used as the response matrix (Y), while the dataset encapsulating the number of flower visits by bees and hummingbirds was used as the predictor matrix (X). The procedure to ascertain whether a locus could be classified as potentially adaptive was based on the methodology described in Capblancq et al. (2018), assuming a false discovery rate of 1%. Then, we evaluated the effects of pollinators, geographic distance, and genetic ancestry on the adaptive divergence of E. chilensis and E. litoralis. We implemented a distance-based RDA (db-RDA) using the subset of candidate loci retrieved in the previous step as the response matrix. As predictor matrices, we used the number of flower visits, the principal coordinates of neighborhood matrices based on sample georeferences (PCNMs; Legendre and Legendre, 2012), and the genetic ancestry dataset ascertained through population structure analysis. We used the LEA (Frichot and François, 2015) and vegan (Oksanen, 2010) packages for population structure and db-RDA analyses, respectively.

3. Results 3.1. Cacti distribution

Species distribution models (SDMs) confirmed that Eriosyce chilensis is restricted to a narrow coastal area, with an extent of occurrence of approximately 2.28 km² along 10.6 km of coastline. In contrast, E. litoralis has a broader geographic range, covering an estimated area of 104 km². E. mutabilis occupies a geographic area of about 41 km². Our results support the co-occurrence of Neoporteria members E. chilensis, E. albidiflora, and E. litoralis, as well as the Horridocactus member E. mutabilis (Fig. S3). The predicted distribution area of E. chilensis significantly overlaps with the ranges of E. mutabilis (79% overlap) and E. litoralis (87% overlap; Fig. S3). Despite some proximal occurrences of E. mutabilis and E. litoralis, our SDMs suggest that these species do not coexist outside the overlapping area shared with E. chilensis (Fig. S3).

3.2. Evolutionary relationships and genetic flow

Phylogenetic inferences strongly support that Eriosyce chilensis/E. albidiflora are sister taxa to E. litoralis, while E. mutabilis is more distantly nested within the Horridocactus clade (Fig. S1). Ancestral reconstruction in Eriosyce showed asymmetric transitions between bee and hummingbird pollination (Fig. S4), averaging 4.503 changes across 1000 trees, with hummingbird-to-bee transitions (3.163) outnumbering bee-to-hummingbird transitions (1.34). The analysis indicated that bee pollination accounted for 80.2% of the total evolutionary time, whereas hummingbird pollination represented 19.8%. E. chilensis and E. albidiflora exhibit bee pollination and are nested within the Neoporteria clade, which had evolved hummingbird pollination prior to the origin of E. chilensis/E. albidiflora.

All three core taxa in the study–Eriosyce chilensis (which includes E. albidiflora), E. litoralis, and E. mutabilis–were recovered as monophyletic groups according to SNP analysis that included several individuals per taxon (Fig. 2a). The two varieties of E. chilensis are placed in a mixed subclade (Fig. 2a); however, they form a differentiated cluster in PCA analyses from Horridocactus (Fig. 2b), as well as E. mutabilis (Fig. 2c), and E. litoralis (Fig. 2d). Time-calibrated phylogeny based on Sanger sequences suggests that the divergence time between E. chilensis and E. litoralis ranges from 0.25 to 0.97 million years ago (Fig. S1).

Gene flow among species, a potential key factor associated with phenotypic convergence, was first assessed using PCA on approximately 3000 unlinked SNPs with no missing data. This analysis revealed distinct genetic clusters for each taxon, suggesting limited gene flow between them (Fig. 2b–d). Indeed, we observed that two allopatric E. litoralis accessions (1075 and 1083; Fig. 2b), located 80 km north of their overlapping range with E. chilensis, form a group apart (Fig. 2a and b). All natural accessions, including three artificial F1 hybrids, appeared unique, as none closely resembled the hybrids (Fig. 2c and d). However, we found two potentially admixed accessions of E. chilensis (samples 1596 and 1597; Fig. 2c and d), which showed genomic similarities to hand-crossed F1 hybrids between E. litoralis × E. mutabilis (Fig. 2c), suggesting subtle gene flow that disrupts the tendency to form distinct genetic clusters.

Our analysis with D-statistics identified notable levels of gene flow (Table 1). Eriosyce chilensis and E. albidiflora share more alleles with the Horridocactus member E. mutabilis than with sympatric members of their sister species E. litoralis, as evidenced by negative D-statistics in Tests 1–4. Similarly, admixture was detected between E. mutabilis and E. chilensis/E. albidiflora, with Tests 3 and 4 indicating gene flow involving allopatric populations of E. litoralis. Analyses of admixture between E. litoralis and hybrids (both natural and hand-made) showed significant gene flow favoring E. litoralis, as highlighted in Tests 5 and 6.

3.3. Pollination niche and chromaticity diagrams

We observed that the narrow-distributed Neoporteria taxa, Eriosyce chilensis/E. albidiflora, exhibit a high degree of pollination niche resemblance to the Horridocactus member E. mutabilis, with niche overlaps of 60% and 62%, respectively (Table 2). In contrast, compared to their sister species E. litoralis, E. chilensis and E. albidiflora showed overlaps of 25% and 27%, respectively (Table 1). The pollination niches of E. litoralis and E. mutabilis exhibit the lowest similarity (9% overlap; Table 1). These results strongly suggest a convergence in pollination niches between E. chilensis (including E. albidiflora) with E. mutabilis.

Unidimensional analyses of niche axes showed that the flower morphospace of Eriosyce mutabilis is more similar to that of E. chilensis (76% overlap; P = 0.095) and E. albidiflora (85% overlap, P = 0.357; Fig. 3a) than to E. litoralis (< 2% overlap; P < 0.001; Table 2). Indeed, the differences between E. chilensis and E. litoralis were smaller than those observed between E. litoralis-E. mutabilis (7% overlap, P < 0.001; Table 2). Floral phenology showed that E. chilensis/E. albidiflora flower later in the spring, aligning more closely with E. mutabilis (Table S1), diverging from its sister species E. litoralis whose flowering peaks occur in winter (Fig. 3b). Despite this transition, E. chilensis/E. albidiflora have greater overlap with E. litoralis (41% and 10%, respectively) than with E. mutabilis (Table 2).

Regarding pollinators, we observed that Eriosyce litoralis is predominantly visited by the giant hummingbird, Patagona gigas (Fig. 4a and d). In contrast, E. chilensis, E. albidiflora and E. mutabilis are exclusively visited by bees (Fig. 4b and d). Niche overlap based on the number of visits showed a greater overlap between the Horridocactus species E. mutabilis with E. chilensis (93% overlap) and E. albidiflora (94% overlap; Table 2). Eriosyce chilensis and E. albidiflora showed strong differentiation from their sister species (18% overlap; Table 2). These results are consistent with PERMANOVA analysis (Fig. 4e), which shows significant differences in pollination assemblages among taxa (df = 3, SS = 37.04, F = 49.56, P < 0.001), with E. litoralis differentiating from the remaining sympatric cacti (Fig. 4e).

The distribution of flower colors in the chromaticity diagram of bees revealed high overlap of flower color loci between Eriosyce mutabilis and E. albidiflora, as well as between E. chilensis and E. litoralis (Fig. 5a), along with low chromatic contrast values between them (13% and 38% overlap, respectively; Fig. 5b). Tepal colors of E. albidiflora showed little to no overlap with those of E. chilensis and E. litoralis (0% and 0.3% overlap, respectively; Fig. 5c). Similarly, E. mutabilis flower colors exhibited no overlap with those of E. chilensis and E. litoralis (overlap 0%, Fig. 5a), with dS values high enough as to predict their discrimination by bees (Fig. 5b). In the chromaticity diagram of birds, color loci showed low overlap between flowers of E. albidiflora and E. mutabilis (overlap 3%; Fig. 5a), with dS values high enough as to predict their discrimination by birds (Fig. 5b). Tepal colors of E. chilensis and E. litoralis showed an overlap of 19% and low dS values (Fig. 5b), indicating that these flowers appear similar in color to birds (Fig. 5a). Flower colors of E. mutabilis showed little to no overlap with those of E. chilensis and E. litoralis (0% and 0.8%, respectively; Fig. 5a).

3.4. Adaptive value

We observed that seed production in both Eriosyce chilensis (PL = −1.1, CI_95% = −19.3 to 28.2) and E. albidiflora (PL = −1.2, CI_95%: −2.5 to 0.4) is not pollen limited, which contrasts with E. litoralis and E. mutabilis that exhibit pollen limitation (Table S2). Upon further evaluation of fruit set and seed production, we found that E. chilensis/E. albidiflora have a higher reproductive output (Fig. S5) compared to the other species (Table 3). Seed germination did not show significant differences among the studied species (Table 3; Fig. S5), although artificial F1 hybrids displayed significantly lower germination rates (Table 3).

We detected a significant influence of flower visitors as selective pressure agents on genomic divergence between Eriosyce chilensis and E. litoralis (1356 SNPs, RDA analysis). The distance-based RDA (db-RDA) revealed that the frequency of floral visitors had a significant effect on the adaptive divergence of this species pair, explaining 33.9% of the candidate SNPs (Table S3). Geographic distance was the next significant factor (32.4% of the candidate SNPs; Table S3), followed by genetic ancestry (27.0%; Table S3). The combined effects of pollinators, genetic ancestry, and geographic distance accounted for 19.5% of the adaptive divergence (Table S3).

4. Discussion 4.1. Pollination niches shifts

We observed a simultaneous convergence and divergence in pollination niches among sympatric Eriosyce cacti inhabiting a narrow coastal strip. Eriosyce litoralis, a widely distributed sister species of E. chilensis/E. albidiflora, experiences pollen limitation and relies on the unpredictable pollination services of the giant hummingbird (Patagona gigas) within its contact zone with E. mutabilis. Our findings suggest that the unreliability of hummingbird pollination and the similarity in flower traits between E. chilensis/E. albidiflora and E. mutabilis might reinforce shifts toward a common advertising display for bee pollination, offering these micro-endemic taxa a greater adaptive advantage while simultaneously promoting reproductive isolation from ornithophilous flowers. The mechanisms of pollination niche shifts identified in this study raise important questions about their potential role as a case of sympatric speciation in a continental but island-like habitat. Specifically, these findings underscore the possibility that E. mutabilis may act as a model species for pollinators, shaping the pollination strategies of sympatric taxa via shared visual and ecological cues.

The open funnel-shaped flowers of Eriosyce mutabilis closely resemble those of E. chilensis/E. albidiflora. Additionally, the flowering peak of E. chilensis/E. albidiflora shifts to spring, coinciding with the blooming period of E. mutabilis. Flowering phenology might be influenced by selective pressure from pollinators, especially bees, which are more active and abundant during warmer days (Heinrich and Esch, 1994). Genomic analyses show that the frequency of floral visitors has a significant effect on adaptive divergence between E. chilensis (including E. albidiflora) and E. litoralis. A cactus using bees as primary pollinators might be an advantageous reproductive strategy at the studied site. The most frequent flower visitors of E. chilensis, E. albidiflora and E. mutabilis are Anthrenoides sp., Chilicola sp., and Liphanthus sp., all are small bee species, suggesting similar pollination strategies among these cacti (Fig. S6). In contrast, E. litoralis shows a unique pollination pattern dominated by the giant hummingbird, aligning with its ornithophilous pollination syndrome. This suggests strong pre-pollination reproductive barriers (floral isolation) between E. litoralis and E. chilensis, which is consistent with the Grant–Stebbins model of pollinator-driven speciation (van der Niet et al., 2014). In bee and hummingbird pollination strategies, plants typically produce specific food resources—pollen and nectar, respectively—tailored to each group's ethological preferences (Willmer, 2011). For bees, an early spring presence of flowers with pollen, such as those provided by E. chilensis/E. albidiflora, is crucial as it becomes the first primary food source (Guerrero et al., 2019b).

The pollination niche agonistic displacement of Eriosyce chilensis/E. albidiflora toward E. mutabilis likely illustrates a common advertising display (Roy and Widmer, 1999; Jersáková et al., 2012). The distribution of color loci in the bee-perceptual space indicates low discrimination between E. chilensis–E. litoralis and E. albidiflora–E. mutabilis flower pairs, aligning with the higher gene flow observed among these species. This independent evidence further supports the hypothesis of speciation driven by bee-mediated differential perception and an ecological disjunction from the hummingbird-pollinated lineage, even though some gene flow between E. litoralis and E. chilensis may still occur, likely facilitated by small bees (e.g., Anthrenoides and Liphanthus, Andrenidae) capable of entering the tubular corolla of E. litoralis. Meanwhile, for E. albidiflora and E. mutabilis, larger bees such as Trichothurgus dubius (Megachilidae) may promote pollen exchange, particularly given the open flower morphology of these cacti.

The total absence of hummingbird visits to Eriosyce chilensis/E. albidiflora strongly supports that their pollination system is exclusively bee mediated. This is further reinforced by the advantage in seed production observed in E. chilensis/E. albidiflora, which exhibit reduced pollen limitation and copious seed production. A higher visitation rate and more efficient pollen deposition on these taxa might lead to increased fertilization of ovules, thereby generating the observed displacement in pollination niches (Pauw, 2013). Our study addresses a specific aspect of the eco-evolutionary processes that could explain this shift in the pollination niche. Seed germination is another crucial process influencing plant population dynamics and fitness, as it represents a critical stage in the life cycle, particularly in arid environments (Arroyo-Cosultchi et al., 2016). In our study, seed germination rates were generally consistent across taxa, indicating that post-reproductive barriers are almost absent at the intra-specific level. However, hand-made F1 hybrids (E. litoralis × E. mutabilis) exhibited significantly lower germination rates, suggesting some level of selection against hybrids and indicating that a hybrid origin of E. chilensis is unlikely. Furthermore, the contrasting pollinator assemblages between E. chilensis and E. litoralis left distinct genomic signatures, indicating disruptive selection driven by either bees or hummingbirds. Such reproductive shift, facilitated by assortative mating, are a clear but scarce example of how new species can emerge in the absence of geographic barriers (Servedio et al., 2011; Wenzell et al., 2023), driven by ecological factors such as pollinator preferences (Schluter, 2009; Schluter and Conte, 2009).

4.2. The origin of Eriosyce chilensis: an example of sympatric speciation?

The divergence of Eriosyce chilensis/E. albidiflora into an open-flowered ecotype, as showed by our study is characterized by distinct pollinator assemblages and phenological shifts, results in strong reproductive isolation from its sister species, E. litoralis. This process, likely initiated by the ethological exclusion of open flowers without nectar from hummingbird foraging behavior and further reinforced by bee within a shared geographical area. Co-occurrence is a key mechanism for sympatric speciation (Stebbins, 1970; van der Niet et al., 2014).

In our study, the spatial coexistence of Eriosyce chilensis/E. albidiflora with E. litoralis arise the question of sympatric speciation (Schliewen et al., 1994; Savolainen et al., 2006). The divergence observed between E. chilensis/E. albidiflora and E. litoralis fulfills all the criteria necessary for sympatric speciation, including reproductive isolation and ecological divergence, thereby supporting the hypothesis that sympatric speciation has occurred in this system (Coyne and Orr, 2004; Hopkins, 2013; Richards et al., 2019).

Firstly, there is geographical co-occurrence among the species, a finding supported by both direct observations and species distribution models. Our field observations indicate that they share the same habitat, specifically rock crevices in exposed marine terraces. Notably, ancestral state reconstruction reveals that the shift from hummingbird to bee pollination in Eriosyce chilensis/E. albidiflora coincides with their divergence from E. litoralis. This evolutionary change represents a reversion to an ancestral bee pollination system within Eriosyce. Such a reversal may have occurred due to genetic changes that retained the necessary genetic information to express bee pollination traits. This evolutionary history likely facilitated adaptation to a bee-mediated pollination niche, reinforcing reproductive isolation and contributing to speciation, emphasizing how shifts in pollination modes drive reproductive isolation in sympatric taxa. Consistently, the second criterion is supported by trait reproductive divergence driving reproductive isolation. Third, genetic analyses support the notion of a recent divergence (mid-Pleistocene divergence), with E. chilensis emerging as a distinct lineage within the shared habitat of crevices in rocky outcrops. This is further supported by observed genomic differences and limited gene flow between E. chilensis and E. litoralis, despite their physical proximity. Additionally, the recent uplift of marine terraces that created the habitat for E. chilensis approximately 500, 000 years ago (Rodríguez et al., 2013) matches the estimated divergence time between E. chilensis and E. litoralis, suggesting that this novel habitat promoted ecological conditions for new evolutionary trajectories.

Demonstrating that parental and derived species co-occurred at the time of origin poses a significant challenge in most sympatric speciation studies. Several factors in our study suggest an unlikely history of allopatry. First, both species (Eriosyce chilensis/E. albidiflora and E. litoralis) share specialized microhabitats in the crevices of rocky outcrops on marine terraces and coastal cliffs. Sandy bays to the north and south of the study area help circumscribe the species within the region, establishing an island-like habitat. This setting, encircled by sandy bays and towering mountains, forms a formidable barrier to dispersal, making it unlikely that the populations were outside each other's mating range, particularly considering the high movement capacity of the giant hummingbird Patagona gigas (Fernández et al., 2011; González-Gómez and Medrano, 2018). These factors collectively support a model of sympatric speciation driven by pollinator-mediated reproductive isolation rather than a scenario aligned with an allopatric speciation model.

4.3. Is floral mimicry possible in the origin of Eriosyce chilensis?

Floral mimicry is a key evolutionary mechanism in angiosperms that intertwines sexual selection with the sensory ecology and behavior of pollinators (Johnson and Schiestl, 2016). Generally, floral adaptive evolution by mimicry is triggered by the avoidance of pollen limitation, enhancing the fitness of mimic species compared to their non-mimic progenitors. This process fosters prezygotic reproductive isolation and evolutionary divergence (Schaefer and Ruxton, 2009). It is important to note that floral mimicry encompasses a broad spectrum, ranging from highly precise examples of specialized, specific interactions to more generalized models based on the functional relationships between interacting plant and animal species (Johnson and Schiestl, 2016). This explanation fits well with the studied system, where E. chilensis/E. albidiflora likely is the mimic species and E. mutabilis the model. Specifically, our system fits well with Müllerian mimicry, which involves both the mimic (E. chilensis/E. albidiflora) and the model species (E. mutabilis) providing rewards to pollinators, creating mutual benefits (Benitez-Vieyra et al., 2007; Johnson and Schiestl, 2016). In contrast, Batesian mimicry involves unrewarding flowers of the mimic species resembling those of a rewarding model species, effectively deceiving pollinators (Johnson and Schiestl, 2016). Although our results suggest a common advertising display between E. chilensis/E. albidiflora and E. mutabilis, determining floral rewards is necessary to strongly support the floral mimicry hypothesis. In addition, some authors suggest that experimental manipulations, such as removing mimic species, are required to demonstrate the fitness advantage of the mimic species (Johnson and Schiestl, 2016). However, the unique and threatened nature of the studied cacti (Faundez et al., 2013) made such experiments unfeasible.

Despite not having direct quantification of floral rewards or experimental manipulations to determine floral mimicry, our study supports the three conditions for floral mimicry (de Jager and Ellis, 2012). First, there is a specialized bee pollination niche of the model species (Eriosyce mutabilis). Second, receivers (pollinators) have a visual bias and might confuse the mimic (E. chilensis/E. albidiflora) with the model species, as they are perceptually similar. Third, the sexual selection exerted by receivers on the mimic favors a shift towards model species. In this case, the pollination niche of Eriosyce chilensis/E. albidiflora provides an adaptive advantage, as the pollination shift may result in an escape from pollen limitation, increasing seed production.

5. Conclusions

Our study provides compelling evidence that the geographical co-occurrence of Eriosyce chilensis/E. albidiflora with E. litoralis within the same habitat highlights the role of pollinator-mediated selection in driving their divergence instead of geographic isolation. In this system, we find that the shift from hummingbird to bee pollination in E. chilensis/E. albidiflora toward E. mutabilis, where similar floral traits and resources attract bees, enhances reproductive success and offers adaptive advantages to these micro-endemic taxa. In this eco-evolutionary context, floral mimicry and sympatric speciation may have concurrently shaped the evolutionary trajectory of sympatric Eriosyce cacti in a narrow coastal habitat. For instance, sympatric speciation requires mechanisms that initiate and maintain reproductive isolation without geographic separation, and for mimicry to evolve, shifts in the pollination niche must involve an imitator and a model that provide adaptive benefits. We acknowledge that further evidence is needed to fully unravel this undoubtedly complex evolutionary mechanism. However, the hypothesis of sympatric speciation through floral mimicry appears plausible, as suggested by ecological, spatial, and genomic evidence. Future research should incorporate functional genomics to uncover the genetic basis of floral divergence and mimicry and conduct behavioral studies of pollinators to further elucidate the mechanisms underlying reproductive isolation.

Acknowledgements

The authors thank S. Figueroa, G. Guerrero, and M. Rosas for help with fieldwork, and D. Craven, T. Hernández-Hernández, I. Larridon, J.J. Wiens, and Y. Yuan for providing comments on previous drafts of this manuscript. This work was supported by the Fondo Nacional de Desarrollo Científico y Tecnológico [1160583 and 1211441 to P.C.G.; 1240877 to G.O.C.], the Comisión Nacional de Investigación Científica y Tecnológica PIA [REDII 170031 to P.C.G. and G.O.C.], ANID PIA/BASAL [FB210006] to the Instituto de Ecología y Biodiversidad (IEB), ANID PIA/BASAL [PFB210018] to the Cape Horn International Center (CHIC). A.V.M. acknowledges the support of ANID/BASAL FB210006 by the Institute of Ecology and Biodiversity (IEB) with counterpart contributions from the Anglo-American Foundation. B.M.M. is grateful to the ANID Scholarship 20210673.

CRediT authorship contribution statement

Pablo C. Guerrero: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. Jaime Martínez-Harms: Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis. Mary T.K. Arroyo: Writing – original draft, Investigation. Deren Eaton: Writing – original draft, Software, Methodology, Investigation, Formal analysis. Beatriz M. Meriño: Methodology, Investigation, Formal analysis. Antonio Varas-Myrik: Writing – original draft, Formal analysis. Heidy M. Villalobos-Barrantes: Writing – original draft, Methodology, Formal analysis. Gastón O. Carvallo: Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Conceptualization.

Data availability

All data are available in https://github.com/eaton-lab/eriosyce-rad, https://github.com/pabloguerrero-cmd/Eriosyce-pollination and GenBank.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.03.001.