Recent and ancient evolutionary events shaped plant elemental composition of edaphic endemics: a phylogeny-wide analysis of Iberian gypsum plants.

The analysis of plant elemental composition and the underlying factors affecting its variation are a current hot topic in ecology. Ecological adaptation to atypical soils may shift plant elemental composition. However, no previous studies have evaluated its relevance against other factors such as phylogeny, climate or individual soil conditions. We evaluated the effect of the phylogeny, environment (climate, soil), and affinity to gypsum soils on the elemental composition of 83 taxa typical of Iberian gypsum ecosystems. We used a new statistical procedure (multiple phylogenetic variance decomposition, MPVD) to decompose total explained variance by different factors across all nodes in the phylogenetic tree of target species (covering 120 million years of Angiosperm evolution). Our results highlight the relevance of phylogeny on the elemental composition of plants both at early (with the development of key preadaptive traits) and recent divergence times (diversification of the Iberian gypsum flora concurrent with Iberian gypsum deposit accumulation). Despite the predominant phylogenetic effect, plant adaptation to gypsum soils had a strong impact on the elemental composition of plants, particularly on sulphur concentrations, while climate and soil effects were smaller. Accordingly, we detected a convergent evolution of gypsum specialists from different lineages on increased sulphur and magnesium foliar concentrations.

Introduction

Plant life, in all its diversity, is the result of the combination of up to 30 chemical elements (Ågren, 2008). The ultimate way in which such elements combine and give rise to the elemental composition of plants depends on processes of nutrient uptake, use, storage and translocation (Baxter, 2009). Although some of these processes are conserved across the phylogeny, many are species specific, leading to an ‘elemental fingerprint’ for each taxon (i.e. the ‘ionome’, sensu Lahner et al., 2003 or ‘elementome’ sensu Li et al., 2008). The elemental composition of plant taxa is recently gaining recognition as a fundamental concept in ecology, as it encompasses each taxon nutritional requirements, which are ultimately the result of its function and life strategy (Peñuelas et al., 2019). The concentrations of the different elements can be considered as axes of variation within a multidimensional space, encompassing functional information and shaping the ‘stochiometric niche’ (sensu González et al., 2017), or ‘biogeochemical niche’ (sensu Peñuelas et al., 2019). Identifying the different factors that determine species elemental composition is, therefore, a current key goal in ecology (Jeyasingh et al., 2017; Peñuelas et al., 2019).

There is ample evidence that evolutionary history is a strong determinant of plant elemental composition (Thompson et al., 1997; Broadley et al., 2004; Watanabe et al., 2007; Neugebauer et al., 2018). The phylogenetic signal seems to be generally stronger in macroelements than microelements, which are mostly determined by environmental factors (Zhao et al., 2016; de la Riva et al., 2018). Previous studies have evaluated phylogenetic effects on plant elemental composition by partitioning elemental variance across different taxonomic levels (i.e. order, family, genus, species; e.g. Broadley et al., 2003; Watanabe et al., 2007). For instance, Hao et al. (2015) detected stronger effects at the subfamily than at the genus level. While indicative of phylogenetic relatedness, taxonomic levels do not correspond to phylogenetic distances or divergence times (Magallón & Castillo, 2009; Massoni et al., 2015). Consequently, new analytical methods that partition elemental variability across the divergence time of species are required to have an integrated overview of phylogenetic effects on plant elemental composition.

The environmental conditions in which individual plants thrive have also an impact on their elemental composition (Han et al., 2011; Zhang et al., 2012; Sardans et al., 2016). For example, both mean annual precipitation (MAP) and mean annual temperature (MAT) correlate with plant elemental concentrations (Zhang et al., 2012; Sardans et al., 2016). Soil is considered one of the main drivers of plant elemental composition (Marschner, 2012). Contrastingly, most previous studies have found only a limited effect of soil properties on the elemental composition of plants growing in the wild (Thompson et al., 1997; Zhao et al., 2016). The resolution level of soil data included in these analyses was likely too coarse to account for meaningful plant–soil interactions (Thompson et al., 1997). Although sampling of soil adjacent to plant individuals included in elemental analyses may improve the evaluation of plant–soil interaction effects on plant elemental composition (Stein et al., 2017), studies incorporating this approach are still scarce (but please refer to Salmerón‐Sánchez et al., 2014; Stein et al., 2017; Pillon et al., 2019).

Species adaptation to certain environmental conditions can also strongly alter the elemental composition of plants (Huang & Salt, 2016). This is particularly true for plants growing on atypical substrates, which have to cope with soil nutrient imbalances (van der Ent et al., 2018; Matinzadeh et al., 2019; Merlo et al., 2019). Some species tolerate atypical soils by excluding phytotoxic elements or nutrients found in excess to keep elemental homeostasis (van der Ent et al., 2018; Matinzadeh et al., 2019; Merlo et al., 2019). However, specialisation to atypical substrates often involves shifts in plant elemental composition (Verboom et al., 2017), frequently related to extreme accumulation of excess elements (Pillon et al., 2010; van der Ent et al., 2018; Merlo et al., 2019). For example, species adapted to serpentine soils are commonly metal hyperaccumulators, reaching markedly high concentrations of Ni, Zn, Cd, Co, Mn, Al or Pb depending on soil pH (Stein et al., 2017; van der Ent et al., 2018). Halophytes accumulate several orders of magnitude higher Na concentrations than co‐occurring species (Matinzadeh et al., 2019). Similarly, some species exclusive to calcareous soils show markedly higher concentrations of Ca and Mg compared with neighbouring plants (Hao et al., 2015). The flora of atypical soils is, therefore, a perfect system to evaluate the relevance of ecological adaptations on the elemental composition of species. However, to our knowledge, no previous attempts have been made to estimate the relative contribution of such adaptations on the elemental composition of plants. Ascertaining if they are the consequence of certain pre‐adaptations of soil specialists (implying a strong phylogenetic signal), or simply an expression of micro‐evolutionary processes with low phylogenetic signal, remains a critical issue to unveil how soil specialisation evolved.

Extending over 100 million hectares worldwide, gypsum soils are amongst the most widespread atypical substrates of the world (Eswaran & Gong, 1991). They occur on all continents in areas where arid and semiarid conditions prevent gypsum from being leached (Eswaran & Gong, 1991). The high concentration of gypsum in soil leads to special physical and chemical conditions that pose serious restrictions to plant life, limiting the development of agriculture and conditioning the livelihood of millions of people (Verheye & Boyadgiev, 1997; Palacio & Escudero, 2014; Escudero et al., 2015). The moderate solubility of gypsum (c. 2.4 g l−1) results in abnormally high Ca and sulphate concentrations, toxic for some plants (Ernst, 1998). Such high Ca and sulphate concentrations decrease nutrient availability in the soil due to the saturation of the soil solution with Ca2+ and sulphate ions (FAO, 1990), which leads to overall low nutrient retention ability (Casby‐Horton et al., 2015). High soil Ca levels enhance sulphate uptake by plants, while a low N supply may impair the N : S balance for protein synthesis, further exacerbating excess sulphate accumulation in plants (Rennenberg, 1984). Similar to other atypical soils, the extremely restrictive conditions of gypsum soils contrast with the highly diversified flora they sustain, rich in endemic and specialised species (Mota et al., 2011; Musarella et al., 2018; Ochoterena et al., 2020), which is a conservation priority of international concern (European Community, 1992).

Depending on their affinity to gypsum, plants growing on gypsum soils are classified as gypsophiles (plants that grow only on gypsum soils) or gypsovags (plants that grow both on and off gypsum) (Meyer, 1986). Gypsophiles can further be segregated into: (1) wide gypsophiles, species occurring on most gypsum outcrops within a given region, which are considered specialised to gypsum and putatively belong to old gypsophilic lineages; and (2) narrow gypsophiles, locally distributed species that mostly belong to young gypsum lineages (Palacio et al., 2007; Muller et al., 2017). Previous studies on gypsum species from Spain, the Chihuahuan desert and Turkey indicate that widely distributed gypsophiles tend to show higher Ca, S and Mg foliar concentrations (elements found in excess in gypsum soils) than their neighbour gypsovags and narrow gysophiles (Duvigneaud & Denayer‐De Smet, 1968; Palacio et al., 2007; Bolukbasi et al., 2016; Muller et al., 2017). For Chihuahuan desert plants, differences remained when phylogenetic effects were removed from the analyses (Muller et al., 2017). This ability could be related to plant ecological adaptation to gypsum soils (Palacio et al., 2014). However, studies specifically addressing the effect of evolutionary history on the elemental composition of gypsum plants are lacking.

The aims of this study were to:

Evaluate the relevance of phylogeny, affinity for gypsum soils and environmental factors (soil and climate) on the elemental composition of species from gypsum ecosystems in Iberia.

Assess the relationship between the ability to accumulate high foliar Ca, Mg and S concentrations and species specialisation to gypsum.

The following hypotheses were tested:

Phylogeny will be the most important factor explaining the variability on the elemental composition of plants from gypsum habitats in Iberia, but environmental and ecological (i.e. gypsum affinity) factors will also play a relevant role. In particular, we expect affinity to gypsum soils to explain an important proportion of the variability in S, Ca and Mg concentrations.

Phylogenetic effects will vary across the divergence time of taxa and trends will be different among elements.

Widely distributed Iberian gypsophyle species will accumulate more Ca, Mg and S in leaves than closely related gypsovags and narrow gypsophiles, independent of their phylogenetic origin.

We tested these hypotheses on a broad dataset including elemental concentrations of several key elements in 83 taxa. Our approach combined multivariate (including 11 elements), univariate (on 15 elements) and phylogenetic statistical tools to evaluate the effects on the elemental composition of plants as a whole and by each individual element separately. To this end, we used a new statistical procedure that allows evaluation of phylogenetic effects and their significance across the divergence time of taxa (multiple phylogenetic variance decomposition, MPVD).

Materials and Methods

Study species and sites

We investigated 83 taxa representing 18 families and 10 orders (Table 1). The affinity for gypsum was assigned based on plant distribution on gypsum soils in two ways: (1) a categorical approach including 15 widely and 27 narrowly distributed gypsophiles and 41 gypsovags (Meyer, 1986; Palacio et al., 2007), and (2) a continuous gypsophily index proposed by Mota et al. (2011) (Supporting Information Methods S1). The selected taxa were good representatives of the gypsum flora of Iberia, comprising 56% of its gypsophile taxa (Mota et al., 2011) plus highly dominant gypsovag species. Whenever possible, the selection included congeneric or confamiliar representatives of the three groups of gypsum affinity considered.

Table 1

Details of taxa and sampling locations included in the study.

Order, family and taxon name	Gypsum affinity 1		Life cycle 2	Locality 3	Coordinates	n 4
	Class	GI
Apiales
Apiaceae
*Ferula communis L.	GV	2	P	Alfajarín (Z)	41°36′34.2″N, 00°41′23.1″W	5
*Ferula loscosii (Lange) Willk.	W	3.75	P	Candasnos (Hu)	41°25′07.3″N, 00°06′49.3″E	2
Asterales
Asteraceae
*Achillea santolinoides Lag.	N	3.5	P	Villena (A)	38°36′40.45″N, 00°53′45.86″W	5
*Artemisia herba‐alba Asso	GV	2	P	Villamayor (Z)	41°41′30.1″N, 00°44′48.8″W	5
*Centaurea hyssopifolia Vahl	N	4.63	P	Valdemoro (M)	40°11′19.4″N, 03°36′05.5″W	5
*Launaea fragilis (Asso) Pau	GV	2	P	Villamayor (Z)	41°42′02.5″N, 00°44′47.7″W	5
*Launaea pumila (Cav.) Kuntze	W	3.22	P	Villamayor (Z)	41°42′34.1″N, 00°44′03.3″W	5
*Santolina viscosa Lag.	N	3.92	P	Venta de los Yesos (Al)	37°04′57.98″N, 02°17′23.70″W	5
Campanulaceae
*Campanula fastigiata Dufour ex Schult	W	4.58	A	Villafranca de Ebro (Z)	41°35′34.2″N, 00°31′32.4″W	5
Brassicales
Brassicaceae
*Boleum asperum (Pers.) Desv.	GV	3.03	P	Monegrillo (Z)	41°37′37.9″N, 00°27′3.30″W	4
*Brassica repanda subsp. gypsicola Gómez‐Campo	N	3.45	P	Yebra (Gu)	40°20′51.7″N, 02°56′23.4″W	5
*Erucastrum nasturtiifolium (Poir.) O.E.Schulz	GV	2	P	Ibars de Noguera (L)	41°50′31.3″N, 00°35′19.0″E	5
*Lepidium cardamines L.	N	3.44	A	Orusco de Tajuña (M)	40°16′06.3″N, 30°08′09.4″W	5
*Lepidium subulatum L.	W	4.91	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.30″W	5
				Venta de Yesos (Al)	37°04′57.98″N, 02°17′23.70″W	4
*Matthiola fruticulosa (L.) Maire	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Vella pseudocytisus subsp. paui Gómez‐Campo	GV	3.47	P	Villastar (Te)	40°14′35.0″N, 01°08′31.0″W	8
*Vella pseudocytisus L. subsp. pseudocytisus	N	3.23	P	Aranjuez (M)	40°02′03.4″N, 03°33′07.7″W	5
Resedaceae
*Reseda barrelieri Bertol. ex Müll.Arg.	GV	2	P	Petrel (A)	38°26′44.9″N, 00°46′50.4″W	5
*Reseda lutea L.	GV	2	P	Villamayor (Z)	41°42′00.8″N, 00°43′21.2″W	19
*Reseda phyteuma L.	GV	2	P	Carchelejo (J)	37°38′24.2″N, 03°6′54.02″W	5
*Reseda stricta Pers.	W	3.97	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Reseda suffruticosa Loefl.	N	4.05	P	Morata de Tajuña (M)	40°12′24.5″N, 03°24′54.0″W	5
*Reseda undata L.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
Caryophyllales
Amaranthaceae
*Bassia scoparia (L.) A.J.Scott	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Salsola vermiculata L.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
Caryophyllaceae
*Arenaria aggregata subsp. cavanillesiana (Font Quer & Rivas Goday) Greuter & Burdet ex Lop	N	3.33	P	San Pedro Palmiches (Cu)	40°25′54.1″N, 02°23′50.7″W	5
*Gypsophila struthium L. subsp. hispanica (Willk.) G. López	W	4.69	P	Villamayor (Z)	41°42′34.1″N, 00°44′03.3″W	10
				Almunia de San Juan (Hu)	41°56′16.7″N, 0°15′31.7″E	5
*Gypsophila tomentosa L.	W	3.34	P	Rueda de Jalón (Z)	41°38′17.9″N, 01°01′10.5″W	5
*Gypsophila bermejoi G.López	N	3.77	P	Yebra (Gu)	40°20′51.7″N, 02°56′23.4″W	5
*Herniaria fruticosa L.	W	4.05	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	9
Frankeniaceae
*Frankenia pulverulenta L.	GV	2	A	Monegrillo (Z)	41°36′54.8″N, 00°29′25.2″W	5
*Frankenia thymifolia Desf.	W	3.28	P	Villafranca de Ebro (Z)	41°35′44.2″N, 00°31′25.4″W	5
Ericales
Primulaceae
*Coris hispanica Lange.	N	4.39	P	Sorbas (Al)	37°4′26.53″N, 20°5′32.61″W	5
*Coris monspeliensis L.	GV	2	P	Alfajarín (Z)	41°36′34.2″N, 00°41′23.1″W	5
Fabales
Fabaceae
*Astragalus guttatus Banks & Sol.	N	3	P	Quesada (J)	37°37′9.10″N, 03°8′6.98″W	1
*Astragalus oxyglottis M.Bieb.	N	3.64	P	Alicún de Ortega (J)	37°36′24.43″N, 03°6′51.73″W	1
*Hedysarum boveanum Basiner subsp. europaeum Guitt. & Kerguelen	GV	2	P	Villamayor (Z)	41°42′0.8″N, 0°43′21.2″W	5
*Hedysarum boveanum Basiner subsp. palentinum Valdés	N	4.14	P	Tudela de Duero (Va)	41°34′12.4″N, 04°33′43.0″W	4
*Ononis fruticosa L.	GV	2	P	Paracuellos de Jiloca (Z)	41°19′24.3″N, 01°36′11.7″W	5
*Ononis rotundifolia L.	GV	2	P	Plan (HU)	42°35′0.6″N, 00°20′0.2″E	5
*Ononis tridentata L. subsp. tridentata Devesa & G. López	W	4.43	P	Villamayor (Z)	41°42′34.1″N, 00°44′03.3″W	10
*Ononis tridentata L. subsp. crassifolia (Dufour ex Boiss.) Nyman	N	4.53	P	Escúzar (Gr)	37°3′33.5″N, 03°44′53.9″W	5
*Ononis tridentata subvar. edentula (Webb ex Willk.) O. Bolós & Vigo 5	W	4.16	P	Orcheta (A)	38°34′43.9″N, 00°14′59.3″W	5
Lamiales
Lamiaceae
*Rosmarinus officinalis L.	GV	2	P	Peñaflor (Z)	41°45′52.8″N, 00°46′26.1″W	5
*Salvia officinalis subsp. lavandulifolia (Vahl) Gams	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Sideritis ilicifolia Willd.	GV	2	P	Almunia de San Juan (Hu)	41°56′23.2″N, 00°15′51.1″E	5
*Teucrium balthazaris Sennen	N	4.39	P	Limaria (Al)	37°24′19.44″N, 02°4′5.85″W	5
*Teucrium capitatum L.	GV	2	P	Azanuy (Hu)	41°59′25.5″N, 00°17′22.1″E	5
*Teucrium lepicephalum Pau	N	4.5	P	Orcheta (A)	38°34′43.9″N, 00°14′59.3″W	5
*Teucrium libanitis Schreb.	N	4.88	P	Elda (A)	38°28′24.2″N, 00°50′47.6″W	5
*Teucrium pumilum Loefl. ex L.	N	4.45	P	Yebra (Gu)	40°20′51.7″N, 02°56′23.4″W	5
*Teucrium turredanum Losa & Rivas Goday	N	4.67	P	Sorbas (Al)	37°4′26.53″N, 2°5′32.61″W	5
*Thymus lacaitae Pau	N	4.54	P	Valdemoro (M)	40°11′19.4″N, 03°36′05.5″W	5
*Thymus vulgaris L.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Thymus zygis L.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
Plantaginaceae
*Chaenorhinum rupestre (Guss.) Speta	W	4.79	A	Villafranca de Ebro (Z)	41°35′44.2″N, 00°31′25.4″W	7
*Chaenorhinum grandiflorum (Coss.) Willk. subsp. carthaginense (Pau) Benedí	GV	3.11	A	Almagro (Ab)	37°21′13.83″N, 01°54′53.32″W	2
				La Malahá (Gr)	37°2′38.26″N, 3°49′38.15″W	2
*Chaenorhinum grandiflorum (Coss.) Willk. subsp. grandiflorum	N	4.23	A	Venta de los Yesos (Al)	37°4′57.98″N, 02°17′23.7″W	1
*Chaenorhinum reyesii (C. Vicioso & Pau) Benedí	W	4.83	A	Villafranca de Ebro (Z)	41°35′44.2″N, 00°31′25.4″W	10
*Chaenorhinum rubrifolium (Robert & Castagne ex DC.) Fourr.	GV	2	A	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	1
*Plantago albicans L.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
Malpighiales
Euphorbiaceae
*Euphorbia minuta Loscos & Pardo subsp. moleroi P. Monts.	N	3	P	Almunia de San Juan (Hu)	41°56′23.2″N, 00°15′51.1″E	5
Linaceae
*Linum castroviejoi Mart. Labarga, Pedrol & Muñoz Garm.	N		P	Ibars de Noguera (L)	41°50′31.3″N, 00°35′19.0″E	5
*Linum strictum L.	GV	2	A	Almagro (Ab)	37°21′15.3″N, 01°52′53.3″W	1
				Monegrillo (Z)	41°36′54.8″N, 00°29′25.2″W	5
*Linum suffruticosum L.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
Malvales
Cistaceae
*Fumana ericifolia Wallr.	GV	2	P	Peñaflor (Z)	41°46′6.6″N, 00°45′17.5″W	5
*Fumana thymifolia (L.) Spach	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Helianthemum alypoides Losa & Rivas Goday	N	4.85	P	Sorbas (Al)	37°4′26.53″N, 02°5′32.61″W	5
*Helianthemum apenninum subsp. stoechadifolium (Brot.) Samp.	GV	2	P	Baltanás (Pa)	41°55′37.4″N, 04°16′24.7″W	3
*Helianthemum marifolium subsp. conquense Borja & Rivas Goday ex G.López	N	4.43	P	Yebra (Gu)	40°20′51.7″N, 02°56′23.4″W	5
*Helianthemum marifolium (L.) Mill. subsp. marifolium	GV	2	P	Barbastro (Hu)	41°59′29.3″N, 00°04′49.2″E	5
*Helianthemum salicifolium (L.) Mill.	GV	2	A	Monegrillo (Z)	41°36′54.8″N, 0°29′25.2″E	5
*Helianthemum squamatum (L.) Dum.Cours.	W	4.87	P	Venta de los Yesos (Al)	37°4′57.98″N, 02°17′23.70″W	5
				Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	10
*Helianthemum syriacum (Jacq.) Dum.Cours.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	15
Thymelaceae
*Thymelaea tinctoria (Pourr.) Endl.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	10
Poales
Cyperaceae
*Schoenus nigricans L.	GV	2	P	Azanuy (Hu)	41°59′25.5″N, 00°17′22.1″E	3
Poaceae
*Agropyron cristatum (L.) Gaertn. subsp. pectinatum (M. Bieb.) Tzvelev.	W	3.66	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Brachypodium retusum (Pers.) P.Beauv.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Koeleria vallesiana (Honck.) Bertol. ex Schult.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Koeleria vallesiana subsp. castellana (Boiss. & Reut.) Domin	N	3.38	P	Yebra (Gu)	40°20′51.7″N, 02°56′23.4″W	5
*Lygeum spartum Loefl. ex L.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Stipa barbata Desf.	GV	2	P	Villamayor (Z)	41°42′34.1″N, 00°44′3.3″W	5
*Vulpia gypsophila (Hack.) Nyman	N	3.92	P	Valdemoro (M)	40°11′19.4″N, 03°36′05.5″W	1
				Orusco de Tajuña (M)	40°16′06.3″N, 03°08′09.4″W	1

Plant nomenclature followed the Taxonomic Name Resolution Service (TNRS) (Internet). iPlant Collaborative. v.4.0 (accessed 8 May 2020). Available from: http://tnrs.iplantcollaborative.org or Flora Ibérica (Castroviejo, 1986–2012) for local taxa not included in TNRS, unless otherwise indicated.

Gypsum affinity either as a class (GV, gypsovag; N, narrow gypsophile; W, wide gypsophile) or as the Gypsum Index proposed by Mota et al. (2011).

P, perennial; A, annual.

Geographic origin: A, Alicante; Ab, Albacete; Al, Almería; Cu, Cuenca; Hu, Huesca; Gu, Guadalajara; Gr, Granada; J, Jaén; L, Lérida; M, Madrid; Pa, Palencia; Te, Teruel; Va, Valladolid; Z, Zaragoza.

Number of individuals sampled for leaf elemental composition analyses.

Sensu Bolòs & Vigo (1984).

Taxa were collected from at least one population growing on gypsum. Gypsum soils were considered as those having more than 10% gypsum, which is higher than the threshold of 5% considered for gypsum soils (FAO, 1990). In fact, 93.5% of the soils included in the study contained over 20% gypsum. Sampling sites were spread across the Iberian Peninsula, covering most of the gypsum outcrops of the region (Fig. 1). They were mainly composed of semiarid Mediterranean shrublands and dominated by subshrubs, many of them restricted to gypsum soils (please refer to Rivas Goday et al., 1957 for a description of the Iberian gypsum vegetation). In general, collection sites had semiarid and dry Mediterranean climates, with intense summer droughts and cold winters (Fig. 1; Methods S1).

Fig. 1

Map showing the distribution of sampling sites (white dots) across the gypsum outcrops of the Iberian Peninsula (hatched shade), plus climatic diagrams of representative locations of the four main regions studied showing average monthly temperature (black lines), average monthly maximum and minimum temperatures (T) (top and low dotted lines, respectively) and average monthly precipitation (P) (grey bars) over 30 yr of records (1982–2012). Gypsum outcrops were drawn from information in Escavy et al. (2012). Meteorological data were obtained from climate‐data.org (https://es.climate‐data.org/) (accessed 6 June 2022). [Colour figure can be viewed at wileyonlinelibrary.com]

Plant and soil sampling

For elemental analyses, 2 g leaf samples were collected from five individuals per taxon (n = 5 replicates in most cases; please refer to Table 1; Methods S1). Leaf elemental composition is broadly used as a good representative of the elemental composition of plants (Watanabe et al., 2007). When plants were small or threatened, composite samples were collected by pooling leaves from different individuals. In some cases, the number of replicates had to be less due to plant scarcity, such as in Astragalus oxyglottis or A. guttatus, in which only one composite sample was gathered (Table 1). We collected leaf samples when plant growth had stopped and leaf nutrient composition was more stable, targeting fully expanded mature leaves and avoiding senescing, young or damaged leaves (Palacio et al., 2007). Leaf samples were rinsed in tap water, dried at 50°C to a constant weight and subsequently milled in a ball mill (MM200; Restch GmbH, Haan, Germany) to a fine powder before elemental analyses.

Soil samples of the first 0–20 cm of the soil were collected at each sampling site after removing the biological soil crust and O horizons. To have a better estimate of leaf‐soil relationships, 269 soil samples were collected paired to leaf samples, by collecting soil close to individuals sampled for leaf material. These samples were included in the analyses individually. When paired soil samples were not available, one to five different soil samples were collected to represent each plant collection site. Unpaired soil samples were processed and analysed individually and the average values per site were included in statistical analyses. Soil samples were air dried for at least 2 months and subsequently sieved through a 2‐mm sieve before chemical analyses.

Chemical analyses

Plant samples were dissolved in HNO3 : H2O2 (8 : 2) using Microwave Acid Digestion (speedwave MWS‐3+; Berghof, Eningen, Germany). The filtered extract solution was used to measure 13 elements Al, Ca, Cu, K, Li, Mg, Mn, Mo, Na, P, S, Ti and Zn using inductively coupled plasma‐optical emission spectrometry (ICP‐OES; Varian ICP 720‐ES, Agilent Technologies Inc., Mulgrave, Vic., Australia) with a detection limit of 0.025 μg g−1. Total nitrogen and carbon (N and C) concentrations were determined using an elemental analyser (TruSpec CN; Leco, St. Joseph, MI, USA). As it was not possible to obtain N and C concentrations for seven species, we ended with two data sets for multivariate leaf composition analyses: CN− dataset (83 spp., nine elements, n = 458) and CN+ data set (76 spp., 11 elements and n = 417).

Gypsum content in soils was measured by gravimetry according to Artieda et al. (2006). Carbonate content was determined by Bernard calcimetry, soil texture was estimated using a particle laser analyser (Mastersizer 2000 Hydro G; Malvern Panalytical Ltd, Malvern, UK) and soil pH and conductivity were measured with a pH/conductivity meter (Orio Star A215; Thermo Scientific, Waltham, MA, USA) by diluting samples with distilled water to 1 : 2.5 (w/v) and 1 : 5 (w/v), respectively. A subsample of each soil sample was finely ground, dissolved in HCl : HNO3 (9 : 3) by microwave digestion and analysed for the same elements and plant samples as described above. All elemental analyses were performed by Estación Experimental del Zaidín (EEZ‐CSIC) Analytical Services, Spain.

Phylogenetic reconstruction

We built a custom phylogenetic tree for the study species using R packages v.phylomaker (Jin & Qian, 2019), phytools (Revell, 2012) and adephylo (Jombart & Dray, 2008). First, we detected which of our target species were already included in the Angiosperm dated phylogeny provided by v.phylomaker. For the target species not represented in the phylogeny, we selected, based on the available bibliography, a close relative in the v.phylomaker phylogeny for replacement (Methods S1; Table S1). We performed ancestral state reconstruction using the contMap function to generate a visual representation of the desired element composition across the phylogeny and tested its phylogenetic signal with Bloomberg's K and Pagel's λ using phylosyg function (both functions were from phytools, Revell, 2012).

Statistical analyses

All statistical analyses were run in R v.3.6.3 (R Core Team, 2020). A full description of the statistical methods is found in Methods S1. For multivariate analyses we used the CN+ dataset. Data were transformed to Centre Log‐ratio coordinates (Aitchison, 1982). Univariate approaches used raw data from 15 elements. Variance partitioning among: (1) the evolutionary levels of phylogeny (i.e. divergence time); (2) species affinity to gypsum; and (3) environmental factors (first axis of a principal components analysis (PCA) with soil features (Table S2), mean annual temperature (MAT) and precipitation (MAP)), were performed with PERMANOVA, when the whole elemental composition at the individual level was considered, and generalised linear models (GLMs) fitted to a gamma distribution, for concentrations of individual elements. We implemented a new procedure, called multiple phylogenetic variance decomposition (MPVD), to evaluate total variance explained by phylogeny and other potential factors along evolutionary time and test for significance (Methods S1). Mantel tests (Mantel, 1967) were performed to test for correlations among the phylogenetic and the elemental composition distance matrices.

To obtain a deeper insight into the effects of ecological and environmental factors on the elemental composition of plants, we used generalised linear mixed models (GLMM) fitted with Markov chain Monte Carlo techniques (MCMC) using the mcmcglmm package (Hadfield, 2010). This Bayesian approach allowed informing models with the phylogenetic relationships of study species, which were incorporated as the random term in models. Models included also species affinity to gypsum and environmental factors (first axis of a PCA with soil features (Table S2), MAT and MAP) as fixed factors, while elemental concentrations of individual elements were the response variables. We ran 500 000 MCMC iterations, with a burn‐in period of 1000 iterations. Model convergence was evaluated by running models repeated times (Methods S1). The following priors were used in accordance to Verdú et al. (2012) and de Villemereuil & Nakagawa (2014), both for the random effects (G) and residual variance (R): (V = 1, nu = 0.02). A redundancy analysis (RDA) was run to evaluate the relationships among elemental concentrations in taxa with different affinities for gypsum soils.

Results

The analysis of the coefficients of variation (CV) of different elements indicated that C, N and P were the least variable elements (with a minimum of 11% for C), whereas Na, Cu, Li, Ti, Al, S and Mg were the most variable ones (with a maximum of 394.2% observed in Na; Table 2). K, Ca, Mn, Zn and Mo showed intermediate CVs close to 60–70% (Table 2). Some elements became less variable when groups of plants with different affinities for gypsum were analysed separately. Such was the case of S, which was markedly less variable within wide gypsophiles than in the other groups (Table 2).

Table 2

Average concentrations (mg g−1, except for Cu, Li, Mn, Mo, Ti and Zn for which data are in μg g−1) and coefficients of variation (CV, %) of different elements measured in the leaves of plants with different affinity for gypsum substrates (i.e. gypsovags, narrow gypsophiles and wide gypsophiles).

Elements	Gypsovag	CV	Narrow	CV	Wide	CV	CV (all)
Al	0.31 ± 0.02	115.3	0.34 ± 0.03	94.2	0.35 ± 0.04	119.3	111.3
C	430.65 ± 3.05	10.4	433.59 ± 5.35	12.3	390.19 ± 4.52	11.7	11.8
Ca	24.71 ± 1.23	75.4	25.53 ± 1.59	66.9	37.85 ± 2.08	58.6	70.8
Cu	17.22 ± 1.68	145.9	12.55 ± 0.69	42.7	19.53 ± 2.45	120.2	133.5
K	13.72 ± 0.61	66.6	10.58 ± 0.47	47.5	14.39 ± 0.88	65	64.7
Li	10.43 ± 1.18	116.8	18.79 ± 3.77	118.8	10.84 ± 1.38	81.8	118.8
Mg	4.06 ± 0.23	85.2	3.88 ± 0.29	79.6	7.88 ± 0.69	93.1	99.6
Mn	47.01 ± 1.87	60.2	37.74 ± 2.34	66.7	48.06 ± 2.76	61.4	62.6
Mo	5.09 ± 0.71	78.9	6.62 ± 0.95	49.6	5.61 ± 1.63	91.8	73.6
N	24.45 ± 0.71	42.5	23.03 ± 1.06	45.9	23.01 ± 0.84	36.8	42.1
Na	1.46 ± 0.38	390.2	0.73 ± 0.18	257.6	0.69 ± 0.20	307.7	394.2
P	1.28 ± 0.04	51.1	1.09 ± 0.06	62.7	1.23 ± 0.05	42.4	52.1
S	8.27 ± 0.45	81.4	11.78 ± 1.23	112.1	23.39 ± 1.47	67.1	99.6
Ti	6.61 ± 0.39	69.7	12.83 ± 1.73	122.9	8.17 ± 1.18	106.5	117.5
Zn	33.20 ± 1.35	61.4	38.68 ± 2.17	60.5	36.78 ± 2.71	78.8	66.6

The CV of all elements in all species analysed together is also shown in the last column. Concentration data are means ± SE. n = 228, 116 and 114 for gypsovags, narrow and wide gypsophiles, respectively except for Cu (n = 223, 61, 92); Li (n = 107, 35, 41); Mo (n = 32, 12, 10); N and C (n = 215, 99, 102) and Ti (n = 141, 83 and 54). Data correspond to the same dataset used in generalised linear mixed models (GLMMs).

Phylogeny explained most of the variability in leaf chemical composition among species, but affinity to gypsum soils and environmental conditions had also a significant effect

PERMANOVA showed that the effects of phylogeny, ecological adaptation to gypsum (gypsovag/narrow gypsophile/wide gypsophile) and environmental conditions (soil, MAT and MAP) on the elemental composition of leaves of Iberian gypsum plants were highly significant (Table 3). However, while phylogeny (i.e. the taxonomic family) accounted for 36.4% of the total variance explained, the affinity for gypsum soils explained 2.1%, MAT explained only 0.9%, and MAP and soil features accounted for < 0.9% of the explained variance each (Table 3). Multilevel pairwise comparisons indicated that, while all groups of gypsum affinity differed in their elemental composition, differences were larger among wide gypsophiles and the rest (Table S3).

Table 3

Output of PERMANOVA indicating the minimum estimated explained variance (TVE) of the foliar elemental composition of gypsum Iberian plants for each independent variable.

Factor	df	SS	F	P‐value	TVE (%)
Phylogeny (family)	16	895.3	16.8	< 0.001	36.4
Gypsum affinity	2	52.3	7.9	< 0.001	2.1
MAT	1	22.6	6.8	< 0.001	0.9
MAP	1	20.0	6.0	< 0.001	0.8
Soil (PC1)	1	16.1	4.9	< 0.001	0.7

Type III sum of squares are used. The model explained 46.8% of the variance. Please refer to Materials and Methods for further details on the analyses. df, degrees of freedom; MAT, mean annual temperature; MAP, mean annual precipitation; SS, sum of squares.

The results of MPVD for the elemental composition of Iberian gypsum plants showed that the proportion explained by phylogeny tended to increase at relatively recent divergence times (Fig. 2; shown by the null model), reaching up to 75.7% at the species level. However, the percentage of the variance explained by phylogeny was significantly higher than that explained by the null model across most of the divergence time of study species (from 130 up to 0.5 million years, Myr). The high percentage of variance for all elements explained by phylogeny at shallow phylogenetic times (c. 10–3 million years ago, Ma) supports the importance of recent evolutionary events for the elemental composition of Iberian gypsum plants (Fig. 2). However, the accumulated phylogenetic signal remained highly explanatory at the genus level (R 2 = 0.58). Furthermore, both at the family (R 2 = 0.36) and order level (R 2 = 0.26), the accumulated effect of the phylogeny remained marked. It is noteworthy that the affinity for gypsum soils became more explanatory when the R 2 of the phylogenetic distance decreased, reaching a maximum R 2 of 0.06 when the largest phylogenetic distances (equivalent to no phylogenetic sorting) were included in the model (Fig. 2).

Fig. 2

Results of multiple phylogenetic variance decomposition (MPVD) showing, for each element, how the variation explained by affinity to gypsum soils and phylogeny changed depending on the phylogenetic grouping applied in generalised linear models (please refer to Supporting Information Methods S1). For the entire elemental composition (last plot), PERMANOVA was applied to extract the minimum variation explained by each factor. The R 2 of each explanatory variable was plotted according to the divergence time derived from each of the 74 nodes used to generate the phylogenetic grouping. A line linking serial points across divergence time was plotted for each explanatory variable. Average results for null models were represented (dashed blue line) plus their 5–95% confidence interval (grey shade). To help interpretation, usual taxonomic classification levels (gen, genus; fam, family; ord, order; sp, species) are plotted using the mean values of our dataset. The effect of other factors such as climate (mean annual temperature (MAT) and mean annual precipitation (MAP)) and soil was too low to be noticeable in figures and therefore is not shown. Myr, million years. [Colour figure can be viewed at wileyonlinelibrary.com]

The analysis of the variance partitioning for each element separately indicated that phylogeny was the most important explanatory factor, having a particularly large effect on the concentrations of Ca, Mg, N and Na (Fig. 2; Table 4). The effect of phylogeny was less for Al, and intermediate for the rest of elements studied (Fig. 2; Table 4). In accordance with our second hypothesis, the analysis of different elements by MPVD highlighted the existence of different patterns between elements in the effect of the phylogeny across the evolution of plants. For example, phylogeny explained a large proportion of the variance in C, Ca, Mg and S concentrations since early divergence times (c. 120–110 Ma), much earlier than the separation of orders and probably at the separation of the main groups of Angiosperms (Fig. 2). In N, the effect of phylogenetic differences became dominant c. 100 Ma, roughly at the time of diversification of different Angiosperm orders (Fig. 2). In metals such as Al and Mn, phylogenetic effects on elemental composition became dominant much later, c. 30 Ma, or were not significant (such as in Zn) (Fig. 2). By contrast, other elements, such as K and P, showed a more or less constant increase in the accumulated proportion of variance explained along the phylogeny (Fig. 2).

Table 4

Results of generalised linear models (GLM) with a gamma distribution and log link function showing the total variance explained (%) by the phylogeny (taxonomic family), affinity to gypsum (gypsovag/wide gypsophile/narrow gypsophile), soil conditions (principal component (PC1) after a principal components analysis (PCA) with soil variables, please refer to Supporting Information Methods S1 for further details) and climate (mean annual temperature, MAT, and mean annual precipitation, MAP).

Element	Family	Gypsum affinity	PC1 soil	MAT	MAP	Residuals
Al	29.81	0.16	0.01	0.58	1.70	67.74
C	56.06	1.72	0.74	0.17	2.48	38.83
Ca	46.13	2.25	1.28	0.07	0.37	49.90
K	56.74	0.30	0.52	0.19	0.10	42.15
Mg	48.35	5.95	0.22	0.08	0.81	44.59
Mn	33.42	2.39	0.35	1.01	0.46	62.37
N	62.60	0.05	1.04	0.11	0.61	35.59
Na	46.20	1.79	1.59	1.77	0.09	48.57
P	34.34	1.40	0.18	0.05	1.23	62.79
S	48.18	9.58	0.66	0.15	0.01	41.41
Zn	30.91	2.56	0.44	0.03	3.35	62.71

Residual variance not explained by our models is also shown.

In agreement with the results from PERMANOVA, affinity to gypsum soils was, in most cases, the second factor explaining a higher percentage of the variance when elements were analysed individually (Table 4). Its effect was significant for Mg and S (Table 5). The highest contribution to the variance across taxa was observed for S, in which affinity for gypsum soils accounted for c. 10% of the explained variance (Table 4). This factor also had a marked effect on C and Mg variability. The rest of factors analysed explained comparatively smaller proportions of the total variance (Table 4). However, MAP had a significant positive effect on plant Mg, S and N concentrations and negative effects on K and C contents (Table 5), and explained 3.4% of the variance of C concentrations (Table 4). MAT had a significant negative effect on K and P concentrations (Table 5), explaining < 0.6% of their variance (Table 4). Soil features (illustrated by the first PCA axis) had a minor effect on the variability of data (Table 4). Nevertheless, plants growing in soils with clayish and loamy textures, with higher K and lower Ca contents (i.e. higher PC1 values) showed higher K, Mg and N and lower Ca and P concentrations in their leaves (Table 5).

Table 5

Results of phylogenetically informed Bayesian linear model (generalised linear mixed models (GLMM) fitted with Markov chain Monte Carlo techniques (MCMC)) on the effect of specificity to gypsum soils on the elemental composition of plants.

	Gypsum affinity (gypsovag – narrow)		Gypsum affinity (gypsovag – wide)		PC1 (soil)		MAT		MAP
	CI	P	CI	P	CI	P	CI	P	CI	P
Al	−0.02, 0.02	0.621	−0.02, 0.02	0.842	−0.01, 0.03	0.182	−4 × 10−3, 4 × 10−3	0.95	−5 × 10−5, 1 × 104	0.429
C	−1.79, 0.50	0.249	−2.04, 0.74	0.337	−2.38, 0.18	0.094	−0.27, 0.21	0.802	−0.02, −4 × 10−3	0.002
Ca	−0.84, 0.67	0.896	−0.24, 1.62	0.108	−1.29, −0.24	< 0.001	−0.16, 0.06	0.367	−1 × 10−3, 4 × 10−3	0.327
Cu	−0.01, 0.02	0.966	−0.02, 0.02	0.982	−0.01, 0.01	0.884	−2 × 10−3, 2 × 10−3	0.998	−5 × 10−5, 3 × 10−5	0.908
K	−0.45, 0.13	0.291	−0.29, 0.40	0.655	0.06, 0.62	0.012	−0.16, −0.05	< 0.001	−4 × 10−3, −1 × 10−3	< 0.001
Li	−0.02, 0.02	0.996	−0.03, 0.03	0.964	−0.02, 0.02	0.928	−3 × 10−3, 3 × 10−3	0.978	−9 × 10−5, 9 × 10−5	0.978
Mg	−0.12, 0.29	0.421	0.17, 0.67	< 0.001 *	0.21, 0.53	< 0.001	−0.02, 0.05	0.35	2 × 10−3, 3 × 10−3	< 0.001
Mn	−0.01, 0.01	0.972	−0.02, 0.01	0.962	−0.01, 0.01	0.754	−2 × 10−3, 1 × 10−3	0.609	−4 × 10−5, −4 × 10−5	0.970
Mo	−0.11, 0.10	0.990	−0.19, 0.17	0.988	−0.08, 0.07	0.962	−0.03, 0.03	0.96	−4 × 10−4, 5 × 10−4	0.972
N	−0.24, 0.43	0.683	−0.40, 0.32	0.858	0.24, 0.89	0.002	−0.13, 0.01	0.09	2 × 10−3, 0.01	< 0.001
Na	−0.10, 0.07	0.756	−0.16, 0.10	0.635	−0.07, 0.10	0.77	−0.01, 0.02	0.733	−5 × 10−4, 3 × 10−4	0.589
P	−0.05, 0.01	0.13	−0.05, 0.03	0.806	−0.06, −0.01	0.024	−0.01, −8 × 10−4	0.016	−2 × 10−4, 4 × 10−5	0.166
S	−0.13, 0.90	0.112	0.90, 2.35	< 0.001 *	−0.04, 0.56	0.104	−0.09, 0.04	0.449	2 × 10−3, 0.01	< 0.001
Ti	−0.02, 0.02	0.880	−0.02, 0.02	0.976	−0.02, 0.02	0.92	−2 × 10−3, 3 × 10−3	0.978	−7 × 10−5, 8 × 10−5	0.962
Zn	−0.01, 0.01	0.986	−0.01, 0.02	0.902	−0.01, −0.01	0.834	−2 × 10−3, 1 × 10−3	0.876	−3 × 10−5, 4 × 10−5	0.776

Models included the first principal component (PC1) of a principal components analysis (PCA) of soil features, mean annual precipitation (MAP) and mean annual temperature (MAT) as covariates. In the results shown, gypsovags were considered the reference category for comparison of the gypsum affinity effect, however, models were also run with narrow gypsophiles as the reference category to verify potential significant differences between wide and narrow gypsophiles. When present, these differences were indicated by asterisks. The credibility interval (CI) for the probability that the estimate is higher than 0 is shown along with mcmcglmm P‐values. Significant effects are highlighted in bold. Negative or positive CI indicate the sign of factor effects. Model convergence was evaluated by running models several times. Please refer to Supporting Information Methods S1 for more details. n = 458 except for Cu (n = 376), Li (n = 183), Mo (n = 54), Ti (n = 278), N and C (n = 416).

Significant effects were also observed when wide and narrow gypsophiles were compared.

Ecological adaptation to gypsum shifts plant elemental composition

We detected a significant phylogenetic signal for the studied species when we analysed the main elements related to plant specialisation to gypsum soils, namely S (Bloomberg's K = 0.1, P = 0.001, Pagel's λ = 0.83, P < 0.001), Ca (Bloomberg's K = 0.08, P = 0.001, Pagel's λ = 0.59, P < 0.001) and Mg (Bloomberg's K = 0.04, P = 0.0026, Pagel's λ = 0.58, P < 0.001) (please refer to Fig. 3 for S and Figs S1, S2 for Ca and Mg, respectively). Foliar Ca accumulation was markedly lower in monocots and also in Lamiaceae and Asteraceae (Fig. S1). Nevertheless, a general trend was observed for most Iberian gypsum plants to accumulate Ca, as suggested by the recurrently high values shown across the phylogeny. In contrast to Ca, S and Mg showed more evident phylogenetic patterns. For instance, there was a pattern in Brassicales to accumulate more S in leaves, irrespective of plant affinity to gypsum (Fig. 3). The phylogenetic signal for S was also shown at the genus level (i.e. Ononis or Gypsophila). Some isolated species (such as the gypsophiles Campanula fastigiata or Vulpia gypsophila) showed markedly high S. By contrast, most monocots and the families Lamiaceae, Linaceae, Asteraceae and Apiaceae showed lower levels of foliar S accumulation. Mg showed also a clear phylogenetic pattern (Fig. S2), with the Caryophyllales displaying a strong phylogenetic signal of accumulation. Similarly, Fabaceae and Resedaceae also showed particularly high levels of Mg in their leaves.

Fig. 3

Leaf S accumulation traced across the phylogeny of the 83 investigated species. Circles of different colours indicate plant affinity to gypsum (gypsovag/narrow gypsophile/wide gypsophile). Tip values correspond to observed mean leaf S content (expressed as mg g−1 dry weight), which is reconstructed for ancestral values across the tree. Bloomberg's K = 0.1 (P‐value = 0.001), Pagel's λ = 0.83 (P < 0.001). [Colour figure can be viewed at wileyonlinelibrary.com]

Despite these general trends, Bayesian models showed that widely distributed gypsophiles had significantly higher foliar S and Mg concentrations compared with the rest, when controlling for the effect of phylogeny (Tables 2, 5). Narrow gypsophiles showed similar trends to gypsovags (Tables 2, 5). The foliar concentrations of the rest of elements showed no significant variation in relation to plant affinity to gypsum. These results persisted when the continuous ‘gypsophily index’ was included in the analyses to account for plant affinity to gypsum soils instead of the gypsophile wide, narrow and gypsovag categories (Methods S1; Table S4). When the elemental composition of plants was evaluated as a whole using RDA (Fig. 4), gypsum affinity was a significant factor explaining 4.3% of the variance. Wide gypsophiles showed higher foliar S, Mg and Ca concentrations compared with the rest of groups, gypsovags showed higher concentrations of C, N, P, K, Mn and Na in their leaves and narrow gypsophiles overlapped largely with gypsovags, reaching intermediate positions for elements such as S (Fig. 4).

Fig. 4

Redundancy analysis (RDA) biplot showing differences in the elemental composition of Iberian gypsum species with different affinity for gypsum soils (black, gypsovags; dark grey, narrow gypsophiles; light grey, wide gypsophiles). Results are based on Euclidean distances from Centre Log‐ratio coordinates (clr‐transformed data). Please refer to Materials and Methods for further details on the analyses.

The results of Mantel tests showed that the elemental composition of wide gypsophiles was less correlated with the phylogeny than that of narrow gypsophiles or gypsovags (Table 6). When all species were analysed together, we observed significant correlations between leaf elemental composition and the phylogenetic distance for elements such as C, Ca, K, Mn, Mg, N, Na, P or S and also for the elemental composition as a whole (Table 6). Concentrations of metals such as Al and Zn did not show any significant correlation with phylogenetic distances. These results generally held in gypsovags and narrow gypsum endemics, when species were grouped by their affinity to gypsum, but not in wide gypsophiles (Table 6). These results showed that the elemental composition of wide gypsophiles is less conditioned by the phylogeny than that of narrow gypsum endemics or gypsovags.

Table 6

Results of Mantel correlations (r M) and Mantel tests applied to a phylogenetic distance matrix and distance matrices of each chemical element concentrations.

Element(s) tested		All species		Gypsovag		Narrow		Wide
		r M	P	r M	P	r M	P	r M	P
Mean	Al	0.03	0.12	0.00	0.449	0.03	0.317	0.03	0.367
	C	0.09	< 0.001	0.11	0.036	0.14	0.036	0.10	0.158
	Ca	0.16	< 0.001	0.18	0.001	0.17	0.013	0.13	0.083
	K	0.09	0.005	0.12	0.018	0.16	0.018	0.00	0.557
	Mg	0.15	< 0.001	0.21	0.001	0.21	0.006	0.04	0.288
	Mn	0.11	< 0.001	0.15	0.006	0.13	0.048	0.02	0.368
	N	0.15	< 0.001	0.11	0.058	0.21	0.005	0.03	0.270
	Na	0.09	0.016	0.12	0.036	0.16	0.029	−0.12	0.913
	P	0.11	0.002	0.13	0.008	0.18	0.010	0.01	0.437
	S	0.11	0.003	0.09	0.063	0.17	0.015	0.02	0.375
	Zn	0.03	0.14	0.10	0.032	0.04	0.306	0.06	0.240
clr	All	0.27	< 0.001	0.34	< 0.001	0.34	< 0.001	0.10	0.130

For distance matrices accounting for more than one element, centred log‐ratio (clr) transformed data were used for Euclidean distance calculation. Significant effects are highlighted in bold. Results are shown for all species together and for groups of plants with different affinity for gypsum soils (gypsovags, narrow gypsophiles and wide gypsophiles), separately. Please refer to Supporting Information Methods S1 for further details on calculations.

Discussion

Phylogeny is a key determining factor of plant elemental composition but specialisation to atypical soils can also be critical

In agreement with our first hypothesis (please refer also to Broadley et al., 2004; Watanabe et al., 2007; Neugebauer et al., 2018), the elemental composition of Iberian gypsum plants was largely determined by its evolutionary history. However, our results indicated that the explanatory power of affinity to gypsum soils was also high, whereas the effect of climatic and soil factors was minor.

The phylogenetic effect on the whole elemental composition of plants varied depending on the divergence time used for phylogenetic grouping, being larger at shallower evolutionary times (up to 76% of the variation), mainly as a result of an increasing number of levels in phylogenetic clustering (as shown by our null models). These results agree with Yang et al. (2017) who found that species identity explained the largest proportion of the total variance in exchangeable cations concentrations. However, at divergence times between 32 and 5 Myr (Oligocene and Miocene), more than 40–60% of the variation (significantly above the null model) was still retained by phylogeny, a time period when diversification of most Mediterranean plant genera occurred (Postigo Mijarra et al., 2009; Vargas et al., 2018) and when most of the Iberian gypsum was formed (Escavy et al., 2012). These findings highlight the existence of recent evolutionary changes in the elemental composition of plants and may explain why Iberian gypsophiles are punctual tips in the phylogeny and not members of entire tribes, genera or families (Mota et al., 2011). Nevertheless, our results also showed a significant effect of ancient phylogenetic events on the elemental composition of species, highlighting the importance of accounting for phylogeny at different time scales (Buira et al., 2021). Divergence times as old as 120 Myr (c. Order diversification) accounted for 25% of the elemental variation in study species. This ancestral phylogenetic signal suggests the presence of adaptive constraints in some lineages and points to a strong conservation of certain metabolic processes acquired during early stages of the evolution of Angiosperms (Watanabe et al., 2007). Current data do not allow ascertaining if gypsum or other factors acted as selective forces at these early divergence times. However, elements such as S, Mg and Ca show a significant increase in the variation explained by phylogeny at 120 Myr, being the three elements critical for plant adaptation to gypsum soils. Dominant in the surface or not, gypsum existed on Earth far before that time. Although most Iberian gypsum outcrops are relatively recent, a significant fraction (30%) are from the Mesozoic (250 Ma) or older (Escavy et al., 2012). It, therefore, seems feasible that a first adaptive evolution to gypsum of certain dicotyledonous lineages occurred c. 120 Ma on gypsum materials. In any case, we cannot dismiss that, alternatively, adaptation to other factors (for example calcareous substrates) could have shaped the leaf elemental composition of plants, particularly towards the accumulation of exchangeable cations, which could have later served as useful pre‐adaptations to adapt to gypsum. Any of these processes might help to understand the unbalanced occurrence of gypsophiles in favour of certain families or orders (Moore et al., 2014). The geologic history of gypsum outcrops in the Iberian Peninsula is, therefore, compatible with a recent convergent evolution of gypsophilic lineages from certain families and orders bearing ancient adaptations or pre‐adaptations, developed in early stages of Angiosperm evolution, to cope with high Ca, Mg and S in the soil (Moore et al., 2014).

In agreement with our second hypothesis, different patterns were observed for phylogenetic effects among elements. Our results showed high phylogenetic effects for C, Ca, K, Mg, N, Na and S. Several previous studies have reported strong phylogenetic effects for most of these elements (Broadley et al., 2003; Kerkhoff et al., 2006; Willey & Wilkins, 2006; Watanabe et al., 2007; White et al., 2015, 2017). Particularly well known are the strong phylogenetic signals of Ca and Mg (Broadley et al., 2003; White et al., 2015), which our MPVD traced as earlier than order diversification, probably at the divergence or large Angiosperm groups (Fig. 2). Accordingly, eudicots showed higher Ca content than monocots (Broadley et al., 2003). We also detected higher Ca concentrations in members of Caryophyllales (marked in Frankenia thymifolia and Gypsophila sp.), Brassicaceae and Ononis sp., and lower concentrations in Asteraceae and Lamiaceae (Fig. S1). These results agree with the putative ability of Caryophyllales and Brassicaceae to accumulate Ca, whereas Asteraceae tended to show lower Ca concentrations (Broadley et al., 2003). The notable accumulation of Mg observed in Caryophyllaceae (Fig. S2) is consistent with the ability of Caryophyllales to accumulate Mg in the cell vacuole (White et al., 2018). The phylogenetic signal for Na, with an abrupt change at the order–family level, agrees with the marked differences previously reported for Caryophyllales (White et al., 2017). Similarly, the strong phylogenetic signal detected for S, with higher accumulation in Brassicales (Fig. 3), is consistent with the high S requirements of certain lineages such as Brassicaceae (Willey & Wilkins, 2006; Neugebauer et al., 2018). The marked accumulation of S observed in isolated clades such as Gypsophila sp. and Ononis or isolated taxa such as V. gypsophila (Fig. 3), however, may be indicative of a recent adaptation to soils rich in gypsum, and requires further study. In agreement with previous studies (Zhang et al., 2012; Zhao et al., 2016, but please refer to Metali et al., 2012), the concentrations of some metals such as Al, Mn and Zn were less phylogenetic structured.

Although phylogenetic effects were strong, our results highlighted the important effect of ecological adaptation to atypical soils on the elemental composition of plants. Affinity to gypsum soils accounted for up to 6% of the total variance in the elemental composition of study species, an effect that increased up to 10% for S, the element identified as most discriminating for gypsum plants (Merlo et al., 2019). The rest of ecological and environmental factors analysed, although significant, showed a much smaller impact on the elemental concentrations of Iberian gypsum plants. MAP and MAT were the environmental factors that had the largest effect. However, climatic factors explained only c. 1.7% of the total variance of our elemental data, which is markedly lower than the effect reported in previous studies (Han et al., 2011; Zhang et al., 2012). As our study is restricted to Iberian gypsum ecosystems, these discrepancies can be explained by the much lower variation in climatic conditions compared with those in previous studies (cf. Han et al., 2011). Noteworthy, the restrictive role of gypsum soils on vegetation only appears under arid and semiarid conditions (Escudero et al., 2015), limiting the potential role of climate on the elemental composition of gypsum plants. The small effect observed for soil properties on the elemental composition of Iberian gypsum plants contrasts with the results of previous studies in which soil had a marked effect on the elemental composition of plants (Ågren & Weih, 2012; Shtangeeva et al., 2020). These discrepancies cannot be explained by a low resolution in our dataset, but are likely to be due to the low variability in the soil conditions included in the study. All soils in our analyses had high gypsum content, with limited variation in terms of nutrient content, texture, pH, Ca or S concentrations. Furthermore, the low effect of soil conditions on the elemental composition of plants may be due to the consideration of total soil elemental composition rather than extractable or phyto‐available fractions. Stein et al. (2017) detected a poor correlation among total elemental concentrations in soil and plants, whereas relationships were markedly improved if the exchangeable fraction of soil elements was considered.

Foliar concentrations of macronutrients or elements that are required by plants in larger amounts (such as C, N and P) were less variable than microelements (such as Na, Cu, Li, Ti and Al; Han et al., 2011; Bai et al., 2019). Nevertheless, S and Mg, both macronutrients, showed relatively large variability when all plants were analysed together, probably due to the different ability to accumulate S and Mg in plants with a contrasting affinity to gypsum. For example, wide gypsophiles such as Gypsophila tomentosa accumulated almost two orders of magnitude more S and Mg than the average values for gypsovags. Interestingly, the variability in S decreased when wide gypsophiles were analysed separately, probably as a result of the convergent evolution towards S accumulation in plants specialised to gypsum soils (please refer to subsequent paragraphs).

Ecological adaptation to gypsum soils shifted leaf elemental composition of Iberian gypsophiles towards increased foliar S and Mg accumulation

Our results, independently of the approach, showed that plant specialisation to gypsum soils is related to an increased ability of plants to accumulate higher concentrations of S and Mg in their leaves. Muller et al. (2017) also detected clear patterns of S and Ca accumulation that were independent of the phylogeny in widely distributed gypsophiles from the Chihuahuan desert. The marked accumulation of S, Mg and Ca in widely distributed Iberian gypsophiles is well documented (Duvigneaud & Denayer‐De Smet, 1968; Palacio et al., 2007; Merlo et al., 2019), but most of these previous studies did not take phylogenetic effects into account. Although average Ca concentrations were markedly high in wide gypsophiles, the effect of plant affinity to gypsum on foliar Ca concentrations was not significant, contrasting with the effects on S and Mg. The lack of significant differences in the Ca concentrations of Iberian gypsum plants when phylogeny is accounted for, agrees with the hypothesis that gypsophilic floras may have evolved from calcicole lineages (Heiden et al., 2022), with a preadaptation to tolerate high Ca concentrations in the soil partly through Ca accumulation. This notion seems highly plausible for the Iberian Peninsula, due to the relatively high Ca concentrations of most taxa (Fig. S1), the intermingled distribution of gypsum and calcareous deposits (Rodríguez‐Fernández, 2004), their close physical and chemical characteristics (FAO, 1990) and the widespread occurrence of calcicole gypsovags (Braun‐Blanquet & Bolòs, 1957). In addition to Ca, Mg and S, Palacio et al. (2007) also reported a marked accumulation of N and P in wide gypsophiles, a pattern that was not confirmed by our results and could be attributed to unaccounted phylogenetic effects.

Phylogeny determined the baseline elemental concentrations of plants, but specialisation to gypsum significantly increased S and Mg concentrations of wide gypsophiles compared with their gypsovag and narrow gypsophile relatives. Our results, consequently, point to a convergence in the accumulation of S and Mg in wide gypsophiles, not observed in gypsovags or narrow gypsophiles. This convergence is exemplified by the marked S accumulation in species belonging to lineages with reported moderate S concentrations (Willey & Wilkins, 2006), such as C. fastigiata (Campanulaceae), Ononis tridentata (Fabaceae) or Gypsophila sp. (Caryophyllaceae; Fig. 3).

According to our results, S and Mg accumulation in Iberian gypsum specialists seems to have evolved several times along the evolution of Angiosperms, involving different possible ‘metabolic solutions’. For example, widespread gypsophilic lineages from Iberia and the Chihuahuan desert show gypsum crystals in their leaves (Palacio et al., 2014; Muller et al., 2017), suggesting that the accumulation of S as gypsum inside the cell vacuole could be de‐toxifying mechanisms for wide gypsophiles. Nevertheless, some lineages with high S concentrations that were extensively represented amongst gypsophiles, such as the Brassicaceae, do not accumulate gypsum crystals (Palacio et al., 2014; Muller et al., 2017). These species could accumulate excess S into organic molecules, such as glucosinolates, thiosulfinates (typical of Allium sp.) or flavone sulphates (Ernst, 1998). In a preliminary analysis of the glucosinolate content of Brassicaceae with different affinity for gypsum soils, Tuominen et al. (2019) detected higher glucosinolate content in coexisting wide gypsophile lineages compared with gypsovags. Mechanisms for Mg accumulation in gypsum plants remain unexplored. Most free Mg in plants accumulates in the cell vacuole, so that cytosolic Mg levels are kept low (Waters, 2011). Plants growing in soils with a high S supply form magnesium sulphate crystals in their tissues (He et al., 2014), which might be a feasible Mg‐accumulating option for gypsophiles showing gypsum biomineralisation.

Our study does not permit a conclusion on the adaptive values of S and Mg accumulation in gypsum specialist plants. Ascertaining if the accumulation of these elements offers a selective advantage to gypsum specialists, or if it is merely a by‐product of other adaptations to succeed on gypseous soils, would require estimates on the fitness of plants, which is out of the scope of the present study. Romao & Escudero (2005) reported a lower survival and growth of the wide gypsophile Helianthemum squamatum when cultivated out of gypsum. More recently, Cera et al. (2021) showed that gypsophiles cultivated on gypsum and limestone soils completed their life cycle with no differences in plant fitness, but wide gypsophiles maintained a physiological dependence on S and Mg, evidenced by increased accumulation, irrespective of the substrate. These results point to a physiological dependence of wide gypsophiles on S and Mg accumulation, which deserves further study.

Conclusions

Phylogeny is a strong determinant of the elemental composition of Iberian gypsum plants. Markedly, the MPVD procedure highlighted that this effect involved processes that occurred at two very different divergence times. Recent evolutionary events played a prominent role in shaping actual plant elemental composition, but ancient relationships at the base of the phylogeny also explain a marked proportion of the total variation in the elemental concentrations of these species. Interestingly, phylogenetic effects were different for different chemical elements across the evolution of study species. On top of the phylogenetic background, our results highlight the relevance of ecological adaptations on the elemental composition of plants. Adaptation to gypsum soils shifted plant elemental composition towards convergent foliar S and Mg accumulation in Iberian widespread gypsophilic lineages. These results are of paramount importance for understanding the evolution of plant life on gypsum and might serve as an example for further research in other atypical substrates.

Author contributions

SP and GM‐M designed the study. GM‐M, SP, PT, AC, JFM, MEM, ES‐S, FM‐H, AJM‐F, AMS and FJP‐G conducted fieldwork. MP‐SS ran chemical analyses. SP and PT analysed the data with input from AC, ALL, AMS and AE. All authors discussed the results. SP and PT wrote the manuscript with contributions from all authors.

Fig. S1 Leaf Ca accumulation traced across the phylogeny of the 83 investigated species indicating their affinity to gypsum.

Fig. S2 Leaf Mg accumulation traced across the phylogeny of the 83 investigated species indicating their affinity to gypsum.

Methods S1 Details on methodology for statistical analyses.

Table S1 List of the 50 species included in the initial tree constructed with v.phylomaker.

Table S2 Scores of different variables plus eigenvalues and their contribution to the three first principal components analysis components summarising the variability in soil features measured.

Table S3 Results of multilevel pairwise comparisons among groups of plants with different affinity for gypsum soils after PERMANOVA.

Table S4 Results of phylogenetically informed Bayesian linear models on the effect of gypsophily index on the elemental composition of plants.

Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.

Click here for additional data file.