Polyploidy promotes divergent evolution across the leaf economics spectrum and plant edaphic niche in the Dianthus broteri complex

Journal of Ecology. 2022;110:605–618. | 605wileyonlinelibrary.com/journal/jec Received: 6 September 2021 | Accepted: 22 November 2021 DOI: 10.1111/1365-2745.13823 RESEARCH ARTICLE Polyploidy promotes divergent evolution across the leaf economics spectrum and plant edaphic niche in the Dianthus broteri complex Javier López- Jurado | Enrique Mateos- Naranjo | Francisco Balao This is an open access article under the terms of the Creat ive Commo ns Attri butio n- NonCo mmerc ial- NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is noncommercial and no modifications or adaptations are made. © 2021 The Authors. Journal of Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society. Departamento de Biología Vegetal y Ecología, Facultad de Biología, Universidad de Sevilla, Sevilla, Spain Correspondence Javier López- Jurado Email: [email protected] Funding information Áridos La Melera S.L., Grant/Award Number: 2234/0724; Fifth Research Plan from the University of Seville, Grant/Award Number: USE- 19606- X; Sixth Research and Transfer Plan from the University of Seville, Grant/Award Number: USE- 21358- Z; Ministerio de Ciencia e Innovación, Grant/Award Number: PGC2018- 098358- B- I00 Handling Editor: Stephen Bonser Abstract 1. The evolution of the leaf economics spectrum (LES) is known to be constrained by genetic relatedness but also promoted at small geographical and phylogenetic scales. In those cases, we hypothesized that polyploidy would play a prominent role as an outstanding source of functional divergence and adaptive potential. 2. We registered leaflevel nutrient, water- and lighteconomy- related traits from the LES as well as edaphic properties in the four cytotypes of the autopolyploid Dianthus broteri complex (2×, 4×, 6× and 12×). We analysed the effect of ploidy level on the integration of the LES network, checked if concerted evolution occurred between LES and soil niche and tested the influence of phylogeny on the variables. Alternative evolutionary models for both sets of traits were compared. 3. We found higher divergence of polyploids (especially 6× and 12×) compared to diploids in the LES and soil niche, but these traits are not coevolving. 6× and 12× showed opposite ecological strategies regarding resource use and higher uncoupling of the LES network. Early divergence of traits prevailed in both LES and edaphic niche (supported by better fitted evolutionary models with one optimum per cytotype), but postpolyploidization processes played an important role for the photochemical behaviour. 4. Synthesis. Our results indicated shifts in ecological strategies across Dianthus broteri cytotypes and suggested a powerful role of polyploidy in overcoming constraints for the evolution of plant functional traits. KEYWORDS adaptive divergence, ecological strategies, functional traits, leaf economics spectrum, polyploidy, soil niche 1 | INTRODUCTION Plant functional traits are those impacting their demographic fitness, that is, fecundity, growth, development or survival (Violle et al., 2007; Volaire et al., 2020). Notably, major functional traits in vascular plants are concentrated and coordinated in a multidimensional space explaining plant form and function at a global scale (Díaz et al., 2016). Within the global plant trait variation, there is a fully differentiated axis corresponding to the leaf economics spectrum (LES; Wright et al., 2004), which accounts for the return on investments in 606 | Journal of Ecology LÓPEZJURADO Et AL. leaves regarding carbon and nutrient acquisition and processing. The coordination among LES traits, especially on a mass basis, causes this spectrum to be defined by a tradeoff between plant longevity and growth potential (Wright et al., 2004). According to this, ecological strategies would be reflected by the positioning of a given plant species along the LES (Reich, 2014), varying from nutrientacquisitive strategies, also known as fast (with high photosynthetic capacity and shortlived leaves), to conservative ones (i.e. slow, with low photosynthetic capacity and longlived leaves). Beyond nutrientassociated traits, this spectrum seems to be coordinated with the essential resources of water and light (Reich, 2014). Previous studies have already demonstrated the correlation between plant nutrient and water economy (e.g. Pérez- Ramos et al., 2012; Prieto et al., 2018), being the species with slow strategies also conservative about water and therefore more tolerant to drought (Hallik et al., 2009). Furthermore, the coordination of light economy with the purely LES traits has also been indirectly addressed by linking shade tolerance with less thick and dense leaf blades, probably to achieve larger harvesting surfaces in lightlimited environments (e.g. Hallik et al., 2009; Lusk et al., 2008). On the basis of these findings and the results obtained by Hallik et al. (2012) about chlorophyll fluorescence characteristics in woody and herbaceous species under different growth irradiance conditions, we hypothesized that light collection at the photosystem level would be enhanced for less lightexposed leaves while PSII photochemical efficiency would be diminished as a result of photoinhibition. In any case, the joint study of variables related to nutrient, water and light use are needed to delve into the interdependence of functional traits that assure plant success in a given environment, following the respective recommendations of Poorter et al. (2014) and Reich (2014). However, globalscale trait covariation in the LES is challenged by nonglobal scale predictions (Mason et al., 2016). Locally, natural selection can generate trait correlations that are relatively easy to escape (soft constraint) by drivers of phenotypic variation, altering the LES tradeoff. These drivers include plasticity, ecological filtering, sexual reproduction, adaptation, genetic drift or ontogenetic mechanisms (Messier et al., 2017) and may lead to a diversity of functional trait associations within phylogenetically closely related groups (Carvalho et al., 2020; Niinemets, 2015) and even at the intraindividual level (between leaf cohorts; Puglielli, 2019; Puglielli & Varone, 2018). Noticeably, a significant relationship occurs between environmental conditions and LES traits, as the influence of soil properties (Lambers & Poorter, 1992; Lambers et al., 2008), but also the interaction between soil and climate (Luo et al., 2015, 2016) has been addressed on plant functional traits and resourceuse strategies. This supports the interpretation of trait– environment associations as a result of adaptive evolution (Reich et al., 2003), highlighting the role of resource availability or climatic conditions in functional divergence along the LES (Brouillette et al., 2014; Sartori et al., 2019). Apart from natural selection, harder constraints (i.e. those that are almost impossible to escape) to novel trait associations shape the expression of the LES. These are of biophysical and genetic nature (Messier et al., 2017), being the former the most restrictive ones and consistent within and across ecological scales. Among biophysical constraints, the link between photosynthetic capacity and leaf life span is of paramount importance since there are unavoidable limitations to construct longlasting leaves with high productivity. Furthermore, the evolution of the LES is expected to be limited to some extent by a genetic component, although this topic deserves further investigation (Donovan et al., 2011). Specifically, the referenced study suggested that both selection and genetic constraints need to be measured jointly to test their relative roles on the LES evolution. To date, it is known that a few pleiotropic genes maintain the strong correlations between LES traits (Vasseur et al., 2012), which would have a major role in the phenotypic integration within this complex network (Ji et al., 2020). In the light of this evidence, the release of genetic constraints would be the baseline for promoting functional divergence in the LES. It would cause evolutionary novelties such as adaptive phenotypic plasticity or epigenetic variation (Kalisz & Kramer, 2008). Thus, polyploidy (i.e. the state of having more than two paired sets of chromosomes) may play a pivotal role on LES evolution over low phylogenetic scales given that it is considered an outstanding driver of these sources of adaptive modifications (Chen, 2007; Hahn et al., 2012). Allopolyploids (those originated from interspecific hybridization) are prone to produce numerous combinations of traits due to their diverse origin (Seehausen, 2004), but genome duplication itself and subsequent adaptation have also demonstrated the capacity to cause phenotypic changes in autopolyploids (those originated from a single species; Martínez et al., 2018; Visger et al., 2016; Wei et al., 2020). Moreover, autopolyploidy per se can lead to novelty by modifying the phenotypic covariance matrix, meaning the disintegration or uncoupling of those traits (Oswald & Nuismer, 2011). This evolutionary innovation in polyploids can be caused either by relaxed selection on duplicate genes due to their higher functional redundancy (Flagel & Wendel, 2009) or increased variation in dosageregulated gene expression levels (Osborn et al., 2003). Using a recently radiated autopolyploid complex, we can not only assess the effect of polyploidization and rapid adaptive processes on functional divergence, but we also would provide new insights into the neglected ecological context of noncrop polyploid systems (Soltis et al., 2016). An approach like the one used in this work would facilitate underlying explanations to several widely studied ecological processes in polyploids, such as competitive interactions (see e.g. Čertner et al., 2019), fitness advantage (see e.g. Wei et al., 2019) or niche differentiation (see e.g. Rojas- Andrés et al., 2020). Considering this background, we studied LES variation across the autopolyploid complex Dianthus broteri (Caryophyllaceae). This group of perennial herbaceous plants consists of four different ploidy levels (2×, 4×, 6× and 12×; Balao et al., 2009) that diverged about 0.9– 2.1 Mya (Balao et al., 2010) and have different evolutionary histories (the 4× cytotype has two independent origins while 6× and 12× have single origins; Balao et al., 2010). Moreover, D. broteri has a Mediterranean distribution (exclusively in the Iberian | 607 Journal of Ecology LÓPEZJURADO Et AL. Peninsula), where its ploidies do not coexist and occupy disjunct geographical ranges with markedly different environmental conditions. Comparatively, diploids and tetraploids inhabit temperate and broad niches with nutrientrich soils, whereas hexaploids and dodecaploids are distributed in warm and arid habitats with infertile soils (López- Jurado et al., 2019). The intercytotype differentiation pattern in this complex also encompasses global methylation levels (Alonso et al., 2016) and photochemical responses (López- Jurado et al., 2020). Interestingly, morphological features (vegetative and reproductive organs) not only showed divergence among D. broteri cytotypes, but also a loss of phenotypic integration with increasing ploidy (i.e. lower covariation among functionally related traits; Balao et al., 2011). The shared conclusions from the mentioned studies suggest adaptations of polyploids (particularly 6× and 12×) to specific biotic and abiotic factors. Besides that, a focus on the 12× has been established caused by unique properties as its floral scent emission (Balao et al., 2011), high tolerance to recurrent drought (López- Jurado et al., 2016) and alleged key role in ecosystem functioning (López- Jurado et al., 2019). Based on the mentioned previous results obtained for this complex, our starting hypotheses were that (a) D. broteri has a variety of LES trait combinations and hence ecological strategies, which would be acquisitive in loworder ploidy levels (2× and 4×) and conservative in highorder ones (6× and 12×; in response to the stressful niches that they inhabit) and (b) polyploidy is enhancing the uncoupling of functional traits in this complex, especially in the 6× and 12× cytotypes. In this study, we aimed to uncover polyploidymediated shifts in the LES and the edaphic niche, as well as resource (nutrient, water and light) acquisition and resourceuse strategies that enable the D. broteri polyploid complex to occupy heterogeneous environments. 2 | MATERIALS AND METHODS 2.1 | Characterization of study sites During July 2017, three natural populations of each Dianthus broteri cytotype encompassing their corresponding geographical distributions were sampled (between 10:00 and 14:00 hr; Figure S1). All these locations were characterized by considerably low D. broteri density, ranging between 30 and 100 individuals, which were accompanied mostly by xerophytic Mediterranean shrubland. In each location, the main climatic conditions were recorded (Table S1), besides performing measurement and sampling procedures described in the following sections. As previously shown in López- Jurado et al. (2019), populations of lower cytotypes (2× and 4×) had milder climatic conditions compared to those of higher ones (6× and 12×), and showed intermediate values in relative humidity and the lowest temperatures. On average, the niche of hexaploids was the most humid during sampling time (summer) and dodecaploids were exposed to the warmest and driest conditions (Table S1). 2.2 | Soil composition and properties In each population, soil samples (n = 3) from the upper 15 cm were collected throughout the area occupied by D. broteri individuals, placed in plastic containers and conserved at 4°C until their processing in the laboratory. Soil granulometric composition was determined and chemical characterization considered, for example, soil electrical conductivity (EC), which is an indirect indicator of the amount of water and watersoluble nutrients available for plant uptake (Smith & Doran, 1997). Organic matter content (OM; by the Walkley– Black method), total carbonates (CaCO3), active limestone and nitrogen (N, by the Kjeldahl method) were also included in the analysis. Finally, we quantified assimilable phosphorus (P) by ultraviolet– visible spectrophotometry, as well as exchangeable potassium (K), calcium (Ca), magnesium (Mg), sodium (Na), iron (Fe), manganese (Mn), zinc (Zn) and copper (Cu) by atomic absorption spectroscopy (Fernandez- Hernandez et al., 2010; Usman et al., 2017). 2.3 | Leaf nutrient concentrations and isotopic composition Completely developed mature distal leaves (n = 5 per individual, n = 50 per population) were collected from the middle canopy layer of plants. Leaves from the same population were blended, dried at 60°C for 48 hr and finely ground using a mixer mill (Retsch GmbH— MM 301). Carbon and nitrogen contents of each leaf pool replicate (n = 3 per population) were quantified using a LECO TruSpec CN elemental analyzer. Additionally, Inductively Coupled Plasma Optical Emission Spectrometry (ICPOES) analysis was carried out in a Thermo Scientific iCAP 6500 spectrometer to determine levels of macro- and micronutrients in leaf pools (n = 3 per population; Table S3). Among the numerous leaf nutrients, we focused on N, P and K levels because of their tight association with plant productivity (Niinemets & Kull, 2005; Sardans & Peñuelas, 2015). N and P influence photosynthesis by having structural and functional roles in the Calvin cycle and the thylakoids (Evans, 1989; Vance et al., 2003) and K is involved in a broad range of plant physiological functions (e.g. enzyme activation and regulation of osmotic potential or stomatal opening; Wang et al., 2013). The calculation of C and N isotopic signature on dry and pulverized leaf samples (n = 3 per population) was conducted using a Thermo Scientific Delta V spectrometer. The carbon and nitrogen isotope ratios (13C/12C and 15N/14N) were expressed in delta (δ) notation, defined as the parts per thousand (‰) relative to a standard material (VPDB). The analytical precision for the measurements was 0.2‰. δ15N was used as a proxy for nitrogen fluxes, assimilation and allocation (Kalcsits et al., 2014) supplementary to leaf N content. Furthermore, wateruse strategies were approached by δ13C, which was employed as a proxy for iWUE (intrinsic wateruse efficiency; Pérez- Harguindeguy et al., 2013) because it relies on the ratio of 608 | Journal of Ecology LÓPEZJURADO Et AL. intercellular to ambient CO2 concentrations (Ci/Ca) and therefore is affected by transpiration at the stomata (gs) and net assimilation (AN; Farquhar et al., 1989). 2.4 | Leaf traits and physiological measurements During field sampling, we collected 10 completely developed mature leaves in each population (1 leaf per individual) and measured their areas in the field by digital image analysis using ‘Easy Leaf Area’ software (Easlon & Bloom, 2014). Dry weight (DW) was determined after drying the leaves at 60°C until constant weight was reached and leaf mass per area (LMA) was then calculated as the ratio between DW and leaf area (g DW/m2). LMA was selected as a crucial variable to account for leaf constructions costs as well as plant growth rate (low LMA indicates thin leaf blade or not very dense tissue, or both, which contribute to short leaf life spans and fast growth rates; Lambers & Poorter, 1992; Westoby et al., 2002; Wright et al., 2004), also being highly variable within functional groups and responsive to environmental conditions (Poorter et al., 2009). Furthermore, the maximum quantum efficiency of PSII photochemistry (Fv/Fm) and the chlorophyll a fast kinetics were registered in 30 min darkadapted, completely developed leaves (n = 25 per protocol) from randomly selected individuals of each population, using FMS- 2 (Hansatech Instruments Ltd., UK) and FluorPen FP100 (Photo System Instruments, Czech Republic) portable pulsemodulated fluorimeters, respectively. For the chlorophyll a transient data, we used the preprogrammed OJIP test protocols of FluorPen FP100 and calculated derived parameters following Strasser et al. (2010). In particular, we analysed the informative energy transduction fluxes on a leaf crosssection basis (absorbed, ABS/CS, trapped, TR/CS, electron transport, ET/CS, and dissipated, DI/CS, energy fluxes) since significant differences in the density of reaction centres (RC/CS) among cytotypes were previously found (oneway analysis of variance, ANOVA: F3,246 = 5.91; p < 0.05). Fv/Fm and OJIP testderived total performance index (PITotal) were considered based on the broad view they offer about the overall photochemical activity and their sensitivity to stressful conditions, as those of the sampling time in summer. 2.5 | Edaphic divergence of polyploids After obtaining a detailed description of the edaphic niche characteristics in the 12 sampled D. broteri populations, we estimated their phylogenetic relationships from a previously performed phylogram based on amplified fragment length polymorphism (AFLP) data (see Balao et al., 2010, 2011). Polytomies were resolved in random order using the multi2di function from the ‘ape’ package in r (Paradis et al., 2004). Using the ‘adephylo’ r package (Jombart, Balloux, et al., 2010), we then performed a phylogenetic principal component analysis (pPCA; Jombart, Pavoine, et al., 2010), which can detect sets of soil traits presenting a similar phylogenetic signal and hence informs about the general influence of phylogeny on the measured variables. Through this method, we attempted to identify sources of intercytotype variation due to phylogenetic relatedness as well as local divergence (i.e. adaptation or phenotypic plasticity). On the one hand, axes with the most positive eigenvalues reflected positive phylogenetic autocorrelation, corresponding to early diverged traits (global structures sensu Jombart, Pavoine, et al., 2010). At the microevolutionary scale studied in the D. broteri complex, early divergence is tightly related to polyploid formation (Balao et al., 2010), although it is not completely possible to separate direct effects due to polyploidization from rapid selection promoting the polyploid establishment without experimenting with synthetic polyploids. On the other hand, axes with the most negative eigenvalues mirrored the most recently diverged traits (lower phylogenetic autocorrelation than null expected or local structures), related to postpolyploidization processes such as adaptation and phenotypic plasticity. 2.6 | Divergence and integration in the LES covariation network For the analyses involving LES and its relationship with polyploidy, we used the previously introduced meaningful variables to describe ecological strategies from nutrient (LMA, leaf N, P, K and δ15N) to water use (δ13C), but also the photochemical behaviour (ABS/CS, TR/CS, ET/CS, DI/CS, Fv/Fm and PITotal). Thus, we carried out a pPCA, as previously described, using the 12 selected functional variables to investigate the effects of phylogeny and downstream selection on the LES. The influence of ploidy level on both sets of variables (functional and edaphic) was analysed using oneway ANOVAs in the two most explanatory axes from each pPCA. In case of significant results, posthoc Tukey's HSD tests (α = 0.05) were performed to identify intercytotype pairwise differences. The overall covariation network in LES traits considering the influence of phylogeny was analysed using the ‘corrr’ package in r (Kuhn et al., 2020). We also quantified the magnitude of integration in the LES network, as a measure of coordination between functional traits, for each D. broteri cytotype separately based on the index of phenotypic integration (INT; Cheverud et al., 1989; Wagner, 1984), following the approach of Damián et al. (2020). We first constructed one correlation matrix per cytotype using Pearson's r and compared them by Mantel tests. We calculated INT values by accounting for the variance among eigenvalues (λi) of each correlation matrix and then corrected the values as INTc = (Var (λi) − (number of traits − 1)/ number of individuals per cytotype). We transformed each INTc to percentage of the maximum integration, which equals the total number of traits in the correlation matrix. Finally, we obtained the 95% confidence intervals for each INTc by the bootstrapping method (R = 10,000). High INTc values were indicative of high variance among eigenvalues and therefore strong phenotypic integration (i.e. high trait coordination), whereas low INTc values meant that phenotypic integration was weak in a particular cytotype (i.e. low | 609 Journal of Ecology LÓPEZJURADO Et AL. trait coordination). Less overlapping confidence intervals among cytotypes were considered as evidence of variation in phenotypic integration. 2.7 | LES and edaphic niche evolution: Divergence in nutrient uptake and use efficiency Complementary to both pPCAs (i.e. soil and LES), we tested for phylogenetic signal on the most explanatory axes with global structures using the Blomberg's K (Blomberg et al., 2003) in the ‘phytools’ r package (Revell, 2012). Given the cytotype relatedness pattern, the global structures can be interpreted as a consequence of niche/functional conservatism (following a Brownian motion model of gradual drift; BM) or divergence due to polyploidization (following an Ornstein– Uhlenbeck, OU, model of stabilizing selection with different optima; Butler & King, 2004). To investigate these alternatives, we compared the relative fit of several evolutionary models in case of significant phylogenetic signal: BM as well as an OU model with a single optimum (OU1) and an OU model with four optima (corresponding to one optimum per ploidy; OU4). Models were performed with the ‘mvMORPH’ r package (Clavel et al., 2015) and compared by the weights of their Akaike information criterion values (AIC). In case of pPCA axes with local structures, we checked them by fitting a null model of trait evolution (NP) through the simulation of a starshaped phylogenetic tree. Furthermore, to investigate the degree to which edaphic niche characteristics predicts LES variables (and vice versa; i.e. correlated evolution), we employed procrustean randomization tests (Peres- Neto & Jackson, 2001) through the functions procrustes and protest in the ‘vegan’ r package. We compared the pPCA scores obtained for each population in the two axes that explained more variation (PC1 and PC2 for the soilpPCA; PC1 and PC11 for the LES- pPCA), using 9,999 permutations. A significant high correlation (p < 0.05) between PCs indicated interdependence of datasets and correlated evolution of characteristics. Finally, we also investigated a wide variety of element contents in leaves. Macronutrient uptake and use efficiency in the four D. broteri cytotypes were determined, for each location, by accounting for the differences between leaves and soils in the concentration of N, K and P (in mg/g of dry matter). The influence of ploidy level on these latter variables was analysed using generalized linear models (oneway ANOVAs) with normal or gamma distribution of errors. Significant results were followed by posthoc Tukey's HSD tests (α = 0.05) to identify intercytotype pairwise differences. All statistical analyses were performed in R software ver. 4.0.0 (R Core Team, 2020). 3 | RESULTS D. broteri populations showed a pattern of intercytotype differentiation in the edaphic space, as represented in the two most explanatory pPCA axes (soil- PC1 and PC2). Both axes explained more than half of the total variance in the edaphic space (37.7% and 20.1%, respectively) and showed positive eigenvalues, so local structures (related to recently diverged traits) had little influence on the variation of soil niche characteristics. Moreover, variable contributions on these axes were generally low (always <0.5 in absolute values), revealing high levels of correlation between edaphic variables. Despite this result, the first axis (soil- PC1) clearly mirrored an alkalinity gradient and the second one (soil- PC2) a nutrient content gradient in soils. In this edaphic space, we noticed niche displacement and expansion of the polyploid cytotypes, being 4×, 6× and 12× niches clearly different and wider than the 2× one. Thus, diploids had the most restricted niche, which almost completely overlapped the tetraploid space (Figure 1B). Congruently, the two lower cytotypes (2× and 4×) did not present significant differences in any pPCA axis (Figure 1C,D) and they were associated with milder edaphic conditions. Their niches were characterized by clay, semiacid soils with high EC as well as elevated contents of essential macronutrients (N, P and Ca) and several important micronutrients (Cu and Mn; Figure 1A,B). Oppositely, the 6× cytotype was located in a significantly divergent position of soil- PC1 compared to the rest of ploidies (F3,8 = 6.18; p < 0.05; Figure 1C). In this axis, the niche space of hexaploids corresponded to silty and alkaline soils, with pH ≥ 7 and high levels of total carbonates and active limestone. Lastly, the 12× cytotype appeared in an extreme position of soil- PC2, marginally different from the other cytotypes (F3,8 = 2.37; p = 0.07; Figure 1D). The niche space of dodecaploids in soil- PC2 was related to nutrientpoor soils (except for Fe), with reduced concentrations of all the essential macronutrients (N, P, K, Ca and Mg), organic matter and electrical conductivity. D. broteri also demonstrated a notable variation of the LES traits, which was mainly explained by two pPCA axes accounting for 47.7% and 25.0% of the total inertia, respectively (Figure 2A). The first most explanatory axis had a positive eigenvalue but the second one had a negative eigenvalue, indicating that both global and local structures affect to LES covariation in D. broteri. Therefore, in the first axis (LESPC1), we found the earliest diverged traits, which were those from the classical LES. LESPC1 reflected intercytotype divergence, being the 12× significantly differentiated from the rest of cytotypes (F3,8 = 16.78; p < 0.05; Figure 2C). The dodecaploids were characterized in this axis by high levels of leaf macronutrients (N, P and K) and high PITotal as well as low LMA and δ13C (Figure 2A,B). In this same axis, the other three ploidies overlapped and appeared in more central positions of the pPCA, with LES traits showed by the 2× cytotype encompassing both 4× and 6× (Figure 2B). Nonetheless, hexaploids were occupying the most opposite location in PC1 compared to dodecaploids. The second represented axis (LESPC11) gathered recently diverged traits and was mainly characterized by photochemical parameters (TR/CS, ET/CS, DI/CS and Fv/Fm). The scores of LESPC11 did not differ among cytotypes (F3,8 = 0.19; p > 0.05), being this axis mostly defined by intracytotype variation, specifically within the 4× (Figure 2D). Interestingly, the photosystems efficiency (Fv/Fm) 610 | Journal of Ecology LÓPEZJURADO Et AL. and performance (PITotal) showed opposite patterns for all the ploidies considering their reverse contribution signs in both axes (Figure 2A). Overall, the welldescribed covariation relationships for functional traits from the LES were maintained in D. broteri. A strong positive relationship was found between the essential macronutrients N and P (r = 0.72; p < 0.05), which were negatively related to δ13C (Figure 3). This isotopic composition was, in turn, negatively related to δ15N but positively to LMA. However, photochemical variables were separated from the classical LES traits in this covariation network but integrated between them (Figure 3). Thus, the total performance index (PITotal) presented tight positive correlations with the absorbed (ABS/CS; r = 0.62; p < 0.05) and dissipated fluxes (DI/CS; r = 0.90; p < 0.05), but the relationship was markedly negative with Fv/Fm (r = −0.85; p < 0.05). This latter parameter of PSII efficiency was more related to fluxes of primary photochemistry (TR/CS and ET/CS; Figure 3). Furthermore, the direction of fundamental, globally described, leaf trait– trait relationships (those between LMA, N and P; Wright et al., 2004) were maintained within cytotypes. However, we found differences when comparing the coordination of functional traits among ploidies considering the complete network. The 2× cytotype had the highest integration of LES traits, which decreased with polyploidization, being 6× and 12× cytotypes FIGURE 1 Phylogenetic principal component analysis (pPCA) obtained for the sampled populations of Dianthus broteri, using soilrelated variables regarding nutrient composition and physicochemical properties. Only the two axes with greater inertia explained are represented. Panel (A) shows the distribution of the selected variables loading on the axes. Panel (B) presents the edaphic niches of the four cytotypes (ellipses encompass occurrence points) in the two main axes. Panels (C) and (D) are a breakdown of (B) representing each axis separately (soil- PC1 and PC2, respectively). These show the edaphic range of the ploidies as violin plots, where white circles correspond to the median, thick black bars to the interquartile range and thin bars to the 95% confidence interval. Different letters indicate cytotypes that are significantly different from each other (oneway ANOVA: p < 0.05, Tukey's HSD test: α = 0.05). Letters with asterisks designate marginally significant differences (Tukey's HSD test: 0.05 < α < 0.10) | 611 Journal of Ecology LÓPEZJURADO Et AL. those with the lowest INTc and percentages of maximum integration (%INTmax; Table 1). Additionally, the Mantel tests revealed a greatly dissimilar structure of LES traits (lowest Pearson's r) between diploids and the other cytotypes (4×, 6× and 12×; Table S2). In contrast, these three cytotypes showed similar LES traits covariation matrices, although the Mantel tests were not significant (p > 0.05; Table S2). In both soil and LES pPCAs, we only found significant phylogenetic signal (K ≥ 1, p < 0.05) for the most explanatory axes (PC1s). The variation in these axes better fitted the OU4 model (i.e. the one considering a scenario with four optima, corresponding to the four cytotypes), showing AIC weights above 75 (virtually 100 for LESPC1; Table 2). These results supported the conservatism of the edaphic niche and functional traits at population level (i.e. within ploidies), but the divergence among D. broteri cytotypes. Nevertheless, both sets of characteristics were not correlated according to the results of the procrustean tests. The two pairs of pPCA scores did not show a pattern of multivariate concerted evolution, based on their low and nonsignificant correlation (r = 0.36; p > 0.05). In contrast to the abovementioned divergence in soil nutrient contents among D. broteri cytotypes, their leaves showed a clearly distinct pattern, marked by significantly higher contents of essential macronutrients in the 12× cytotype (C, N, P and K; p < 0.05; Table S3) and a greater variation in the 2× cytotype (as indicated its overall higher standard errors; Table S3). Thus, nutrient acquisition and accumulation rates and efficiencies differed among ploidies, especially affecting 6× and 12×. With regard to the variables accounting for the difference in macronutrient concentrations between leaf FIGURE 2 Phylogenetic principal component analysis (pPCA) obtained for the sampled populations of Dianthus broteri, using leaf economics spectrum (LES) and other functional variables. Only the two axes with greater inertia explained are represented. Panel (A) shows the distribution of the selected variables loading on the axes and panel (B) presents the functional space of the four cytotypes (ellipses encompass occurrence points) in the two main axes. Panels (C) and (D) are a breakdown of (B) representing each axis separately (LESPC1 and PC11, respectively). These show the functional range of the ploidies as violin plots, where white circles correspond to the median, thick black bars to the interquartile range and thin bars to the 95% confidence interval. Different letters indicate cytotypes that are significantly different from each other (oneway ANOVA: p < 0.05, Tukey's HSD test: α = 0.05) 612 | Journal of Ecology LÓPEZJURADO Et AL. and soil, dodecaploids presented again the overall most elevated values (Figure 4A– C). Thus, the 12× cytotype had significantly greater differences in N than the rest of ploidies (ANOVA: F3,8 = 6.66; p < 0.05; Figure 4A). Dodecaploids also exhibited marginally higher differences in K (F3,8 = 2.23; p = 0.08; Figure 4B) and P (F3,8 = 2.39; p = 0.07; Figure 4C) compared to 2× and 6× cytotypes, respectively. We also detected large differences in several micronutrients such as Na, which was significantly higher in the 6× cytotype than in the remaining ones (F3,8 = 5.17; p < 0.05; Table S3). 4 | DISCUSSION 4.1 | Soil niche divergence and expansion associated with polyploidy Most of the variation in the edaphic niche of D. broteri would be explained by processes associated with the shortterm effects of polyploidization as suggested the global structures found in the two main soilpPCA axes. The immediate phenotypic novelties derived from genome duplication, such as those found in D. broteri synthetic FIGURE 3 Correlation network of functional variables regarding nutrient, water and light use in the sampled populations of Dianthus broteri. Tightly related variables appear closer and lines connecting them are wider and have darker colours (red for positive relationships and blue for negative ones). Only correlations with Pearson's r > 0.3 are drawn TABLE 1 Summary of the integration analysis in the leaf economics spectrum traits for each Dianthus broteri cytotype Ploidy NINTc Lower CI Upper CI %INTmax 2×64 3.04 ± 1.00 2.29 3.33 25.33 4×66 2.95 ± 0.90 2.49 3.10 24.58 6×64 2.39 ± 0.93 1.83 2.52 19.92 12×81 2.70 ± 0.88 2.41 2.76 22.50 Note: Values represent mean ± SD. Abbreviations: %INTmax, percentage of the total maximum integration; CI, confidence interval; INTc, corrected phenotypic integration index; N, total number of sampled individuals. TABLE 2 Performance of models for the evolution of leaf economic traits and soil niche properties in Dianthus broteri. Only pPCA axes with significant phylogenetic signal are represented. Models are ranked according to their Akaike information criterion (AIC) weight LESPC1 Soil- PC1 Model AIC Weight Model AIC Weight OU4 39.2 99.1 OU4 48.3 77.9 BM 49.3 0.6 BM 51.8 13.8 OU1 51.3 0.2 OU1 53.8 5.1 NP 55.1 0.0 NP 54.7 3.2 FIGURE 4 Violin plots representing differences in the main macronutrients between leaves from the four Dianthus broteri cytotypes and soils from their respective niches. Panel (A) corresponds to nitrogen, panel (B) to potassium and panel (C) to phosphorus. For each element, higher values represent more efficient acquisition rates from soil to leaves. In these plots, white circles depict the median, thick black bars correspond to the interquartile range and thin bars show the 95% confidence interval. Different letters indicate cytotypes that are significantly different from each other (oneway ANOVA: p < 0.05, Tukey's HSD test: α = 0.05). Letters with asterisks designate marginally significant differences (Tukey's HSD test: 0.05 < α < 0.10) | 613 Journal of Ecology LÓPEZJURADO Et AL. neotetraploids for several functional traits similar to the ones explored in this study (gasexchange and photochemical variables; Domínguez- Delgado et al., 2021), would be the main drivers of the divergence in soil niche. However, an early postpolyploidization character displacement caused by competition for niche space (Beans, 2014) or niche filling cannot be ruled out. The highest cytotypes (6× and 12×) would have been displaced to stressful environmental edges with much poorer soils in terms of macronutrient content and potential availability than the lower cytotypes (López- Jurado et al., 2019). Furthermore, the divergence pattern was accompanied by an expansion of the edaphic niche in polyploids, which could be linked with amplified ecological tolerances (Molina- Henao & Hopkins, 2019; te Beest et al., 2012) and/or enhanced competitive abilities (Rey et al., 2017). 4.2 | Divergent positioning along the LES associated with polyploidy Although soil variables are supposed to be the strongest environmental predictors of leaf traits (Maire et al., 2015), the results of the procrustean tests linking edaphic niche characteristics and LES variables suggested that they are not coevolving. Nevertheless, functional traits related to nutrient (N, P and K) and water economy (δ13C) explained by LESPC1 also showed an early divergence (global structure) associated with polyploidization. Leal- Bertioli et al. (2017) and Wei et al. (2020) also demonstrated large effects of polyploidy per se on leaf morphology and anatomy by finding analogous results in wild and induced polyploids. Moreover, Münzbergová and Haisel (2019) found a similar pattern for the photosynthetic pigment contents in autotetraploids. In agreement with soil- PC1, each cytotype reached a different LESPC1 optimum given the highest performance of OU4 models. Thus, divergent selection was the major force driving the evolution of D. broteri polyploids in functional traits and soil niche characteristics, which agreed with a scenario of intercytotypic competitive interactions (Laport et al., 2013; Rey et al., 2017). In contrast to its most restricted edaphic niche, the 2× cytotype had the overall greatest variation in leaf nutrient concentrations. This fact could be explained by larger pools of genetic variation in this cytotype than the other ploidies (Balao et al., 2010). Notwithstanding the higher genetic variation of the diploids, they showed reduced leaf nutrient contents, which were similar to those of 4× and 6× cytotypes. However, we would have expected the leaf nutrients to be increased in the polyploids based on the greater need for these in taxa with larger genomes (Leitch & Leitch, 2008) and the complex polygenic uptake mechanisms for essential macronutrients (N, P and K; Lea & Azevedo, 2006; Vance et al., 2003; Véry & Sentenac, 2003), which would be facilitated by postpolyploidization processes (Cheng et al., 2018). This finding could be explained by the phenotypic diploidization phenomenon, an evolutionary tendency of returning to diploidlike phenotype which has been recently demonstrated for photosynthetic characteristics of established D. broteri 4× individuals (Domínguez- Delgado et al., 2021). Oppositely, the highestorder cytotype (12×) showed, as expected, an elevated leaf accumulation of K, P and N. On the basis of the low- P and N soils where the 12× ploidy is distributed, their elevated content in leaves would be useful for dealing with uncertainties in availability (Clarkson et al., 2005; Ostertag, 2010; Ruiz- Navarro et al., 2016). Moreover, dodecaploids would benefit from a higher K content because this element has a key role in numerous stress response mechanisms (e.g. drought, cold or salinity; Wang et al., 2013). In the D. broteri complex, LESPC1 reflected the acquisitiveconservative continuum (fast– slow, sensu Reich, 2014), which was in accordance with wateruse strategies (through δ13C) as Prieto et al. (2018) described in other Mediterranean species. Meanwhile, the highest ploidy (12×) occupied an extreme of the continuum, the second highest ploidy (6×) was located at the opposite extreme, and 2× and 4× filled more intermediate spaces. The 6× cytotype showed a nutrient and waterconservative or slow strategy (i.e. high LMA, low leaf nutrient accumulation plus high δ13C), whereas the 12× cytotype presented a nutrient and wateracquisitive or rapid strategy (i.e. low LMA, high leaf nutrient accumulation plus low δ13C). The 12× cytotype showed the lowest intrinsic wateruse efficiency (iWUE; estimated from leaf δ13C), which would explain its elevated leaf nutrient concentrations since the low iWUE is offset with high nutrient uptake efficiency through increased stomatal conductance (gs; Cramer et al., 2009). Nevertheless, the other highly droughtexposed cytotype (6×) would prioritize water conservation (causing high iWUE) through the alleged reduction of gs. The functioning at lower internal CO2 concentrations (Ci) would help to reduce water loss (Goud et al., 2019). Additionally, the higher LMA in hexaploids would allow a better adjustment of leaf water potential with small changes in water content (optimizing water uptake from dry soils; Niinemets, 2001) and, at the same time, the alleviation of extreme heat stress (Leigh et al., 2012) and the adaptation to high light intensities (Poorter et al., 2009). Congruently, high- LMA species have been extensively associated with semiarid Mediterranean environments as the one inhabited by the 6× cytotype (e.g. Costa- Saura et al., 2016; de la Riva et al., 2016). However, increased LMA entails generally low photosynthetic rates (AN), since it would limit CO2 diffusion through the mesophyll (Tomás et al., 2013), as well as N investment in photosynthetic proteins because of the higher cell wall mass fraction (Onoda et al., 2017). Interestingly, the rapid or acquisitive strategy characterizing the 12× cytotype cannot be maintained indefinitely in a perennial plant (Onoda et al., 2017). Therefore, this functional strategy should be promoted in periods of opportunity, limiting the growing season to when water and nutrients become sufficiently available (Querejeta et al., 2018). These adaptations would enable D. broteri dodecaploids to cope with the remarkably seasonal Mediterranean climate under which they live, with highly variable precipitation events (Deitch et al., 2017) and long dry spells (Sánchez et al., 2011). Supporting this hypothesis, we have found stomatal closure and increased iWUE during summer in the 12× cytotype (López- Jurado et al., 2019).