Ecology

Study area and field observations
We recorded the territorial calls of male savanna nightjars in eight areas of Taiwan, from north to south: Jinshan District (16 nightjars), Taichung City (7 nightjars), Hualien City (9 nightjars), Yuanlin City (8 nightjars), Beigang Township (8 nightjars), Chiayi City (14 nightjars), Taitung City (4 nightjars), and Hengchun Township (1 nightjar) (Fig. 1) during April-June of 2018. Sound recordings were collected in the downtown of each area using a Denon Portable IC Recorder (DN-F20R, sampling rate = 48 kHz, 16 bit, wav format) equipped with a Sennheiser ME67 unidirectional microphone. We made recordings between 19:00 to 24:00 in good weather during 1–2 consecutive nights for each area. If we recorded in the same area for more than one consecutive nights, the recording range of the second night was at least 1000 m away from the previous recording range. The maximum territory size of the savannah nightjar observed by Chan48 in urban areas of Taiwan was 83,424 m2 with a radius of about 163 m. Therefore, we are confident that we avoided recording the same territorial male twice. Territorial calls from an individual were recorded until the individual stopped calling or flew out of our recording range. When emitting territorial calls, nightjars either perched on some artificial structure, such as antennas or fences on roof tops, or they flew around the tops of buildings. Because the loud territorial calls recorded in these two situations demonstrate the same time–frequency patterns on spectrograms and sound the same when listened to, we did not differentiate between them while recording. We always attempted to record at the closest possible distance to the calling individual by moving closer to the individual at the street level.
After measuring the individual’s calls, we immediately took three samples of the maximum noise levels (maximum hold function, C-weighting function on Sound Level Meter TES-1350, TES, Taiwan) near the calling location and close to a road intersection within the calling range following the procedures detailed in Shieh et al.49. During each noise sampling, the sound level meter was held horizontally at a height of about 1.5 m and turned 360° clockwise to measure the maximum noise level from all directions within a time period of about 30 s. We then averaged the three samples of maximum noise measurements to arrive at our value of ‘ambient noise levels’ for each individual. Therefore, ambient noise levels were measured at a height of 1.5 m, which was not only the height at which anthropogenic sources such as traffic and human activity are the main sources of noise but also the height at which we held the microphone to record the nightjar calls emitted at a height of usually more than 10 m.
Playback-recording experiments
Seven artificial calls were generated using the frequency shift function (+ 300 Hz, + 200 Hz, + 100 Hz, 0 Hz, − 100 Hz, − 200 Hz, − 300 Hz) under the Frequency Domain Transformations tool in Avisoft-SASLab Pro software v5.2.12 from a source call with good quality (see Supplementary Table S6 for descriptions of acoustic measurements) and after band-pass filters (2–7 kHz) for noise removal. These seven artificial frequency-shifted calls represented seven different artificially created individuals and were identified by their frequency shift values (+ 300 Hz, + 200 Hz, + 100 Hz, 0 Hz, − 100 Hz, − 200 Hz, − 300 Hz). We then copied each frequency-shifted call 10 times with equal silent intervals and same amplitude as a group, and then we merged the seven groups of frequency-shifted calls into one sound file with alternated orders following a Latin square design. A total of seven merged sound files were obtained, and each was broadcast once in three sites of different urban noise levels: high, medium and low. The playback-recording experiments were conducted at three sites near or on the campus of Kaohsiung Medical University (22.648 N, 120.310 E). We also took 10 samples of the maximum noise level (maximum hold function, C-weighting function on Sound Level Meter TES-1350, TES, Taiwan) at the recording sites during the playback periods to obtain the ambient noise levels for each site. The first site was on the traffic roadside near a road intersection and had the highest ambient noise levels with a mean of 83.7 dB (n = 10) and a range of 79.3–89.6 dB. The second site was on the sport field of the campus about 30–50 m away from the traffic road and had the medium ambient noise levels with a mean of 74.6 (n = 10) and a range of 73.1–76.3 dB. The third site was on a 4th floor roof garden of the campus and had the lowest ambient noise levels with a mean of 71.4 dB (n = 10) and a range of 69.9–74 dB.
We used a Denon Portable IC Recorder (DN-F20R) connected to a speaker (Sony SRS-X11) for the playback experiments. The speaker was placed 1.5 m above the ground on a tripod. Playback-recording experiments were conducted between 17:00 and 18:30 h, a period with high traffic levels. The seven merged sound files were played back at a standardised volume with a sound-pressure level of 82.5 dB at 1 m from the speaker, and about 62.5 dB at 10 m from the speaker. The sound-pressure level of the broadcasting sound which we received at 10 m is about the same amplitude that we recorded a nightjar call at a distance of 28.8 m away from its calling spot with an amplitude of 97.7 dB. We recorded with a Denon Portable IC Recorder (DN-F20R) connected to a shotgun microphone (Sennheiser ME67) that was placed on a stand 1.5 m above the ground and oriented toward the broadcasting speaker at a distance of 10 m. All the recordings were set at the same recording levels and same settings (sampling rate = 48 kHz, 16 bit).
Sound analyses
We selected high-quality calls with clear acoustic structures on spectrograms and thus excluded any recordings with low-quality calls, which were, for example, emitted when the calling nightjar was too far away, or its calls were overlapped with calls from neighboring nightjars.
We also excluded any recordings where the individual only uttered one call (two calls being the minimum needed for inclusion). This left us with a total of 1925 calls from 67 individuals for our analyses, with a mean of 28.7 ± 3.1 calls/individual and a range of 2–97 calls/individual. The band-pass filters were set from 2.0–6.8 kHz and adjusted for each individual to reduce noise components. Using the recordings, we produced spectrograms with the following spectrogram parameters in the software: sampling frequency = 22.05 kHz, FFT = 512, hamming window, frequency resolution = 43 Hz, and time resolution = 2.9 ms. We quantified 30 acoustic variables (Table 1) from the spectrograms using the Automatic Parameter Measurements setup in the Avisoft-SASLab Pro software v5.2. We marked a call with the section label manually on the spectrogram by eye, and then two time-based parameters were measured (duration of the call and the temporal distance from the start to the location of the maximum amplitude) (see the manual50 of the Avisoft-SASLab Pro software for details). To automatically measure parameters other than the time-based on the labelled section, we specified four spectrum-based parameters (peak frequency, quartile 25%, quartile 50% and quartile 75%) to be measured at seven locations of the labelled section (start, end, maximum amplitude of the call, minimum parameter of entire call, maximum parameter of entire call, mean parameter of entire call, relative standard deviation of entire call); thus, 24 frequency-based parameters and four frequency-modulation-based parameters were automatically measured based on the labelled sections on the spectrograms.
For the playback-recording experiments, the recordings were first analyzed using the same settings as the above except for two differences. We adjusted the band-pass filters to 2.0–7.3 kHz because of the high frequency shift value (up to 300 Hz). Furthermore, we used the duration of the source call as the duration of all received calls; that is, we fixed the duration of the call section while marking, and the other 29 variables were automatically measured on the marked section by the software to reduce any possible human measurement errors51.
Statistical analyses
To examine possible geographic variation of the calls, we used individuals as our sample units (n = 67), and the averaged measurements from the calls of each individual were analyzed using a PCA. The PCA was performed on the 30 acoustic variables after a normalised transformation of each variable to a mean of zero and unit standard deviation (software PRIMER 6, version 6.1.5). We retained the five components with eigenvalues greater than one and interpreted each component based on its correlations with the original variables. However, because there is a sharp decrease of the eigenvalue and of the explained variance to smaller values from PC2 to PC3 (Supplementary Table S2), only PC1 and PC2 were used for examining the geographic patterns of the calls. The 95% confidence ellipses of the groups of individuals from different geographic areas are shown on the plot of PC1 against PC2 (Fig. 3). If the 95% confidence ellipses of the eight geographic areas overlapped, we can treat all the sampled individuals as one population and pool them for further analysis.
The following statistical tests were performed on the untransformed data of each acoustic variable using JMP Pro 14.2.0. First, we treated individuals as our sample units, and we used the averaged measurements from the calls of each individual. We performed Spearman rank tests to examine the relationships between ambient noise levels and each acoustic variable for 65 individuals (because noise measurements were not taken for two individuals). To distinguish variables’ relationship to ambient noise levels, we then classified the 30 acoustic variables into two categories: (1) noise-related variables had a significant relationship with ambient noise levels, and (2) noise-unrelated variables did not. Acoustic variables which had a significant (positive or negative) correlation with ambient noise levels using Spearman rank tests were classified as noise-related variables; all remaining variables lacking such a significant correlation were classified as noise-unrelated variables.
To determine variables which can distinguish the calls of different individuals, we treated calls as sample units and individuals as groups. We then performed a Kruskal–Wallis test to select variables which were more likely to encode individual information, that is, with significant individual differences. We only included those variables with a significant P-value into the DFA (see details below). To investigate how ambient noise levels affected the transmission efficacy of vocal individuality, we used DFA which calculates the accuracy of correctly assigning a particular call to a particular individual. The purpose of a DFA is thus to discriminate sampled individuals based on all the possible acoustic measurements which encode information about individual identity; therefore, we did not perform a variable selection process. Specifically, we used the accuracy value (1—misclassification rate) calculated from the DFA to assess the transmission efficacy of vocal individuality information. Higher accuracy values are assumed to indicate higher transmission efficacy of vocal individuality information from calls recorded through various ambient noise levels.
To compare possible differences in the transmission efficacy of vocal individuality for different sets of acoustic variables, we then calculated a separate DFA for three datasets: (1) all the 30 acoustic variables, (2) the noise-related variables, and (3) the noise-unrelated variables. For each dataset, a separate DFA was used to calculate one overall accuracy value and 67 individual accuracy values. The overall accuracy value (1—misclassification rate) describes the ability of the DFA to correctly assign the 1925 calls to the 67 recorded individuals. The individual accuracy values then describe the ability of the DFA to correctly assign the calls of one particular individual to that individual. To account for small sample sizes of calls by some individuals, the overall misclassification rates were bootstrapped with fractional weights option (number of bootstrap samples, n = 20,000), and the overall accuracy values of the three datasets and the associated 95% biased-corrected confidence intervals (BCI) were obtained. Sixty-seven individual accuracy values (corresponding to the 67 sampled nightjars) were also calculated for each set of variables.
To compare the individual accuracy values using the noise-related variables with the individual accuracy values using the noise-unrelated variables, we used the Wilcoxon signed rank test. To examine how noise levels affected the transmission efficacy through different acoustic structures (noise-related vs. noise-unrelated), we investigated the correlation between the individual accuracy value and the ambient noise level associated with that particular individual by using Spearman rank tests. The differences of the individual accuracy values, which was taken as the accuracy value of the DFA using the noise-unrelated variables minus the accuracy value of the DFA using the noise-related variables from the same individual, were correlated with the ambient noise levels by using Spearman rank tests to examine the similarity in trends.
For the playback-recording experiments, seven sound found files were played in each site (high, medium or low urban noise levels), and seven corresponding recordings were received as seven samples for each site. In each recording sample, although 10 calls of each individual (identified by frequency shift values) were played, the number of received calls might be less than 10 because the calls were overlapped with other unexpected sounds and thus excluded for measurements. Therefore, we first averaged the measurements of the possible received calls for each individual in a sample and then used the averaged measurements as the measurements of the sample. Thus, we obtained 49 samples (7 samples × 7 individuals) for each site for analyzing the overall accuracy of vocal individuality and variable accuracy for the playback-recording experiments. Since we excluded the duration variable, only 29 variables were used for DFA because the duration was fixed in measurements for all individuals. Furthermore, we only used six main frequency-based variables (PFSTART, PFEND, PFMIN, PFMAXA, PFMAX, PFMEAN) for further comparison on transmission accuracy between noise-related variables and noise-unrelated variables because the artificial calls were generated by transforming only the frequency domain. Thus, for each site, we calculated a separate DFA for three datasets: (1) all the 29 acoustic variables, (2) the three noise-related variables (PFSTART, PFEND, PFMIN), and (3) the three noise-unrelated variables (PFMAXA, PFMAX, PFMEAN). For the DFA, the overall misclassification rates were also bootstrapped (number of bootstrap samples, n = 20,000), and the overall accuracy values of the three datasets and the associated 95% biased-corrected confidence intervals (BCI) were obtained. Furthermore, for the six main frequency-based variables, an accuracy value for each variable was calculated as 1 − (|R – B| /B), in which R indicated the measurement value of the received call and B indicated the measurement value of the broadcast call. In each site for each variable, the accuracy values of the seven samples of the same individual were averaged as the variable accuracy value for the individual. In each site, we then investigated the differences of variable accuracy values among variables by Friedman rank tests (with individual as block) and Wilcoxon signed rank test for paired comparison between variables. We set the significance level at 0.05 for all tests and report the two-tailed probability values. More

1.
Parekh, P., Dutkiewicz, S., Follows, M. J. & Ito, T. Atmospheric carbon dioxide in a less dusty world. Geophys. Res. Lett. 33, L03610 (2006).
2.
Longhurst, A. R. In Ecological Geography of the Sea (Second Edition) (ed Longhurst, A. R.) 19– 34 (Academic Press, Burlington, 2007).

3.
Sverdrup, H. U. On conditions for the vernal blooming of phytoplankton. J. du Cons. Int. pour l’ Exploration de. la Mer. 18, 287–295 (1953).
Article  Google Scholar 

4.
Gran, H. H. & Braarud, T. A quantitative study of the phytoplankton in the Bay of Fundy and the Gulf of Maine (including observations on hydrography, chemistry and turbidity). J. Biol. Board Can. 1, 279–467 (1935).
CAS  Article  Google Scholar 

5.
Gran, H. H. Phytoplankton. methods and problems. J. du Cons. Int. Pour l’Exoploration de. la Mer. 7, 343–358 (1932).
Article  Google Scholar 

6.
Atkins, W. R. G. The chemistry of sea-water in relation to the productivity of the sea. Sci. Prog. 7, 298–312 (1932).
Google Scholar 

7.
Uchida, T. et al. Southern Ocean phytoplankton blooms observed by biogeochemical floats. J. Geophys. Res.: Oceans 124 (2019).

8.
Banse, K. In Primary Productivity and Biogeochemical Cycles in the Sea (eds Falkowski, P. G., Woodhead, A. & Vivirito, K.) vol. 43, 409–440 (Springer US, Boston, MA, 1992).

9.
Evans, G. T. & Parslow, J. S. A model of annual plankton cycles. Biol. Oceanogr. 3, 327–347 (1985).
Google Scholar 

10.
Behrenfeld, M. J. & Boss, E. S. Student’s tutorial on bloom hypotheses in the context of phytoplankton annual cycles. Glob. Change Biol. 24, 55–77 (2018).
ADS  Article  Google Scholar 

11.
Behrenfeld, M. J. Abandoning Sverdrup’s critical depth hypothesis on phytoplankton blooms. Ecology 91, 977–989 (2010).
PubMed  Article  Google Scholar 

12.
Boss, E. & Behrenfeld, M. In situ evaluation of the initiation of the North Atlantic phytoplankton bloom. Geophys. Res. Lett. https://doi.org/10.1029/2010GL044174. L18603 (2010).

13.
Mignot, A., Ferrari, R. & Claustre, H. Floats with bio-optical sensors reveal what processes trigger the North Atlantic bloom. Nat. Commun. 9, 190 (2018).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

14.
Westberry, T. K. et al. Annual cycles of phytoplankton biomass in the subarctic Atlantic and Pacific Ocean. Glob. Biogeochem. Cycles 30, 175–190 (2016).
ADS  CAS  Article  Google Scholar 

15.
Behrenfeld, M. J. et al. Annual boom-bust cycles of polar phytoplankton biomass revealed by space-based lidar. Nat. Geosci. 10, 118–122 (2017).
ADS  CAS  Article  Google Scholar 

16.
Arrigo, K. R., van Dijken, G. L. & Bushinsky, S. Primary production in the Southern Ocean, 1997-2006. J. Geophys. Res.: Oceans https://doi.org/10.1029/2007JC004551. C08004 (2008).

17.
Buitenhuis, E. T., Hashioka, T. & Quéré, C. L. Combined constraints on global ocean primary production using observations and models. Glob. Biogeochem. Cycles 27, 847–858 (2013).
ADS  CAS  Article  Google Scholar 

18.
Behrenfeld, M. J. Climate-mediated dance of the plankton. Nat. Clim. Change 4, 880–887 (2014).
ADS  Article  Google Scholar 

19.
Behrenfeld, M. J. & Milligan, A. J. Photophysiological expressions of iron stress in phytoplankton. Annu. Rev. Mar. Sci. 5, 217–246 (2013).
Article  Google Scholar 

20.
Banse, K. Grazing and zooplankton production as key controls of phytoplankton production in the open ocean. Oceanography 7, 13–20 (1994).
Article  Google Scholar 

21.
Behrenfeld, M. J., Doney, S. C., Lima, I., Boss, E. S. & Siegel, D. A. Annual cycles of ecological disturbance and recovery underlying the subarctic atlantic spring plankton bloom. Glob. Biogeochem. Cycles 27, 526–540 (2013).
ADS  CAS  Article  Google Scholar 

22.
Zehr, J. P. & Ward, B. B. Nitrogen cycling in the ocean: new perspectives on processes and paradigms. Appl. Environ. Microbiol. 68, 1015–1024 (2002).
CAS  PubMed  PubMed Central  Article  Google Scholar 

23.
Moore, C. M. et al. Processes and patterns of oceanic nutrient limitation. Nat. Geosci. http://www.nature.com/doifinder/10.1038/ngeo1765 (2013).

24.
Wong, A. P. S. & Riser, S. C. Profiling float observations of the upper ocean under sea ice off the wilkes land coast of Antarctica. J. Phys. Oceanogr. 41, 1102–1115 (2011).
ADS  Article  Google Scholar 

25.
Johnson, K. S. et al. Biogeochemical sensor performance in the SOCCOM profiling float array. J. Geophys. Res.: Oceans 122, 6416–6436 (2017).
ADS  Article  Google Scholar 

26.
Arrigo, K. R. et al. Massive phytoplankton blooms under Arctic sea ice. Science 336, 1408–1408 (2012).
ADS  CAS  PubMed  Article  Google Scholar 

27.
Arrigo, K. R. et al. Early spring phytoplankton dynamics in the western Antarctic Peninsula. J. Geophys. Res.: Oceans 122, 9350–9369 (2017).
ADS  Article  Google Scholar 

28.
Geider, R. J., Platt, T. & Raven, J. Size dependence of growth and photosynthesis in diatoms: a synthesis. Mar. Ecol. Prog. Ser. 30, 93–104 (1986).
ADS  CAS  Article  Google Scholar 

29.
Martin, J., Gordon, R. M. & Fitzwater, S. E. Iron in Antarctic waters. Nature 345, 156–158 (1990).

30.
Boyd, P. W. et al. Mesoscale iron enrichment experiments 1993-2005: Synthesis and future directions. Science 315, 612–617 (2007).
ADS  CAS  PubMed  Article  Google Scholar 

31.
Sarmiento, J. L. et al. Response of ocean ecosystems to climate warming. Glob. Biogeochem. Cycles 18, GB3003 (2004).
ADS  Article  CAS  Google Scholar 

32.
Riebesell, U., Körtzinger, A. & Oschlies, A. Sensitivities of marine carbon fluxes to ocean change. Proc. Natl Acad. Sci. USA 106, 20602–20609 (2009).
ADS  CAS  PubMed  Article  Google Scholar 

33.
Polovina, J. J., Howell, E. A. & Abecassis, M. Ocean’s least productive waters are expanding. Geophys. Res. Lett. 35, L03618 (2008).
ADS  Article  Google Scholar 

34.
Smith, M. J., Tittensor, D. P., Lyutsarev, V. & Murphy, E. Inferred support for disturbance-recovery hypothesis of north atlantic phytoplankton blooms. J. Geophys. Res.: Oceans 120, 7067–7090 (2015).
ADS  Article  Google Scholar 

35.
Yang, B. et al. Phytoplankton phenology in the North Atlantic: insights from profiling float measurements. Front. Mar. Sci. 7, 139 (2020).
ADS  Article  Google Scholar 

36.
Gruber, N. et al. Oceanic sources, sinks, and transport of atmospheric CO2. Glob. Biogeochem. Cycles 23 (2009).

37.
Gruber, N., Landschützer, P. & Lovenduski, N. S. The variable Southern Ocean carbon sink. Annu. Rev. Mar. Sci. 11, 159–186 (2019).
ADS  Article  Google Scholar 

38.
Sarmiento, J. L., Gruber, N., Brzezinski, M. A. & Dunne, J. P. High-latitude controls of thermocline nutrients and low latitude biological productivity. Nature 427, 56–60 (2004).
ADS  CAS  PubMed  Article  Google Scholar 

39.
Pahlow, M. & Prowe, F. Model of optimal current feeding in zooplankton. Mar. Ecol. Prog. Ser. 403, 129–144 (2010).
ADS  Article  Google Scholar 

40.
Johnson, K. S. et al. SOCCOM float data — Snapshot 2019-03-12. In Southern Ocean Carbon and Climate Observations and Modeling (SOCCOM) Float Data Archive (UC San Diego Library Digital Collections, 2019).

41.
Haëntjens, N., Boss, E. & Talley, L. D. Revisiting ocean color algorithms for chlorophyll a and particulate organic carbon in the Southern Ocean using biogeochemical floats. J. Geophys. Res.: Oceans https://doi.org/10.1002/2017JC012844 (2017).

42.
Graff, J. R. et al. Analytical phytoplankton carbon measurements spanning diverse ecosystems. Deep Sea Res. Part I: Oceanographic Res. Pap. 102, 16–25 (2015).
CAS  Article  Google Scholar 

43.
Behrenfeld, M. J., Doney, S. C., Lima, I., Boss, E. S. & Siegel, D. A. Reply to a comment by Stephen M. Chiswell on: Annual cycles of ecological disturbance and recovery underlying the subarctic Atlantic spring plankton bloom- by M. J. Behrenfeld et al. (2013). Glob. Biogeochemical Cycles 27, 1294–1296 (2013).
ADS  CAS  Article  Google Scholar 

44.
de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A. & Iudicone, D. Mixed layer depth over the global ocean: an examination of profile data and a profile-based climatology. J. Geophys. Res.: Oceans https://doi.org/10.1029/2004JC002378. C12003 (2004).

45.
Tagliabue, A. et al. A global compilation of dissolved iron measurements: focus on distributions and processes in the Southern Ocean. Biogeosciences 9, 2333–2349 (2012).
ADS  CAS  Article  Google Scholar 

46.
Westberry, T., Behrenfeld, M. J., Siegel, D. A. & Boss, E. Carbon-based primary productivity modeling with vertically resolved photoacclimation. Glob. Biogeochemical Cycles 22, GB2024 (2008).
ADS  Google Scholar 

47.
Ricchiazzi, P., Yang, S., Gautier, C. & Sowle, D. SBDART: a research and teaching software tool for plane-parallel radiative transfer in the Earth’s atmosphere. Bull. Am. Meteorological Soc. 79, 2101–2114 (1998).
ADS  Article  Google Scholar 

48.
Banse, K. Rates of phytoplankton cell division in the field and in iron enrichment experiments. Limnol. Oceanogr. 36, 1886–1898 (1991).
ADS  CAS  Article  Google Scholar 

49.
Behrenfeld, M. J., Boss, E., Siegel, D. A. & Shea, D. M. Carbon-based ocean productivity and phytoplankton physiology from space. Glob. Biogeochem. Cycles 19, GB1006 (2005).
ADS  Article  CAS  Google Scholar 

50.
Geider, R. J. Light and temperature dependence of the carbon to chlorophyll ratio in microalgae an cyanobacteria: Implications for physiology and growth of phytoplankton. N. Phytologist 106, 1–34 (1987).
CAS  Article  Google Scholar 

51.
Arteaga, L., Pahlow, M. & Oschlies, A. Global patterns of phytoplankton nutrient and light colimitation inferred from an optimality-based model. Glob. Biogeochem. Cycles 28, 648–661 (2014).
ADS  CAS  Article  Google Scholar 

52.
Geider, R. J. & LaRoche, J. The role of iron in phytoplankton photosynthesis, and the potential for iron-limitation of primary productivity in the sea. Photosynthesis Res. 39, 275–301 (1994).
CAS  Article  Google Scholar 

53.
Orsi, A. H., Whitworth, T. & Nowlin, W. D. On the meridional extent and fronts of the Antarctic Circumpolar Current. Deep Sea Res. Part I: Oceanographic Res. Pap. 42, 641–673 (1995).
ADS  Article  Google Scholar 

54.
Roemmich, D. & Gilson, J. The 2004–2008 mean and annual cycle of temperature, salinity, and steric height in the global ocean from the argo program. Prog. Oceanogr. 82, 81–100 (2009).
ADS  Article  Google Scholar 

55.
Bushinsky, S. M., Gray, A. R., Johnson, K. S. & Sarmiento, J. L. Oxygen in the Southern Ocean from argo floats: determination of processes driving air-sea fluxes. J. Geophys. Res.: Oceans 122, 8661–8682 (2017).
ADS  CAS  Article  Google Scholar 

56.
Arteaga, L. A., Pahlow, M., Bushinsky, S. M. & Sarmiento, J. L. Nutrient controls on export production in the Southern Ocean. Glob. Biogeochem. Cycles https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019GB006236 (2019), More

1.
Kasting, J. What caused the rise of atmospheric O2?. Chem. Geol. 362, 13–25 (2013).
ADS  CAS  Article  Google Scholar 
2.
Des Marais, D. J. Isotopic evolution of the biogeochemical carbon cycle during the Proterozoic Eon. Org. Geochem. 27(5–6), 185–193 (1997).
CAS  Article  Google Scholar 

3.
Karhu, J. A. & Holland, H. D. Carbon isotopes and the rise of atmospheric oxygen. Geology 24(10), 867–870 (1996).
ADS  CAS  Article  Google Scholar 

4.
Schidlowski, M. Carbon isotopes as biogeochemical recorders of life over 3.8 Ga of Earth history: evolution of a concept. Precambrian Res. 106(1–2), 117–134 (2001).
ADS  CAS  Article  Google Scholar 

5.
Aharon, P. Redox stratification and anoxia of the early Precambrian oceans: implications for carbon isotope excursions and oxidation events. Precambrian Res. 137, 207–222 (2005).
CAS  Google Scholar 

6.
Krissansen-Totton, J., Buick, R. & Catling, D. C. A statistical analysis of the carbon isotope record from the Archean to Phanerozoic and implications for the rise of oxygen. Am. J. Sci. 315(4), 275–316 (2015).
ADS  CAS  Article  Google Scholar 

7.
Martin, A. P., Condon, D. J., Prave, A. R. & Lepland, A. A review of temporal constraints for the Palaeoproterozoic large, positive carbonate carbon isotope excursion (the Lomagundi-Jatuli Event). Earth Sci. Rev. 127, 242–261 (2013).
ADS  CAS  Article  Google Scholar 

8.
Bekker, A. et al. Fractionation between inorganic and organic carbon during the Lomagundi (2.22–2.1 Ga) carbon isotope excursion. Earth Planet. Sci. Lett. 271(1–4), 278–291 (2008).
ADS  CAS  Article  Google Scholar 

9.
Maheshwari, A. et al. Global nature of the Paleoproterozoic Lomagundi carbon isotope excursion: a review of occurrences in Brazil, India, and Uruguay. Precambrian Res. 182(4), 274–299 (2010).
ADS  CAS  Article  Google Scholar 

10.
Melezhik, V. A., Huhma, H., Condon, D. J., Fallick, A. E. & Whitehouse, M. J. Temporal constraints on the Paleoproterozoic Lomagundi-Jatuli carbon isotopic event. Geology 35(7), 655–658 (2007).
ADS  CAS  Article  Google Scholar 

11.
Frauenstein, F., Veizer, J., Beukes, N., Van Niekerk, H. S. & Coetzee, L. L. Transvaal supergroup carbonates: implications for paleoproterozoic δ18O and δ13C records. Precambr. Res. 175, 149–160 (2009).
ADS  CAS  Article  Google Scholar 

12.
Hayes, J. M. & Waldbauer, J. R. The carbon cycle and associated redox processes through time. Philos. Trans. R. Soc. B 361, 931–950 (2006).
CAS  Article  Google Scholar 

13.
Frimmel, H. E. On the reliability of stable carbon isotopes for Neoproterozoic chemostratigraphic correlation. Precambrian Res. 182, 239–253 (2010).
ADS  CAS  Article  Google Scholar 

14.
Shields, G. A., Brasier, M. D., Stille, P. & Dorjnamjaa, D. I. Factors contributing to high δ13C values in Cryogenian limestones of western Mongolia. Earth Planet. Sci. Lett. 196(3–4), 99–111 (2002).
ADS  CAS  Article  Google Scholar 

15.
De PaulaSantos, G. M., Caetano-filho, S., Babinski, M. & Enzweiler, J. Rare elements of carbonate rocks from the Bambui Group, southern Sao Francisco Basin, Brasil, and their significance as paleoenvironmental proxies. Precambrian Res. 305, 327–340 (2017).
Article  CAS  Google Scholar 

16.
Klaebe, R. M., Kennedy, M. J., Jarrett, A. J. M. & Brocks, J. J. Local paleoenvironmental controls on the carbon-isotope record defining the Bitter Springs Anomaly. Geobiology 15(1), 65–80 (2017).
CAS  PubMed  Article  Google Scholar 

17.
Melezhik, V. A., Fallick, A. E., Medvedev, P. V. & Makarikhin, V. V. Extreme 13Ccarb enrichment in ca. 2.0 Ga magnesite-stromatolite-dolomite-red beds’ association in a global context: a case for the world-wide signal enhanced by a local environment. Earth-Sci. Rev. 48(1–2), 71–120 (1999).
ADS  CAS  Article  Google Scholar 

18.
Blättler, C. L. et al. Two-billion-year-old evaporites capture Earth’s great oxidation. Science 360(6386), 320–323 (2018).
PubMed  Article  CAS  Google Scholar 

19.
Hodgskiss, M. S., Crockford, P. W., Peng, Y., Wing, B. A. & Horner, T. J. A productivity collapse to end Earth’s Great Oxidation. Proc. Natl. Acad. Sci. 116(35), 17207–17212 (2019).
ADS  CAS  PubMed  Article  Google Scholar 

20.
Partin, C. A. et al. Uranium in iron formations and the rise of atmospheric oxygen. Chem. Geol. 362, 82–90 (2013).
ADS  CAS  Article  Google Scholar 

21.
Kanzaki, Y. & Murakami, T. Estimates of atmospheric O2 in the Paleoproterozoic from paleosols. Geochim. Cosmochim. Acta 174, 263–290 (2016).
ADS  CAS  Article  Google Scholar 

22.
Sheen, A. I. et al. A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia. Geochim. Cosmochim. Acta 227, 75–95 (2018).
ADS  CAS  Article  Google Scholar 

23.
Galili, N. et al. The geologic history of seawater oxygen isotopes from marine iron oxides. Science 365(6452), 469–473 (2019).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

24.
Knauth, L. P. Temperature and salinity of the Precambrian ocean: implications for the course of microbial evolution. Palaeogeogr. Palaeoclimatol. Palaeoecol. 219, 53–69 (2005).
Article  Google Scholar 

25.
Tartèse, R., Chaussidon, M., Gurenko, A., Delarue, F. & Robert, F. Warm Archaean oceans reconstructed from oxygen isotope composition of early-life remnants. Geochem. Perspect. Lett. 3, 55–65 (2017).
Article  Google Scholar 

26.
Kasting, J. Methane and climate during the Precambrian era. Precambr. Res. 137, 119–129 (2005).
ADS  CAS  Article  Google Scholar 

27.
Kasting, J. Early Earth: faint young Sun redux. Nature 464(7289), 687 (2010).
ADS  CAS  PubMed  Article  Google Scholar 

28.
Zinke, J., Reijmer, J. J. & Thomassin, B. Systems tracts sedimentology in the lagoon of Mayotte associated with the Holocene transgression. Sed. Geol. 160, 57–79 (2003).
CAS  Article  Google Scholar 

29.
Feuillet, N. MAYOBS1 Cruise, RV Marion Dufresne (Institut de Physique du Globe de Paris, 2019), https://doi.org/https://doi.org/10.17600/18001217

30.
Leboulanger, C. et al. Microbial diversity and cyanobacterial production in Dziani Dzaha crater lake, a unique tropical thalassohaline environment. PLoS ONE 12, e0168879 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

31.
Milesi, V. et al. Formation of Mg-smectite during lacustrine carbonates early diagenesis: study case of the volcanic crater lake Dziani Dzaha (Mayotte – Indian Ocean). Sedimentology (2018).

32.
Gérard, E. et al. Key role of alphaproteobacteria and cyanobacteria in the formation of stromatolites of Lake Dziani Dzaha (Mayotte, Western Indian Ocean). Front. Microbiol. 9, 1–20 (2018).
Article  Google Scholar 

33.
Cellamare, M. et al. Characterization of phototrophic microorganisms and description of new cyanobacteria isolated from the saline-alkaline crater-lake Dziani Dzaha (Mayotte, Indian Ocean). FEMS Microbiol. Ecol. 94(8), 1–25 (2018).
Article  CAS  Google Scholar 

34.
Hugoni, M. et al. Spatiotemporal variations in microbial diversity across the three domains of life in a tropical thalassohaline lake (Dziani Dzaha, Mayotte Island). Molecular Ecology (2018).

35.
Marty, B., Avice, G., Bekaert, D. V. & Broadley, M. W. Salinity of the Archaean oceans from analysis of fluid inclusions in quartz. Compte Rendus Geosci. 350(4), 154–163 (2018).
ADS  Article  Google Scholar 

36.
Hay, W. W. et al. Evaporites and the salinity of the ocean during the Phanerozoic: implications for climate ocean circulation and life. Palaeogeogr. Palaeoclimatol. Palaeoecol. 240(1–2), 3–46 (2006).
Article  Google Scholar 

37.
Marin-Carbonne, J., Chaussidon, M. & Robert, F. Micrometer-scale chemical and isotopic criteria (O and Si) on the origin and history of Precambrian cherts: Implications for paleo-temperature reconstructions. Geochim. Cosmochim. Acta 92, 129–147 (2012).
ADS  CAS  Article  Google Scholar 

38.
Marin-Carbonne, J., Robert, F. & Chaussidon, M. The silicon and oxygen isotope compositions of Precambrian cherts: a record of oceanic paleo-temperatures?. Precambr. Res. 247, 223–234 (2014).
ADS  CAS  Article  Google Scholar 

39.
Halevy, I. & Bachan, A. The geologic history of seawater pH. Science 355, 1069–1071 (2017).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

40.
Isson, T. T. & Planavsky, N. J. Reverse weathering as a long-term stabilizer of marine pH and planetary climate. Nature 560(7719), 471–475 (2018).
ADS  CAS  PubMed  Article  Google Scholar 

41.
Krissansen-Totton, J., Arney, G. N. & Catling, D. C. Constraining the climate and ocean pH of the early Earth with a geological carbon cycle model. Proc. Natl. Acad. Sci. 115(16), 4105–4110 (2018).
ADS  CAS  PubMed  Article  Google Scholar 

42.
Stüeken, E. E., Buick, R. & Schauer, A. J. Nitrogen isotope evidence for alkaline lakes on late Archean continents. Earth Planet. Sci. Lett. 411, 1–10 (2015).
ADS  Article  CAS  Google Scholar 

43.
Bartley, J. K. & Kah, L. C. Marine carbon reservoir, Corg-Ccarb coupling, and the evolution of the Proterozoic carbon cycle. Geology 32(2), 129–132 (2004).
ADS  CAS  Article  Google Scholar 

44.
Halevy, I., Alesker, M., Schuster, E. M., Popovitz-Biro, R. & Feldman, Y. A key role for green rust in the Precambrian oceans and the genesis of iron formations. Nat. Geosci. 10(2), 135–139 (2017).
ADS  CAS  Article  Google Scholar 

45.
Fakhraee, M., Hancisse, O., Canfield, D. E., Crowe, S. A. & Katsev, S. Proterozoic seawater sulfate scarcity and the evolution of ocean-atmosphere chemistry. Nat. Geosci. 12(5), 375–380 (2019).
ADS  CAS  Article  Google Scholar 

46.
Poulton, S. W. & Canfield, D. E. Ferruginous conditions: a dominant feature of the ocean through Earth’s history. Elements 7(2), 107–112 (2011).
CAS  Article  Google Scholar 

47.
Reinhard, C. T., Lalonde, S. V. & Lyons, T. W. Oxidative sulfide dissolution on the early Earth. Chem. Geol. 362, 44–55 (2013).
ADS  CAS  Article  Google Scholar 

48.
Och, L. M. & Shields-Zhou, G. A. The Neoproterozoic oxygenation event: environmental perturbations and biogeochemical cycling. Earth Sci. Rev. 110(1–4), 25–57 (2012).
ADS  Google Scholar 

49.
Planavsky, N. J., Bekker, A., Hofmann, A., Owens, J. D. & Lyons, T. W. Sulfur record of rising and falling marine oxygen and sulfate levels during the Lomagundi event. Proc. Natl. Acad. Sci. 109(45), 18300–18305 (2012).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

50.
Knoll, A. H., Bergmann, K. D. & Strauss, J. V. Life: the first two billion years. Philos. Trans. R. Soc. B Biol. Sci. 371, 1–13 (2016).
Article  Google Scholar 

51.
Butterfield, N. J. Early evolution of the Eukaryota. Palaeontology 58, 5–17 (2014).
Article  Google Scholar 

52.
Butterfield, N. J. Oxygen, animals and oceanic ventilation: an alternative view. Geobiology 7(1), 1–7 (2009).
MathSciNet  CAS  PubMed  Article  Google Scholar 

53.
Lenton, T. M., Boyle, R. A., Poulton, S. W., Shields-Zhou, G. A. & Butterfield, N. J. Co-evolution of eukaryotes and ocean oxygenation in the Neoproterozoic era. Nat. Geosci. 7(4), 257 (2014).
ADS  CAS  Article  Google Scholar 

54.
Peters, S. E., Husson, J. M. & Wilcots, J. The rise and fall of stromatolites in shallow marine environments. Geology 45(6), 487–490 (2017).
ADS  Article  Google Scholar 

55.
Gu, B., Schelske, C. L. & Hodell, D. A. Extreme 13C enrichments in a shallow hypereutrophic lake: implications for carbon cycling. Limnol. Oceanogr. 49, 1152–1159 (2004).
ADS  CAS  Article  Google Scholar 

56.
Zhu, Z., Chen, J. A. & Zeng, Y. Abnormal positive δ13C values of carbonates in lake Caohai, southwest China, and their possible relation to lower temperature. Quatern. Int. 288, 85–93 (2013).
Article  Google Scholar 

57.
Birgel, D. et al. Methanogenesis produces strong 13C enrichment in stromatolites of Lagoa Salgada, Brazil: a modern analogue for Palaeo- /Neoproterozoic stromatolites?. Geobiology 13, 245–266 (2015).
CAS  PubMed  Article  Google Scholar 

58.
Valero-Garcés, B. L., Delgado-Huertas, A., Ratto, N. & Navas, A. Large 13C enrichment in primary carbonates from Andean Altiplano lakes, northwest Argentina. Earth Planet. Sci. Lett. 171(2), 253–266 (1999).
ADS  Article  Google Scholar 

59.
Anoop, A. et al. Palaeoenvironmental implications of evaporative gaylussite crystals from Lonar Lake, central India. J. Quat. Sci. 28(4), 349–359 (2013).
Article  Google Scholar 

60.
Talbot, M. R. & Kelts, K. Primary and diagenetic carbonates in the anoxic sediments of Lake Bosumtwi, Ghana. Geology 14(11), 912–916 (1996).
ADS  Article  Google Scholar 

61.
Saba, V. S., Friedrichs, M. A., Antoine, D., Armstrong, R. A., Asanuma, I., Behrenfeld, M. J., Ciotti, A. M., Dowell, M., Hoepffner, N., Hyde, K. J. & Ishizaka, J. An evaluation of ocean color model estimates of marine primary productivity in coastal and pelagic regions across the globe. Biogeosciences. 2011, 489-503

62.
Lambrecht, N. et al. Biogeochemical and physical controls on methane fluxes from two ferruginous meromictic lakes. Geobiology 18(1), 54–69 (2020).
CAS  PubMed  Article  PubMed Central  Google Scholar 

63.
Bastviken, D., Cole, J., Pace, M. & Tranvik, L. Methane emissions from lakes: Dependence of lake characteristics, two regional assessments, and a global estimate. Global Biogeochem. Cycles 18(4), (2004)

64.
Bižić, M. et al. Aquatic and terrestrial cyanobacteria produce methane. Sci. Adv. 6(3), eaax5343 (2020).
ADS  PubMed  PubMed Central  Article  Google Scholar 

65.
Caetano-Filho, S., Sansjofre, P., Ader, M., Paula-Santos, G.M., Guacaneme, C., Babinski, M., Bedoya-Rueda, C., Kuchenbecker, M., Reis, H. L. & Trindade R. I. A large epeiric methanogenic Bambuì sea in the core of Gondwana supercontinent? Geosci. Front. (2020)

66.
Karl, D. M. & Knauer, G. A. Microbial production and particle flux in the upper 350 m of the Black Sea. Deep Sea Res. Part A Oceanogr. Res. Papers 38, S921–S942 (1991).
ADS  Article  Google Scholar 

67.
Katsev, S. & Crowe, S. A. Organic carbon burial efficiencies in sediments: the power law of mineralization revisited. Geology 43(7), 607–610 (2015).
ADS  CAS  Article  Google Scholar 

68.
Cowie, G. L., Hedges, J. I., Prahl, F. G. & De Lange, G. J. Elemental and major biochemical changes across an oxidation front in a relict turbidite: an oxygen effect. Geochim. Cosmochim. Acta 59(1), 33–46 (1995).
ADS  CAS  Article  Google Scholar 

69.
Logan, G. A., Hayes, J. M., Hieshima, G. B. & Summons, R. E. Terminal Proterozoic reorganization of biogeochemical cycles. Nature 376(6535), 53–56 (1995).
ADS  CAS  PubMed  Article  Google Scholar 

70.
Kuntz, L. B., Laakso, T. A., Schrag, D. P. & Crowe, S. A. Modeling the carbon cycle in Lake Matano. Geobiology 13(5), 454–461 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

71.
Laakso, T. A. & Schrag, D. P. Methane in the Precambrian atmosphere. Earth Planet. Sci. Lett. 522, 48–54 (2019).
ADS  CAS  Article  Google Scholar 

72.
Lambert, M. & Fréchette, J. L. Analytical techniques for measuring fluxes of CO2 and CH4 from hydroelectric reservoirs and natural water bodies. In Greenhouse Gas Emissions—Fluxes and Processes, Springer, Berlin, Heidelberg, 37–60 (2005).

73.
Abril, G. et al. Carbon dioxide and methane emissions and the carbon budget of a 10-year old tropical reservoir (Petit Saut, French Guiana). Global Biogeochem. Cycles 19(4), 1–16 (2005).
Article  CAS  Google Scholar 

74.
Assayag, N., Rivé, K., Ader, M., Jézéquel, D. & Agrinier, P. Improved method for isotopic and quantitative analysis of dissolved inorganic carbon in natural water samples. Rapid Commun. Mass Spectrom. 20(15), 2243–2251 (2006).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

75.
Lebeau, O., Busigny, V., Chaduteau, C. & Ader, M. Organic matter removal for the analysis of carbon and oxygen isotope compositions of siderite. Chem. Geol. 372, 54–61 (2014).
ADS  CAS  Article  Google Scholar 

76.
Galès, A., Triplet, S., Geoffroy, T., Roques, C., Carré, C., Le Floc’h, E., Lanfranchi, M., Simier, M., d’Orbcastel, E. R., Przybyla, C. & Fouilland, E. Control of the pH for marine microalgae polycultures: A key point of CO2 fixation improvement in intensive cultures. J. CO2 Util. 38, 187–193 (2020)

77.
Falkowski, P. G. & Raven, J. A. Aquatic Photosynthesis (Blackwell Science, Oxford, 1997).
Google Scholar 

78.
Silsbe, G. M. & Malkin, S. Y. Package “phytotools”: Phytoplankton Production Tools. CRAN library repository. https://cran.r-project.org/package=phytotools (2015).

79.
Eilers, P. H. C. & Peeters, J. C. H. A model for the relationship between light intensity and the rate of photosynthesis in phytoplankton. Ecol. Model. 42(3–4), 199–215 (1988).
Article  Google Scholar 

80.
Kirk, J. T. O. Light and Photosynthesis in Aquatic Environments 3rd edn. (Cambridge University Press, UK, 2010).
Google Scholar 

81.
Berner, R. A. Early Diagenesis: A Theoretical Approach (Princeton University Press, Princeton, 1980).
Google Scholar 

82.
Milesi, V. P. et al. Early diagenesis of lacustrine carbonates in volcanic settings: the role of magmatic CO2 (Lake Dziani Dzaha, Mayotte, Indian Ocean). ACS Earth Space Chem. 4(3), 363–378 (2020).
CAS  Article  Google Scholar  More

Affiliations

Laboratoire des Sciences de l’Environnement Marin (LEMAR), UMR 6539 UBO/CNRS/IRD/IFREMER, Institut Universitaire Européen de la Mer (IUEM), Brest, France
Hugo Berthelot, Stéphane L’Helguen, Jean-Francois Maguer & Nicolas Cassar

Division of Biology and Paleo Environment, Lamont-Doherty Earth Observatory, PO Box 1000, 61 Route 9W, Palisades, NY, 10964, USA
Solange Duhamel

Division of Earth and Ocean Sciences, Nicholas School of the Environment, Duke University, Durham, NC, 27708, USA
Seaver Wang & Nicolas Cassar

NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Code 616, Greenbelt, MD, USA
Ivona Cetinić

GESTAR/Universities Space Research Association, Columbia, MD, USA
Ivona Cetinić

Authors
Hugo Berthelot

Solange Duhamel

Stéphane L’Helguen

Jean-Francois Maguer

Seaver Wang

Ivona Cetinić

Nicolas Cassar

Corresponding authors
Correspondence to Hugo Berthelot or Nicolas Cassar. More

Our method is based on the assumption that the information structure of a system captures every effective interaction among its agents and thereby reflects the condition of the system. The abstract nature of information theory allows us to analyse systems independently of their specific instantiation. We only rely on the availability of longitudinal data reflecting the strength of the system’s individual components in a very broad sense. Hence, in general, our method can be applied to any complex system. The only condition is that for a given period of time and for every component of the system, a time series of quantitative data reflecting the outcome of interactions exists. Such time series could exemplarily be biomass of a plant species, number of individuals of an animal species, or sales of a company. The data type can differ among the components of the system, i.e. be heterogeneous.
In a first step, networks of information transfer are inferred via pairwise estimation of transfer entropy9 among all agents. Considering these networks and, in particular, their development over time, offers insights into functional interactions.
In the second step, potential, connectedness and resilience are computed solely using the networks of information transfer (see Supplementary A for a review of the adaptive cycle and its defining variables). Here, we utilize capacity and ascendency as being defined by Ulanowicz in the context of ascendency theory11. Note that Ulanowicz also used information theory to define capacity as an entropy of flows and ascendency as mutual information between inflow and outflow. While the first one is a measure of the average indeterminacy in the fluxes of the network, the latter quantifies the efficiency the system has in making use of its capacity12. However, being rooted in systems ecology, Ulanowicz always considered flows of physical quantities, such as energy or resource fluxes. In contrast, we will derive the quantities from networks of information transfer, abstracting from the physical representation of the interaction. Thus, potential is the capacity of the network of information transfer, and connectedness the corresponding ascendency.
The challenging part of our approach is to find an appropriate measure of resilience. There exist various conceptions and following definitions of resilience13. For our purposes, Holling’s definition of resilience fits best, namely to define resilience as ”the magnitude of disturbance that can be absorbed before the system changes the variables and processes that control behavior” (Ref.1, p. 28). There have been various approaches to make this notion measurable, however, all of them either depending on the specific system under observation14,15,16,17 or requiring deep knowledge of the system dynamics18. Resilience has also been studied from a network perspective (see e.g.19,20). Since we are modeling complex systems as networks of information transfer, our definition is inspired by a common concept in spectral graph theory. We use the so-called graph Laplace operator, which captures vulnerability of a network with respect to perturbations of its topological structure.
Taken together, the development of these three variables displays the system’s course through the adaptive cycle, helping to better understand system maturation and to evaluate its current condition. We will now provide a detailed description of our method, its implementation, and its application in the three case studies presented in this paper.
Step 1: estimation of networks of information transfer
Let (mathscr {V}) be a collection of variables, quantifying the state of agents defining a system. Let (I = (i_1, dots , i_N)) and (J = (j_1, dots , j_N)) be two sets of samples of states for the components I and J, say. For example, I and J can be identified with abundances of two interacting species at time points (1, dots , N). We consider the time series I and J as realisations of two approximately stationary discrete Markov processes. This allows us to compute Schreiber’s transfer entropy9, serving as a measure of their effective interaction. Transfer entropy from J to I is defined as

$$begin{aligned} T_{J rightarrow I} = sum _{n = 1}^{N-1} p left( i_{n+1},i_{n}, j_{n} right) cdot log left( frac{p left( i_{n+1}|i_{n}, j_{n} right) }{p left( i_{n+1} | i_{n} right) } right) . end{aligned}$$

(T_{J rightarrow I}) quantifies the average reduction in uncertainty about the future of I given the past of J. In other words, how much additional information do we gain about the next state of I, knowing not only the past of I itself, but the past of J as well. In the literature, a multitude of studies on the interpretation of transfer entropy in general and in specific contexts can be found21,22,23.
As the probabilities occurring in the definition of transfer entropy are in general not known, we have to estimate transfer entropy on the basis of given realizations of the random variables, i.e. the data given as samples of the time series. Typically, we do not use all available samples to estimate transfer entropy at time t but samples falling within a certain window of time preceding time t. The size (w_t) of this windows can either be fixed, or depend on the time t, e.g. (w_t = t). In the first case, the window is “shifted” going along with t to guarantee transfer entropy always being estimated on the same number of samples. In the second case, the window starts at the beginning of the time series and is extended with increasing t. In this case, the full history of the time series is considered for estimating transfer entropy. The choice of the window size depends on the system under consideration. In any case, it should be at least as large as the assumed order of the underlying Markov process. We then compute the information transfer from J to I at time t estimating transfer entropy over the period (t-w_t+1,dots ,t). To be precise,

$$begin{aligned} T_{J rightarrow I}^t = sum _{n = t-w_t+1}^{t} p left( i_{n+1},i_{n}, j_{n} right) cdot log left( frac{p left( i_{n+1}|i_{n}, j_{n} right) }{p left( i_{n+1} | i_{n} right) } right) . end{aligned}$$

Depending on the size of (w_t) and the data being available, it can be useful to increase the number of data points falling within every window by interpolation. For our calculations, we used the Matlab function pchip. Interpolation stabilizes the estimation in case of small window sizes. At the same time, interpolating too many points reduces stochasticity in the time series due to the deterministic component being introduced by the interpolation model. Thus, there is a trade-off between stochasticity and stability which has to be taken into account.
We estimated (T_{J rightarrow I}^t) using the Kraskov-Stögbauer-Grassberger (KSG) estimator TransferEntropyCalculatorKraskov as being provided with the JIDT toolkit24. For all our calculations, the function has been called using the data ((j_{t-w_t+1},dots ,j_t)) and ((i_{t-w_t+1},dots ,i_t)) in the mode computeAverageLocalOfObservations. For all other parameters of the estimation procedure, we used the default values (k=k_{tau}=l=l_{tau}=delay=1). Other choices of these parameters can be reasonable depending on the specific system to be analysed. To distinguish actual interactions from random noise, we tested all estimates via hundred-fold bootstrapping using the function computeSignificance(100) incorporated in the toolkit. Tests passing below a given significance level have been accepted and thus lead to an edge between the corresponding components with the estimated transfer entropy defining the corresponding weight. Estimating and testing for all pairs of components at fixed time t, we finally obtained a weighted, directed graph

$$begin{aligned} G^t = left( mathscr {V},{T_{J rightarrow I}^t|(J,I) in mathscr {V} times mathscr {V} } right) end{aligned}$$

as being our inferred model of interaction at time t. Given time series of abundances of length N for each component, this results in a sequence of interaction networks for time points (w_1,dots ,N).
Summarizing, the first step infers models of interaction among the given variables in form of a series of networks capturing the interaction patterns and their strengths. These network models can then be used in the second step to actually determine the position of the system within the adaptive cycle.
Step 2: determining potential, connectedness, and resilience
As mentioned before, our definitions of potential and connectedness are based on Ulanowicz’s notions of capacity and ascendency11. Ulanowicz provides further information on the theoretical background of these measures. Let

$$begin{aligned} T^t = sum _{(J,I) in mathscr {V} times mathscr {V}} T_{J rightarrow I}^t end{aligned}$$

be the total transfer of the system at time t. We further introduce the following shorthand notation

$$begin{aligned} T_{J}^{text {out},t} = sum _{I in mathscr {V}} T_{J rightarrow I}^t qquad text{ and } qquad T_{I}^{text {in},t} = sum _{J in mathscr {V}} T_{J rightarrow I}^t. end{aligned}$$

Define

$$begin{aligned} P^t = – sum _{(J,I) in mathscr {V} times mathscr {V}} T_{J rightarrow I}^t cdot log left( frac{T_{J rightarrow I}^t}{T^t} right) qquad hbox {as the system’s} ,potential ,hbox {at time} ,t end{aligned}$$

and

$$begin{aligned} C^t = sum _{(J,I) in mathscr {V} times mathscr {V}} T_{J rightarrow I}^t cdot log left( frac{T_{J rightarrow I}^tT^t}{T_{J}^{text {out,t}}T_{I}^{text {in,t}}} right) qquad ,hbox {as its} ,connectedness ,hbox {at time}, t. end{aligned}$$

Being essentially a sum over the indeterminacy in each transfer within the system, potential can be interpreted as a measure of the system’s power for evolution and its ability to develop. Recall that development of the system as a whole necessarily relies on communication, i.e., transfer of information among its components. In contrast, connectedness measures the degree of internal coherence of the system by contrasting information leaving one component with information arriving at another component.
In order to define resilience, we need to capture vulnerability of the system with respect to unforeseen perturbation. In terms of graph theory, this can be achieved by studying the eigenvalues of a certain matrix, being associated with the graph. Indeed, the smallest non-trivial eigenvalue of the so-called graph Laplacian of an undirected graph quantifies the vulnerability of the graph with respect to disturbance of the topology of the graph25,26. In our case, we need to transfer this idea to the case of the directed graphs (G^t).
Thus, given (G^t=left( mathscr {V},T^t_{J rightarrow I}|(J,I) in mathscr {V} times mathscr {V}right)) be a non-empty, weighted, directed graph with vertex set (mathscr {V}). Let further (c > 0) be a constant. Let (D_{out}) and (D_{in}) be the diagonal matrix of out-degrees and in-degrees, respectively, and A the weighted adjacency matrix. We then define the following Laplace type operators of (G^t):

$$begin{aligned} L_{out} = c cdot D^{-frac{1}{2}}_{out} left( D_{out}- A right) D^{-frac{1}{2}}_{out}, quad hbox { and }quad L_{in} = c cdot D^{-frac{1}{2}}_{in} left( D_{in} – A right) D^{-frac{1}{2}}_{in}, end{aligned}$$

following the convention that (D^{-frac{1}{2}}_*(u,u) = 0) for (D_*(u,u) = 0). Note that, for the sake of readability, we omitted the superscript t in these definitions. For all case studies presented in this paper, we used

$$begin{aligned} c = frac{1}{max { T^t_{J rightarrow I}|(J,I) in mathscr {V} times mathscr {V} }} end{aligned}$$

as standardization constant.
Since A is no longer symmetric, the spectrum of (L_{out}) and (L_{in}) is complex in general. Nevertheless, the distance of the spectrum to the imaginary axis in the complex plane still determines the stability of the graph. Therefore, we define resilience of the graph G as the smallest, non-trivial absolute value of the real parts of all eigenvalues of its two Laplacian matrices, i.e.

$$begin{aligned} R^t = min left{ |mathfrak {R}sigma | :sigma in {{,mathrm{Spec},}}(L_{out}) cup {{,mathrm{Spec},}}(L_{in}), sigma ne 0right} . end{aligned}$$

See Supplementary E for a more detailed explanation motivating this definition as well as for an alternative definition of the involved Laplacian matrices.
Our definitions of the three systemic variables are summarized in Table 1. In addition, Table 2 displays basic information concerning the data sets and parameters of our case studies. Note that, for visualization purposes, Figs. 3, 4, and 5 show a smoothed version of the estimated variables as being obtained by applying the R functions smooth.spline and splinefun.
Table 1 Summary of the definitions of the three systemic variables.
Full size table

Table 2 Data and parameters of the presented case studies.
Full size table

Figure 2

Schematic representation of our quantification method. In the first step, time series of abundance data (a) are transferred into networks of information transfer (b). In the second step, the three systemic variables (c) are computed on basis of the networks. The figure depicts the window shifting method.

Full size image

We created the R package QtAC in order to enable a straightforward application of our method27. The package comprises all functions required to compute a system’s course through the adaptive cycle and to visualize the results.
Figure 2 illustrates the key idea of our approach. Figure 2a shows randomly generated abundances of five components (A,B,C,D,E). To estimate the position of this small system within the adaptive cycle at time t and (t+1), we estimate transfer entropy for all pairs of components based on the samples within the window ((t-w+1, dots , t)) and ((t-w+2, dots , t+1)), respectively, and test for significance. This results in two inferred interaction networks shown in Fig. 2b. Using these networks, we can compute potential, connectedness and resilience at these two points in time. Figure 2c depicts the shift the system has made in the coordinate system spanned by the three characteristic variables.
The decrease in resilience from t to (t+1) mainly follows from loosing the edge (Drightarrow C) at (t+1). With component D being connected with the rest of the system by one edge, only, the system becomes more vulnerable, since perturbation of the edge (Drightarrow E) would fully decouple D from the rest of the system. Similarly, the loss of this edge also leads to a decrease of potential. Heuristically, the more edges a system has, the more potential there is to change from one state to another. Note that the even distribution of weights also added to the system’s potential, as for example edges (Arightarrow E) and (Arightarrow C) both loose weight. The moderate decline of connectedness follows from the loss of the edge (Drightarrow C) as well as from the smaller capacity of the edges (Arightarrow C) and (Arightarrow E), decreasing the overall total edge weight. More

Reagents and equipment
Para-nitrophenol, analytical grade; triethanolamine, analytical grade; calcium acetate (Ca(CO2CH3)2), analytical grade, China National Pharmaceutical Group Co., Ltd. Potassium chromate (K2CrO4), analytical grade; hydrochloric acid (HCl), analytical grade, Beijing Chemical Works. Calcium chloride (CaCl2·2H2O), analytical grade, Shanghai Macklin Biochemical Co., Ltd. Sodium hydroxide (NaOH), analytical grade, Xilong Chemical Industry Incorporated Co., Ltd. pH meter, Thermo Fisher Scientific, USA. ICP-OES 730, Agilent. X ray diffractometer (XRD), Bruker D8, Germany.
Soil conditioner: Developed by Tianjin Cement Industry Design and Research Institute Co., Ltd. Specific preparation method: The raw materials were K-feldspar (produced in Inner Mongolia, with ({mathrm{KAlSi}}_{3}{mathrm{O}}_{8}), ({mathrm{NaAlSi}}_{3}{mathrm{O}}_{8}), and SiO2 being the main components), CaCO3, and CaMg(CO3)2, which were crushed and ball-milled in order to get the sizes that could pass through an 80 μm sieve, before being mixed in appropriate ratios. Next, these were sintered by using an alumina crucible placed in a box furnace at 1270 °C for 60 min, then naturally cooled inside, and finally ground to approximately 80 μm to obtain the MSCs.
Soils samples
14 typical acid soils were listed in Table 1, which are from three provinces in China: Hunan (No. 1–5), Sichuan (No. 6–10), and Guangdong (No. 11–14). Soils were sampled in 0–20 cm, with all vegetation residues in the surface layer removed, then were placed indoors, naturally air-dried, and passed through a 2 mm sieve. A portion of each sample was used for testing and analysis, and the remained was used for the culture test.
Table 1 Information of the 14 types of acid soil samples in China.
Full size table

The evaluation of relationship between LR from schematics of ΔpH and soil acidity by using SMP-DB method
Mclean’s improved SMP-DB method was used to calculate LR. The calculation principle for the double buffer method is shown in Fig. 1, while the specific operating method and calculation principles of the experiment are as follows.
Figure 1

Computation of soil LR from the double buffer schematics and the relationships of the resulting similar triangles. (Notes: d = acidity of the soil neutralized by the buffer solution when the soil–buffer solution was at the ideal pH (6.5); d1 = acidity of the soil neutralized by the buffer solution when pH of the soil–buffer solution reduced from 7.5 to 1; d2 = acidity of the soil neutralized by the buffer solution when pH of the soil–buffer solution reduced from 6.0 to 2; pH1 = pH of the soil–buffer solution after addition of the SMP buffer solution; and pH2 = pH of the soil–buffer solution after addition of HCl.).

Full size image

Based on the schematics of the double buffer method and in accordance with the isosceles triangle principle, the proportional relationship was established as shown in Eq. (1):

$$frac{d-{d}_{2}}{{d}_{1}-{d}_{2}}=frac{6.5-{pH}_{2}}{{pH}_{1}-{pH}_{2}}$$
(1)

Equation (2) was obtained after conversion of Eq. (1):

$$d={d}_{2}+left({d}_{1}-{d}_{2}right)frac{6.5-{pH}_{2}}{{pH}_{1}-{pH}_{2}}$$
(2)

According to Fig. 1, Eqs. (3) and (4) were obtained by making (Delta {pH}_{1}=7.5-{pH}_{1}) and (Delta {pH}_{2}=6.0-{pH}_{2}), respectively. ({left(frac{Delta x}{Delta y}right)}^{^circ }) is the milligram equivalent (meq) of H+ that must be depleted to increase the pH of the buffer solution by one unit. It is determined based on the standard curve of the buffer solution.

$$frac{{d}_{1}}{Delta {pH}_{1}}={left(frac{Delta x}{Delta y}right)}^{^circ }$$
(3)

$$frac{{d}_{2}}{Delta {pH}_{2}}={left(frac{Delta x}{Delta y}right)}^{^circ }$$
(4)

Equations (5) and (6) were obtained after conversion of Eqs. (3) and (4):

$${d}_{1}=Delta {pH}_{1}times {left(frac{Delta x}{Delta y}right)}^{^circ }$$
(5)

$${d}_{2}=Delta {pH}_{2}times {left(frac{Delta x}{Delta y}right)}^{^circ }$$
(6)

These were substituted into Eq. (2) and then integrated to obtain Eq. (7):

$$mathrm{d}=Delta {pH}_{2}times {left(frac{Delta x}{Delta y}right)}^{^circ }+left(Delta {pH}_{1}-Delta {pH}_{2}right)times {left(frac{Delta x}{Delta y}right)}^{^circ }times frac{6.5-{pH}_{2}}{{pH}_{1}-{pH}_{2}}$$
(7)

Equation (7), intended for theoretical calculations, was derived from the mathematical relationships among the various parameters. (d) is the equivalent acidity of the soil neutralized by the buffer solution when the soil–buffer solution was at the ideal pH (6.5). The LR required to neutralize 5 g acid soil to pH 6.5 could then be extrapolated based on the measured data and the aforementioned equation.
Since the molecular weight of 1 mol of CaCO3 is 100, there is a clear conversion relationship between lime and CaCO3. For calculation convenience, LR is commonly expressed as the mass of CaCO3 needed to deplete the H+ present in 100 g of soil. In other words, ({mathrm{L}}_{R}= {text{meq CaCO}}_{3}/100 ,mathrm{g}=20mathrm{d}). After comparing the results obtained via the SMP-DB method and the Ca(OH)2-titrated acidity method, Mclean found that Eq. (8), a revision of the earlier equation, had a better correlation with the actual situation.

$${mathrm{L}}_{R}= {text{meq CaCO}}_{3}/100 ,mathrm{g}=1.69left(20dright)-0.86$$
(8)

Under normal circumstances, the weight of a 20 cm thick layer of ploughed soil would be 2250 t per hectare. For one hectare of soil, the LR of CaCO3 could be calculated using Eq. (9), with the unit being tons per hectare. This amount is expressed in meq CaCO3/100 g soil; to obtain approximate rates in metric tons per hectare (0–20 cm), it can be multiplied by 1.125.

$${L}_{R}= {text{meq CaCO}}_{3}/100 ,mathrm{g}times 1.125=left[1.69left(20dright)-0.86right]times 1.125=38.03d-0.97$$
(9)

SMP-DB buffer performance
Preparation of buffer solution27
800 mL distilled water was poured into a 1 L beaker, and then 1.8 g of para-nitrophenol, 2.5 mL of triethanolamine, 3.0 g of K2CrO4, 2.0 g of Ca(CO2CH3)2, and 53.1 g of CaCl2·2H2O were added into the beaker; the mixture was then stirred and mixed. NaOH 40% (w/w) or HCl 50% (v/v) was used to adjust the pH to 7.5 before the buffer solution was transferred to a 1 L volumetric flask. The beaker was rinsed for 3 times with 50 mL of distilled water, and the rinses were transferred to the volumetric flask. The eventual constant volume was 1 L.
Test method for buffer standard curve titration
1 mL of 0.05 M HCl was added into a 50 mL beaker with 10 mL of buffer solution in it and the pH of the solution was measured after stirring for 1 min. This operation was repeated 8 times and the corresponding pH was recorded to plot a titration curve.
The standard curve of the prepared buffer solution is shown in Fig. 2. It has a pH of 5–8 and its standard curve is linear, which could be fitted using a proportional function. The fitting yielded the linear equation y = −7.19x + 8.05, r2 = 0.998, which was highly significant. The buffering performance of the buffer solution was calculated based on the fitting equation, and the specific calculation process is stated below.
Figure 2

XRD pattern of the MSCs.

Full size image

Two points (x1, y1), (x2, y2) were selected and substituted into the fitting equation to obtain Eqs. (10) and (11):

$${y}_{1}=-7.19{x}_{1}+8.05$$
(10)

$${y}_{2}=-7.19{x}_{2}+8.05$$
(11)

The two equations were subtracted to obtain Eq. (12):

$${mathrm{y}}_{1}-{mathrm{y}}_{2}=-7.19left({x}_{1}-{x}_{2}right)$$
(12)

Let (Delta y={y}_{1}-{y}_{2}, Delta x=-left({x}_{1}-{x}_{2}right)), then Eq. (13) would be established.

$$frac{Delta y}{Delta x}=frac{{y}_{1}-{y}_{2}}{-left({x}_{1}-{x}_{2}right)}=7.19$$
(13)

When the pH increased by one unit, ∆ y = 1 was substituted into Eq. (13) to obtain Eq. (14):

$${left(frac{Delta x}{Delta y}right)}^{^circ }=frac{-left({x}_{1}-{x}_{2}right)}{{y}_{1}-{y}_{2}}=frac{1}{7.19}=0.139$$
(14)

In other words, 0.139 meq of H+ must be depleted per unit increase in the pH of the buffer solution.
During the process of the buffer titration test, the pH of the buffer was adjusted from 7.5 to 6.0. It was shaken again to determine pH2. Titration was performed using 0.05 M HCl (1 mL of 0.05 M HCl = 0.05 meq HCl). When the pH of the buffer solution was adjusted from 7.5 to 6.0, the reduction of 1.5 units required 0.139 × 1.5 = 0.2085 meq HCl, which converted to 4.2 mL of 0.05 M HCl.
Test procedure for buffer titration28
Soil pH was measured using a glass electrode pH meter. 5.00 g soil sample was weighed and placed in a 50 mL beaker. Deionized water was added at a 1:1 water-to-soil ratio, and the beaker was shaken for 10 min at 250 r min−1. After standing for 30 min, the pH (suspension) was measured. 10.00 mL of the SMP buffer solution was added and the mixture was shaken again for 10 min and then allowed to stand for 30 min. The suspension’s pH was measured to obtain pH1.
After measuring the pH and pH1, 4.2 mL of 0.05 M HCl was added to the suspension. This was the equivalent amount needed to adjust the buffer solution’s pH from 7.5 to 6.0 and was calculated according to the buffer solution’s standard curve. The mixture was shaken again for 10 min, and stand for 30 min before the pH of the soil suspension (pH2) was measured. The steps were repeated for 3 times.
ICP-OES measurement
The MSC main elemental contents were determined using ICP-OES. The operating parameters are presented in Table 2.
Table 2 Operating parameters of the ICP-OES spectrometer.
Full size table

0.200 g of each MSC was weighed and put into a 30 mL platinum crucible, and 1.500 g of molten agent was added (the mass ratio of sodium carbonate to sodium tetraborate was 2:1). After the molten agent and samples were mixed, the crucible was placed in a muffle furnace and its temperature was raised to 950 °C for 60 min to melt the contents. The crucible was taken out of the furnace after cooling and the sample inside was leached using 70 mL of HCl (3 + 7) to reach a constant volume of 100 mL. This solution was directly used to determine the Ca, Mg, Ba, Ti, and Mn content. Next, the solution was diluted 10 times to determine the high concentrations of K, Al, Si, Na, and Fe content.
XRD measurement
The MSC samples were ground to 0.045 mm (300 mesh) using an agate mortar and uniformly distributed inside sample frames. These were then pressed, flattened, and compacted using glass slides before being placed on the sample stage of the XRD sample chamber for analysis. The powder XRD patterns were obtained using a Bruker D8 Advance powder diffractometer working at 40 kV and 40 mA, using monochromatized Cu Kα radiation (λ = 0.154056 nm). The measurement was performed in the range angle 2θ = 15°–70°. Before the XRD test, all the samples were ground to 80 μm. The MDI Jade 5.0 software package (USA Materials Data Inc.) was used for qualitative analysis of the XRD spectra being tested.
Soil culture experiment
The soil samples were mixed with the MSCs and cultured for 30 days. Changes in the soil pH values were used to calculate the MSCs’ pH adjustment capacity and MSCR. The test treatments involved the addition of 0, 0.2%, 0.4%, 0.8%, 1.2%, or 1.6% of MSCs (total six levels) to the 14 soil samples and have 3 replications. The specific operating steps were as follows: 50 g of each soil sample and MSCs were mixed in each plastic cup uniformly, and water was added to 60% of the field moisture capacity. The cups were then sealed with plastic film to prevent excessive evaporation. The soil pH was measured after 30 days (1:1 water-to-soil ratio, measured after 10 min of shaking)29. The measurements were repeated twice. The measured pH and actual MSCR were subjected to regression analysis, and the regression equation was used to determine the MSCR required to neutralize the soil pH to 6.5 (the ideal pH for this study). More

Study sites and experimental design
The study sites are part of the NutNet experiment (Supplementary Data 1; http://nutnet.org/)27. Plots at each site are 5 × 5 m separated by at least 1 m. All sites included in the analyses presented here included unmanipulated plots and fertilized plots with nitrogen (N), phosphorus (P), and potassium and micronutrients (K) added in combination (NPK+). N, P, and K were applied annually before the beginning of the growing season at rates of 10 gm−2 y−1. N was supplied as time-release urea ((NH2)2CO) or ammonium nitrate (NH4NO3). P was supplied as triple super phosphate (Ca(H2PO4)2), and K as potassium sulfate (K2SO4). In addition, a micronutrient mix (Fe, S, Mg, Mn, Cu, Zn, B, and Mo) was applied at 100 gm−2 y−1 to the K-addition plots, once at the start of the experiment, but not in subsequent years to avoid toxicity. Treatments were randomly assigned to the 25 m2 plots and were replicated in three blocks at most sites (some sites had fewer/more blocks or were fully randomized). Sampling was done in 1 m2 subplots and followed a standardized protocol at all sites27.
Site selection
Data were retrieved on 1 May 2020. To keep a constant number of communities per site and treatment, we used three blocks per site, excluding additional blocks from sites that had more than three (Supplementary Data 1). Sites spanned a broad envelope of seasonal variation in precipitation and temperature (Supplementary Fig. 1), and represent a wide range of grassland types, including alpine, desert and semiarid grasslands, prairies, old fields, pastures, savanna, tundra, and shrub-steppe (Supplementary Data 1).
Stability and asynchrony measurements are sensitive to taxonomic inconsistencies. We adjusted the taxonomy to ensure consistent naming over time within sites. This was usually done by aggregating taxa at the genus level when individuals were not identified to species in all years. Taxa are however referred to as “species”.
We selected sites that had a minimum of 4 years, and up to 9 years of posttreatment data. Treatment application started at most sites in 2008, but some sites started later resulting in a lower number of sites with increasing duration of the study, from 42 sites with 4 years of posttreatment duration to 15 sites with 9 years of duration (Supplementary Data 1). Longer time series currently exist, but for a limited number of sites within our selection criteria.
Primary productivity and cover
We used aboveground live biomass as a measure of primary productivity, which is an effective estimator of aboveground net primary production in herbaceous vegetation36. Primary productivity was estimated annually by clipping at ground level all aboveground live biomass from two 0.1 m2 (10 × 100 cm) quadrats per subplot. For shrubs and subshrubs, leaves and current year’s woody growth were collected. Biomass was dried to constant mass at 60 °C and weighed to the nearest 0.01 g. Areal percent cover of each species was measured concurrently with primary productivity in one 1 × 1 m subplot, in which no destructive sampling occurred. Cover was visually estimated annually to the nearest percent independently for each species, so that total summed cover can exceed 100% for multilayer canopies. Cover and primary productivity were estimated twice during the year at some sites with strongly seasonal communities. This allowed to assemble a complete list of species and to follow management procedures typical of those sites. For those sites, the maximum cover of each species and total biomass were used in the analyses.
Diversity, asynchrony, and stability across spatial scales
We quantified local scale and larger-scale diversity indices across the three replicated 1-m2 subplots for each site, treatment and duration period using cover data37,38. In our analysis, we treated each subplot as a “community” and the collective subplots as the “larger scale” sensu Whittaker28. Local scale diversity indices (species richness, species evenness, Shannon, and Simpson) were measured for each community, and averaged across the three communities for each treatment at each site resulting in one single value per treatment and site. Species richness is the average number of plant species. Shannon is the average of Shannon–Weaver indices39. Species evenness is the average of the ratio of the Shannon–Weaver index and the natural logarithm of average species richness (i.e., Pielou’s evenness40). Simpson is the average of inverse Simpson indices41. Due to strong correlation between species richness and other common local diversity indices (Shannon: r = 0.90 (95% confidence intervals (CIs) = 0.87–0.92), Simpson: r = 0.88 (0.86–0.91), Pielou’s evenness: r = 0.62 (0.55–0.68), with d.f. = 324 for each), we used species richness as a single, general proxy for those variables in our models. Results using these diversity indices did not differ quantitatively from those presented in the main text using species richness (Supplementary Fig. 5), suggesting that fertilization modulate diversity effects largely through species richness. Following theoretical models15,16, we quantified abundance-based gamma diversity as the inverse Simpson index over the three subplots for each treatment at each site and abundance-based beta diversity, as the multiplicative partitioning of abundance-based gamma diversity: abundance-based beta equals the abundance-based gamma over Simpson28,42, resulting in one single beta diversity value per treatment and site. We used abundance-based beta diversity index because it is directly linked to ecosystem stability in theoretical models15,16, and thus directly comparable to theories. We used the R functions “diversity”, “specnumber”, and “vegdist” from the vegan package43 to calculate Shannon–Weaver, Simpson, and species richness indices within and across replicated plots.
Stability at multiple scales was determined both without detrending and after detrending data. For each species within communities, we detrended by using species-level linear models of percent cover over years. We used the residuals from each regression as detrended standard deviations to calculate detrended stability17. Results using detrended stability did not differ quantitatively from those presented in the main text without detrending. Stability was defined by the temporal invariability of biomass (for alpha and gamma stability) or cover (for species stability and species asynchrony), calculated as the ratio of temporal mean to standard deviation14,17. Gamma stability represents the temporal invariability of the total biomass of three plots with the same treatment, alpha stability represents the temporal invariability of community biomass averaged across three plots per treatment and per site, and species stability represents the temporal invariability of species cover averaged across all species and the three plots per treatment14. The mathematical formula are:

$${mathrm{Species}},{mathrm{stability}} = frac{{sum _{i,k}m_{i,k}}}{{sum _{i,k}sqrt {w_{ii,kk}} }},$$
(1)

$${mathrm{Alpha}},{mathrm{stability}} = frac{{sum _kmu _k}}{{sum _ksqrt {v_{kk}} }},$$
(2)

$${mathrm{Gamma}},{mathrm{stability}} = frac{{sum _kmu _k}}{{sqrt {sum _{k,l}nu _{kl}} }},$$
(3)

where mi,k and wii,kk denote the temporal mean and variance of the cover of species i in subplot k; μk and vkk denote the temporal mean and variance of community biomass in subplot k, and vkl denotes the covariance in community biomass between subplot k and l. We then define species asynchrony as the variance-weighted correlation across species, and spatial asynchrony as the variance-weighted correlation across plots:

$${mathrm{Species}},{mathrm{asynchrony}} = frac{{sum _{i,k}sqrt {w_{ii,kk}} }}{{sum _ksqrt {sum _{ij,kl}w_{ij,kl}} }},$$
(4)

$${mathrm{Spatial}},{mathrm{asynchrony}} = frac{{sum _ksqrt {v_{kk}} }}{{sqrt {sum _{k,l}nu _{kl}} }},$$
(5)

where wij,kl denotes the covariance in species cover between species i in subplot k and species j in subplot l.
These two asynchrony indices quantify the incoherence in the temporal dynamics of species cover and community biomass, respectively, which serve as scaling factors to link stability metrics across scales14 (Fig. 1). To improve normality, stability, and asynchrony measures were logarithm transformed before analyses. We used the R function “var.partition” to calculate asynchrony and stability across spatial scales14.
Climate data
Precipitation and temperature seasonality were estimated for each site, using the long-term coefficient of variation of precipitation (MAP_VAR) and temperature (MAT_VAR), respectively, derived from the WorldClim Global Climate database (version 1.4; http://www.worldclim.org/)44.
Analyses
All analyses were conducted in R 4.0.2 (ref. 45) with N = 42 for each analysis unless specified. First, we used analysis of variance to determine the effect of fertilization, and period of experimental duration on biodiversity and stability at the two scales investigated. Models including an autocorrelation structure with a first-order autoregressive model (AR(1)), where observations are expected to be correlated from 1 year to the next, gave substantial improvement in model fit when compared with models lacking autocorrelation structure. Second, we used bivariate analyses and linear models to test the effect of fertilization and period of experimental duration on biodiversity–stability relationships at the two scales investigated. Again, models including an autocorrelation structure gave substantial improvement in model fit (Supplementary Table 1)46,47,48. We ran similar models based on nutrient-induced changes in diversity, stability, and asynchrony. For each site, relative changes in biodiversity, stability, and asynchrony at the two scales considered were calculated, as the natural logarithm of the ratio between the variable in the fertilized and unmanipulated plots (Supplementary Fig. 9). Because plant diversity, asynchronous dynamics, and temporal stability may be jointly controlled by interannual climate variability22, we ran similar analyses on the residuals of models that included the coefficient of variation among years for each of temperature and precipitation. Results of our analyses controlling for interannual climate variability did not differ qualitatively from the results presented in the text (Supplementary Fig. 4). In addition, to test for temporal trends in stability and diversity responses to fertilization, we used data on overlapping intervals of four consecutive years. Results of our analyses using temporal trends did not differ qualitatively from the results presented in the text (Supplementary Fig. 6). Inference was based on 95% CIs.
Second, we used SEM29 with linear models, to evaluate multiple hypothesis related to key predictions from theories (Table 1). The path model shown in Fig. 1e was evaluated for each treatment (control and fertilized), and we ran separate SEMs for each period of experimental duration (from 4 to 9 years of duration). We generated a summary SEM by performing a meta-analysis of the standardized coefficients across all durations for each treatment. We then tested whether the path coefficients for each model differed by treatment by testing for a model-wide interaction with the “treatment” factor. A positive interaction for a given path implied that effects of one variable on the other are significantly different between fertilized and unfertilized treatments. We used the R functions “psem” to fit separate piecewise SEMs49 for each duration and combined the path coefficients from those models, using the “metagen” function50.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

Credit: Tad Theimer

Joseph (Joe) Connell altered both what and how ecologists study. Tree by tree, coral by coral, barnacle by barnacle, he saw patterns and processes across diverse ecosystems. Simply and with incontrovertible evidence, he demonstrated that interactions such as competition and predation could determine where species lived.
Before his classic experiments on Scotland’s rocky shores, field ecology was mainly descriptive, focusing on physical conditions such as temperature or moisture in determining where species lived. Connell, who died last month aged 96, inspired thousands of ecologists to test their hypotheses by manipulating conditions in the field.
Connell established long-term studies of coral reefs at Heron Island in the Great Barrier Reef and of tropical rainforests in Queensland, Australia, that spanned more than three and five decades, respectively. Monitoring revealed the dynamic nature of plant and animal communities that had long been considered stable. He discovered that natural variability in biological interactions and physical factors maintains diversity in these and other endangered ecosystems.
Born in 1923, Connell grew up just outside Pittsburgh, Pennsylvania. When the United States entered the Second World War in 1941, he enlisted in the Army Air Corps and was trained in meteorology. Later, conducting weather surveillance in the Azores — the Portuguese Atlantic archipelago — in support of army operations in Europe, he spent his free time birdwatching and identifying trees. Meeting army recruits who worked as wildlife managers, he realized it was possible to have a career as a biologist. After the war, and a degree in meteorology at the University of Chicago, Illinois, he headed to the University of California, Berkeley, for a master’s in zoology.
Connell produced what he described as a dull, unsatisfying thesis on brush rabbits (Sylvilagus bachmani) in the Berkeley Hills. Discouraged by the difficulties of conducting a population study (he trapped only 40 rabbits in 2 years), he adopted a rule of thumb — never again to study anything bigger than his thumb. As a doctoral student at the University of Glasgow, UK, he gleefully discovered what Charles Darwin had found a century before: that thousands of barnacles could easily be studied on the seashore, no traps required.
Connell realized that he could test his hypotheses about what factors determined where on the shore certain species lived by removing, adding or transplanting barnacles and their snail predators. Classic papers ensued, inspiring other ecologists to rethink distribution patterns, and, importantly, to test their ideas with controlled field experiments.
After a postdoc at the Woods Hole Oceanographic Institution in Massachusetts, Connell joined the faculty at the University of California, Santa Barbara, where he remained for the rest of his career. He was curious about processes that affected distribution and abundance, and those that might keep biodiversity high. Shifting to species that live for hundreds or thousands of years on coral reefs and in rainforests, he set up his Australian long-term monitoring studies in 1962 and 1963. Both recorded the demography and interactions of organisms in permanent plots, tracking community dynamics and the impact of disturbances, ranging from fallen trees to cyclones.
Visiting Connell’s sites with him in the 1970s and 1990s, we were impressed with his foresight and inspired by his insights. On the reef, he explained, physical disturbance by large waves associated with recurring cyclones intermittently reduced the cover of dominant species such as staghorn coral (Acropora aspera). This prompted recolonization by a diverse assemblage of weaker competitors such as encrusting or mound-like species. Connell coined the term ‘intermediate disturbance hypothesis’ to describe this process.
We strolled through the larger of his rainforest plots, avoiding stinging trees, biting flies, ticks and leeches, and relishing the richness — more than 300 tree species and about 100,000 individual plants. Connell outlined another hypothesis, that forests are more diverse when rarer species such as the conifer Sundacarpus amarus are favoured over more common ones such as the flowering tree Planchonella sp. Patterns of seedling establishment, growth or survival depend on that difference in frequency. Because common species grow more densely than rare ones, they are more vulnerable to specialist herbivores or pathogens.
This pattern of density-dependent predation or infection thins out common species, enabling a richer mix to coexist. It is a central component of the Janzen–Connell hypothesis (independently proposed by US ecologist Daniel Janzen in 1970), which predicts that seedlings are more likely to die under the canopies of their parent trees than farther away, ensuring diversity.
Connell was unfailingly kind, generous and devoted to his family. He never lost his profound curiosity about the natural world or his delight in exploring ideas with students and colleagues. He loved to be challenged and, if proven wrong, he gladly moved on to a new hypothesis or question. He sought truth, not fame. Moreover, he empowered everyone around him to think critically by focusing on ideas and evidence, not personalities. Fortunately for the world, his way of exploring science proved powerful, infectious, fun and enduringly productive. More