Floral evolution by simplification in Monanthotaxis (Annonaceae) and hypotheses for pollination system shifts

Simplification by reduction has occurred many times independently in the floral evolution of angiosperms. These reductions have often been attributed to changes in reproductive biology. In the angiosperm plant family Annonaceae, most species have flowers with six petals, and many stamens and carpels. In the genus Monanthotaxis several deviations from this pattern have been observed, including flowers that contain three petals and three stamens only. New DNA sequences were generated for 42 specimens of Monanthotaxis. Five chloroplast markers and two nuclear markers for 72 out of 94 species of Monanthotaxis were used to reconstruct a phylogeny of the genus, which revealed several well-supported, morphologically distinct clades. The evolution of four quantitative and two qualitative floral characters was mapped onto this phylogeny, demonstrating a reduction in flower size and number of flower parts in Monanthotaxis. A large variation in stamen forms and numbers, strong correlations between petal size, stamen and carpel number, combined with a non-gradual mode of evolution and the sympatric co-occurrence of Monanthotaxis species from different clades suggest that the high diversity in the African rainforest of this genus is caused by switches in pollination systems.

two perianth and stamen whorls 5 . Both increases in size and complexity, as well as simplification occurred later in evolution and can now be observed in extant angiosperms 12,13 . A reduction in the number of whorls of tepals and stamens has occurred in most modern angiosperms. The ancestral state of many other floral characters is uncertain for angiosperms, which is likely caused by the variability of many of those states in the early diverging angiosperms, complicating the inference of ancestral states 13 .
Among early diverging angiosperms, the Annonaceae are second in size after the Lauraceae, containing ca 2430 species 14 , and are amongst the dominant plant families in tropical forests worldwide 15 . The flowers of most species of Annonaceae contain three perianth whorls, viz. a single whorl of three sepals and two whorls of three petals. In the majority of species the flowers have a large number of stamens in multiple whorls, with the whorled pattern often becoming irregular with increasing numbers of stamens 16 .
Several deviations from the general floral pattern in Annonaceae were observed in the genus Monanthotaxis during an ongoing taxonomic revision [17][18][19] . This genus consists of lianas and is endemic to Africa and Madagascar. With 94 species it is one of the species-rich genera of Annonaceae and the second largest genus for Africa. In contrast to most other genera in the family, Monanthotaxis displays a huge variation in floral characters (Fig. 1). The flowers of most species are small relative to the range of sizes observed across Annonaceae. Species such as M. tripetala are amongst the smallest-flowered species in the family. Petals commonly occur in either one or two whorls. In some species, however, the petals form a single whorl at the floral base, whilst the outer petals overtop the inner petals at the apex in flower buds and give the appearance of two petal whorls in bud stage. In some species (e.g. M. bidaultii, M. tripetala), reduction of floral parts had led to a decrease in size, or even absence, of the inner whorl of petals. The number of stamens greatly varies from a few stamens in a single whorl (e.g. M. bidaultii, M. heterantha), up to 120 stamens in many whorls (e.g. M. gracilis). Several species possess an outer whorl of staminodes 18 . In most other genera of Annonaceae these characters vary hardly.
As most species have a moderate number of small flowers high up in the canopy of tropical rainforests, almost nothing is known about the ecology and pollination biology of Monanthotaxis. However, well-sampled phylogenies are a helpful tool to test hypotheses about evolution 20 , facilitated by observations on floral morphology made from herbarium specimens. The exact generic delimitation of Monanthotaxis and allied genera had long been unclear 21 . Recently, the generic boundaries of the genus have been clarified with a well-supported phylogeny containing about 40% of the species of Monanthotaxis. This resulted in the transfer of the genera Exellia, Gilbertiella and most of the African species of Friesodielsia to the genus Monanthotaxis, rendering Monanthotaxis morphologically well delimited and monophyletic 19 .
In this study, we increased the taxon sampling of Monanthotaxis to over 75% of the species diversity. Four quantitative and two qualitative characters, which are elements of floral reductions, were scored for these species. We traced changes in these characters over the phylogeny to produce and tested hypotheses about the evolution and diversification of the genus. More specifically, we examined the evolutionary mode of the floral characters, i.e. whether these have evolved following a gradual mode or pulsed mode of evolution, and the correlations between floral characters to test the possibility that small flowers in some Monanthotaxis species have evolved in concert with a reduction in number of parts. This is done to create testable hypotheses about the diversification of the genus.

Results
Phylogenetic analyses. The dataset used for the phylogenetic analyses comprised 6649 bp of sequence data that has been gathered for 88 specimens, including 80 specimens of Monanthotaxis representing 77% of the species of this genus. The best partitioning scheme found by Partitionfinder divided the nuclear and chloroplast markers in separate partitions and rendered the HKY + gamma substitution model best fitting in both cases. The analyses using CodeML did not find evidence for positive selection in the matK gene, while an indication for positive selection was found in 5 sites in the rbcL gene. Removing those five sites slightly improved the support values for a few nodes (Supplementary information Figs S1 and S2). Excluding the nine species for which less than half of the DNA sequences were available did not change the principle topology, but did improve the support values of seven nodes ( Fig. 2 and Supplementary information Figs S1 and S3). As the absence of these nine species would result in the exclusion of some of the relevant morphological character variation, all subsequent analyses were done using all species, but without the five codons of rbcL with possible positive selection.
The genus Monanthotaxis was highly supported to be monophyletic in both the Bayesian and maximum likelihood analyses, with Monanthotaxis obovata highly supported as sister to the rest of the genus. Most backbone nodes are not well supported, in contrast to more recent nodes, notably clades A, B, C, E and H ( Fig. 2 and Supplementary information Fig. S1).
Analyzing chloroplast and nuclear loci separately revealed a well-supported incongruent pattern within clade A (Supplementary information Figs S4 and S5). The nuclear data reveal a sistergroup relationship between M. ferruginea and M. bokoli, whereas the chloroplast data show that M. ferruginea is sister to M. orophila with strong support. All three species involved are closely related and differ only slightly in floral morphological characters. Therefore, we consider the effect of this incongruence on the ancestral state analyses negligible and performed all analyses after concatenation of the chloroplast and nuclear data.

Morphological characters and analyses. The results clearly show that simplification has taken place in
Monanthotaxis with strong correlations between number of stamens, petal size and number of carpels (Table 1). Furthermore, the outgroups and clades J and K of Monanthotaxis have many more stamens and larger petal sizes than the remaining Monanthotaxis species. The number of stamens ranged from 105 to 346 in the outgroups, and from 125 in Monanthotaxis gracilis to as few as three in M. bidaultii. Maximum petal size in Monanthotaxis ranged from 2.2 mm in M. tripetala to 50 mm in M. hirsuta, while reaching 130 mm in the outgroup species Dasymaschalon dasymaschalum (Fig. 3). Due to the low support for some of the backbone nodes, it is unclear ranged from 1 to 16 in Monanthotaxis (Fig. 4). Staminodes were present in species scattered across five different clades of Monanthotaxis and have arisen or disappeared multiple times. Unisexual flowers were confined to a single clade (Fig. 5).
The ancestor of Monanthotaxis most probably had bisexual flowers and lacked staminodes ( Fig. 5 and Supplementary Table S1). The ancestral state inferences of the continuous traits should be interpreted with caution, but suggest that the ancestor of Monanthotaxis had petals of 10 mm long, 31 stamens, 21 carpels and 3 ovules per carpel (Figs 3 and 4 and Supplementary Table S1). Strong to very strong correlations were found between number of stamen, petal size and the number of ovules per carpel, indicating that the simplification of flowers has occurred jointly for different parts of the flowers. The number of carpels is only strongly correlated with the number of stamens when considering Bayes factors, and less strong after the likelihood ratio test (in only 76% of the input trees). The correlation between petal size and number of carpels was strong considering Bayes factors and weak between number of ovules and number of carpels, while the likelihood ratio tests found no significant correlations between these variables in the majority of input trees. The discrete tests of correlation did not retrieve evidence for a correlation between presence of staminodes and flower sexuality ( Table 1).
The tests for mode of evolution significantly favoured a pulsed mode for the number of stamens and the number of carpels, indicated by kappa estimates of 0.01 and 0.1 respectively. No evidence for a significant deviation from a gradual mode of evolution was found for the number of ovules per carpel, while evidence for a deviation from a gradual mode of evolution was found using likelihood ratio tests for the petal length. The kappa estimate for petal length of 1.9 indicates that long branches contribute more to the evolution of petal length than would be expected under a gradual mode of evolution. However the Bayes factors did not indicate a deviation from a gradual mode of evolution for the petal length, suggesting that the signal is very weak if present ( Table 2).

Discussion
In this study the majority of the Monanthotaxis species were included in phylogenetic analyses and data of the floral characters were acquired for all species of which flowers are known. The results demonstrate that a reduction in flower size and number of floral parts has occurred in most species of Monanthotaxis. The methods also proved useful for the generation of new insights into the evolution of a group for which hardly any ecological information is available. In a previous study the monophyly and generic delimitation of Monanthotaxis were resolved 19 , based on a sample of 40% of the species of Monanthotaxis that still lacked some crucial species, such as those with unisexual flowers. In this study, the monophyly of Monanthotaxis is reaffirmed and although some nodes along the backbone of the tree receive weak support, most nodes within the genus are well supported (Fig. 2). Morphologically similar species are consistently closely related. As a result, the majority of the clades are easily recognizable by the combination of a few morphological characters. For example, the species of the well-supported clade A are distinct by the extra-axillary inflorescences, rounded flower buds, six petals in two whorls, and 30 or fewer stamens. The well-supported clade C is distinguishable by axillary inflorescences, ovate flower buds and petals in a single whorl at the base, with the outer petals overlapping the inner petals distally.  There is a trend across the phylogeny towards a reduction of stamen number and petal length in Monanthotaxis (Fig. 3). Indeed, all outgroup species have more than 100 stamens per flower (e.g. Fig. 1a), while the only clades of Monanthotaxis with more than 40 stamens per flower are the early diverging clades J and K. There are some other genera of Annonaceae which show only a few stamens, such as Orophea, Bocagea and Miliusa 22 . Within the tribe Uvarieae, Monanthotaxis is the only genus in which species with a low number of stamens occur. Consequently, the low number of stamens in most Monanthotaxis species is considered a reduction. Note, however, that the ancestral state for the genus Monanthotaxis was inferred as 31 stamen and petals of 10 mm, which is fewer than the number of stamens and a smaller petal size than the species in clades J and K, suggesting an increase has occurred in those clades. The species of these two clades were previously assigned to the genus Friesodielsia based on their superficial resemblance with that genus 19 .
The very strong correlation between number of stamens and petal length, and between number of stamens and number of carpels within the genus Monanthotaxis is consistent with the general notion that the number of floral structures depends on available space 23 . However, the reduction in flower size apparently did involve all floral parts of the species with unisexual flowers. Despite the consistent small petals, most of the unisexual flowers contain many carpels. Monanthotaxis wieringae, for example, has a maximum petal length of 5 mm, but has around 130 carpels. This discrepancy probably explains the weak correlation between carpel number and petal length in most of the likelihood ratio tests. However, the male flowers do follow the general reduction trend with fewer stamens in species with smaller flowers.
Pollen-ovule ratio is an indicator of the breeding system in angiosperms 24 . A high ratio indicates cross-pollination, while low pollen-ovule ratios generally are found in self-pollinating species. In view of this, the strong correlation between number of ovules per carpel and number of stamen, and between number of carpels and number of stamen, may suggest that no differences in breeding system exist within the genus Monanthotaxis. This is a tentative hypothesis, since the diversity in stamen forms (Fig. 1b,c,e) as well as the presence of species with unisexual flowers and bisexual flowers could indicate that different breeding systems occur in Monanthotaxis 25 . Moreover, it is questionable whether number of stamen can directly be linked to the number of pollen grains as there is considerable variation in anther cell size within the genus Monanthotaxis.
The presence of non-functional staminodes is often thought to be an intermediate stage in the reduction of stamens 10 . In the genus Monanthotaxis the staminodes do not have an apparent function; they form an outer whorl of reduced stamens which are often only visible as very tiny appendages less than half a millimeter in length, such as in M. zenkeri 18 . The fact that the presence of staminodes occurs on multiple independent branches on the phylogeny could indicate that the process of reduction is still ongoing in Monanthotaxis. Monanthotaxis whytei (Fig. 1d) is the only species on which an ontogenetic study has been performed 26 . First an outer whorl of 6 staminodes is formed; this is followed by an inner whorl of 9 staminodes then a further inner whorl of 9 stamens. The outer whorl of staminodes already stops developing early in the floral development and is not visible in the mature flowers, while the whorl of 9 staminodes is only visible as very small appendages (Fig. 1d). This could indicate that the reduction of stamens in the genus Monanthotaxis takes place centripetally, i.e. from the outer to the inner whorls and that the stamen whorl in species with only one whorl is homologous with the inner whorl in species with multiple whorls. It is interesting to note that the presence of staminodes in species with unisexual flowers (clade I, Fig. 1b,f,g) only occurs in male flowers and not in the female flowers a characteristic rarely reported in plants with unisexual flowers 10 .
Reduction trends can follow a gradual mode of evolution or the slightly controversial 27,28 pulsed mode of evolution, in which rapid change is followed by periods of stasis or little changes. Genome reductions and duplications, for example, by their nature have been demonstrated to follow a pulsed mode of evolution 29,30 , but this mode of evolution has rarely been established for morphological characters. The strong indication for a pulsed mode of evolution of the number of carpels and stamens may suggest that sudden events, such as changes in environment or pollination shifts, have triggered the floral reduction in Monanthotaxis. It must be noted that a limitation of the test used in this study is the simultaneous inference of the mode of evolution and of the amount of morphological change associated with speciation events 28 . Currently, methods are being developed which can disentangle these questions 31 . Alternatively, the inference of a pulsed mode of evolution for the number of stamens and carpels could be explained by the reduction of entire whorls. A gradual loss of whorls will be discernable as saltational changes in numbers of individual stamens and carpels, given the high number of floral parts per whorl. The irregular pattern of stamens and carpels in some species hampered the observation of the number of whorls in this study. Ontogenetic studies are needed for Monanthotaxis species with higher numbers of stamens and carpels to assess the exact number of whorls and assess the mode of evolution of those characteristics.   In general, we demonstrate an evolutionary reduction of flowers with a high number of floral parts to reduced flowers with as few as three petals and three stamens in the genus Monanthotaxis (Fig. 1b). The question remains whether a shift in pollination regime has triggered this reduction. In most species of Annonaceae the flowers are pollinated by beetles. Self-pollination is generally prevented by protogyny, viz. flowers entering the pistillate phase first, with a subsequent non-sexual phase and finally the staminate phase 32 . This is also the case in the genera Desmos, Dasymaschalon and Friesodielsia, the genera most closely related with Monanthotaxis for which pollination studies exist 32 . Protogyny also occurs in Monanthotaxis as in herbarium specimens only flowers in the pistillate phase or staminate phase have been observed. Flowers with overlapping staminate and pistillate functional phases have not been observed. The pollinators of Monanthotaxis however are unknown. Despite great differences in petal morphology in the tribe Uvarieae, most species are reported to be pollinated by beetles. However, pollination by cockroaches as well as bees has been reported in the genus Uvaria and flies and thrips have been reported as pollinators in some other genera outside the Uvarieae 32 . Small unisexual flowers, as those of most Monanthotaxis species in clade I, have been correlated with fly as well as wind pollination 25 . However, without further studies the question remains if the great diversity in stamen forms and number in Monanthotaxis also reflect a diversity in pollination syndromes in this genus.
Monanthotaxis is the second-most diverse genus of Annonaceae in Africa, especially in central Africa as many as 10 species can occur sympatrically. The majority of these sympatric species belong to different subclades, with most pairs of sister species showing allopatric distributions. Therefore, the diversification of Monanthotaxis may have been promoted by floral adaptations of different subclades. This enabled them to occupy different ecological microniches and consequently species of different clades are able to co-occur. Dispersal limitation could have promoted subsequent diversification. The pattern of allopatric and sympatric distributions can only be tested directly if more accurate distribution data become available. Recent expeditions have shown that the known distributions of species are larger than previously thought and these also uncovered many undescribed species 18 showing that the genus is undersampled in some areas of Africa. Finally, follow-up studies on the ontogeny, pollination and ecology are needed to understand more about the diversification of the genus Monanthotaxis and to answer questions on biodiversity of tropical rain forests.

Methods
Taxon sampling. For each species, subspecies and variety of Monanthotaxis at least one representative voucher specimen was selected and six species have been included twice in this study as these belong to allopatric populations with some morphological differences. Eight outgroup species were selected based on the phylogeny of Guo et al. 19 including all genera in the Dasymaschalon clade and representative genera of the tribe Uvarieae. 262 DNA sequences were newly generated for this study and 299 sequences were taken from previous studies 19,33 . All voucher specimens were identified by the first author in the framework of a taxonomic revision of the genus Monanthotaxis. Sequences for a total of 72 out of the 94 species (74 of 96 taxa when including varieties) of Monanthotaxis were successfully produced. This includes nine recently described new species 17,18 and six undescribed species (Hoekstra unpublished). An initial exploration of specimens from Madagascar indicated that there are at least 12 undescribed species from that region. Of the 22 taxa not included 8 species were only very recently recognized as new species to science. The other taxa are only known from very old collections and DNA extraction and amplification failed or no permission was obtained from the herbaria to sample those old collections.
DNA extraction, amplification and sequencing. DNA sequences of five chloroplast markers (matK, ndhF, psbA-trnH, rbcL and trnL-F) and two nuclear markers (ETS and ITS1-5.8S-ITS2) were generated for 42 specimens and downloaded from genbank for 46 specimens. DNA extraction, amplification and sequencing followed the protocol that has been previously described 34,35 and the PCR reaction for the nuclear markers as in Guo et al. 19 with the modification that 5 µl of 5x Betaine was added to each 25 µl PCR, because the ITS and ETS sequences have a high GC content (55-60%). The PCR products were sequenced on an ABI 3730 Sequencing platform by BaseClear (Leiden, The Netherlands). For voucher information and GenBank accession numbers see Supplementary Table S2.

Phylogenetic analyses. Trace files were checked and assembled in Sequencer 5.4 (Gene Codes
Corporation, MI U.S.A.) or BioEdit 7.2.5 36 . Sequences were aligned using the L-ins-I option in MAFFT v 7.245 37 and verified manually in Mesquite v 3.11 38 . Ambiguously aligned regions (in ITS and ETS) as well as microsatellites were excluded from the analyses 33 . Some species had an inversion of 14 positions in the psbA-trnH spacer. Following Pirie et al. 39 , this fragment was inverted when necessary, to retain any phylogenetic information. Gaps were coded following the method of Simmons and Ochoterena 40 using FastGap 1.2 41 .
Phylogenetic trees were inferred using maximum parsimony analyses, maximum likelihood analyses and Bayesian Inference. With the program Partitionfinder v. 1.1.1 42 we inferred the best partitioning scheme as well as the best substitution models for those partitions. Individual loci, and separate codon positions for the coding genes, were defined as the data blocks. All possible partitioning schemes and substitution models were fitted on the dataset using the greedy algorithm of PhyML in Partitionfinder and the best partitioning scheme was selected based on the Bayesian Information Criterion.
Maximum parsimony analyses were performed in PAUP* version 4.0a151 43 with the heuristic search option and tree bisection-reconnection (TBR) branch swapping, with 1000 random addition sequence replicates, saving 50 trees per replicate. Character states were treated as unordered and of equal weight 44 . To assess clade support nonparametric bootstrap analyses were performed with 1000 bootstrap replicates and 100 random addition sequence replicates per bootstrap replicate, saving 50 trees per replicate. Bayesian inference was conducted with MrBayes v3.2.6 45 on the CIPRES gateway server 46 . The best fitting model of sequence evolution as defined by Partitionfinder was applied to the DNA sequence data and the restriction site model with the setting "coding = variable" was applied to the gap-coding data. It has been shown that MCMC analyses of closely related species may get biased towards excessively long branch-length estimates 47 and therefore following Guo et al. 19 , the temperature parameter was set to 0.08 and the mean branch length prior was set to 0.01 (brlenspr = unconstrained: exponential (100.0)). These settings improved mixing of the chains compared to the default values of the program. Four different chains were run for 10 million generations sampling every 1000 th generation. Convergence was assessed using the R-package RWTY version 1.0.1 48 as this package includes an estimation of the effective sample size (ESS) for the tree topology parameter 49 as well as ESS values for other parameters and diagnostic plots.
Maximum likelihood analyses were run on the CIPRES Gateway server with RAxML version 8 50 . The GTR-gamma model of substitutions was used for the DNA sequence data and the bingamma model for the gap-coded data. To assess clade support a rapid bootstrap analysis was performed on the best-scoring tree with 1000 nonparametric bootstrap iterations.
In many angiosperms positive selection takes place on the active sites of the rbcL-gene and removing those sites could improve phylogenetic resolution 51 . This has been found in a family-wide study of Annonaceae 52 . We tested for signs of positive selection in the genes rbcL and matK, using the branch-site model of positive selection 53 in the program Codeml of the PAML4.8 package 54 . A simplified version without branch lengths of the most likely tree from the maximum likelihood analyses was used as a backbone phylogeny. We assessed the presence of positive selection for the branch subtending the Monanthotaxis crown node, for each gene separately. These models were compared using a likelihood ratio test against a null model in which the value ω was fixed. When the likelihood ratio test significantly demonstrated the presence of positive selection (P < 0.05), the Bayes empirical Bayes (BEB) 53 was used to calculate the posterior probabilities for site classes and to identify the sites under positive selection. These codons were deleted from the alignment and subsequently all phylogenetic inferences were redone with the new dataset without codons under positive selection.
Bayesian, maximum likelihood and maximum parsimony analyses were rerun separately on the nuclear and chloroplast loci to check for incongruences. Nine specimens, for which more than half of the DNA sequences were missing (3000 bp), were excluded from the dataset and all phylogenetic inferences were rerun to test if support values for any nodes improved. The support values from all different runs were summarized and compared using a customized Python script with the library Dendropy version 4.1.0 55 . Maximum parsimony and maximum likelihood bootstrap values from 85 to 100 were considered as strong node support, values from 75 to 84 as moderate support and values from 50 to 74 as weak support. Posterior probabilities higher than 0.95 were considered as supported, below 0.95 as unsupported.
Morphological character sampling. Six morphological characters were selected to be used as a proxy for floral reduction. Four characters are quantitative, i.e. maximum petal length, number of stamens, number of carpels and number of ovules per carpel, and two characters are qualitative, i.e. flower sexuality and presence of staminodes. For species of Monanthotaxis all these characters were (re-)measured from herbarium material as in published literature sometimes errors were encountered. For the outgroup species data was taken from literature [56][57][58] and missing data (including number of stamen for all species) was measured from herbarium material.
Quantitative variables were transformed to logarithmic scale. However, variation occurs in the number of stamens and carpels within species and even within specimens, with variation increasing with the number of stamens and carpels per species. Species with generally nine stamens often may have a few flowers with eight stamens and sometimes even twelve stamens, while in a species with around 100 stamens the number can vary from 80 to 120 stamens. To take the increasing absolute value of variation with increasing number of stamens and or carpels into account, the average of the logarithm of the minimum and maximum number of stamens and carpels per species was used in the analyses.
The logarithm of the average of the maximum and minimum number of ovules per carpel was taken for species with a variable number of ovules per carpels. Species with a single ovule per carpel rarely have carpels with two ovules, and vice versa; these species were scored as having the predominantly occurring number of ovules.
The morphological characters could be observed for almost all species, with the exception of Monanthotaxis sterilis for which neither flowers nor fruits are known. Other species for which some of the characters could not be Character mapping and analyses. Morphological characters were optimized and analyzed on phylogenetic trees using BayesTraits version 2.0 59 . One thousand randomly chosen trees from the Bayesian analyses, after discarding the burnin and outgroup species, were selected to account for phylogenetic uncertainty.
Throughout the analyses in BayesTraits, model testing was used to test whether more complex models better fitted the data than simpler models. Three different model tests were performed: likelihood ratio tests were used for Maximum Likelihood (ML) analyses for nested models, AIC relative likelihoods were used for Maximum Likelihood (ML) analyses of non-nested models and Bayes factors for Markov Chain Monte Carlo (MCMC) methods. The likelihood ratio tests were performed for each of the one thousand input trees based on the complex and simpler models and for each tree was tested if the p-value of the likelihood ratio statistic was below 0.05. For the AIC-test, the AIC values were calculated for each of the one thousand input trees of the complex and simpler model. Next, the relative likelihood of the simpler model was calculated and values less than 0.05 were considered evidence that the complex model better fitted the data than the simpler model. For the MCMC analyses Bayes factors were calculated from the marginal likelihoods using a stepping stone sampler 60 , as marginal likelihoods estimated by a stepping stone sampler have been shown to be more robust than the harmonic mean marginal likelihoods 61 . One hundred stones were used to calculate an estimate of the marginal likelihood and each stone was run for 10,000 iterations. Bayes factors higher than 2 were considered positive evidence, higher than 5 strong evidence and higher than 10 as very strong evidence for selection of the more complex model compared to the simpler model with less parameters.
To infer the ancestral states of the discrete characters, a continuous-time Markov model using MCMC methods 62 was applied. Hyper-priors were used to seed the mean of the exponential priors from a uniform distribution ranging from 0 to 100. Model tests using Bayes factors were applied to test whether different transition rates from one state to another state better fitted the data than both rates set equal. For the character of uni-or bisexual flowers, the Bayes factors did not reject the hypothesis that the transition rates are equal, this was accommodated in the prior settings. While equal rates were significantly rejected (P < 0.05) for the ancestral states of the staminodes presence and those were inferred using reversible-jump Markov Chain Monte Carlo 63 .
A two-step process was run to infer the ancestral states of the continuous characters. First an MCMC with a continuous random walk model 64 was run to estimate a distribution of models from the available data and subsequently these models were used to infer the ancestral states. Bayes factors of initial analyses indicated that the trait 'number of stamens' does not evolve along the phylogeny following a Brownian motion model of evolution. Therefore, the lambda parameter was set to be estimated in analyses with the number of stamens character. Finally, after discarding the burnin, the mean and standard deviation of the ancestral state inferences were calculated for each node.
For both the discrete and continuous characters five runs using MCMC methods were performed for each trait to infer the ancestral nodes. The runs consisted of 10 million generations with a burnin of 100,000 generations and sampling every 10,000 generations. Convergence was checked by calculating the ESS-values with the R-package coda version 0.19-1 and by verifying if the acceptance rate for proposed changes to the chain lies between 20 and 40%.
To test for correlated evolution between 2 character traits, both 5 MCMC and 2 ML analyses were run for each model. The binary traits were first run with the discrete independent model in BayesTraits in which the traits are assumed to have evolved independently. Then the traits were run with the discrete dependent model which assumes that the traits are correlated. For the continuous characters they were first run with the continuous random walk model and subsequently those analyses were rerun, but with the correlation between the two traits being forced to be zero. Both Bayes Factors and likelihood ratio tests were performed to test if the models that assume correlation between the traits better fitted the data than the models without such correlation.
We tested whether the four continuous morphological characters evolved along the phylogeny following a pulsed or a gradual mode of trait evolution. In the same analyses we also tested if morphological changes concentrate around speciation events. First, analyses were run with the kappa scaling parameter set to be estimated and subsequently run with the kappa parameter fixed to 0 and additionally run with the kappa parameter fixed to 1. A kappa close to 0 indicates a pulsed mode of evolution for the involved trait, while a kappa of 1 indicates a gradual mode of evolution. Values lower than 1 indicate that shorter branches contribute more to the character evolution, while values higher than 1 indicate that longer branches contribute more. Model testing as described for the correlation tests were used to test if the kappa parameter significantly differed from 1 and/or from 0.