Correlated evolution between climate and suites of traits along a fast–slow continuum in the radiation of Protea

Abstract Evolutionary radiations are responsible for much of Earth's diversity, yet the causes of these radiations are often elusive. Determining the relative roles of adaptation and geographic isolation in diversification is vital to understanding the causes of any radiation, and whether a radiation may be labeled as “adaptive” or not. Across many groups of plants, trait–climate relationships suggest that traits are an important indicator of how plants adapt to different climates. In particular, analyses of plant functional traits in global databases suggest that there is an “economics spectrum” along which combinations of functional traits covary along a fast–slow continuum. We examine evolutionary associations among traits and between trait and climate variables on a strongly supported phylogeny in the iconic plant genus Protea to identify correlated evolution of functional traits and the climatic‐niches that species occupy. Results indicate that trait diversification in Protea has climate associations along two axes of variation: correlated evolution of plant size with temperature and leaf investment with rainfall. Evidence suggests that traits and climatic‐niches evolve in similar ways, although some of these associations are inconsistent with global patterns on a broader phylogenetic scale. When combined with previous experimental work suggesting that trait–climate associations are adaptive in Protea, the results presented here suggest that trait diversification in this radiation is adaptive.

Tracing the simultaneous evolution of individual traits and climate allows us to assess the role that adaptation to climate played in generating trait diversity. If adaptation is important in trait diversification, then evolutionary changes in functional phenotypic traits, those with presumed effects on survival, growth, and reproduction in the context of the abiotic environment (Violle et al., 2007), should be associated with changes in some aspect of the climatic-niche (although changes in biotic associations may also play a role). In plants, there is a wealth of interest in so-called economics spectra along a fast-slow continuum, such as the worldwide leaf economics spectrum (Wright et al., 2004), the wood economics spectrum (Chave et al., 2009), and the wholeplant spectrum (e.g., Edwards, Chatelet, Sack, & Donoghue, 2014;Reich, 2014). A recent global analysis identified two main dimensions of variation in plant traits: one along the leaf economics spectrum and another related to plant size (Díaz et al., 2016). Different suites of traits may be related to different aspects of the climatic-niche, and the global relationships that exist between traits and climate variables such as precipitation and temperature need to be tested within individual regions and lineages (Messier, Lechowicz, McGill, Violle, & Enquist, 2017;Moles et al., 2014;Wright et al., 2005). Trait-climate associations along the branches of a phylogeny suggest that adaptation plays a role in trait diversification. There is also substantial evidence for integrated trait evolution (phenotypic integration through evolutionary time) in plants (e.g., the worldwide leaf economics spectrum, Wright et al. 2004), so patterns of covariation within traits and within climate variables also need to be identified.
Trait-climate associations can be observed at many spatial and temporal scales. For example, statistical associations between fieldmeasured traits and environmental parameters provide evidence that contemporary trait differences are associated with important physiological and ecological functions both among and within distantly related genera in South Africa (Mitchell et al., 2015). By themselves, however, such associations do not provide evidence that differentiation in those functions played an important role in evolutionary diversification among species. Moreover, in field-based studies, observed differences could be the result of phenotypic plasticity rather than differential adaptation, a distinction that can be made only with common-garden studies (e.g., Givnish & Montgomery, 2014;Mason & Donovan, 2015). Evidence for the adaptive nature of trait-climate associations is strengthened by analyses incorporating evolutionary history. Specifically, assessing associations across a phylogeny requires the use of methods incorporating evolutionary relationships among taxa (such as phylogenetically independent contrasts, Felsenstein, 1985, or phylogenetic generalized least squares, Martins & Hansen, 1997). These approaches provide measures of trait-trait or trait-climate associations while controlling for phylogeny, that is, "evolutionary associations." We search for evolutionary trait-climate associations in a speciesrich and morphologically diverse plant lineage located in a biodiversity hotspot characterized by multiple climatic gradients. The genus Protea L. (Proteaceae) is a diverse group with 112 known evergreen plant species displaying diversity in growth form (ranging from subshrubs to shrubs and small trees), leaf shape and size, and inflorescence architecture (Rebelo, 2001;Valente et al., 2010). The age of the group is uncertain, but best estimates place the crown age at 5-18 my (Sauquet et al., 2009). The genus has its center of diversity in the Cape Floristic Region (CFR) of South Africa (Valente et al., 2010), a biodiversity hotspot characterized by high levels of species diversity (over 9,000 plant species) and endemism (about 70%, Goldblatt & Manning, 2002) in addition to being particularly threatened by human impacts (Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000). Protea is one of the dominant members of the fynbos community, although its range extends into northern and eastern portions of South Africa, Lesotho, Kenya, and central Africa (Rebelo, 2001;Rourke, 1980;Valente et al., 2010).
The extraordinary plant diversity in South Africa has been attributed to several different factors: the topographical complexity of multiple mountain ranges and "sky islands," sharp changes in soil types, soils that are extremely low in nutrients, steep gradients in temperature and in rainfall amount and seasonality, and the onset of the present-day climate dated at the Miocene-Pliocene boundary some 10 million years ago (Campbell, 1983;Linder, 2003;Verboom, Bergh, Haiden, Hoffmann, & Britton, 2015;Verboom et al., 2009). This diversity is largely accounted for by high species diversity in just 33 evolutionary radiations (Linder, 2003;Linder & Hardy, 2004;Schnitzler et al., 2011). Although an abundance of work has addressed the extent to which climatic and environmental heterogeneity has driven speciation and radiation throughout the region (e.g., Lamont, He, & Downes, 2013;Latimer, Silander, Rebelo, & Midgley, 2009;Linder, 2003;Verboom et al., 2009Verboom et al., , 2015, few systems have linked evidence from population genetics, common-garden experiments, and evolutionary analyses in this framework. Previous studies in Protea have documented contemporary associations between morphological traits and the environment across the genus (Mitchell et al., 2015). Experimental work within a smaller clade (the white proteas) and a single species (Protea repens) have found evidence for genetic differentiation of traits related to the environment consistent with adaptive differentiation in physiology and fitness (Carlson, Adams, & Holsinger, 2015;Carlson, Holsinger, & Prunier, 2011;Prunier, Holsinger, & Carlson, 2012). Here, we expand on these results and ask whether the diversification of traits is correlated with climatic factors across the genus in South Africa. More specifically, we ask: 1. Are there correlations between suites of fast-slow spectrum traits and species' climatic-niches across Protea? How do these compare with global patterns?

2.
Have traits and climatic-niche evolved in similar ways?
3. Are patterns of integrated evolution consistent with patterns of covariation observed in global datasets?
The answers to these questions allow us to identify suites of traits that have correlated evolution with suites of climatic-niche variables. In answering these questions, we also address how phylogenetic uncertainty and variation in both traits and climatic-niche values affect our conclusions.

| Incorporation of uncertainty
Existing methods for analyzing trait-by-trait or trait-by-climate associations at evolutionary timescales suffer from two important limitations: (1) They often assume that traits are uniform within species and (2) they usually ignore uncertainty in phylogenetic estimates (but see Huelsenbeck, Rannala, and Masly (2000)). These limitations may be especially important in radiations where diversification has occurred quickly, resulting in soft polytomies and uncertainty in species relationships. The use of Bayesian posterior tree samples or bootstrap replicates can account for some phylogenetic uncertainty (Huelsenbeck & Rannala, 2003), and very recent methods have begun to incorporate this uncertainty into estimates of correlated trait evolution (Caetano & Harmon, 2017), but the role of intraspecific trait variation has typically been neglected. The magnitude of intraspecific trait variation is often quite large, especially in some often-used plant functional traits (Auger & Shipley, 2013;Carlson et al., 2015;Donovan, Mason, Bowsher, Goolsby, & Ishibashi, 2014). This variation may have large impacts on comparative studies (Garamszegi & Møller, 2010), motivating the modeling of trait variances as well as means in evolutionary studies (Kostikova, Silvestro, Pearman, & Salamin, 2016).
These analyses, however, are difficult to carry out on large datasets.
We incorporate intraspecific variation in traits, climatic-niche, and phylogeny throughout our analyses.

| Trait measurements
We measured a suite of traits on plants from 58 different Protea species in the field from 2011-2013, including leaf and whole-plant traits associated with the leaf economics spectrum, as well as those found to be important in previous common-garden work (Carlson et al., 2011Prunier et al., 2012). We incorporated additional traits measured by Carlson et al. (2011Carlson et al. ( ) from 2008Carlson et al. ( -2009 in 133 species × site combinations that covered most of the range of Protea (Figure 1, average No. of observations per species = 26, range = 1-203). There were 1520 observations, of which <10 percent are complete, by design (See Table S1 for full data). For most populations (species × site combinations), we sampled eight plants for trait measurements, including height and canopy area (estimated from measured orthogonal dimensions of the plant using the formula for an ellipse) and sampled one of the most recently fully expanded leaves per plant. For shrub-like species, we also harvested wood samples from the previous year's growth for two plants per population. We measured leaf fresh weights and scanned leaves for analysis of length, width, and area in ImageJ (National Institutes of Health, Bethesda, MD, USA). Leaves were then dried and reweighed for dry weights. We did not include saturated weights in results presented here because preliminary analyses from 2011 found a correlation of 0.994 between saturated and fresh weights. For one to two leaves per population, we made stomatal peels on the adaxial side using clear nail varnish and tape that were later analyzed under a light microscope to estimate stomatal size and density. Four leaves per population were analyzed for carbon and nitrogen isotopes at the Stable Light Isotope Laboratory in the Archaeology Department at the University of Cape Town.
Wood density was estimated using a water displacement method as dry mass/wet volume (Cornelissen et al., 2003). We combined these F I G U R E 1 Map of individuals sampled for trait data for this study. Colors correspond to biomes as defined by Mucina and Rutherford (2006). Voucher and latitude/longitude data can be found in Table S1 -32.5 plant height (cm), plant canopy area (cm 2 ), leaf mass per area (lma, g/ cm 2 ; dry leaf mass/fresh leaf area), leaf freshwater content (fwc; [leaf fresh weight-leaf dry weight/leaf dry weight]), leaf length-to-width ratio (lwr; leaf length/leaf width), leafarea (cm 2 ), stomatal density (sd, stomates/cm 2 ), leaf nitrogen per unit mass (nmass; mg/g), leaf 13 C/ 12 C ratio (δ 13 C), leaf carbon-to-nitrogen ratio (cnratio), and wood density (wood, g/cm 3 ), Table 1. We natural-log-transformed all traits except for δ13C prior to analysis.
Instead of using species means to characterize traits in comparative analyses, we used Bayesian linear models to estimate distributions of trait values within each species, thereby obtaining a distribution that specifies the probability of particular values for each species and is consistent with the observed data. Specifically, we used the stan_glmer (or stan_glm if only one population) function in the R package "rstanarm" (Gabry & Goodrich, 2016) Tables S2 and S3, respectively. See Figure S1 for density plots of 100 randomly selected samples per species per phenotypic trait.

| Characterizing climatic-niches
We used Maximum Entropy Modeling (MaxEnt, Phillips and Dudik 2008) to characterize the climatic-niche for all 58 Protea species in our dataset. Latitude and longitude occurrence data for each species were extracted from the Protea Atlas database <http://www.
proteaatlas.org.za/>. The occurrence data included 94,715 individual georeferenced records across our species. Climate variables for each georeferenced point were extracted from the South African atlas of agrohydrology and climatology layers (Schulze et al., 2007) at the resolution of 1 by 1 min, or 1.55 by 1.85 km. Sites were grouped together if they were in the same grid cell, and all species observed in that grid cell were recorded. We retained ten variables that capture climatic gradients in the CFR and were used previously in the literature (Table 1): mean annual temperature (mat), average daily minimum temperature in July (t min ), average daily maximum temperature in January (t max ), elevation (elev), mean annual precipitation (map), interannual rainfall variability, measured as the coefficient of variation of mean annual rainfall across years (rflcv), temperature variability measured as maximum-minimum temperature at a site (t var ), mean annual potential evapotranspiration (pet), and two more direct measures of drought: the number of days receiving >2 mm of rain in the three driest months, hence lower values reflect more drought (rfl2 mm), and mean summer rainfall, December-February (summer_rain). Although many of these are highly correlated, we retain them to identify suites of climate variables that are associated with trait variation. All analyses were performed using MaxEnt in R through the package "dismo" (Hijmans, Phillips, Leathwick, & Elith, 2013). For each species in the dataset, we used 90% of the data for training and left 10% for testing.
Random pseudoabsences were taken from the extent of South Africa, because species of Protea are found in most South African biomes.
Details on the MaxEnt sampling and model results can be found in Table S4.
Similar to the trait values, we also incorporated uncertainty into estimates for climate niche variables to obtain mean values and a set of likely values for each species. We used the raw probabilities generated from the MaxEnt models to generate histograms for each species and square-root transformed climate variables (except for t min , which was left untransformed) using custom scripts from S.D. Smith (see Evans et al. 2009). These give a distribution of the predicted occupancy profile for each species in each climate variable independent of the other variables. From each of these distributions, we randomly sampled 10,000 observations and then calculated the 95% HPD intervals for each variable and each species in the R package "coda" (Plummer, Best, Cowles, & Vines, 2006). From these 9,500 samples, we calculated means per variable per species (Table S3) and also randomly selected 100 observations for use in downstream analyses (Table S2). Figure S1 has density plots of 100 randomly selected samples per species per niche trait.

| Accounting for uncertainty in phylogeny and traits
It is often difficult to estimate phylogenetic relationships in rapid radiations (Knowles & Chan, 2008). In earlier work, we used an anchored phylogenomics approach (Lemmon, Emme, & Lemmon, 2012) to sequence almost 500 nuclear genes conserved across all angiosperms (Buddenhagen et al., In Review)  We thus have a two-by-two

| Correlated evolution between traits and climate
We tested for correlations between traits and climate niche variables taking into account evolution using the program BayesTraits (Pagel & Meade, 2007) through R using the wrapper package "btw" (Griffin, 2015

| Evolutionary models
Comparisons of models of evolution can help to determine whether phylogenetic signal in traits is the result of phylogenetic niche con- similar to those of their ancestors (Peterson 1999). If, however, phylogenetic similarity in traits is not driven by climate, we may find a lack of niche conservatism.
To determine the evolutionary model that best fits each trait or climatic-niche variable, we used the fit continuous function in the R package "geiger" (Harmon, Weir, Brock, Glor, & Challenger, 2008) to compare the fit under Brownian motion (BM), Ornstein-Uhlenbeck (OU), and white noise models, setting the bounds for the alpha parameter in the OU model from 0 to 1,000. Model fit was evaluated using AICc scores and Akaike weights. We identify the best fit model using AICc scores using the "best" tree and mean trait/climate value, but also report results incorporating uncertainty across all four datasets (1 × 1, 1 × 100, 100 × 1, and 100 × 100).

| Correlated evolution among traits and among climate niche
A priori we expect correlations among traits and among climate variables in our dataset. We estimated coefficients taking evolution into account using BayesTraits as above for the 55 trait-trait and 45 climate-climate comparisons. To account for multiple comparisons, we used a Bayes factor cutoff of 9.8 for trait-trait comparisons and 9.6 for climate-climate comparisons. We built separate distance-based dendrograms for traits and climate variables using the correlation matrices from the 1 × 1 BayesTraits analysis to identify clusters for later analysis.

| Supplemental analyses
We performed two additional sets of analyses to ensure that our results are insensitive to important modeling choices. First, to assess the influence of log-transforming trait and climate data we performed all analyses on the 1 × 1 datasets using untransformed data and found results qualitatively similar to those we obtained using log-transformed data. Second, to assess the influence of the method of phylogenetic inference we used, we also performed all analyses on the 1 × 1 datasets using the best trees identified using SVDquartets (Chifman & Kubatko, 2014) and RAxML (Stamatakis, 2014) and found results qualitatively similar to those we obtained using ASTRAL-II.

| Trait-climate correlated evolution
We detected correlated evolution between morphological and climatic-niche traits in only a relatively small number of cases in two main clusters of strongly supported evolutionary associations. These include (1) plant size (leafarea, height, wood density) and its positive association with temperature (mat, t max , t min ) and negative relationships with rflcv and elev and (2) leaf composition (fwc, δ13C, lwr, nmass, cnratio, lma) with precipitation (summer_rain, map), where higher investment in leaves is found in drier areas (Figure 3). In addition, stomatal density has a strongly supported positive association with elevation, while wood density has associations with both temperature (positive) and rainfall (negative) variables. Estimated correlation coefficients for the BayesTraits analyses ranged from −0.513 to 0.627 in the 1 × 1 analyses, but only eight of 110 were individually very strongly supported. Nine additional correlations were strongly supported, 17 were weakly supported, and the remaining 76 lacked substantial support. The most strongly supported evolutionary correlations were between plant size (height) and variables related to elevation or temperature (mat, elevation, and t min ), where taller plant are found in warmer areas (Figure 3, Table S5).
Eight comparisons had a Bayes factor that exceeded the threshold of 9.7 to hold across all comparisons (Figure 3).

F I G U R E 3
Trait-climate associations for evolutionary analyses. Correlations are either positive (blue) or negative (magenta), vary in strength (size of circle), and have different levels of support indicated by transparency of circle color (weak support: logBF > 2, p < .10, most transpaent; moderate support: logBF > 5, p < .05, medium transparency; strong support: logBF > 10, p < .01, darkest circles). Correlations not supported at any level have no fill (completely transparent). Asterisks indicate correlations that meet the criterion threshold correcting for multiple comparisons. Dendrograms for the evolutionary analyses are based on distance matrices, and order is preserved in the contemporary data to more easily make visual comparisons

| Models of evolution
Morphological and climatic-niche traits appear to evolve in largely similar ways when analyses are based on the "best" tree and on trait or climate variable means for each species. All analyses of morphological traits for the 1 × 1 analysis (best tree and mean trait value) favored the OU model that indicates evolution around an optimum, or phylogenetic inertia in traits with limits on trait variance ( Table 2, Table   S6). Values for the alpha parameter are often quite high, although the OU model is still preferred over white noise (for which alpha is equal to infinity). Thus, the morphological traits of descendants are largely similar to those of their ancestors, and the range of morphological trait variation in the genus is somewhat constrained. Analyses of climaticniche traits also favored OU models for climate variables except for t var (which followed a white noise model) when analyses are based on the "best" tree and on species' climate means ( Table 2). Incorporating uncertainty (in the 1 × 100, 100 × 1, and 100 × 100 analyses) resulted in largely the same results for traits, although some replicates included support for BM or white noise models (particularly in canopy area, leafarea, and wood density; Figure S2, Table S6). For climate niche, incorporation of uncertainty in niche trait values resulted more often in white noise, especially with climatic-niche uncertainty, which may be an artifact of the MaxEnt based resampling scheme ( Figure S2).

| Correlated evolution along axes of variation
BayesTraits analyses reveal a wide range of pairwise correlations among morphological traits and among climatic variables (Figure 4).
We identified two major axes of morphological traits with strong patterns of covariation based on BayesTraits correlations: (1) general size (leafarea, height, canopy area) and (2) leaf composition (leafarea, cnratio, lma, δ13C, lwr, nmass). Leaf area is a component of both of these, while wood density also has weak to moderate connections with both axes. Stomatal density (sd) appears to evolve somewhat independently of the other traits ( Figure 4).
Eight of these met the Bayes factor threshold of 9.8 adjusting for multiple comparisons (Figure 4).
We identified three major suites of climatic variables with strong patterns of covariation: (1) temperature (mat, t min , and t max ), (2) rainfall (summer_rain, map, rfl2 mm), and (3) a group combining temperature, rainfall, and variability (pet, rflcv, elev, t var ; Figure 4). There is also an association between the temperature and rainfall axes of variation.
Estimated Bayes Traits pairwise correlations for climate variables ranged from −0.905 to 0.865 in the 1 × 1 analyses. Twenty-one of the 45 pairwise correlations were very strongly supported, four were strongly supported, six were weakly supported, and 14 lacked support ( Figure 4, Table S5). The strongest correlations were those between elevation and t min (corr = −0.905, logB = 94.5) and map and rfl2 mm (corr = 0.865, logBF = 78.8). Twenty-one of these met the Bayes factor threshold of 9.6 adjusting for multiple comparisons (Figure 4).
Incorporating trait uncertainty across all 210 pairwise BayesTraits associations (100 × 100 datasets) resulted in qualitatively similar outcomes (Figure 5a). The median values of the correlations for the 100 × 100 dataset were similar to the values from 1 × 1 analyses, although they tended to be less extreme. Notably, the upper and lower (97.5th and 2.5th percentiles) indicate that point estimates are very imprecise. In contrast, using other species trees has a minor influence on point estimates using other species trees, and the interval of the estimates is fairly narrow, indicating that phylogenetic uncertainty is not contributing heavily to differences in these estimates (Figure 5b).

| Correlated trait-climate evolution indicative of adaptive radiation
Our results provide several examples where morphological traits are correlated with climatic-niche variables when taking into account evolutionary history. Although these data are collected in the field and patterns could be due to phenotypic plasticity, we find similar traitclimate relationships in this and previous studies within a subset of the region (Mitchell et al., 2015) and in experimental common-garden studies (Carlson et al., 2011;Prunier et al., 2012). The overall patterns found here involve the same two dimensions of variation in plant traits found in the global analysis by Díaz et al. (2016): plant size and leaf economics. Here, we find a positive association between plant and leaf size and temperature, and a negative association between leaf construction cost (e.g., LMA) and rainfall amount.
The finding of correlated trait-climate evolution along two main axes is consistent with previous work in Protea, suggesting that traitclimate associations are adaptive. For example, Carlson et al. (2011) showed in the white proteas (a strongly supported clade of six species: Mitchell et al., 2017) that trait-performance differences in two experimental gardens are consistent with patterns expected from trait-climate associations in wild populations and that in-garden survival is associated with trait differences in the predicted way. That work provided direct evidence that many trait-climate associations are adaptive. Similarly, selection gradient analyses using an estimate of lifetime seed production as a proxy for fitness demonstrated differential selection in two wild populations consistent with predictions of trait-climate associations in P. repens .
Although the results here do not provide direct evidence for adaptation, they suggest that adaptive differences identified within species and among closely related species may be extrapolated across the entire genus.
Plants commonly tend to be taller and have larger leaves with dense wood in warm areas and to be shorter with smaller leaves and less dense wood at high elevations (see Kdimer, 1999), reflecting both hydraulic and allometric constraints to woody plant height (Givnish, Wong, Stuart-Williams, Holloway-Phillips, & Farquhar, 2014 with soil depth, where lowland areas generally have deeper soil than high elevation sites (Campbell and Werger 1988). Our associations between plant size and temperature are consistent with the finding that mean annual temperature is generally a better predictor of plant traits than mean annual precipitation (Moles et al., 2014). In contrast, Moles et al., 2009 found that precipitation during the wettest month was most strongly associated with plant height, while Thuiller, Lavorel, Midgley, Lavergne, and Rebelo (2004) found no association between climate and plant height and found that leaf size was correlated with precipitation, not temperature, in the related genus Leucadendron. Yates, Anthony Verboom, Rebelo, and Cramer (2010) concluded that leaf size in Proteaceae may reflect the ability to dissipate heat during warm periods and quickly transpire during wet periods, although our patterns of leaf size and temperature are opposite to their findings.
The other broad pattern that we detected is that "slow" trait values on the leaf economics spectrum are associated with lower levels of rainfall. Plants in areas with low rainfall invest more in their physical construction, with higher values of LMA and wood density, higher C:N ratios, lower nmass, and lower leaf water contents. These results are consistent with the hypothesis that higher values of LMA are associated with optimizing plant growth in drier or nutrient-poor areas (Givnish, 1979;Reich et al. 2003, Wright et al., 2005, but they may be inconsistent with previous results, suggesting that temperature is more strongly associated with higher carbon investment in leaf construction (Moles et al., 2014). Lamont, Groom, and Cowling (2002) tested these associations more explicitly using Proteaceae species in both South Africa and Australia along precipitation gradients. They found that LMA was inversely correlated with rainfall and that this relationship is perhaps more important than soil nutrient status.
Wood density is a component of both suites of traits, and its positive association with temperature and negative association with precipitation are consistent with global patterns (Swenson & Enquist, 2007) and may reflect an adaptation to decrease the chance of xylem cavitation during periods of drought stress and the lack of water availability associated with hotter and drier downslope regions in the CFR (Hacke, Sperry, Pockman, Davis, & McCulloh, 2001).
Similar associations between the evolution of morphological traits and shifts in climatic-niche have also been seen in Pelargonium in the CFR (Jones et al., 2013), but the trait-climate associations found in Protea and Pelargonium can differ (Mitchell et al., 2015). For example, larger plants are associated with more drought in Pelargonium, while the opposite is true in Protea. These differences may be due to the occupation of unique microhabitats or differing growth forms and life history strategies. Similarly, the lack of detectable associations between leaf mass per area and temperature variables is not at odds with the weak negative association in evergreen plants on the global scale, while the negative relationship with summer rain here is consistent with the global relationships across evergreens (Wright et al., 2005).
Precipitation may be a more important driver than temperature for LMA in evergreen plants in the CFR, where rainfall amount is more variable and results in periods of drought stress (as opposed to tropical regions).

| Integrated trait evolution
Analysis of both contemporary and evolutionary trait associations provides evidence for covariation of morphological traits and climate variables in Protea. Individual morphological traits and climatic-niche variables are used to indicate significant aspects of the "whole organism" and "n-dimensional hypervolume" climaticniche variation (Hutchinson, 1957;Reich, 2014;Reich et al., 2003).
Correlations among morphological traits are one way of measuring the degree to which phenotypes are integrated, whether as a result of shared function, development, or genetics (Pigliucci, 2003). Patterns of integration may be especially interesting when analyzing traits across organ systems or those involved with different physiological processes in plants. For instance, hydraulic capacity in plant stems can be linked to both hydraulic and photosynthetic rates in leaves (Brodribb, Holbrook, & Gutiérrez, 2002;Sack, Cowan, Jaikumar, & Holbrook, 2003). Similarities or differences between contemporary and evolutionary associations between morphological traits may provide clues to the causes of phenotypic integration. For example, the strong negative associations between height and canopy area in both sets of analyses could reflect fundamental biophysical constraints. In contrast, the negative association between lma and δ 13 C (sclerophyllous leaves with higher lma have more negative δ 13 C values (less water use efficient)) may reflect a functional association between strategies enhancing resource conservation and those enhancing water use efficiency.  Reich, 2014;Wright et al., 2004). Each of these syndromes involves suites of correlated traits, combined with trade-offs among them.

Some patterns of trait integration in
In Protea, for example, we find a positive evolutionary association between cnratio and lma and a negative association between nmass and lma, consistent with worldwide patterns in the LES. Similarly, we find (1) a positive evolutionary association between leafarea and plant height and then between plant height and canopy area, and (2) a positive association between leafarea and cnratio and a corresponding negative association between leafarea and nmass, suggesting that there are integrated traits related to overall investment and plant size (Figure 4a). At the evolutionary scale, wood density is weakly or moderately linked to attributes of both axes (plant size and leaf investment), while stomatal density evolves largely independently. Wood density may be bridging the gap between the suites of trait, as it can be tied to plant size and physical support associated with temperature, and the need for cavitation resistance associated with low rainfall.
Not surprisingly, several climatic-niche variables also covary. As the climatic-niches descendants occupy are, at best, weakly correlated with those of their progenitors, these associations probably reflect intrinsic physical climate correlations rather than correlated niche evolution. For example, the positive network of mat, t max , and t min , and their negative associations with elevation are expected due to adiabatic cooling and would likely be detected in any random sample of geographic locations. In contrast, the associations between seasonality, precipitation, and temperature reflect niche hypervolumes characteristic of the CFR and perhaps of other Mediterranean climate regions around the world characterized by cool, wet winters and hot, dry summers.

| Incorporating uncertainty in comparative analyses
Two sources of uncertainty should be recognized in any comparative analysis: (1) uncertainty about trait values for the species that arises because of trait variation within species and (2) uncertainty about species relationships that arises because phylogenetic relationships are imperfectly estimated. In our case, incorporating these sources of uncertainty did not qualitatively affect our results in the BayesTraits analyses. The median values incorporating both intraspecific trait variation and phylogenetic uncertainty (i.e., the 100 × 100 analyses) were very similar to those from the 1 × 1 analyses based on the species mean value and the "best" phylogenetic tree. Nonetheless, the range of values generated in the 100 × 100 analyses indicates that results using any single point estimate, such as a mean, must be interpreted with considerable caution. This is particularly evident in the analysis of evolutionary models for climatic-niche traits, where an Ornstein-Uhlenbeck model was supported in the 1 × 1 analysis for all but t var , while white noise models were often favored in the 100 × 100 analyses, although this may be due to the nonnormal sampling from the MaxEnt output.

| CONCLUSIONS
The correlated evolution of traits and climatic-niches in Protea, combined with previous experimental work in the group, suggests that broadly adaptive processes played an important role in divergence of many morphological traits during this radiation. Moreover, our results provide evidence for two different suites of trait-climate associations, possibly reflecting a combination of fundamental biophysical or functional relationships among traits and intrinsic physical properties of climate variables. In particular, we identified that the same two dimensions of plant variation found at the global scale are associated with different aspects of the climate, namely: (1) a positive association between plant size and temperature, and (2) a "faster" leaf investment strategy with higher amounts of rainfall. These are supported by consistent findings between evolutionary and contemporary associations.
Phenotypic traits are relatively conserved, which may be accounted for by climatic-niche conservatism. Overall, we find substantial evidence for broadly adaptive coevolution among traits and climate even when incorporating both uncertainty in phylogenetic relationships and to within-species variation in trait values. Together with past work demonstrating that these traits are adaptive and heritable, we find evidence that the radiation of Protea is indeed adaptive. Future work on trait and physiological differentiation in closely related cooccurring species of Protea may provide more robust evidence for the mechanisms underlying these phenotype-climate associations.