1. The evolution of hummingbird pollination is common across angiosperm lineages throughout the Americas, presenting an opportunity to examine convergence in both traits and environments to better understand how complex phenotypes arise. We examine multiple independent shifts from bee to hummingbird pollination in the Neotropical spiral gingers (Costus) and use our data to address several common explanations for the prevalence of bee to bird pollination transitions.

2. We use floral traits of species with observed pollinators to predict pollinators of unobserved species and reconstruct ancestral pollination states on a well-resolved phylogeny. We examine whether independent transitions evolve towards the same phenotypic optimum and whether shifts to hummingbird pollination are associated with high elevation or climatic niche.

3. Traits predicting hummingbird pollination include small flower size, brightly-colored floral bracts, and the absence of nectar guides. We find many shifts to hummingbird pollination and no reversals, a single shared phenotypic optimum across hummingbird flowers, and no association between pollination and elevation or climatic niche.

4. Costus presents surprising findings compared to other plant clades. Hummingbird flowers are consistently smaller than bee flowers and primary flower colors are not predictive of pollinators. Moreover, hummingbird pollination shows no association with high elevation, as found in other tropical plants. 

Pollinator observations We collated both published and unpublished pollinator observation data (Table S1). For each species in the phylogenetic study, we recorded whether it had been observed, the total number of visits observed (if known), the proportion of observed visits that were by hummingbirds, the proportion of hummingbird visits that were by hermit hummingbirds, the proportion of visits matching the assigned pollination syndrome, and the data source. Pollinator visits were defined as flower visits in which there was likely contact between the visitor and the anthers and stigma. Thus, illegitimate visits in which nectar or pollen were removed without causing pollination were excluded. When a published source noted that a species was observed being pollinated by either hummingbirds or bees but did not present quantitative data, we recorded it as a single observation. We first asked whether the pollination syndromes assigned in taxonomic treatments predict the most frequent pollinators. This type of analysis was previously done using observations of 11 species (Kay & Schemske, 2003), but here we were able to include 28 Costus species. We also asked whether the proportion of visits matching the most frequent pollinator type differed between bee and hummingbird pollinated taxa using both ANOVA and phylogenetic ANOVA (aov.phylo; Harmon et al., 2008).

Floral traits -- We gathered a dataset of morphological and color traits from live plants, taxonomic publications, and well-documented photographs (Figure 1; Table S2). Live plants were measured in the UCSC Costus greenhouses and field with digital calipers. The flower was first removed intact from the inflorescence. Corolla measures were taken from the dorsal petal, and then the corolla was removed to reveal the floral tube, which consists of a single fleshy stamen and a labellum. Stamen exsertion was measured as the distance that the stamen tip protruded beyond the labellum and was negative when the stamen was inserted into the floral tube. Labellum and stamen size measures were taken from dissected flowers in which the parts were laid flat. The gynoecium was removed from the rest of the flower and the style was laid flat and measured from the ovary to the base of the stigma. 

We combined these with compatible measurements from monographs and protologues of Costus (Maas, 1972, 1977; Garcia-Mendoza, 1991; Maas & Maas-van de Kamer, 1997; Maas-van de Kamer et al., 2016) for species sampled in the Vargas et al. (2020) phylogenetic study. When taxonomic publications described a trait range for a species, we took the midpoint value.

In order to accommodate various data sources, we used categorical assignments of color based on human perception. Color descriptions were simplified into a series of categorical traits with a limited number of flower colors, including white, red, orange, yellow, green, and purple. We scored the color of various flower parts, including the exposed portion of the floral bract, the corolla (dorsal petal), the primary labellum color (avoiding nectar guides), the main portion of the stamen, and the stamen tip. We scored two binary characters for nectar guides -- whether or not the labellum had red stripes, and whether or not the labellum had yellow stripes. Red stripes typically occur around the periphery of the labellum, whereas yellow stripes are typically restricted to the center of the labellum and align with the path bees follow when entering the flower.

Our study omitted certain traits that may be important components of syndromes but are not noted in taxonomic treatments based on herbarium specimens or observable from photos. For example, we did not include nectar reward, floral scent, quantitative color reflectance, or differences in the texture or stiffness of flowers, which all may be important in pollinator attraction and efficiency.

Assembling the trait dataset -- We assessed whether we could combine data from measurements of fresh flowers with measurements from the taxonomic literature by amassing data for 21 species from both sources. For continuous traits, we analyzed Pearson correlations and linear model slopes and intercepts between data sources. For all traits except stamen and labellum lengths, correlations were high (> 0.8). Stamen and labella lengths had been measured from the apex of the ovary on fresh material but from the apex of the fused corolla tube in the monographs. Therefore, we corrected the monograph values of stamen and labellum length by adding corolla tube length. In addition, stamen exertion was not reported in the taxonomic literature, so we calculated it as stamen length minus labellum length. For categorical color traits, there were no mismatches in the data between sources. 

We combined data sources to arrive at species means for continuous traits and modes for categorical traits. In the greenhouse, when possible, we measured multiple flowers from a plant, multiple plants from a particular site, and plants from multiple sites representing a particular species. We combined the data by first averaging across flowers per individual, then across individuals for a site, and then across sites for a species. To combine data sources, we prioritized data from fresh material and high quality photographs and filled in missing data from the taxonomic literature, where possible. The combined dataset comprised 52 species, with 43 measured from fresh material and/or photographs and data for 9 taken from taxonomic literature.

Estimating pollinators and importance of floral traits using machine-learning -- We used machine learning algorithms (Random Forest Analysis, RF) to characterize pollination syndromes following Dellinger et al. (2019, 2021).  First, we trained and validated RF models on 28 Costus taxa with empirically documented pollinators.  We calculated 1000 RFs, each consisting of 500 trees with four characters tried at each split. To quantify the robustness of the RF predictions, we calculated the percentage of instances where a taxon was correctly classified in the training models (N=1000).  The importance of each floral character (N=17) for predicting pollinators was ranked by the mean decrease in Gini index over all 1000 RFs (standard deviation was also reported). This index is a measure of how important a floral character is for estimating pollinators across all the trees that make up a random forest. Finally, we estimated pollinators for 21 taxa missing pollinator observations using the 1000 RFs above.  

Reconstructing ancestral states and biases in direction of pollination shifts -- Using all taxa with observed and inferred pollinators (N=52), ancestral states of pollinator syndrome were estimated using maximum-likelihood methods under three separate models: equal rates, all rates different, and unidirectional rate from bee to bird with no reversals. Models were performed with the fitDiscrete function in geiger (Harmon et al., 2008). Model selection was then used to determine the optimal model using the aic.w function in phytools (Revell 2012). Using the optimal model, ancestral states were then estimated using 1000 iterations of Bayesian stochastic character mapping using the make.simmap function in the phytools package. Finally, to ensure that our results were not driven by inferred pollination syndrome (which differed somewhat from syndromes assigned in taxonomic treatments), this analysis was repeated using 1) syndromes assigned in taxonomic treatments and 2) only those tip taxa with pollinator observations (N=28).

Examining convergence in pollination syndromes – Because many of the 17 floral traits were highly correlated, we first reduced the dimensions of our trait data set. We estimated multivariate phenotypes using factor analysis of mixed data in the R package FactoMineR (Lê et al., 2008). We first imputed 34 missing data across 10 continuous traits in 52 taxa using the regularized algorithm and five components in the missMDA package (Josse & Husson, 2016). a. No data were missing for the 7 categorical traits. We then performed factor analysis on  the complete data set of 17 floral traits for 52 taxa using the FAMD function and visualized the output with the factoextra package (Kassambara & Mundt, 2020). We output the coordinates in the first ten dimensions for each species for use in further analyses. 

We visualized the evolution of floral traits by projecting the phylogeny on species’ values for dimensions 1 and 2 using Phytools (Revell 2012). We then asked whether floral dimension 1-10 were predicted by pollination syndrome using both ANOVA and phylogenetic ANOVA (aov.phylo; Harmon et al., 2008). Given the significant relationship between dimension 1 and pollination syndrome, we then assessed patterns of convergent evolution of dimension 1 under an Ornstein-Uhlenbeck process using the l1ou R package (Khabbazian et al., 2016; for a similar approach, see Smith & Kriebel, 2018). We note that l1ou does not require an a priori designation of where regime shifts (i.e., shifts in adaptive optima) occurred and instead estimates both the phylogenetic placement and magnitude of shift (function estimate_shift_configuration), and whether shifts are towards one or multiple regimes (function estimate_convergent_regimes). Convergence is inferred as independent shifts to the same regime.

Assessing whether pollination syndrome varies by elevation and climate -- Bees are thought to be less efficient pollinators at high elevations due to cool wet conditions whereas hummingbirds are efficient across all elevations (Cruden, 1972; Armbruster & McCormick, 1990; Dellinger et al., 2021). Therefore, we hypothesized that bee syndrome taxa should occupy lower elevations that are warmer and drier on average, and have lower variance in these attributes, than hummingbird syndrome taxa. To estimate the median elevation of each taxon, we used a set of previously published cleaned occurrence records (Vargas et al., 2020) and extracted elevation associated with each unique occurrence based on latitude and longitude (R package ‘elevatr’ function ‘get_elev_raster’; Hollister et al., 2020). This resulted in 3772 unique occurrences for 47 taxa (average 80 per taxa, range 4-500).  Climate niche estimates were used directly from Vargas et al., 2020, where principal component analysis summarized four climate variables for all unique occurrences: mean annual temperature, mean annual precipitation, temperature seasonality, and precipitation seasonality. Climate PC1 primarily captured variation in precipitation and the seasonality of temperature and precipitation, whereas PC2 primarily captured variation in temperature. Niche position of each species was then estimated as the mean of PC1 and PC2. To test whether median elevation, climate PC1 and climate PC2 were predicted by pollination syndrome we used phylogenetic ANOVA in three separate models (aov.phylo; Harmon et al., 2008). To test whether the variance in median elevation, climate PC1 and climate PC2 differed by syndrome we used three separate Levene’s tests (leveneTest function, car package in R; Fox & Weisberg, 2019). Ideally, we would account for shared ancestry in the levene’s test, however, this has yet to be implemented in R to our knowledge.