Microbial physiology conserves phosphorus across long-term ecosystem development

One of the paradigms of terrestrial biogeochemistry is that the long-term trajectory of ecosystem development—from primary succession through to retrogression—is paralleled by a decline in the availability of soil phosphorus (P)^1,2,3. This decline is due to a combination of progressive soil weathering and increasing sequestration of P in organic compounds that become tightly sorbed to soil minerals^1,4,5,6 and results in an increased tendency towards P-limited productivity^3,7,8, a proliferation of specialized plant taxa with highly specific P requirements^9,10 and an increasing degree of microbial control over belowground P cycling¹¹. In late successional and retrogressive ecosystems, a large proportion of soil P is bound within the soil microbial biomass^11,12,13. Thus, physiological traits that affect the efficiency of intracellular allocation and retention of P by the soil microbial biomass could exert strong control over rates of P supply for plant uptake and, by extension, plant community productivity and composition, particularly in later stages of ecosystem development. However, our understanding of the physiological strategies used by soil microbes for P efficiency is limited¹⁴, particularly relative to those of plants, the study of which is a major goal of contemporary plant ecophysiology¹⁵.

One strategy for high P efficiency observed widely across plants is that of phospholipid substitution, in which PO₄–bearing membrane lipids (phospholipids) are replaced by non-phosphorous membrane lipids such as galactolipids and sulfolipids^16,17. Comparable lipid remodelling has been observed for marine microorganisms¹⁸, and recent studies have shown that soil fertility can influence soil microbial lipid composition^19,20, with phospholipid substitution occurring in direct response to microbial P limitation²¹. Another strategy for nutrient efficiency found in plants is a high allocation of ‘excess’ photosynthate (carbon (C)) to C-rich organic structures^22,23. Such structures may serve structural and protective purposes, as in the cases of woody tissues (thick, protective bark and seed capsules) and defensive compounds (for example, tannins and polyphenols), or storage purposes as in the case of starch-based lignotubers²². Analogous storage compounds have been observed in ex situ microbial populations and include triacylglycerols (TAGs), which are found across both eukaryotes and prokaryotes, polyhydroxybutyrate (PHB, exclusive to prokaryotes), glycogen, wax esters and cyanophycin²⁴. Recent studies of natural soil communities indicate that C storage represents a major pathway of soil microbial growth²⁵ that varies with soil fertility and microbial P demand²⁶. The influence of C storage on microbial nutrient efficiency is revealed by modelling^24,27, which indicates that C storage enhances microbial nutrient acquisition and retention by reducing stoichiometric imbalances between microorganisms and their substrates and by providing an energy source that enables rapid uptake of scarce resources^24,25.

To investigate the physiological strategies deployed by soil microorganisms to cope with declining soil P throughout long-term ecosystem development, we studied an ~700 kyr successional chronosequence at Cooloola in eastern Australia (26° S, 153° E)²⁸. The Cooloola sandmass comprises six to ten Holocene and late Pleistocene dune systems with estimated ages ranging from 0–0.04 thousand years (ka) to 350–730 ka. The dune systems were formed by marine and/or aeolian deposition of siliceous sands from the inner continental shelf during the last ~700 kyr^28,29. Ecosystems associated with each dune system are based on equivalent quartz-based substrates, such that Cooloola represents a true ecological chronosequence and a pristine example of ecosystem development in a largely unmodified landscape³⁰. We collected soils from numerous sites on dune crests across the chronosequence (Fig. 1a and Extended Data Table 1), quantified microbial phospholipid substitution, intracellular C storage (TAG and PHB), resource limitation and C/nitrogen (N)/P stoichiometry, and applied a stoichiometric model to evaluate how microbial physiology might influence net P mineralization and consequently primary productivity throughout ecosystem development.

a, Map of the Cooloola dune field showing 24 sampling sites representing the 6 major dune systems at Cooloola (DS1–DS6), which comprise an ~700 kyr chronosequence of soil formation (P for 0–30 cm soil depth). b, Soil total P as a function of dune age. c, The contribution of MBP to soil total P as a function of dune age (white circles, 0–5 cm depth; grey circles, 5–30 cm depth); P values indicate significance of dune age or the dune age × depth interaction based on F values from Wald tests (two-sided; n = 48); grey shaded bands indicate 95% confidence bands. Credit: a, © Stadia maps © OpenMapTiles © OpenStreetMap © Stamen design.

Full size image

Dynamics of soil P and ecosystem productivity at Cooloola are mostly typical of long-term pedogenesis and primary succession worldwide^1,3,11,31,32. Soil total P levels in the 0–30 cm depth decline from 35–39 mg kg⁻¹ in the youngest dunes, which support early successional plant species (for example, Allocasuarina littoralis, Alphitonia excelsa, Banksia intermedia) to 9–10 mg kg⁻¹ in the oldest dunes, which support retrogressive ‘wallum’ heath vegetation (Table 1, Fig. 1b and Extended Data Fig. 1). Phosphorus extracted by the Bray-2 extractant (which corresponds to PO₄³⁻–P in solution plus weakly sorbed PO₄³⁻–P, and is therefore a rough index of plant-available P) follows a similar trend, ranging from 6.5 to <0.1 mg kg⁻¹ across the sequence (Table 1). Microbial biomass P (MBP) increases throughout succession, from 0.2 mg kg⁻¹ in young dunes, peaking at 4.7 mg kg⁻¹ in middle-aged dunes, and subsequently declines to 2.4 mg kg⁻¹ in the oldest dunes. Across all but the youngest dunes, microbial P exceeds Bray-2 P 6–97-fold and comprises 6.9–52.0% (dune system averages) of total soil P (Fig. 1c), suggesting strong potential for plant–microbial competition for P and extensive microbial control over P cycling and availability¹¹. Thus, any ecophysiological strategies that shape the soil microbial P efficiency could influence ecosystem productivity and composition throughout much of primary succession and retrogression.

We found substantial microbiome-level phospholipid substitution and intracellular C storage allocation at Cooloola (Figs. 2 and 3), with the extents of both strategies reaching their maxima in retrogressive ecosystems and increasing with declining soil total P. Phospholipids were replaced extensively with non-phosphorous betaine lipids in older dunes, with the phospholipid contribution to total membrane lipid declining from 80% in dune system one to 60% in dune system six, and the estimated contribution of phospholipids to total microbial P declining from 36.0% to 6.8%. Microbial allocation to TAG (soil TAG content scaled to soil microbial polar lipid fatty acid (PLFA) content), the primary intracellular C storage compound for eukaryotic microorganisms (for example, fungi), tended to be greatest in old dune systems, particularly in subsurface soils, with soils of dune system five having particularly high TAG allocation. Likewise, bacterial PHB allocation increased sharply throughout ecosystem development.

a,b, Contributions of phospholipids (triangles) and non-phosphorous betaine lipids (circles) to total soil membrane lipid content in 0–5 cm soils as functions of dune age (a) and soil total P content (b) (n = 24). c,d, Contribution of phospholipid-P to total soil MBP as functions of dune age (c) and soil total P content (d) (n = 48). P values indicate significance of predictor variables based on F values from Wald tests (two-sided); grey shaded bands indicate 95% confidence bands.

Full size image

a,b, Allocation of C to TAG as a function of dune age (a) and soil total P content (b), scaled to soil microbial PLFA content (n = 48). c,d, Extent of allocation of C to PHB as a function of dune age (c) and soil total P content (d), scaled to soil bacterial PLFA content (n = 36). P values indicate significance of predictor variables based on F values from Wald tests (two-sided); grey shaded bands indicate 95% confidence bands. Allocation corresponds to soil TAG/PHB content for a given value of PLFA; that is, plots reflect trends after the significant, positive effect of notionally ‘structural’ microbial biomass (as indexed by PLFA content) has been taken into account.

Full size image

These findings indicate that at least two of the physiological strategies for nutrient efficiency exhibited by soil microorganisms—phospholipid substitution and C storage—parallel those of plants^16,17,23 and follow similar dynamics to those of plants across long-term ecosystem development. At Cooloola and other long-term chronosequences worldwide, plant P content declines as ecosystems mature (Extended Data Fig. 2)^11,32,33, and plant taxa associated with retrogressive ecosystems in Western Australia exhibit superior phospholipid substitution and P efficiency relative to other taxa^34,35. Likewise, the belowground C storage organs (lignotubers) of resprouting Eucalpytus species tend to be enlarged in response to declining soil nutrient availability at Cooloola and elsewhere^36,37. These organs enable rapid resprouting and uptake of nutrient pulses in post-fire environments^22,38,39, such that their fitness benefits to plants are comparable to the theoretical advantages conferred by microbial C storage compounds under fluctuating resource regimes²⁷. The coordination of microbial physiological traits that simultaneously lower cellular P and increase cellular C—and thereby increase biomass C/P ratios—in low-P environments is analogous to the ‘sclerophylly’ characteristic of vegetation adapted to infertile soils²². Mean microbial biomass C/N/P values at Cooloola (89/14/1 atomic; 0–30 cm depth) are roughly double the global mean (estimates range from 42/6/1 to 60/7/1)^40,41. Our work therefore provides a valuable basis for integrating soil microorganisms—the biogeochemical ‘engines’ of ecosystems—into existing theoretical frameworks for the ecology of strongly weathered, nutrient-depleted landscapes^22,42. In this context, it will be necessary to understand not just how microbial and plant P management strategies are similar, but also how they are likely to diverge in many ways given that plants and microbes tend strongly towards nutrient and C limitation, respectively.

The strong physiological responses of soil microorganisms to soil aging and ecosystem development were not paralleled by equally strong variation in microbial biomass, metabolic activity or resource (C, N or P) limitation (Fig. 4 and Extended Data Fig. 3), consistent with the theoretical expectation that microbial strategies for nutrient efficiency buffer against stoichiometric imbalance between microorganisms and their substrates²⁷. Microbial biomass C (Fig. 4a) and respiration (Fig. 4b) increased sharply from dune system one to two, but thereafter declined only gradually, at much slower rates than soil P. Positive linear relationships between soil microbial biomass C and both soil C/N and soil C/P ratios (Fig. 4c,d) along with strong respiratory responses to C addition (Fig. 4e) suggest that microbial C limitation was present and largely constant across the chronosequence. This was further supported by enzyme vector-length analysis (based on potential activities of β-glucosidase, N-acetyl-β-D-glucosaminidase and acid phosphomonoesterase + phosphodiesterase at 37 °C; Extended Data Fig. 3a). Likewise, moderate respiratory responses to P but not N (Fig. 4e), combined with enzyme vector-angle analysis (Extended Data Fig. 3b; but note the limitations of this approach as outlined in Methods), indicated that microorganisms tended towards relative P limitation versus N limitation across most of the chronosequence, with only modest (1.46-fold) variation in the strength of relative P limitation despite the ~fivefold gradient of surface soil total P and the ~138-fold gradient of Bray-2-extractable P. Thus, we suggest that microbial phospholipid substitution and intracellular C storage could contribute to the long-term stability of biogeochemical functioning in retrogressive ecosystems at Cooloola and elsewhere²⁶. For example, we estimate that substitution of phospholipids with betaine lipids supports 16–41% of the surface soil microbial biomass (less storage C) in the terminal dune system at Cooloola (Fig. 5e and equation (1)), with corresponding implications for the cycling and retention of nutrients and C in soil¹¹.

a,b, Trends in the dynamics of surface soil MBC (n = 24) (a) and surface soil CO₂ respiration (n = 118) (b). Inset (b): respiration per unit MBC, qCO₂ (n = 24) throughout pedogenesis. c,d, Soil microbial biomass expressed as functions of soil total C/N (c) and C/P (d) ratios (n = 48). e, Soil respiratory responses (log_e[treatment/control]) to factorial additions of C, N and P (n = 4 for all dune system × treatment combination; red arrows provide pairwise mean comparisons based on Tukey adjustments); P values indicate significance of predictor variables in linear models based on F values from Wald tests (two-sided); grey shaded bands indicate 95% confidence bands.

Full size image

a–d, Correlations between MBC/MBP ratios as per direct measurement and calculated microbial C/P ratios if the nutrient conservation strategies are discounted (a), MBN/MBP ratios as per direct measurement and calculated microbial N/P if phospholipid substitution is discounted (b), soil total C/P ratios and MBC/MBP ratios calculated with and without nutrient conservation strategies (c) and soil total N/total P ratios and MBN/MBP ratios calculated with and without nutrient conservation strategies (d); n = 21 (only surface soil data are shown, and three of the replicates for dune system one had no detectable MBP). e,f, The modelled implications of microbial nutrient efficiency strategies for the notional proportion of surface soil structural microbial biomass that is supported by phospholipid substitution (e) and the potential for microbially mediated mineralization of P with and without phospholipid substitution and intracellular C storage (f), as revealed by comparison of measured soil C/P ratios (C/P_soil) with modelled threshold elemental ratios for C/P_soil (TER_C/P) with and without nutrient efficiency strategies: values of natural log-transformed ratios of C/P_soil/TER_C/P > 0 indicate net microbial P mineralization, while values <0 indicate net P immobilization; model assumes fixed values of the uptake coefficient α (9.05) and microbial C use efficiency (0.3).

Full size image

At face value, the concurrence of persistent microbial C limitation and extensive intracellular C storage appears paradoxical; however, there are at least two reasons why these are not necessarily contradictory outcomes. First, this stored C could represent the ‘reserve’ storage mode, which is constitutive (determined by genes rather than environmental conditions) and has been documented within soil microbiomes that are demonstrably C limited^25,26. Microbial taxa that allocate high proportions of intracellular C to reserve storage tend to dominate in low-fertility soils, leading to high levels of C storage on a microbiome level²⁶. Second, previous evidence shows that P-limited taxa can form a subset of an apparently overarchingly C-limited soil microbiome²⁶. A subset of the bacterial assemblage, in particular, could be P limited and consequently allocating an increasing amount of surplus C to PHB as pedogenesis advances. Additional work is warranted to tease out the role of taxonomic turnover and reserve versus surplus C storage in shaping the ecophysiological responses of soil microbial communities to long-term pedogenesis.

Although the tendency towards ecosystem retrogression in the absence of major disturbance has been observed across several biomes, the trajectories of microbial biomass and nutrient limitation throughout long-term ecosystem development vary on the basis of factors such as climate, parent material and chronosequence duration³. We suggest that our findings are particularly applicable to ecosystems based on low-fertility parent materials and those at low latitudes where mineral weathering and podzolization proceed rapidly (for example, Extended Data Fig. 4), such that evolution has optimized microbial physiology to low P availability. However, evidence from other chronosequences with comparatively P-rich parent material indicates that buffering of soil microorganisms against the soil P depletion that drives ecosystem retrogression could be widespread. For example, soil microbial biomass remains essentially stable across the latter stages of the retrogressive Jurien Bay chronosequence (120,000–2,000,000 years) in Western Australia, across which soil P levels decline threefold and microorganisms are P limited¹³. Soil microbial biomass is similarly uniform across the retrogressive Waitutu⁴³ and Franz Josef⁴⁴ chronosequences in New Zealand, where microbial P demand increases sharply in older chronosequence stages^11,43,44. We suggest that this buffering has a physiological basis and anticipate that physiology will emerge as a major driver of the diverse structural and functional trajectories of plant and microbial communities throughout long-term ecosystem development.

The importance of soil microbial phospholipid substitution and intracellular C storage for ecosystem function across long-term ecosystem development is supported by stoichiometric modelling (Fig. 5). Phospholipid substitution and intracellular C storage increased microbial C/P fourfold, on average, across the chronosequence (mean surface microbial C/P = 89 ± 15 versus 22 ± 9 as calculated without strategies; Fig. 5a) and microbial N/P almost twofold (surface microbial N/P = 11.7 ± 1.4 versus 6.3 ± 1.2 as calculated without substitution; Fig. 5b). These effect sizes were essentially constant across the chronosequence, and P efficiency strategies made no obvious contribution to the homoeostasis of microbial C/P and N/P ratios at Cooloola (evidenced by orthogonality of microbial versus substrate stoichiometry; Fig. 5c,d). However, the modelled potential for net microbial P mineralization increased threefold because of phospholipid substitution and intracellular C storage (Fig. 5f). Assuming moderate preferential uptake of P (α = 9, the mean modelled value across all dunes), microbial strategies for P efficiency bring microbes into the domain of C limitation and, thus, P cycling into the domain of net mineralization (Fig. 5f). Thus, we suggest that microbial P efficiency strategies sustain ecosystem function by pushing microbes close to absolute C limitation and, to an extent, stabilizing the strength of relative P versus N limitation for a given C supply (Fig. 4c,d), resulting in gradual but reliable recycling of P. These conditions should likewise ease competition between plants and soil microorganisms for P in P-depleted ecosystems. Such competition is increasingly recognized as an important driver of ecosystem responses to elevated atmospheric CO₂ concentrations^45,46. Thus, our insights into the physiological mechanisms that shape microbial P efficiency can support efforts to predict ecosystem C balance under future emissions scenarios.

Contemporary models of retrogression assert that declines in ecosystem productivity and function are driven by increasing P limitation as pedogenesis advances^3,31. Our study demonstrates that such models will be improved by considering soil microbial physiology. The proliferation of physiological strategies such as phospholipid substitution and intracellular C storage—which is probably underpinned by intraspecific plasticity and microbial taxonomic turnover across the chronosequence²⁶—could stabilize microbial productivity and function in the face of soil nutrient depletion. Moreover, by facilitating the fine-tuning of intracellular C and P content, these strategies could enhance the extent of microbial control over, and general temporal stability of, ecosystem P cycling throughout retrogression, enabling microbes to compete with plants for P while promoting gradual P mineralization and ultimately facilitating overall ecosystem productivity. By combining published vegetation biomass data for Cooloola⁴⁷ with P allometry data for relevant plant species, we estimate that more P is saved through the reduced contribution of phospholipids to total microbial P in Cooloola’s oldest dunes (relative to that in dune system one; 4.0 ± 0.5 kg P saved per hectare for 0–30 cm) than exists in the entire aboveground plant biomass on those dunes (~2.4 kg ha⁻¹). The former value underestimates the total quantity of P saved by microorganisms across the entire soil profile, as our ancillary data indicate meaningful, albeit somewhat reduced, quantities (1.8 mg kg^–1 of soil) of MBP at depths >30 cm in the terminal dune system at Cooloola. While there remain some uncertainties regarding the respective depth distributions of microbial biomass and P acquisitive roots, we therefore speculate that physiological microbial ‘P savings’ could contribute to the productivity and persistence of vegetation associated with well-weathered, nutrient-depleted soils. These communities, which include Mediterranean-climate shrublands and lowland tropical forests, contain some of the highest levels of floristic diversity and endemism worldwide^42,48,49. Thus, the ecophysiology of soil microorganisms is a potentially major driver of the long-term dynamics of terrestrial ecosystems and the biogeochemical function of the most biodiverse landscapes on Earth.

Study site and sample collection

The Cooloola sandmass is situated within the Great Sandy Region of coastal southeast Queensland, Australia (26° S, 153° E) and represents one of several similar dune fields that occur on Australia’s east coast from far north Queensland to New South Wales²⁸. The contemporary climate at Cooloola region is subtropical, with 1,343 mm of annual precipitation, a mean minimum temperature of 18.8 °C and a mean maximum temperature of 24.3 °C (ref. ⁵⁰). Vegetation at Cooloola is diverse and varies across the dune sequence but comprises primarily sclerophyllous shrubland, woodland or forest dominated by Myrtaceae and Proteaceae^47,51.

The sandmass contains as many as ten distinct stabilized parabolic dune systems, ranging in age from 0–0.04 ka to 350–730 ka (ref. ³⁰), with each system associated with a distinct period of active dune building via marine and/or aeolian deposition of sand^28,29. The sands at Cooloola are primarily quartz, containing <2% heavy minerals, feldspar, rutile, zircon, monzanite and ilmenite²⁸. Extensive weathering of the soils derived from these parent materials has led to substantial losses of mineral-derived nutrients, most notably P, throughout pedogenesis. We adhered to the six well-recognized dune systems described in foundational surveys of Cooloola⁵¹. We note here that chronosequence studies vary widely in terms of chronosequence length from, for example, 6 kyr in Arjeplog, Sweden, to 4.1 Myr across the Hawaiian islands, but generally have between 6–10 ‘stages’, with older stages becoming increasingly difficult to differentiate. Thus, chronosequences are often subject to a trade-off between temporal resolution number of stages per year) and long-term insight. At 350–730 kyr in duration, with ~6 stages, Cooloola optimizes this trade-off and thus provides one of the best systems for study of pedogenesis and primary succession worldwide. Nevertheless, the challenge of both extrapolating and interpolating our results should be considered when interpreting our findings.

We conducted the most rigorous sampling campaign possible within a single chronosequence by selecting replicate dunes for each dune system in a manner that minimizes non-independence between dune replicates. To the extent permitted by geomorphology, we spatially interspersed the replicate dunes of the six dune systems and/or maximized the distance among replicates of a given dune system (Fig. 1a).

Approximate sampling locations were selected on the basis of a combination of coordinates from historical sampling campaigns^30,47,51, existing mapping of the six geomorphologically recognized soil landscapes or dune systems³⁰ and Google Earth satellite imagery. Final specific locations were established in the field to ensure samples were collected from positions on trailing-arm dune crests close to dune apices (the youngest parts of the stabilized dunes) as these positions reduce the risks of upper soils being reworked by wind or overlain by sand from water erosion³⁰. Coordinates for each sampling location/replicate dune are provided in Extended Data Table 1.

To collect soil samples, overlying plant litter was first collected (and retained) using a 45 × 45 cm quadrat, and a 5 cm polyvinyl chloride tube was then gently inserted into the soil surface using a mallet. This 0–5 cm soil core was carefully extracted and transferred to a clip-seal bag. A 25 cm tube was then used to capture the 5–30 cm soil depth (Extended Data Fig. 4). Three to five sets of 0–5 cm and 5–30 cm cores were collected per site and bulked together to form two composite soil depth samples. At four sites representing one replicate each of dune systems 1, 2, 4 and 6, we collected an additional composite core from the 30–40 cm depth to make a preliminary evaluation of the persistence of MBP at depths >30 cm. Coring tubes were cleaned with water and then sterilized with 70% EtOH at each replicate dune. Root samples were subsequently obtained during the soil sieving stage (see the following) and were washed with ultrapure water before oven drying at 65 °C alongside the litter samples.

Samples were collected over three contiguous days in the late austral spring of 2022–2023 (22–24 November 2022), with dune systems stratified randomly across the sampling period to eliminate potential biases associated with sampling and storage. Weather conditions were consistent (dry and warm) across the sampling period. After sampling, samples were stored in plastic clip-seal bags with a small amount of ventilation and kept in the dark until they were brought to the laboratory for immediate processing and analysis, which occurred 3–5 days after sampling.

Laboratory analyses

Soil samples were either kept fresh under ambient conditions for rapid extraction and/or analysis of key biological properties (with prior sieving at 2 mm unless otherwise specified) or sieved at 2 mm, air-dried and finely ground before total elemental analyses. Soil gravimetric moisture content was determined immediately on return to the laboratory by oven drying soils at 105 °C until constant weight. Oven-dried litter and root samples were finely ground before elemental analyses.

Soil total C and N were determined via dry combustion (Elementar Vario MACRO cube). Soil total P was determined via molybdenum-blue colorimetry⁵² after digestion in concentrated H₂SO₄ with hydrogen peroxide^53,54. The same method was used for total P content of litter and roots. We used the Bray-2 extraction to derive an estimate of soil P that is likely to be potentially available for plant uptake. The Bray-2 P extraction quantifies PO₄³⁻–P in solution (immediately available P) plus PO₄³⁻–P that is weakly sorbed to soil minerals, and could be rendered available by root exudates, and involves a 60-second shaking extraction of fresh soil samples with a dilute NH₄-HCl solution, filtration with Whatman 42 filter paper and subsequent PO₄³⁻–P determination with molybdenum-blue colorimetry^53,55. For all PO₄³⁻–P determinations (including lipid P as described in the following), molybdenum-blue colorimetric analyses were carried out in 96-well plates using a Multi-Mode Microplate Reader (BioTek, Synergy 2) with linear calibration curves (0–1 ug ml^–1 P) constructed for each plate using the appropriate matrix²⁶. Before plating, sample pH was adjusted manually with H₂SO₄ and NaOH according to the endpoint of p-nitro-phenol (pH 7.5). Every sample was analysed in triplicate, and each plate included multiple checking standards. Reagent blanks were carried through all analyses, and where greater than zero, blank values were deducted from sample values.

Soil MBC, MBN and MBP were quantified on the basis of the conventional vacuum-infiltration CHCl₃ fumigation methods^56,57,58. For MBC and MBN, calculations invoked the standard correction factors of 2.64 and 2.22, while for MBP we derived a new correction factor (kP) based on the strong quadratic relationship between soil total P and empirically determined kP values as reported previously⁵⁹ (kP = 1.0238x2 – 1.5768x + 0.8649, where x = soil total P (g kg^–1 soil); R² = 0.85; n = 6). Modelled kP values ranged from 0.77 to 0.85 across our 48 samples. To avoid biasing relationships between final MBP values and soil total P values, we uniformly applied the overall mean kP factor of 0.83. MBP calculations further incorporated an empirically determined correction for PO₄³⁻ sorption during extraction of fumigated soils as per the original published protocol⁵⁶. Here a unique sorption correction value was determined for every individual soil sample (n = 48 + 4 samples from the 30–40 cm depth) based on spikes of ~10 µg PO₄³⁻–P per gram of soil. Rates of P recovery differed between the younger two dune systems (mean recovery = 39 ± 3%) and the older four dune systems (mean recovery = 91 ± 9%), but no other strong trends in PO₄³⁻ sorption were evident across the chronosequence.

The activities of soil extracellular enzymes (β-glucosidase (BG), N-acetyl-β-D-glucosaminidase (NAG), acid phosphomonoesterase (AP) and phosphodiesterase (PD)) were assayed using the p-nitrophenol spectrophotometric technique adapted for microplate analysis^60,61,62,63. During the assays, soils were incubated at 37 °C at the prescribed, notionally optimal pH (pH 6.0 for BG, pH 5.0 for NAG, pH 6.5 for AP and pH 8.0 for PD), such that measured enzyme activities represent potential, rather than actual, activities in soil. We chose not to measure leucine aminopeptidase (LAP) activity because the low and effectively uniform soil pH values at Cooloola (surface soil pH 5.96 ± 0.04) suggest that LAP probably makes a small (10–20%) and constant contribution to microbial N acquisition and NAG + LAP activity. Thus, we used activities of BG, NAG and AP + PD as indices of microbial C, N and P demand, respectively⁶⁴. With respect to enzyme analyses, there are a number of limitations that should be considered. First, we cannot discount a contribution of enzymes from plant roots. This would probably increase the signal of increasing belowground P limitation with dune age as species-level plant P content declines sharply in older dunes at Cooloola³², reflecting increasing P deficiency. Second, we did not measure the activity of phytases, which degrade inositol phosphates, a class of organic P compounds present in large quantities in many soils. However, recent evidence indicates that higher-order inositol phosphates could be largely absent from coastal dune soils in Australia⁶⁵. Third, certain assay conditions are somewhat arbitrary and based on convention (for example, the incubation temperature of 37 °C and the buffer pH values). This aids inter-site comparisons but may lead to biases if the respective enzymes have markedly different temperature responses or pH optima that change across the chronosequence. Finally, some conceptual issues surround the use of eco-enzymes as proxies for microbial resource demand⁶⁶, most notably that enzymes that notionally liberate one element (for example, P) also render C available. Such issues should be considered carefully when interpreting our results.

Soil lipids were extracted quantitatively using the Bligh and Dyer extraction⁶⁷ optimized for Australian soils⁶⁸. In brief, ~1 g fresh soil was extracted twice with a 2/1/0.8 mixture of methanol, chloroform (CHCl₃) and citric acid buffer (0.15 M, pH 4.0) in an ultrasonic bath. The organic CHCl₃ phase was subsequently separated, collected and stored at –20 °C until subsequent analysis. This phase was separated into three portions: one for determination of total lipid P, one for quantification of fatty acid methyl esters (FAMES) via gas chromatography-mass spectrometry (GC-MS) and one for quantification of intact lipids via liquid chromatography-mass spectrometry (LC-MS).

Total lipid P was determined by initially evaporating ~1 ml CHCl₃ under ambient conditions, digesting the dry lipid pellet in H₂SO₄ and hydrogen peroxide, and subsequently quantifying P with molybdenum-blue spectrophotometry as described in the preceding, with 5–6 analytical replicates per sample and checking standard P concentration = 0.05 µg ml⁻¹. To evaluate our lipid P measurements, we examined the relationship between surface soil lipid P content and (1) soil phospholipid content as determined via LC-MS and (2) soil polar lipid fatty acid content as determined by GC-MS (see the following for details of LC-MS and GC-MS analyses). As expected, assuming lipid P measurements accurately reflect the phosphate headgroups of membrane phospholipids, these relationships were strongly positive and linear, with R² = 0.87 for LC-MS phospholipids and R² = 0.73 for PLFAs.

For determination of FAMES, 1 ml of CHCl₃ extract was subjected to silica solid-phase extraction fractionation (Discovery DSC-SPE Tube bed wt 100 mg, volume 1 ml) to separate polar (methanol-eluted) membrane lipids from comparatively non-polar (CHCl₃-eluted) C storage lipids (TAGs). Methanol and CHCl₃ eluates were dried under nitrogen gas before and subsequently subjected to mild alkaline methanolysis (incubation at 37 °C for 30 minutes with 0.2 M methanolic KOH and CHCl₃). The methanolysed samples were again dried under nitrogen and taken up in hexane that included 25 μg ml⁻¹ methyl heptadecanoate-d₃₃ (Sigma-Aldrich) to serve as an internal standard. FAMES were then quantified via GC-MS (GCMSQP2010Plus, Shimadzu). The GC-MS was fitted with an arylene-modified 5% diphenyl–95% dimethyl polysiloxane stationary phase (30 m long × 0.25 mm internal diameter × 0.25 mm film thickness; Rxi-5SilMS, Restek) with helium as a carrier gas. We identified FAMES with reference to standards in commercial FAMES mixtures²⁶, as well as retention indices based on even-numbered alkanes from C8 through to C36 and EI mass spectra from the National Institute of Standards and Technology library. We followed standard taxonomic assignments for FAMES⁶⁹. PLFAs were used subsequently as an index of structural microbial biomass for analyses of microbial intracellular C storage allocation (see the following). Soil TAG content (mg C kg^–1 of soil) was calculated as the sum of all non-polar fatty acids from C13 to C24.

Determination of intact soil membrane lipids followed an established protocol^19,20,68. Aliquots (1 ml) of the CHCl₃ phase of initial soil lipid extracts were dried under nitrogen and reconstituted in a 2/1/1 mixture of isopropanol, acetonitrile and ultrapure water, with sphingosyl phosphatidylethanolamine (D17/1/12/0) as an internal standard. Samples were then subjected to nano-LC-MSⁿ following previous publications^20,68. In brief, samples were injected into a nano-scale analytical column (150 mm long × 75 μm internal diameter column, packed with Acclaim PepMap RSLC C18, 2 μm, 100 Å) and separated by gradient elution. Mobile phase A comprised 60% acetonitrile, 40% water, 10 mM ammonium formate and 0.1% formic acid, while mobile phase B comprised 90% isopropanol, 10% acetonitrile, 10 mM ammonium formate and 0.1% formic acid. Following nano-electrospray ionization, lipids were detected in positive mode Auto-MSⁿ and negative mode Auto-MSⁿ by mass spectrometer (AmaZon, Bruker Daltonics). After LC-MS, lipids were identified to class level on the basis of commercial standards, with reference to diagnostic neutral losses and product ions in either positive or negative mode, and subsequently to the level of sum composition. Quantification of each lipid species was achieved by dividing the lipid peak area by that of the internal standard, and subsequently dividing by an empirically determined response factor that captures the relationship between peak area and concentration for the lipid in question^19,20.

Polyhydroxybutyrate was analysed via GC-MS following a direct acid extraction–methylation protocol developed for similar sandy soils in eastern Australia²⁶. In this procedure, PHB is extracted, depolymerized and converted to its methylated monomer in a single step^70,71. Approximately 200 mg of fresh soil was digested at 100 °C for 3 hours with 800 μl CHCl₃ and 800 μl 5% H₂SO₄ in methanol containing 0.05 mg ml⁻¹ sodium benzoate as an internal standard. The digests were then phase separated with ultrapure water, and the organic phases were dried with solid Na₂SO₄ before analysis via GC-MS (with GC-MS instrument specifications as per FAMES analyses). PHB was quantified for all surface soil samples; however, due to instrumental constraints, PHB in subsoils (5–30 cm) was quantified only for dune systems one, four and six.

Microbial respiratory responses to factorial additions of C, N and P were tested for all replicates of dune systems two, four and six, as representatives of the young, middle-aged and old dune systems. Details of our incubation set-up and measurements are provided in the Supplementary Methods.

Calculations, statistical analyses and stoichiometric modelling

Enzyme vector analysis

Enzyme vector analyses followed previous studies, with the caveat that we measured a slightly different suite of enzymes than that used previously⁷².

Enzyme vector length—a unitless index of microbial C limitation—was calculated as follows:

Enzyme vector angle—an index of relative microbial N versus P limitation—was calculated as the inverse tangent, expressed in degrees, of the coordinates x (BG/BG + AP + PD) and y (BG/BG + NAG). Values of enzyme vector angle greater than 45° are associated with P limitation, and values less than 45° are associated with N limitation. Additional details of vector analyses are provided in the Supplementary Methods.

Statistical analysis

Statistical analyses were carried out in R (version 4.4.0). We tested the predictive significance dune age (the mean of dune age range estimates³⁰ on a natural log scale) and, separately, soil total P content for response variables using linear models or linear mixed-effects models using the ‘nlme’ package. Mixed-effects models were used when measurements were made at two soil depths (0–5 cm and 5–30 cm depth samples) to account for the non-independence. Accordingly, sampling site was used as a random intercept term in such models. Significances of dune system, sampling depth and the dune system × sampling depth interaction were evaluated with Wald tests. In the case of intracellular C storage compounds, where our primary interest was in allocation as opposed to absolute soil content, our linear mixed-effects models included an additional covariate term for microbial PLFA (in the case of TAG) and bacterial PLFA (in the case of PHB), in line with previous work on microbial storage allometry²⁶. As such, estimated marginal means presented reflect soil intracellular storage content at a constant value of microbial or bacterial PLFA.

To evaluate soil microbial respiratory responses to experimental C, N and P additions, we first calculated a set of response ratios to capture the effect of experimental treatments relative to untreated controls. Response ratios were calculated as the natural log of the ratio of experimental CO₂ efflux to corresponding control efflux. Next, we evaluated the significance of dune system, experimental treatment and the dune system × experimental treatment as predictors of response ratios across the 7-day incubation period using a linear mixed-effect model. Because CO₂ efflux was measured on multiple days (incubation days 1, 3, 5 and 7), we also included a random intercept term for incubation sample (sample jar). Statistical P values were obtained via Wald tests, and treatment means were compared using Tukey tests in the ‘emmeans’ package.

Stoichiometric modelling

Our stoichiometric modelling approach was based on previous modelling work^27,73 and underpinned by various arguments with which we link flows of C and P between soil and microbial biomass. Details of our modelling approach are provided in the Supplementary Methods.

Estimates of ecosystem-scale implications for microbial phospholipid substitution

Estimates of (1) the amount of microbial biomass that is supported by betaine lipid as a substitute for phospholipid and (2) microbial P savings that result from phospholipid substitution were based on the preservation of ratios of non-phospholipid to total membrane lipid from dune system one and lipid P to total microbial P, respectively. Aboveground vegetation P stock for dune system six was estimated on the basis of previously published vegetation biomass estimates⁴⁷ in combination with established allometric scaling relationships for similar taxa (conspecifics or congenerics) in similar environments^74,75,76. Details of these estimates and calculations are provided in the Supplementary Methods.