Philippa Kaur

Cliffs present a unique flora that has been little studied until now mainly because of the inherent difficulties to access this unique environment, as shown in Fig. 2. The techniques currently used to access plants on steep slopes and cliffs (e.g., abseiling, helicopter) are generally dangerous, costly and time consuming. Using a small aerial manipulator to sample plants on the cliffs can represent many advantages, including safety and portability, as well as the capability of reaching otherwise inaccessible locations easily, quickly and at low cost.Figure 2Examples of the cliff habitats of some critically endangered species on the Kauaʻi Island along with the count of known individuals as of February 2022.Full size imageHowever, several technical challenges make it difficult to develop suitable aerial manipulators for this task. Indeed, the sampling of plants on cliffs necessarily leads to significant collision risks, as well as contact forces and moments during sampling that can destabilize the drone. The samples collected would also need to be accessed from the side of the aerial platform22. Any weight (e.g., sampling tool, collected samples) located horizontally away from the center of mass of the drone creates large additional demands on the propulsion system of most drones. To collect specific plant parts in windy conditions (e.g., scion, flowers, seeds, etc.), precise and fast motion is required even in degraded Global Navigation Satellite System (GNSS) coverage near the cliffs. The great diversity of plant species and morphology found on cliffs, as well as the variety of targeted sections of plant, also represent a major design challenge. Finally, to maximize the adoption of this tool, it is also desirable that scientists with minimal training could use this platform. The next sections describe how these challenges were addressed through the development of the Mamba.Suspended sampling platformThere are a multitude of configurations that could have been explored to sample plants on cliffs. Some drones have manipulators rigidly attached to their structure20,23. However, these manipulators tend to have a limited reach to keep the center of mass within the propeller footprint and to minimize the inertia of the system. This could result in a high collision risk with the propellers in the uneven terrain found on cliffs. The contact forces created during the sampling operation also generate destabilizing moments through manipulators rigidly attached to the drone. To address these challenges, concepts involving a compliant manipulator operated from specialized drones were also explored10. Alternatively, some aerial manipulators were also passively suspended under the drone through a long rod21,24. This keeps the drone above potential obstacles within the environment, significantly reducing the operator’s mental demand and stress while also reducing the disturbances transmitted to the drone to a downward force aligned with the rod and yaw torque. To maintain these advantages while providing better precision, some projects have developed cable suspended platforms equipped with thrusters25,26. As these platforms do not have to counter gravity, the thrusters can be positioned to fight external disturbances more efficiently (e.g., wind, contact forces, drone movements). Existing systems however only stabilize the suspended platform close to its equilibrium point.The chosen concept for the Mamba, illustrated at Fig. 3, consists of a suspended platform that can stabilize itself far from its natural equilibrium to provide a large workspace. The lifting drone in this system stays safely away and above from steep cliff faces, while supporting the platform and providing rough positioning in space through better GNSS coverage. The platform is suspended 10 m below the lifting drone using four attachment points to prevent pitch and roll motions. The cable also acts as a low pass filter, isolating the platform from the fast drone movements required to fight wind disturbances. The suspended platform design can then focus on fast and precise positioning, while also being tolerant to contacts during sampling. To do so, four pairs of bidirectional actuators are used to control the motion in the plane of the pendulum (i.e., x and y translation, as well as yaw). Two pairs of actuators are installed in the x-direction to provide sufficient force to reach plants as far as 4 m from the equilibrium position. This corresponds to roughly 3.3 m from the tip of the lifting drone’s propellers.Figure 3(a) General concept of the Mamba and lifting drone during transit and sampling on cliffs. (b) Side view of the Mamba showing the components and cable installations. (c) Top view showing the antagonist thrusters configuration. (d) Close-up of the sampling tool and 2 degrees of freedom (DOF) wrist specifically designed to sample small fragile plants.Full size imageSince the Mamba is self-powered and has its own communication system, the lifting drone function is simply to lift the platform and hold it in place. This made it possible to select amongst the many commercially available products to accelerate the development of the Mamba. The DJI M300 was chosen as it comes equipped with a 360° optical obstacle avoidance vision system, an IP45 rating, and a flight time of 20 min with the Mamba attached (3.3 kg). It also advertised a four constellation GNSS receiver for better coverage around buildings, structures, and cliffs.Precise control in windsWinds under 20 km/h represent a gentle breeze on the Beaufort scale. At this level, the wind only moves the leaves, and not the branches, which allows for ideal sampling conditions. According to historical weather data from 2020, daily maximum winds are less than 20 km/h for 40 to 70% of the year, depending on the exact location on Kauaʻi Island (i.e., Lihuʻe International airport, as reported by the National Oceanic and Atmospheric Administration, and the Makaha Ridge Weather Station, as reported in the MesoWest database). This also implies that Kauaʻi experiences stronger winds on certain days which would make precise sampling difficult. Wind conditions are also more challenging near cliff faces, with increased turbulence and vertical airflow along the cliff.To allow operations on most days, while providing precise positioning and fast rejection of wind disturbances, the actuators of the Mamba are oriented in the horizontal plane. This allows the actuator forces to directly affect the motion of the suspended platform. Each actuator of the Mamba consists of a pair of brushless DC motors and 23 cm propellers capable of producing 7 N of force. The motors are installed in opposite directions, are always idling at their minimum rotation speed, and are commanded to only create force in their preferred direction. This antagonistic configuration avoids the low-velocity dead zone of a brushless motor during thrust reversal. This makes it possible to quickly revert the direction of the thrust and nearly triples the bandwidth of the actuators to approximately 2.5 Hz27. This configuration, however, comes at the expense of added mass and components.The Mamba is equipped with a flight controller that includes a control system, and a state estimator. To avoid degraded GNSS coverage issues, the state estimator only uses data from a high accuracy inertial measurement unit (IMU) to estimate the attitude of the platform. This provides the relative position of the platform with respect to the drone and is sufficient for teleoperation. Three separated proportional-derivative controllers are used for each of the DOF controlled by the actuators. This control system also provides attitude-hold assistance (i.e., pitch and roll, which correspond to x and y displacements, as well as yaw). This implies that if the user does not send any commands, the suspended platform maintains its current state.Figure 4 illustrates the stabilization accuracy of the Mamba when moving along a representative trajectory when suspended indoors from a 5.7 m cable (limited by ceiling height). This experiment confirmed that the sampling tool can maintain a position at a horizontal reach of 2.25 m with a precision of about 5 cm for 30 s. As the horizontal reach and precision are limited by the cable angular displacements (e.g., component of weight acting on the pendulum, IMU angular resolution), the resulting workspace when operating with a 10 m long cable would reach a radius of 4 m with a positioning accuracy of about 9 cm. To account for potential external disturbances like wind, the sampling tool was designed with an opening of 15 cm. This creates some margin for the pilot to align the target with the sampling mechanism. Field trials detailed below demonstrated that the Mamba actuators and controller could maintain a sufficiently stable position to sample plants in winds During the sampling phase, wind speed averaged 15.7 km/h with a standard deviation of 6.8 km/h, while wind gusts reached an average of 20.1 km/h with a standard deviation of 6.5 km/h. The maximum average wind speed recorded during sampling was 28 km/h with gusts up to 37 km/h. This represents a lower bound of the system performance, as no failure resulted from the wind conditions experienced during the trials. The a ttached Supplementary Video also demonstrates the stability of the system.Figure 4Representative motion of the sampling tool within its workspace based only on feedback from a high accuracy IMU and recorded using a motion capture system. The natural equilibrium point is at (0,0). The experiment starts with a 90° rotation around the z axis, followed by a forward movement along the x-axis of the Mamba and a lateral movement along its y-axis. The system then maintains this position for 30 s without any user inputs. Produced in MATLAB R2021a.Full size imageTeleoperated sampling of cliffs habitatsPlants growing on Kauaʻi cliffs exhibit a wide morphological variety. For this project, targets ranged from small herbaceous plants such as Euphorbia eleanoriae (plants  More

Cole, J., Findlay, S. & Pace, M. Bacterial production in fresh and saltwater ecosystems: a cross-system overview. Mar. Ecol. Prog. Ser. 43, 1–10 (1988).ADS 
Article 

Google Scholar 
Azam, F. et al. The Ecological Role of Water-Column Microbes in the Sea. Mar. Ecol. Prog. Ser. 10, 257–263 (1983).ADS 
Article 

Google Scholar 
Cotner, J. B. & Biddanda, B. A. Small players, large role: Microbial influence on biogeochemical processes in pelagic aquatic ecosystems. Ecosystems. 5, 105–121 (2002).CAS 
Article 

Google Scholar 
Falkowski, P. G., Fenchel, T. & Delong, E. F. The microbial engines that drive earth’s biogeochemical cycles. Science. 320, 1034–1039 (2008).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Newton, R. J., Jones, S. E., Eiler, A., McMahon, K. D. & Bertilsson, S. A Guide to the Natural History of Freshwater Lake Bacteria. Microbiol. Mol. Biol. Rev. 75, 14–49 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Coleman, M. L. et al. Genomic islands and the ecology and evolution of Prochlorococcus. Science. 311, 1768–1770 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Franzosa, E. A. et al. Sequencing and beyond: Integrating molecular ‘omics’ for microbial community profiling. Nat. Rev. Microbiol. 13, 360–372 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hanson, C., Fuhrman, J., Horner-Devine, M. & Martiny, J. Beyond biogeographic patterns: processes shaping the microbial landscape. Nat. Rev. Microbiol. 10, 497–506 (2012).CAS 
PubMed 
Article 

Google Scholar 
Dai, A. & Trenberth, K. E. Estimates of freshwater discharge from continents: Latitudinal and seasonal variations. J. Hydrometeorol. 3, 660–687 (2002).ADS 
Article 

Google Scholar 
White, W. R. World water: resources, usage and the role of man-made reservoirs. Report No. FR/R0012. Fundation for Water Research, (2010).Clark, E. A., Sheffield, J., van Vliet, M. T. H., Nijssen, B. & Lettenmaier, D. P. Continental runoff into the oceans (1950–2008). J. Hydrometeorol. 16, 1502–1520 (2015).ADS 
Article 

Google Scholar 
Stevaux, J. C., Paes, R. J., Franco, A. A., Mário, M. L. & Fujita, R. H. Morphodynamics in the confluence of large regulated rivers: The case of Paraná and Paranapanema Rivers. Lat. Am. J. Sedimentol. Basin Anal. 16, 101–109 (2009).
Google Scholar 
Brêda, J. P. L. F. et al. Climate change impacts on South American water balance from a continental-scale hydrological model driven by CMIP5 projections. Clim. Change 159, 503–522 (2020).ADS 
Article 

Google Scholar 
Llames, M. E. & Zagarese, H. E. Lakes and Reservoirs of South America. In Encyclopedia of Inland Waters vol.2 (ed. Linkens, G. E.). (Oxford: Elsevier, 2009).Cabrera, A. L. & Willink, A. Biogeografia De America Latina 2da edn (Organización de los Estados Americanos, 1980).Morrone, J. J. Biogeografía de América Latina y el Caribe 1st edn. (Nature, 2001).Morrone, J. J. Biogeographical regionalisation of the neotropical region. Zootaxa 3782, 1–110 (2014).PubMed 
Article 

Google Scholar 
Antonelli, A. et al. Amazonia is the primary source of Neotropical biodiversity. Proc. Natl. Acad. Sci. USA 115, 6034–6039 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Sarmento, H. New paradigms in tropical limnology: The importance of the microbial food web. Hydrobiologia 686, 1–14 (2012).Article 

Google Scholar 
Meerhoff, M. et al. Environmental Warming in Shallow Lakes. A Review of Potential Changes in Community Structure as Evidenced from Space-for-Time Substitution Approaches. Adv. Ecol. Res. 46, 259–349 (2012).Article 

Google Scholar 
Herlemann, D. P. R. et al. Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea. ISME J. 5, 1571–1579 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Metz, S. & Huber, P. et al. A georeferenced rRNA amplicon database of aquatic microbiomes from South America (Dataset), Zenodo, https://doi.org/10.5281/zenodo.6802178 (2022).Callahan, B. J., Sankaran, K., Fukuyama, J. A., McMurdie, P. J. & Holmes, S. P. Bioconductor Workflow for Microbiome Data Analysis: from raw reads to community analyses. F1000 Research 5, 1492 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17, 10 (2011).Article 

Google Scholar 
Edgar, R. C. UNOISE2: improved error-correction for Illumina 16S and ITS amplicon sequencing. Preprint at https://www.biorxiv.org/content/10.1101/081257v1 (2016).Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).CAS 
PubMed 
Article 

Google Scholar 
Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596 (2013).CAS 
PubMed 
Article 

Google Scholar 
Griffith, G. E., Omernik, J. M. & Azevedo, S. H. Ecological classification of the Western Hemisphere http://ecologicalregions.info/htm/ecoregions.htm (1998).Salcedo, J. C. R. South America: Argentina, Bolivia, and Peru https://www.worldwildlife.org/ecoregions/nt1002 Accessed (2018).Vidal, J. Geografía del Perú: las ocho regiones naturales, la regionalización transversal, la microregionalización 9th edn (PEISA, 1987).Paruelo, J. M., Beltran, A., Jobbagy, E., Sala, O. E. & Golluscio, R. A. The climate of Patagonia: General patterns and controls on biotic processes. Ecol. Austral 8, 85–101 (1998).
Google Scholar 
Iriondo, M. Quaternary lakes of Argentina. Palaeogeogr. Palaeoclimatol. Palaeoecol. 70, 81–88 (1989).Article 

Google Scholar 
Soto, D. & Campos, H. in Ecología de los bosques templados de Chile vol. 1 (eds. Khalin, J. M. & Villagrán, C.) (Editorial Universitaria, 1995).Modenutti, B. et al. Structure and dynamic of food webs in Andean North Patagonian freshwater systems: Organic matter, light and nutrient relationships. Ecol. Austral 20, 95–114 (2010).
Google Scholar 
Modenutti, B. E. et al. Structure and dynamics of food webs in Andean lakes. Lakes Reserv. Res. Manag. 3, 179–186 (1998).Article 

Google Scholar 
Quirós, R. & Drago, E. The environmental state of Argentinean lakes: An overview. Lakes Reserv. Res. Manag. 4, 55–64 (1999).Article 

Google Scholar 
Morris, D. P. et al. The attenuation of solar UV radiation in lakes and the role of dissolved organic carbon. Limnol. Oceanogr. 40, 1381–1391 (1995).ADS 
CAS 
Article 

Google Scholar 
Bastidas Navarro, M., Balseiro, E. & Modenutti, B. Bacterial Community Structure in Patagonian Andean Lakes Above and Below Timberline: From Community Composition to Community Function. Microb. Ecol. 68, 528–541 (2014).PubMed 
Article 

Google Scholar 
Modenutti, B. et al. Environmental changes affecting light climate in oligotrophic mountain lakes: The deep chlorophyll maxima as a sensitive variable. Aquat. Sci. 75, 361–371 (2013).CAS 
Article 

Google Scholar 
Bastidas Navarro, M., Martyniuk, N., Balseiro, E. & Modenutti, B. Effect of glacial lake outburst floods on the light climate in an Andean Patagonian lake: implications for planktonic phototrophs. Hydrobiologia 816, 39–48 (2018).CAS 
Article 

Google Scholar 
Sioli, H. Hydrochemistry and Geology in the Brazilian Amazon Region. Amazoniana 1, 267–277 (1968).
Google Scholar 
Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
Salati, E. & Vose, P. B. Amazon Basin: A system in equilibrium. Science. 225, 129–138 (1984).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Melack, J. M. & Forsberg, B. R. In The Biogeochemistry of the Amazon Basin Vol. 1 (eds. MacCLain, M. E., Victoria, R. & Richey, J. E.). (Oxford Scholarship Online, 2001).Junk, W. J., Bayley, P. B. & Sparks, R. E. The flood pulse concept in river-floodplain systems. Can. J. Fish. Aquat. Sci. 106, 110–127 (1989).
Google Scholar 
Ratter, J. A., Ribeiro, J. F. & Bridgewater, S. The Brazilian cerrado vegetation and threats to its biodiversity. Ann. Bot. 80, 223–230 (1997).Article 

Google Scholar 
Haridasan, M. Nutritional adaptations of native plants of the cerrado biome in acid soils. Braz. J. Plant Physiol. 20, 183–195 (2008).Article 

Google Scholar 
Vasconcelos, V., de Carvalho Júnior, O. A., de Souza Martins, É. & Couto Júnior, A. F. in World Geomorphological Landscapes. Vol. 1 (eds. Vieira, B., Salgado, A. & Santos, L.) (Springer, 2015).Bichsel, D. et al. Water quality of rural ponds in the extensive agricultural landscape of the Cerrado (Brazil). Limnology 17, 239–246 (2016).CAS 
Article 

Google Scholar 
Cunha, D. G. F., Calijuri, M., do, C. & Lamparelli, M. C. A trophic state index for tropical/subtropical reservoirs (TSItsr). Ecol. Eng. 60, 126–134 (2013).Article 

Google Scholar 
Morellato, L. P. C. & Haddad, C. F. B. Introduction: The Brazilian atlantic forest. Biotropica 32, 786–792 (2000).Article 

Google Scholar 
Galindo-Leal, C. & Câmara, I. de G. The Atlantic Forest of South America: Biodiversity status, threats, and outlook 1st edn (Island Press, 2003).Joly, C. A., Metzger, J. P. & Tabarelli, M. Experiences from the Brazilian Atlantic Forest: Ecological findings and conservation initiatives. New Phytologist 204, 459–473 (2014).PubMed 
Article 

Google Scholar 
Caliman, A. et al. Temporal coherence among tropical coastal lagoons: A search for patterns and mechanisms. Brazilian J. Biol. 70, 803–814 (2010).CAS 
Article 

Google Scholar 
Junger, P. C. et al. Salinity Drives the Virioplankton Abundance but Not Production in Tropical Coastal Lagoons. Microb. Ecol. 75, 52–63 (2018).CAS 
PubMed 
Article 

Google Scholar 
Depetris, P. J., Probst, J. L., Pasquini, A. I. & Gaiero, D. M. The geochemical characteristics of the Paraná River suspended sediment load: An initial assessment. Hydrol. Process. 17, 1267–1277 (2003).ADS 
Article 

Google Scholar 
Orfeo, O. & Stevaux, J. Hydraulic and morphological characteristics of middle and upper reaches of the Paraná River (Argentina and Brazil). Geomorphology 44, 309–322 (2002).ADS 
Article 

Google Scholar 
Neiff, J. J. Large rivers of South America: toward the new approach. Verh. Internat. Verein. Limnol 26, 167–180 (1996).
Google Scholar 
Unrein, F. Changes in phytoplankton community along a transversal section of the Lower Paraná floodplain, Argentina. Hydrobiologia 468, 123–134 (2002).Article 

Google Scholar 
Devercelli, M. Changes in phytoplankton morpho-functional groups induced by extreme hydroclimatic events in the Middle Paraná river (Argentina). Hydrobiologia 639, 5–19 (2010).CAS 
Article 

Google Scholar 
Huber, P. et al. Environmental heterogeneity determines the ecological processes that govern bacterial metacommunity assembly in a floodplain river system. ISME J. 14, 2951–2966 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Olson, D. M. et al. Terrestrial ecoregions of the world: A new map of life on Earth. Bioscience 51, 933–938 (2001).Article 

Google Scholar 
Peel, M. C., Finlayson, B. L. & McMahon, T. A. Updated world map of the Köppen-Geiger climate classification. Hydrol. Earth Syst. Sci. 11, 1633–1644 (2007).ADS 
Article 

Google Scholar 
Conde, D., Arocena, R. & Recursos, R.-G. L. acuáticos superficiales de Uruguay: ambientes, algunas problemáticas y desafios para la gestión. Ambios 10, 1–7 (2003).
Google Scholar 
Martin, L. & Suguio, K. Variation of coastal dynamics during the last 7000 years recorded in beach-ridge plains associated with river mouths: example from the central Brazilian coast. Palaeogeogr. Palaeoclimatol. Palaeoecol. 99, 119–140 (1992).Article 

Google Scholar 
Alonso, C. et al. Environmental dynamics as a structuring factor for microbial carbon utilization in a subtropical coastal lagoon. Front. Microbiol. 4, 1664–302X (2013).Article 
CAS 

Google Scholar 
Amaral, V., Graeber, D., Calliari, D. & Alonso, C. Strong linkages between DOM optical properties and main clades of aquatic bacteria. Limnol. Oceanogr. 61, 906–918 (2016).ADS 
Article 

Google Scholar 
Rennella, A. M. M., Quiro, R. & Quirós, R. The effects of hydrology on plankton biomass in shallow lakes of the Pampa Plain. Hydrobiologia 556, 181–191 (2006).Article 

Google Scholar 
Diaz, M., Pedrozo, F. & Baccala, N. Summer classification of Southern Hemisphere temperate lakes (Patagonia, Argentina). Lakes Reserv. Res. Manag. 5, 213–229 (2000).Article 

Google Scholar 
Izaguirre, I. et al. Influence of fish introduction and water level decrease on lakes of the arid Patagonian plateaus with importance for biodiversity conservation. Glob. Ecol. Conserv. 14, e00391 (2018).Article 

Google Scholar 
Porcel, S., Saad, J. F., Sabio y García, C. A. & Izaguirre, I. Microbial planktonic communities in lakes from a Patagonian basaltic plateau: influence of the water level decrease. Aquat. Sci. 81, 51 (2019).Article 
CAS 

Google Scholar 
Bernal, M. C. et al. Spatial variation of picoplankton communities along a cascade reservoir system in Patagonia, Argentina. J. Limnol. 80, 84–99 (2021).
Google Scholar 
Leinonen, R. et al. The European nucleotide archive. Nucleic Acids Res. 39, 44–47 (2011).Article 
CAS 

Google Scholar 
ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA217932 (2013).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA302313 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA294718 (2022).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA309832 (2016).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA326475 (2016).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB48609 (2022).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA289691 (2015).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA414894 (2018).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA323673 (2016).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA356055 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA310230 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA390178 (2019).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA411849 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA725228 (2021).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA292014 (2015).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA310230 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA411849 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA316315 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA406945 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA515842 (2019).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA310230 (2017).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA321235 (2016).ENA European Nucleotide Archive https://identifiers.org/ena.embl:SAMN07998328 (2015).ENA European Nucleotide Archive https://identifiers.org/ena.embl:SAMN07998330 (2015).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB36116 (2020).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB29989 (2019).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA788397 (2021).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB48353 (2022).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB37379 (2020).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB46122 (2021).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB40710 (2020).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB40864 (2020).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJEB40854 (2020).ENA European Nucleotide Archive https://identifiers.org/ena.embl:PRJNA268541 (2015). More

Incentivizing farmers to shift from conventional to regenerative practices could help fulfil the United Nations Food Systems commitments to transform food supply chains — as well as reducing carbon emissions (see L. A. Schulte et al. Nature Sustain. 5, 384–388; 2022).
Competing Interests
The authors declare no competing interests. More

Galloway, J. N., Townsend, W. H., Erisman, J. W., Bekunda, M. & Cai, Z. Transformation of the nitrogen cycle: Recent trends, questions and potential solutions. Science 320, 889–892 (2008).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Phoenix, G. K. et al. Impacts of atmospheric nitrogen deposition: Responses of multiple plant and soil parameters across contrasting ecosystems in long-term field experiments. Glob. Change Biol. 18, 1197–1215 (2012).ADS 
Article 

Google Scholar 
Pons, T. L., van der, Werf, A. & Lambers, H. Photosynthetic nitrogen use efficiency of inherently slow- and fast-growing species: Possible explanations for observed differences. In A Whole Plant Perspective on Carbon-Nitrogen Interactions (eds Roy, J., Garnier, E.) 61–77 (SPB Academic Publishing, The Hague, 1994).Ai, Z. M., Xue, S., Wang, G. L. & Liu, G. B. Responses of Non-structural carbohydrates and C:N: P stoichiometry of Bothriochloa ischaemum to nitrogen addition on the Loess Plateau, China. J. Plant Growth Regul. 36, 714–722 (2017).CAS 
Article 

Google Scholar 
Marklein, A. R. & Houlton, B. Z. Nitrogen inputs accelerate phosphorus cycling rates across a wide variety of terrestrial ecosystems. New Phytol. 193, 696–704 (2012).CAS 
PubMed 
Article 

Google Scholar 
Dietze, M. C. et al. Nonstructural carbon in woody plants. Annu. Rev. Plant Biol. 65, 667–687 (2014).CAS 
PubMed 
Article 

Google Scholar 
Hartmann, H. & Trumbore, S. Understanding the roles of nonstructural carbohydrates in forest trees -from what we can measure to what we want to know. New Phytol. 211, 386–403 (2016).CAS 
PubMed 
Article 

Google Scholar 
Yang, Q. P. et al. Different responses of non-structural carbohydrates in above-ground tissues/organs and root to extreme drought and re-watering in Chinese fir (Cunninghamia lanceolata) saplings. Trees 30, 1863–1871 (2016).CAS 
Article 

Google Scholar 
Peng, Z. T. et al. Non-structural carbohydrates regulated by nitrogen and phosphorus fertilization varied with organs and fertilizer levels in Moringa oleifera Seedlings. J. Plant Growth Regul. 40, 1777–1786 (2021).CAS 
Article 

Google Scholar 
Nardini, A. et al. Rooting depth, water relations and non-structural carbohydrate dynamics in three woody angiosperms deferentially affected by an extreme summer drought. Plant Cell Environ. 39, 618–627 (2016).CAS 
PubMed 
Article 

Google Scholar 
Wang, F. C. et al. Effects of experimental nitrogen addition on nutrients and nonstructural carbohydrates of dominant understory plants in a Chinese Fir plantation. Forests 10, 155 (2019).Article 

Google Scholar 
Elser, J. J. et al. Growth rate-stoichiometry couplings in diverse biota. Ecol. Lett. 6, 936–943 (2003).Article 

Google Scholar 
Jinm, X. M. et al. Ecological stoichiometry and biomass response of Agropyron michnoi under simulated N deposition in a sandy grassland, China. J. Arid Land. 12, 741–751 (2020).Article 

Google Scholar 
Jing, H. et al. Nitrogen addition changes the stoichiometry and growth rate of different organsin pinus tabuliformis seedlings. Front. Plant Sci. 8, 1922 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Zhan, S. X., Wang, Y., Zhu, Z. C., Li, W. H. & Bai, Y. F. Nitrogen enrichment alters plant N: P stoichiometry and intensifies phosphorus limitation in a steppe ecosystem. Environ. Exp. Bot. 134, 21–32 (2017).CAS 
Article 

Google Scholar 
Stiling, P. & Cornelissen, T. How does elevated carbon dioxide (CO2) affect plant–herbivore interactions? A field experiment and meta-analysis of CO2 -mediated changes on plant chemistry and herbivore performance. Glob. Change Biol. 13, 1823–1842 (2007).ADS 
Article 

Google Scholar 
Wang, X. G. et al. Responses of C:N: P stoichiometry of plants from a Hulunbuir grassland to salt stress, drought and nitrogen addition. Phyton-Int. J. Exp. Bot. 87, 123–132 (2018).
Google Scholar 
Liu, X. J. et al. Enhanced nitrogen deposition over China. Nature 494, 459–462 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Huang, W. J., Houlton, B. Z., Marklein, A. R., Liu, J. X. & Zhou, G. Y. Plant stoichiometric responses to elevated CO2 vary with nitrogen and phosphorus inputs: Evidence from a global-scale meta-analysis. Sci. Rep. 5, 18225 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Wang, X. et al. Effects of nutrient addition on nitrogen, phosphorus and non-structural carbohydrates concentrations in leaves of dominant plant species in a semiarid steppe. Chin. J. Ecol. 33, 1795–1802 (2014).
Google Scholar 
Yang, D. X., Song, L. & Jin, G. Z. The soil C:N: P stoichiometry is more sensitive than the leaf C:N: P stoichiometry to nitrogen addition: A four-year nitrogen addition experiment in a Pinus koraiensis plantation. Plant Soil 442, 183–198 (2019).CAS 
Article 

Google Scholar 
Chong, P. F., Zhan, J., Li, Y. & Jia, X. Y. Carbon dioxide and precipitation alter Reaumuria soongorica root morphology by regulating the levels of soluble sugars and phytohormones. Acta Physiol. Plant 41, 184 (2019).Article 
CAS 

Google Scholar 
Ma, X. Z. & Wang, X. P. Biomass partitioning and allometric relations of the Reaumuria soongorica shrub in Alxa steppe desert in NW China. For. Ecol. Manag. 468, 118–178 (2020).Article 

Google Scholar 
He, F. L., Bao, A. K., Wang, S. M. & Jin, H. X. NaCl stimulates growth and alleviates drought stress in the salt-secreting xerophyte Reaumuria soongorica. Environ. Exp. Bot. 162, 433–443 (2019).CAS 
Article 

Google Scholar 
Xu, D. H. et al. Photosynthetic parameters and carbon reserves of a resurrection plant Reaumuria soongorica during dehydration and rehydration. Plant Growth Regul. 60, 183–190 (2010).CAS 
Article 

Google Scholar 
Zhang, H. et al. miRNA–mRNA integrated analysis reveals roles for miRNAs in a typical halophyte, Reaumuria soongorica, during seed germination under salt stress. Plants 9, 351 (2020).CAS 
PubMed Central 
Article 

Google Scholar 
Bai, Y. M., Li, Y., Shan, L. S., Su, M. & Zhang, W. T. Effects of precipitation change and nitrogen addition on root morphological characteristics of Reaumuria soongorica. Arid Zone Res. 37, 1284–1292 (2020).
Google Scholar 
Hedwall, P. O., Nordin, A., Strengbom, J., Brunet, J. & Olsson, B. Does background nitrogen deposition affect the response of boreal vegetation to fertilization?. Oecologia 173, 615–624 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Wang, G., Fahey, T. J., Xue, S. & Liu, F. Root morphology and architecture respond to N addition in Pinus tabuliformis, west China. Oecologia 171, 583–590 (2013).ADS 
PubMed 
Article 

Google Scholar 
Grechi, I. et al. Effect of light and nitrogen supply on internal C: N balance and control of root-to-shoot biomass allocation in grapevine. Environ. Exp. Bot. 59, 139–149 (2007).CAS 
Article 

Google Scholar 
Xiao, L., Liu, G., Li, P. & Xue, S. Nitrogen addition has a stronger effect on stoichiometries of non-structural carbohydrates, nitrogen and phosphorus in Bothriochloa ischaemum than elevated CO2. Plant Growth Regul. 83, 325–334 (2017).CAS 
Article 

Google Scholar 
Quentin, A. G. et al. Non-structural carbohydrates in woody plants compared among laboratories. Tree Physiol. 35, 1146–1165 (2015).CAS 
PubMed 

Google Scholar 
White, L. M. Carbohydrate reserves of grasses: A review. J. Range Manag. 26, 13–18 (1973).CAS 
Article 

Google Scholar 
Millard, P., Sommerkorn, M. & Grelet, G. A. Environmental change and carbon limitation in trees: A biochemical, ecophysiological and ecosystem appraisal. New Phytol. 175, 11–28 (2007).CAS 
PubMed 
Article 

Google Scholar 
Zhang, T., Cao, Y., Chen, Y. M. & Liu, G. B. Non-structural carbohydrate dynamics in Robinia pseudoacacia saplings under three levels of continuous drought stress. Trees 29, 1837–1849 (2015).CAS 
Article 

Google Scholar 
Chapin, F. S., Schulze, E. D. & Mooney, H. A. The ecology and economics of storage in plants. Annu. Rev. Ecol. Syst. 21, 423–447 (1990).Article 

Google Scholar 
Sardans, J., Rivas-Ubach, A. & Peñuelas, J. The C:N: P stoichiometry of organisms and ecosystems in a changing world: A review and perspectives. Perspect. Plant Ecol. Evolut. Syst. 14, 33–47 (2012).Article 

Google Scholar 
Xia, J. Y. & Wan, S. Q. Global response patterns of terrestrial plant species to nitrogen addition. New Phytol. 179, 428–439 (2008).CAS 
PubMed 
Article 

Google Scholar 
Yang, Y. H., Luo, Y. Q., Lu, M., Schädel, C. & Han, W. X. Terrestrial C: N stoichiometry in response to elevated CO2 and N addition: a synthesis of two meta-analyses. Plant Soil 343, 393–400 (2011).CAS 
Article 

Google Scholar 
Mayor, J. R., Wright, S. J. & Turner, B. L. Species-specific responses of foliar nutrients to long-term nitrogen and phosphorus additions in a lowland tropical forest. J. Ecol. 102, 36–44 (2014).CAS 
Article 

Google Scholar 
Koerselman, A. & Meuleman, A. F. The vegetation ratio: A new tool to detect the nature of nutrient limitation. J. Appl. Ecol. 33, 1441–1450 (1996).Article 

Google Scholar 
Gusewell, S. N: P ratios in terrestrial plants: variation and functional significance. New Phytol. 164, 243–266 (2004).PubMed 
Article 

Google Scholar 
Wang, S., Shan, L. S., Li, Y., Zhang, Z. Z. & Ma, J. Effect of Precipitation on the Stoichiometric Characteristics of Carbon, Nitrogen and Phosphorus of Reaumuria soongarica and Salsola passerina. Acta Bot. Boreal. Occident. Sin. 40, 0335–0344 (2020).ADS 

Google Scholar 
Niu, D. C., Li, Q., Jiang, S. G., Chang, P. J. & Fu, H. Seasonal variations of leaf C:N: P stoichiometry of six shrubs in desert of China’s Alxa Plateau. Chin. J. Plant Ecol. 37, 317–325 (2013).Article 

Google Scholar 
Kleyer, M. & Minden, V. Why functional ecology should consider all plant organs: An allocation-based perspective. Basic Appl. Ecol. 16, 1–9 (2015).Article 

Google Scholar 
Yemm, E. & Willis, A. The estimation of carbohydrates in plant extracts by anthrone. Biochem. J. 57, 508–514 (1954).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bao, S. D. Soil and Agricultural Chemistry Analysis 3rd edn. (China Agriculture Press, Beijing, 2000).
Google Scholar  More

Historical land use and cover changesExisting databases differed significantly in representing historical LUCC in China (Fig. 1). Generally, datasets agree on the direction of change in cropland area until 1980 in Liu and Tian18, Ramankutty19, Houghton20, and this study (Fig. 1b, c), while the magnitude of change varied greatly. Specifically, the total cropland expansion in China was comparable between our new data set and the LUH2-GCB from 1900 onwards (56 vs 60 Mha, Fig. 1b), but cropland area changes since 1980 diverged considerably (−14 vs 41 Mha, Fig. 1c). The differences were also evident across space and more distinct during the period of 1980 to 2019 (Fig. 2a–d), in which the cropland coverage was mainly declining in our reconstructed data but increasing in LUH2-GCB (Fig. 2b, d). We found that the distinct changes are derived from the abrupt cropland increases in the FAO data reported from China, upon which LUH2-GCB was based (see Supplementary Information 3).Fig. 1: Temporal, net changes of cropland and forest from 1900 (unit: Mha).Panel a–c: cropland; panel d–f: forest; the bar charts indicate the total accumulated areas (b, e) from 1900 and (c, f) from 1980 until the last available year; LUH2-GCB was the latest version of LUH2 data used in Global Carbon Budget assessments projects (LUH2 used in MsTMIP and TRENDY were showed in Supplementary Figs. S7 and S10); Houghton data were derived from Houghton and Nassikas20 and the data in 1900 were interpolated from 1850 and 1950; Liu&Tian and Ramankutty data were derived from the works of Liu and Tian16 and Ramankutty and Foley18; the open circles indicate the changes of cropland and forest areas derived from inventory-based benchmark data; details of the benchmark data for cropland and forest were presented in Yu et al.11 and Supplementary Information 1.2 of this study, respectively; error bars: one standard deviation from the mean.Full size imageFig. 2: Spatial distribution of the fractional coverage changes of cropland and forest in China (unit: %).Panels a–d: cropland; panel e–h: forest; panels a, b, e, and f indicate the results derived from this study; data in panels c, d, g, and h were from LUH2-GCB; panels a, c, e, and g show the changes from 1900 to 1980, whereas panels b, d, f, and h show the changes from 1980 to 2019; negative and positive values indicate coverage reduction and increment, respectively.Full size imageThe problems of cropland area expansion reported to FAO are likely caused by changes in the underlying database, in which the Chinese Agricultural Yearbook (CAY) was used prior to 1996, the China Land and Resources Statistical Yearbook (LRSY) from 1996 to 2007, and the National Land and Resources Bulletin (NLRB) after 2007 (Supplementary Information 3).These three datasets are not consistent with each other because surveying methods were distinct. For example, cropland area in CAY before 1982 used an extrapolation method (i.e. “production-to-acreage” approach) due to limited field survey data11. Specifically, the extrapolation method was widely adopted for convenience and for taxation purposes in the early period, such as in the framework of the first benchmark cropland survey conducted in 1953. Such methods assumed that low-productivity cropland occupied an area of 1/3–1/8 of a predetermined, “standard-productivity” cropland21, which greatly underestimates the acreages of low productivity cropland. Biases accumulated in this method persisted until the satellite era (1980s), while the 1953 surveying data were used as the baseline for CAY to update cropland area on an annual basis.Besides the survey method, policies also contributed to a bias of reported cropland area. To tackle rising food demands, cropland expansion was highly encouraged by the government before the 1980s, implementing an incentive policy to allow new tax-free cropland without reporting to the government for the first 3–5 years22,23. Even after the initial reporting free period, these newly cultivated croplands continued to be unreported due to political incentives to show increasing crop yield to the local authorities23,24.When the first comprehensive and systematic survey (i.e. the second national cropland survey conducted during 1985–1996) was completed, the cropland area was found to be larger than previously reported in CAY11. Similarly, the shift from the use of LRSY to NLRB also introduced a spurious cropland area increment from 2007 to 2010 as small, fragmented croplands were identified by better technologies adopted in NLRB, which had remained undetected previously (Supplementary Fig. S10).Thus, LUH2-GCB has inherited spurious temporal signals of abrupt cropland increment in FAO from the 1980s to 2010 (Fig. 1a and Supplementary Fig. S10). Therefore, if the areas of other land cover types (e.g. forest) are indirectly constrained from cropland area change, cropland area biases were mirrored in the area change of other land use types. This is the case for the LUH2-GCB and for Liu and Tian’s previous land use gridded datasets. Our new database, rebuilt from Yu et al.11, corrected these problems in temporal dynamics by assimilating multiple data sources (Fig. 1a). More specifically, we retrospectively reconstructed information about cropland and forest areas year by year, using tabular data from official agencies (Supplementary Information 1 and Supplementary Data 1). To further reduce the aforementioned biases, we used the most recent and authoritative record of provincial cropland and forest areas available as the benchmark, and then spatialized the cropland and forest distributions using gridded maps as ancillary data (Supplementary Information 1). The area changes were also validated using inventory-based benchmark data (Fig. 1a, d, details were presented in Yu et al.11 and Supplementary Information 1.2).Changes in forest area in China also varied dramatically among databases. Based on Ramankutty and Foley19 and LUH2-GCB, a net forest loss was found from 1900 to the last available year, at 33–108 Mha whereas Liu and Tian18 and Houghton and Nassikas20 reported a net increase of 15 Mha (1900–2005) and 70 Mha (1900–2015) in forest area, respectively (Fig. 1d, e).By assimilating multiple source records, reports, and national surveys, however, our newly reconstructed and intensively validated database (Supplementary Figs. S4, S5, and S8) with corrected biases suggests that the forest area increased by 58 Mha from 1900 to 2019 (Fig. 1e). In particular, our data suggest that there is a surprisingly large underestimation of forest expansion in all other databases (38–102 Mha) after 1980 (Fig. 1f). We performed spatial analyses and show that widespread forest expansion in our reconstructed data was represented as a forest decline in LUH2-GCB during the period 1980–2019 (Fig. 2f, h). These existing biases in the dataset during the last four decades can be simply removed using recently available and spatially explicit forest products (Supplementary Table S2).Bias in forest change might be explained by two reasons. First, gridded datasets inherited and transferred errors from the use of FAO-based cropland dataset in developing global land use databases such as HYDE and thus LUH2-GCB8. Second, the FAO forest area reported is an important reference data used in these databases. The FAO forest area is reported based on a “land use” definition, which underestimated gross “land cover” change signals between reported years (Supplementary Information 1.3). Specifically, the FAO forest area describes lands that have been forested and will continue to be used for forestry (e.g. cut-over area, fired-over area, unestablished afforestation land) (Supplementary Table S5). This approach overestimates forest area by including lands used for reforestation where no forest was yet created. Thus, for example, the FAO statistics reported a 157.2 Mha forest area in 1990 (Supplementary Fig. S7), which is ~30 Mha higher than officially released data.More importantly, newly established forests were underestimated in such an accounting approach. The forest area expansion in China reported in the FAO statistics was 61 Mha from 1990 to 2019, which is 30 Mha lower than the officially released data16. Our reconstructed dataset, in agreement with officially released forest area, uses a “land cover” definition that characterizes the distribution of annually established forests. Therefore, the FAO statistics – a data set with definition specified to describe the area of land use – should be used with caution for constraining the temporal evolution of forest cover distribution in gridded data reconstruction, and the modeling community should be alerted to treat the LUCC data appropriately.Nonetheless, the FAO and the related LUH2 products were the dominant LUCC forcing data used in multiple studies3,25, including various process-model-based intercomparison projects (e.g. MsTMIP, LUMIP, NMIP, TRENDY), annually released Global Carbon Budget reports2,26, and IPCC reports5, implying a potential bias of these assessments for the China region. In contrast, changes in forest area from our database were independently developed (Supplementary Information 1.2), intensively calibrated, and validated using officially released national forest inventories (NFIs, see Supplementary Figs. S4 and S5), which can help to reduce the potential bias of C balance assessment in China. More specifically, the total forest area and PF area in our database were compared with historical NFIs released by the National Forestry and Grassland Administration at provincial level since 1949 (Supplementary Figs. S4 and S5), which supports the reliability of our reconstructed data.Historical carbon stock changesTo illustrate the bias in the C balance of China when using previous LUCC dataset, we performed simulations with the DLEM model for the period 1900–2019 at a resolution of 0.5 × 0.5 degree forced by our new LUCC dataset. We validated the distribution and changes of C stock using published studies and previously reported inventory-based estimations (Supplementary Information 6 and 7). The model could capture well C dynamics in China using inventory-based forest C stock changes at both provincial and national levels as the validation data set (Supplementary Fig. S14).Our results show that the total C stock decreased by 6.9 ± 0.6 Pg from 1900 to 1980 and increased by 8.9 ± 0.8 Pg C from 1980 to 2019 (Fig. 3, derived from experiment S1 in Supplementary Table S10). Such a large C stock increment since the 1980s, which is dominated by vegetation biomass C accumulation, was not captured in the MsTMIP and TRENDY projects driven by different versions of the LUH2 data (Fig. 3). This is attributed to the fast expansion of forest area(s) that was not captured by this land use forcing (Fig. 1).Fig. 3: Temporal changes of carbon storage from 1900 to 2010s in China.Panel a–c indicate vegetation carbon, soil organic carbon, and total ecosystem carbon, respectively. Results derived from experiment designed to have all environmental factors vary historically from 1900 to the 2010s, for model design details of this study see Supplementary Information 8); pink color: MsTMIP (1900–2010); blue color: TRENDY (1900–2019); dark color: this study (1900–2019); the shade areas represent the ranges of 1 standard deviation; unit: Pg C.Full size imageWe found that the large-scale forest expansion in China alone has caused a substantial C accumulation since 1980 (0.21 ± 0.006 Pg C per year, Table 1). In contrast, the forest C sink of the TRENDY models is negligible (−0.02 ± 0.05 Pg C per year, Table 1). A moderate C source (0.10 ± 0.08 Pg C per year, Table 1) was even found in the MsTMIP models, since these models were driven by continuous forest area loss and cropland expansion since the 1980s (Supplementary Fig. S7).Table 1 Comparison of reported carbon fluxes from various biomes in ChinaFull size tableA recent atmospheric inversion-based study reported that China’s land ecosystems were a large CO2 sink of −1.11 ± 0.38 Pg C per year27, which seems to be ecologically implausible and critically sensitive to the assimilation of the CO2 record from one station28. The compilation of previous studies from inventory- and satellite-based estimation, atmospheric inversion, and process-based models suggested that the Chinese C sink was much smaller (−0.18– −0.45 Pg C per year; Table 1). Our model-simulated terrestrial sink (~−0.28 ± 0.06 Pg C per year) was in this range (Table 1).While our simulated C balance in different categories or biomes is close to previous estimations, three major differences are observed (Table 1). First, because the LUCC data used in previous global models suffered from biases as shown above, the national C sink was generally underestimated in these simulations (Table 1). Second, our estimation of the forest sink is around two to three times larger than the previous one during 1949–199829. This was mainly because forest area was underestimated by over 33% (53 Mha) in the previous study29 compared to the national forest inventory (NFI)16. This underestimation may stem from exclusion of economic and bamboo forests. The third major difference is the role of grassland soils in C balance during the period 1980–2000. China’s grassland soils were previously reported as a minor sink of −0.007–−0.022 Pg C per year from the 1980s to the 2000s (Table 1), while our simulations suggest that grassland soils were a C source of 0.062–0.066 Pg C per year. This discrepancy lies in the approaches used and the accounting boundaries between studies (i.e. whether the transitions of grassland were considered), in which LUCC impacts were represented differently. For example, impervious surfaces (part of urbanized area) expanded into ~15 Mha of natural lands in China from 1978 to 201730, which further drove redistribution of cropland into marginal lands with the majority converted from grassland, causing wind erosion, habitat loss, and more water and fertilizer consumption31. Earlier studies using a static grassland map exclude the C stock loss in the land-use transition32. Thus, the distinct roles of grassland soils (i.e. sink vs source) derived from our simulations and earlier studies are not contradictory but are due to differences in accounting boundaries.LUCC impacts on carbon stock changesOur DLEM simulation indicates that LUCC induced a C loss of 5.1 ± 0.7 Pg C from 1900 to 2010s (Fig. 4a), which is substantially lower than that from MsTMIP (13.8 ± 7.7 Pg C, 1900–2010) and TRENDY (9.4 ± 3.3 Pg C, 1900–2019; Fig. 4e, f and Supplementary Fig. S18d, g). From 1980 onward, LUCC increased C storage by 4.3 ± 0.7 Pg C, with the major contribution from vegetation biomass C increment in the southwestern and northeastern regions (Fig. 4d and Supplementary Fig. S19a). Nonetheless, this C increase in biomass was not captured in MsTMIP and TRENDY models (Fig. 4e, f and Supplementary Fig. S19d, g), which simulated that LUCC continued to reduce C stock by 7.5 ± 1.6 and 5.3 ± 2.3 Pg C during the period 1980 to the 2010s, respectively (Fig. 4 and Supplementary Fig. S20).Fig. 4: Spatial distribution of LUCC impacts on ecosystem carbon storage.Panel a–c: LUCC impacts for period of 1900–2019; panel d–f: LUCC impacts for period of 1980–2019 (d–f). Panels a and d are from this study; data in panels b and e are from MsTMIP; data in panels c and f are from TRENDY; negative and positive values indicate sink and source, respectively; green and yellow bar stacked in the insert indicate LUCC impacts on vegetation and soil organic carbon in Pg C; spatial map unit: g C m−2; error bars: one standard deviation from the mean of LUCC impacts on total carbon storage.Full size imageTo confirm that such discrepancy was induced by LUCC data but not the DLEM model, we set up additional DLEM simulations using the LUH2-GCB database (Supplementary Information 8). The simulated C losses induced by LUCC when DLEM was driven with LUH2-GCB were 6.5 ± 0.4 and 11.4 ± 0.6 Pg C during the periods of 1980–2019 and 1900–2019, which are close to MsTMIP and TRENDY simulations. These results confirm that the LUCC forcing database is the major contributor to the difference between our simulations and the MsTMIP and TRENDY projects. An earlier study reported that global LUCC-induced C emissions are substantially underestimated due to underrepresented tree harvesting and land clearing from shifting cultivation33. Our simulation revealed that regional LUCC-induced C emission could also be overestimated in China due to a bias in the LUCC data.There are also disputes over whether the LUCC induced a C sink in China since the 1990s or not (Supplementary Table S8). By using an updated LUCC database, our simulations revealed that LUCC was a strong C sink in China, and that its magnitude was larger than previous estimates since the 1990s (Supplementary Table S8). Our results using an improved LUCC forcing data can facilitate narrowing down the well-known, large uncertainty in LUCC-induced C change at regional scale.Attributions of different factors on C stock changes since 1980By using the DLEM model with factorial simulations (see Supplementary Information 8 for details), we examined the direct and interactive contributions of different drivers to terrestrial C stock change in China for the period 1980–2019, including LUCC, climate, forest management, N deposition, and CO2 fertilization (see Methods, Fig. 5). Note that historical C stock change is not equivalent to the sum of factorial attributions as the baseline conditions differ (see Supplementary Information 8).Fig. 5: Attributions of different environmental factors on carbon stock change in China from 1980 to 2019.Panels a–c indicate attributions of impacts on the changes of vegetation carbon, soil organic carbon, and total ecosystem carbon, respectively; CLM: climate; CO2: rising atmospheric CO2 concentration; Ndep: N deposition; Man: forest management; Nfer: N fertilizer and manure application.Full size imageOverall, 81.9% (6.5 Pg C) of the terrestrial C sink during this period was attributed to direct impacts of all major factors, while the interactive effect contributed 18.1% (1.43 Pg C; Fig. 5c). Among all the factors examined, LUCC was the dominant driver accounting for 50.3% (3.96 Pg C) of the total C increment during the period 1980–2019 (Fig. 5c), which was largely attributed to biomass C accumulation (70.0%; Fig. 5a, c). Tian et al.13 reported that LUCC’s contribution to the sink in China was at 0.05 Pg C yr−1 since the 1980s – an amount that is only about 30% of our simulations. The discrepancy is attributed to the different representation of forest expansion in model simulations, which was 65 Mha from 1980 to 2005 in our database but only ~14 Mha in Tian et al.13. Similarly, the increase in the global land sink during the recent period (1998–2012) was also mainly attributed to LUCC (i.e. decreased tropical forest area loss and increased afforestation in northern temperate regions), instead of CO2 or climate change34.Climate change enhanced biomass C stocks by 1.63 Pg but caused a soil C loss of 0.30 Pg, thus contributing to land sink of 1.41 Pg C (18.0% of the total with all factors) since 1980 (Fig. 5). Other global change factors, such as N fertilizer application, atmospheric N deposition, and rising CO2, had a relatively minor contribution (0.1–9.54%) to the terrestrial C sink. Therefore, conversely to previous studies13,35,36,37, we showed that LUCC was the dominant driver of the recent land C sink in China, and other factors including climate change, rising CO2, and N deposition, contributed much less (0.1–18.0%) to the C stock increment in China (Fig. 5c). Tian et al.13 pointed out that LUCC effects in China should not be ignored and that the CO2 fertilization effect might be overestimated in Piao et al.38.Our simulations confirm these statements, and further show that LUCC was actually the largest contributor to land sink in China since 1980 (Fig. 5). In those studies which did not account for the influence of LUCC separately, the effects of other global change factors may have been overestimated by including LUCC impacts. For example, Chen et al.39 and He et al.37 attributed China’s C sink into different components including climate change, leaf area index (LAI) change, rising CO2, and N deposition. Such partition inevitably masked the separate contribution from LUCC, because LAI changes are closely related to land-cover changes. Thus, the accurate representation of the LUCC should be prioritized in future modeling attribution studies.Carbon stock changes in each land cover type since 1980The contribution of the establishment of young and new forest plantations to C sink has received increasing attention3,40,41,42. Our simulation (experiment S1, see Methods section) revealed that the increase in terrestrial C stock was dominantly contributed by biomass C accumulation (76.3%) (Fig. 5), in which the natural and planted forests accounted for 65% (2.9 Pg C) and 35% (1.6 Pg C) during the last four decades. We examined the LUCC effect (i.e. the largest contributor of C stock increment in Fig. 5) on the C stock of different biomes and confirmed that forest was the major contributor of the net C accumulation in China since 1980, while other biomes, including cropland, grassland, shrubland, and wetland, were relatively stable, varying from −0.3 to 0.3 Pg C during the same period (Fig. 6). A recent study documented that forest expansion was essential for a large C sink in southern China during 2002–2017, where newly-established and existing forests contributed to 32% and 34% of land C sink in the region43. In comparison to the large biomass C increase since 1980 (3.0 Pg C, Fig. 6a), the SOC increase was much lower (0.7 Pg C) during the concurrent period, although SOC changes in each biome varied greatly (–3.4–8.6 Pg C; Fig. 6b) due to area change from land conversions. The biome-level analyses further revealed that the LUCC-induced C stock increment was dominantly contributed from forest and by area expansion, while C storage in grassland and shrubland was reduced by LUCC (Fig. 6).Fig. 6: LUCC-induced carbon storage changes by land cover types based on model simulations during 1980–2019.Panel a–c indicate vegetation carbon, soil organic carbon, and total ecosystem carbon, respectively; the widths of the red blocks indicate the estimation ranges of net changes in model simulations; purple error bars indicate one standard deviation of multiple model runs; negative and positive changes indicate carbon loss and gain, respectively.Full size imageThis study highlights the dominant role of LUCC in determining the terrestrial C sink in China. Because of inaccurate representations of land cover change in China, previous estimates of the terrestrial C sink have been strongly underestimated. In contrast, forest expansion and cropland abandonment have been overestimated in the U.S., resulting in an underestimated C emission since 19807. Hence, we highlighted that the global LUCC database should be further improved, which could potentially narrow down the C imbalance reported in global C budget accounting2. In contrast to the previous studies, we showed that the contributions of factors including rising CO2, N deposition, and climate change to the land C sink in China were much smaller than LUCC over the past four decades (1980-present time). Thus, reforestation projects could represent important climate change mitigation pathways, with co-benefits for biodiversity33. To achieve the ‘C neutrality’ goal as the Chinese government declared, future climate policy should be directed to improve land management, especially forest ecosystems.Implications for future LUCC data improvementsThis study provides a novel reconstruction of recent land use change in China and assesses its implications in quantifying for terrestrial C storage dynamics. The improved dataset more accurately depicts the spatiotemporal dynamics of LUCC in China because the historically contradictory surveying records were identified, which helped to correct the biased temporal signals. Specifically, the improved surveying methods and the socioeconomic factors have greatly shaped the LUCC signals. We advocate that these impacts should be considered in the reconstruction of the national and global LUCC dataset, especially in the areas that have been intensively disturbed by human activities as is the case of China. These endeavours will be worthwhile, as demonstrated by the large impact that these bias corrections have on China’s C dynamic assessments since 1900. Thus, accurate delineation of LUCC forcing should be stressed in global simulations, including C budget accounting, biodiversity assessments, and ecosystem services evaluations. More

Intra-bull semen quality variationTo understand variation in bull semen quality, we assessed 1271 ejaculates from 79 different bulls (11 different breeds) housed at Rockhampton stud farm, in the state of Queensland, Australia, over a period of 5 years (2014–2018). The raw data, together with the semen analysis and when the samples for each individual bull were collected is available in Supplementary 1. The climate in this area (23.3786° S, 150.5089° E) is considered sub-tropical, ranging from 16 °C in winter to over 30 °C in summer. A comprehensive semen analysis was undertaken, including sperm morphology and motility. To determine the variation in semen quality, we plotted the percentage of sperm normal forms for each bull that had 5 or more ejaculates taken annually. Morphology was used as a measure of sperm quality, as Söderquist et al.17 demonstrated that sperm motility is heavily influenced by the collection/collectors and, therefore potentially unreliable and irreproducible. This resulted in the analysis of 1178 ejaculates from 50 bulls, with an average of 23 ejaculates per bull. The percentage of sperm normal forms as a box and whiskers plot for each bull is given (Fig. 1). As shown, many bulls demonstrated extremely high variation between ejaculates, with several males ranging from  70% (considered an outright “pass” in terms of cryopreservation potential) of normal sperm morphology. On the contrary, some bulls appeared to produce consistent semen samples across the year.Figure 1Changes in sperm normal forms. Semen samples were taken from bulls via electroejaculation and the percentage sperm normal forms were counted. The data show a box and whiskers plot consisting of 50 bulls, each of which had at least 5 different ejaculates across a minimum one month. Each box and whiskers plot represents an individual Bull showing the median, upper and lower quartile range. Outliers are represented by individual dots.Full size imageTo determine the amplitude and the proportion of bulls demonstrating variation in the number of normal sperm forms, we measured the difference between the maximum and the minimum values recorded for each animal. From this analysis we found that: 9 (18%) bulls showed less than 20% variation in normal forms; 15 (30%) bulls had between 20–40% variation; 13 (26%) bulls were between 40–60% and for 13 (26%) bulls this number was over 60%. These data have major implications when interpreting semen analysis, since a bull could be classified as either fertile or infertile depending on which ejaculate was considered. This data also sheds light into why correlations between the vBBSE parameters such as morphology and the bull fertility are so variable.Seasonal effect on semen qualitySeveral sources of environmental influence have been suggested to affect bull sperm quality. These include feed availability (i.e., higher conception rates in rainy seasons)27, excessive protein intake28, day length29, thermal heat stress and age30,31. To better understand the dynamics of semen quality variation within our samples, we plotted sample “pass” and “fail” cryopreservation criteria against the month of collection. A raw bull semen sample is classified as “pass” when motility is above 60% and normal forms greater than 70%. When samples were between 30 and 60% motility and 50–70% normal forms, they were classified as a “compensatory” (or qualified) pass (q-pass). The compensatory pass relied on there being the ability to have at least 10 million motile normal forms of spermatozoa in each straw to allow for conception. An outright failure was given to any sample with less than 30% progressive motility or 50% normal forms. This allowed each ejaculate to be placed into a binary “pass” or “fail”.The data for the percentage of total males that “failed” within each month (1271 ejaculates) is shown (Fig. 2A). Clearly, there is a seasonal pattern, with over 90% pass rate in winter (June–August) that fell to 50% or lower in summer (Dec-Feb). Considering that all bulls were greater than 4 years old, housed on the same stud farm and received the same dietary supplement we found no relationship in terms of “pass” or “fail” rates to these parameters. Thus, the data clearly suggested that Temperature/Temperature-Humidity or day length were responsible for the increased failure rates seen during Summer. Therefore, to understand if there was any causal relationship, we correlated either the average monthly temperature (Fig. 2B) or daylight (Fig. 2C) with monthly failure rates. The data showed a correlation with monthly temperature (r2 = 0.55; and temperature-humidity index – see further modelling below) but not with daylight hours (r2 = 0.05). Combined, these data suggest that temperature was the most likely reason for increased failure rates during the warm/hot months.Figure 2Seasonal variation in the semen quality of 1271 bull semen ejaculates. Semen samples were taken from bulls via electroejaculation and a full semen analysis was undertaken. Each sample was then classified as a pass or fail as described in Materials and Methods. (A) The percentage failure rate for each month is shown for all bulls. The number above each column indicate how many semen ejaculates were processed that month. (B). Scatter plot showing the average monthly temperature of Rockhampton and the percentage of samples that fail/month. Line of best fit indicates and r2 = 0.55. (C) Scatter plot showing the average daily sunlight in Rockhampton and the percentage of samples that fail/month. Line of best fit indicates and r2 = 0.04.Full size imageChanges in normal sperm forms categorised by breedThe present study investigated 11 different breeds of cattle, and we reasoned that maybe one, or more breed(s) contributed to failure rates more than others. Therefore, we plotted the percentage of normal forms for every ejaculate against the breed (Fig. 3). All breeds showed similar variation except for the Belmont Red, Boran and Wagyu. However, a relatively small number of bulls from the Belmont Red and Boran breeds were assessed in this study, therefore, it is unclear if they are indeed more resistant to heat. In the case of the Wagyu, it is worth mentioning that only one animal exhibited poor sperm morphology in several ejaculates (Fig. 3 circled) during winter. A close inspection of the records showed that during this time the animal had a fever episode, with body temperature reaching 39.4 °C, and that the sperm morphology returned to normal in approximately 70 days.Figure 3Variation in Semen quality as judged by Bull breed. Semen sample was collected and analysed for sperm morphology. The animals were then separated according to breed and the percentage normal forms for each ejaculate are shown.Full size imageSome bulls are heat-sensitive, whilst others are heat-tolerantAnalysis of the present data clearly illustrated that some bulls showed marked variation in terms of their semen quality throughout the year (Fig. 1). Meanwhile, others demonstrated much less variation, and were reasonably consistent. To further clarify these differences, we closely analysed the percentage of sperm morphology from two bulls, both of whom had several ejaculates were taken throughout the year, including during and after summer (Fig. 4). There was a clear pattern, and evidence of two types of bulls. Prior to the summer season, bull 1 (Fig. 4, red), designated here as “heat-sensitive”, exhibited  > 70% normal forms of spermatozoa. This value decreases dramatically, reaching its lowest point (10%) mid-January, before undergoing a recovery by April ( > 70%). In contrast, bull 2 (Fig. 4, green) showed a consistent semen profile throughout the year. The data suggest this bull was more “heat-tolerant”.Figure 4Identification of Heat-Sensitive and Heat-Tolerant bulls. The percentage normal sperm morphology from two bulls, both Droughtmasters, which had several ejaculates taken over the course of the year were plotted against the month in which the semen sample was taken. The first bull (red) is an example of a heat-sensitive bull. The second bull (Green) an example of heat-tolerant response.Full size imageTo further explore the concept of “heat-tolerant” and “heat-sensitive” bulls, we subjected 20 Wagyu bulls to a single event of controlled heat stress (40 °C, 12 h). This experiment was performed during Winter, at Singleton (New South Wales, Australia, 32.5695° S, 151.1788° E), where the average temperature was 17 °C and never exceeded 18 °C. Prior to the heat stress event, baseline semen samples were taken from each animal. After heat stress, semen samples were taken every week for 11 weeks. During the experiment, two bulls were removed from the program due to infection and sickness whilst a 3rd bull was removed as it refused to co-operate with electroejaculation procedure. From the remaining bulls, we were able to reproduce the heat-sensitive and heat-tolerant bull phenomenon. The raw data from this work is given in Supplementary 1, and an example of the data is shown (Fig. 5). For 14 bulls, we found no difference in terms of their baseline samples, which were between 70–90% normal forms. This is consistent with the Wagyu bull characteristics and their heat-tolerance (Fig. 5, yellow, green, blue lines). Within these “heat-tolerant” bulls, there was a variation of 16–22% sperm normal forms. For the other three bulls, two of them showed a decline in sperm quality, which began 2–3 weeks after the heat event, dropping from a baseline of 85% and 90% normal forms to 55% and 59%, respectively (30–31% variation in normal forms; Fig. 5, grey and orange line). The third bull showed a greater degree of heat-sensitivity. Starting at 77% morphologically normal sperm, the spermiogram of this bull illustrated a rapid decrease in normal forms in a short time (2 weeks), reaching around 40% after 4–5 weeks. Sperm morphology remained at this level (37% variation in normal form) for four weeks, before recovery. These data show that under experimental condition, the phenomenon of heat-sensitive and heat-tolerant animals can be reproduced. Further, it appears that there are degrees of heat-sensitivity.Figure 5Heating of Wagyu bulls to identify heat-sensitive and heat-tolerant effect. Twenty Wagyu bulls all 3 years of age and over were heated to 40 °C for 12 h in an insulated barn. Before heating, bassline samples were taken (week 1). After heating, electroejaculation was used to collect semen every week for 11 weeks. For every sample, sperm morphology was counted by a qualified theriogenologist. The data show the percentage normal morphology for 5 bulls. The light blue line indicates a heat-sensitive bulls, whose morphology was affected by heat, then returned back to baseline. The orange and grey line represent two related bulls (same father) who also produced less than 70% normal forms. The yellow, green and dark blue lines represent three heat-tolerant bulls, whose semen profile did not drop below the 70% normal spermatozoa threshold.Full size imageEnvironmental heat stress leads to poor sperm quality 17 days laterSimilar to previous reports, we noted that sperm quality does not begin to deteriorate until 2–3 weeks after the heat stress event of the bulls32. Based on the timing of spermatogenesis, this is consistent with reports that meiotic cells are more susceptible to heat stress following a heating event, with poor quality spermatozoa appearing in the ejaculate around 2–3 weeks later. To better understand the relationship between a “heat-event” and the production of poor-quality spermatozoa, we modelled both maximum temperature and maximum temperature humidity index (THI) and their relationship to the proportion of morphologically normal spermatozoa. The THI is an index representing the effect of humidity on the heat stress of an animal. THI was obtained using the following formula:$$mathrm{THI}=0.8* frac{{T}_{max}}{100}+frac{left(humidity*left({t}_{max}-14.4right)right)}{1}+46.4$$where Tmax = maximum temperature, (oF), and H = relative humidity.We plotted the correlation between semen quality and Tmax on the day, and every day prior (up to 40 days) to semen collection (Fig. 6). This modelling demonstrated that poor semen quality was due to maximum daytime temperature 17 days prior (Fig. 6a, arrow). Notably, 1 day of heat-stress appears to be sufficient to cause poor sperm quality, since if we take the average of 2 (Fig. 6b) or 3-day maximal temperatures prior to collection (Fig. 6c) the correlation patterns were similar. Supplementary 3 shows further modelling for Tmax and THI using between 1 and 5-day average temperatures prior to semen collection.Figure 6Bull semen quality (as percentage sperm normal forms) is related to the temperature that occurred 17–19 days ago. Correlation between sperm quality and maximum Temperature (Tmax). The Y axis is the Pearson correlation coefficient and X axis represents the number of days before the day the sperm sample was taken. (a) Uses one day of Tmax data whilst (b) averages two and (c) averages three consecutive days of Tmax data. The arrow shows the best correlation between Tmax and poor sperm quality, which occurs around 17–19 days before the semen sample is collected.Full size imageUnderstanding the temperatures at which heat-sensitive bulls failTo determine the Tmax at which bulls in the paddock begin to produce poor quality spermatozoa, we modelled data using both parameters measured at 17 days prior to the heat event, and plotted samples from 12 heat-sensitive bulls (6 Brahmans, 4 Drought Masters and 2 Santa Gertrudis). The relationship between sperm morphology and Tmax 17 days prior to heat even was plotted, with a spline smoothing cure to show the mean quality as a function of Tmax (Fig. 7a). As the temperature increase, so the quality of sperm morphology decreases as expected. To gain further clarity, we next fitted a nominal logistic regression analysis to model the proportion of spermatozoa that would either pass, Q-pass or fail sperm cryopreservation criteria as a function of Tmax 17 days prior. Tmax effect was highly significant for both outcome categories, with both p  More

Acidification of CPHistorical geochemical data suggest that the water chemistry of Cinder Pool (CP) has been relatively stable from the time of first reported geochemical data in 1947 until autumn 2018, followed by pronounced acidification between winter and spring 2019 (Supplementary Data 1, Fig. 1a, b). Images and documentation dating to even earlier (1927) reveal the presence of cinders covering ~50% of the spring surface at that time, a temperature near boiling (91.5 °C), and a description of having high sulfate and chloride levels (although data was not provided), suggesting that its chemistry has been generally stable since its discovery1. Spring pH ranged between ~3.6 and 4.5 in 22 yearly measurements spanning 71 years (1947–2018; multiple measurements in the same year were averaged to represent each year) (Fig. 1b), while the pH has been subsequently measured after 2018 as low as 2.5 (Fig. 1b). A single pH measurement of 2.5 was also recorded in a 2003 publication27, although other measurements in 2003, 2000, and 2001 were more consistent with the long-term average (i.e., pH 4.2–4.3; Supplementary Data 1). Scrutiny of chemical data accompanying the pH 2.5 measurement in 2003 indicates a SO42− concentration (~48 mg L−1) that is considerably lower than would be expected for CP, even when the pH is much higher (SO42− = 80 mg L−1; pH = 4.2–4.3). Considering that sulfuric acid is the predominant buffer of pH in these systems7,28, the pH 2.5 reading in 2003 is considered questionable. Nevertheless, the 2018 shift in pH towards more acidic conditions was accompanied by a notable change in the appearance of CP. Prior to autumn 2018, the spring waters were cloudy gray with the considerable suspension of kaolinite clay particles20 and black cinders10. However, between autumn 2018 and spring 2019, the spring waters visibly turned blue-green and contained colloidal S° particles that were also deposited along the pool shelves, while the pool also lacked its characteristic black cinders (Fig. 1a). The spring has maintained this appearance since spring 2019 until at least July 2022.Fig. 1: Historical geochemistry of Cinder Pool (CP).a Top panel shows the visual change in the appearance of CP in 2016 (left) and 2020 (right). Scale bars in the bottom right are ∼1 m. b Measurements of pH (n = 21; black line) and sulfate (SO42−) concentrations (n = 12; red line) in CP waters between 1947 and 2021. Years with multiple measurements were averaged to represent the entire year. c Paired measurements of SO42− and chloride (Cl−) concentrations (n = 12) between 1947 and 2021 in the context of the same measurements for 488 YNP springs derived from previous studies. Paired points for CP are colored based on the year they were recorded (averaged for multiple measurements/year as described above). End member fluid compositions as described in the manuscript text are indicated based on the abbreviations: MO meteoric only, HO hydrothermal only, MG meteoric plus gas, HB hydrothermal plus boiling, HBG hydrothermal plus boiling plus gas. Points for 2016, 2018, 2019, 2020, and 2021 are indicated by “16”, “18”, “19”, “20”, and “21”, respectively.Full size imageThe source of fluids in YNP hot springs can be broadly defined by concentrations of sulfate (SO42−) and chloride (Cl−)2,7. These indicators have been previously used to define the source of YNP springs as either (1) hydrothermal only (HO) waters that have moderate concentrations of SO42− (~30 mg L−1 depending on the depth of boiling; described below) but high concentrations of Cl− (~300 mg L−1), (2) meteoric-only (MO) waters containing lower concentrations of both solutes, or (3) MO waters infused with gas (MG) that have lower Cl− concentrations and higher SO42− concentrations (Fig. 1c). Subsequent boiling and/or evaporation of HO waters can concentrate Cl− and SO42− to higher concentrations (termed hydrothermal plus boiling; HB), while additional gas input into HO or HB waters can lead to particularly high concentrations of both Cl− and SO42− (hydrothermal + boiling + gas; HBG)7 (Fig. 1c). Geochemical data from surveys spanning 1947 to 2018 suggest that CP was largely sourced by hydrothermal (HO) waters that have undergone boiling and/or evaporation (HB) during this time frame (Fig. 1c).HO and HB waters are typically circumneutral7, while CP (which is also sourced by HB waters) has maintained a moderately acidic pH of ~4 until autumn 2018 (Fig. 1b). Several other low pH HB waters have been previously observed within the NGB7. The moderately acidic pH in CP (prior to 2018) has been attributed to the hydrolysis of molten S° that occurs at depths of >18 m that leads to the formation of S2O32– 11. Oxygen (O2)-dependent oxidation of S2O32−, catalyzed by trace iron sulfide in the cinders, forms SxO62− that can then react with sulfide to yield S2O32− and S° 11. Alternatively, SxO62− can be disproportionated to form S2O32− and SO42− 11. The relative rates of these reactions in CP prior to 2018 are not known although similar concentrations of S2O32− measured between 1995 and 1997 suggest that rates of S° hydrolysis and rates of S2O32− formation have been relatively constant over yearly time scales11. The consumption of O2 by reaction with S2O32− and the consumption of sulfide involving reactions with SxO62− would limit the amount of sulfuric acid that could be formed, thereby maintaining a less acidic pH than other sulfuric acid buffered acidic springs in YNP7.Between November 2018 and March 2019, the pH of CP markedly decreased to 2.8 in 2019, 2.7 in 2020, and 2.6 in 2021. This coincided with a marked increase in SO42− concentrations of ~3–5 fold above historical ranges (Fig. 1b), while Cl− concentrations fluctuated without clear trends during this time (Supplementary Fig. 1c). Thus, CP transitioned from an HB water type to an HBG water type between autumn 2018 and spring 2019 and has remained this way since (Fig. 1c). This is interpreted to reflect a substantial increase in H2S/S° oxidation that results in the formation of SO42− and H+ (sulfuric acid). Several observations suggest a fundamental restructuring of CP’s unique sulfur cycling due to dramatic physical and chemical changes at this time. As described in more detail below, the molten S° layer was detected at a depth of 18 m in 2016. However, in 2020 and 2021 there was no evidence of molten S° at ~18 to 20 m depth as previously documented, and sampling equipment could be freely dropped to a depth of 22 m (length of the cable) without interruption. In the absence of the molten S° at depth, the S° hydrolysis product S2O32−, and the cinders that catalyze SxO62− formation from S2O32− and H2S, it is possible that such reactions that previously competed for H2S or O2 (i.e., those involving S2O32− and SxO62−) are no longer taking place in CP. This in turn would allow for sulfur compounds (H2S and S°) to now be oxidized, thereby contributing to spring acidification.Alternative scenarios underlying the dramatic changes in CP waters also warrant consideration, and the three most logical are presented below. First, it is possible that the waters sourcing CP may have shifted either via replacement of the primary source or by altered mixing of multiple water sources. Water isotope values (δ2H and δ18O) can be used to further deconvolute the sources of hydrothermal waters because distinctive isotope values are associated with distinct water sources and the various influences upon them including meteoric water recharge, boiling (and/or evaporation), and water–rock interactions7,29. The water isotope values measured among the measured depths in CP in 2020 were near the range of water isotope values observed in CP across multiple months in 201613 (depth-resolved water isotope measurements were not made in 2016). The 2020 CP water isotope values were slightly right-shifted relative to those of 2016, suggesting a minor increase in the evaporation and concentration of CP water isotopes between 2016 and 20207 (Supplementary Fig. 2). These data thus do not support the hypothesis that the source of waters in CP dramatically shifted between 2016 and 2020, consistent with the SO42− and Cl− measurements indicating that the primary change to CP waters was increased input or availability of H2S for oxidation.A second alternative explanation is that a change in the water level of CP could potentially alter residence times which could allow for more oxidation of sulfur compounds in the spring and increased acidification. Such a scenario would also likely result in increased evaporation and concentration of solutes. However, the minimal increase in water isotope values (Supplementary Fig. 2) and similar Cl− concentrations (Supplementary Fig. 1c) accompanying a ~3–5 fold increase in SO42− concentration pre- and post-acidification (Fig. 1b) argue that increased residence time was of minimal importance in acidification.A third possible explanation is that a change in the plumbing system of CP is now delivering more vapor phase gas that contributes H2S and acidity when oxidized. Such a scenario could be consistent with increased surface deformation, subsurface gas accumulation, and seismic activity that has been taking place near NGB just prior to these changes21, and the transition from HB-type to HBG-type waters in CP. Sulfur species isotope analyses would help deconvolute the sources of SO42− in CP, but samples for sulfur isotopic analyses were not collected prior to acidification. Thus, it is unclear if this process may also be contributing to the acidification of CP. Regardless, the disappearance of the molten S° cap either by consumption or displacement would in effect make H2S more available for oxidation, similar to increased vapor phase input. The acidification of hot springs involves the oxidation of H2S by O230. More specifically, partial oxidation of H2S at acidic pH (90% amino acid identity to other homologs from UYS MAGs), but that was only present on unbinned contig sequences. Proteins are grouped based on their functionalities and associations in complexes. TetH (tetrathionate hydrolase), SQO sulfide:quinone oxidoreductase, SOR sulfur oxygenase reductase, SoxABCD Sulfolobus oxidase, SoxM Sulfolobus oxidase, CbsAB cytochrome b 558/566, SoxLN cytochrome ba complex, DoxBCE Desulfurolobus oxidase, DoxAD/TQOab Desulfurolobus oxidase/thiosulfate-quinone oxidoreductase, HdrAB1C1B2C2 (heterodisulfide reductase), DsrE3 DsrE3 sulfurtransferase, Dld dihydrolipoamide dehydrogenase, LplA lipoate-protein ligase A, LbpA lipoate binding protein A/glycine cleavage system H protein, TusA tRNA 2-thiouridine synthesizing protein A, SreABC sulfur reductase, SAOR sulfite:acceptor oxidoreductase, HcaLS [NiFe]-hydrogenase group 1 g. SoxEFGHI and FoxABCDEFGH (ferrous iron oxidation) gene sets were also investigated, but not identified in any of the MAGs and not shown here for brevity. A complete description of the enzymes/proteins found in individual UYS MAGs is provided in Supplementary Data 4.Full size imageTo assess the potential role of the UYS in sulfur biogeochemical cycling, the metabolic functional potentials of these populations were evaluated in greater detail based on their reconstructed genomes (Fig. 5, Supplementary Data 3). The UYS encoded the capacity for autotrophy via full complements of enzymes involved in the 3-hydroxypropionate/4-hydroxybutyrate cycle (3HP-4HB) (Supplementary Data 4), consistent with the general potential for autotrophy in most other Sulfolobales36. Consistently, the SoxM subunit that has been suggested as a marker for (facultatively) heterotrophic growth of Sulfolobales37 was absent in all UYS MAGs (Fig. 5, Supplementary Data 4). Given that all known Acidilobus and Vulcanisaeta spp. are characterized heterotrophs without known autotrophic capacity38,39, the UYS are likely the sole primary producers in the CP surface and subsurface waters, consistent with their considerable dominance in CP water communities over time.Also consistent with almost all other Sulfolobales36, the UYS universally encode the ability to reduce O2 via terminal cytochrome oxidases, although not via Sulfolobus oxidase (SoxABCD) complexes that are common among many Sulfolobales36 but rather via Desulfurolobus oxidase complexes (DoxBCE) (Fig. 5, Supplementary Data 4). An additional terminal oxidase complex (CbsAB-SoxLN) was encoded in the 2020 CP MAGs along with several other UYS MAGs from other YNP springs, although homologs of CbsAB-SoxLN were not present in the 2016 CP MAGs or several others recovered from sediments of other hot springs (Fig. 5). Thus, a potentially important metabolic difference between the pre- and post-acidification (2016 and 2020, respectively) CP Sulfolobales was the ability to use different terminal cytochrome oxidase compliments for aerobic respiration. The capacity to use multiple terminal oxidases has been suggested as an adaptation to varying oxygen tensions/availabilities37,40 that likely substantively differed between the low ORP 2016 CP waters and the high ORP 2020 CP waters (Fig. 2c). Consequently, these data point to the ecological succession of UYS strains within CP that are, at least in part, related to strain-level differences in aerobic respiration capacities.A defining feature of most cultured Sulfolobales is the ability to grow chemolithoautotrophically by coupling the oxidation of sulfur compounds (e.g., S0) to aerobic respiration37. The slow kinetics associated with abiotic oxidation of S0 with O2 at temperatures More

Ahmad, M., Saleem, M. A. & Sayyed, A. H. Efficacy of insecticide mixtures against pyrethroid-and organophosphate-resistant populations of Spodoptera litura (Lepidoptera: Noctuidae). Pest. Manag. Sci. 65, 266–274 (2009).CAS 

Google Scholar 
Shekhawat, S. S., Shafiq, A. M. & Basri, M. Effect of host plants on life table parameters of Spodoptera litura. Ind. J. Pure Appl. Biosci. 6, 324–332 (2018).
Google Scholar 
Sang, S. et al. Cross-resistance and baseline susceptibility of Spodoptera litura (Fabricius) (Lepidoptera: Noctuidae) to cyantraniliprole in the south of China. Pest Manag. Sci. 72, 922–928 (2016).CAS 
PubMed 

Google Scholar 
Ortega, D. S., Bacca, T., Silva, A. P. N., Canal, N. A. & Haddi, K. Control failure and insecticides resistance in populations of Rhyzopertha dominica (Coleoptera: Bostrichidae) from Colombia. J. Stored Prod. Res. 92, 101802 (2021).CAS 

Google Scholar 
Li, L. Pest biological control: Goals throughout my life. Annu. Rev. Entomol. 67, 1–10 (2022).PubMed 

Google Scholar 
Razaq, M., Shah, F. M., Ahmad, S. & Afzal, M. in Pest management for agronomic crops. Agronomic Crops (ed. Hasanuzzaman M.) 365–384 (Springer, 2019).Shah, F. M. & Razaq, M. in From agriculture to sustainable agriculture: Prospects for improving pest management in industrial revolution 4.0. Handbook of Smart Materials, Technologies, and Devices: Applications of Industry 4.0. Cham. (ed. Cham) 1–18 (Springer, 2020).Razaq, M. & Shah, F. M. in Biopesticides for management of arthropod pests and weeds. Biopesticides. Biopesticides Voulme 2: Advances in Bioinoculants 7–18 (Elsevier, 2022).Kishinevsky, M., Keasar, T. & Bar-Massada, A. Parasitoid abundance on plants: Effects of host abundance, plant species, and plant flowering state. Arthropod-Plant Interact. 11, 155–161 (2017).
Google Scholar 
Islam, Y. et al. Age-stage, two-sex life table and predation parameters of Harmonia axyridis Pallas (Coleoptera: Coccinellidae), reared on Acyrthosiphon pisum (Harris) (Hemiptera: Aphididae), at four different temperatures. Crop Prot. 2, 106029 (2022).
Google Scholar 
Furlong, M. J. & Zalucki, M. P. Climate change and biological control: The consequences of increasing temperatures on host–parasitoid interactions. Curr. Opin. Insect Sci. 20, 39–44 (2017).PubMed 

Google Scholar 
Islam, Y. et al. Functional response of Harmonia axyridis preying on Acyrthosiphon pisum nymphs: The effect of temperature. Sci. Rep. 11, 1–13 (2021).
Google Scholar 
Keva, O. et al. Increasing temperature and productivity change biomass, trophic pyramids and community-level omega-3 fatty acid content in subarctic lake food webs. Glo. Change Bio. 27, 282–296 (2021).ADS 
CAS 

Google Scholar 
Chi, H. et al. Age-stage, two-sex life table: An introduction to theory, data analysis, and application. Entomol. Gen. 40, 103–124 (2020).
Google Scholar 
Guedes, C. Preferência alimentar e estratégias de alimentação em Coccinellidae (Coleoptera). Oecol. Aust. 17, 59–80 (2013).
Google Scholar 
Hodek, I. & Honêk, A. Ecology of coccinellidae. Vol. 54 464 (Kulver Academic Publisher, 2013).Sutherland, A. M. & Parrella, M. P. Mycophagy in Coccinellidae: Review and synthesis. Biol. Control 51, 284–293 (2009).
Google Scholar 
Hagen, K. & Ks, H. The significance of predaceous Coccinellidae in biological and integrated control of insects. Entomophaga 7, 25–44 (1974).
Google Scholar 
Jawad, D. S., Rashid, Y. D. & Hamzah, A. G. in IOP Conference Series: Earth and Environmental Science. 012029 (IOP Publishing).Kumari, S., Suroshe, S. S., Kumar, D., Budhlakoti, N. & Yana, V. Foraging behaviour of Scymnus coccivora Ayyar against cotton mealybug Phenacoccus solenopsis Tinsley. Saudi J. Biol. Sci. 28, 3799–3805 (2021).PubMed 
PubMed Central 

Google Scholar 
Alloush, A. A. Developmental duration and predation rate of the coccidophagous coccinellid Rhyzobius lophanthae (Blaisdell) (Coleoptera: Coccinellidae) on Aspidiotus nerii Bouche. Bull. Entomol. Res. 109, 612–616 (2019).PubMed 

Google Scholar 
Koch, R., Hutchison, W., Venette, R. & Heimpel, G. Susceptibility of immature monarch butterfly, Danaus plexippus (Lepidoptera: Nymphalidae: Danainae), to predation by Harmonia axyridis (Coleoptera: Coccinellidae). Biol. Control 28, 265–270 (2003).
Google Scholar 
Islam, Y., Shah, F. M., Güncan, A., DeLong, J. P. & Zhou, X. Functional response of Harmonia axyridis to the larvae of Spodoptera litura: The combined effect of temperatures and prey instars. Front. Plant Sci. 13, 849574 (2022).PubMed 
PubMed Central 

Google Scholar 
Dixon, A. F. G. & Dixon, A. E. Insect predator-prey dynamics: ladybird beetles and biological control. (Cambridge University Press, 2000).Thompson, S. Nutrition and culture of entomophagous insects. Annu. Rev. Entomol. 44, 561–592 (1999).CAS 
PubMed 

Google Scholar 
Chaudhary, D. D., Kumar, B. & Mishra, G. Functional response in Coccinellid beetles (Coleoptera: Coccinellidae) is modified by prey-density experience. Can. Entomol. 154, 55068 (2022).
Google Scholar 
Castro, C., Almeida, L. & Penteado, S. The impact of temperature on biological aspects and life table of Harmonia axyridis (Pallas) (Coleoptera: Coccinellidae). Fla. Entomol. 94, 923–932 (2011).
Google Scholar 
Noman, Q. M., Shah, F. M., Mahmood, K. & Razaq, M. Population dynamics of Tephritid fruit flies in citrus and mango orchards of Multan, Southern Punjab, Pakistan. Pakistan J. Zool. 54, 325–330 (2021).
Google Scholar 
Eliopoulos, P. & Stathas, G. Life tables of Habrobracon hebetor (Hymenoptera: Braconidae) parasitizing Anagasta kuehniella and Plodia interpunctella (Lepidoptera: Pyralidae): Effect of host density. J. Econ. Entomol. 101, 982–988 (2008).CAS 
PubMed 

Google Scholar 
Yu, J.-Z., Chi, H. & Chen, B.-H. Comparison of the life tables and predation rates of Harmonia dimidiata (F.) (Coleoptera: Coccinellidae) fed on Aphis gossypii Glover (Hemiptera: Aphididae) at different temperatures. Biol. Control 64, 1–9 (2013).
Google Scholar 
Roy, H. E. & Ten Brown, P. M. years of invasion: Harmonia axyridis (Pallas)(Coleoptera: Coccinellidae) in Britain. Ecol. Entomol. 40, 336–348 (2015).PubMed 
PubMed Central 

Google Scholar 
Koch, R. The multicolored Asian lady beetle, Harmonia axyridis: a review of its biology, uses in biological control, and non-target impacts. J. Insect. Sci. 3, 5689 (2003).
Google Scholar 
de Castro-Guedes, C. F., de Almeida, L. M., do Rocio, C. P. S. & Moura, M. O. Effect of different diets on biology, reproductive variables and life and fertility tables of Harmonia axyridis (Pallas) (Coleoptera, Coccinellidae). Rev. Bras. Entomol. 60, 260–266 (2016).
Google Scholar 
Abdel-Salam, A. & Abdel-Baky, N. Life table and biological studies of Harmonia axyridis Pallas (Col., Coccinellidae) reared on the grain moth eggs of Sitotroga cerealella Olivier (Lep., Gelechiidae). J. Appl. Entomol. 125, 455–462 (2001).
Google Scholar 
Islam, Y. et al. Temperature-dependent functional response of Harmonia axyridis (Coleoptera: Coccinellidae) on the eggs of Spodoptera litura (Lepidoptera: Noctuidae) in laboratory. Insects 11, 583 (2020).PubMed Central 

Google Scholar 
Di, N. et al. Predatory ability of Harmonia axyridis (Coleoptera: Coccinellidae) and Orius sauteri (Hemiptera: Anthocoridae) for suppression of fall armyworm Spodoptera frugiperda (Lepidoptera: Noctuidae). Insects 12, 1063 (2021).PubMed 
PubMed Central 

Google Scholar 
Saljoqi, A.-U.-R., Khan, J. & Ali, G. Rearing of Spodoptera litura (Fabricius) on different artificial diets and its parasitization with Trichogramma chilonis (Ishii). Pak. J. Zool. 47, 1104 (2015).
Google Scholar 
Brown, P. M. et al. The global spread of Harmonia axyridis (Coleoptera: Coccinellidae): Distribution, dispersal and routes of invasion. Biocontrol 56, 623–641 (2011).
Google Scholar 
Chi, H. TWOSEX-MSChart: a computer program for the age-stage, two-sex life table analysis. Available from http://140.120.197.173/ecology/Download/TWOSEX-MSChart-B100000.rar. (2022).Chi, H. Life-table analysis incorporating both sexes and variable development rates among individuals. Environ. Entomol. 17, 26–34 (1988).
Google Scholar 
Chi, H. & Liu, H. Two new methods for the study of insect population ecology. Bull. Inst. Zool. Acad. Sin. 24, 225–240 (1985).
Google Scholar 
Goodman, D. Optimal life histories, optimal notation, and the value of reproductive value. Am. Nat. 119, 803–823 (1982).MathSciNet 

Google Scholar 
Chi, H. & Su, H.-Y. Age-stage, two-sex life tables of Aphidius gifuensis (Ashmead) (Hymenoptera: Braconidae) and its host Myzus persicae (Sulzer) (Homoptera: Aphididae) with mathematical proof of the relationship between female fecundity and the net reproductive rate. Environ. Entomol. 35, 10–21 (2006).
Google Scholar 
Tuan, S.J., Lee, C.C., Chi, H. Population and damage projection of Spodoptera litura (F.) on peanuts (Arachis hypogaea L.) under different conditions using the age-stage, two-sex life table. Pest Manag. Sci. 70, 805–813 (2014a).CAS 
PubMed 

Google Scholar 
Tuan, S.J., Lee, C.C., Chi, H. Erratum: Population and damage projection of Spodoptera litura (F.) on peanuts (Arachis hypogaea L.) under different conditions using the age-stage, two-sex life table. Pest Manag. Sci. 70, 1936 (2014b).CAS 

Google Scholar 
Chi, H. & Yang, T.-C. Two-sex life table and predation rate of Propylaea japonica Thunberg (Coleoptera: Coccinellidae) fed on Myzus persicae (Sulzer)(Homoptera: Aphididae). Environ. Entomol. 32, 327–333 (2003).
Google Scholar 
Chi, H. CONSUME-MSChart: a computer program for consumption rate analysis based on the age stage, two-sex life table analysis. http://140.120.197.173/ecology/Download/CONSUME-MSChart.rar. (2022).Akca, I., Ayvaz, T., Yazici, E., Smith, C. L. & Chi, H. Demography and population projection of Aphis fabae (Hemiptera: Aphididae): With additional comments on life table research criteria. J. Econ. Entomol. 108, 1466–1478 (2015).PubMed 

Google Scholar 
Akköprü, E. P., Atlıhan, R., Okut, H. & Chi, H. Demographic assessment of plant cultivar resistance to insect pests: A case study of the dusky-veined walnut aphid (Hemiptera: Callaphididae) on five walnut cultivars. J. Econ. Entomol. 108, 378–387 (2015).
Google Scholar 
Huang, Y. B. & Chi, H. Age-stage, two-sex life tables of Bactrocera cucurbitae (Coquillett) (Diptera: Tephritidae) with a discussion on the problem of applying female age-specific life tables to insect populations. Insect Sci. 19, 263–273 (2012).
Google Scholar 
Wei, M. et al. Demography of Cacopsylla chinensis (Hemiptera: Psyllidae) reared on four cultivars of Pyrus bretschneideri (Rosales: Rosaceae) and P. communis pears with estimations of confidence intervals of specific life table statistics. J. Econ. Entomol. 113, 2343–2353 (2020).PubMed 

Google Scholar 
Huang, H.-W., Chi, H. & Smith, C. L. Linking demography and consumption of Henosepilachna vigintioctopunctata (Coleoptera: Coccinellidae) fed on Solanum photeinocarpum (Solanales: Solanaceae): with a new method to project the uncertainty of population growth and consumption. J. Econ. Entomol. 111, 1–9 (2018).PubMed 

Google Scholar 
Chi, H.Timing of control based on the stage structure of pest populations: a simulation approach. J. Econ. Entomol. 83,
1143–1150 (1990).
Google Scholar 
Chi, H. TIMING-MSChart: a computer program for the population projection based on age-stage, two-sex life table. (http://140.120.197.173/Ecology/Download/TIMING-MSChart.rar). (2022).Mignault, M.-P., Roy, M. & Brodeur, J. Soybean aphid predators in Quebec and the suitability of Aphis glycines as prey for three Coccinellidae. BioControl 51, 89–106 (2006).
Google Scholar 
Brown, M. Intraguild responses of aphid predators on apple to the invasion of an exotic species, Harmonia axyridis. BioControl 48, 141–153 (2003).
Google Scholar 
Pervez, A., Chandra, S. & Kumar, R. Effect of dietary history on intraguild predation and cannibalism of ladybirds’ eggs. Int. J. Trop. Insect Sci. 41, 2637–2642 (2021).
Google Scholar 
Lundgren, J. G. Nutritional aspects of non-prey foods in the life histories of predaceous Coccinellidae. Biol. Control 51, 294–305 (2009).
Google Scholar 
Yu, J.Z. et al. Demography and mass-rearing Harmonia dimidiata (Coleoptera: Coccinellidae) using Aphis gossypii (Hemiptera: Aphididae) and eggs of Bactrocera dorsalis (Diptera: Tephritidae). J. Econ. Entomol. 111, 595–602 (2018).PubMed 

Google Scholar 
De Oliveira, R. T., dos Santos-Cividanes, T. M., Cividanes, F. J. & da Conceic, L. Harmonia axyridis Pallas (Coleoptera: Coccinellidae): Biological aspects and thermal requirements. Adv. Entomol. 2014, 5589 (2014).
Google Scholar 
Ali, S. et al. Using a two-sex life table tool to calculate the fitness of Orius strigicollis as a predator of Pectinophora gossypiella. Insects 11, 275 (2020).PubMed Central 

Google Scholar 
Merene, Y. Population dynamics and damages of onion thrips (Thripstabaci)(Thysanoptera: Thripidae) on onion in Northeastern Ethiopia. J. Entomol. Nematol. 7, 1–4 (2015).
Google Scholar 
Mou, D. F., Lee, C. C., Smith, C. & Chi, H. Using viable eggs to accurately determine the demographic and predation potential of Harmonia dimidiata (Coleoptera: Coccinellidae). J. Appl. Entomol. 139, 579–591 (2015).
Google Scholar 
Farhadi, R., Allahyari, H. & Chi, H. Life table and predation capacity of Hippodamia variegata (Coleoptera: Coccinellidae) feeding on Aphis fabae (Hemiptera: Aphididae). Biol. Control 59, 83–89 (2011).
Google Scholar 
Hance, T., van Baaren, J., Vernon, P. & Boivin, G. Impact of extreme temperatures on parasitoids in a climate change perspective. Annu. Rev. Entomol. 52, 107–126 (2007).CAS 
PubMed 

Google Scholar 
Ma, X., Zhu, J., Yan, W. & Zhao, C. Projections of desertification trends in Central Asia under global warming scenarios. Sci. Total Environ. 781, 146777 (2021).ADS 
CAS 
PubMed 

Google Scholar  More