Information

19.3: Adaptive Evolution - Biology

19.3: Adaptive Evolution - Biology


We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

19.3: Adaptive Evolution

Speciation and adaptive evolution reshape antioxidant enzymatic system diversity across the phylum Nematoda

Nematodes have evolved to survive in diverse ecological niches and can be a serious burden on agricultural economy, veterinary medicine, and public health. Antioxidant enzymes in parasitic nematodes play a critical role in defending against host oxidative stress. However, the features of the evolution of antioxidant enzymes in the phylum Nematoda remain elusive.

Results

Here, we systematically investigated the evolution and gene expression of antioxidant enzymes in the genomes of 59 nematodes and transcriptomes of 20 nematodes. Catalase has been independently lost in several orders, suggesting that it is unnecessary for some nematodes. Unlike in mammals, phospholipid hydroperoxide glutathione peroxidase is widely distributed in nematodes, among which it has evolved independently. We found that superoxide dismutase (SOD) has been present throughout nematode evolutionary process, and the extracellular isoform (SOD3) is diverged from the corresponding enzyme in mammals and has undergone duplication and differentiation in several nematodes. Moreover, the evolution of intracellular and extracellular SOD isoforms in filaria strongly indicates that extracellular SOD3 originated from intracellular SOD1 and underwent rapid evolution to form the diversity of extracellular SOD3. We identify a novel putative metal-independent extracellular SOD presenting independently in Steinernema and Strongyloididae lineage that featured a high expression level in Strongyloides larvae. Sequence divergence of SOD3 between parasitic nematodes and their closest free-living nematode, the specifically high expression in the parasitic female stage, and presence in excretory-secretory proteome of Strongyloides suggest that SOD3 may be related with parasitism.

Conclusions

This study advances our understanding of the complex evolution of antioxidant enzymes across Nematoda and provides targets for controlling parasitic nematode diseases.


Introduction

In December 2019, an outbreak of pneumonia cases in the city of Wuhan, China, was linked to a novel coronavirus. Evolutionary analysis identified this new virus to humans as a severe acute respiratory syndrome-related virus [1], in the Sarbecovirus subgenus of the Betacoronavirus genus, sister lineage to the original SARS virus subsequently named Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) to reflect this relationship to SARS-CoV [2]. This novel lineage represents the seventh known human-infecting member of the Coronaviridae. The initial outbreak of human cases of the virus was connected to the Huanan Seafood Wholesale Market in Wuhan [3], and while related viruses have been found in horseshoe bats [4] and pangolins [5], their divergence represents decades of evolution [6] leaving the direct origin of the pandemic unknown. In addition to elucidating the transmission route from animals to humans, key questions for assessing future risk of emergence are: (i) what is the extent of evolution, if any, required for a bat virus to transmit to humans, and (ii) what subsequent evolution will occur once the virus is established within the human population?

The first SARS virus outbreak in 2002/2003, causing approximately 8,000 infections, and its reemergence in late 2003, causing 4 infections, were linked to Himalayan palm civets and raccoon dogs in marketplaces in Guangdong Province [7,8]. Later, it became clear that while these animals may have been conduits for spillover to humans, they were not true viral reservoirs [9]. Extensive surveillance work subsequently identified related viruses circulating in horseshoe bats in China, some of which can replicate in human cells [10,11]. The bat viruses most closely related to SARS-CoV (hereafter referred to as SARS-CoV-1 for clarity), can bind to human angiotensin-converting enzyme (hACE2, the receptor SARS-CoV-2 also uses for cell entry), while the addition of a host protease is required to cleave the Spike protein before it can bind hACE2 for the more divergent bat viruses tested (S1 Fig) [12]. Two of the key changes which appear to have been required to generate this highly capable pathogen, the specific receptor binding domain sequence and the inserted furin cleavage site, can all be traced to bat coronaviruses [6,13–15]. Collectively, these results demonstrate that, unlike most other RNA viruses which acquire adaptations after switching to a new host species [16,17] for efficient replication and spreading as successfully as exhibited by SARS-CoV-2, the Sarbecoviruses—which already transmit frequently among bat species [18]—can exploit the generalist properties of their ACE2 binding ability, facilitating successful infection of non-bat species, including humans. A main difference between SARS-CoV-1 and SARS-CoV-2 is the increased binding affinity for hACE2 in the latter [19] permitting more efficient use of human cells and the upper respiratory tract, and on average lower severity but, paradoxically—due to the higher number of infections—higher disease burden.

There is intense interest in the mutations emerging in the SARS-CoV-2 pandemic [20–22]. Although the vast majority of observed genomic change is expected to be “neutral” [23,24], mutations with functional significance to the virus will likely arise, as they have in many other viral epidemics and pandemics [25]. In SARS-CoV-2, amino acid replacements in the Spike protein may reduce the efficacy of vaccines, replacements in proteases and polymerases will result in acquired drug resistance, and other mutations could change the biology of the virus, e.g., enhancing its transmissibility as demonstrated for the replacement D614G [26], contributing to adaption to humans, a new host species.

A main way to begin to understand the functional impact of mutations is to characterise the selective regime they are under. Mutations which are under positive selection are of particular interest as they are more likely to reflect a functional change. However, identifying mutations under positive selection from frequency data alone can be misleading, as allele frequencies in viral pandemics are significantly driven by biased sampling, founder effects, and superspreading events [22]. Exponentially growing populations can increase in average fitness [27] however, they are also expected to exhibit elevated genetic drift, with deleterious mutations “surfing” expansion waves [28].

Here, we investigate the evolutionary history of bat Sarbecoviruses that shaped the emergence and rapid spread of SARS-CoV-2 by contrasting results from comprehensive searches for signatures of positive selection (a measure of molecular adaptation) in the virus circulating in humans since the Coronavirus Disease 2019 (COVID-19) outbreak began, to signatures of historic selection acting on related bat viruses. We use an array of selection detection methods (summarised in S2 Fig) on a set of (i) 133,741 SARS-CoV-2 genomes sampled from December 2019 to October 2020 and on (ii) 69 Sarbecovirus genomes (including a representative of SARS-CoV-1 and SARS-CoV-2 sequences) separated into phylogenetically congruent regions based on detected recombination patterns [6]. Finally, we detect a shift in CpG representation on the Sarbecovirus tree, associated with the phylogenetic clade SARS-CoV-2 is found in.


Introduction

Current climate change projections suggest that the average global temperature in 2100 will be higher than the average of 1986–2005 by 0.3 °C–4.8 °C [1], coupled with an increase in temperature fluctuations in certain areas [2]. Therefore, it is now more important than ever to understand how temperature changes affect biological systems, from individuals to whole ecosystems. At the level of individual organisms, temperature affects functional traits in the form of the thermal performance curve (TPC). Typically, this TPC, especially when the trait is a rate (e.g., respiration rate, photosynthesis, growth), takes the shape of a negatively skewed unimodal curve (Fig 1) [3, 4]. The curve increases (approximately) exponentially to a maximum (Tpk) and then also decreases exponentially, with the fall being steeper than the rise. Understanding how various aspects of the shape of this TPC adapt to a changing thermal environment is crucial for predicting how rapidly organisms can respond to climate change.

(A) Tpk (K) is the temperature at which the curve peaks, reaching a maximum height that is equal to Bpk (in units of trait performance). E and ED (eV) control how steeply the TPC rises and falls, respectively. B0 (in units of trait performance) is the trait performance normalised at a reference temperature (Tref) below the peak. In addition, Wop (K), the operational niche width of the TPC, can also be calculated a posteriori as the difference between Tpk and the temperature at the rise of the TPC where B(T) = 0.5 ⋅ Bpk. We note that we use Wop instead of a metric that captures the entire TPC width because previous studies have shown that species generally experience temperatures below Tpk [6, 7]. Thus, Wop is a measure of the thermal sensitivity of the trait near typically experienced temperatures. (B) TPCs of individual- and population-level traits (such as rmax) are usually well described by the Sharpe-Schoolfield model. The raw data for panel B are available at https://doi.org/10.6084/m9.figshare.12816140.v1. TPC, thermal performance curve.

According to the Metabolic Theory of Ecology (MTE) as well as a large body of physiological research, the shape of the TPC is expected to reflect the effects of temperature on the kinetics of a single rate-limiting enzyme involved in key metabolic reactions [5, 8–11]. Under this assumption, the rise in trait values up to Tpk can be mechanistically described using the Boltzmann-Arrhenius equation: (1)

Here, B is the value of a biological trait, B0 is a normalisation constant—that includes the effect of cell or body size—which gives the trait value at a reference temperature (Tref), T is temperature (in K), k is the Boltzmann constant (8.617⋅10 −5 eV ⋅ K −1 ), and E (eV) is the thermal sensitivity of the trait at the rising component of the TPC up to Tpk. Because Tpk tends to be higher than the mean environmental temperature [6, 7, 12], E represents the thermal sensitivity within the organism’s typically experienced thermal range.

Early MTE studies argued that, because of strong thermodynamic constraints, adaptation will predominantly result in changes in B0, whereas E will remain almost constant across traits (e.g., respiration rate, rmax), species, and environments around a range of 0.6–0.7 eV [8–10]. The latter assumption is referred to in the literature as universal temperature dependence (UTD). This restricted range of values that E can take is centered on the putative mean activation energy of respiration (≈ 0.65 eV). A notable exception to the UTD is photosynthesis rate, which is expected to have a lower E value of ≈ 0.32 eV, reflecting the activation energy of the rate-limiting steps of photosynthesis [13].

The existence of a UTD has been strongly debated. From a theoretical standpoint, critics of the UTD have argued that the Boltzmann-Arrhenius model is too simple to mechanistically describe the complex physiological mechanisms of diverse organisms [3, 14–16] and is inadequate for describing TPCs emerging from the interaction of multiple factors, and not just the effects of temperature on enzyme kinetics. That is, the E calculated by fitting the Boltzmann-Arrhenius model to biological traits is an emergent property that does not directly reflect the activation energy of a single rate-limiting enzyme. For example, a fixed thermal sensitivity for net photosynthesis rate is not realistic because it depends on the rate of gross photosynthesis as well as photorespiration, which is in turn determined not only by temperature but also by the availability of CO2 in relation to O2 [17].

Indeed, there is now overwhelming empirical evidence for variation in E (thermal sensitivity) far exceeding the narrow range of 0.6–0.7 eV, with such variation being, to an extent, taxonomically structured [12, 18–23]. Furthermore, the distribution of E values across species is typically not Gaussian but right-skewed. If we assume that E is nearly constant across species—and therefore that variation in E is mainly due to measurement error—such skewness could be the outcome of the proximity of the E distribution to its lower boundary (0 eV). In that case, however, we would expect a high density of E values close to 0 eV, but such a pattern has not been observed [18]. Both the deviations from the MTE expectation of a heavily restricted range for E and the shape of its distribution have been argued to be partly driven by adaptation to local environmental factors by multiple studies. These include selection on prey to have lower thermal sensitivity than predators (the “thermal life-dinner principle”) [18], adaptation to temperature fluctuations within and/or across generations [3, 21, 24–26], and adaptive increases in carbon allocation or use efficiency due to warming [27–30].

In general then, adaptive changes in the TPCs of underlying (fitness-related) traits are expected to influence the TPCs of higher-order traits such as rmax, resulting in deviations from a UTD. Therefore, understanding how the thermal sensitivity of rmax and its distribution evolves is particularly important, as it may also yield useful insights about the evolution of the TPCs of underlying physiological traits (e.g., respiration rate, photosynthesis rate, and carbon allocation efficiency). Indeed, systematic shifts in the thermal sensitivity of fundamental physiological traits have been documented [27, 31–33], albeit not through comparative analyses of large datasets.

In particular, phylogenetic heritability—the extent to which closely related species have more similar trait values than species chosen at random—can provide key insights regarding the evolution of thermal sensitivity. A phylogenetic heritability of 1 indicates that the evolution of the trait across the phylogenetic tree is indistinguishable from a random walk (Brownian motion) in the parameter space. Note that this does not necessarily indicate that the trait evolves neutrally, as it may be under selection towards a nonstationary optimum that itself performs a random walk [34]. In contrast, a phylogenetic heritability of 0 indicates that trait values are independent of the phylogeny. This is the case either because (i) the trait is practically invariant across species and any variation is due to measurement error, or (ii) the evolution of the trait is very fast and with frequent convergence (i.e., independent evolution of similar trait values by different lineages). It is worth clarifying that rapid trait evolution that does not result in convergence (e.g., when major clades are extremely separated in the parameter space) will not lead to a complete absence of phylogenetic heritability. Phylogenetic heritabilities between 0 and 1 reflect deviations from Brownian motion (e.g., due to occasional patterns of evolutionary convergence). Among phytoplankton, measures of thermal sensitivity of rmax (E and Wop) have previously been shown to exhibit intermediate phylogenetic heritability [35]. This indicates that, among phytoplankton, thermal sensitivity is not constant but evolves along the phylogeny, albeit not as a purely random walk in trait space, reflecting either thermodynamically constrained evolution or rapid evolution in response to selection.

To understand (i) how variation in thermal sensitivity accumulates across multiple autotroph and heterotroph groups and (ii) whether its distribution is shaped by environmental selection, here we conduct a thorough investigation of the evolutionary patterns of thermal sensitivity, focusing particularly on rmax. Using a phylogenetic comparative approach, we test the following hypotheses:

1. Thermal sensitivity does not evolve across species and any variation is noise-like

In this scenario, thermodynamic constraints would force E to be distributed around a mean of 0.65 eV (or 0.32 eV in the case of photosynthesis), with deviations from the mean being mostly due to measurement error. Depending on the magnitude of the error, the E distribution would either be approximately Gaussian (little measurement error) or non-Gaussian with a high density near 0 eV (substantial measurement error). This hypothesis agrees with the UTD concept of early MTE studies. If this hypothesis holds, thermal sensitivity would have 0 phylogenetic heritability and would not vary systematically across different environments.

2. Thermal sensitivity evolves gradually across species but tends to revert to a key central value, without ever moving very far from it

This hypothesis is also consistent with the UTD assumption, as it is a relaxed variant of hypothesis 1. Here, small deviations from the central tendency of 0.65 eV (or 0.32 eV) are possible, as they would reflect adaptation of species’ enzymes to certain ecological lifestyles or niches. Therefore, thermal sensitivity would be weakly phylogenetically heritable. Thermodynamic constraints would prevent large deviations from the central tendency.

3. Thermal sensitivity evolves in other ways

This is an “umbrella” hypothesis that encompasses multiple subhypotheses that do not invoke the UTD assumption. For example, a key central tendency (thermodynamic constraint) may still exist, but its influence would be very weak, allowing for a wide exploration of the parameter space away from it. In this case, changes in thermal sensitivity could be the outcome of adaptation to different thermal environments. Another subhypothesis is that clades differ systematically in the rate at which thermal sensitivity evolves, due to the occasional emergence of evolutionary innovations. Thus, clades with high evolutionary rates would be able to better explore the parameter space of thermal sensitivity (i.e., through large changes in E and Wop values), compared to low-rate clades in which thermal sensitivity would evolve more gradually. A third subhypothesis is that evolution may favour individuals (and metabolic variants) that are relatively insensitive to temperature fluctuations. In that case, the central tendency of E would not be stationary but moving towards lower values with evolutionary time. It is worth clarifying that these 3 subhypotheses are not necessarily mutually exclusive.


Methods

Taxon sampling and pollination syndrome classification

Ethanol-preserved flowers of 30 Merianieae species, covering the major clades and morphological diversity of the tribe 25 , were used for this study (Supplementary Table 13). Our material stems from six different Latin American countries and has been collected on various sampling trips between 2002 and 2015. Due to difficulties associated with fieldwork (research permits, species occurrence in remote and isolated places), we were unable to increase sample sizes. Fifteen out of 30 species were only represented by a single specimen, the other 15 species were represented by eight specimens on average (Supplementary Table 13). Only fully anthetic and relatively undamaged flowers were used in our study (see paragraph on Estimation of missing landmarks). For 14 species, pollinators are documented and include bees (seven sp.), passerines (three sp.) and mixed assemblages of hummingbirds, bats, rodents and flowerpiercers (five sp. 25 ). For the 16 species with unknown pollinators, the syndrome classification of Dellinger et al. 25 , based on an extensive dataset of 61 floral traits not included in this study, was used. As none of the traits used for the delimitation of syndromes was used in this study, we avoid problems of circularity. Also, syndrome classification of 25 was based on rigorous field studies and objective statistical classification methods which yielded highly precise syndrome predictions. Hence, we are convinced that the risk of misclassification of species in this study is very low. In total, pollination syndromes are represented by 16 buzz-bee, eight mixed-vertebrate, and six passerine syndrome species in this study. All species have tubular anthers, releasing pollen only by a small apical pore 25 . Marked differences in pollen expulsion mechanisms differentiate the three pollination syndromes 25,28 . Stamen appendages are the key for activating pollen expulsion in the buzz-bee and passerine syndrome, while they have lost their function in the mixed-vertebrate syndrome 25 .

Phylogeny, dating and estimation of ancestral pollination syndromes

To analyse floral shape evolution across Merianieae, we inferred a Bayesian phylogeny for Merianieae using BEAST2 (v2.5.0) 52 , as implemented through the CIPRES portal 53 . We determined the best partition scheme with PartitionFinder 2 54 , using each locus as a separate probable partition, and in the case of the three coding genes, also allowing for each of the three codon positions to be considered a partition. A seven partition scheme was found to be the best fit for the data (each locus as an independent partition, and in the case of ndhF, first codon position separate from second and third position). We assigned each partition the GTR + Γ + i model of sequence evolution and unlinked the partitions. We set rate variation across branches as uncorrelated and log-normally distributed, and with tree prior set to the Yule process. Based on previous analyses across the Melastomataceae, calibrated with fossils across the Myrtales, we fixed the age of the Merianieae at 29.25 MY (Michelangeli et al., unpublished). We ran three independent analyses of 60 million generations each, sampling every 20,000 generations with a 20 % burn-in. Convergence was assessed using Tracer v.1.6 55 , and runs were considered satisfactory with effective sample size (ESS) values greater than 200. We combined the stable posterior distributions of the independent runs using LogCombiner v2.5.0 56 and a maximum clade credibility tree summarized with TreeAnnotator v2.5.0 57 .

We then pruned this tree to only include the 30 species present in this study (drop.tips PHYTOOLS 58 ). We reconstructed pollination syndromes using ML methods (ace APE 59 Supplementary Table 14) and stochastic character mapping to show that bee-pollination is ancestral (make.simmap PHYTOOLS Fig. 2a) using the ‘equal-rates’ model (lower AIC than ‘all-rates-different’). Reconstructions of pollination syndromes were later used to paint branches using OU-models (see Flower shape evolution).

High-Resolution X-ray Computed Tomography (HRX-CT), 3D-models

We prepared 137 ethanol-preserved flowers of 30 species (one to 29 flowers per species, four on average Supplementary Table 13 for exact numbers of specimens per species) for HRX-CT scanning by putting them into a contrasting agent for four weeks (1% PTA–70% EtOH, Supplementary Tables 13, 23). We then mounted fully contrasted flowers in plastic cups (Semadeni Plastics Group) and stabilized them by acrylic-pillow foam to prevent movement during the scanning process. We HRX-CT scanned the samples using the Xradia MicroXCT-200 system. We reconstructed three-D image stacks from the raw scan data (XMReconstructor XRadia Inc.) and deposited tiff-stacks on the public repository https://phaidra.univie.ac.at/ from where they can be downloaded.

Landmark placement

We used the imaging software AMIRA 5.5.0 to create 3D-models of the image stacks. We calculated isosurface models of each flower to place landmarks on. In total, we selected 37 landmarks under the criteria of homology and repeatability (ability to accurately locate homologous landmark positions in different specimens) to capture patterns of floral shape variation in the three different pollination syndromes (Fig. 1). We placed landmarks as follows: five on the typical notch on the petal tips, one at the base of the style (on top of the syncarpous ovary, not visible in Fig. 1), ten on the stamen appendage tips, ten on the base of the stamen appendages, ten on the anther pores, and one on the stigma. All landmarks were placed by one of us (S. A.) in order to minimize variation due to potential observer inconsistencies.

Statistics and reproducibility

Assessment of landmark quality. We performed all data analyses in R 60 . In order to assure accurate landmark placement and to minimize observer error, we performed a precision test at the beginning of the landmarking process for two specimens (one passerine and one hummingbird/bat pollinated) following the methodology adopted by ref. 61 . We landmarked ten replicates of the two specimens and 10 additional specimens stemming from different pollination syndromes and Procrustes fitted those three datasets separately. In optimal landmark configurations, error in replicated samples should be close to 0 and at least one magnitude smaller than in non-replicated samples. To calculate the error around each single landmark, we compared the mean distance of each landmark (of the 10 replicates and the 10 independent samples, respectively) to the consensus. Using T-and F-tests, we compared the mean replicate distances to the mean distances of the non-replicates at each landmark. All landmarks placed in both replication sets were significantly less variable than in the non-replicate placements both using T- and F-tests and observer errors (mean distance of landmarks to consensus) were more than one magnitude smaller in replicates than in the non-replicate set (set1-replicate: 0.00139, set2-replicate: 0.00117, non-replicate: 0.0689). Thus, selected landmarks were accurate enough to proceed with further landmark placement.

Estimation of missing landmarks. In 72 of the 137 specimens used for analyses, all landmarks could be placed accurately without problems. The remaining 65 specimens showed minor damages due to handling and transport or damage by herbivores or pollen thieves (e.g. broken tip of one petal, broken style tip, broken stamen or stamen tip chewed up by pollen robbing Trigona bees) so that one to maximally ten landmarks could not be placed. Most geometric morphometric analyses require the placement of exactly the same number of homologous landmarks in all specimens and are intolerant of missing data 62 . Our dataset includes a number of rare taxa collected at sites with difficult access from six different Latin American countries and excluding those from our analyses would have greatly reduced the breadth (in terms of taxonomic and morphological diversity) of our study. Since we aimed at capturing the actual 3-dimensional floral architecture of flowers, a study like ours also could not make use of flowers from herbarium specimens. We thus chose to estimate missing landmarks for the 65 specimens in questions, following methods developed by Arbour and Brown 62 . For these specimens, we estimated the missing landmarks by four different landmark estimation techniques (Bayesian PCA (BPCA), mean substitution (MS), thin-plate spline interpolation (TPS) and least-squares regression (REG)) using the R-package LOST (see ref. 61 for a thorough comparison of estimation techniques J. Arbour provided updated R scripts to run TPS in 3D, currently not implemented in LOST). To improve estimation accuracy, we only estimated missing landmarks from specimens most similar to the specimen for which landmarks should be estimated 63 . Thus, we divided the dataset of the 72 intact specimens into six subsets for estimation (first column Supplementary Table 15). For each of the subsets, we performed a test run by randomly removing one to ten landmarks in one intact individual 50 times and estimating the missing landmarks. We Procrustes fitted each estimated set and performed a PCA. We used the function protest() from the R-package ‘vegan’ to compare PCA-coordinates (first two axes) of the estimated subset and the intact subset to test if the estimation procedure significantly altered relative morphospace occupation patterns. In addition, we used T- and F-tests to test for significant alteration of each landmark position between the estimated and the intact set in all 50 runs. All estimation techniques gave PCA results that were significantly correlated to the respective intact subset but the four techniques differed in the quality of single landmark estimation (Supplementary Table 16) with MS and REG performing worst. We chose TPS as method to estimate landmarks in all 65 specimens. In order to keep possible errors due to missing data small, we estimated each specimen with missing data separately with its respective subset.

Procrustes fitting and shape space calculation. We performed generalized Procrustes superimposition of landmarks in GEOMORPH 64 to remove variation in position, orientation and size. For each species with more than one specimen present (15 species), we calculated the mean shape. For the other 15 species, which only were represented by a single specimen, we directly used the Procrustes fitted coordinates in subsequent analyses. We visualized shape space by Principal Component Analyses (PCA). In addition, we calculated phylomorphospaces using the phylomorphospace function in PHYTOOLS 58 . To visualize shape change along PC1 and PC2, we used wireframes based on codes from http://rgriff23.github.io/2017/11/10/ plotting-shape-changes-geomorph.html (last accessed 22 November 2018).

Testing hypotheses on modularity using the CR coefficient. We used the covariance ratio (CR) as a metric to test the five modularity hypotheses as it generates robust results even with small and variable sample sizes 24 . The CR-metric determines the degree of modularity between pre-defined modules (from our Hypotheses 1–5) and estimates whether they are significantly more modular than when landmarks are randomly re-assigned to modules (null-hypothesis of random trait association). The CR-coefficient ranges between 0 and positive values, smaller values indicate less covariation between partitions of data and hence modularity. We tested the five modularity hypotheses for each pollination syndrome separately but on joint Procrustes fitted landmark coordinates using the function test.modularity (GEOMORPH). Thousand random permutations were used to evaluate the statistical significance of the observed CR-coefficient.

Evaluating the strength of modularity within and among syndromes. Summary measures of trait correlation are sensitive to various attributes of the data and hence cannot be readily compared between different groups 24,33,65 such as, for instance, the three different pollination syndromes considered here. Adams and Collyer 33 proposed the z-score as a standardized test statistic for the rPLS (Partial Least Squares correlation coefficient) where the rPLS is scaled by its permutation-based sampling distribution (effect size of the rPLS is calculated as standard deviates for the permuted samples). Calculating the effect size of the difference between two rPLS effect sizes allows for direct comparison of the strength of morphological integration across datasets 33 . We extended this approach for the CR-coefficient using the formulas provided by Adams and Collyer 33 in order to statistically evaluate the strengths of modularity between the three different pollination syndromes. We performed two-sample tests to assess if levels of modularity differed significantly between pollination syndromes.

Assessing floral modularity across Merianieae. In order to understand if detected floral modules represent relatively independent units also in an evolutionary context, we tested the five different modularity hypotheses across the Merianieae phylogeny. We calculated the CR-coefficient for all species together while accounting for phylogenetic relatedness using the function phylo.modularity (GEOMORPH).

Selecting the best-fit hypothesis of floral modularity. The approaches outlined above allow for detection of modularity and an evaluation of the strength of modularity between the different pollination syndromes. However, they do not permit conclusions on which modularity hypothesis fits the data best. We thus used the maximum-likelihood approach proposed by Goswami and Finarelli 34 to assess the fit of the five competing hypotheses. First, vector congruence coefficient correlation matrices were calculated on the Procrustes fitted landmark coordinates for each pollination syndrome separately, resulting in three 37 × 37 element matrices 66 using the dotcorr function (PALEOMORPH 67 ). We then ran the function EMMLi (EMMLi 34 ) to detect the best fitting model for each pollination syndrome by comparing the finite-sample corrected Akaike Information Criterion (AICc). EMMLi allows for complex models with different correlation coefficients between and within hypothesized modules, so that a total of 15 different models were tested, including a model of no modularity. The same procedure was repeated for all species together to assess the best-fit modularity hypotheses across Merianieae.

Assessing the rate of morphological evolution. In order to understand whether different floral modules evolve at different rates (i.e. whether some traits respond to changes in pollinator selection regimes more quickly than others), we calculated multivariate net evolutionary rates under Brownian motion for each module of Hypothesis 4 and Hypothesis 5 24 . We used the function compare.multi.evol.rates (GEOMORPH).

Flower shape evolution. We calculated phylogenetic signal in flower shape on the landmark data by the Kmult statistic, which is an extension of Blomberg’s Kappa statistic and designed for multivariate data 68 . We then assessed the fit of four different evolutionary models (Brownian motion (BM), Lambda, Early Burst (EB), Ornstein-Uhlenbeck (OU)) to the landmark data using the newly developed penalized likelihood framework for highly multivariate datasets (fit_t_pl in RPANDA 35 ). Based on the clear clustering of the three different pollination syndromes in shape space as assessed by PCA, we used PC1 and PC2 to visualize flower shape change on the phylogeny by constructing a traitgram (PHYTOOLS). We then modelled trait evolution (PC1–2) under an Ornstein-Uhlenbeck (OU) process 69 to screen for different phenotypic optima within Merianieae using the l1ou R-package 36 . We used a LASSO (Least Absolute Shrinkage and Selection Operator) procedure 70 to estimate shifts in phenotypic optima from the data without an a-priori definition of where regime shifts may have occurred (estimate_shift_configuration function, estimated shifts-model). Convergence of these shifts was then evaluated using the estimate_convergent_regimes function (L1OU). We evaluated model fit using the phylogenetic Bayesian information criterion (pBIC) and calculated weights (aicw from GEIGER 71 ).

Finally, we reconstructed morphospace evolution through time on PC1 and PC2 using the evomorphospace function (EVOMAP 72 ). We did ancestral character estimation for PC1 and PC2 (ace, method REML, APE) and painted pollination syndromes onto branches according to the estimation of ancestral pollination syndromes (Fig. 2).

Assessing the robustness of the data. Since our dataset is limited in size (ca. 10% of Merianieae, 15 species only represented by one specimen), we worked towards carefully assessing the robustness of our results. First, we randomly rarefied the landmark dataset 100 times to only include one specimen per species. This rarefaction helps understand the impact of intraspecific variability (i.e. calculation of mean shape or representing each species by one specimen only). Second, we randomly down sampled the landmark dataset 100 times to 50% of species per pollination syndrome (hence, eight buzz-bee, four mixed-vertebrate and four passerine) to understand how a reduction in species numbers affects our results. Again, we only included a single specimen per species in each down sampled dataset. Note that we included four (instead of three) species in the passerine syndrome since this was the minimum number required in assessments of the best-fit modularity hypothesis.

We tested all hypotheses on modularity on these two additional datasets following the methods described above. We calculated CR- and p-values, z-scores and significant differences in strength of modularity between syndromes for each of the 100 runs. We summarized results by calculating average CR-, p- and z-scores (Supplementary Tables 3, 4) and by reporting the proportion of times a hypothesis of modularity was significant (Supplementary Tables 5, 6). Also, we assessed the best-fit modularity hypothesis for the subsampled datasets and summarized these results by counting how often a specific hypothesis resulted as best fit (Supplementary Table 8). We also used the rarefied datasets to compare rates of morphological evolution among modules (Supplementary Table 9 reporting averages) and tested which hypothesis of evolution fits the landmark data best (Supplementary Table 10). We further used rarefied datasets to estimate regime shifts under an OU-process of floral shape evolution. We summarized these results by calculating the proportion of times a species was included in a regime shift (Supplementary Table 11). Overall, both the rarefaction and the down sampling results are congruent with results obtained from the original data and support the view that the buzz-bee syndrome is most modular and that functional modularity better explains floral shape evolution in Merianieae than developmental modularity.

Reporting summary

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