Phenotypic plasticity triggers rapid morphological convergence

Phenotypic convergence, the independent evolution of similar traits, is ubiquitous in nature, happening at all levels of biological organizations and in most kinds of living beings. Uncovering its mechanisms remains a fundamental goal in biology. Evolutionary theory considers that convergence emerges through independent genetic changes selected over long periods of time. We show in this study that convergence can also arise through phenotypic plasticity. We illustrate this idea by investigating how plasticity drives Moricandia arvensis, a mustard species displaying within-individual polyphenism in flowers, across the morphological space of the entire Brassicaceae family. By compiling the multidimensional floral phenotype, the phylogenetic relationships, and the pollination niche of over 3000 Brassicaceae species, we demonstrated that Moricandia arvensis exhibits a plastic-mediated within-individual floral disparity greater than that found not only between species but also between higher taxonomical levels such as genera and tribes. As a consequence of this divergence, M. arvensis moves outside the morphospace region occupied by its ancestors and close relatives, crosses into a new region where it encounters a different pollination niche and converges phenotypically with distant Brassicaceae lineages. Our study suggests that, by inducing phenotypes that explore simultaneously different regions of the morphological space, plasticity triggers rapid phenotypic convergence.


305
We used the Gower distance, the number of mismatched traits over the number of shared traits.

306
This dissimilarity index is preferable to the raw Euclidean distance when there are discrete and 307 continuous traits co-occurring in the same dataset (52).

308
We reduced the dimensionality of this phenotypic matrix by projecting it in a two-dimensional 309 space. For this, to ensure an accurate description of the distribution of the species in the 310 morphospace, we first run a principal coordinate analysis (PCoA), a technique providing a 311 Euclidean representation of a set of objects whose relationship is measured by any dissimilarity 312 index. We corrected for negative eigenvalues using the Cailliez procedure (25). Afterwards, we 313 used this metric configuration as the initial configuration to run a non-metric multidimensional 314 scaling (NMDS) algorithm (25), a method that will further optimise the sample distribution so as 315 more variation in species composition is represented by fewer ordination axes. Unlike methods 316 that attempt to maximise the variance or correspondence between objects in an ordination,

317
NMDS attempts to represent, as closely as possible, the pairwise dissimilarity between objects in 318 a low-dimensional space. NMDS is a rank-based approach, where the original distance data is 319 substituted with ranks, preserving the ordering relationships among species (25). Objects that are 320 ordinated closer to one another are likely to be more similar than those further apart (53). This 321 method is more robust than distance-based methods when the original matrix includes variables 322 of contrasting nature. However, NMDS is an iterative algorithm that can fail to find the optimal 323 solution. We decreased the potential effect of falling in local optima by running the analysis with 324 5000 random starts and iterating each run 1 x 10 6 times (54). The NMDS was run using a 325 11 monotone regression minimizing the Kruskal's stress-1 (55, 56), and compared each solution 326 using Procrustes analysis, retaining that with the lowest residual. Because many species did not 327 share trait states, a condition complicating ordination, we used stepacross dissimilarities, a 328 function that replaces dissimilarities with shortest paths stepping across intermediate sites while

329
regarding dissimilarities above a threshold as missing data (57). Furthermore, we used weak tie 330 treatment, allowing equal observed dissimilarities to have different fitted values. The scores of the 331 species in the final ordination configuration were obtained using weighted averaging. We checked 332 if the reduction in dimensionality maintained the between-species relationship by checking the 333 stress of the resulting ordination and finding goodness of fit measure for points in nonmetric 334 multidimensional scaling (54). Both PCoA and NMDS ordinations were done using the R package 335 vegan (58) and ecodist (59). It is important to note that, although the transfer function from 336 observed dissimilarities to ordination distances is non-metric, the resulting NMDS configuration is

337
Euclidean and rotation-invariant (60). were calculated using the function dispRity of the R package dispRity using the command 350 centroid (62). We checked whether the disparity between spring and summer M. arvensis 351 phenotypes was significantly different from the disparities of each of these sets of species using 352 Z-score tests.

354
Family-wide phylogeny. We retrieved 80 phylogenetic trees from the literature and from the 355 online repositories TreeBase (Table S9). All trees were downloaded in nexus format. The 356 taxonomy of the species included in each tree was checked and updated using the species 357 checklist with accepted names provided by Brassibase (https://brassibase.cos.uni-heidelberg.de/) 358 (7, 23, 63). All trees were converted to TreeMan format (64) and concatenated into a single

359
TreeMen file that was then converted into a multiPhylo class. Afterward, we estimated a 360 supertree from this set of trees. Because trees did not share the same taxa, we used the Matrix 12 representation parsimony method (65). To make this supertree more accurate, it was re-362 constructed using as backbone phylogeny the tree provided by Walden et al. (7). We removed 363 from the supertree those species without information on floral phenotype, resulting in a tree with 364 1876 taxa. Because the original trees used to assemble this supertree where very 365 heterogeneous, this supertree was not dated. We finally rooted the supertree using several 366 species belonging to the sister families Capparaceae and Cleomaceae (66). All phylogenetic 367 manipulations were performed using the R libraries treebase (67), ape (68), treeman (64), 368 phangorn (69) and phytools (70).

369
We tested whether the position of the Brassicaceae species in the morphospace was

384
We counted the number of intersections between lineages as a measurement of the 385 disorder of the phylomorphospace and evidence of the mode of evolution of the phenotypes (11).

386
For this, we used R codes provided in Ref 11. We compared the observed number of crossings    Table S10). In those species studied by us 415 (coded as UNIGEN data origin in the Supplementary Data 3), we conducted flower visitor counts 416 in 1-16 populations per plant species. We visited the populations during the blooming peak, 417 always at the same phenological stage and between 11:00 am and 5:00 pm. In these visits, we 418 recorded the insects visiting the flowers for two hours without differentiating between individual 419 plants. Insects were identified in the field, and some specimens were captured for further 420 identification in the laboratory. We only recorded those insects contacting anthers or stigma and 421 doing legitimate visits at least during part of their foraging at flowers. We did not record those

454
Carlo approach to find the best division of populations into modules. A maximum of 10 10 MCMC 455 steps with a tolerance level = 10 -10 was used in 100 iterations, retaining the iterations with the 456 highest likelihood value as the optimal modular configuration. We tested whether our network was 457 significantly more modular than random networks by running the same algorithm in 100 random 458 networks, with the same linkage density as the empirical one (81). Modularity significance was 459 tested for each iteration by comparing the empirical versus the random modularity indices using a 460 Z-score test (80). After testing the modularity of our network, we determined the number of

515
The second approach to measure convergence was based on comparing the angles 516 formed by two tested clades from their most recent common ancestor with the expected angle 517 according to null evolutionary models (38). Under the "state case", search.conv computes the 518 mean angle over all possible combinations of species pairs using one species per state. Each 519 individual angle is divided by the patristic distance between the species. Significance is assessed 520 by contrasting this value with a family of 1,000 random angles obtained by shuffling the state 521 across the species (38). These analyses were performed using the R package RRphylo (90).

522
The third approach to measure convergence used the Wheatleaf metric (40). This index 523 generates phenotypic (Euclidean) distances from any number of traits across species and 524 penalizes them by phylogenetic distance before investigating similarity (in order to weight close 525 phenotypic similarity higher for distantly related species). It uses an a priori designation of 526 convergent species, which are defined as species belonging to a niche for which the traits are 527 hypothesized to converge. The method then calculates a ratio of the mean (penalized) distances 528 between all species to the mean (penalized) distances between allegedly convergent species.

529
The index detects if convergent species diverge more in phenotypic space from the non-   Jorquera, and Iván Rodríguez Arós for helping us during several phases of the study. We also 539 thank all contributors to the pollinator database (Table S10) Table S8). (C)

810
Floral disparity to the nearest ancestor, according to the supertree and 18 time-calibrated 811 phylogenies (phylogeny codes in Table S8). We show the disparity between the two M. arvensis        Table S7. Outcome of the analyses testing for morphological convergence between the Moricandia clade and the rest of clades included in each time-calibrated phylogeny. Clade size is the number of species within the Moricandia clade. θ real is the mean angle over all possible combinations of pairs of species taking one species per clade. θ ace is the mean angle between ancestral states between each pairs of clades. dist mrca is the patristic distance (sum of brach length) between the most recent common ancestors of each pair of clade. We indicate the congervent clades and the pollination niches of each species included in the convergent clades.