This is the accepted version of the article: Fernández-Martínez, Marcos; Pearse, I.; Sardans i Galobart, Jordi; [et al.]. «Nutrient scarcity as a selective pressure for mast seeding». Nature Plants, Vol. 5 (Dec. 2019), p. 1222–1228. DOI 10.1038/s41477-019-0549-y This version is avaible at https://ddd.uab.cat/record/216950 under the terms of the license 1 Nutrient scarcity as a selective pressure for mast seeding *1 2 3,4 5,6 7 2 Authors: M. Fernández-Martínez , I. Pearse , J. Sardans , F. Sayol , W.D. Koenig , J.M. 8 9 10,11 12 13 4 3 LaMontagne , M. Bogdziewicz , A. Collalti , A. Hacket-Pain , G. Vacchiano , J.M. Espelta , 3,4 1 4 J. Peñuelas , I.A. Janssens 5 Affiliations: 1 6 Centre of Excellence PLECO (Plants and Ecosystems), Department of Biology, University of 7 Antwerp, 2610 Wilrijk, Belgium. 2 8 U.S. Geological Survey, Fort Collins Science Center, Fort Collins, Colorado 80526 USA 3 9 CSIC, Global Ecology Unit, CREAF-CSIC-UAB, Bellaterra 08193 Barcelona, Catalonia, Spain 4 10 CREAF, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, Spain 5 11 Department of Biological and Environmental Sciences, University of Gothenburg, SE 405 30 12 Gothenburg, Sweden 6 13 Gothenburg Global Biodiversity Centre, SE 405 30 Gothenburg, Sweden 7 14 Lab of Ornithology, Cornell University, Ithaca, NY, USA 8 15 Department of Biological Sciences, DePaul University, Chicago, IL, USA 9 16 Department of Systematic Zoology, Faculty of Biology, Adam Mickiewicz University in Poznań, 17 Poland 10 18 Division on Impacts on Agriculture, Forests and Ecosystem Services, Foundation Euro- 19 Mediterranean Centre on Climate Change (CMCC-IAFES), 01100, Viterbo, Italy 11 20 Institute for Agriculture and Forestry Systems in the Mediterranean, National Research 21 Council of Italy (CNR-ISAFOM), 87036, Rende (CS), Italy 12 22 Department of Geography and Planning, School of Environmental Sciences,University of 23 Liverpool, Liverpool, UK 13 24 Università di Milano, DISAA. Via Celoria 2, 20133, Milan - Italy 25 * Correspondence to: M. Fernández-Martínez, marcos.fernandez-martinez@uantwerpen.be 26 Length of the main text: 2495 words 27 Length of the methods: 2356 words 28 Figures and tables: 3 figures, 1 table 29 Length of the captions: 387 words 30 References: 55 (+17 from Methods) 31 Introductory paragraph 32 Mast seeding is one of the most intriguing reproductive traits in nature. Despite its 33 potential drawbacks in terms of fitness, the widespread existence of this phenomenon 34 suggests that it should have evolutionary advantages under certain circumstances. 35 Using a global dataset of seed production time series for 219 plant species from all the 36 continents, we tested whether masting behaviour appears predominantly in species 37 with low foliar N and P concentrations, when controlling for local climate and 38 productivity. Here we show that masting intensity is higher in species with low foliar N 39 and P concentrations and especially imbalanced N:P ratios, and that the evolutionary 40 history of masting behaviour has been linked to that of nutrient economy. Our results 41 support the hypothesis that masting is stronger in species growing under limiting 42 conditions and suggest that this reproductive behaviour might have evolved as an 43 adaptation to nutrient limitations and imbalances. 44 Main text 45 Mast seeding, often called masting, has long intrigued biologists as one of the most 46 bizarre reproductive behaviours found in nature1,2. This behaviour consists of the 47 synchronous production of highly variable seed crops over time3. Masting has often 48 been considered an evolutionary paradox because organisms that skip reproductive 49 attempts should have lower fitness than those that reproduce at every opportunity4. 50 Nonetheless, the fact that this reproductive behaviour is found in different lineages 51 suggests that masting behaviour should be beneficial, at least under certain scenarios5. 52 The most widely accepted hypotheses explaining the selective advantages of masting 53 are all related to economies of scale6,7. Briefly, these hypotheses state that, in terms of 54 fitness, it is more efficient for plants to produce large number of seeds every few to 55 several years than producing a constant amount every year. This general mechanism 56 includes the predator satiation hypothesis2,8–10, where predators are starved during 57 years of null or low reproduction and satiated during high reproduction mast years, 58 leaving large numbers of seeds intact. Another example is the pollination efficiency 59 hypothesis6,11,12, stating that, particularly for wind-pollinated plants, saturating the 60 atmosphere with pollen in a given year is more efficient than producing regular 61 amounts of pollen each year in order to ensure pollination. Given that masting is 62 present in only a modest percentage of plant species13, such economies of scale are 63 apparently advantageous only under certain circumstances. What those circumstances 64 are remains, so far, under debate. 65 The environmental stress hypothesis14 suggests that masting behaviour should be 66 stronger under unfavourable growing conditions or limitation of resources, conditions 67 under which economies of scale should be more beneficial3,12,15. This is because plants 68 growing in unfavourable environments presumably experience more difficulties in 69 acquiring the required resources to reproduce, as suggested by the resource 70 accumulation hypothesis16,17. According to this hypothesis, plants growing under 71 favourable conditions will be able to accumulate the required amount of resources 72 every year and, therefore, present a regular pattern in seed production, without 73 exhibiting any underlying negative temporal autocorrelation that could indicate resource 74 depletion after reproduction16. The opposite is true for plants growing in unfavourable 75 conditions, which will exhibit high interannual variability and negative temporal 76 autocorrelation in seed production due to potential resource depletion after seeding. 77 However, there is no current empirical evidence suggesting that species with higher 78 interannual variability in fruit production are more likely to exhibit negative temporal 79 autocorrelation than species that produce seeds more regularly5. On the other hand, 80 weather variability has been found to be a key factor driving interannual variability in 81 fruit production in many plant species12,18–21. Temporal patterns in weather events (i.e., 82 temporal variability and autocorrelation) could, therefore, potentially shape the temporal 83 patterns of fruit production22. 84 Foliar nutrient concentrations play a key role in plant ecophysiology and ecosystem 85 functioning: photosynthetic rates are linked to foliar nitrogen (N) and phosphorus (P) 86 concentrations23–25. Together with carbon, they are the basis of ecological 87 stoichiometry26,27 and are fundamental parts of the elementome or the biogeochemical 88 niche28, useful for inferring ecological traits from the elemental composition of 89 organisms29. Carbon (C), and especially N and P have been suggested to be potential 90 resources determining seed production and masting behaviour15,30–32 because seeds 91 and fruits are enriched with N and P compared to vegetative tissues33. Low foliar 92 concentrations of N and P would imply lower photosynthetic rates23 that would in turn 93 result in reduced accumulation of C, in addition to low accumulation of N and P. Thus, 94 based on the resource budget model, plant species with low foliar N and P 95 concentrations would be expected to exhibit stronger masting behaviour in terms of 96 variability and negative temporal autocorrelation. 97 Both foliar N and P have been reported to be positively correlated with investment in 98 plant reproduction34–36. Additionally, the interaction between foliar N and P, (i.e., the 99 N:P ratio) has been suggested to be an important indicator of plant nutritional status, to 100 promote vegetative productivity37–40 and to relate to reproductive effort36 and ecological 101 strategies28. The lower N:P ratio of seeds compared to vegetative tissues suggest that 102 nutrient imbalances may constrain reproduction more than growth38. Previous studies 103 have reported that plants and shoots producing flowers often exhibit lower N:P values 104 than those that do not38,41, highlighting the role of nutrient imbalances in plant 105 reproduction. However, no study has yet been focused on the evolutionary linkage 106 between plant reproduction and foliar nutrient stoichiometry over a large set of species. 107 Here we hypothesize that mast seeding behaviour evolved predominantly in plant 108 species with low foliar N and P concentrations and high N:P ratios. To test this 109 hypothesis, we first explored the relationship between temporal variability and 110 autocorrelation in reproductive effort using a global dataset of 219 plant species from 111 all the continents to see if higher variability in seed production is related to higher 112 potential resource depletion (negative temporal autocorrelation)5. We then fitted 113 phylogenetic linear models to test whether masting intensity, here defined as the 114 combination of temporal variability and autocorrelation using the consecutive disparity 115 (D) index42,43, was related to foliar N, P and N:P ratio, while controlling for the 116 evolutionary history of the species and other potentially influencing factors such as 117 local climate and productivity of the regions where species were sampled. Finally, we 118 tested the model of evolution of masting behaviour by means of ancestral character 119 reconstructions and explored whether it has evolved in concert with foliar N and P 120 economies. 121 Results 122 Redefining masting intensity: temporal variability, potential resource depletion and their 123 evolutionary relationships 124 Our results indicated that temporal variability (here calculated as the proportional 125 variability index44 PV, see Methods) in seed production is evolutionary conserved 126 (Figure S1a), showing a strong phylogenetic signal45. In contrast, lag-1 temporal 127 autocorrelation (AR1, indicating potential resource depletion when it takes negative 128 values) is not phylogenetically conserved. Additionally, temporal variability and 129 autocorrelation are not evolutionary correlated. Hence, species exhibiting higher 130 temporal variability do not necessarily exhibit any particular pattern of temporal 131 autocorrelation (Figure S1b). Masting intensity, defined as D in this study (see 132 Methods), accounted for both features of masting behaviour (Figure 1), temporal 133 variability and potential resource depletion (negative AR1 coefficients), hence defining 134 masting behaviour more broadly than the coefficient of variation (CV) alone, as it is 135 typically assessed42. Masting intensity was also preserved phylogenetically (Figure 1). 136 Our results suggest that most species exhibit low or intermediate masting intensity, 137 while only a few exhibit strong masting behaviour (see histogram in Figure 1a). 138 Estimating controls and the mode of evolution of masting behaviour 139 Phylogenetic models indicated a statistically significant negative interaction between 140 foliar N and P explaining variability in masting intensity across species (Table 1, Figure 141 2). Model results indicated that when foliar P is low, increasing foliar N increases 142 masting intensity as N:P increases. The inverse situation (high masting intensity at very 143 low N:P ratios) is also possible according to model results, despite being a less likely 144 scenario (we only had 19 species with N:P < 8 in our dataset). The highest masting 145 intensity was found in species with high N:P while the lowest was found in species with 146 high foliar N and P concentrations. Species with low foliar N and P concentrations 147 showed intermediate values of masting intensity. Temporal variability of the local 148 climate or productivity of the region from where species were sampled did not explain 149 variation in masting intensity. Mean annual precipitation, however, was positively 150 related to masting intensity (Table 1). Temporal variability in seed production was 151 negatively related to mean annual temperature and positively related to temporal 152 variability of annual precipitation. Additionally, species with higher foliar P exhibited 153 lower temporal variability. When assessing temporal variability using the CV, we found 154 that species from regions with higher temporal variability in productivity also presented 155 higher CV of seed production (β=0.22 ± 0.08, P=0.004, λ=0.58). Potential resource 156 depletion was more likely to occur in species living in climates exhibiting negative 157 temporal autocorrelation in annual precipitation (Table 1). Interestingly, the model also 158 showed a positive interaction between foliar N and P, indicating that potential resource 159 depletion after seeding was more likely to occur in species with low foliar N and P. 160 Our analyses also revealed the evolutionary processes linking masting behaviour and 161 foliar nutrient concentrations (Table S1). We found that adaptive processes - Ornstein- 162 Uhlenbeck (OU) models based on ancestral character reconstructions (Figure S2) - 163 are the models that better explain current patterns in masting intensity, linked to foliar 164 nutrient concentrations. Adaptive OU models explaining evolutionary history of masting 165 intensity with N, P and N×P (a measure of total N and P foliar concentration) assumed 166 that species evolved toward different optimum values depending on whether they 167 present high or low masting intensity (masting and non-masting species) (OUM model). 168 However, a single optimum for all species was equally likely for foliar N, P and N×P 169 (less than 2 units of AICc between models were found, OU1, see Methods). 170 Nonetheless, given that OUM models with different optimum values performed equally 171 well for foliar N, P and N×P (Table S1), OUM models are used to illustrate our results. 172 The model that best fitted the data for foliar N:P was an adaptive OU model with 173 different optimum and phenotypic variation for high and low masting intensity species 174 (OUMV model). 175 Results from OU models indicated that masting species had, on average, 9.5 ± 0.4% 176 and 18.3 ± 0.5% lower foliar N and P concentrations, respectively, than non-masting 177 species (Figure 3, Table S2). Foliar N:P ratio was 11.4 ± 0.2% (mean ± standard error 178 of the mean) higher in masting species, while N×P, combined availability of N and P, 179 was 28.7 ± 1.0% lower than in non-masting species. These results were consistent 180 when using five different thresholds of masting intensity (see Methods) to classify 181 species as masting or non-masting species (Table S2). Like masting intensity, both 182 foliar N and P concentrations were preserved throughout the phylogeny (Figure S3). 183 Discussion 184 The role of foliar nutrients in seed production 185 Our results suggest that masting intensity co-evolved with species-specific optimal 186 foliar N and P concentrations and that species with lower N and especially lower P, and 187 therefore high N:P, exhibit higher masting intensity (Figure 3). Hence, our analyses 188 supported our initial hypothesis stating that masting behaviour evolved predominantly 189 in plant species with low foliar N and P concentrations and high N:P ratios. This 190 observed evolutionary pattern may have originated because of different underlying 191 mechanisms driven by environmental and physiological constraints. One potential 192 mechanism explaining these findings could be the physiological role of foliar N and P 193 concentrations in plants. Foliar N is well known for being the primary limiting nutrient for 194 reproduction34,35 and vegetative growth46,47. Elevated foliar P concentrations, and low 195 N:P ratios, have been shown to allow larger seed crops in multiple species15,36,41,48, 196 coinciding with fast growing species according to the growth rate hypothesis26. P is also 197 essential to maintain water-use efficiency and growth, particularly during drought 198 conditions49. Both nutrients are essential elements of ribosomes and, therefore, play a 199 major role in organismal metabolism50. The production of seeds, which are structures 200 enriched with N and P33, may potentially benefit from high concentrations of nutrients 201 even more than vegetative tissues. 202 Low foliar N and P concentrations and high N:P imbalance are often indicative of 203 unfavourable environments for plant growth15,51,52, such as infertile soils. These poor 204 growing conditions are often related to dry or cold climates where decomposition of 205 organic matter is constrained, rates of nutrient mineralisation and weathering are 206 low49,53, or disturbances (e.g., wildfires) are frequent54. However, differences in foliar N 207 and P concentrations may not always reflect differences in nutrient availability across 208 species, as coexisting species may have different elemental compositions (the 209 biogeochemical niche hypothesis28). Nonetheless, plants with high foliar N and P 210 concentrations must either be growing in environments without nutrient limitations or 211 have developed mechanisms that allow them to sustain high foliar nutrient 212 concentrations even if they are growing under unfavourable conditions. In any case, 213 these higher concentrations of nutrients should confer a competitive advantage in 214 terms of C acquisition, because higher N and P concentrations are, on average, linked 215 to higher photosynthetic rates23. 216 Evolutionary history of masting intensity and foliar nutrients 217 Being capable to invest more in reproduction does not explain masting behaviour by 218 itself, because equal average seed crops over time could lead to different reproductive 219 behaviours15. The necessary link between the ecological stoichiometry and masting 220 theoretical backgrounds lies in the resource accumulation hypothesis16,17. Plant species 221 with lower or imbalanced availability of N and P may present more difficulties in 222 acquiring the necessary amount of C, N and P to successfully produce seeds regularly, 223 thereby mechanistically producing a reproductive behaviour aligned with high masting 224 intensity: high interannual variability and negative temporal autocorrelation (i.e., 225 potential resource depletion). The combination of low and imbalanced nutrient 226 availability, causing high variability in seed production and potential resource depletion 227 after crops, and environmental variability that synchronises the reproduction of 228 individuals through weather cues7,19 would finally shape the reproductive behaviour of 229 masting species. In fact, our results also revealed that mean annual precipitation, and 230 its temporal variability and autocorrelation, are related to the reproductive behaviour of 231 plant species (Table 1). These results indicate that even nutrient-rich species can 232 present masting behaviour if they grow under climates with highly fluctuating weather 233 conditions. Therefore, we suggest that the interaction between weather conditions and 234 the availability of nutrients, both conditioning photosynthetic rates, are the triggering 235 factor of the common nonlinear (often exponential-like) response22 between seed 236 production and weather variability in masting species8,18,21. 237 As a result of the environmental and physiological constraints, species may have been 238 selected to exhibit distinctive reproductive behaviours in order to increase their fitness. 239 Nutrient-rich species may not have developed a resource-conservative masting-like 240 reproduction strategy because of their capacity to produce abundant seeds with 241 regularity, avoiding losing reproduction attempts. Instead, because the investment in 242 reproduction in terms of C and nutrients should be proportionally more expensive for 243 nutrient-poor or nutrient-imbalanced species, the pressure to exhibit more cost-efficient 244 reproduction may have selected such species to produce fewer but larger reproductive 245 events in order to take advantage of one or more economies of scale6,7,11. As a side 246 effect of these massive reproductive events, negative temporal autocorrelation would 247 also have appeared in masting species because of potential resource depletion 248 (Figure S1). Hence, these particular traits would have been preserved throughout 249 evolutionary history (Figure 1) because foliar functional traits and masting intensity 250 have co-evolved (Figure 3, Table S2). Nonetheless, our results do not discard other 251 potential selective pressures that may have triggered the evolution of masting 252 behaviour. Some species may have been selected to mast to improve their pollination 253 efficiency55, to escape seed predation from voracious predators8,9 or because of 254 interspecific competition in different ways (e.g., seedling establishment). Also, a 255 selective pressure towards reproducing more constantly could happen in animal- 256 pollinated species, where a more constant production of flowers would favour 257 populations of pollinators and, hence, pollination. Further research on reproductive 258 behaviour of early plants, such as bryophytes, and taxa from different realms using 259 concepts from the masting literature may facilitate better understanding of the evolution 260 of different reproductive behaviours in nature. 261 Methods 262 Datasets 263 Masting database 264 Data on interannual reproductive effort (seed or fruit production) were compiled from 265 Web of Science searches, scanning the literature cited of published papers to look for 266 more records of reproductive effort, contacting managers of wildlife surveys, forestry 267 districts, and regional seed surveys, and soliciting datasets in the Ecolog listserv 268 (https://listserv.umd.edu/archives/ecolog-l.html). See Pearse et al.14 for more 269 information on data collection methods and characteristics of the dataset. 270 We only included records when: i) data were available for more than four consecutive 271 years per species at a given site (with clear geographical coordinates), ii) records could 272 be clearly assigned to plants of a particular species, iii) records were not measured in 273 such an indirect way their accuracy could be jeopardised (e.g., anticipated correlations 274 with gamete abundance), iv) data represented seed or fruit production, or inflorescence 275 set only for those cases where inflorescences are strongly linked to seed or fruit 276 production, and v) records from iteroparous perennial plants whose seed set could not 277 be explained by changes in population size. We did not include records of pollen 278 production nor from records from agricultural settings. 279 Overall, the reproductive effort dataset contained 1084 records of reproductive effort 280 including 363 plant species (trees, shrubs, vines, grasses and herbs) from 205 studies, 281 ranging from 1900 to 2014 and covering the six majorly vegetated continents (Africa 282 (17), North America (466), Europe (280), Japan (68), New Zealand (67), Central 283 America (118) and other regions (68)). On average, records were 11.9 years long, 284 although 131 studies had more than 20 years of data14. 285 For each site and species, we calculated the proportional variability index in seed 286 production (PV)44 as a measure of temporal variability. The PV index overcomes 287 several statistical and mathematical issues of the CV index44. The PV index was 288 calculated as: ∑ ( ) 289 Where z is calculated as: ( ) ( ) 290 Where ―z‖ represents the list of individual values (seed production per year) from which 291 to calculate the pairwise comparisons and ―n‖ indicates the number of values in a 292 variable. We also calculated the lag-1 temporal autocorrelation (AR1) as a measure of 293 potential resource depletion using the acf function in R, and the consecutive disparity 294 index (D)42,43 as a measure of masting intensity. We defined masting intensity as D 295 because, like the PV index, this index also overcomes several statistical and 296 mathematical limitations of the coefficient of variation (CV = standard deviation × mean- 1 297 ), including dependence on the mean, dependence on the length of the time series, 298 and bias associated with non-Gaussian data commonly used to describe masting13,14. 299 More importantly, we used D because it combines two of the main features describing 300 masting behaviour: temporal variability and lag-1 autocorrelation (AR1) in seed 301 production42. D is defined as: ∑| | 302 Where pi is the series value (seed production in our case) at time i, n is the series 303 length and k is a constant (often one unit) to avoid numerical indetermination in time 304 series with 0 values. D is high when temporal variability is high and lag-1 temporal 305 autocorrelation is negative (i.e., strong masting behaviour, showing potential resource 306 depletion after large seed crops20,31). Conversely, D is low when temporal variability is 307 low and temporal autocorrelation is close to zero or positive (i.e., describing a constant 308 pattern of seed production and hence, no masting behaviour). Given that the CV is still 309 the most widely used index to assess temporal variability in masting studies, we also 310 calculated the CV of seed production per species and site to compare its results with 311 those shown by the PV index. When multiple records from the same species were 312 available, we calculated the average masting intensity (D), temporal variability (PV and 313 CV) and the potential resource depletion coefficient (AR1) per species. 314 Climate, productivity and foliar nutrients 315 We extracted mean annual temperature (MAT, °C) and precipitation (MAP, mm yr-1) 316 and their seasonality (MATs [standard deviation of monthly values × 100] and MAPs 317 [CV]) for each location in our masting database from the WorldClim2 database56. The 318 climate database contains long-term means (1950-2000), calculated on a 30 arc- 319 second grid. Data for variability and autocorrelation for temperature and precipitation 320 was extracted from the CRU TS v3.25 dataset57. To estimate site productivity we used 321 a remotely-sensed gross primary production (GPP) database58. For annual 322 temperature and precipitation of the CRU TS and the GPP global databases we 323 calculated D, PV and AR1 indices for each pixel. We then extracted climate (i.e., MATD, 324 MAPD, MATPV, MAPPV, MATAR1 and MAPAR1) and productivity (i.e., GPPD, GPPPV and 325 GPPAR1) data for each site and species in our masting database and calculated the 326 average per species. We used these variables in our statistical analyses to control for 327 site-specific differences in temporal variability and autocorrelation patterns of climate 328 and productivity. 329 Data for foliar concentration of N and P for the species in our database were gathered 330 from the TRY trait database (https://www.try-db.org/TryWeb/Home.php)59, the BIEN 331 database (http://bien.nceas.ucsb.edu/bien/)60, the ICP Forests database on foliar 332 elemental concentration61 (www.icp-forests.net), and the Catalan Forest Inventory 333 (http://www.creaf.uab.cat/iefc/). To estimate an average value per species, we first 334 calculated the mean value per species and database. We then merged all databases 335 and calculated the mean values per species. Species names in our database were 336 checked and corrected using The Plant List database in the R package Taxonstand62. 337 Phylogenetic analyses were performed using the plant phylogeny provided by Qian & 338 Jin (2016). Out of the 363 species in the masting database, 219 species (~60%) names 339 matched those in the phylogenetic tree, and therefore, only those were used for further 340 analyses. Analyses using foliar N and P data were restricted to the 168 species (~46%) 341 for which we could find data. 342 Data analyses 343 Evolutionary link between masting intensity, temporal variability and autocorrelation 344 To explore how well the phylogenetic ancestry can explain masting behaviour, we first 345 estimated the phylogenetic signal (i.e. tendency for related species to resemble each 346 other more than they resemble species drawn at random from the phylogenetic tree) of 347 masting intensity (D) , PV and AR1 using the phylosig function in the R64 package 348 phytools65. Phylogenetic signal was assessed by the lambda (λ) metric, which varies 349 from 0 (where phylogenetic and trait similarity are totally independent) to 1 (where the 350 traits are completely explained by shared ancestry). We then used continuous mapped 351 phylogenetic trees (contMap function in R package phytools) to visualise their 352 phylogenetic signal. Finally, we explored the evolutionary relationship between 353 temporal variability, temporal autocorrelation and masting intensity using pairwise 354 correlations correcting for the phylogeny. Relationships were shown using 355 phylomorphospace plots 65, which depicts each species as a data point in a trait space, 356 together with the phylogenetic relationship of each species-point. 357 Controls of masting intensity and its mode of evolution 358 We first tested whether masting intensity (D) was related to climate, productivity, foliar 359 N and P concentrations and their interaction. To do so, we fitted phylogenetic linear 360 models in which the response variable was masting intensity (D) and the predictors 361 were foliar N and P concentrations and their interaction, MAT, MAP, MATs, MAPs, 362 MATD and MAPD for climate and GPP and GPPD for productivity. Phylogenetic models 363 were fitted optimising lambda (i.e., the strength of phylogenetic signal) and using the 364 phylolm function in the R phylolm package66. The final model was achieved by 365 removing the least significant terms from the full model, in a step-by-step process, until 366 all variable estimates were significant. The same models were fitted for the PV and 367 AR1, but changing the predictors from D to PV or AR1 respectively (e.g., MAT PV 368 instead of MATD when predicting PV). Because the CV has been widely used to assess 369 temporal variability, we also fitted a model using CV as the predictor variable to 370 compare its results with those of the model using the PV index. Masting intensity and 371 PV were transformed to natural logarithms to normalise model residuals. We used the 372 package visreg67 to visualise model results. 373 Evolution of masting intensity and foliar N and P and their interaction 374 To test the hypothesis that masting behaviour has evolved as an adaptation to nutrient 375 imbalances and low foliar N and P concentrations, we performed three step analysis. 376 First, we classified species as masting and non-masting. Second, we reconstructed the 377 ancestral state between the two types of reproducting behaviour, and third, we fitted 378 different evolutionary models to test whether foliar concentrations of N, P and N:P ratio 379 and N×P (N times P, as a measure of the overall availability of nutrients) evolved under 380 the reconstructed discrete selective regimes (masting or non-masting). 381 To define masting behaviour, and as a test for robustness of our results, we classified 382 species as subsets that represent masting (high temporal variability and strong 383 potential resource depletion) or non-masting (low temporal variability and no resource 384 depletion) behaviours based on 5 different thresholds of masting intensity (D), selecting 385 only the lower and upper bounds for the analyses and discarding the intermediate 386 species. The selected percentile thresholds were: i) from 0% to 33% non-masting 387 (N=38 species) and from 66% to 100% masting (N=58); ii) from 0% to 25% non- 388 masting (N=32) and from 75% to 100% masting (N=43); iii) from 0% to 20% non- 389 masting (N=27) and from 80% to 100% masting (N=34); iv) from 0% to 15% non- 390 masting (N=22) and from 85% to 100% masting (N=26) and v) from 0% to 10% non- 391 masting (N=16) and from 90% to 100% masting (N=17). 392 To reconstruct ancestral states of masting behaviour, for each of the five classifications 393 we performed stochastic character mappings68, which reconstructs the state of the 394 ancestors of a phylogeny based on the observed traits of the current species and the 395 phylogenetic structure. Ancestral reconstructions were performed using the 396 make.simmap function in the phytools R package65, simulating 1000 stochastic 397 ancestral reconstructions, specifying equal rates of transition amongst the character 398 states and using the ―mcmc‖ method. 399 Once we had the 1000 stochastic character mappings for each masting classification, 400 we performed generalised Ornstein–Uhlenbeck Hansen models to test whether the 401 inferred evolutionary trajectories in foliar concentrations of N, P, N:P and N×P 402 (hereafter ―continuous traits‖) were associated with the two alternative masting 403 behaviour strategies (hereafter ―ancestral states‖) and whether they followed an 404 adaptive (Ornstein–Uhlenbeck: OU) or random (Brownian motion—BM) model of 405 evolution69–71. To do so, we used the OU function from OUwie R package72. We fitted 406 five different models using the 1000 ancestral reconstructions mentioned above for 407 each classification. The five different models represent different types of underlying 408 evolutionary processes, being: i) single-state BM models (BM1), where evolutionary 409 rates for the continuous traits are equal for all ancestral states, ii) BM models with 410 different evolutionary rates for each ancestral state (BMS), iii) OU models with a single 411 optimal value for the continuous traits for all ancestral states (OU1), iv) OU models with 412 different optimal values but a single alpha (the strength of the pull towards the optimal 413 values of the trait) and single theta (the rate of phenotypic variation around the optimal 414 value) for each state (OUM), and v) OU models that assumed different optimal values 415 with multiple rates of phenotypic evolution (theta) for each state (OUMV). 416 Models containing negative eigenvalues (non-sound models) were deleted when 417 summarising our results29. Different evolutionary models were compared using second- 418 order Akaike information criterion (AICc) amongst all sound models. Those models with 419 the lowest AICc were considered to be those that fitted the data best. For models 420 assuming different optimal values of foliar N, P, N:P and N×P for masting and non- 421 masting species (i.e., OUM and OUMV models) we calculated the geometric mean of 422 the percentage differences of each model. Statistical differences in optimal values 423 estimated between subsets of masting and non-masting species by the OU models 424 with different state means (OUM and OUMV) were tested using paired t-tests. Given 425 that results pointed out to the same direction (see Supplementary Information) for all 426 masting intensity thresholds, we only show those from the intermediate (0% to 20% for 427 non-masting and 80% to 100% for masting). These analyses used the 168 species for 428 which we had data for masting intensity and foliar N and P concentrations. All statistical 429 analyses were perform with R statistical software version 3.5.164. 430 Data availability 431 The authors declare that the data supporting the findings of this study will be made 432 openly available upon acceptance at USGS ScienceBase. Any use of trade, firm, or 433 product names is for descriptive purposes only and does not imply endorsement by the 434 US Government. 435 436 References: 437 1. Salisbury, E. The reproductive capacity of plants. (Bell, 1942). 438 2. Silvertown, J. W. The evolutionary ecology of mast seeding in trees. Biol. J. Linn. 439 Soc. 14, 235–250 (1980). 440 3. Kelly, D. & Sork, V. L. Mast seeding in perennial plants: why, how, where? Annu. 441 Rev. Ecol. Syst. 33, 427–447 (2002). 442 4. Waller, D. M. Models of mast fruiting in trees. J. Theor. Biol. 80, 223–232 (1979). 443 5. Pearse, I. S., LaMontagne, J. M., Lordon, M. & Koenig, W. D. Biogeography and 444 phylogeny of masting: Do global patterns fit functional hypotheses? J. Ecol. 445 (Accepted article) 446 6. Norton, D. A. & Kelly, D. Mast Seeding Over 33 Years by Dacrydium 447 cupressinum Lamb. (rimu) (Podocarpaceae) in New Zealand: The Importance of 448 Economies of Scale. Funct. Ecol. 2, 399–408 (1988). 449 7. Pearse, I. S., Koenig, W. D. & Kelly, D. Mechanisms of mast seeding: resources, 450 weather, cues, and selection. New Phytol. 212, 546–562 (2016). 451 8. Espelta, J. M., Cortés, P., Molowny-horas, R., Sánchez-humanes, B. & Retana, 452 J. Masting Mediated by Summer Drought Reduces Acorn Predation in 453 Mediterranean Oak Forests. Ecology 89, 805–817 (2008). 454 9. Espelta, J. M. et al. Beyond predator satiation: Masting but also the effects of 455 rainfall stochasticity on weevils drive acorn predation. Ecosphere 8, e01836 456 (2017). 457 10. Donaldson, J. S. Mast-seeding in the cycad genus Encephalartos: a test of the 458 predator satiation hypothesis. Oecologia 94, 262–271 (1993). 459 11. Smith, C. C., Hamrick, J. L. & Kramer, C. L. The Advantage of Mast Years for 460 Wind Pollination. The American Naturalist 136, 154 (1990). 461 12. Fernández-Martínez, M., Belmonte, J., Maria Espelta, J. & Espelta, J. M. 462 Masting in oaks: Disentangling the effect of flowering phenology, airborne pollen 463 load and drought. Acta Oecologica 43, 51–59 (2012). 464 13. Herrera, C., Jordano, P., Guitián, J. & Traveset, A. Annual variability in seed 465 production by woody plants and the masting concept: reassessment of principles 466 and relationship to pollination and seed dispersal. Am. Nat. 152, 576–594 467 (1998). 468 14. Pearse, I. S., LaMontagne, J. M. & Koenig, W. D. Inter-annual variation in seed 469 production has increased over time (1900–2014). Proc. R. Soc. B Biol. Sci. 284, 470 20171666 (2017). 471 15. Fernández-Martínez, M., Vicca, S., Janssens, I. A., Espelta, J. M. & Peñuelas, J. 472 The role of nutrients, productivity and climate in determining tree fruit production 473 in European forests. New Phytol. 213, 669–679 (2016). 474 16. Isagi, Y., Sugimura, K., Sumida, a. & Ito, H. How Does Masting Happen and 475 Synchronize? J. Theor. Biol. 187, 231–239 (1997). 476 17. Ye, X. & Sakai, K. A new modified resource budget model for nonlinear 477 dynamics in citrus production. Chaos, Solitons & Fractals 87, 51–60 (2016). 478 18. Pérez-Ramos, I. M., Ourcival, J. M., Limousin, J. M. & Rambal, S. Mast seeding 479 under increasing drought: results from a long-term data set and from a rainfall 480 exclusion experiment. Ecology 91, 3057–68 (2010). 481 19. Fernández-Martínez, M., Vicca, S., Janssens, I. A., Espelta, J. M. & Peñuelas, J. 482 The North Atlantic Oscillation synchronises fruit production in western European 483 forests. Ecography (Cop.). 39, 864–874 (2016). 484 20. Sork, V. L., Bramble, J. & Sexton, O. Ecology of mast-fruiting in three species of 485 North American deciduous oaks. Ecology 74, 528–541 (1993). 486 21. Perez-Ramos, I. M., Padilla-Díaz, C. M., Koenig, W. D. & Marañón, T. 487 Environmental drivers of mast-seeding in Mediterranean oak species: Does leaf 488 habit matter? J. Ecol. 103, 691–700 (2015). 489 22. Fernández-Martínez, M., Bogdziewicz, M., Espelta, J. M. & Peñuelas, J. Nature 490 beyond Linearity: Meteorological Variability and Jensen’s Inequality Can Explain 491 Mast Seeding Behavior. Front. Ecol. Evol. 5, 1–8 (2017). 492 23. Wright, I. J. et al. The worldwide leaf economics spectrum. Nature 428, 821–7 493 (2004). 494 24. Domingues, T. F. et al. Co-limitation of photosynthetic capacity by nitrogen and 495 phosphorus in West Africa woodlands. Plant, Cell Environ. 33, 959–980 (2010). 496 25. Field, C., Merino, J. & Mooney, H. a. Compromises between water-use efficiency 497 and nitrogen-use efficiency in five species of California evergreens. Oecologia 498 60, 384–389 (1983). 499 26. Elser, J. J. et al. Nutritional constraints in terrestrial and freshwater food webs. 500 Nature 408, 578–580 (2000). 501 27. Sardans, J. et al. Factors influencing the foliar elemental composition and 502 stoichiometry in forest trees in Spain. Perspect. Plant Ecol. , Evol. Syst. 18, 52– 503 69 (2016). 504 28. Peñuelas, J. et al. The bioelements, the elementome and the ―biogeochemical 505 niche‖. Ecology 1–20 (2019). 506 29. Fernández-Martínez, M. et al. Nutrient-rich plants emit a less intense blend of 507 volatile isoprenoids. New Phytol. 220, 773–784 (2018). 508 30. Camarero, J. J., Albuixech, J., López-Lozano, R., Casterad, M. A. & Montserrat- 509 Martí, G. An increase in canopy cover leads to masting in Quercus ilex. Trees 510 24, 909–918 (2010). 511 31. Fernández-Martínez, M., Garbulsky, M., Peñuelas, J., Peguero, G. & Espelta, J. 512 M. Temporal trends in the enhanced vegetation index and spring weather predict 513 seed production in Mediterranean oaks. Plant Ecol. 216, 1061–1072 (2015). 514 32. Sala, A., Hopping, K., McIntire, E. J. B., Delzon, S. & Crone, E. E. Masting in 515 whitebark pine (Pinus albicaulis) depletes stored nutrients. New Phytol. 196, 516 189–99 (2012). 517 33. Reekie, E. G. & Bazzaz, F. A. Reproductive effort in plants. II Does carbon 518 reflect the allocation of other resources? Am. Nat. 129, 897–906 (1987). 519 34. Miyazaki, Y. et al. Nitrogen as a key regulator of flowering in Fagus crenata: 520 understanding the physiological mechanism of masting by gene expression 521 analysis. Ecol. Lett. 17, 1299–1309 (2014). 522 35. Han, Q., Kabeya, D., Iio, A., Inagaki, Y. & Kakubari, Y. Nitrogen storage 523 dynamics are affected by masting events in Fagus crenata. Oecologia 174, 679– 524 687 (2013). 525 36. Fujita, Y. et al. Low investment in sexual reproduction threatens plants adapted 526 to phosphorus limitation. Nature 505, 82–6 (2014). 527 37. Sardans, J., Rivas-Ubach, A. & Peñuelas, J. The C:N:P stoichiometry of 528 organisms and ecosystems in a changing world: A review and perspectives. 529 Perspect. Plant Ecol. Evol. Syst. 14, 33–47 (2012). 530 38. Güsewell, S. N:P ratios in terrestrial plants: Variation and functional significance. 531 New Phytologist 164, 243–266 (2004). 532 39. McGroddy, M. E., Daufresne, T. & Hedin, O. L. Scaling of C : N : P Stoichiometry 533 in Forests Worldwide : Implications of Terrestrial Redfield- Type Ratios. Ecology 534 85, 2390–2401 (2004). 535 40. Sardans, J. et al. Foliar and soil concentrations and stoichiometry of nitrogen 536 and phosphorous across European Pinus sylvestris forests: relationships with 537 climate, N deposition and tree growth. Funct. Ecol. 30, 676–689 (2016). 538 41. Eckstein, R. L. & Karlsson, P. S. Above-ground growth and nutrient use by 539 plants in a subarctic environment: effects of habitat, life-form and species. Oikos 540 79, 311–324 (1997). 541 42. Fernández-Martínez, M., Vicca, S., Janssens, I. A., Martín-Vide, J. & Peñuelas, 542 J. The consecutive disparity index, D, as measure of temporal variability in 543 ecological studies. Ecosphere 9, e02527 (2018). 544 43. Martín-Vide, J. Notes per a la definició d’un índex de «desordre» en 545 pluviometria. Soc. Catalana Geogr. 7, 89–96 (1986). 546 44. Heath, J. P. Quantifying temporal variability in population abundances. Oikos 547 115, 573–581 (2006). 548 45. Losos, J. B. Phylogenetic niche conservatism, phylogenetic signal and the 549 relationship between phylogenetic relatedness and ecological similarity among 550 species. Ecol. Lett. 11, 995–1003 (2008). 551 46. Elser, J. J. et al. Global analysis of nitrogen and phosphorus limitation of primary 552 producers in freshwater, marine and terrestrial ecosystems. Ecol. Lett. 10, 1135– 553 42 (2007). 554 47. LeBauer, D. S. & Treseder, K. K. Nitrogen limitation of net primary productivity in 555 terrestrial ecosystems is globally distributed. Ecology 89, 371–379 (2008). 556 48. Erel, R. Flowering and fruit set of olive trees. Am. Soc. 133, 639–647 (2008). 557 49. Sardans, J. & Peñuelas, J. Drought changes phosphorus and potassium 558 accumulation patterns in an evergreen Mediterranean forest. Funct. Ecol. 21, 559 191–201 (2007). 560 50. Elser, J. J. et al. Growth rate-stoichiometry couplings in diverse biota. Ecol. Lett. 561 6, 936–943 (2003). 562 51. Fernández-Martínez, M. et al. Nutrient availability as the key regulator of global 563 forest carbon balance. Nat. Clim. Chang. 4, 471–476 (2014). 564 52. Nogueira, C. et al. Leaf nutrients, not specific leaf area, are consistent indicators 565 of elevated nutrient inputs. Nat. Ecol. Evol. 3, 400–406 (2019). 566 53. Smaill, S. J., Clinton, P. W., Allen, R. B. & Davis, M. R. Climate cues and 567 resources interact to determine seed production by a masting species. J. Ecol. 568 99, 870–877 (2011). 569 54. Saura-Mas, S. & Lloret, F. Linking post-fire regenerative strategy and leaf 570 nutrient content in Mediterranean woody plants. Perspect. Plant Ecol. Evol. Syst. 571 11, 219–229 (2009). 572 55. Kelly, D., Hart, D. & Allen, R. B. Evaluating the wind pollination benefits of mast 573 seeding. Ecology 82, 117–126 (2001). 574 56. Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate 575 surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017). 576 57. Harris, I., Jones, P. D. D., Osborn, T. J. J. & Lister, D. H. H. Updated high- 577 resolution grids of monthly climatic observations - the CRU TS3.10 Dataset. Int. 578 J. Climatol. 34, 623 – 642 (2013). 579 58. Zhang, Y. et al. A global moderate resolution dataset of gross primary production 580 of vegetation for 2000–2016. Sci. Data 4, 170165 (2017). 581 59. Kattge, J. et al. TRY - a global database of plant traits. Glob. Chang. Biol. 17, 582 2905–2935 (2011). 583 60. Maitner, B. S. et al. The bien r package: A tool to access the 584 Botanical Information and Ecology Network (BIEN) database. Methods Ecol. 585 Evol. 2017, 1–7 (2017). 586 61. Rautio, P., Fürst, A., Stefan, K., Raitio, H. & Bartels, U. Sampling and Analysis of 587 Needles and Leaves. ICP Forests Manual XII, 19 (2010). 588 62. Cayuela, L. & Oksanen, J. Taxonstand: Taxonomic Standardization of Plant 589 Species Names. (2016). 590 63. Qian, H. & Jin, Y. An updated megaphylogeny of plants, a tool for generating 591 plant phylogenies and an analysis of phylogenetic community structure. J. Plant 592 Ecol. 9, 233–239 (2016). 593 64. R Core Team. R: A Lenguage and Environment for Stasitical Computing. (2018). 594 65. Revell, L. J. phytools: An R package for phylogenetic comparative biology (and 595 other things). Methods Ecol. Evol. 3, 217–223 (2012). 596 66. Tung Ho, L. S. & Ané, C. A linear-time algorithm for gaussian and non-gaussian 597 trait evolution models. Syst. Biol. 63, 397–408 (2014). 598 67. Breheny, P. & Burchett, W. Visualization of Regression Models Using visreg, R 599 package version 2.2-0. (2015). 600 68. Bollback, J. P. SIMMAP: Stochastic character mapping of discrete traits on 601 phylogenies. BMC Bioinformatics 7, 1–7 (2006). 602 69. O’Meara, B. C., Ané, C., Sanderson, M. J. & Wainwright, P. C. Testing for 603 different rates of continuous trait evolution using likelihood. Evolution (N. Y). 60, 604 922–933 (2006). 605 70. Thomas, G. H., Freckleton, R. P. & Székely, T. Comparative analyses of the 606 influence of developmental mode on phenotypic diversification rates in 607 shorebirds. Proc. R. Soc. London B 273, 1619–24 (2006). 608 71. Beaulieu, J. M., Jhwueng, D. C., Boettiger, C. & O’Meara, B. C. Modeling 609 stabilizing selection: Expanding the Ornstein-Uhlenbeck model of adaptive 610 evolution. Evolution (N. Y). 66, 2369–2383 (2012). 611 72. Beaulieu, J. M. & O’Meara, B. OUwie: Analysis of Evolutionary Rates in an OU 612 Framework. (2016). 613 614 Correspondence and requests for materials should be addressed to: 615 Marcos Fernández-Martínez: marcos.fernandez-martinez@uantwerpen.be 616 Acknowledgements 617 This research was supported by the Spanish Government project CGL2016-79835-P 618 (FERTWARM), the European Research Council Synergy grant ERC-2013-726 SyG- 619 610028 IMBALANCE-P, and the Catalan Government project SGR 2017-1005. 620 M.F.M.is a postdoctoral fellow of the Research Foundation – Flanders (FWO). MB was 621 supported by (Polish) NSF grants Sonatina 2017/24/C/NZ8/00151 and Uwertura 622 2018/28/U/NZ8/00003. The research was also supported by NSF grant DEB-1745496 623 to J.M.L. and DEB-1256394 to W.D.K. 624 625 Author Contributions 626 M.F-M., I.P., and I.A.J. conceived the paper. M.F-M and F.S analyzed the data. M.F- 627 M., J.S., J.P., I.P., W.K., and J.M.L, provided data. All authors contributed substantially 628 to the writing and discussion of the paper. 629 630 Figure captions 631 Figure 1: Masting behaviour intensity per species and its relationship with 632 potential resource depletion (negative AR1) and temporal variability (PV) of 633 reproductive effort. Panel (a) shows the reconstructed evolution of masting intensity 634 based on the disparity (D) index (see Methods) and 219 species. The inset graph 635 shows the distribution of values of masting intensity. Panels (b) and (c) are 636 phylomorphospace plots showing the evolutionary correlation between masting 637 intensity (ln(D), natural logarithm of disparity), potential resource depletion (AR1, 638 negative values indicate potential resource depletion may happen) and temporal 639 variability (PV). Plant silhouettes in panel (a) were drawn by FS and are available at 640 PhyloPic (http://phylopic.org). 641 Figure 2: 3D graph showing the interaction between foliar N and P on masting 642 intensity. Combined high foliar N and P concentrations decreased masting intensity, 643 although the highest masting intensity was found in species with the highest N to P 644 imbalances. Response surface of the negative interaction between N and P (P = 0.007) 645 was estimated from a phylogenetic linear model (Table 1). See Methods for further 646 information about the models. 647 Figure 3: Different optimum values of foliar N and P for subsets of masting and 648 non-masting species. Panel a shows model results for foliar N and P concentrations 649 for a subset of species identified as masting (high masting intensity) and non-masting 650 based on percentiles 20th and 80th of the D distribution (see Methods). Panel b shows 651 model results for foliar N:P and N×P (overall nutrient availability) for masting and non- 652 masting species. For panels a and b, 1000 and 927 stochastic character mapping 653 simulations were respectively used (see Table S2). For foliar N:P results from OUMV 654 models were used, while for N, P and N×P, OUM models were used (see Table S1 655 and Table S2). See Methods for further information. 656 Table 1: Model summary of the phylogenetic linear models based on 168 657 species. Masting intensity was defined as the D index, PV indicates temporal 658 variability and AR1 potential resource depletion coefficient (i.e., negative values 659 indicate potential resource depletion). Beta (β) shows the standardised model 660 coefficients, followed by the standard error of the mean (s.e.m), the P-value (P) and 661 lambda (λ) indicating the phylogenetic signal in the models. See Methods for further 662 information. Acronyms: mean annual precipitation (MAP), foliar nitrogen (N), foliar 663 phosphorus (P), mean annual temperature (MAT), temporal autocorrelation at lag 1 664 (AR1). 665 666 Figure 1 667 668 Figure 2 669 670 Figure 3 671 672 673 Table 1 β s.e.m P λ Masting intensity (D) MAP 0.22 0.09 0.016 N 0.04 0.10 0.708 P -0.04 0.10 0.674 N×P -0.19 0.07 0.007 0.47 Temporal variability (PV) MAPPV 0.16 0.08 0.042 MAT -0.18 0.09 0.040 P -0.16 0.08 0.044 0.41 Potential resource depletion (AR1) MAPAR1 0.27 0.08 0.001 N -0.09 0.10 0.001 P -0.05 0.10 0.008 N×P 0.22 0.07 0.002 0.00 674 675 Supplementary Material 676 Figure S1: Evolutionary relationship between potential resource depletion coefficient 677 (AR1) and temporal variability (PV) in seed production shown in a continuous trait 678 phylogenetic reconstruction (a) and a phylomorphospace plot (b). Potential resource 679 depletion and variability in seed production were not evolutionary correlated. Negative 680 values of AR1 indicate that potential resource depletion may happen, see Methods. 681 N=219 species. 682 683 Figure S2: Phylogenetic tree including the subset of low (non-masting) and high 684 masting intensity (masting) species used to perform the generalised Ornstein- 685 Uhlenbeck model results presented in the main text (20th – 80th percentile thresholds 686 for non-masting and masting species, Figure 3, Table S1 and S2). The phylogenetic 687 tree includes the estimated probability that ancestor nodes were masting or non- 688 masting species (large circles) as pie charts. Small circles indicate the current category 689 of the species. The ancestral character reconstruction was performed using 1000 690 stochastic character-mapped trees (see Methods for further information). 691 692 693 Figure S3: Evolutionary relationship between foliar N and P shown in a continuous trait 694 phylogenetic reconstruction (a) and a phylomorphospace plot (b). Foliar N and P 695 concentrations were evolutionary correlated. N=168 species. 696 697 Table S1: Mean differences (ΔAICc, second-order Akaike information criterion) 698 between each of the model’s AICc and the model with the lowest AICc. Evolutionary 699 models were Brownian motion (BM1, BMS) and generalised Ornstein-Uhlenbeck- 700 based Hansen (OU1, OUM, OUMV), fitting ―masting‖ and ―non-masting‖ species-state 701 and foliar nutrient concentrations (N: nitrogen, P: phosphorus, N:P: ratio N-to-P and, 702 N×P: N times P (overall nutrient availability). Average AICc values were calculated 703 using the subset of models in which none of them presented negative eigenvalues 704 (sound models, column N). Non-masting and masting columns indicate the number of 705 species used in each category depending on the percentile of masting intensity used to 706 classify species as non-masting (i.e., higher than e.g., 33%) and masting (i.e., lower 707 than e.g., 66%). Models with ΔAICc lower than 2 (indicating equal performance) were 708 highlighted. See Methods for further information. BM1 BMS OU1 OUM OUMV N Non-masting Masting 33 - 66% 38 58 N 36.36 30.79 0.00 1.35 3.19 1000 P 42.52 46.34 0.51 1.39 0.00 1000 N:P 123.81 98.59 18.25 20.23 0.00 903 N×P 31.99 35.48 0.00 0.69 0.25 1000 25 - 75% 32 43 N 35.51 26.82 0.00 1.11 2.43 1000 P 33.48 36.88 0.00 0.77 1.04 1000 N:P 112.98 85.34 21.92 23.69 0.00 906 N×P 27.60 30.71 0.00 0.50 1.34 1000 20 -80% 27 34 N 29.05 17.30 0.00 1.11 1.28 1000 P 29.46 32.96 0.00 0.54 1.35 1000 N:P 97.95 72.61 19.86 21.29 0.00 927 N×P 23.34 26.12 0.00 0.46 2.07 1000 15 - 85% 22 26 N 25.21 18.42 0.00 0.68 2.34 1000 P 23.22 27.27 1.08 0.27 0.00 998 N:P 97.78 68.37 30.33 31.11 0.00 890 N×P 15.37 18.94 0.25 0.00 0.96 1000 10 - 90% 16 17 N 24.60 19.25 0.00 2.15 4.33 999 P 14.88 17.99 0.00 1.78 3.96 999 N:P 76.40 52.25 29.94 31.03 0.00 914 N×P 11.35 13.51 0.00 1.89 4.31 1000 709 Table S2: Estimated foliar nitrogen (N) and phosphorus (P) concentrations, N:P and N×P (overall nutrient availability) optimal values for 710 masting and non-masting species using OUMV and OUM models (see Methods for further information about the models), chosen based on the 711 lowest ΔAICc estimating different state means for masting and non-masting species (Table S1). Masting and non-masting species were 712 classified depending on the percentile of masting intensity (e.g., masting for higher than 66%, non-masting for lower than 33%, see subheaders 713 within the table). Columns 2.5%, 50 and 97.5% indicate, for masting and non-masting species, the percentiles of the optimal values based on 714 the sound models (without negative eigenvalues, N column) used. M>N% indicate the percentage of models in which masting species 715 presented average higher N, P, N:P or N×P optimal values than non-masting species. ΔM-N, followed by s.e.m (standard error of the mean), 716 indicate the paired (across simulations) difference between optimal values in masting and non-masting species. P (t-test) shows the P-value of 717 the paired t-test testing for differences in the mean optimal values of masting and non-masting species. ΔM-N%, followed by s.e.m., indicates 718 the average percentual difference (geometric, paired differences) in mean optimal values between masting and non-masting species. 719 720 Table S2: Masting Non-masting 2.5% 50% 97.5% 2.5% 50% 97.5% M>N% ΔM-N s.e.m P (t-test) ΔM-N% s.e.m N Model 33 - 66% N 1.59 1.78 1.92 1.70 1.92 2.13 20.5% -0.15 0.01 <0.001 -7.6% 0.4% 1000 OUM P 0.15 0.17 0.18 0.17 0.20 0.21 6.2% -0.03 0.00 <0.001 -14.7% 0.3% 1000 OUMV N:P 10.57 11.15 11.44 10.46 10.79 11.28 83.5% 0.31 0.01 <0.001 2.9% 0.1% 903 OUMV N×P 0.26 0.31 0.38 0.33 0.42 0.50 9.3% -0.10 0.00 <0.001 -24.3% 0.7% 1000 OUMV 25 - 75% N 1.57 1.73 1.91 1.74 1.94 2.11 11.5% -0.20 0.01 <0.001 -10.5% 0.3% 1000 OUM P 0.14 0.16 0.18 0.17 0.19 0.22 10.1% -0.03 0.00 <0.001 -15.4% 0.4% 1000 OUM N:P 10.77 11.38 11.93 10.43 10.79 11.17 92.4% 0.58 0.01 <0.001 5.3% 0.1% 906 OUMV N×P 0.23 0.31 0.39 0.33 0.42 0.51 9.9% -0.11 0.00 <0.001 -27.2% 0.8% 1000 OUM 20 - 80% N 1.55 1.76 2.02 1.68 1.97 2.18 23.3% -0.19 0.01 <0.001 -9.5% 0.4% 1000 OUM P 0.14 0.16 0.19 0.18 0.20 0.23 9.0% -0.04 0.00 <0.001 -18.3% 0.5% 1000 OUM N:P 10.85 11.74 12.53 10.06 10.51 10.88 97.7% 1.20 0.02 <0.001 11.4% 0.2% 927 OUMV N×P 0.22 0.32 0.42 0.33 0.45 0.54 11.9% -0.13 0.00 <0.001 -28.7% 1.0% 1000 OUM 15 - 85% N 1.54 1.71 1.95 1.77 2.03 2.22 10.00% -0.29 0.01 <0.001 -14.3% 0.4% 1000 OUM P 0.14 0.16 0.18 0.20 0.22 0.24 0.9% -0.06 0.00 <0.001 -26.1% 0.3% 998 OUMV N:P 10.44 11.88 12.80 9.28 9.71 9.99 99.7% 2.12 0.02 <0.001 21.7% 0.2% 890 OUMV N×P 0.21 0.30 0.41 0.38 0.49 0.59 3.80% -0.19 0.00 <0.001 -38.7% 1.0% 1000 OUM 10 - 80% N 1.67 1.87 2.05 1.64 1.85 2.05 54.2% 0.02 0.01 <0.001 1.3% 0.4% 999 OUM P 0.14 0.18 0.21 0.18 0.20 0.24 12.0% -0.03 0.00 <0.001 -14.8% 0.5% 999 OUM N:P 9.90 12.52 13.96 9.49 9.85 10.22 96.6% 2.46 0.04 <0.001 24.5% 0.3% 914 OUMV N×P 0.23 0.36 0.48 0.33 0.45 0.56 20.6% -0.09 0.00 <0.001 -22.2% 2.0% 1000 OUM 721