Phenotypic plasticity guides Moricandia arvensis divergence and convergence across the Brassicaceae floral morphospace

Phenotypic plasticity guides Moricandia arvensis divergence and convergence across the Brassicaceae floral morphospace JoseM.Gomez 1,2 , Adela Gonzalez-Megıas 2,3 , Eduardo Narbona 4 , Luis Navarro 5 , Francisco Perfectti 2,6 and Cristina Armas 1 1 Estacion Experimental de Zonas  Aridas (EEZA-CSIC), E-04120 Almerıa, Spain; 2 Research Unit Modeling Nature, Universidad de Granada, E-18071 Granada, Spain; 3 Departamento de Zoologıa, Universidad de Granada, E-18071 Granada, Spain; 4 Departamento de Biologıa Molecular e Ingenierıa Bioquımica, Universidad Pablo de Olavide, E-41013 Sevilla, Spain; 5 Departamento de Biologıa Vegetal y Ciencias del Suelo, Universidad de Vigo, E-36310 Vigo, Spain; 6 Departamento de Genetica, Universidad de Granada, E-18071 Granada, Spain Authors for correspondence: Jos eM.G  omez Email: [email protected] Adela Gonz alez-Meg ıas Email: [email protected] Received: 19 July 2021 Accepted: 10 October 2021 New Phytologist (2021) doi: 10.1111/nph.17807 Key words: Brassicaceae, floral morphospace, Moricandia, phenotypic convergence, phenotypic divergence, phenotypic plasticity, pollination niches. Summary Many flowers exhibit phenotypic plasticity. By inducing the production of several phenotypes, plasticity may favour the rapid exploration of different regions of the floral morphospace. We investigated how plasticity drives Moricandia arvensis, a species displaying withinindividual floral polyphenism, across the floral morphospace of the entire Brassicaceae family. We compiled the multidimensional floral phenotype, the phylogenetic relationships, and the pollination niche of over 3000 species to construct a family-wide floral morphospace. We assessed the disparity between the two M. arvensis floral morphs (as the distance between the phenotypic spaces occupied by each morph) and compared it with the family-wide disparity. We measured floral divergence by comparing disparity with the most common ancestor, and estimated the convergence of each floral morph with other species belonging to the same pollination niches. Moricandia arvensis exhibits a plasticity-mediated floral disparity greater than that found between species, genera and tribes. The novel phenotype of M. arvensis moves outside the region occupied by its ancestors and relatives, crosses into a new region where it encounters a different pollination niche, and converges with distant Brassicaceae lineages. Our study suggests that phenotypic plasticity favours floral divergence and rapid appearance of convergent flowers, a process which facilitates the evolution of generalist pollination systems. Introduction The expansion of lineages across the phenotypic space is a complex process that may result from the ecological opportunities that emerge when said lineages are exposed to new ecological niches (Schluter, 2000; Nosil, 2012). When this occurs, divergent selection on some phenotypes results in phenotypic diversification among lineages, boosting morphological disparity (the phenotypic space occupied by the lineages), triggering a morphological radiation, and eventually filling the phenotypic space (Sim~oes et al., 2016; Walden et al., 2020). Because morphological spaces fill up as lineages diversify (Pie & Weitz, 2005), unoccupied regions become rare in highly diversified lineages (Winemiller et al., 2015). Under these circumstances, entering into a new region usually entails sharing it with other species exploiting the same ecological niche (Losos, 2011; Winemiller et al., 2015; Stayton, 2020). In this situation, independent lineages tend to evolve similar phenotypes through convergent evolution (Losos, 2011; Pearce, 2011). In diversified lineages occupying filled morphospaces, divergent and convergent evolution are ineludibly connected (Winemiller et al., 2015; Pigot et al., 2020), and both processes contribute significantly to shaping the geometry of the morphospace occupation (Pearce, 2011; Stayton, 2020). The evolution of flowers in angiosperms is largely driven by interaction with pollinators (Lloyd & Barrett, 1996; Fenster et al., 2004; Endress, 2011; Van der Niet & Johnson, 2012; Gomez et al., 2014; Peter & Johnson, 2014; van der Niet et al., 2014). Shifts between pollination niches explain floral divergence and are associated with diversification of flowers in multiple groups of plants (Whittall & Hodges, 2007; Smith et al., 2008; Lagomarsino et al., 2016; Serrano-Serrano et al., 2017; Kriebel et al., 2020; Dellinger et al., 2021). On the other hand, floral convergence has also been tied to pollinator shifts in some plant groups (Smith & Kriebel, 2018). In fact, selection exerted by functionally equivalent pollinators is considered to be the main driver of the convergent evolution of floral traits in different plant species, the so-called pollination syndromes (Faegri & van der Pijl, 1980; Dellinger, 2020; Phillips et al., 2020; Wessinger & Hileman, 2020). All of this suggests that pollinators are an important motor shaping the pattern and geometry of occupation of floral phenotypic space (Ollerton & Watts, 2000; Ollerton Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation New Phytologist (2021) 1 www.newphytologist.com This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. Research et al., 2009; Chartier et al., 2014, 2017; Smith & Kriebel, 2018; Hervıas-Parejo et al., 2019). Phenotypic plasticity might elicit the emergence of novel phenotypes with new adaptive possibilities, which may be beneficial in some contexts (West-Eberhard, 2003; Sultan, 2015). It is increasingly acknowledged that many floral traits exhibit phenotypic plasticity in response to diverse environmental stimuli and to cope with many biotic and abiotic stressors (Harder & Johnson, 2005; Campbell et al., 2019; Rusman et al., 2019a,b). By inducing the production of contrasting phenotypes, plasticity may favour the exploration of different regions of the floral phenotypic space almost simultaneously (Gomez et al., 2020). This creates an opportunity for plastic species to display derived floral morphs (phenotypes) that may diverge from their linages/ancestors and converge with species already located in other regions of the floral morphospace. We explored this idea by investigating how plasticity drives the Purple Mistress, Moricandia arvensis, across the morphological space of the entire Brassicaceae family. This mustard species exhibits extreme polyphenism in flowers (Gomez et al., 2020). Individuals of M.arvensis display two flower morphs in response to seasonal changes in temperature, radiation, and water availability: large, cross-shaped lilac flowers in spring but small, rounded, white flowers in summer (Fig. 1; Gomez et al., 2020). By studying the multidimensional floral phenotypes, the phylogenetic relationships, and the pollination niches of over 3000 Brassicaceae species, we demonstrate that phenotypic plasticity causes this mustard species to produce a phenotype that diverges from its ancestors and close relatives, crosses into a new region of the ecological space, and converges morphologically with distant Brassicaceae lineages belonging to the same pollination niche. This finding may have striking implications for understanding how pollination systems may evolve in generalist species. Materials and Methods Floral traits Brassicaceae is one of the largest angiosperm families, with almost 4000 species grouped in 351 genera and 51 tribes (Kiefer et al., 2014; Koch et al., 2018; Walden et al., 2020). We recorded from the literature the following 31 floral traits in 3140 Brassicaceae plant species belonging to 330 genera and 51 tribes: plant height; flower display size; inflorescence architecture; presence of apetalous flowers; number of symmetry axes of the corolla; orientation of dominant symmetry axis of the corolla; corolla with overlapped petals; corolla with multilobed petals; corolla with visible sepals; petal length; sepal length; asymmetric petals; petal limb length; length of long stamens; length of short stamens; stamen dimorphism; tetradynamous condition; visible anthers; exserted stamens; number of stamens; concealed nectaries; petal carotenoids; petal anthocyanins; presence of bullseyes; presence of veins in the petals; coloured sepals; relative attractiveness of petals vs sepals; petal hue; petal colour; sepal hue; sepal colour (see Supporting Information Notes S1; Dataset S1 for details). All of these traits have been proven to be important for interactions with pollinators (Notes S1). Family-wide floral morphospace Using the original multidimensional trait-species matrix, we built a floral morphospace. For this, we reduced the highdimensional matrix of floral traits to a two-dimensional space using an ordination technique (Legendre & Legendre, 2012). Because the floral traits included in this study were quantitative, semi-quantitative and qualitative, we used ordination techniques based on dissimilarity values. We first constructed a pairwise square distance matrix of length equal to the number of species included in the analysis (n= 3140). We used the Gower distance, the number of mismatched traits over the number of shared traits. This dissimilarity index is preferable to the raw Euclidean distance when discrete and continuous traits co-occur in the same dataset (Legendre & Legendre, 2012; Guillerme & Cooper, 2018). To ensure an accurate description of the distribution of the species in the morphospace, we first run a principal coordinate analysis (PCoA), a technique providing a lowdimensional Euclidean representation of a set of objects whose relationship is measured by any dissimilarity index. We corrected for negative eigenvalues using the Cailliez procedure (Legendre & Legendre, 2012). Afterwards, we used this metric 10 cm25 cm (b) (a) Fig. 1 Moricandia arvensis life habit. Detailed view of plants and flowers of M. arvensis in spring (a) and summer (b). New Phytologist (2021) www.newphytologist.com Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation Research New Phytologist 2 configuration as the initial configuration to run a nonmetric multidimensional scaling (NMDS) algorithm (Legendre & Legendre, 2012), a method that further optimises the sample distribution such that more variation in species composition is represented by fewer ordination axes. Unlike methods that attempt to maximise the variance or correspondence between objects in an ordination, NMDS attempts to represent, as closely as possible, the pairwise dissimilarity between objects in a low-dimensional space (Legendre & Legendre, 2012). Objects that are ordinated closer to one another are likely to be more similar than those further apart (Legendre et al., 2005). This method is more robust than distance-based methods when the original matrix includes variables of contrasting nature. Nonmetric multidimensional scaling is an iterative algorithm that can fail to find the optimal solution. We decreased the potential effect of falling in local optima by running the analysis with 5000 random starts and iterating each run 1 910 6 times (Mair et al., 2016). The NMDS used a monotone regression minimizing Kruskal’s stress-1 (Kruskal, 1964a,b) and compared each solution using Procrustes analysis, retaining that with the lowest residual. Because many species did not share trait states, a condition complicating ordination, we used stepacross dissimilarities, a function that replaces dissimilarities with the shortest paths stepping across intermediate sites, while regarding dissimilarities above a threshold as missing data (De’ath, 1999). We used weak tie treatment, allowing equal observed dissimilarities to have different fitted values. The scores of the species in the final ordination configuration were obtained using weighted averaging. We checked whether the reduction in dimensionality maintained the between-species relationship by checking the stress of the resulting ordination and calculating the goodness of fit for points in nonmetric multidimensional scaling (Mair et al., 2016). All ordinations were done using the R packages VEGAN (Oksanen et al., 2013) and ECODIST (Goslee & Urban, 2007). It is important to note that, although the transfer function from observed dissimilarities to ordination distances is nonmetric, the resulting NMDS configuration is Euclidean and rotationinvariant (Oksanen, 2020). Morphological disparity The morphological disparities were calculated using indices related to the distance between elements (Guillerme et al., 2020a, b). We first determined the disparity between the spring and summer phenotype of M. arvensis as their Euclidean distance in the floral morphospace (Guillerme et al., 2020b). Afterwards, we calculated the pairwise disparities between all species included in our study using the same approach. To find out how intense the plasticity-mediated disparity in M. arvensis is, we compared its value with the disparity values observed at five different taxonomic levels in our dataset: within species (Dataset S2), between Moricandia species, between species of the same genera, between species of different genera and between species of different tribes (Dataset S1). Disparity values were calculated using the function dispRity of the R package DISPRITY using the command centroid (Guillerme, 2018). The statistical differences between the M. arvensis plasticity-mediated disparity and the disparity of each taxonomic level were tested using Z-score tests. Morphological divergence of the plastic phenotypes Divergence in floral phenotype was estimated by calculating the disparity of the two floral morphs of M. arvensis to its ancestors. We determined the floral phenotype of the most recent common ancestor (MRCA) by projecting a recent time-calibrated phylogeny of the genus Moricandia (Perfectti et al., 2017) onto the morphospace. This phylogenetic analysis included all known species belonging to the genus Moricandia (eight species), as well as the sister genus Rytidocarpus and the close genus Eruca, totalling 15 species, and was made using a nuclear (the internal transcribed spacers of the ribosomal DNA) and two plastidial regions (parts of the NADH dehydrogenase subunit F gene and the trnT-trnF region). Because disparity analyses are sensitive to the tree topology and the inferred branch lengths (Guillerme et al., 2020a), we also calculated the divergence using three recent time-scaled phylogenies that included information on M. arvensis and other Moricandia (Gaynor et al., 2018; Smith & Brown, 2018; Huang et al., 2020). The taxonomy of the species included in each tree was checked and updated using the species checklist, with accepted names provided by Brassibase (https://brassibase. cos.uni-heidelberg.de/) (Koch et al., 2012, 2018; Kiefer et al., 2014; Walden et al., 2020). Once we inferred the coordinates of the MRCA in the morphospace according to each of the four phylogenies, we calculated the disparity of all Moricandia species and the two plastic phenotypes of M. arvensis to the MRCA. Afterwards, we calculated the divergence of the two plastic phenotypes from the direct ancestor of M. arvensis. In addition, we calculated the divergence from the direct ancestors of the rest of the Brassicaceae species included in these four phylogenies. All floral divergences were calculated using the command ancestral.dist of the function dispRity in the R package DISPRITY (Guillerme, 2018). Morphospatial variation in pollination niches We compiled a massive database of 21 212 records, comprising 456 031 visits by over 800 animal species from 19 taxonomical orders and 276 families to 554 Brassicaceae species belonging to 39 tribes (Dataset S3). Information for the database came from the literature, personal observations, online repositories and personal communications of several colleagues (Dataset S3; Table S1). In those species studied by us (coded as UNIGEN data origin in Dataset S3), we conducted flower visitor counts in 1–16 populations per plant species. We visited the populations during the blooming peak, always at the same phenological stage and between 11:00 h and 17:00 h. During these visits, we recorded the insects visiting the flowers for 2 h without differentiating between individual plants. Insects were identified in the field, and some specimens were captured for further identification in the laboratory. We only recorded those insects contacting anthers or stigma and making legitimate visits for at least part of their foraging on flowers. We did not record petal-eating or Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation New Phytologist (2021) www.newphytologist.com New Phytologist Research 3 nectar-thieving insects which did not make a legitimate visit. The information obtained from the literature and online repositories (coded as LITERATURE data origin in the Dataset S3) includes records done during ecological studies, taxonomical studies and naturalistic studies. The reference of every record is included in the dataset. We thereby grouped all pollinators visiting the Brassicaceae species into 43 functional groups using criteria of similarity in body length, proboscis length, morphological match with the flower, foraging behaviour, and feeding habits (Fenster et al., 2004; Gomez et al., 2015a, 2016; Table S2). We tested the autocorrelation across the morphospace in the abundance of each of these functional groups using a multivariate Mantel test. The Pearson correlation method was used, and the statistical significance was found after performing 999 bootstrap replicates (Legendre & Legendre, 2012). The test was performed in the R package VEGAN (Oksanen et al., 2013). We determined the occurrence of different pollination niches in the set of plant species included in this study using bipartite modularity, a complex-network metric. Modularity has proven to be a good proxy of interaction niches both in ecological networks, those including coexisting species or populations, as well as in clade-oriented networks, those including species with information coming from disparate and contrasting sources (Gomez et al., 2010). We constructed a weighted bipartite network, including pollinator data of four M. arvensis populations during the spring and summer flowering. In this network, we pooled the data from the different individuals in a population and did not consider the time difference involved in sampling across different species. We removed all plant species with <20 visits. We subsequently determined the modularity level in this weighted bipartite network using the QUANBIMOalgorithm (Dormann & Strauss, 2014). This method uses a simulated annealing Monte-Carlo approach to find the best division of populations into modules. A maximum of 10 10 Markov chain Monte Carlo (MCMC) steps with a tolerance level =10 –10 was used in 100 iterations, retaining the iterations with the highest likelihood value as the optimal modular configuration. We tested whether our network was significantly more modular than random networks by running the same algorithm in 100 random networks, with the same linkage density as the empirical one (Guimera & Amaral, 2005). Modularity significance was tested for each iteration by comparing the empirical vs the random modularity indices using a Z-score test (Dormann & Strauss, 2014). After testing the modularity of our network, we determined the number of modules (Newman, 2004). We subsequently identified the pollinator functional groups defining each module and the plant species ascribed to each module. Modularity analysis was performed using the R package BIPARTITE v.2.0 (Dormann et al., 2008). We quantified the niche overlap between all pairs of Brassicaceae species using the Czekanowski index of resource utilization, an index that measures the area of intersection of the resource utilization histograms of each species pair (Feinsinger et al., 1981). This index was calculated using the function niche.overlap in the R package SPAA (Zhang & Zhang, 2013). Estimation of ancestral values of pollination niches The ancestral states of the pollination niche were inferred for the Moricandia lineage by simulated stochastic character mapping of discrete traits with Bayesian posterior probability distribution (Huelsenbeck et al., 2003; Bollback, 2006). Three models of character evolution (ER, equal Rates; SYM, symmetric; and ARD, all rates different) were evaluated using the fitDiscrete function of the R package GEIGER (Harmon et al., 2008). The best model was selected using the Akaike Information Criterion (AIC), and was used for stochastic character mapping. The posterior distribution of the transition rate matrix was determined using an MCMC simulation, and the stochastic mapping was simulated 100 times. Stochastic character mapping was performed using the make.simmap function and a plot of posterior probabilities was mapped using the describe.simmap function in R package PHYTOOLS (Revell, 2012). This analysis was performed using the Perfectti et al. (2017) phylogeny, the most accurate phylogeny to date for Moricandia. Morphological convergence We used three different approaches to detect morphological convergence, two based on the comparison of phenotypic and phylogenetic distances (Arbuckle et al., 2014; Stayton, 2015) and the other based on comparing the angles formed by two tested clades from their most recent common ancestor with the expected angle according to null evolutionary models (Castiglione et al., 2019). Because all of these analyses are sensitive to the number of tips in the phylogeny and the inferred branch lengths, we tested for the occurrence of morphological convergence using the three independent, time-calibrated phylogenies described above (Gaynor et al., 2018; Smith & Brown, 2018; Huang et al., 2020). Under the first approach, four distance-based measures of convergence (C1–C4) were calculated between two lineages relative to their distance at the point in evolutionary history where the two lineages were maximally dissimilar (Stayton, 2015). C1 specifically measures the proportion of phenotypic distance closed by evolution, ranging from 0 to 1 (where 1 indicates complete convergence). To calculate C1, ancestral states are reconstructed (via a Brownian motion model of evolution) for two or more putatively convergent lineages, back to their most recent common ancestor. The maximum phenotypic distance between any pair of ancestors (Dmax) is calculated, and compared with the phenotypic distance between the current putatively convergent taxa (Dtip). The greater the difference between Dmax and Dtip, the higher the index. C2 is the raw value of the difference between the maximum and extant distance between the two lineages. C3 is C2 scaled by the total evolution (sum of squared ancestor-to-descendant changes) between the two lineages. C4 is C2 scaled by the total evolution in the whole clade. These four measures quantify incomplete convergence in multidimensional space. The significance of these four metrics was found by running 1000 simulations for each comparison using Brownian motion on a variance–covariance matrix based on data-derived parameters, with convergence measures for each simulation New Phytologist (2021) www.newphytologist.com Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation Research New Phytologist 4 calculated to determine whether the observed Cvalue is greater than that expected by chance. A priori focal groups forming the basis of convergence tests were the two M. arvensis niches. These analyses were performed using the R package CONVEVOL (Stayton, 2018). The second approach to measuring convergence used the Wheatsheaf metric (Arbuckle et al., 2014). This index generates phenotypic (Euclidean) distances from any number of traits across species and penalizes them by phylogenetic distance before investigating similarity (in order to weight close phenotypic similarity higher for distantly related species). It uses an a priori designation of convergent species, which are defined as species belonging to a niche for which the traits are hypothesized to converge. The method then calculates a ratio of the mean (penalized) distances between all species to the mean (penalized) distances between allegedly convergent species. The index detects whether convergent species diverge more in phenotypic space from the nonconvergent species and show a tighter clustering to each other (Arbuckle et al., 2014). The significance of this index was found by comparing the empirical values of the index with a distribution of simulated indices obtained by running 5000 bootstrap simulations. These analyses were performed using the R package WINDEX (Arbuckle & Minter, 2015). The third approach to measuring convergence was based on comparing the angles formed by two tested clades from their most recent common ancestor with the expected angle according to null evolutionary models (Castiglione et al., 2019). Under the ‘state case’, search.conv computes the mean angle over all possible combinations of species pairs using one species per state. Each individual angle is divided by the patristic distance (the sum of the lengths of the branches that link two nodes in a phylogenetic tree) between the species. Significance is assessed by contrasting this value with a family of 1000 random angles obtained by shuffling the state across the species (Castiglione et al., 2019). These analyses were performed using the R package RRPHYLO (Raia et al., 2019). Once we tested the occurrence of convergence in the M. arvensis pollination niches, we assessed whether plasticity caused the evolution of morphological convergence in any of the two morphs of this species. We first assessed the convergence region of the Moricandia lineage, the region that includes the species converging morphologically to the Moricandia lineage. To do this, we found the clades that had significantly smaller angles than the expected angle according to null evolutionary models. These analyses were performed using the R package RRPHYLO (Raia et al., 2019). Afterwards, we checked whether any of the two M. arvensis floral phenotypes entered the region of the phylomorphospace defined by their pollination niches. We used the C5 index, a frequency-based measure that quantifies and reports the number of convergent events where lineages evolve into a specific region of morphospace (crossing it from outside) (Stayton, 2015, 2018). C5 sums the number of times through the evolution of a clade that lineages cross into a convergent region of a phylomorphospace (the phylogenetic connections between taxa represented graphically in a plot of morphological space). This analysis was performed using the RpackageCONVEVOL (Stayton, 2018). Results Floral morphospace of Brassicaceae The resulting morphospace (Fig. 2a) was significantly correlated with the initial PCoA configuration (r= 0.40, P< 0.0001; Mantel test) and was a good representation of the original relationship among the species (R 2 >0.95, Stress =0.2; Fig. 2b). The distribution of the species across the morphospace was significantly associated with different pollination traits (Fig. S1; Table S3). Species in the central region were mostly medium-sized plants bearing a moderate to high number of small, polysymmetric white flowers with short corolla tubes, exposed nectaries and visible sepals (Figs a, S1). Species in the bottom right corner were small or prostrate, bearing minute flowers, often apetalous and with just two or four stamens, whereas species located in the bottom left corner were medium-sized plants with asymmetric flowers arranged in corymbous inflorescences. Plants with yellow flowers were located in the right region of the morphospace. By contrast, large plants with strongly tetradynamous androecium and large, veined, dissymmetrical to asymmetrical, pink to blue flowers with concealed nectaries, long corolla tubes and bullseyes were located in the upper left region (Figs 2a, S1). Plasticity-mediated floral disparity Moricandia arvensis, when blooming in spring (Fig. 1a), occupies this later peripheral region of the morphospace, close to other Moricandia species (black dots in Fig. 2a). However, during summertime, individuals of M. arvensis are shorter and produce fewer, much smaller flowers with white, unveined and rounded corollas with overlapped petals and green sepals that are mostly arranged alone without forming inflorescences (Fig. 1b). Due to this radical phenotypic change, the summer phenotype of M. arvensis was located in a different, more central position of the floral morphospace (Fig. 2a), far away from the region occupied by any Moricandia species. As a consequence of this jump, the morphological disparity between the spring and summer phenotypes of M. arvensis, calculated as their distance in the morphospace, was very high (0.264). In fact, it was much higher than the average pairwise disparities among all studied Brassicaceae species (0.155 0.090, mean SE, 4912 545 pairwise disparities) and almost 50% of the largest observed disparity (0.55). This outcome suggests that phenotypic plasticity prompts M. arvensis to explore two distant regions of the Brassicaceae floral morphospace almost simultaneously. Although several polymorphic species showed considerable values of between-morph disparity, they were significantly smaller than the disparity between spring and summer floral phenotypes of M. arvensis (Z-score =5.06, P< 0.0001; Table S4). The plasticity-mediated disparity of M.arvensis was significantly higher than the disparity existing between Moricandia species (0.057 0.033, mean 1 SE, Z-score =6.27, P< 0.0001; Table S5) and between species belonging to the same genus (0.069 0.055, Z-score =3.51, P< 0.0002) (Fig. 3). It was marginally different from the disparity between species of Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation New Phytologist (2021) www.newphytologist.com New Phytologist Research 5 different genera but the same tribe (0.150 0.085, Z-score = 1.34, P= 0.089) and it was statistically similar to the disparity occurring between species belonging to different tribes (0.167  0.087, Z-score =1.11, P= 0.133; Fig. 3). Plasticity-mediated floral divergence Disparities to the most recent common ancestor and to the direct ancestors were much larger for the summer floral phenotype of M. arvensis than for the spring phenotype and the other species of Moricandia (Table 1). This result was consistent across phylogenies (Fig. 4a,b) Consequently, the spring phenotype did not significantly diverge from either the MRCA of Moricandia (Zscore =0.36, P= 0.36) or from its direct ancestor (Z-score = 1.24, P= 0.108). By contrast, the summer phenotype of M. arvensis diverged significantly both from the Moricandia MRCA (Z-score =2.48, P= 0.007) and from its direct ancestor (Z-score =1.77, P= 0.038). Hence, the summer phenotype explores a region of the floral morphospace located out of its phylogenetic clade range (Fig. 4a). The ancestral disparity of the summer phenotype was even higher than the ancestral disparity of most other Brassicaceae species (Fig. 4b). Morphospatial variation in pollination niches The plant–pollinator network included in this study was significantly modular (Modularity =0.385, P< 0.0001) and identified eight different pollination niches associated with different groups of pollinators (Fig. S2) located in different regions of the morphospace (F=44.4, P<0.001, R 2 =0.39; Adonis test; Fig. 5; Table S6). Because different insects visited M. arvensis in spring and summer (Table S7), this plant species shifted between pollination niches seasonally (Fig. 5). During spring, M. arvensis belonged to a niche where the most frequent pollinators were long-tongued bees, beeflies, and hawkmoths (pollination niche 5 in Fig. 5). This pollination niche was also shared by the other Moricandia species (Fig. 6) and by 60 out of the 500 species we were able to include in this analysis (Fig. 5). By contrast, during summer, M. arvensis belonged to a niche dominated by short-tongued bees (pollination niche 3 in Fig. 5). This niche was shared with 116 other species. This niche shift was substantial, since the overlap between the spring and summer pollinator niches of M. arvensis (Czekanowski overlap index =0.35) was significantly lower than the niche overlap between congeneric species of Brassicaceae (0.57 0.42, Z-score =0.51, P=0.003). Evolution of pollination niche in Moricandia lineage The stochastic character mapping analysis suggests that the ancestral niche of the Moricandia lineage was pollination niche 5 (Fig. 6). This niche seems to be ancestral to the origin of the genus Moricandia, since it is shared with the closest genera, Rytidocarpus and Eruca (Fig. 6). Two shifts have occurred later during the evolution of the Moricandia lineage, both towards pollination niche 3. One has been made by Moricandia foetida, and the other shift is that made by the summer phenotype of M.arvensis (Fig. 6). –0.2 –0.1 0.0 0.1 0.2 –0.2 –0.1 0.0 0.1 0.2 NMDS1 NMDS2 Chorispora greigii Heliophila latisiliqua Solms-laubachia gamosepala Streptanthus hispidus (a) n = 3140 spp Stress = 0.210 R = 0.956 2 Plasticity-mediated disparity 0.96 Linear R 0.80 Observed dissimilarity Observed distance Moricandia arvensis spring M. moricandioides M. nitens M. rytidocarpoides M. sinaica M. spinosa M. suffruticosa M. foetida M. arvensis summer Non metric R (b) Fig. 2 Brassicaceae floral morphospace. (a) Floral morphospace of the Brassicaceae, shown as the projection of 31 traits recorded in 3140 species onto two nonmetric multidimensional scaling (NMDS) axes. The positions of the spring and summer phenotypes of Moricandia arvensis are linked by a thick, dashed line. We have also indicated the position of the other Moricandia species, marked as small black dots. (b) Shepard plot showing the goodness of fit of the NMDS ordination. New Phytologist (2021) www.newphytologist.com Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation Research New Phytologist 6 Plasticity-mediated floral convergence The three methods gave similar results irrespective of the phylogeny used (Table 2). Brassicaceae species belonging to pollination niche 3 converged morphologically according to the difference between the observed and expected angles from their most recent common ancestor, the C1–C3 distance measures of convergence, and the Wheatsheaf metric only for Huang et al. (2020)’s phylogeny (Table 2). Similarly, the species belonging to pollination niche 5 also converged morphologically according to the difference between the observed and expected angles from their most recent common ancestor, and the C1–C3 distance measures of convergence (Table 2). Altogether, these results suggest that floral convergence was frequent among the species belonging to either of the two studied niches, irrespective of the method and the phylogeny used. A convergence region including all species of Moricandia was detected only when using the Smith & Brown (2018) and Gaynor et al. (2018) phylogenies (Table 3). This region extended across the morphospace to include two species of Erysimum also belonging to pollination niche 5 (Table 3). However, it did not include the summer phenotype of M. arvensis in any case (Fig. 7). When using Huang et al. (2020)’s phylogeny, including only species from the tribe Brassiceae, the analysis was unable to detect a convergence region including the entire Moricandia genus (Table 3). The C5 index detected between two and six convergent events towards pollination niche 5, depending on the phylogeny used (blue arrows in Fig. 7), but none was associated with the spring phenotype of M. arvensis. By contrast, the C5 index consistently detected that the summer phenotype of M. arvensis has converged with the species belonging to pollination niche 3 (red arrow in Fig. 7). Discussion Brassicaceae floral diversification Brassicaceae is a plant family with a conserved floral bauplan consisting of flowers with four sepals, four petals, six stamens, and two carpels (Appel & Al-Shehbaz, 2003) that are easily recognizable due to their cross-like appearance (Endress, 1992; Ronse de Craene, 2010). The floral morphospace built here, including >80% of the known Brassicaceae species, has revealed considerable diversity in floral forms due to the considerable variation in morphological, architectural, and colour traits. We presume that the magnitude of floral variation is even larger than that reported here, since information on other floral traits, such as pollen production, intrafloral phenology, morphology and number of nectaries, production and composition of nectar, and amount and composition of floral scents and volatiles, is not available for many of the species included in our study. Brassicaceae flowers seem to have diversified considerably while the number of floral organs has remained relatively invariant (Bowman et al., 1999; Gomez et al., 2016; Nikolov, 2019). The floral morphospace resulting from this diversification process is dense and lacks large empty regions. This suggests that Brassicaceae has explored many of the possible floral morphospace regions during its evolution, including those associated with the loss or gain of some floral organs (Bowman et al., 1999; Lee et al., 2002; Mendez & Gomez, 2006; Hameister et al., 2013; Nikolov, 2019). The macroevolutionary processes shaping this floral diversification are Between Moricandia species Between morphs of the same species Between species of the same genus Between species of different genera Between species o f different tribes Pollination-traits disparity 28 34 Disparity between spring and summer phenotypes of Moricandia arvensis 0.0 0.1 0.2 0.3 0.4 0.5 136 768 4578 031 147 694 Fig. 3 Floral disparity. Magnitude of floral disparity between different taxonomic levels of Brassicaceae species. The horizontal lines in the boxplots denote the median values (50 th percentile); the boxes contain the 25 th –75 th percentiles of the dataset; the whiskers mark the 5 th and 95 th percentiles; and values beyond these upper and lower bounds are considered outliers, marked with dots. The number above each boxplot shows the number of disparities per level. We have compared this value with the disparity between spring and summer phenotypes of Moricandia arvensis (comparisons with boxplots in red are statistically significant at P < 0.05, in orange are marginally significant at P< 0.1, and in grey are not significant). Table 1 Floral divergence: floral disparity of each species of Moricandia from the most recent common ancestor (MRCA) of the genus and from the direct ancestor of each species, according to the Moricandia phylogeny made by Perfectti et al. (2017). Species Disparity to MRCA Disparity to direct ancestor Moricandia foetida 0.0396 0.1399 Moricandia moricandioides 0.0591 0.0373 Moricandia nitens 0.0417 0.0735 Moricandia rytidocarpoides 0.0258 0.1050 Moricandia sinaica 0.0276 0.1055 Moricandia spinosa 0.0196 0.2096 Moricandia suffruticosa 0.0639 0.2070 Moricandia arvensis spring phenotype 0.0808 0.0239 M. arvensis summer phenotype 0.1951 0.2874 Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation New Phytologist (2021) www.newphytologist.com New Phytologist Research 7 difficult to deduce due to the absence of a well-solved familywide phylogeny comprising the species used in this study. Nevertheless, we could infer some of these processes by checking for the position of related species in the morphospace. We found that congeneric species were often widely spread across the morphospace, whereas species from different genera and tribes were located close together in the same regions. This pattern strongly suggests frequent events of floral divergence and convergence during the evolution of the family (Pie & Weitz, 2005; Chartier et al., 2014). In fact, divergent, convergent, and parallel evolutionary patterns seem to be frequent in Brassicaceae (Huang et al., 2016; Zhong et al., 2019; Rellstab et al., 2020; Yang et al., 2020; Bohutınskaet al., 2021a,b). Several biotic and abiotic factors may have triggered the evolution of floral diversity in Brassicaceae. Pollinators constitute an important biotic factor that explains floral diversification in many angiosperms (van der Niet & Johnson, 2012). It has been previously shown that, despite exhibiting a generalist pollination system, different Brassicaceae species are pollinated by different sets of pollinators and exploit different pollination niches (Wilson et al., 1999; Denisow, 2004; Gomez et al., 2014, 2015a,b, 2016; Schlinkert et al., 2015; Gibson-Forty et al., 2020). Pollinators have indeed contributed to the floral divergence of some Brassicaceae species and clades (Gomez et al., 2015b, 2016; Gervasi & Schiestl, 2017; Schiestl et al., 2018; Ramos & Schiestl, 2019). How they have contributed to the floral convergent and divergent evolution at the level of the whole family is difficult to ascertain with the information at hand. But the fact that pollination niches are not randomly distributed across the floral morphospace suggests that pollinators may have played a role in the diversification pattern of Brassicaceae flowers. Our study suggests that, despite the widespread generalization observed in the pollination systems of Brassicaceae, species converging phenotypically interact with similar pollinators. Further research is nevertheless necessary to get an accurate conclusion. Plastic floral divergence in M. arvensis Moricandia is a small and morphologically uniform genus, with most species displaying similar medium-to-large, cross-shaped purple-to-pale lilac flowers (de Bolos, 1946; Tahir & Watts, 2011; Perfectti et al., 2017), with only two narrow endemic species departing slightly from this floral pattern, M. foetida and Moricandia rytidocarpoides (Perfectti et al., 2017). Consequently, the magnitude of floral disparity among these species is small, with all of them located in the same upper-left region of the Brassicaceae floral morphospace. They all were also located very close to the region where the most recent common ancestor of the genus was presumably located. This suggests that during the almost 7 Myr of evolution of the genus Moricandia (Perfectti et al., 2017), floral divergence has been negligible. In contrast to this phylogenetic conservatism, we found that the summer phenotype of M.arvensis was located in a completely different region of the floral morphospace. This is not surprising, considering that the floral forms of this polyphenic species switch between two extremes (Gomez et al., 2020). Plasticity-mediated movements across the morphospace have been reported for several other organisms (Pfennig, 2021), like fish (Muschick et al., 2011; Lofeu et al., 2021), ants (Oettler et al., 2019), nematodes (Susoy et al., 2015) and even plants (Barnett et al., 2018). Our case is remarkable because, in contrast to the findings of (a) –0.2 –0.1 0.0 –0.2 –0.1 0.0 0.1 0.2 M. arvensis summer Moricandia arvensis spring M. nitens M. moricandioides M. spinosa M. sinaica Most common recent ancestor of Moricandia M. rytidocarpoides M. foetida 0.0 0.1 0.2 0.3 0.4 Phylogenies Disparity to direct ancestor (b) Summer phenotype Spring phenotype 12 43 NMDS1 NMDS2 M. suffruticosa Fig. 4 Floral divergence. (a) Phylomorphospace of the Moricandia species, using the phylogeny of Perfectti et al. (2017). (b) Floral disparity between each of the two phenotypes (i.e. spring and summer) and the nearest ancestor, according to the four dated phylogenies used in this study (1, Perfectti et al., 2017; 2, Smith & Brown, 2018; 3, Gaynor et al., 2018; 4, Huang et al., 2020). The horizontal lines in the boxplots denote the median values (50 th percentile) of the disparity of each Brassicaceae species included in each phylogeny and its direct ancestor; the boxes contain the 25 th –75 th percentiles; the whiskers mark the 5 th and 95 th percentiles; and values beyond these upper and lower bounds are considered outliers, marked with dots. New Phytologist (2021) www.newphytologist.com Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation Research New Phytologist 8 previous studies, plasticity is expressed in M. arvensis within individuals (Gomez et al., 2020), meaning that the plasticitymediated movement across the morphospace is made by single individuals as the season changes. Moreover, the magnitude of this plasticity-induced within-individual disparity was even larger than that found between Moricandia species, and even between species belonging to different genera or tribes. More importantly, the plastic flowers of M. arvensis explored a region of the floral morphospace outside its phylogenetic clade range, diverging radically from the ancestral phenotype. Finally, within-individual plasticity is expressed by most populations across the whole geographical range of the species (Gomez et al., 2020; J. M. Gomez et al., unpublished). Altogether, these findings suggest that plasticity mediates the pattern of floral divergence in M. arvensis and contributes to expanding the area of the floral morphospace explored by Moricandia. Studying this process in other crucifers will help establish if it is exclusive of Moricandia or, on the contrary, it is a motor of floral diversification in Brassicaceae. Plasticity-mediated floral divergence has allowed the summer morph of M. arvensis to explore a new region of the ecological space and to exploit a pollination niche that differs markedly from that exploited by its closest relatives. The flowers of most Moricandia species are pollinated by long-tongued bees (Apidae: Anthophorini) and, to a lesser extent, by other long-tongued insects (Dukas & Shmida, 1989; K€uchmeister et al., 1995; Gonzalez-Megıas, 2016; Gomez et al., 2020). This niche is even shared by their closest relatives, such as Rytidocarpus moricandioides or Eruca spp. (Gomez et al., 2016, 2020; Barazani et al., 2019; Shakeel et al., 2019), and was reconstructed as the most likely pollination niche of the Moricandia ancestor. However, the summer floral morph of M. arvensis is pollinated by a different set of insects and belongs to a different pollination niche associated with short-tongued small and medium-size bees (Gomez et al., 2020). Consequently, this morph exploits a different pollination niche that has largely diverged from the ancestral Moricandia niche. We do not know yet whether the niche shift is adaptive, with plants expressing during summer a set of floral traits that attract efficient pollinators and maximize fitness, or is a mere consequence of the change in floral phenotype, with the pollinators visiting summer flowers being those attracted to the type of flower displayed during summer. Interestingly, the summer pollination niche is also exploited by other co-generic –0.2 –0.1 0.0 0.1 0.2 –0.2 –0.1 0.0 0.1 0.2 Niches Main pollinators 26 37 15 48 Niches Main pollinators NMDS1 NMDS2 Fig. 5 Pollination niches. Morphospatial distribution of the eight pollination niches detected in Brassicaceae. Insect silhouettes were drawn by Divulgare (www.divulgare.net) and are reproduced under a Creative Commons license (http://creativecommons.org/licenses/by-nc-sa/3.0). E. pinnatifida E. vesicaria Eruca foleyi M. arvensis M. spinosa M. nitens M. suffruticosa M. moricandioides M. rytidocarpoides Moricandia foetida Rytidocarpus moricandioide s 1.00 1.00 1.00 1.00 1.00 0.97 0.97 0.91 1.00 Fig. 6 Evolution of pollination niches. Estimate of the ancestral pollination niche of the Moricandia lineage using a stochastic character mapping inference analysis. The numbers underneath each ancestral node indicate the posterior Bayesian probability of belonging to pollination niche 5. Pink dots, Pollination niche 5; White dots, Pollination niche 3. Ó2021 The Authors New Phytologist Ó2021 New Phytologist Foundation New Phytologist (2021) www.newphytologist.com New Phytologist Research 9