Population structure in Neotropical plants: Integrating pollination biology, topography and climatic niches

Abstract Animal pollinators mediate gene flow among plant populations, but in contrast to well‐studied topographic and (Pleistocene) environmental isolating barriers, their impact on population genetic differentiation remains largely unexplored. Comparing how these multifarious factors drive microevolutionary histories is, however, crucial for better resolving macroevolutionary patterns of plant diversification. Here we combined genomic analyses with landscape genetics and niche modelling across six related Neotropical plant species (424 individuals across 33 localities) differing in pollination strategy to test the hypothesis that highly mobile (vertebrate) pollinators more effectively link isolated localities than less mobile (bee) pollinators. We found consistently higher genetic differentiation (F ST) among localities of bee‐ than vertebrate‐pollinated species with increasing geographical distance, topographic barriers and historical climatic instability. High admixture among montane populations further suggested relative climatic stability of Neotropical montane forests during the Pleistocene. Overall, our results indicate that pollinators may differentially impact the potential for allopatric speciation, thereby critically influencing diversification histories at macroevolutionary scales.

explicitly evaluate how different pollinators affect population genetic differentiation, particularly so across multiple closely related plant species (and hence considering their shared macroevolutionary history; but see Barbará et al., 2007;Kramer et al., 2011). Such studies are urgently needed not only to better resolve macroevolutionary processes of angiosperm diversification, but also to choose appropriate management strategies in human-altered, fragmented landscapes and under current climate change (Castilla et al., 2017;Hadley et al., 2012;Toon et al., 2014).
Recently, attempts have been made to formalize theoretical predictions and to provide empirical evidence that pollinator mobility and foraging behaviour significantly impact population genetic structure (Wessinger, 2021). Specifically, less mobile pollinators with small foraging ranges are expected to generate more localized mating patterns between neighbouring individuals (Levin, 1979;Wessinger, 2021). This may result in lower heterozygosity and nucleotide diversity (particularly in small populations, Ness et al., 2010). In self-compatible species, increased inbreeding may be observed due to repeated visits to the same (geitonogamy) or closely related neighbouring plant individual(s) (Bezemer et al., 2019;Breed et al., 2015;Schoen & Clegg, 1984;Wessinger, 2021). This smallscale foraging behaviour may lead to strong isolation-by-distance even within populations, and particularly high genetic differentiation among populations (Opedal et al., 2016;Schmidt-Lebuhn et al., 2019;Wessinger, 2021). Small insects, nonflying vertebrates and territorial hummingbirds, with their systematic, repeated visitation of flowers on the same plant(s), fall into this category of "less mobile" pollinators (Medina et al., 2018;Schmidt-Lebuhn et al., 2019). Highly mobile pollinators, on the other hand, are expected to promote outcrossing and gene flow by visiting dispersed individuals over larger geographical distances, hence resulting in low population differentiation and high admixture among localities (Ballesteros-Mejia et al., 2016;Gamba & Muchhala, 2020;Hughes et al., 2007;Whelan et al., 2009). Moths and large bees (i.e., carpenter bees) may be considered relatively mobile pollinators since they forage over distances of several hundred metres and more (Castilla et al., 2017;Wikelski et al., 2010). They often adopt a traplining foraging strategy, where they systematically follow a route to visit flowers across spatially dispersed individuals, thereby transferring mixed pollen loads from different donors (Rhodes et al., 2017;Whelan et al., 2009;Wikelski et al., 2010). Traplining is also found in nonterritorial flying vertebrates (e.g., bats, passerines; Fleming, 1982;Tello-Ramos et al., 2015), and these may be considered as highly mobile pollinators since they often forage over distances of many kilometres and hence may more effectively bridge gaps between isolated populations (Bawa, 1992;Breed et al., 2015;Krauss et al., 2017).
These differences in mobility and foraging strategy of pollinators may further be increased by differences in their susceptibility to abiotic climatic conditions, and their ability to surmount geographical or environmental isolating barriers (Breed et al., 2015;Cruden, 1972;. Mountains, in particular, with their highly heterogeneous, rugged terrain, create dispersal barriers and strong abiotic (e.g., geological, climatic) environmental gradients across relatively small spatial and temporal scales (Cortés et al., 2018;Kisel & Barraclough, 2010;Nevado et al., 2018;Surget-Groba & Kay, 2013). Endothermic vertebrate pollinators may continue foraging across such variable environments since their activity is less affected by, for example, fluctuating, adverse weather conditions common in mountains (Cruden, 1972). Vertebrates may hence pursue their traplining routes across larger areas, and even surmount barriers imposed; that is, by habitat fragmentation or patches of unsuitable habitat (Breed et al., 2015). Ectothermic insect pollinators, on the other hand, are limited in their flower visitation activity to periods of sunny and warm weather (Cruden, 1972), probably further reducing their foraging ranges in mountains.
Besides present-day pollinator-mediated gene flow, a species' current population genetic structure also reflects climatic history and past adaptations to changing abiotic environmental conditions (Cortés et al., 2018;Helmstetter et al., 2020;Ramírez-Barahona & Eguiarte, 2013;Vasconcellos et al., 2019). Pleistocene glacial cycles, in particular, have led to repeated periods of population isolation, range shifts and habitat recolonization across the globe, leaving imprints in current population genetic diversity and differentiation (Valencia et al., 2010, Escobar et al., 2020Ornelas et al., 2013;Ortego et al., 2014). These patterns remain poorly understood in the tropics, with idiosyncratic responses documented among different groups of organisms and geographical regions (Ramírez-Barahona & Eguiarte, 2013;Vasconcellos et al., 2019). Establishing a solid understanding of the geographical and (historical) climatic factors, which may have shaped population structure, is therefore essential when evaluating the relative effect of different pollinators in driving population differentiation across plant species.
Here, we test the hypothesis that more mobile vertebrate pollinators promote higher levels of outcrossing and less population differentiation than less mobile bee pollinators. Using population genomic approaches based on genotype likelihoods, we contrast six plant species from the Neotropical tribe Merianieae (Melastomataceae).
Two of the six species are pollinated by bees, three are pollinated by mixed assemblages of vertebrates (i.e., traplining hummingbirds, bats and rodents) and one is pollinated by nonterritorial passerine birds ( Figure 1; Table S1; Dellinger et al., 2014Dellinger et al., , 2019a. We aim to understand whether mating patterns are indeed more localized (i.e., lower heterozygosity and nucleotide diversity, higher inbreeding, stronger within-population isolation-by-distance) in localities of bee-pollinated species than of vertebrate-pollinated species, and whether flying vertebrate pollinators consistently promote higher levels of gene flow across larger (geographical, topographical, environmental) distances than bees. Given the strong impact of abiotic environmental conditions on the activity of ectothermic bee pollinators, we further expect stronger effects of variation in climatic conditions in bee-than in vertebrate-pollinated species. For each species, we contrast two adjacent populations (<12 km apart) against three to four more distant populations, totalling 424 individuals from 33 localities (Figure 1; Table S1). Since all six species have small, dry, wind-dispersed seeds, gene flow attributable to seed dispersal should be comparable across the species (and potentially limited | 2267 DELLINGER Et aL. to short distances within populations given the dense structure of tropical forests, Kartzinel et al., 2013, Nazareno et al., 2020. Also, all species are self-compatible, but require pollinator visits to effect pollen release (A. S. Dellinger, unpublished;Dellinger et al., 2019a).
Finally, the six species differ in distribution ranges and ecosystems colonized (i.e., lowland rainforests, cloud forests, Figure 1). Our approach of directly addressing these differences through landscape genetics and (historical) niche modelling allows us to objectively evaluate the relative impact of geographical distance, topography, habitat suitability, climatic instability during the Pleistocene and environment on structuring population genetic diversity across a plant group of closely related species with distinct pollination strategies.

| Study group
Merianieae (Melastomataceae) are a Neotropical plant tribe of ~300 species, which has radiated recently in the tropical Andes . Bee pollination is ancestral and common both among lowland rainforest and cloud forest Merianieae, while shifts to vertebrate pollination (three times passerine, three times mixed assemblages of vertebrates) are restricted to cloud forest species (for detailed empirical pollinator observations, see Dellinger et al., 2014;2019b, 2021 and

| Sample collection, GPS coordinates and DNA extraction
We collected leaf and bark material in silica gel and sampled five to six localities per species in Costa Rica (2015) or Ecuador (2016, 2017; Figure 1; Table S1). To ensure comparability between species, we sampled at least two localities in close vicinity (<12 km distance) and the other localities at larger distances (>20 km; Table S2). On average, we sampled 13 (SD 3) individuals per locality across a distance of 1 km (Table S1). We recorded the exact location of each individual sampled to assess mating patterns within localities (Gelmi-Candusso et al., 2017) and calculated the centroid to obtain the average coordinate of each locality. We calculated Euclidean distances between localities based on these averaged coordinates in R (R Core Team, 2020).
We extracted genomic DNA from ~60 mg of dried plant material of 424 individuals using a CTAB protocol with sorbitol washing (Barfuss et al., 2016), RNAse treatment and subsequent clean-up (gDNA Cleanup Kit; Machery-Nagl). We quantified double stranded DNA content using a Qubit 3.0 Fluorometer with the dsDNA HS Assay Kit (Thermofisher) and only used samples with more than 6 ng µl -1 of DNA.

| Genome size estimation and RADseq library preparation
For three species (Ad. adscendens, Ax. costaricensis, M. phlomoides), we estimated genome sizes to select appropriate restriction enzymes for RADseq (restriction site-associated DNA sequencing).
We prepared fresh leaf material for three individuals per species in Otto's buffer (Otto et al., 1981) for propidium iodide flow cytometry (CyFlow ML flow cytometer, Partec; 532-nm/100-mW laser, Cobolt Samba, Cobolt AB) following Temsch et al., (2010). We used Solanum pseudocapsicum (1C = 1.29 pg DNA) as an internal standard to calculate 1C values for each run (three runs per sample) and calculated the mean 1C per species over all runs and individuals. Average genome sizes (1C value) for the three species were estimated as follows: Ad. adscendens 0.344 pg, SD 0.003; Ax. costaricensis 0.587 pg, SD 0.01; M. phlomoides 0.639 pg, SD 0.006. According to the Kew C-value database, these genome sizes are the second, third and fourth estimates for the family (Hanson et al., 2001).
In accordance with the relatively small genome sizes, we selected the restriction enzyme PstI (New England Biolabs) for single-digest RADseq (protocol modified from Paun et al., 2016).
We prepared eight RADseq libraries starting with 240 ng DNA per sample, pooling 72 samples per library, including also some repeats (using 120 ng DNA per sample). We ligated 300 nm P1 barcoded F I G U R E 1 Sampled localities (red dots, grey circles give population ID and clustering as estimated from admixture analyses [ Figure 3]) depicted on maps of habitat suitability (IBR Habitat , 0-unsuitable, 1-highly suitable). Red lines represent least-cost distances between study localities, and white areas on maps represent elevations below/above possibly suitable distribution ranges (for habitat suitability across full the distribution, see Figure S3). Per-species minimum/maximum F ST is given, as well as the pairwise F ST for the two adjacent populations. I-spline plots show the variable most strongly explaining population genetic differentiation as estimated through GDM (Table 2). Black bars represent 100 km, black outlines represent coastlines of Costa Rica and Ecuador, and map inserts represent each species' entire distribution range and the focal study areas. (a) Bee-pollinated Adelobotrys adscendens, lowland rainforests, Costa Rica; (b) bee-pollinated Meriania maxima, montane cloud forests, Ecuador; (c) hummingbird-bat-pollinated Meriania phlomoides, montane cloud forests, Costa Rica; (d) Meriania sanguinea with hummingbird-rodent-pollinated southern populations and hummingbird-bat-pollinated northern population, cloud forest-Páramo ecotone, Ecuador; (e) hummingbird-bat-pollinated Meriania tomentosa, cloud forests, Ecuador; (f) passerine-pollinated Axinaea costaricensis, cloud forests, Costa Rica adapters (150 nm if 120 ng input DNA) to the restricted samples at 16°C overnight. The P1 adapters included both index and inline barcodes, which were different from each other by at least three nucleotide positions. After P1 ligation, we sheared the DNA by sonication in a Bioruptor Pico (Diagenode) using two cycles (45 s "on," 30 s"off"; at 4°C) to obtain DNA fragments of 400 bp on average. After P2 adapter ligation, PCR (polymerase chain reaction) amplification and clean-up steps (using the MiniElute PCR purification Kit, Qiagen), we performed a final size selection for the range 220-850 bp using a 1.5% precasted dye-free cassette and a Pippin Prep

| Identification of RAD loci and variant filtering
We demultiplexed the data to sublibraries (index barcodes) with bamindexdecoder version 1.03 (included in Picard Illumina2Bam, available from http://gq1.github.io/illum ina2b am/), and qualityfiltered and further demultiplexed the reads to individual accessions (inline barcodes) using process_radtags.pl from stacks version 1.47 (Catchen et al., 2013). We removed low-quality reads with an uncalled base and corrected inline barcodes and RAD-tags with one mismatch, retaining only full-length reads (94 bp). We concatenated samples which had been sequenced twice. As there is no reference genome available for Merianieae, we followed Brandrud et al., (2020) to first create a pseudoreference using denovo_map. pl from stacks and followed up with a mapping approach to improve coverage and the recovery of loci. For each of the six species, we selected the 10 samples with the highest number of reads to be used for building a unique pseudoreference. We built several catalogues using different settings and chose the best following the r80 rule of Paris et al., (2017) for parameter optimization. We tested m (minimum number of identical reads required to create a stack) of 4 and 5, subsequently increased M (number of mismatches between stacks within an individual) starting from 1, and n (number of mismatches allowed between stacks between individuals) as n = M or n = M + 1 (Paris et al., 2017). For each setting, we recorded the number of tags retained with data for at least 80% of individuals and chose the settings m = 4, M = 1 and n = 2 which gave the maximum number of reliable polymorphic loci. We extracted the.fasta pseudoreference from the optimized catalogue by including polymorphic RAD loci that were present in at least 30% of samples and contained up to nine single nucleotide polymorphisms (SNPs) using export_sql.pl from stacks. We then mapped the raw reads of each accession to this pseudoreference separately using the mem algorithm of bwa version 0.7.12-5 (Li & Durbin, 2009), flagging shorter split hits as secondary (-M). Next, we sorted the resulting aligned.sam-files by reference coordinates and added read groups in the output.bam-files with picard tools version 2.18.17 (Wysoker et al., 2013). Finally, we performed a realignment around indels using the Genome Analysis Toolkit version 3.8.1 (gatk; McKenna et al., 2010).
The six species compared in this study are not equally closely related. Specifically, Ad. adscendens, the lowland bee-pollinated species, is part of a clade that diverged from core Merianieae (to which the other five species belong) ~25 million years ago .
Whereas the common pseudoreference we used for mapping should be appropriate for the five core Merianieae species, ascertainment bias may affect the inference for Ad. adscendens based on this reference. To address this issue, we tested whether the mapping rates for the six study species to the pseudoreference differ significantly from each other using ANOVA and Tukey tests as post-hoc tests.
Indeed, the mapping rate of Ad. adscendens was significantly lower than of all other species (F = 86.31, df = 5, p < .01; for pairwise comparisons see Table S3). We hence decided to create a separate pseudoreference for Ad. adscendens, based on its own accessions only, following the optimization procedure described above (final settings used are m = 5, M = 6, n = 6). Using this approach, we increased the mapping rate for Ad. adscendens from 31.7% to 35.1% (Table S4; the reference includes only polymorphic loci). We treated Ad. adscendens the same way as all other species in subsequent analyses unless otherwise stated.

| Mating patterns: Population genetic summary statistics
We measured genetic diversity within localities by calculating persite nucleotide diversity (θ π , i.e., average number of pairwise differences between sequences) using the genotype-free method implemented in angsd (Korneliussen et al., 2014). This approach yields more accurate estimates of population genomic parameters from medium to low coverage data like ours (Warmuth & Ellegren, 2019). We generated a theta file for each locality based on the locality-specific site frequency spectra (SFS). We intersected the theta files of all localities within each species in order to only compare SNPs that were shared by all localities following Peñalba et al., (2018). From these intersected sites files, we calculated nucleotide diversity using thetaStat as implemented in angsd (Maas et al., 2018).
We divided θ π estimates by the number of sites and tested for significant differences between localities using Kruskal-Wallis ANOVA and Dunn tests with Bonferroni correction for multiple comparisons (Dinno, 2017). For calculating heterozygosity, we calculated unfolded SFS per individual and divided per-individual heterozygosity by the total number of sites. Again, we used Kruskal-Wallis ANOVA and Dunn tests to evaluate differences among localities.
We calculated per-individual inbreeding coefficients (F) on genotype likelihoods as the degree of deviation from Hardy-Weinberg equilibrium (ngsf version 1.2.0-STD, Vieira et al., 2013). We performed two runs in ngsf, the first run to calculate reliable starting values of F per individual and a second performing a deeper search, allowing for a maximum of 1,500 iterations.
We used two measures of genetic distance to test whether populations were significantly different from each other. First, we converted pcangsd-derived covariance matrices to distance matrices (dist.from.cov, rwc version 1.11; Hanks, 2018). Second, we calculated pairwise genetic distances between all individuals of each species (ngsdist). We then tested for significant differences in genetic diversity between populations of the same species (adonis, vegan; Oksanen et al., 2019). We used pairwise.adonis (with corrected pvalue estimation) as a post-hoc test to assess which populations were significantly different. To test whether certain populations were more disparate (genetically variable) than others, we used betadisper and TukeyHSD as post-hoc test (results presented in the Supporting Information).
To assess population genetic differentiation, we calculated F ST values on genotype likelihoods starting from unfolded, pairwise SFS (core Merianieae) or folded SFS (Ad. adscendens) using realSFS in angsd.
We further calculated the admixture coefficient for each individual by estimating the likelihood of genetic clustering in the data (ngsadmix, Skotte et al., 2013). We randomly selected only one SNP per locus and used 10 random initializations to estimate admixture from K = 1 to K = n + 1 ancestral populations, n being the total number of sampled localities in each species. We compared the rate of change in the log-likelihood of different successive K-values (Evanno et al., 2005) to select the K-value best describing clustering in the data and used bar plots for visualization.

| Testing for the isolating effects of distance, terrain, current and historical climatic suitability and environmental niche
We used landscape genetic approaches to estimate the degree of connectivity among localities, assessing the relative impact of isolating barriers (IBR Terrain , based on topographic complexity), current climatic suitability (IBR Habitat , based on environmental niches), climatic suitability since the Last Glacial Maximum (LGM) (IBI, isolation-byinstability) and the environmental niche (IBE) on population genomic differentiation (F ST ). For each species' entire distribution range, we downloaded occurrence data from GBIF (Table S5; Chamberlain et al., 2021, Figure S1) and pruned the data using custom cleaning techniques (coordinatecleaner, Zizka et al., 2019), flagging records with equal longitude/latitude, zero coordinates, coordinate-country mismatches, records located in country centroids, in the sea or around GBIF headquarters, as well as duplicates and altitudinal outliers. To estimate IBR Terrain , we used raw elevation data (30 arc-sec resolution, earthenv.org/topography) and calculated topographic complexity using the Terrain Ruggedness Index (TRI, tri function, spatialeco, Evans, 2021). For each species, we drew convex hulls around GBIF occurrences with a 0.4° buffer (rgeos, Bivand & Rundel, 2020) and assessed each species' elevational distribution range (Table S5).
Within the convex hulls, we restricted the TRI layer to the plausible altitudinal distribution range of each species by removing areas 100 m above sea level below and above the lowest and highest GBIF occurrence (mask, raster; Table S5). To define a resistance cost surface based on TRI, we chose cost values to represent different proportions of the GBIF occurrences, with 1 (low cost) representing the central 70% of the GBIF occurrences of each species, 2 (moderate cost) representing the adjacent 5% quantiles (10%-15% and 85%-90%), 4 (high cost) representing the next 5% quantiles (5%-10% and 90%-95%), 8 (very high cost) representing the 0%-5% and 95%-100% quantiles, and 16 (extremely high cost) representing TRI values outside of the range covered by actual GBIF occurrences. Next, we calculated a transition object specifying a "knight and one-cell queen move" direction (transition, gdistance, van Etten, 2017) and correcting for map distortion (geoCorrection). Finally, we calculated all pairwise least-cost distances (paths with least-accumulative cost over cost surface) between the sampling localities (costDistance function).
To estimate IBR Habitat , we selected four bioclimatic variables which significantly impact pollinator activity and may hence be important in forming environmental barriers (Cruden, 1972;: Annual Mean Temperature (bio1), Mean Diurnal Temperature Range (bio2), Annual Precipitation (bio12), and Precipitation Seasonality (bio15; WorldClim 2.0, 30 arc-sec, Fick & Hijmans, 2017). We restricted these rasters as we did for IBR Terrain and estimated the abiotic climatic niche of each species using all GBIF occurrences using maxent 3.4.1 (dismo, Hijmans et al., 2020). We created 10,000 pseudoabsence points (the default number of background points used by maxent) from raster cells lacking GBIF occurrences (dismo, Hijmans et al., 2020). We set 80% of the GBIF occurrences of each species as the training data set and 20% to validate the models, using 500 iterations. We fit models 10 times and validated model fit using AUC (area under the receiver operating curve, values >0.75 indicate good fit, Swets, 1988) and TSS (true skill statistic, values range between −1 and +1; −1 to 0 no better than random; values >0.4 to 0.8 acceptable models, see Ornelas et al., 2019, Table S6;sdmtune, Vignali et al., 2020). In addition, we evaluated model performance using fourfold spatial-block cross-validation (blockCV, Valavi et al., 2019).
In this approach, the occurrence data set is split into four spatially distinct blocks, which are used as training data sets separately and subsequently used to test how well a model can be transferred to environmental space not used for calibration (Table S6). After assuring good model performance, we created raster maps of habitat suitability (0 0% suitable, 1 100% suitable) and averaged all 10 models.
To arrive with a resistance cost surface, we subtracted the climatic niche raster from 1 so that 0 means low cost and 1 means high cost.
We then calculated the transition object as for IBR Terrain and calculated all pairwise least-cost distances between the sampled localities.
To estimate IBI, we used the same four bioclimatic variables (bio1, bio2, bio12, bio15, at 2.5 arc-min resolution) as for the current climate, but derived from three different general circulation mod- For estimating habitat stability, we first estimated the current abiotic climatic niche of each species as for IBR Habitat , and then projected current habitat suitability to the LGM and mid-Holocene. We fit each model 10 times and averaged habitat suitability across runs and circulation models. Next, to derive a measure of habitat suitability through time, we overlaid the suitability maps of the three times (current, mid-Holocene, LGM) and averaged suitability of each raster cell. A raster cell providing 100% suitable conditions at all three times hence received a value of 1, while a raster cell providing unsuitable conditions at all three times received a value of 0. Again, we subtracted the suitability raster by 1 to arrive at a resistance cost surface and calculated all pairwise least-cost distances between localities of each species (see above; Figure S3).
Finally, we calculated IBE by extracting the raw values of the four bioclimatic variables (current climate) per population and calculating the Euclidean distances in environmental space between sampled localities using the dist function in R.
Naturally, the resistance surfaces we calculated for IBR Terrain , IBR Habitat and IBI represent refined geographical matrices and may thus be strongly correlated with the Euclidean geographical IBD matrices. We hence first ran separate Mantel tests (with 10,000 permutations) on these four matrices (each matrix standardized by the mean to account for differences in scale) to assess collinearity (García-Rodríguez et al., 2020). At least two of the four matrices were strongly correlated in all species (Table S7). We hence next constructed separate Mantel tests (with 10,000 permutations) to test for the effect of the respective distance matrix on normalized population genetic differentiation (F ST /(1 − F ST )). Then, for each species, we selected the distance measure with the highest R² and used multiple matrix regressions with randomizations (10,000 random permutations, Wang, 2013) to test for the relative impact of the respective distance matrix and IBE on population genetic differentiation using rglMMRR (popgenreport, Adamack & Gruber, 2014).  (Ferrier et al., 2007). For each species, we fit a GDM, specifying F ST as the response variable and IBD plus the resistance matrices and IBE as explanatory variables using the R-package gdm (Fitzpatrick et al., 2021). We plotted I-splines to visualize how magnitudes and rates of genetic differentiation vary with explanatory variables and estimated the (combination of) explanatory variable(s) best explaining genetic differentiation using backward elimination with permutation (see Ferrier et al., 2007, García-Rodríguez et al., 2020 for recent implementation). Briefly, the unique contribution of each predictor variable to total explained deviance is calculated, and the variable contributing least is discarded before fitting a new GDM. We ran 500 permutations until all variables retained in the final model made significant (p < .05), unique contributions to the explained deviance.

| Mating patterns and within-population IBD
Population nucleotide diversity (θ π ) across all sites ranged between 0.0006 and 0.0073 and was variable among species (Figure 2a; Tables S8 and S9). Per-species average estimates were lowest for bee-pollinated Ad. adscendens and highest for hummingbird-batpollinated Meriania tomentosa (Table S8). We detected significant differences in θ π both among adjacent and distant localities in all species (Tables S9 and S10).
Heterozygosity was most variable in bee-pollinated Ad. adscendens ( Figure 2b) and least variable in vertebrate-pollinated M. tomentosa (Figure 2b). We detected significant differences in heterozygosity among localities in all species except the vertebratepollinated Axinaea costaricensis and Meriania sanguinea (Tables S11 F I G U R E 2 Population genetic summary statistics (per-site nucleotide diversity, heterozygosity and disparity) for the six study species, colour coded according to pollination strategy. (a) Population nucleotide diversity differed significantly among species, but not among pollination strategies. (b) Individual heterozygosity differed significantly among species, but not among pollination strategies. (c) Disparity (variance in genetic diversity estimated on covariance matrix) did not differ significantly among pollination strategies. (d) Disparity (variance in genetic diversity estimated through ngsdist) did not differ significantly among pollination strategies. Colours indicate the different pollination strategies. Aa-Adelobotrys adscendens, Ac-Axinaea costaricensis, Mm-Meriania maxima, Mp-Meriania phlomoides, Ms-Meriania sanguinea, Mt-Meriania tomentosa and S12). Overall, only few localities differed significantly in heterozygosity (Table S12).

Per-individual inbreeding coefficients (F) were generally low for
all species ( Figure S4, Table S13). Eleven individuals of Ad. adscendens showed intermediate to high levels of inbreeding (F = 0.05-0.5), followed by M. sanguinea (one individual) while all other species had lower F values (Table S13).
Within localities, there was inconsistent IBD (Table S14)

| Population structure and disparity
Population structuring was variable between species and clustering mostly reflected geographical relationships (Figure 3; Figure   S5). Individuals of lowland bee-pollinated Ad. adscendens grouped into two distinct clusters (South-Western versus Northern Costa Rica, Figures 1a and 3a). While admixture analyses supported two ancestral populations with low admixture, co-ancestry estimates revealed shared co-ancestry between the two geographical clusters ( Figure 3a). In bee-pollinated M. maxima, individuals clustered into three distinct groups ( Figure 3b Figure S5).
In all species, at least some localities were significantly differentiated from each other (Tables S15 and S16), but adjacent localities were generally not differentiated (Figure 1). Genetic differentiation was weakest in passerine-pollinated Ax. costaricensis and highest in bee-pollinated Ad. adscendens. We detected most significant differences in genetic disparity among localities in Ad. adscendens (6/15 comparisons) and M. tomentosa (4-5/10 comparisons) and no differences in M. sanguinea and Ax. costaricensis (Tables S17 and S18).
Compared to the other species, montane bee-pollinated M. maxima showed lowest disparity within localities ( Figure S5).

| IBD, resistance, instability and environment
Pairwise population fixation indices (F ST ) confirmed patterns found in admixture and co-ancestry analyses and indicated moderate to high population structuring, ranging from 0.0178 (adjacent localities 1 and 2 of M. maxima) to 0.661 (distant localities 1 and 6 of Ad. adscendens; Figure 1; Table S2) we also detected a significant correlation with Pleistocene climatic instability (IBI , Table 1; Figures S3, S6 and S7). We did not detect any significant correlation with environmental distances (IBE, Figure 4, Table 1). These patterns were mostly confirmed by multiple matrix regressions, although effects of IBI/ IBR Habitat were not significant in Ad. adscendens and M. maxima, respectively (Table S19).
These results were mostly confirmed by GDM (Table 2; Tables S20 and S21 detail model fit). In bee-pollinated Ad. adscendens, we detected strong effects of IBD, IBR Terrain and IBE on F ST (Figure S8), although only IBD was retained as significant in the final model (Table 2). In bee-pollinated M. maxima, we found significant effects of current habitat suitability (IBR Habitat ) and IBD (Table 2; Figure   S9). F ST in hummingbird-bat-rodent-pollinated M. sanguinea and hummingbird-bat-pollinated M. phlomoides was significantly affected only by geographical distance (IBD, Table 2; Figures S10 and S11), while F ST in hummingbird-bat-pollinated M. tomentosa was most strongly affected by climatic instability since the LGM (IBI) and IBD (Table 2; Figure S12). In Ax. costaricensis, again, no variable significantly affected F ST (Table 2; Figure S13).

| DISCUSS ION
Here we tested the hypothesis that pollination by less mobile insect pollinators results in higher population differentiation and more localized mating patterns than pollination by highly mobile flying vertebrates. We explore this hypothesis across a sample of six closely related Neotropical plant species with disparate pollination systems. Assessing geographical, climatic and topographic isolating barriers, we detected consistently stronger isolating effects (particularly IBD) across localities of bee-than vertebrate-pollinated species (Table 1, Figure 4). These results support the hypothesis that more mobile pollinators (i.e., flying vertebrates) may connect populations more effectively than less mobile (i.e., [small] bee) pollinators (Medina et al., 2018;Wessinger, 2021). Even though large bees (pollinators of Meriania maxima) are generally considered as relatively mobile (Gamba & Muchhala, 2020;Wikelski et al., 2010), their flight activity is strongly reduced under adverse weather conditions in tropical mountains, probably limiting (large-distance) pollen dispersal . Within localities, we also detected IBD more frequently among bee-than vertebrate-pollinated individuals (Table S14). In contrast to the idea that less mobile pollinators would F I G U R E 3 Population structure, co-ancestry and F ST of the six study species. The greyscale bar plots (left) give admixture proportions according to the K-value best describing clustering in the data. In co-ancestry heatmaps (right, lower triangle), darker tones represent higher pairwise relatedness and stronger differentiation from other individuals; F ST values are shown on the upper triangle. Note that (d) and (e) were sampled across much larger distances than the other species. (a) Localities of bee-pollinated Adelobotrys adscendens form two clusters with low admixture proportions; heatmaps indicate relatively high relatedness within clusters, and considerable relatedness among individuals between clusters. (b) Bee-pollinated Meriania maxima with three clusters and very little admixture among clusters. (c) Three ancestral clusters in hummingbird-bat-pollinated Meriania phlomoides, with overall weak clustering, high admixture and high shared co-ancestry across localities. (d) In hummingbird-bat-rodent-pollinated Meriania sanguinea, the five hummingbird-rodent-pollinated localities from southern Ecuador clustered together with low admixture and shared co-ancestry across localities. They were significantly different from locality 4 (northern Ecuador, hummingbird-bat-pollinated, Table S15) and also showed floral adaptations to their different nocturnal (rodent/bat) pollinators (Table S1; Dellinger et al., 2019b). (e) Two ancestral clusters with low admixture in hummingbird-bat-pollinated Meriania tomentosa, with low shared co-ancestry between those clusters. (f) All localities were intermixed without clear clustering (best K = 2) in passerine-pollinated Axinaea costaricensis, with considerable admixture and shared co-ancestry among all localities promote more localized mating patterns, we did not find consistent differences in nucleotide diversity, heterozygosity or disparity according to pollination strategy (Figure 2). This suggests that, at least across small spatial scales (i.e., within localities, between adjacent localities), bees may be equally effective outcrossing pollinators as vertebrates (Castilla et al., 2020;Opedal et al., 2016, Schmidt-Lebuhn et al., 2019. Clearly, within-population genetic diversity may also be affected by factors other than pollinator foraging behaviour. The overall size of a population, for example, is generally positively correlated with genetic diversity (Hague & Routman, 2015).
Furthermore, differences in wind properties (i.e., wind speed, direction) among localities of the same species may impact the dispersal of fruits or seeds (such as, for instance, the tiny, wind-dispersed seeds of Merianieae), leading to pollen-dispersal-independent, TA B L E 1 Results of Mantel's tests on the impact of IBD, IBR Terrain , IBR Habitat , IBI and IBE on normalized population genetic differentiation (F ST ). Significant isolation by distance and/or resistance was observed in all species except Axinaea costaricensis, no IBE, the highest significant R² for each species is highlighted in bold, and significant values are given in italics

F I G U R E 4
Relationship between genetic distance (F ST ) and geographical distance (IBD), topographic barriers (IBR Terrain ), climatic instability (IBI), habitat suitability (IBR Habitat ) and climatic distance (IBE) in the six study species. We detected significant IBD/IBR/IBI in all species except Axinaea costaricensis, and no IBE in any species. Our results suggest that large geographical distances among localities of species pollinated by less mobile bee pollinators (blue) result in larger genetic differentiation than in species pollinated by mixed assemblages of (volant) vertebrates (red, yellow). Correlations of genetic distance and current habitat suitability and climatic variables were (mostly) not significant. Relationships are depicted for each species separately; each dot represents a pairwise population comparison spatial patterns of genetic diversity and clustering of individuals within populations (Mueller-Landau et al., 2008). We also want to emphasize that the set-up of our study with only two bee-vs. four vertebrate-pollinated species from two different biogeographical  Table S20, Figures S4-S10); Imp.-1/-2/-3 give models where, sequentially, the respective explanatory variable(s) contributing least to model fit were removed (indicated by empty cell). Note that models may fail to fit if relationships between F ST and (some) explanatory variables are weak (indicated by "/")  Figure S1). We may, hence, expect partly idiosyncratic responses to isolating barriers and (current and past) climatic habitat suitability. Following the expectation that mountain terrain generates strong physical barriers (Barbará et al., 2007), we did indeed detect significant IBR Terrain in all species except Ax. costaricensis (Table 1). It was, however, never recovered as a factor best explaining population differentiation (Table 2). While rugged mountain terrain may act as effective barrier even across small spatial scales (i.e., populations 3 and 6 in M. phlomoides, Figure 1), admixture was overall high among close (<12 km) localities in all species, suggesting considerable connectivity by both insect and vertebrate pollinators among localities (Figure 3; Castilla et al., 2017;Tello-Ramos et al., 2015).
Ectothermic insect pollinators may be affected more strongly by harsh abiotic climatic conditions (i.e., low temperatures, high precipitation, strong winds) than vertebrate pollinators (Cruden, 1972;. Accordingly, we found some isolating effects related to current climatic conditions ( isolating effects related to mountain topography and current climatic conditions. The southwestern localities are cut off from the northeastern localities by the central Costa Rican mountain range. These mountains feature moist and cool cloud forests, which are generally unhostile to smaller bees, the pollinators of Ad. adscendens (Dellinger et al., 2019c). Along the western coast, occasional dry habitats (i.e., Nicoya peninsula), on the other hand, represent unsuitable habitats for the moisture-adapted plants (Pröhl et al., 2010). In accordance with this, our niche models estimated the "least-cost" path connecting the southern and northern localities of Ad. adscendens through the southeastern lowlands along the Caribbean side of the high Central American mountains (Figure 1a; Patten & Smith-Patten, 2008). In M. maxima, the marked differentiation among localities 3 and 4 further correlates with environmental differences. These two localities were significantly differentiated from each other genetically, albeit only 20 km apart. Our habitat-suitability models indicated that locality 4 lies in a climatically less suitable area (Figure 1).
It is therefore possible that a difference in habitat suitability strongly reduced the probability of bee flight among these localities (given reduced bee activity under adverse montane weather conditions, , hence generating strong population differentiation. Interestingly, we also found individuals of M. maxima in locality 4 to differ morphologically (i.e., smaller flowers, nonrevolute leaf basis; A. S. Dellinger, pers. observ.). Whether these observed phenotypic differences are the result of random genetic drift, or a response to different selection pressures, remains to be investigated.
Among the vertebrate-pollinated species, we did not find isolating effects caused by habitat suitability (IBR Habitat ) or environment (IBE; Figure 1, Table 2). Indeed, cloud forests form a relatively continuous ecosystem, particularly on the eastern slopes of the Andes and Central American mountains (Balslev, 1988;Luebert & Weigend, 2014), and possibly provide continuously suitable habitats for the  (Escobar et al., 2020;Jørgensen & Ulloa Ulloa, 1994).
This demarcation line, part of the Amotape-Huancabamba zone, has acted both as a dispersal barrier for montane species as well as a corridor for lowland species (Trénel et al., 2008;Weigend, 2002).
Understanding how tropical plants reacted to Pleistocene climatic fluctuations, that is whether they retracted into refugia ("dryrefugia" hypothesis) or underwent down-and upslope migrations ("moist-forest" hypothesis), remains a major conundrum (Carnaval et al., 2009, Valencia et al., 2010, Ramírez-Barahona & Eguiarte, 2013, Escobar et al., 2020. If species retracted into small refugia, climatic instability through time should explain population genetic variation . In our study, we found significant associations between genetic differentiation and climatic instability in five species (  Figure S3). These results support the "moist-forest" hy- However, our sampling covered only parts of each species' distribution range and to fully understand whether and where cloud forests or foothill forests may have acted as refugial areas, it will be necessary to sample localities across the full distribution range of each species (e.g., see approach in Escobar et al., 2020;Helmstetter et al., 2020). Such broad-scale sampling will also be relevant to test more explicitly (i.e., using distance-based redundancy analyses; He et al., 2013) whether genetic variation of species differing in pollination strategy responds differently to topographic and climatic factors.
Together, our results are highly valuable in adding to the limited data available on the diverse processes shaping population genetic differentiation of tropical plants, including differences in pollination strategy. While we highlight that a wider sampling across bee-pollinated Merianieae is required to firmly establish the role of pollinators in promoting population differentiation, our result on stronger isolation among bee-than vertebrate-pollinated populations suggests a critical role of pollinator mobility in shaping population-level processes. Extrapolating to a macroevolutionary perspective, pollinator shifts are often invoked as "key innovations" spurring diversification (van der Niet et al., 2014). Potentially, pollinator shifts may also alter a population's susceptibility to isolation and, consequently, its potential for allopatric divergence (Kisel & Barraclough, 2010). Interestingly, in various Neotropical/Andean plant groups, shifts from bee to vertebrate (particularly hummingbird) pollination go hand in hand with increases in diversification rates (e.g., Lagomarsino et al., 2016, Serrano-Serrano et al., 2017. This is somewhat counterintuitive, however, if vertebrates are expected to buffer isolating effects among populations (see Serrano-Serrano et al., 2017 for a discussion on additional factors). Clearly, comparative studies, ideally focusing on small monophyletic plant complexes, and documenting both pollination ecology and population genetics of multiple populations across the landscape (e.g., Opedal et al., 2016), will be essential for resolving the relative contribution of pollinator shifts in spurring or limiting speciation through gene flow (Abrahamczyk et al., 2014;Kartzinel et al., 2013;Kisel et al., 2010). To date, we largely lack pollinator observations from multiple populations of the same plant species in the tropics, and hence know little about the variability of pollinator composition across a species' distribution range (pollinator mosaic; Gowda & Kress, 2013). Within our own data set, we have such information for only a subset of all populations (Table S1). While we documented the same functional groups acting as pollinators in most populations, we documented effective hummingbird and rodent pollinators in southern Ecuadorian populations of M. sanguinea, and hummingbird and bat pollinators in northern Ecuador. Obtaining such natural history information, in addition to population genomic data, will be key to adding a more realistic understanding of the processes governing speciation. This will ultimately help in addressing broad-scale questions such as how the tropics worldwide have become exceptionally biodiverse, but also in tracing similarities and differences in drivers of diversification across different tropical habitats and continents, namely Andean uplift in the Neotropics (Lagomarsino et al., 2016) and aridification in the African tropics .