The resource availability hypothesis predicts that plants adapted to infertile soils have high levels of anti-herbivore leaf defenses. This hypothesis has been mostly explored for secondary metabolites such as phenolics, while it remains underexplored for silica-based defenses. We determined leaf concentrations of total phenols and silicon (Si) in plants growing along the 2-million-year Jurien Bay chronosequence, exhibiting an extreme gradient of soil fertility. We found that nitrogen (N) limitation on young soils led to a greater expression of phenol-based defenses, whereas old, phosphorus (P)-impoverished soils favored silica-based defenses. Both defense types were negatively correlated at the community and individual species level. Our results suggest a tradeoff among these two leaf defense strategies based on the strength and type of nutrient limitation, thereby opening up new perspectives for the resource availability hypothesis and plant defense research. This study also highlights the importance of silica-based defenses under low P supply.

Study area and site description

The 2 Ma Jurien Bay chronosequence is located in south-western Australia, approximately 200 km north of Perth (Fig. S1), and is described in detail in Laliberté et al. (2012) and Turner & Laliberté (2015). The chronosequence, part of the Swan Coastal Plain, comprises a series of dunes parallel to the coast, formed by periodic interglacial sea-level high-stands since the Early Pleistocene/Late Pliocene (Kendrick et al. 1991), with a clear gradient of soil age with increasing distance from the Indian Ocean. The dunes comprise three units: the Quindalup dunes date from the Holocene (up to 7 ka), the Spearwood dunes from the Middle Pleistocene (120 to 500 ka) and the Bassendean dunes from the Early Pleistocene or Late Pliocene (~2 Ma) (McArthur & Bettenay 1974; Playford et al. 1976). The parent material of the dunes is calcareous sand from the nearshore coastal environment (Turner & Laliberté 2015). The climate is Mediterranean, with a mean annual temperature of 19°C, mean annual rainfall 533 mm and potential annual evapotranspiration of 1433 mm, which results in a water balance of –900 mm yr-1 (data from the Jurien Bay Bureau of Meteorology from 1968 to 2015 in Turner et al. (2018).

We selected the same five chronosequence stages as in Hayes et al. (2014); these include both the early and retrogressive phases of long-term ecosystem development. The main soil properties of these five stages can be found in Table 1. Soil total P and carbonate concentrations, cation exchange capacity and pH-CaCl2 continually decrease with increasing soil age. Soil total N concentrations increase from stage 1 to stage 2 during the progressive phase of ecosystem development, then decrease towards the last stages during the retrogressive phase (Laliberté et al. 2012; Turner & Laliberté 2015). Plant growth is most strongly limited by low N availability in the early stages, and by P availability in the advanced stages (Laliberté et al. 2012; Hayes et al. 2014). Previous studies showed that plant-available [Si] is low in the early stages of soil development, increases in stage 4 in the Spearwood dune system, and finally decreases in the oldest stage of soil development, where it is controlled by intense biocycling (Table 1) (de Tombeur et al. 2020b, c).

Site selection

For each chronosequence stage, we randomly selected five plots (10 m × 10 m each) among the 10 plots already characterized for soil and vegetation in previous studies (Hayes et al. 2014; Zemunik et al. 2016) (Fig. S1). The plots were originally selected using a random stratified sampling design (Zemunik et al. 2016). To characterize vegetation, seven 2 m × 2 m subplots were randomly positioned in each plot in which all individuals of all vascular plant species were counted  (Zemunik et al. 2016). The percent canopy cover of each species was estimated, and the relative cover of each species was calculated as a fraction of the total canopy cover over the seven subplots (Zemunik et al. 2016).

Sampling procedure

In the 25 plots selected, we sampled leaves according to two procedures. First, we sampled leaves from one individual plant for each of the 10 most-abundant species of each plot as defined in Zemunik et al. (2016). The number of leaves sampled per individual was adapted according to their mass, but was never less than 10. Occasionally, a species originally included in the 10 most-abundant species was not found on the plot, which resulted in less than 10 species for some plots (Table S1). The 234 species sampled with this first procedure still accounted for 57% to 88% of the total cover of each plot (Table S1). The community-level analyses were performed only on these species. Second, we systematically sampled the species belonging to nine families, even if they were not included in the 10 most-abundant species, in order to study family-level variation in leaf [Si] and [phenols], following the same sampling procedure: Asparagaceae, Cyperaceae, Ericaceae, Fabaceae, Haemodoraceae, Myrtaceae, Poaceae, Restionaceae and Rhamnaceae (Table S2). These families were selected because they were well represented and found at all stages of the chronosequence (Zemunik et al. 2016), and likely had contrasting [Si] based on known phylogenetic patterns (Hodson et al. 2005). In total, 298 leaf samples belonging to 24 families were collected (Tables S1 and S2).

All leaf material was collected over two weeks in November 2018. Leaves were sampled from one healthy mature individual plant per species in each plot; when an individual did not provide sufficient biomass for analysis (e.g., Poaceae spp.), leaf samples from several individuals within the plot were combined.

Leaf analyses

Leaves were washed with distilled water, dried at 70 °C for 48 h and finely ground. Leaf material (0.5 g) was placed in a porcelain crucible and calcinated at 450°C for 24 h. The weight after calcination was used to calculate the ash content. The ash was mixed with 1.6 g lithium-metaborate and 0.4 g of lithium-tetraborate in a graphite crucible and heated at 1000 °C for 5 min (Chao & Sanzolone 1992). The bead was then dissolved in 15% HNO3 and the concentrations of Si, P, calcium (Ca), magnesium (Mg) and potassium (K) were determined by inductively coupled plasma-optical emission spectrometry (Agilent Technologies, 700 series ICP-OES). Phenolic compounds were extracted from a 0.25 g ground sample stirred with 10 mL of 70% acetone for 30 min (Salminen & Karonen 2011; Schaller et al. 2012; Bettaieb Rebey et al. 2020). Total phenols were determined in triplicate as described in Salminen & Karonen (2011) using a Folin-Ciocalteu assay with gallic acid monohydrate as standard (Merckx, Darmstadt, Germany). Total phenol concentrations were expressed as g of gallic acid equivalents (GAE) per kilogram of dry weight.

Soil sampling and analyses

In order to determine how Si availability in soils affected species-level variations in leaf [Si], three soil samples (top 20 cm) were taken in each of the 25 plots, for a total of 75 soil samples. Samples were air-dried and sieved (< 2 mm). The pool of ‘plant-available Si’ was determined by extraction in 0.01 M CaCl2 (Haymsom & Chapman 1975; Sauer et al. 2006). Soil was shaken for 5 h in a 1:10 soil-to-solution ratio, filtered (cellulose filter, pore size < 2 µm, Healthcare Whatman™), acidified with 50 µL of ultrapure 65% HNO3, and stored in darkness at 4°C prior to Si determination by ICP-OES.

Data analyses

To characterize leaf [phenols], [Si], [Ca], [Mg], [K] and [P] in plant communities across the chronosequence, we calculated the mean values of the 10 most-abundant species per plot, weighted or not by their relative canopy cover. The cover-weighted mean (CWM) was calculated as follows (Garnier et al. 2004; Violle et al. 2007):

where ti and RCi are, respectively, the value of the trait t and its relative cover RC for a species i and S is the number of species.

The differences in plant-available [Si], leaf [phenols], [Si], [P], [Ca], [Mg] and [K] across the chronosequence stages were tested by one-way analysis of variance (ANOVA), followed by post-hoc multiple comparison (Fisher’s Least Significant Difference [LSD] tests). When these analyses considered all individuals together (i.e., not the mean and CWM of the plant communities), we treated species and plots as a random factor (mixed-effect models). We tested the relationships between leaf [Si], [phenols] and major soil properties (total P, total N, ratio soil N:P) with linear mixed-effect models, treating plot and species as random factors when all individuals were considered together, and treating chronosequence stage as a random factor when the means and CWM of the 25 plots were considered. We also explored relationships between leaf [Si], [phenols] and foliar nutrient concentrations through Pearson tests of correlation. For the nine plant families selected, we tested the differences in leaf [phenols] and [Si] across the chronosequence stages using mixed-effect models with species and plot as random factors, followed by Fisher’s LSD tests, and we tested the relationships between leaf [Si] and [phenols] with linear mixed-effect models (with plot and species as random factors). Additionally, we studied intraspecific variation for seven taxa, and detailed explanation of the statistical analyses used is presented in the Supplementary Information. Finally, a t-test was performed to examine the differences in leaf [Si] between dicots and monocots. All residuals were visually inspected for heteroscedasticity and appropriate transformations were performed to meet the model assumptions. All analyses were conducted in R using the ‘nlme’ (Pinheiro et al. 2020) and ‘multcomp’ (Hothorn et al. 2008) packages.