Adaptation to hummingbird pollination is associated with reduced diversification in Penstemon

Abstract A striking characteristic of the Western North American flora is the repeated evolution of hummingbird pollination from insect‐pollinated ancestors. This pattern has received extensive attention as an opportunity to study repeated trait evolution as well as potential constraints on evolutionary reversibility, with little attention focused on the impact of these transitions on species diversification rates. Yet traits conferring adaptation to divergent pollinators potentially impact speciation and extinction rates, because pollinators facilitate plant reproduction and specify mating patterns between flowering plants. Here, we examine macroevolutionary processes affecting floral pollination syndrome diversity in the largest North American genus of flowering plants, Penstemon. Within Penstemon, transitions from ancestral bee‐adapted flowers to hummingbird‐adapted flowers have frequently occurred, although hummingbird‐adapted species are rare overall within the genus. We inferred macroevolutionary transition and state‐dependent diversification rates and found that transitions from ancestral bee‐adapted flowers to hummingbird‐adapted flowers are associated with reduced net diversification rate, a finding based on an estimated 17 origins of hummingbird pollination in our sample. Although this finding is congruent with hypotheses that hummingbird adaptation in North American Flora is associated with reduced species diversification rates, it contrasts with studies of neotropical plant families where hummingbird pollination has been associated with increased species diversification. We further used the estimated macroevolutionary rates to predict the expected pattern of floral diversity within Penstemon over time, assuming stable diversification and transition rates. Under these assumptions, we find that hummingbird‐adapted species are expected to remain rare due to their reduced diversification rates. In fact, current floral diversity in the sampled Penstemon lineage, where less than one‐fifth of species are hummingbird adapted, is consistent with predicted levels of diversity under stable macroevolutionary rates.


Hummingbird adaptation is associated with reduced diversification in Penstemon
Supplemental Information 1. Supplemental methods.

DNA sequencing
We extracted DNA following a modified CTAB protocol and prepared MSG libraries using the restriction enzyme NdeI, exactly as described by (Wessinger et al. 2016). The new samples were multiplexed using unique barcoded adapters and sequenced in two batches, on the equivalent of one Illumina HiSeq 2500 lane at the University of Kansas Genome Sequencing Core. We have deposited new sequence data in the NCBI Sequence Read Archive (PRJNA548140). Sequence data for samples from our earlier study (Wessinger et al. 2016) is was previously deposited in the SRA (PRJNA314392).
We generated aligned sequence datasets from raw sequencing reads using ipyrad v0.7.22 (https://github.com/dereneaton/ipyrad). Reads were first demultiplexed according to barcodes, allowing one mismatch in barcode sequences. We then trimmed adapter reads using the cutadapt algorithm implemented in ipyrad (filter_adapters = 2). Reads were additionally trimmed (to a minimum length of 35 bp) if they contained low quality bases and reads were discarded if more than 5 bases had phred quality scores less than 20. We enabled reverse complement clustering by setting datatype to 'gbs'. We then clustered reads de novo within individuals using a minimum clustering threshold of 0.90 (90% similarity), the clustering threshold previously used for phylogenomic analyses in Penstemon (Wessinger et al. 2016). The minimum read depth was set to 6, the maximum number of uncalled bases per consensus sequence was set to 5, and the maximum alleles per consensus sequence per individual was 2 (all included species are presumed diploid). Sequences were then clustered across individuals using a clustering threshold of 0.90. Consensus loci were aligned if there were at least 20 taxa per locus, no more than 50 SNPs and 8 indels per locus, and no more than 8 shared heterozygous sites across samples. We generated two aligned data matrices. The first includes the full set of 120 Penstemon species ('full dataset'). The second includes the 104 crown clade Penstemon species ('crown dataset').
The aligned datasets are relatively sparse due to our low threshold for minimum taxa per locus. This allows homologous loci to be sampled, even in the face of large amounts of missing data. Causes of missing data include mutational disruption of restriction enzyme cut sites and stochasticity due low and uneven sequencing coverage. Even datasets with large amounts of missing data can be phylogenetically informative, particularly as the number of included samples increases (e.g., Huang and Knowles 2014;Eaton et al. 2017).

Phylogeny inference
We used IQ-TREE v1.6.2 (Nguyen et al. 2014) to infer a maximum likelihood (ML) phylogeny for each aligned and concatenated dataset. IQ-TREE performs favorably compared to other fast ML-based phylogenetic software (Zhou et al. 2017). Specifically, IQ-TREE is less prone to become trapped in local optima because it maintains a set of candidate trees during the analysis. We performed ten independent IQ-TREE analyses for each dataset using a GTR + Γ model of nucleotide substitution, 1000 ultrafast bootstrap approximations (Hoang et al. 2017), and otherwise default settings. We used the highestlikelihood tree of the ten replicate analyses for downstream analyses.
In addition, we estimated a species tree for each dataset using SVDquartets (Chifman and Kubatko 2014), a species tree approach that accounts for incomplete lineage sorting, implemented in Paup* v4.0a161 (Swofford 2003). We sampled 100,000 quartets from the data matrix and assessed support using 200 bootstrap replicates. SVDquartets does not currently estimate branch lengths, therefore we optimized branch lengths for the SVDquartets topology using a constrained tree search in IQ-TREE under a GTR + Γ model.
We rooted phylogenetic trees for the full dataset on the branch leading to P. newberryi (Wolfe et al. 2006) and we rooted trees for the crown dataset on the internal branch separating members of subg. Leptostemon from other members, a consistent relationship identified in the full dataset ML tree and in our previous phylogenomic study (Wessinger et al. 2016).

BAMM analysis
For the crown-ML and crown-ML-cropped trees, we estimated character-independent patterns of diversification using Bayesian Analysis of Macroevolutionary Mixtures implemented in BAMM v2.5.0 (Rabosky 2014). BAMM is designed to detect diversification rate heterogeneity across lineages and through time. We set the global sampling fraction to 0.8254 and the expected number of diversification rate shifts to 1. For each tree, we ran two replicate MCMC chains of 1 million generations, sampling every 100 generations, and assessed convergence and stationarity using Tracer (Rambaut et al. 2018). We analyzed results using the R package BAMMtools  to find the best model of shift configurations using Bayes factor model comparison.

Generation of ultrametric trees
We generated ultrametric trees for the crown dataset ML tree and the crown dataset SVDquartets tree using BEAST2 v2.5.1 (Bouckaert et al. 2014) to estimate relative divergence times on a constrained topology. We included all sites in a single partition and used a Γ site model with a GTR substitution model. We used a relaxed normal distribution clock model and a Yule model prior on diversification. For each of these analyses, we performed four MCMC chains, each consisting of 50,000,000 generations, sampling every 50,000 generations. We ran analyses on CIPRES (Miller et al. 2010). We discarded the first 10% of chains as burnin and combined replicate runs. We ensured combined sample reached stationarity using Tracer v1.7.1 (Rambaut et al. 2018) and used TreeAnnotator to exported Maximum Clade Credibility trees with median node heights for downstream analyses.

State-dependent diversification analyses implemented in RevBayes
We tested our prediction that transitions between pollination syndromes in Penstemon are asymmetric and that hummingbird syndrome lineages have reduced diversification rates by fitting state-dependent speciation and extinction (SSE) models. SSE models are designed to detect the effect of a trait on lineage diversification using a likelihood framework that jointly models character evolution and characterdependent birth-death processes. Model parameters include state-specific speciation rates (λi), statespecific extinction rates (μi), and transition rates between states (qij).
The Binary State Speciation and Extinction model (BiSSE) model considers whether a single observed character potentially affects diversification rate. We fit a series of four BiSSE models with states 0 (bee syndrome) and 1 (hummingbird syndrome). Model B1 has no constraints on parameter values, allowing the six BiSSE model parameters to vary (λ0 ≠ λ1, μ0 ≠ μ1, q01 ≠ q10). Model B2 constrains transition rates to be equal (q01 = q10). Model B3 does not allow transitions from hummingbird to bee syndrome (q10 = 0). Finally, model B4 is the null model of a single diversification rate that constrains speciation and extinction rates to be equal (λ0 = λ1, μ0 = μ1).
In addition to these BiSSE analyses, we performed Hidden State Speciation and Extinction (HiSSE) (Beaulieu and O'Meara 2016) analyses. The BiSSE model assumes no additional sources of diversification rate heterogeneity beyond the modeled binary character (Maddison et al. 2007). However, diversification rate variation across the tree may be independent of the observed character and instead be driven by another unmeasured biological character that is nested or overlapping with the observed character. In this case, the BiSSE model will usually be preferred over a "null" model of a single diversification rate regime (Rabosky and Goldberg 2015). To avoid a spurious conclusion of characterdependent diversification, an appropriate null model for comparison to a BiSSE model is a characterindependent model with two diversification rate categories (CID2), specified by an unobserved binary hidden character (Beaulieu and O'Meara 2016). Similarly, a full HiSSE model, where the observed character and a hidden character can each affect diversification rate (four diversification rate categories), should be compared to the character-independent model with four diversification rate categories (CID4).
The CID2 and HiSSE models include an unobserved character, with states A and B, that can additionally impact diversification. Model parameters include state-specific speciation and extinction rates for the four possible character-state combinations (0A, 1A, 0B, 1B) and transition rates between the four states. Here we used four transition rates (qAB, qBA, q01, and q10) that allow transition rates in a given character to be independent of the other character, except that simultaneous transitions in both characters are not possible. We considered three hidden state models. The first is a character-independent model where the hidden character, but not the observed character, affects diversification rate (CID2 model: λ0Α = λ1Α, λ0B = λ1B, μ0Α = μ1Α, μ0B = μ1B). We also fit a full HiSSE model where both the hidden and observed characters can affect diversification (model H1), with no constraints on parameter values. Finally, we fit a characterindependent model where the hidden character has four states (A-D), leading to eight character state combinations (0A, 1A, 0B, 1B, 0C, 1C, 0D, 1D). This CID4 model has 14 transition rates (qAB, qAC, qAD, qBA, qBC, qBD, qCA, qCB, qCD, qDA, qDB, qDC, q01, and q10) that are the natural expansion of the transition rates found in our CID2 and HiSSE models.
We generated posterior distributions of parameter values and sampled ancestral character states using MCMC implemented in RevBayes (Höhna et al. 2016) for the crown-ML and crown-ML-cropped datasets. We set lognormal prior distributions for speciation and extinction rates, where log-speciation and log-extinction were distributed with mean = 0 and standard deviation = 1. We set exponential prior distributions for transition rates with mean = 15 transitions / total treelength. We accounted for incomplete sampling by setting the probability of sampling species in the crown clade to 104/126. We note that we have sampled bee vs. hummingbird syndrome species in the crown clade at similar rates (83.5% of bee syndrome species and 78.2% of hummingbird syndrome species). Each MCMC chain drew 20,000 samples from the posterior distribution and we discarded the first 20% of samples as burn-in. We sampled ancestral character states during the MCMC in order to visualize estimates of the ancestral states by plotting posterior probabilities for each state at each node. We checked stationarity using Tracer to ensure ESS values for all parameters were above 200. The CID4 model did not reach stationarity after 20,000 generations, so we ran 40,000 generations for this model, at which point it reached stationarity.
We compared models by estimating the marginal likelihoods of competing models and comparing their relative fit to the data using Bayes factors. We estimated marginal likelihoods using the stepping-stone algorithm implemented in RevBayes, which estimates the marginal likelihood of a given model using a series of MCMC-like simulations that iteratively sample likelihood space between the prior and posterior distributions. We refer the reader to the RevBayes manual for additional information on the steppingstone algorithm (Höhna et al. 2017). We sampled from 50 steps between the prior and posterior probability distributions. For each step, we sampled 10,000 generations from the posterior distribution after discarding 10,000 as burnin. We calculated Bayes factor as twice the difference in marginal log likelihoods, using Bayes Factor = 1.16 as a significance threshold for model preference (Kass and Raftery 1995), as recommended in the RevBayes manual (Höhna et al. 2017).

FiSSE statistic calculations
We used the non-parametric Fast, intuitive State-dependent Speciation and Extinction (FiSSE) statistic (Rabosky and Goldberg 2017) to test for an association between pollination syndrome and diversification rate in the crown-ML and crown-ML-cropped datasets. FiSSE computes the average inverse equal splits measure for each character state (Λi), a metric that is correlated with, but not equal to, a tip-specific speciation rate. It then compares the observed difference in these values (Λ1 -Λ0) to a distribution of values simulated under a model where diversification rate is uncorrelated with the focal trait. We implemented FiSSE on the crown-ML and crown-ML-cropped datasets with the R FISSE.binary function (Rabosky and Goldberg 2017), using 10,000 simulations to test significance.

Additional information on metrics calculated for simulated trees
We calculated a set of metrics on our simulated trees to assess the degree of tippiness for our derived trait of hummingbird-adapted flowers. These metrics are well-described by Bromham et al. (2015), and we briefly describe them here. Tip age rank sum (TARS) represents the rank sum of terminal branches with the derived state and provides a relative measure of the ages of species with the derived state relative to those with the ancestral state. Number of tips per origin (NoTO) calculates the number of tips per origin, as determined through parsimony-based reconstruction, and provides a metric for the degree of tippiness of trait origins. Sum of sister-clade differences (SSCD) sums, for each node in the tree, the number of times the descendent sister clades differ in the modeled binary trait, providing a metric of phylogenetic clustering.

Simulated equilibrium proportions of bee vs. hummingbird syndrome species in Penstemon
Here we describe our methods for numerically calculating equilibrium proportions of bee vs. hummingbird syndrome species. In practice, this method produces effectively identical results as the equations presented in Appendix 2 of Maddison et al. (2007), yet allows us to calculate the time taken to reach equilibrium proportions.
The estimated parameters for BiSSE analyses are: r0 = diversification rate for bee syndrome species r1 = diversification rate for hummingbird syndrome species q01 = transition rate from bee to hummingbird syndrome q10 = transition rate from hummingbird to bee syndrome These parameters correspond to the rate of change over a unit time of 1.0, the total tree length. For simulations, we converted these rate parameters into approximations of instantaneous rates by considering time intervals of length 0.001. We performed the following transformations to obtain these instantaneous rate parameters: " * = 1 + " 1000 ⁄ ) * = 1 + ) 1000 ⁄ ") * = ") 1000 ⁄ )" * = )" 1000 ⁄ Since the diversification rates ri represent the net rate at which additional species are added (or subtracted), we must add 1 to the instantaneous rates to represent the current number of species.
With these quantities, two matrices are constructed. The first is the diversification matrix R: The species numbers after two time steps is found by: More generally, after t time units: The equilibrium proportions of the two species are then given by the eigenvector of the matrix: = corresponding to the largest eigenvalue of M, where the elements of the eigenvector are rescaled to sum to 1. We calculated these proportions for the posterior distribution of parameter values, after discarding the first 20% of values as burnin. From the distribution of proportions, we calculated the mean, 95% credible intervals, and 99% credible intervals.
We examined the time taken to approach equilibrium proportions by letting 16" = 8 1 0 9, starting with a single ancestral bee syndrome species, and calculating Nt until the proportion of hummingbird syndrome species reaches 90% and 95% of the expected equilibrium proportion. Again, this analysis was performed for all parameter combinations in the post-burnin posterior distributions. From these calculated distributions of times taken to reach 90% and 95% of equilibrium proportions, we calculated means and 95% credible intervals.

Features of aligned sequence datasets
Adding 45 new samples to our phylogenomic dataset (Wessinger et al. 2016) yielded concatenated alignments for the full dataset of 120 species and for the crown dataset of 104 species. The full dataset alignment included 2306 loci present in at least 20 species, with 72% missing data. The crown dataset included 2051 loci in at least 20 species, with 68% missing data.

Analyses on crown-SVD tree
Results from analyses on the crown-SVD tree are consistent with the analyses on the crown-ML presented in the main text. Table S5 reports model marginal likelihoods and confidence intervals for inferred parameters of BiSSE and HiSSE models. Figure S8 shows parameter estimates and ancestral character reconstructions for the full BiSSE model (B1) and the full HiSSE model (H1). The FiSSE statistic finds the inverse equal splits speciation rate for bee syndrome species is greater than that for hummingbird syndrome species (Λ0 = 2.48, Λ1 = 1.91), a difference that is significantly larger than expected if diversification rate is not associated with pollination syndrome (p = 0.015). Figure S9 shows calculated metrics on simulated trees. Figure S10 shows calculated equilibrium proportions of hummingbird syndrome species based on BiSSE model parameters. Tables   Table S1. Character data used in state-dependent diversification analyses on crown clade. Assignments of hummingbird specialist vs. bee syndrome/ despecialized ("partial") species follows (Wilson et al. 2007

Supplemental figures
Fig. S1. Species tree inferred for the full sample of 120 Penstemon species using SVDquartets. Hummingbird specialist species are marked with red dots. Nodes are colored according to bootstrap support values: black is greater than 95%, grey is 75-95%, white is less than 75%. Branch lengths were estimated using a constrained tree search in IQ-TREE under a GTR + Γ model. P . c y a n e u s P . s p e c io s u s P . p a y e t t e n s is P .p e n n e ll ia n u s P .e a to n ii P .l a e v is P .l e io p h y ll u s P .l o n g if lo ru s P .t id e s tr o m ii P .c y a n a n th u s P .le m hi en si s P. ba rb at us P. vir ga tus P.n eom exi can us P.co mar rhen us P.navajo ensis P.nudiflorus P.pseud oputus P.ca ryi P.p enl and ii P. de av er i P. pu tu s P .la br os us P .p a rv u s P .s tr ic tu s P .g la b e r P .h e n ri c k s o n ii P .c a r d in a li s P .f lo r id u s P .a n g u s ti fo li u s P .a r e n ic o la P .g ra n d if lo ru s P .m u rr a y a n u s P .n it id u s P .a c u m in a tu s P .c ya nt ho ph or us P. ha rri ng to ni i    . Character-independent patterns of diversification in the Penstemon crown clade. Panels A and C: visualization of model-averaged net diversification across the tree inferred using BAMM (red = higher diversification rate, blue = lower rate). Panels B and D: model-averaged net diversification rate inferred using BAMM as a function of time. Plots are shown for the crown-ML tree (panels A and B) and the crown-ML-cropped tree (panels C and D). The crown-ML and crown-ML-cropped trees each display a single diversification regime with no lineage-specific shifts in diversification rate, according to Bayes factor model comparison implemented in BAMMtools. In each of these trees, the net diversification rate slows down through time within the single regime. In panel B the slowdown is from 2.5105-2.400; in panel D the slowdown is from 3.1163-3.0193.