Ecology

United Nations. Transforming Our World: The 2030 Agenda for Sustainable Development (United Nations, 2015).
Google Scholar 
Salimpour, S., Khavazi, K., Nadian, H., Besharati, H. & Miransari, M. Enhancing phosphorous availability to canola (Brassica napus L.) using P solubilizing and sulfur oxidizing bacteria. Plant Biol. 6, 629–642 (2010).
Google Scholar 
Ezawa, T., Smith, S. E. & Smith, F. A. P metabolism and transport in AM fungi. Plant Soil 244, 221–230 (2002).CAS 
Article 

Google Scholar 
Halajnia, A., Haghnia, G. H., Fotovat, A. & Khorasani, R. Phosphorus fractions in calcareous soils amended with P fertilizer and cattle manure. Geoderma 150, 209–213 (2009).ADS 
CAS 
Article 

Google Scholar 
Adnan, M. et al. Coupling phosphate-solubilizing bacteria with phosphorus supplements improve maize phosphorus acquisition and growth under lime induced salinity stress. Plants 9, 900 (2020).CAS 
PubMed Central 
Article 

Google Scholar 
Khan, A. A., Jilani, G., Akhtar, M. S., Naqvi, S. M. S. & Rasheed, M. Phosphorus solubilizing bacteria, occurrence, mechanisms and their role in crop production. J. Agric. Biol. Sci. 1, 48–58 (2009).
Google Scholar 
Torrent, J., Barron, V. & Schwertmann, U. Phosphate adsorption and desorption by goethites differing in crystal morphology. Soil Sci. Soc. Am. J. 54, 1007–1012 (1990).ADS 
Article 

Google Scholar 
Rehim, A. Band-application of phosphorus with farm manure improves phosphorus use efficiency, productivity, and net returns of wheat on sandy clay loam soil. Turk. J. Agric. For. 40, 319–326 (2016).CAS 
Article 

Google Scholar 
Bieleski, R. L. Phosphate pools, phosphate transport and phosphate availability. Annu. Rev. Plant Physiol. 24, 225–252 (1973).CAS 
Article 

Google Scholar 
Goldstein, A. H. Recent progress in understanding the molecular genetics and biochemistry of calcium phosphate solubilization by gram negative bacteria. Biol. Agric. Hortic. 12, 185–193 (1995).Article 

Google Scholar 
Lopez-Bucio, J., Vega, O. M., Guevara-Garcıa, A. & Herrera-Estrella, L. Enhanced phosphorus uptake in transgenic tobacco plants that overproduce citrate. Nat. Biotechnol. 18, 450–453 (2000).CAS 
PubMed 
Article 

Google Scholar 
Tilman, D. et al. Forecasting agriculturally driven global environmental change. Science 292, 281–284 (2001).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Sato, S., Solomon, D., Hyl, C., Ketterings, Q. M. & Lehmann, J. Phosphorus speciation in manure and manure-amended soils using XANES spectroscopy. Environ. Sci. Technol. 39, 7485–74919 (2000).ADS 
Article 
CAS 

Google Scholar 
Brady, N. C., Weil, R. R. & Weil, R. R. The Nature and Properties of Soils Vol. 13, 662–710 (Prentice Hall, 2008).
Google Scholar 
Adnan, M. et al. Coupling phosphate solubilizing bacteria with Phosphorus supplements improve maize phosphorus acquisition and growth under lime induced salinity stress. Plants 9, 900 (2020).CAS 
PubMed Central 
Article 

Google Scholar 
Caravaca, F., Alguacil, M. M., Azcon, R., Diaz, G. & Roldan, A. Comparing the effectiveness of mycorrhizal inoculum and amendment with sugar beet, rock phosphate and Aspergillus niger to enhance field performance of the leguminous shrub Dorycnium pentaphyllum L.. Appl. Soil Ecol. 25, 169–180 (2004).Article 

Google Scholar 
Zaidi, A., Khan, M., Ahemad, M. S., Oves, M. & Wani, P. A. Recent advances in plant growth promotion by phosphate-solubilizing microbes. In Microbial Strategies for Crop Improvement (eds Khan, M. S. et al.) 23–50 (Springer, 2009).Chapter 

Google Scholar 
Illmer, P., Barbato, A. & Schinner, F. Solubilization of hardly-soluble AlPO4 with P-solubilizing microorganism. Soil Biol. Biochem. 27, 265–270 (1995).CAS 
Article 

Google Scholar 
Ryan, P. R., Delhaize, E. & Jones, D. L. Function and mechanism of organic anion exudation from plant roots. Annu. Rev. Plant Biol. 52, 527–560 (2001).CAS 
Article 

Google Scholar 
Chen, Y. P. et al. Phosphate solubilizing bacteria from subtropical soil and their tricalcium phosphate solubilizing abilities. Appl. Soil Ecol. 34, 33–41 (2006).Article 

Google Scholar 
Adnan, M. et al. Integration of poultry manure and phosphate solubilizing bacteria improved availability of Ca bound P in calcareous soils. 3 Biotech 9, 368 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
He, Z. & Zhu, J. Microbial utilization and transformation of phosphate adsorbed by variable charged minerals. Soil Biol. Biochem. 30, 917–923 (1988).Article 

Google Scholar 
Kucey, R. M. N. Effect of Penicillium bilajion the solubility and uptake of P and micronutrients from soil by wheat. Can. J. Soil Sci. 68, 261–270 (1988).CAS 
Article 

Google Scholar 
Bünemann, E. K., Bossio, D. A., Smithson, P. C., Frossard, E. & Oberson, A. Microbial community composition and substrate use in a highly weathered soil as affected by crop rotation and P fertilization. Soil Biol. Biochem. 36, 889–901 (2004).Article 
CAS 

Google Scholar 
McGill, W. B. & Cole, C. V. Comparative aspects of cycling of organic C, N, S and P through soil organic matter. Geoderma 26, 267–268 (1981).ADS 
CAS 
Article 

Google Scholar 
Chaiharn, M. & Lumyong, S. Screening and optimization of indole-3-acetic acid production and phosphate solubilization from rhizobacteria aimed at improving plant growth. Curr. Microbiol. 62, 173–181 (2011).CAS 
PubMed 
Article 

Google Scholar 
Kucey, R. M. N., Janzen, H. H. & Legett, M. E. Microbially mediated increases in plant-available phosphorus. Adv. Agron. 42, 198–228 (1989).
Google Scholar 
Rodriguez, H. & Fraga, R. Phosphate solubilizing bacteria and their role in plant growth promotion. Biotechnol. Adv. 17, 319–339 (1999).CAS 
PubMed 
Article 

Google Scholar 
Xiao, Y., Wang, X., Chen, W. & Huang, Q. Isolation and identification of three potassium-solubilizing bacteria from rape rhizospheric soil and their effects on ryegrass. Geomicrobiol. J. 34, 873–880 (2017).CAS 
Article 

Google Scholar 
Sugihara, S., Funakawa, S., Kilasara, M. & Kosaki, T. Dynamics of microbial biomass nitrogen in relation to plant nitrogen uptake during the crop growth period in a dry tropical cropland in Tanzania. Soil Sci. Plant Nutr. 56, 105–114 (2010).CAS 
Article 

Google Scholar 
Jalili, F. et al. Isolation and characterization of ACC deaminase producing fluorescent pseudomonads, to alleviate salinity stress on canola (Brassica napus L.) growth. J. Plant Physiol. 166, 667–674 (2009).CAS 
PubMed 
Article 

Google Scholar 
Tiwari, V. N., Lehri, L. K. & Pathak, A. N. Effect of inoculating crops with phospho-microbes. Exp. Agric. 25, 47–50 (1989).Article 

Google Scholar 
Pal, S. S. Interaction of an acid tolerant strain of phosphate solubilizing bacteria with a few acid tolerant crops. Plant Soil 213, 221–230 (1999).MathSciNet 
Article 

Google Scholar 
Afzal, A., Ashraf, M., Asad, S. A. & Faroog, M. Effect of phosphate solubilizing microorganism on phosphorus uptake, yield and yield traits of wheat (Triticum aestivum L.) in rainfed area. Int. J. Agric. Biol. 7, 207–209 (2005).
Google Scholar 
Bolan, N. S., Naidu, R., Mahimairajaand, S. & Baskaran, S. Influence of low-molecular-weight organic acids on the solubilization of phosphates. Biol. Fertil. Soils 18, 311–319 (1994).CAS 
Article 

Google Scholar 
Mihoub, A., Amin, A. E. E. A. Z., Motaghian, H. R., Saeed, M. F. & Naeem, A. Citric acid (CA)–modified biochar improved available phosphorus concentration and its half-life in a P-fertilized calcareous sandy soil. J. Soil Sci. Plant Nutr. 22(1), 465–474 (2022).CAS 
Article 

Google Scholar 
Adnan, M., Shah, Z., Sharif, M. & Rahman, H. Liming induces carbon dioxide (CO2) emission in PSB inoculated alkaline soil supplemented with different phosphorus sources. Environ. Sci. Pollut. Res. 25(10), 9501–9509 (2018).CAS 
Article 

Google Scholar 
Amin, A. E. E. A. Z. & Mihoub, A. Effect of sulfur-enriched biochar in combination with sulfur-oxidizing bacterium (Thiobacillus spp.) on release and distribution of phosphorus in high calcareous p-fixing soils. J. Soil Sci. Plant Nutr. 21(3), 2041–2047 (2021).CAS 
Article 

Google Scholar 
Tawaraya, K., Hirose, R. & Wagatsuma, T. Inoculation of arbuscularmycorrhizal fungi can substantially reduce phosphate fertilizer application to Alliumfis-tulosum L. and achieve marketable yield underfield condition. Biol. Fertil. Soils 48, 839–843 (2012).Article 

Google Scholar 
Islam, M. T. & Hossain, M. M. Plant probiotics in phosphorus nutrition in crops, with special reference to rice. In Bacteria in Agrobiology, Plant Probiotics (ed. Maheshwari, D. K.) 325–363 (Springer, 2012).Chapter 

Google Scholar 
Amruthesh, K. N., Raj, S. N., Kiran, B., Shetty, H. S. & Reddy, M. S. Growth promotion by plant growth-promoting rhizobacteria in some economically important crop plants. In Sixth International PGPR Workshop, 5–10 October, Calicut, India, 97–103 (2003).Kumar, S. et al. Impacts of nitrogen rate and landscape position on soils and switchgrass root growth parameters. Agron. J. 111, 1046–1059 (2019).CAS 
Article 

Google Scholar 
Mihoub, A. & Boukhalfa-Deraoui, N. Performance of different phosphorus fertilizer types on wheat grown in calcareous sandy soil of El-Menia, Southern Algeria. Asian J. Crop Sci. 6, 383–391 (2014).Article 

Google Scholar 
Piccini, D. & Azcon, R. Effect of phosphate solubilizing bacteria and vesicular-arbuscular mycorrhizal fungi on the utilization of Bayovar rock phosphate by alfalfa plants using a sand-vermiculite medium. Plant Soil 50, 45–50 (1987).Article 

Google Scholar 
Dwivedi, B. S., Singh, V. K. & Dwivedi, V. Application of phosphate rock, with or without Aspergillus awamori inoculation, to meet phosphorus demands of rice–wheat systems in the Indo Gangetic plains of India. Aus. J. Exp. Agric. 44, 1041–1050 (2004).CAS 
Article 

Google Scholar 
Saad, O. A. O. & Hammad, A. M. M. Fertilizing wheat plants with rock phosphate combined with phosphate dissolving bacteria and V.A mycorrhiza as alternate for ca–superphosphate. Ann. Agric. Sci. Cairo 43, 445–460 (1998).
Google Scholar 
Chabot, R. & Antoun, H. Growth promotion of maize and lettuce by phosphate solubilizing Rhizobium leguminosarum. Plant Soil. 184, 311–321 (1996).CAS 
Article 

Google Scholar 
Kundu, B. S. & Gaur, A. C. Rice response to inoculation with N2 fixing and P solubilizing microorganisms. Plant Soil. 79, 227–234 (1984).CAS 
Article 

Google Scholar 
Sharma, G. D., Thakur, R., Raj, S., Kauraw, D. L. & Kulhare, P. S. Impact of integrated nutrient management on yield, nutrient uptake, protein content of wheat (Triticum aestivum) and soil fertility in a typic Haplustert. Bioscan 8, 1159–1164 (2013).CAS 

Google Scholar 
Afzal, A. & Asghari, B. Rhizobium and phosphate solubilizing bacteria improve the yield and phosphorus uptake in wheat (Triticum aestivum). Int. J. Agric. Biol. 10, 85–88 (2008).CAS 

Google Scholar 
Jalili, G. et al. Enhancing crop growth, nutrients availability, economics and beneficial rhizosphere micro flora through organic and bio fertilizers. Ann. Microbiol. 57(2), 177–183 (2007).Article 

Google Scholar 
Sharma, S. N. & Prasad, R. Yield and P uptake by rice and wheat grown in a sequence as influenced by phosphate fertilization with diammonium phosphate and Mussoorie rock phosphate with or without crop residues and phosphate solubilizing bacteria. J. Agric. Sci. 141, 359–369 (2003).CAS 
Article 

Google Scholar 
Vyas, P. & Gulati, A. Organic acid production in vitro and plant growth promotion in maize under controlled environment by phosphate-solubilizing fluorescent Pseudomonas. BMC Microbiol. 9, 174 (2009).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Mukherjee, P. K. & Rai, R. K. Sensitivity of P uptake to change in root growth and soil volume as influenced by VAM, PSB and P levels in wheat and chickpeas. Ann. Agric. Res. 20, 528–530 (1999).
Google Scholar 
Egamberdiyeva, D. Proc. Inst. Microbiol. Tashkent, Uzekistan (2004).Mihoub, A., Daddi Bouhoun, M., Naeem, A. & Saker, M. L. Low-molecular weight organic acids improve plant availability of phosphorus in different textured calcareous soils. Arch. Agron. Soil Sci. 63, 1023–1034 (2017).CAS 
Article 

Google Scholar 
Thakuria, D. et al. Characterization and screening of bacteria from rhizosphere of rice grown in acidic soils of Assam. Curr. Sci. 86, 978–985 (2004).
Google Scholar 
Mamta, P. et al. Stimulatory effect of phosphate solubilizing bacteria on plant growth, stevioside and rebaudioside-A content of Stevia rebaudiana Bertoni. Appl. Soil Ecol. 46, 222–229 (2010).Article 

Google Scholar 
Banik, S. B. K. Solubilization of inorganic phosphate and production of organic acids by micro-organisms isolated in sucrose tricalcium phosphate agar plate. Zentralblat. Bakterol. Parasilenkl. Infektionskr. Hyg. 136, 478–486 (1981).CAS 

Google Scholar 
Stevenson, F. J. Cycles of Soil: Carbon, Nitrogen, Phosphorus, Sulfur, Micro-nutrients (Wiley, 2005).
Google Scholar 
Ekin, Z. Performance of phosphorus solubilizing bacteria for improving growth and yield of sun flower (Helianthus annuus L.) in the presence of phosphorus fertilizer. Afr. J. Biotechnol. 9, 3794–3800 (2010).CAS 

Google Scholar 
Zabihi, H. R., Savaghebi, G. R., Khavazi, K., Ganjali, A. & Miransari, M. Pseudomonas bacteria and phosphorus fertilization, affecting wheat (Triticum aestivum L.) yield and P uptake under green house and field conditions. Acta Physiol. Plant 33, 145–152 (2010).Article 

Google Scholar 
Gulati, A., Rahi, P. & Vyas, P. Characterization of phosphate-solubilizing fluorescent Pseudomonas from the rhizosphere of seabuckthorn growing in the cold deserts of Himalayas. Curr. Microbiol. 56, 73–79 (2007).PubMed 
Article 
CAS 

Google Scholar 
Kloepper, J. W., Lifshitz, R. & Zablotowicz, R. M. Free-living bacterial inocula for enhancing crop productivity. Trends Biotechnol. 7, 39–44 (1989).Article 

Google Scholar 
Satchell, J. E. Ecology and environment in the United Arab Emirates. J. Arid. Environ. 1, 201–226 (1978).ADS 
Article 

Google Scholar 
Biswas, D. R. Nutrient recycling potential of rock phosphate and waste mica enriched compost on crop productivity and changes in soil fertility under potato–soybean cropping sequence in an Inceptisol of Indo-Gangetic Plains of India. Nutr. Cycl. Agroecosyst. 89, 15–30 (2011).Article 

Google Scholar 
Mitra, S. et al. Effect of integrated nutrient management on fiber yield, nutrient uptake and soil fertility in jute (Corchorus olitorius). Indian J. Anim. Sci. 80(9), 801–804 (2010).
Google Scholar 
Laxminarayana, K. Effect of integrated use of inorganic and organic manures on soil properties, yield and nutrient uptake of rice in Ultisols of Mizoram. J. Indian Soc. Soil Sci. 54, 120–123 (2006).
Google Scholar 
Sanyal, S. K. & De Datta, S. K. Chemistry of phosphorus transformations in soil. Adv. Soil Sci. 16, 1–120 (1991).CAS 

Google Scholar 
Briedis, C. et al. Soil organic matter pools and carbon-protection mechanisms in aggregate classes influenced by surface liming in a no-till system. Geoderma 170, 80–88 (2012).ADS 
CAS 
Article 

Google Scholar 
Bronick, C. J. & Lal, R. Soil structure and management: A review. Geoderma 124, 3–22 (2005).ADS 
CAS 
Article 

Google Scholar 
Krieg, N. R. & Holt, J. G. Bergey’s Manual of Systemetic Bacteriology Vol. 1, 984 (Williams & Wilkin, 1984).
Google Scholar 
Holt, J. G. et al. (eds) Bergey’s Manual of Determinative Bacteriology 9th edn, 787 (The Williams & Wilkin, 1994).
Google Scholar 
Gordon, R. E., Haynes, W. C. & Pang, C. N. The Genus Bacillus. Agricultural Handbook. No. 427 283 (Department of Agriculture, 1973).
Google Scholar 
Nautiyal, C. S. An efficient microbiological growth medium for screening phosphate solubilizing microorganisms. FEMS Microbiol. Lett. 170(1), 265–270 (1999).CAS 
PubMed 
Article 

Google Scholar 
Nelson, D. W. & Sommers, L. E. Total carbon, organic carbon, and organic matter. In Methods of Soil Analysis, Part 2 2nd edn, Vol. 14 (ed. Page, A. L.) 961–1010 (Wiley, 1996).
Google Scholar 
Eivazi, F. & Tabatabai, M. Phosphatases in soils. Soil Biol. Biochem. 9, 167–172 (1977).CAS 
Article 

Google Scholar 
Alexander, D. B. & Zuberer, D. A. Use of chrome azurol S reagents to evaluate siderophore production by rhizosphere bacteria. Biol. Fertil. Soils 12, 39–45 (1991).CAS 
Article 

Google Scholar 
Vincet, J. M. A. Manual for the Practical Study of the Root-Nodule Bacteria; IBPH and Book No. 15 (Blackwell Scientific Publication, 1970).
Google Scholar 
Alagawadi, A. R. & Gaur, A. C. Associative effect of Rhizobium and phosphate solubilizing bacteria on the yield and nutrient uptake of chickpea. Plant Soil. 105, 241–246 (1988).Article 

Google Scholar 
Satyaprakash, M., Nikitha, T., Reddi, E. U. B., Sadhana, B. & Vani, S. S. Phosphorous and phosphate solubilising bacteria and their role in plant nutrition. Int. J. Curr. Microbiol. Appl. Sci. 6, 2133–2144 (2017).CAS 
Article 

Google Scholar 
Wu, S. C., Cao, Z. H., Li, Z. G., Cheung, K. C. & Wong, M. H. Effects of biofertilizer containing N-fixer, P and K solubilizers and AM fungi on maize growth: A greenhouse trial. Geoderma 125, 155–166 (2005).ADS 
Article 

Google Scholar 
Thomas, G. W. Soil pH and soil acidity. In Methods of Soil Analysis, Part 3, Chemical Methods Vol. 5 (eds Sparks, D. L. et al.) 475–490 (Wiley, 1996).
Google Scholar 
Rhoades, J. D. Salinity, electrical conductivity and total dissolved solids. In Methods of Soil Analysis, Part 3, Chemical Methods Vol. 5 (eds Sparks, D. L. et al.) 417–435 (Soil Science Society of America, 1996).
Google Scholar 
Bremner, J. M. & Breitenbeck, G. A. A simple method for determination of ammonium in semi-micro Kjeldahl analysis of soil and plant material using a block digestor. Commun. Soil Sci. Plant Anal. 14, 905–913 (1983).CAS 
Article 

Google Scholar 
Ryan, J., Estefan, G. & Rashid, A. Soil and Plant Analysis Laboratory Manual 2nd edn, 172 (The National Agricultural Research Center (NARC), 2001).
Google Scholar 
Olsen, S. R., Cole, C. V., Watanabe, F. S. & Dean, L. A. Estimation of Available Phosphorus in Soils by Extraction with Sodium Bicarbonate (No. 939) (Department of Agriculture Circular, 1954).
Google Scholar 
Loeppert, R. H. & Suarez, D. L. Carbonate and gypsum. In Methods of Soil Analysis, Part 3, Chemical Methods Vol. 9 (eds Sparks, D. L. et al.) 181–197 (Soil Science Society of America, 1996).
Google Scholar 
Bahadur, L., Tiwari, D. D., Mishra, J. & Gupta, B. R. Effect of integrated nutrient management on yield, microbial population and changes in soil properties under rice-wheat cropping system in sodic soil. J. Indian Soc. Soil Sci. 60(4), 326–329 (2012).CAS 

Google Scholar 
Nelson, D. W. & Sommers, L. E. Total carbon, organic carbon, and organic matter. In Methods of Soil Analysis, Part 2 2nd edn, Vol. 9 (eds Sparks, D. L. et al.) 961–1010 (Soil Science Society of America, 1996).
Google Scholar 
Richards, L. A. Diagnosis and improvement of saline and alkali soils. LWW 78(2), 154 (1954).
Google Scholar 
Steel, R. G. D. & Torrie, J. H. Principles and Procedures of Statistics, a Biometrical Approach 195–233 (McGraw Hill, 1996).MATH 

Google Scholar  More

Site descriptionThe experiment was conducted at the Beef Cattle Research Center of the Institute of Animal Science/APTA/SAA, Sertãozinho, São Paulo, Brazil (21°08′16″ S e 47°59′25″ W, average altitude 548 m), during two consecutive years. The climate in this region is Aw according to the Köppen’s classification, characterized as humid tropical, with a rainy season during summer and drought during winter. The meteorological data is reported in Fig. 1. The soil in the experimental area is classified as an Oxisol42. Before the experiment, soil samples were collected for chemical characterization (Table 4), which was performed following the methodology described in Van Raij et al.43. Samples were collected in 18 experimental paddocks, at the depths of 0- to 10- and 10- to 20-cm layers, from 10 distinct sampling points in each paddock, in order to create one composite sample per unit, totaling 36 samples analyzed.Figure 1Meteorological data during the study period, obtained from the meteorological station located at Centro de Pesquisa de Bovinos de Corte, Instituto de Zootecnia/Agência Paulista de Tecnologia dos Agronegócios (APTA)/Secretaria de Agricultura e Abastecimento de São Paulo (SAA), Sertãozinho, São Paulo, Brazil.Full size imageTable 4 Chemical attributes of the soil in the experimental area, before installing the experiment (November 2015).Full size tableThe nitrogen total (Nt) content was determined by the micro-Kjeldahl method44, and the soil nitrogen stocks (SN) were calculated using the following equation below, according to Veldkamp et al.45.$${text{SN }}left[ {{text{Mg ha}}^{ – 1} {text{ at a given depth}}} right], = ,({text{concentration }} times {text{ BD}}, times ,{1}/{1}0),$$ where concentration refers to the Nt concentration at a given depth (g kg−1), BD is the bulk density at a certain depth (average 1.24 kg dm−3), and 1 is the layer thickness (cm).Description of treatments and managementsThe experiment was carried out in a 16-ha area, divided into 18 paddocks of 0.89 ha each (Fig. 2), organized in a randomized blocks design with three replicates and six treatments, namely conventional crop system with grain maize production (CROP), conventional livestock system with beef cattle production in pasture using Marandu grass (LS), and four ICLS for the production of intercropped maize grain with beef cattle pasture. All production systems were sowed in December 2015, under a no-tillage system. The fertilization recommendations in the systems were based on the recommendation presented in the Boletim 10046.Figure 2Localization and representation of the area of the experiment carried out in the study. Google Earth version Pro was used to construct the map (http://www.google.com/earth/index.html).Full size imageIn the CROP system, the maize Pioneer P2830H was cultivated, sowed in a spacing of 75 cm and sowing density of 70 thousand plants. Applications of 32 kg ha−1 of nitrogen (urea), 112 kg ha−1 of P2O5 (single superphosphate) and 64 kg ha−1 of KCl (potassium chloride) were performed. Complementarily, a topdressing fertilization was made using 80 kg ha−1 of nitrogen (urea) and 80 kg ha−1 of KCl. Sowing was carried out for two consecutive years (December 2015 and 2016), providing two harvests of maize grains (May 2016 and 2017), and between one harvest and the other, the soil remained in fallow without any cover crop. The total amount of fertilizer applied in two years was 224 kg ha−1 of nitrogen (urea), 224 kg ha−1 of P2O5 (single superphosphate) and 288 kg ha−1 of KCl (potassium chloride).For the LS treatment, Urochloa brizantha (Hoechst. ex A. Rich) R.D. Webster cv. Marandu (syn. Brachiaria brizantha cv. Marandu) was sowed in a spacing of 37.5 cm, with a density of 5 kg ha−1 of seeds (76% of crop value) for the pasture assemblage. Marandu grass seeds were mixed with the planting fertilizer, applying 32 kg ha−1 of nitrogen (urea), 112 kg ha−1 of P2O5 (as single superphosphate) and 64 kg ha−1 of KCl. Applications of 40 kg ha−1 of nitrogen, 10 kg ha−1 of P2O5 and 40 kg ha−1 of KCl were also performed as topdressing fertilization in October 2016 and March 2017. 90 days after sowing, the pasture was ready to be grazed (March 2016). Three grazing periods were carried out in continuous stocking systems, with the first period between March and April 2016, the second period between August and October 2016 and the third between November 2016 and December 2017. The total amount for 2 years was 112 kg ha−1 of nitrogen (urea), 132 kg ha−1 of P2O5 (single superphosphate) and 144 kg ha−1 of KCl (potassium chloride).The same cultivar, spacing, sowing density and fertilization rates described in the CROP treatment were used in all ICLS, as well as the same density of Marandu grass seeds and topdressing fertilization adopted in the pasture of the LS treatment. The total amount for two years was 192 kg ha−1 of nitrogen (urea), 132 kg ha−1 of P2O5 (single superphosphate) and 224 kg ha−1 of KCl (potassium chloride). In ICLS-1, Marandu grass was sowed in lines simultaneously with maize, while in ICLS-2, the sowing was also simultaneous, but the application of an under-dose of 200 mL of the herbicide Nicosulfuron was used, 20 days after seedlings emergence. In the ICLS-3, Marandu grass seeds were sown the time of topdressing fertilization of maize, thus the grass seeds were mixed with the fertilizer, and sowing was carried out in the interlines of maize, using a minimum cultivator. In ICLS-4, the sowing of Marandu grass was performed simultaneously with maize, but the grass seeds were sowed in both rows and inter-rows of maize, resulting in a spacing of 37.5 cm. In this treatment, the application of 200 mL of the herbicide Nicosulfuron was adopted, 20 days after seedlings emergence.In all ICLS treatments, maize harvest was carried out in May 2016. Ninety days after harvesting the plants, the pastures were ready to be grazed. Therefore, two grazing periods were made in continuous stocking, being the first period between August and October 2016 and the second period between November 2016 and December 2017. The method for animal stocking in treatments LS and ICLS was continuous with a stocking rate (put and take) being defined according to Mott47. Caracu beef cattle with 14 months of age were used at the beginning of the experiment, with an average body weight of 335 ± 30 kg.Estimations of the nutrient balance (NB) and nutrient use efficiency (NUE)In this study, the inputs and outputs of N were assessed at the farm level48,49. The NB was calculated by the equation below19,45,50.$${text{NB}}_{{text{N}}} = {text{ Input}}_{{text{N}}} {-}{text{ Output}}_{{text{N}}}$$As for the NUE, this parameter was evaluated as defined by the EU Nitrogen Expert Panel51, being calculated as the ratio between outputs and inputs of nitrogen.$${text{NUE}}_{{text{N}}} = , left[ {{text{Output}}_{{text{N}}} /{text{ Input}}_{{text{N}}} } right]$$where NB is the nutrient balance, N is nitrogen, Input is the N concentration in the mineral fertilizer (urea), Output is the nitrogen concentration in export (maize grain and animal tissue), and NUE is the use efficiency of the nutrient.The amount of N exported in maize grains, the grain production results (Table 2) were multiplied by the mean value of N, consulted in Crampton and Harris52.In order to estimate the amounts of nutrient exported by the animals in their tissues, the values of live weight gain were considered [kg ha-1 of live weight (PV)] (Table 2), as well as the nitrogen values of the tissue, according to the methodology proposed by Rasmussen et al.21. Those authors reported that for animals weighting less than 452 kg/PV, it represents 2.7%, while heavier animals have a 2.4% nitrogen content representation of their body weight.The inputs and outputs of N in each production system are represented in Figs. 3, 4 and 5. Biological N fixation, atmospheric deposition, denitrification, leaching, rainfall, and volatilization and absorption of ammonia were not considered in the calculation of NB.Figure 3Representation of inputs and outputs of nitrogen and organic residues generated in the crop system.Full size imageFigure 4Representation of inputs and outputs of nitrogen and organic residues generated in the livestock system.Full size imageFigure 5Representation of inputs and outputs of nitrogen and organic residues generated in the integrated systems.Full size imageData for animal tissue, animal excreta, and N concentration in grains were obtained from key manuscripts from the scientific literature in order to estimate the N balance.Calculation of nitrogen quantity and valuation of organic residuesThe amount of N in the organic residues was determined as a function of the system (Figs. 3, 4, 5). The residue considered in the CROP was the straw derived from maize, while for LS it was the litter deposited (LD) in the grass Marandu, and animal manure (feces and urine). The ICLS were considered as the straw, LD, and animal manure.The N concentration in straw and LD was determined following the methods of AOAC (1990). Straw was sampled immediately after maize grain harvest, using a 1-m2 frame in the field. The material was collected in two spots of the plot that were chosen randomly. All straw deposited on the soil was sampled, weighted and dried in an oven with air circulation (60 °C) until constant weight, for the determination of dry matter in kg of straw per hectare (Table 2). The LD in the pasture system (Table 2) was analyzed according to Rezende et al.53.In order to estimate the daily amount of excreta, we considered the stocking rate adopted in the experiment (Table 2) and the values proposed by Haynes and Williams54. According to those authors, adult beef cattle can defecate on average 13 times a day and urinate 10 times a day, totaling a daily amount of 28.35 kg of feces and 19 L of urine.The valuation was calculated based on the mean value of urea for the last 10 years in the fertilizer market55,56,57, namely $0.28 kg−1 ha−1 of urea, and considering the loss of nitrogen by volatilization, which according to Freney et al.58 and Subair et al.59 can reach up to 28%.Statistical analysisThe experiment was assembled in a randomized blocks design. The model adopted for the analysis of all response variables included the block’s and treatments fixed effects (3 blocks and 6 treatments), in addition to the random error. Statistical analysis were carried out by the function “dbc()” of the package “ExpDes.pt” of the software R Development Core Team60, and the mean values were compared by the Tukey’s test at a 5% probability level. More

Rundel, P. W., Graham, E. A., Allen, M. F., Fisher, J. C. & Harmon, T. C. New Phytol. 182, 589–607 (2009).Article 

Google Scholar 
Gibb, R., Browning, E., Glover‐Kapfer, P. & Jones, K. E. Methods Ecol. Evol. 10, 169–185 (2019).Article 

Google Scholar 
O’Connell, A. F. (ed) Camera Traps in Animal Ecology: Methods and Analyses. Vol. 271 (Springer, 2011).Hale, R. C., Seeley, M. E., Guardia, M. J. L., Mai, L. & Zeng, E. Y. J. Geophys. Res. Oceans 125, e2018JC014719 (2020).Article 

Google Scholar 
Widmer, R., Oswald-Krapf, H., Sinha-Khetriwal, D., Schnellmann, M. & Böni, H. Environ. Impact Assess. Rev. 25, 436–458 (2005).Article 

Google Scholar 
Hwang, S.-W. et al. Science 337, 1640–1644 (2012).CAS 
Article 

Google Scholar 
Ashammakhi, N. et al. Adv. Funct. Mater. 31, 2104149 (2021).Boutry, C. M. et al. Nat. Biomed. Eng. 3, 47–57 (2019).CAS 
Article 

Google Scholar 
Boutry, C. M. et al. Nat. Electron. 1, 314–321 (2018).Article 

Google Scholar 
Hori, K., Inami, A., Kan, T. & Onoe, H. In Proc. 21st International Conference on Solid-State Sensors, Actuators and Microsystems (Transducers) 863–866 (IEEE, Orlando, 2021).Dincer, C. et al. Adv. Mater. 31, 1806739 (2019).Article 

Google Scholar 
Kocer, B. B. et al. In Proc. Aerial Robotic Systems Physically Interacting with the Environment (AIRPHARO) 1–8 (IEEE, Biograd na Moru, 2021).Pandolfi, C. & Izzo, D. Bioinspir. Biomim. 8, 025003 (2013).Article 

Google Scholar 
Wiesemüller, F., Miriyev, A. & Kovac, M. In Proc. Aerial Robotic Systems Physically Interacting with the Environment (AIRPHARO) 1–6 (IEEE, Biograd na Moru, 2021).Boutry, C. M. et al. Sens. Actuators A Phys. 189, 344–355 (2013).CAS 
Article 

Google Scholar 
Tsang, M., Armutlulu, A., Martinez, A. W., Allen, S. A. B. & Allen, M. G. Microsyst. Nanoeng. 1, 15024 (2015).CAS 
Article 

Google Scholar 
Lee, G. et al. Adv. Energy Mater. 7, 1700157 (2017).Article 

Google Scholar 
Dagdeviren, C. et al. Small 9, 3398–3404 (2013).CAS 
Article 

Google Scholar 
Sadasivuni, K. K. et al. J. Mater. Sci. Mater. Electron. 30, 951–974 (2019).CAS 
Article 

Google Scholar 
Luvisi, A., Panattoni, A. & Materazzi, A. Comput. Electron. Agric. 123, 135–141 (2016).Article 

Google Scholar 
Yin, L. et al. Adv. Mater. 26, 3879–3884 (2014).CAS 
Article 

Google Scholar 
Demetillo, A. T., Japitana, M. V. & Taboada, E. B. Sustain. Environ. Res. 29, 12 (2019).CAS 
Article 

Google Scholar 
Salvatore, G. A. et al. Adv. Funct. Mater. 27, 1702390 (2017).Article 

Google Scholar 
Farinha, A., Zufferey, R., Zheng, P., Armanini, S. F. & Kovac, M. IEEE Robot. Autom. Lett. 5, 6623–6630 (2020).Article 

Google Scholar 
Miriyev, A. & Kovač, M. Nat. Mach. Intell. 2, 658–660 (2020).Article 

Google Scholar 
Kang, S.-K., Koo, J., Lee, Y. K. & Rogers, J. A. Acc. Chem. Res. 51, 988–998 (2018).CAS 
Article 

Google Scholar 
Goel, V., Luthra, P., Kapur, G. S. & Ramakumar, S. S. V. J. Polym. Environ. 29, 3079–3104 (2021).CAS 
Article 

Google Scholar  More

Brondizio, E. S., Settele, J., Díaz, S. & Ngo, H. T. Global Assessment Report on Biodiversity and Ecosystem Services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. (IPBES, 2019).Tingley, M. W., Monahan, W. B., Beissinger, S. R. & Moritz, C. Proc. Natl Acad. Sci. USA 106(Suppl 2), 19637–19643 (2009).CAS 
Article 

Google Scholar 
Schloss, C. A., Nuñez, T. A. & Lawler, J. J. Proc. Natl Acad. Sci. USA 109, 8606–8611 (2012).CAS 
Article 

Google Scholar 
Senior, R. A., Hill, J. K. & Edwards, D. P. Nat. Clim. Chang. 9, 623–626 (2019).Article 

Google Scholar 
Viana, D. S. & Chase, J. M. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-022-01814-y (2022).Article 

Google Scholar 
Sauer, J. R. et al. Condor 119, 576–593 (2017).Article 

Google Scholar 
Nowak, L., Schleuning, M., Bender, I. M. A., Kissling, W. D. & Fritz, S. A. Divers. Distrib. https://doi.org/10.1111/ddi.13518 (2022).Article 

Google Scholar 
Allen, C. D. et al. For. Ecol. Manage. 259, 660–684 (2010).Article 

Google Scholar 
Janis, C. M., Damuth, J. & Theodor, J. M. Proc. Natl Acad. Sci. USA 97, 7899–7904 (2000).CAS 
Article 

Google Scholar 
Stuart-Smith, R. D., Mellin, C., Bates, A. E. & Edgar, G. J. Nat. Ecol. Evol. 5, 656–662 (2021).Article 

Google Scholar 
Watanabe, Y. Y. Ecol. Lett. 19, 907–914 (2016).Article 

Google Scholar 
Bladon, A. J. et al. J. Anim. Ecol. 89, 2440–2450 (2020).Article 

Google Scholar 
Claramunt, S., Hong, M. & Bravo, A. Biotropica https://doi.org/10.1111/btp.13109 (2022).Article 

Google Scholar 
Zurell, D., Gallien, L., Graham, C. H. & Zimmermann, N. E. J. Biogeogr. 45, 1459–1468 (2018).Article 

Google Scholar 
Bowler, D. E., Heldbjerg, H., Fox, A. D., O’Hara, R. B. & Böhning-Gaese, K. J. Anim. Ecol. 87, 1034–1045 (2018).Article 

Google Scholar 
Warren, D. L., Cardillo, M., Rosauer, D. F. & Bolnick, D. I. Trends Ecol. Evol. 29, 572–580 (2014).Article 

Google Scholar 
Gómez, C., Tenorio, E. A., Montoya, P. & Cadena, C. D. Proc. R. Soc. Lond. B. Biol. Sci. 283, 20152458 (2016).
Google Scholar 
Amano, T., Lamming, J. D. L. & Sutherland, W. J. Bioscience 66, 393–400 (2016).Article 

Google Scholar 
Rosenberg, K. V. et al. Science 366, 120–124 (2019).CAS 
Article 

Google Scholar 
Howard, C. et al. Divers. Distrib. 26, 1442–1455 (2020).Article 

Google Scholar  More

Banks, N. C., Paini, D. R., Bayliss, K. L. & Hodda, M. The role of global trade and transport network topology in the human-mediated dispersal of alien species. Ecol. Lett. 18, 188–199 (2015).
Google Scholar 
Epanchin-Niell, R. et al. Controlling invasive species in complex social landscapes. Front. Ecol. Environ. 8, 210–216 (2009).
Google Scholar 
Charles, H. & Dukes, J. S. in Biological Invasions (ed. Nentwig, W.) 217–237 (Springer, 2007). https://doi.org/10.1007/978-3-540-36920-2_13Gallardo, B., Clavero, M., Sánchez, M. & Vilà, M. Global ecological impacts of invasive species in aquatic ecosystems. Glob. Change Biol. 22, 151–163 (2016).
Google Scholar 
Diagne, C. et al. High and rising economic costs of biological invasions worldwide. Nature 592, 571–576 (2021).CAS 

Google Scholar 
Sardain, A., Sardain, E. & Leung, B. Global forecasts of shipping traffic and biological invasions to 2050. Nat. Sustain. 2, 274–282 (2019).
Google Scholar 
Epanchin-Niell, R. S. & Hastings, A. Controlling established invaders: integrating economics and spread dynamics to determine optimal management. Ecol. Lett. 13, 528–541 (2010).
Google Scholar 
Chades, I. et al. General rules for managing and surveying networks of pests, diseases, and endangered species. Proc. Natl. Acad. Sci. USA 108, 8323–8328 (2011).CAS 

Google Scholar 
Epanchin-Niell, R. S. & Wilen, J. E. Optimal spatial control of biological invasions. J. Environ. Econ. Manag. 63, 260–270 (2012).
Google Scholar 
Epanchin-Niell, R. S. & Wilen, J. E. Individual and cooperative management of invasive species in human-mediated landscapes. Am. J. Agric. Econ. 97, 180–198 (2015).
Google Scholar 
Aadland, D., Sims, C. & Finnoff, D. Spatial dynamics of optimal management in bioeconomic systems. Comput. Econ. 45, 545–577 (2015).
Google Scholar 
Baker, C. M. Target the source: optimal spatiotemporal resource allocation for invasive species control. Conserv. Lett. 10, 41–48 (2017).
Google Scholar 
Bushaj, S., Büyüktahtakın, İ. E., Yemshanov, D. & Haight, R. G. Optimizing surveillance and management of emerald ash borer in urban environments. Nat. Res. Model. 34, e12267 (2021).
Google Scholar 
Fischer, S. M., Beck, M., Herborg, L.-M. & Lewis, M. A. Managing aquatic invasions: optimal locations and operating times for watercraft inspection stations. J. Environ. Manag. 283, 111923 (2021).
Google Scholar 
Büyüktahtakın, İ. E. & Haight, R. G. A review of operations research models in invasive species management: state of the art, challenges, and future directions. Ann. Oper. Res. 271, 357–403 (2018).
Google Scholar 
Epanchin-Niell, R. S. Economics of invasive species policy and management. Biol. Invasions 19, 3333–3354 (2017).
Google Scholar 
Bodin, Ö. et al. Improving network approaches to the study of complex social–ecological interdependencies. Nat. Sustain. 2, 551–559 (2019).CAS 

Google Scholar 
Nowzari, C., Precaido, V. M. & Pappas, G. J. Analysis and control of epidemics: a survey of spreading processes on complex networks. IEEE Control Syst. 36, 26–46 (2016).
Google Scholar 
Newman, M. E. J. Spread of epidemic disease on networks. Phys. Rev. E 66, 016128 (2002).CAS 

Google Scholar 
Kempe, D., Kleinberg, J. & Tardos, E. Maximizing the spread of influence through a social network. In Proc. 9th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining 137–146 (ACM Press, 2003).Pastor-Satorras, R. & Vespignani, A. Immunization of complex networks. Phys. Rev. E 65, 036104 (2002).
Google Scholar 
Pastor-Satorras, R., Castellano, C., Van Mieghem, P. & Vespignani, A. Epidemic processes in complex networks. Rev. Mod. Phys. 87, 925–979 (2015).
Google Scholar 
Holme, P., Kim, B. J., Yoon, C. N. & Han, S. K. Attack vulnerability of complex networks. Phys. Rev. E 65, 056109 (2002).
Google Scholar 
Muirhead, J. R. & Macisaac, H. J. Development of inland lakes as hubs in an invasion network. J. Appl. Ecol. 42, 80–90 (2005).
Google Scholar 
de la Fuente, B., Saura, S. & Beck, P. S. Predicting the spread of an invasive tree pest: the pine wood nematode in southern europe. J. Appl. Ecol. 55, 2374–2385 (2018).
Google Scholar 
Minor, E. S. & Urban, D. L. A graph-theory framework for evaluating landscape connectivity and conservation planning. Conserv. Biol. 22, 297–307 (2008).
Google Scholar 
Morel-Journel, T., Assa, C. R., Mailleret, L. & Vercken, E. Its all about connections: hubs and invasion in habitat networks. Ecol. Lett. 22, 313–321 (2019).
Google Scholar 
Perry, G. L. W., Moloney, K. A. & Etherington, T. R. Using network connectivity to prioritise sites for the control of invasive species. J. Appl. Ecol. 54, 1238–1250 (2017).
Google Scholar 
Kvistad, J. T., Chadderton, W. L. & Bossenbroek, J. M. Network centrality as a potential method for prioritizing ports for aquatic invasive species surveillance and response in the Laurentian Great Lakes. Manag. Biol. Invasions 10, 403 (2019).
Google Scholar 
Haight, R. G., Kinsley, A. C., Kao, S.-Y., Yemshanov, D. & Phelps, N. B. Optimizing the location of watercraft inspection stations to slow the spread of aquatic invasive species. Biol. Invasions 23, 3907–3919 (2021).
Google Scholar 
McEachran, M. C. et al. Stable isotopes indicate that zebra mussels (Dreissena polymorpha) increase dependence of lake food webs on littoral energy sources. Freshw, Biol. 64, 183–196 (2019).CAS 

Google Scholar 
Karatayev, A. Y., Burlakova, L. E. & Padilla, D. K. in Invasive Aquatic Species of Europe. Distribution, Impacts and Management (eds Leppäkoski, E. et al.) 433–446 (Springer, 2002).Prescott, T. H., Claudi, R. & Prescott, K. L. Impact of Dreissenid mussels on the infrastructure of dams and hydroelectric power plants. In Quagga and Zebra Mussels (eds Nalepa, T. F. & Schloesser, D. W.) 243–258 (CRC Press, 2013).Invasive Species of Aquatic Plants and Wild Animals in Minnesota: Annual Report for 2020 (Minnesota Department of Natural Resources, 2020).Kanankege, K. S., Alkhamis, M. A., Phelps, N. B. & Perez, A. M. A probability co-kriging model to account for reporting bias and recognize areas at high risk for zebra mussels and eurasian watermilfoil invasions in Minnesota. Front. Vet. Sci. 4, 231 (2018).
Google Scholar 
Mallez, S. & McCartney, M. Dispersal mechanisms for zebra mussels: population genetics supports clustered invasions over spread from hub lakes in Minnesota. Biol. Invasions 20, 2461–2484 (2018).
Google Scholar 
Kao, S.-Y. Z. et al. Network connectivity of Minnesota waterbodies and implications for aquatic invasive species prevention. Biol. Invasions 23, 3231–3242 (2021).
Google Scholar 
Kleinberg, J. M. Authoritative sources in a hyperlinked environment. In Proc. 9th Annual ACM-SIAM Symposium on Discrete Algorithms 668–677 (1998).McDonald-Madden, E. et al. Using food-web theory to conserve ecosystems. Nat. Commun. 7, 10245 (2016).CAS 

Google Scholar 
Bossenbroek, J. M., Kraft, C. E. & Nekola, J. C. Prediction of long-distance dispersal using gravity models: zebra mussel invasion of inland lakes. Ecol. Appl. 11, 1778–1788 (2001).
Google Scholar 
Leung, B., Bossenbroek, J. M. & Lodge, D. M. Boats, pathways, and aquatic biological invasions: estimating dispersal potential with gravity models. Biol. Invasions 8, 241–254 (2006).
Google Scholar 
Beger, M. et al. Integrating regional conservation priorities for multiple objectives into national policy. Nat. Commun. 6, 8208 (2015).Runting, R. K. et al. Larger gains from improved management over sparing–sharing for tropical forests. Nat. Sustain. 2, 53–61 (2019).
Google Scholar 
Kinsley, A. C. et al. AIS Explorer: prioritization for watercraft inspections. A decision-support tool for aquatic invasive species management. J. Environ. Manage. 314, 115037 (2022).
Google Scholar 
Vander Zanden, M. J. & Olden, J. D. A management framework for preventing the secondary spread of aquatic invasive species. Can. J. Fish. Aquat. Sci. 65, 1512–1522 (2008).
Google Scholar 
Kanankege, K. S. et al. Lessons learned from the stakeholder engagement in research: application of spatial analytical tools in one health problems. Front. Vet. Sci. 7, 254 (2020).
Google Scholar 
Kroetz, K. & Sanchirico, J. The bioeconomics of spatial-dynamic systems in natural resource management. Annu. Rev. Resour. Econ. 7, 189–207 (2015).
Google Scholar 
Cade, B. S. & Noon, B. R. A gentle introduction to quantile regression for ecologists. Front. Ecol. Environ. 1, 412–420 (2003).
Google Scholar 
Koenker, R. in Asymptotic Statistics (eds Mandl, P. & Hušková, M.) 349–359 (Springer, 1994).Ashander, J. Analysis code and data for ‘Guiding large-scale management of invasive species using network metrics’. figshare https://doi.org/10.6084/m9.figshare.14402447 (2021). More

Climate and land cover dataOur study of tropical dry season dynamics required climatic variables with high temporal resolution (i.e., daily) and full coverage of tropic regions. To reduce uncertainties associated with the choice of precipitation (P) and evapotranspiration (Ep or E) datasets, we used an ensemble of eight precipitation products, three reanalysis-based products for Ep, and one satellite-based land E product. These precipitation datasets were derived four gauge-based or satellite observation (CHIRPS58, GPCC59, CPC-U60 and PERSIANN-CDR61), three reanalyses (ERA-562, MERRA-263, and PGF64) and a multi-source weighted ensemble product (MSWEP v2.865). The potential evapotranspiration (Ep) was calculated using the FAO Penman–Monteith equation66 (Eqs. (1, 2)), which requires meteorological inputs of wind speed, net radiation, air temperature, specific humidity, and surface pressure. We derived these meteorological variables from the three reanalysis products (ERA-5, MERRA-2, and GLDAS-2.067). Since PGF reanalysis lacked upward short- and long-wave radiation output and thus net radiation, we used available meteorological outputs from GLDAS-2.0 instead, which was forced entirely with the PGF input data.$${Ep}=frac{0.408cdot triangle cdot left({R}_{n}-Gright)+gamma cdot frac{900}{T+273}cdot {u}_{2}cdot left({e}_{s}-{e}_{a}right)}{triangle +{{{{{rm{gamma }}}}}}cdot left(1+0.34cdot {u}_{2}right)}$$
(1)
$${VPD}={e}_{s}-{e}_{a}=0.6108cdot {e}^{frac{17.27cdot T}{T+237.3}}cdot left(1-frac{{RH}}{100}right)$$
(2)
Where Ep is the potential evapotranspiration (mm day−1). Rn is net radiation at the surface (MJ m−2 day−1), T is mean daily air temperature at 2 m height (°C), ({u}_{2}) is wind speed at 2 m height (m s−1), ((,{e}_{s}-{e}_{a})) is the vapor pressure deficit of the air (kPa), ({RH}) is the relative air humidity near surface (%), ∆ is the slope of the saturation vapor pressure-temperature relationship (kPa °C−1), γ is the psychrometric constant (kPa °C−1), G is the soil heat flux (MJ m−2 day−1, is often ignored for daily time steps G ≈ 0).We derived the daily evapotranspiration data from the Global Land Evaporation Amsterdam Model (GLEAM v3.3a68), which is a set of algorithms dedicated to developing terrestrial evaporation and root-zone soil moisture data. GLEAM fully assimilated the satellite-based soil moisture estimates from ESA CCI, microwave L-band vegetation optical depth (VOD), reanalysis-based temperature and radiation, and multi-source precipitation forcings. The direct assimilation of observed soil moisture allowed us to detect true soil moisture dynamic and its impacts on evapotranspiration. Besides, the incorporation of VOD, which is closely linked to vegetation water content69,70, allowed us to detect the effect of water stress, heat stress, and vegetation phenological constraints on evaporation. Other observation-driven ET products from remote-sensing physical estimation and flux-tower are not included due to their low temporal resolution (i.e., monthly)71 or short duration72,73. ET outputs of reanalysis products are not considered in our analysis, because the assimilation systems lack explicit representation of inter-annual variability of vegetation activities and thus may not fully capture hydrological response to vegetation changes62,63,67.We used land cover maps for the year 2001 from the Moderate-Resolution Imaging Spectroradiometer (MODIS, MCD12C1 C574) based on the IGBP classification scheme to exclude water-dominated and sparely-vegetated pixels (like Sahara, Arabian Peninsula). All climate and land cover datasets mentioned above were remapped to a common 0.25° × 0.25° grid and unified to daily resolution. The main characteristics of the datasets mentioned above are summarized in Supplementary Table 1.Outputs of CMIP6 simulationsTo understand how modeled dry season changes compare with observed changes, we analyzed outputs from the “historical” (1983-2014) runs of 34 coupled models participating in the 6th Coupled Model Inter-comparison Project75 (CMIP6, Supplementary Table 3). We used these models because they offered daily outputs of all climatic variables needed for our analysis, including precipitation, latent heat (convert to E), and multiple meteorological variables for Ep (air temperature, surface specific humidity, wind speed, and net radiation). All outputs were remapped to a common 1.0° × 1.0° grid and unified to daily resolution.Defining dry season length and timingFor each grid cell and each dry season definition (P  More

The management of introduced species, whether kudzu or zebra mussels, is costly and complex. Now, a paper reports a workable, effective solution that harnesses network analyses of ecological phenomena.Invasive species can pose severe economic and environmental problems, costing more than US$1 trillion worldwide since 1970 (ref. 1). Yet managing this human-driven issue is difficult in itself. The regions involved can be vast — entire continents or countries, for instance — while budgets are typically limited. As well, the sites potentially affected and management options can be numerous. Real systems (for example, all the lakes in the United States) can have thousands of locations that could potentially be infested. By contrast, considering just 40 locations means dealing theoretically with over 1 trillion unique combinations (240) where management could be applied (for instance, to reduce the number of invasive species leaving infested areas or entering uninfested ones). Given these constraints, a key problem is how and where to deploy control measures such as invasive-species removal. While sophisticated optimization approaches exist2, which use mathematical rules to exclude most suboptimal combinations and quickly zoom in to which locations should be managed to minimize new invasions, these algorithms are generally unfeasible for very large systems. Now, writing in Nature Sustainability, Ashander et al.3 demonstrate that simpler network metrics revealing linkages between patches can provide solutions that are often comparable to the more complex optimization algorithms. More

Climate driversTo explore the impact of climate on forest resilience (see the following sections), we used monthly averaged total precipitation, 2-m air temperature, evapotranspiration deficit and surface solar radiation downwards acquired from the ERA5-Land reanalysis product at 0.1° spatial resolution for the 2000–2020 period (https://cds.climate.copernicus.eu/cdsapp#!/home). Evapotranspiration deficit was quantified as the total precipitation minus evapotranspiration. In this study, we referred to climate regions as defined by the Köppen–Geiger world map of climate classification51 (http://koeppen-geiger.vu-wien.ac.at/present.htm). The original 31 climatic zones were merged into major zones and only those characterized by vegetation cover were included in our study (tropical, arid, temperate and boreal; Extended Data Fig. 8).Vegetation dynamicsNDVI data acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument aboard the Terra satellite was used to derive changes in global vegetation for the period 2000–2020. We used cloud-free spatial composites provided at 16-day temporal resolution and 0.05° spatial resolution (MOD13C1 Version 6; https://lpdaac.usgs.gov/products/mod13c1v006/) and retained only pixels with good and marginal overall quality. The MODIS-derived NDVI dataset represents a state-of-the-art product of vegetation state whose retrieval algorithm is constantly improved52, and being derived from a unique platform and sensor, it is temporally and spatially consistent. Vegetation dynamics were analysed in terms of kNDVI, a nonlinear generalization of the NDVI based on ref. 22 and derived as follows:$$text{kNDVI=}tanh left({text{NDVI}}^{2}right)$$
(1)
kNDVI has recently been proposed as a strong proxy for ecosystem productivity that shows high correlations with both plot level measurements of primary productivity and satellite retrievals of sun-induced fluorescence22. In addition, kNDVI has been documented to be more closely related to primary productivity, to be resistant to saturation, bias and complex phenological cycles, and to show enhanced robustness to noise and stability across spatial and temporal scales compared to alternative products (for example, NDVI and near-infrared reflectance of vegetation). For these reasons, it has been retained in this study as the preferred metric to describe the state of the forest ecosystem.To obtain an accurate estimate of resilience indicators, vegetation time series need to be stationary without seasonal periodic patterns or long-term trends53. To this aim, vegetation anomalies were obtained from kNDVI data by first subtracting the multi-year 16-day sample mean and then removing linear trends from the resulting time series. Missing data, due for instance to snow cover affecting the retrieval of reflectance properties, have been gap-filled by climatological kNDVI values. The time series of kNDVI-based vegetation anomalies was used to derive resilience indicators and assess their spatial and temporal variations (see next sections).Interannual changes in vegetation were assessed in terms of growing-season-averaged kNDVI. To this end, a climatological growing season that spanned months with at least 75% of days in the greenness phase was derived from the Vegetation Index and Phenology satellite-based product54 (https://vip.arizona.edu/) and acquired for the 2000–2016 period at 0.05° spatial resolution. In addition, forest cover (FC) fraction was derived from the annual land-cover maps of the European Space Agency’s Climate Change Initiative (https://www.esa-landcover-cci.org/)55 over the 2000–2018 period at 300-m spatial resolution. FC was retrieved by summing the fraction of broadleaved deciduous, broadleaved evergreen, needle leaf deciduous and needle leaf evergreen forest. FC was resampled to 0.05° to match the kNDVI spatial resolution.Spatial patterns of slowness and its dependence on environmental factorsIn this study, we quantified the resilience of forest ecosystems—their ability to recover from external perturbations—by the use of the 1-lag TAC (refs. 3,4,5). Such an indicator was initially computed on the whole time series of vegetation anomalies (2000–2020) for forest pixels with less than 50% missing data in the original NDVI and FC greater than 0.05 and referred to in the text as long-term TAC. This analysis was used to assess the spatial patterns of the forest slowness mediated by environmental factors that affect plant growth rates and capacity to recover from perturbations. The long-term TAC was explored both in the geographic and climate space (Extended Data Fig. 1). In the climate space, long-term TAC was binned in a 50 × 50 grid as a function of average annual precipitation and temperature, both computed over the 2000–2020 period, using the average as an aggregation metric weighted by the areal extents of each record. We retained only bins with at least 50 records.To explore the potential drivers of long-term TAC, we developed an RF regression model23 and predicted the observed long-term TAC (response variable) based on a set of environmental features (predictors). The use of machine learning in general and of RF in particular, being nonparametric and nonlinear data-driven methods, does not require a priori assumptions about the functional form relating the key drivers and the response functions. The environmental variables include vegetation properties (FC and growing-season-averaged kNDVI) and climate variables (total precipitation, 2-m air temperature, evapotranspiration deficit and surface solar radiation downwards). Each of the climate variables was expressed in terms of average, coefficient of variation and 1-lag autocorrelation and resampled to 0.05° spatial resolution to match the spatial resolution of kNDVI. All environmental variables were computed annually and then averaged over time, except the autocorrelation that was computed directly for the whole period, analogously to the long-term TAC. This resulted in a set of 14 predictors representing the forest density, the background climate, the climate variability and its TAC in the observational period (Extended Data Table 1). The RF model was developed by splitting the observed long-term TAC into two separate samples: 60% of records were used for model calibration, and the remaining 40% were used to validate model performances in terms of coefficient of determination (R2), mean squared error and percentage bias (PBIAS). Each record refers to a 0.05° pixel. The RF implemented here uses 100 regression trees, whose depth and number of predictors to sample at each node were identified using Bayesian optimization. The general model formulation is as follows:$$text{TAC},=,fleft(Xright)+{varepsilon }_{{rm{f}}}$$
(2)
in which f is the RF regression model, X are the environmental predictors and εf are the residuals. We found that the model explains 87% of the spatial variance (R2) of the observed long-term TAC with a mean squared error of 0.007 and an average overestimation of 0.058 (PBIAS; Extended Data Fig. 2a). By definition, machine learning methods are not based on the mechanistic representation of the phenomena and therefore cannot provide direct information on the underlying processes influencing the system response to drivers. However, some model-agnostic methods can be applied to gain insights into the outputs of RF models. Here we used variable importance metrics to quantify and rank how individual environmental factors influence TAC (Extended Data Fig. 2b). Furthermore, using partial dependence plots derived from the machine learning algorithm RF, we explored the ecosystem response function (TAC) across gradients of vegetation and climate features (Supplementary Discussion 1 and Extended Data Fig. 2c–f).CSD indicatorsTo explore the temporal variation in forest resilience, we used CSD indicators, here quantified in terms of temporal changes in TAC retrieved for two consecutive and independent periods ranging from 2000 to 2010 and from 2011 to 2020, and assessed the significance of the change in the sampled mean aggregated for different climate regions through a two-sided t-test (Fig. 1c). This analysis was complemented by the computation of TAC on the annual scale over a 2-year lagged temporal window (3-year window size) to track the temporal changes in CSD. This resulted in a time series of TAC with an annual time step.We point out that temporal dynamics of annual TAC are driven by two processes: the changes in the resilience of the system that affect the velocity of the recovery from external perturbations and the confounding effects of the changes in autocorrelation of the climate drivers (Xac) that directly affect the autocorrelation of NDVI. Given the specific goals of this study, we factored out the second process from the total TAC signal to avoid that an increasing autocorrelation in the drivers would affect our analysis and conclusions about the resilience and the potential increase in instability56. For this purpose, we disentangled the temporal changes in TAC due to variations in autocorrelation in the climate drivers (({rm{TAC}}| {X}_{{rm{ac}}})) by adopting the space-for-time analogy and applied the RF model (f) at an annual time step (t) in a set of factorial simulations as follows:$${text{TAC}}^{t},{rm{| }},{X}_{{rm{ac}}}=fleft({X}^{t}right)-fleft({X}_{-{rm{ac}}}^{t},{X}_{{rm{ac}}}^{2000}right)$$
(3)
The first term on the right side of equation (3) is the RF model simulation obtained by accounting for the dynamics of all predictors, and the second term is the RF model simulation generated by considering all predictors dynamic except the factors of autocorrelation in climate that are kept constant to their first-year value (year 2000). For such runs, we used predictors computed on an annual scale over a 2-year lagged temporal window, consistently to the TAC time series. We found that the direct effects of autocorrelation in climate have led to a positive trend of TAC in dry zones (due to the increasing autocorrelation of the drivers in these regions) and to an opposite effect in temperate humid forests (Supplementary Fig. 3). To remove these confounding effects, the estimated term ({{rm{TAC}}}^{t}| {X}_{{rm{ac}}}) is factored out from the TACt by subtraction to derive an enhanced estimate of annual resilience that is independent of autocorrelation in climate (Extended Data Fig. 3).Long-term linear trends computed on the resulting enhanced TAC time series (δTAC) represent our reference CSD indicator used in this study to explore the changes in forest resilience. δTAC was quantified for each grid cell (Fig. 1a) and represented in the climate space following the methodology previously described (Fig. 1b). We then assessed the significance of the trends at bin level by applying a two-sided t-test for the sampled trend distributions within each bin. This significance test is independent from the structural temporal dependencies originating from the use of a 2-year lagged temporal window to compute the TAC time series.Following an analogous approach described in equation (3), we disentangled the effect of the variation in forest density, background climate and climate variability on temporal changes in TAC (Fig. 1d,e). We recognize that other environmental factors not explicitly accounted for in our RF model could play a role in modulating the temporal variations in TAC. However, given the comprehensiveness of the suite of predictors used in equation (2) (Extended Data Table 1), it seems plausible that residuals mostly reflect the intrinsic forest resilience, the component intimately connected to the short-term responses of forests to perturbations, which is not directly related to climate variability. Forest ecosystem evolutionary processes could also play a role, but longer time series would be required to reliably capture these dynamics. Furthermore, abrupt declines (ADs) in the vegetation state and following recoveries, similarly to those potentially originating from forest disturbances (for example, wildfires and insect outbreaks), could influence the TAC changes. However, such occurrences, being distributed across the globe throughout the whole period, are expected to only marginally affect the resulting trend in TAC time series.Sensitivity analysisTo assess the robustness of our results with respect to the modelling choices described above, we performed a series of sensitivity analyses for the difference in TAC retrieved for the two independent periods (2000–2010 and 2011–2020). To this aim, we tested their dependence on: the quality flag of the NDVI data used for the analyses (good, good and marginal); the gap-filling procedure tested on different periods (year and growing season); the inclusion or exclusion of forest areas affected by ADs; the threshold on the maximum percentage of missing NDVI data allowed at the pixel level (20%, 50% and 80%); the threshold on the minimum percentage of FC allowed at the pixel level (5%, 50% and 90%); and the pixel spatial resolution used for the analyses (0.05°, 0.25° and 1°). In addition, we tested the sensitivity of the trend in total TAC signal on the moving temporal window length used to calculate autocorrelation at lag 1. Results obtained for the different configurations were compared in terms of frequency distributions, separately for climate regions (Extended Data Fig. 4), and further explored in the climate space (Extended Data Figs. 5 and 6). Outcomes of the sensitivity analysis are discussed in Supplementary Discussion 2.Interplay between GPP and CSDResilience and GPP interact with each other through mutual causal links. On one hand, a reduction in forest resilience makes the system more sensitive to perturbations with potential consequent losses in GPP (ref. 26). On the other hand, a reduction in GPP may lead to a decline in resilience according to the carbon starvation hypothesis, and may be associated with increasing hydraulic failure46. To explore the link between forest resilience and primary productivity, we quantified the correlation between TAC and GPP. Estimates of GPP were derived from the FluxCom Model Tree Ensemble for the 2001–2019 period at 8-daily temporal resolution and 0.0833° spatial resolution and generated using ecosystem GPP fluxes from the FLUXNET network and MODIS remote sensing data as predictor variables36 (http://www.fluxcom.org/). Annual maps of GPP were quantified and resampled to 0.05° to match the temporal and spatial resolution of TAC time series. The Spearman rank correlation (ρ) was then computed between annual GPP and TAC over a 1° spatial moving window to better sample the empirical distribution of the two variables (Fig. 2d). The significance of ρ(GPP,TAC) was assessed over the climate space separately for each bin (Fig. 2e), similarly to the approach used to test the significance of δTAC. Furthermore, we explored the relationships between the trend in GPP (δGPP) and the trend in TAC (δTAC) by clustering the globe according to the directions of the long-term trajectories of the above-mentioned variables (Fig. 2f).Disentangling the impact of forest managementTo characterize TAC on different forest types and disentangle the potential effects originating from forest management, results were separately analysed for intact forests and managed forests. Intact forests were considered those forest pixels constituting the Intact Forest Landscapes57 dataset (https://intactforests.org/). Intact Forest Landscapes identifies the forest extents with no sign of significant human activity over the period 2000–2016 based on Landsat time series. The remaining forests pixels—not labelled as intact—were considered as managed forests (Extended Data Fig. 8). The resulting forest type map is consistent with those used for United Nations Framework Convention on Climate Change reporting58, although with more conservative estimates of intact forests in the boreal zone due to the masking based on FC and percentage of missing data applied in this study.We analysed the differences in long-term TAC (computed for the whole 2000–2020 period) between managed and intact forests by masking out the potential effect of climate background. To this aim, we compared the climate spaces generated separately for managed and intact forests by extracting only those bins that are covered by both forest classes. The resulting distributions—one for each forest class—have the same sample size, and each pair of elements shares the same climate background. Potential confounding environmental effects on average recovery rates are, therefore, minimized. We then applied a two-sided t-test for analysing the significance of the difference in the sampled means (Fig. 2a). An analogous approach was used to test the differences in δTAC and ρ(GPP,TAC) between managed and intact forests (Fig. 2b,c).Early-warning signals of abrupt forest declinesWhen forest ecosystems are subject to an extended and progressive degradation, the loss of resilience can lead to an AD (refs. 3,4,5). Such abrupt changes can trigger a regime shift (tipping point) depending on the capacity of the system to recover from the perturbations (Supplementary Methods 1 and 2). We investigated the potential of changes in TAC as early-warning signals of ADs in intact forests over the 2010–2020 period. To this aim, we quantified at the pixel level ADs as the events occurring on a certain year when the corresponding growing-season average kNDVI was more than n-times local standard deviation below the local mean. Local mean and standard deviation (σ) were computed over the 10-year antecedent temporal window (undisturbed) period and n varies between 1 and 6 with higher values reflecting more severe changes in the state of the system. For each pixel and for each fixed n value, we recorded only the first AD occurrence, thus imposing a univocal record for each abrupt change in the state of the system.We then explored whether the retrieved ADs were statistically associated with antecedent high values of δTAC. To avoid confusion with the attribution of causality, for each AD that occurred at time t (over the 2010–2020 period), we derived the δTAC over the temporal window 2000 − (t − 1). The resulting trend in δTAC is therefore antecedent and independent of the changes in vegetation associated with the AD. Then, for each pixel with an AD at time t, we also extracted randomly one of the undisturbed (with no AD) adjacent pixels and retrieved δTAC over the same temporal window. This analysis produced two distributions of δTAC associated with pixels with and without ADs (AD and no AD, respectively). The two distributions have the same size and each pair of elements shares similar background climate. We calculated the probability of occurrence of AD conditional on the trend in δTAC (({rm{AD}}| delta {rm{TAC}})) as the frequency of ADs for which (delta {rm{TAC}}left(mathrm{AD}right)| > delta {rm{TAC}}left(mathrm{no; AD}right)), and the significance of the difference in the two sampled means (AD and no AD) was evaluated through a two-sided t-test. Probability and significance were assessed for different climate regions and severity of ADs (Fig. 3a). High statistically significant probabilities suggest that the AD is following the drifting towards a critical resilience threshold plausibly associated with changes in environmental drivers.We complemented the aforementioned analyses by retrieving the tolerance and proximity to AD, hereafter determined for a 3σ severity. We first quantified the TAC that proceeded the occurrence of an AD and followed a progressive loss of resilience as captured by positive δTAC. This value, hereafter referred to as abrupt decline temporal autocorrelation (TACAD), reflects the TAC threshold over which we observed an abrupt change in the forest state (Fig. 3b). The tolerance to AD was quantified as the difference between the local TACAD and the TAC value averaged over the 2000–2009 period to characterize the pre-disturbance conditions. The tolerance metric was explored across a gradient of aridity index59 (Fig. 3c).TACAD can be directly retrieved only on those forest pixels that have already experienced an AD. As a considerable fraction of undisturbed forests could potentially be close to their critical TAC threshold, or even have already passed it, it is important to determine their TACAD. To this aim, we developed an RF regression model that expresses the TACAD as a function of the set X of environmental variables used in model f (equation (2)) but excluding the autocorrelation in climate drivers (Xreduced) already disentangled in the TAC signal. The general formulation is as follows:$${{rm{TAC}}}_{{rm{AD}}}=gleft({X}_{text{reduced}}right)+{varepsilon }_{{rm{g}}}$$
(4)
in which g is the RF regression model, Xreduced are the environmental predictors and εg are the residuals. Implementation, calibration and validation of g follow the same rationale described before for the f model. We found that the RF model explains 50% of the variance (R2) of the observed TACAD, with a mean squared error of 0.019 and an average underestimation of 0.86 (PBIAS).The RF model was then used to predict the TACAD over the whole domain of intact forests and served as input to quantify the proximity to AD of undisturbed forest pixels at the end of the observational period (year 2020). Here we defined the proximity metric as the difference between the value of TAC in 2020 and TACAD. Proximity takes negative or zero values when TACAD has already been reached (({{{rm{TAC}}}^{2020}ge {rm{TAC}}}_{{rm{AD}}})) and positive values when there are still margins before reaching the critical threshold (({{{rm{TAC}}}^{2020} < {rm{TAC}}}_{{rm{AD}}})). Together (delta {rm{TAC}} > 0) and ({{{rm{TAC}}}^{2020}ge {rm{TAC}}}_{{rm{AD}}}) therefore represent the most critical conditions, as they indicate that the critical resilience threshold for AD has already been reached and the ecosystem is continuing to lose its capacity to respond to external perturbations. We finally quantified the amount of GPP potentially exposed to such critical conditions by linearly extrapolating the GPP for the year 2020 (available GPP data stop in 2019) and overlaying it on the map of critical conditions (proximity to ({rm{AD}} < 0) and (delta {rm{TAC}} > 0)).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this paper. More