Within-individual phenotypic plasticity in flowers fosters pollination niche shift

Field sampling design

To determine if there was within-individual plasticity in floral traits between spring and summer conditions, 50 plants of each of four populations from SE Spain (Supplementary Table 1) were marked at the onset of the flowering period in late February–early March 2018. The phenotype of two flowers per individual was quantified (see below). We revisited each population during summer (June 2018) and the same floral traits were quantified in the summer flowers of those plants still flowering (117 plants; Supplementary Table 1).

Experimental design

We performed an experiment testing the effect of temperature and photoperiod in floral plasticity. It included three treatments: (1) Treatment 1, where 30 plants flowered first in conditions mimicking the spring temperature and photoperiod of Mediterranean Spain (day/night = 10/14 h, temperature = 20/10 °C, average daily temperature = 14.2 °C; see Supplementary Table 1), and afterwards in conditions mimicking a mild summer (day/night = 16/8 h, temperature = 30/20 °C, average daily temperature = 23.8 °C). (2) Treatment 2, where 30 plants flowered first in spring conditions and afterwards in hot summer conditions (day/night = 16/8 h, temperature = 35/25 °C, average daily temperature = 28.8 °C). (3) Treatment 3 (control) where 15 plants flowered first in spring conditions, and afterwards they flowered again in spring conditions. For all treatments, we removed flowers before starting the second round of flowering.

We experimentally tested the occurrence of reverse plasticity by performing a Treatment 4 in which 15 plants from Treatment 1 that flowered both during spring and summer conditions were again submitted to a period mimicking spring conditions (Supplementary Table 1).

Floral traits

We measured, both in field and experimental conditions, three floral traits during spring (or under experimental spring conditions) and in summer (or under experimental hot summer conditions). These traits were corolla size, corolla shape and corolla colour.

Corolla size of each studied flower was estimated by means of two traits: (1) corolla diameter, estimated as the distance in mm between the edge of two opposite petals. (2) Corolla tube length, the distance in mm between the corolla tube aperture and the base of the sepals. These variables were measured by using a digital calliper with ±0.1 mm of error.

Corolla shape variation was studied using geometric morphometric tools based on a landmark-based methodology⁴³. For this, in each of the two selected flowers per individual plant studied in each of the four populations, we took a digital photo of the front view and planar position. We defined 32 co-planar landmarks covering the corolla shape and using midrib, primary and secondary veins and petal extremes and connections^21,44. From the two-dimensional coordinates of landmarks, we extracted shape information and computed the generalized orthogonal least-squares Procrustes averages using the generalized procrustes analysis (GPA) superposition method. Due to the intrinsic symmetry pattern exhibited by Brassicaceae flowers, we did the analyses considering both the symmetric and asymmetric components of the shape^45,46,47. We performed a principal component analysis (PCA) on the GPA-aligned specimens, and afterwards, we did a canonical variate analysis (CVA) to explore the difference in shape between season and populations^43,47. Geometric morphometric analyses were performed in the R packages ‘geomorph’⁴⁸, ‘Morpho’⁴⁷ and ‘shapes’^49,50.

To explore the relative position of the corolla shape of spring and summer flowers in the morphospace created by the species most related phylogenetically with M. arvensis, we performed a phylomorphospace. This analysis creates a plot of the main principal dimensions (the three first principal components in this case) of a tangent space for the Procrustes shape variables of the pool of species considered in the analysis and superimposed the phylogenetic tree relating this species in this plot^51,52. By doing this, this analysis reveals how the shape evolves. To perform this analysis, we collected information on the corolla shape of 72 additional species belonging to the Brassicaceae tribe Brassiceae, the tribe to which M. arvensis belongs (Supplementary Table 3). We followed the same procedure as with M. arvensis, using the same number of landmarks and computing the generalized orthogonal least-squares Procrustes averages using GPA superposition method. In this analysis, we kept separate the spring and summer flowers of M. arvensis. The phylogenetic relationship between these 72 species was obtained by making a supertree using Brassicaceae trees hosted in the repository TreeBASE Web (TreeBase.org)⁵³. We first downloaded individual phylogenetic trees from TreeBASE. Second, we concatenated all these individual trees and made a skeleton supertree. Finally, we pruned this supertree, keeping only the species included in the geometric morphometric analysis, and insert the two ‘pseudospecies’ of M. arvensis (spring and summer) as sister species. Afterwards, we projected the value of the three first components of each species on a 3D phylogenetically explicit plot. The phylogenetic analysis was performed in the R packages ‘treeman’⁵⁴, ‘phangorn’⁵⁵, ‘phytools’⁵⁶ and ‘treebase’⁵³, whereas the phylomorphospace analysis was performed in the R packages ‘geomorph’⁴⁸.

The corolla colour of M. arvensis is produced by the accumulation of flavonoids^57,58. Anthocyanin and non-anthocyanin flavonoids present in the petals of M. arvensis were analysed by ultra-performance liquid chromatography (UPLC) (ACQUITY System I-Class, Waters) coupled with quadrupole time-of-flight mass spectrometry (SYNAPT G2 HDMS Q-TOF, Waters). Analytical separation of flavonoids was performed on an Acquity HSST33 analytical column (150 mm × 2.1 mm internal diameter, 1.8 μm). A mobile phase with a gradient programme combining deionized water with 0.5% of acetic acid as solvent A and acetonitrile with 0.5% of acetic acid as solvent B was used. The initial conditions were 95% A and 5% B and a linear gradient was then established to reach 95% (v/v) of B. The total run time was 15 min and the post-delay time was 5 min. The mobile phase flow rate was 0.4 mL min⁻¹. After chromatographic separation, a high-resolution mass spectrometry analysis was carried out in positive electrospray ionization (ESI+). The ionization source parameters using high-purity nitrogen were set at 600 L h⁻¹ for desolvation gas flow and 30 L h⁻¹ for cone gas flow. Spectra were recorded over the mass/charge (m/z) range of 50–1500. Data were recorded and processed using MassLynx software. The flavonoids present in the petal extracts were characterized according to their retention times, mass spectra and molecular formula, and compared with published data when available. We calculated the relative abundance of each compound in both lilac and white petal samples (N = 5 and 2, respectively) using peak intensities.

Quantification of flavonoids present in flowers of M. arvensis was performed spectrophotometrically. Two flowers of each plant used in field and experimental studies were analysed in each blooming period. We collected the four petals of a flower. Flavonoids were extracted in 1.5 ml of MeOH:HCl (99:1% v-v) and stored at −80 °C in the dark, following the procedure described in ref. ³⁴. Two replicas of 200 μL for each sample were measured in a Multiskan GO microplate spectrophotometer (Thermo Fisher Scientific Inc., MA, USA). Main flavonoid classes present in the petals of M. arvensis are anthocyanins (cyanidin derivatives) and flavonols (kaempferol, quercetin and isorhamnetin derivatives; Supplementary Table 4)^57,58. Thus, total anthocyanins and flavonols were quantified as absorbance at 520 and 350 nm, respectively. Their concentrations were calculated using five-point calibration curves of cyanidin-3-glucoside chloride (Sigma-Aldrich, Steinheim, Germany) and kaempferol-3-glucoside standards (Extrasynthese, Genay, France) and expressed as cyanidin-3-glucoside and kaempferol-3-glucoside equivalents in fresh weight (mg g⁻¹ FW), respectively.

Objective quantification of petal colour of lilac and white petals of M. arvensis was performed by measuring their UV–Vis spectral reflectance. A petal of a flower of each colour morph (N = 10) were measured with a Jaz portable spectrometer (Ocean Optics Inc., Dunedin, FL, USA) equipped with a deuterium–tungsten halogen light source (200–2000 nm) and a black metal probe holder (6 mm diameter opening at 45°). Reflectance, relative to a white standard (WS-1-SL), was analysed with SpectraSuite v.10.7.1 software (Ocean Optics). To maximize the amount of light used in reflectance measurements and to reduce occasionally erratic reflectance values at individual nm, we set an integration time of 2 s and smoothing boxcar width of 12, respectively⁵⁹.

Foliar traits

We measured, both in field and experimental conditions, five leaf traits during spring (or under experimental spring conditions) and in summer (or under experimental hot summer conditions). These traits were the specific leaf area (SLA, m² kg⁻¹), the leaf dry matter content (LDMC, mg g⁻¹), the carbon-to-nitrogen content of leaves (C:N ratio), the isotopic signature of ¹³C in leaves (δ¹³C, ‰), and the CO₂ compensation point and the slope of the A–C_i curve.

SLA and LDMC were measured following standard protocols⁶⁰. For SLA and LDMC we collected three fully expanded and mature leaves without any visible damage (e.g., herbivory, pathogen attack) from the base, midsection and apical part of outer stems (that is, leaves were not shaded by other leaves) and at random aspects. Leaves were rehydrated overnight in the dark and subsequently weighted and scanned. Leaf area was measured using the Midebmp software (Almería, Spain). Leaves were dried in the oven at 60 °C and weighted after 72 h. From these measurements, we calculated the SLA as the one-sided area of the fully rehydrated fresh leaf divided by its dry mass, while the LDMC is the ratio between the leaf dry mass and the fully rehydrated fresh mass.

Carbon isotopic signature (δ¹³C), as well as the C and N relative content in leaves, were analysed in a couple of fully expanded leaves per plant without any visible damage. Oven-dry leaves were ground in a ball mill MM400 (Retsch GmbH, Haan, Germany) at 3000 rpm for 1 min to obtain a fine powder, which was stored in Eppendorf tubes. We wrapped 0.003 g of each sample in tin capsules D1008 (Elemental Microanalysis, United Kingdom). Leaf δ¹³C and leaf C and N relative content (in mass percentage) were determined at the Stable Isotope Analysis Lab—Centro de Instrumentación Científica (CIC) of the University of Granada (Spain) with a GC IsoLink—MS—Delta V continuous flow mass spectrometer (MS) system that includes a ISQ-QD single quadrupole MS and a gas chromatographer Trace 1310 (Thermo Fisher Scientific™, Spain). The isotopic abundance was expressed in parts per thousand (‰) as

where R_sample and R_standard are the molar ratios of heavy (¹³C) to light (¹²C) stable isotopes of the sample (R_sample) and an international standard (R_standard). MS precision was 0.15‰ for carbon, based on replicate analyses of standard reference materials.

We measured responses of CO₂ assimilation rate (A) versus calculated substomatal or intercellular CO₂ concentration (C_i) (henceforth, A–C_i curves) to determine the instantaneous photosynthetic metabolism of plants of the intermediate C₃–C₄ species M. arvensis on plants grown under the two experimental conditions (N= 22 plants, spring and hot summer conditions). Gas exchange measurements were performed on one to two mature, fully expanded leaves per plant and experimental condition using a LICOR 6400 (LI-COR Biosciences, Lincoln, USA) and following the standard recommendations to correct leakage errors^61,62,63. Cuvette conditions were maintained at a constant photosynthetic photon flux density (PPFD) of 1500 µmol m⁻² s⁻¹, a vapour pressure deficit (VPD) that ranged from 1.0 to <2.4 kPa and a cuvette temperature that was either 20 or 30 °C (see below). A–C_i curves were measured at a series of ambient (i.e., reference) CO₂ concentrations (C_a) starting at 400 ppm. The C_a was then lowered stepwise from 400 to 40 ppm and then increased again to 1600 ppm, all on 21 steps: 400, 300, 200, 150, 100, 80, 60, 50, 40, 40, 40, 60, 100, 200, 400, 600, 800, 1000, 1200, 1400, 1600 ppm of CO₂.

The CO₂ compensation point (substomatal CO₂ concentration value at which A is zero) and initial slope of each A–C_i curve were calculated performing a linear regression of 10 different C_a values (all below 300 ppm). C_a instead of C_i values were used in the calculations as we had some technical issues with the LICOR in one of the two measurement phases. This fact might overestimate absolute results for CO₂ compensation points, but the relative differences (if any) between the two experimental conditions (spring and summer) remain. Measurements of gas exchange were made at two different LICOR cuvette temperatures (20 and 30 °C) for plants growing in spring and hot summer experimental conditions, respectively. Although the CO₂ compensation point can increase with cuvette temperature⁶¹, we measured the compensation points of six plants grown in spring chamber conditions and at these two contrasting cuvette temperatures, and results did not differ (CO₂ compensation point ± s.e.m.: 29.72 ± 4.03 vs. 31.42 ± 4.03 ppm at 20 and 30 °C bulk temperature of the LICOR chamber, respectively; GLM results, F_1,4 = 0.09; P-value = 0.7796; N= 6 individuals). Moreover, mean CO₂ compensation points of plants grown in spring chamber conditions were significantly higher (not lower) than values for the same plants grown in hot-summer conditions (see main text; mean ± s.e.m.: 25.0 ± 0.21 ppm in spring versus 14.7 ± 0.16 ppm in hot summer conditions).

To compare the CO₂ compensation point of M. arvensis with pure C₃ and C₄ species, we also calculated the CO₂ compensation point for four individuals of the C₃ species Moricandia moricandioides. In addition, we have compiled from the literature information on the CO₂ compensation point of several C₄ species belonging to the genus Cleome^64,65,66. This genus is the closest C₄ plant to Moricandia⁶⁷.

Reproduction and mating system

We collected 10 ripe fruits belonging to spring flowers and 10 ripe fruits belonging to summer flowers from each of the studied plants. These fruits were taken to the lab, where we determined under magnifying glasses (×60) the total number of ovules produced per flower and the number of ripe seeds per fruit^68.

We experimentally determined the ability of plants of producing seeds by self-fertilization during each season. To avoid undesirable side effects of local conditions, this experiment was performed in controlled conditions, in the experimental chambers, at the same time. We selected in March 2019, 10 plants growing in spring conditions and 24 plants growing in summer conditions at the onset of flowering. Each newly opened flower was marked and assigned to each of the following crosses: (1) self-pollination, where the flower was hand-pollinated with own pollen; (2) outcrossing, where the flower was pollinated with pollen from one different conspecific individual from the same season. Consequently, the number of flowers per treatment varied across individual, but was always higher than two, with a mean ± 1 s.e.m. equal to 6.1 ± 0.8 flowers per plant. Because all plants were located inside the experimental chambers, they were excluded from the visit of any pollinator. At the end of the fruiting period, we counted all flowers setting fruits, collected the fruits that were taken to the lab for counting their seeds as explained in the previous section.

Spring and summer pollinators of M. arvensis

We identified the pollinators visiting the flowers of the studied populations both during spring and summer. For this, we conducted flower visitor counts in each population × season combination. We were able to obtain information in three populations (Malaha, Quesada, Tabernas). We visited each of these populations both during spring and during summer, always between 11:00 a.m. and 5:00 p.m. In these visits, we recorded the insects visiting the flowers of the plants for two hours without differentiating between individual plants. Each survey was done at least by two researchers simultaneously, sampling each population for at least 8–9 h person⁻¹. We only recorded those insects contacting anthers or stigma and doing legitimate visits at least during part of their foraging at flowers. We did not record those insects only eating petals or thieving nectar without doing any legitimate visit. Previous studies using the same methodology carried out with similar Brassicaceae species and performing rarefaction analysis indicate that a sample of 130–150 insects provides an accurate estimate of the diversity of the pollinator assemblages^21,33,69,70. To ensure that our sampling was representative, we recorded 471 insects in Malaha population, 334 insects in Quesada population and 300 insects in Tabernas population (Supplementary Table 13).

Pollinators of the relative species

We compiled information on 308,096 insect visits belonging to 38 functional groups and more than 5000 morphospecies visiting 114 plant species belonging to the same tribe than M. arvensis, the Brassiceae (Supplementary Data 3). We use our data and data from the literature. In those species studied in the field (74 species), we conducted flower visitor counts in 1–16 populations per plant species. We visited the populations during the peak of the bloom, always at the same phenological stage and between 11:00 a.m. and 5:00 p.m. In these visits, we recorded the insects visiting the flowers of the plants for two hours 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 doing legitimate visits at least during part of their foraging at flowers. We did not record those insects only eating petals or thieving nectar without doing any legitimate visit. In addition, we included information on pollinators obtained from the literature to supplement our data (Supplementary Data 3 and 6). In this case, we use the information provided in the primary literature in terms of pollinator species and abundance at flowers. In this case, the plant species included in the network do not coexist, implying that this is a clade-oriented network rather than an ecological network³².

Determination of pollination niches

In plant species with highly diverse pollination systems, like those included in this study, many pollinator species interact with the flowers in a similar manner, have similar effectiveness and exert similar selective pressures and are thus indistinguishable for the plant^10,33. These pollinators are thus grouped into functional groups, that are the relevant interaction units in generalized systems^10,21,33,40. We thereby grouped all pollinators visiting M. arvensis and the other Brassicaceae species using criteria of similarity in body length, proboscis length, morphological match with the flower, foraging behaviour, and feeding habits²¹. Supplementary Data 5 describes the 38 functional groups used in this study and Supplementary Data 4 shows the distribution of these functional groups among the studied Brassiceae taxa.

We determined the occurrence of different pollination niches in our studied populations and seasons using bipartite modularity, a complex-network metric. Modularity has proven to be a good proxy of interaction niches both in ecological networks, those included coexisting species or population, as well as in clade-oriented network, those including species with information coming from disparate and contrasting sources³². We constructed a weighted bipartite network, including pollinator data of four populations both during spring and summer flowering period. 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 by using the QuanBiMo algorithm⁷¹. This method uses a simulated annealing Monte-Carlo approach to find the best division of populations into modules. A maximum of 10¹⁰ MCMC steps with a tolerance level = 10⁻¹⁰ were 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⁷². Modularity significance was tested for each iteration by comparing the empirical versus the random modularity indices using a z-score test⁷¹. After testing the modularity of our network, we determined the number of modules⁷³. We subsequently identified the pollinator functional groups defining each module and the plant species that were ascribed to each module. Modularity analysis was performed using R package bipartite 2.0⁷⁴.

Pollinator preference experiments

We carried out two experiments, one in late June 2019 and the other in late February/early March 2020, where 10 plants displaying spring-type flowers and 10 plants displaying summer-type flowers were offered at the same time to pollinators in a natural M. arvensis population. Plants were experimentally grown in chambers under each of the levels of the Treatment 2 conditions (spring conditions and hot summer conditions) and taken to the field in pots. Plants with each type of flowers were randomly distributed in a 5 × 4 grid with plants separated one metre. We performed ten 2-h trials in 2019 and nine 2-h trials in 2020, where the abundance and identity of the insects visiting the flowers of each experimental plant were recorded. Each trial was done by two researchers simultaneously, totalling 40 h of observation in 2019 and 36 h in 2020. Trials were performed between 11:00 and 13:00 local time in 2019 and between 12:00 and 14:00 local time in 2020. We counted each day the number of open flowers per plant, and, to control per between-plant differences in number of open flowers, insect abundance was expressed in number of visits per flower per hour. Afterwards, we followed the approach explained in the previous section to build up the bipartite interaction network, determine modularity and obtain the different modules, using only those plants having more than 10 visits per flower per hour (18 plants in the 2019 experiment and 20 plants in the 2020 experiment).

De novo transcriptome analysis

In five plants from Treatment 1 (labelled MAR-70, MAR-81, MAR-83, MAR-98, and MAR-120), we sampled one starting-to-open flower bud during period 1 (experimental spring conditions) and another flower bud during period 2 (experimental mild summer conditions). We introduced the flower buds in liquid N₂ and conserved them at −80 °C until RNA extraction. We used the RNeasy Plant Mini Kit (Qiagen) for total-RNA extraction. We checked the quality and quantity of the extracted RNA with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The 10 RNA samples (five individuals, two conditions) were sequenced by Macrogen Inc. (Seoul, South Korea) after treatment with RiboZero to remove ribosomic RNA⁷⁵. Libraries were produced using the TruSeq Stranded Total RNA LT Sample Prep Kit (Plant) and sequenced in an Illumina Novaseq6000 platform run (paired-end 150 bp) for a minimum of 40M reads/sample. See Supplementary Table 12 for information on the total number of bases reads and the total number of reads. RNA-Seq raw reads were submitted to Sequence Read Archive with the project accession number PRJNA604514. We retrieved 171,210 trinity genes (Supplementary Table 12 and Supplementary Data 1), but only 47,440 passed the criteria for inclusion in the analyses (at least 10 counts per million in spring or summer conditions).

We checked the quality of the raw sequences with FastQC⁷⁶. After that, we used cutadapt (v. 1.15)⁷⁷ and sickle (v. 1.33)⁷⁸ to quality trim and remove adaptors. To produce a reference transcriptome, we concatenated the libraries (only paired sequences) and submitted them to RNA-Seq de novo assembly using Trinity v2.8.4⁷⁹. A total of more than 372 Mbp were assembled in 424,981 transcripts (171,210 trinity ‘isogenes’), with a median contig length of 499 bp, an average length of 876.97 bp, and 10% of the assembled nucleotides appearing in transcripts of more than 4257 bp. We used trinity add-on scripts to compute contigs statistics with the help of Bowtie2 (v. 2.3.4.1)⁸⁰. A total of 97.24% of reads aligned to the reference, with >80% aligned more than one time. We performed a BLASTX search of the assembled transcriptome to the SWISSPROT database and found that 10,696 proteins of this database mapped in more than 90% of their sequences.

Analysis of DEG

We estimated the abundance of transcripts and genes using the Trinity abundance_estimates_to_matrix.pl script. Since the objective was to compare the within-individual expression of genes, we produced a matrix of raw gene counts as the input for the following analyses. We normalized (TMM)⁸¹ and filtered (conserving genes with at least 10 transcripts per gene and treatment) the matrix using edgeR⁸². A total of 47,870 genes passed this criterion. We analysed this matrix using a design of repeat samples with two conditions (spring as control and summer) and fit a negative binomial generalized log-linear model to the read counts for each gene. We selected as DEG the genes with false discovery rate adjusted-P values < 0.05.

Annotation and GO-enrichment analysis

We used Trinotate v3⁸³ to annotate the 47,870 genes using the uniport_sprot and Pfam-A databases, and mapping the longest isoform of each of those genes. Additionally, we used sma3s⁸⁴ with similar results. For the gene ontology (GO) enrichment analyses, we used the Bioconductor package goseq⁸⁵. As this package requires gene lengths, we obtained them using the Trinity script TPM_weighted_gene_length.py. We compared the GO terms of DEG with those of the 47,870 genes that passed the previous filtering.

Statistical analyses of floral integration

Floral integration was computed as the relative variance of eigenvalues of the covariance matrix of the five floral traits included in this study⁸⁶. One aspect of integration is that variation is concentrated in one or a few of the available dimensions⁸⁷. As a consequence, there will be one or a few large and many small eigenvalues for the covariance matrix of integrated data, whereas eigenvalues of the covariance matrix will be more homogeneous for data lacking integration. To control for among-species differences in sampling size, we re-scaled the relative variance of eigenvalues by the total variance and number of dimensions⁸⁷. By doing this, corolla shape integration ranges between 0 and 1 and can be interpreted as the percentage of integration regarding the maximum possible integration. This index is thus directly comparable to other integration indices found using different approaches.

Statistical analyses of within-individual floral plasticity

The magnitude and significance of the within-individual plasticity were calculated for each floral and foliar trait by random slope mixed models, an analysis that fit individual-level reaction norms and thus assesses their variation in a single step⁸⁸. We considered in all these models as explanatory variable the average daily temperature of each population or treatment during each season (spring and summer) or experimental condition (Supplementary Table 1). We first determined the overall population-level average effect of temperature on each response variable including temperature as a fixed effect function. Second, we quantified how much variation there is among individuals around the average population-level function by including in a second model the individual as a random effect (random intercepts). Finally, we quantified the variation around the average responses in the slopes of the individual reaction norms by adding a random regression term to the model. The temperature was included in all these models mean-centred⁸⁸. Because we compared phenotypes between two environments, we always fitted linear models. The magnitude and significance of the population-level plasticity were found by estimating the coefficient of the fixed effect function. The among-individual differences in the magnitude of each trait were statistically tested by comparing the log-likelihoods of the first and second models and performing a likelihood ratio test (LRT)⁸⁸. Similarly, the among-individual differences in the slopes of their reaction norms (the occurrence of G × E interaction) were statically tested by comparing the log-likelihoods of the second and third models and performing a LRT. The goodness-of-fit of these three models were compared by means of their Akaike Information Criteria (AIC)⁸⁸.

Bayesian generalized multivariate multilevel models

Plasticity integration was determined by quantifying the correlations between the plasticities of pair of traits. For this, we used Bayesian generalized multivariate multilevel models⁸⁹. The structure of these models was similar to that explained previously for random regression, including individuals as random intercepts and slopes. As the dependent variable, we used a composite variable including all the five floral traits analysed independently in the previous analyses. We scaled and centred both the independent and dependent variables and used weakly informative prior with normal distribution and centred on zero⁹⁰. We ran four chains with 2000 iteration each, burning 1000 samples per chain. In total, we analysed 400 post-warmup samples. We got in all cases the potential scale reduction factor on split chains to be 1 or very close to 1 at convergence. Significance was obtained for each effect by means of the posterior distribution of the 95% credible interval of its mean estimate. Bayesian analyses were performed in the R packages ‘brm’^89,91.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com