Ecology

Brown MV, Philip GK, Bunge JA, Smith MC, Bissett A, Lauro FM, et al. Microbial community structure in the North Pacific ocean. ISME J. 2009;3:1374–86.CAS 
Article 

Google Scholar 
Chénard C, Wijaya W, Vaulot D, dos Santos AL, Martin P, Kaur A, et al. Temporal and spatial dynamics of Bacteria, Archaea and protists in equatorial coastal waters. Sci Rep. 2019;9:1–13.Article 

Google Scholar 
Yeh YC, McNichol J, Needham DM, Fichot EB, Berdjeb L, Fuhrman JA. Comprehensive single‐PCR 16S and 18S rRNA community analysis validated with mock communities, and estimation of sequencing bias against 18S. Environ Microbiol. 2021;23:2340–3250.Article 

Google Scholar 
Needham DM, Fuhrman JA. Pronounced daily succession of phytoplankton, archaea and bacteria following a spring bloom. Nat. Microbiol. 2016;1:16005.CAS 
Article 

Google Scholar 
Parada AE, Needham DM, Fuhrman JA. Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ Microbiol. 2016;18:1403–14.CAS 
Article 

Google Scholar 
Needham DM, Fichot EB, Wang E, Berdjeb L, Cram JA, Fichot CG, et al. Dynamics and interactions of highly resolved marine plankton via automated high-frequency sampling. ISME J. 2018;12:2417.CAS 
Article 

Google Scholar 
Steinberg DK, Carlson CA, Bates NR, Johnson RJ, Michaels AF, Knap AH. Overview of the US JGOFS Bermuda Atlantic Time-series Study (BATS): a decade-scale look at ocean biology and biogeochemistry. Deep Sea Res Part II Top Stud Oceanogr. 2001;48:1405–47.CAS 
Article 

Google Scholar 
Karl DM, Church MJ. Microbial oceanography and the Hawaii Ocean Time-series programme. Nat Rev Microbiol. 2014;12:699–713.CAS 
Article 

Google Scholar 
Mestre M, Höfer J, Sala MM, Gasol JM. Seasonal variation of bacterial diversity along the marine particulate matter continuum. Front Microbiol. 2020;11:1590.Article 

Google Scholar 
Cram JA, Chow C-ET, Sachdeva R, Needham DM, Parada AE, Steele JA, et al. Seasonal and interannual variability of the marine bacterioplankton community throughout the water column over ten years. ISME J. 2015;9:563–80.Article 

Google Scholar 
Berelson WM. The flushing of two deep‐sea basins, southern California borderland. Limnol Oceanogr. 1991;36:1150–66.CAS 
Article 

Google Scholar 
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.CAS 
Article 

Google Scholar 
Traving SJ, Kellogg CT, Ross T, McLaughlin R, Kieft B, Ho GY, et al. Prokaryotic responses to a warm temperature anomaly in northeast subarctic Pacific waters. Commun Biology. 2021;4:1–12.Article 

Google Scholar 
Butler TM, Wilhelm A-C, Dwyer AC, Webb PN, Baldwin AL, Techtmann SM. Microbial community dynamics during lake ice freezing. Scient Rep. 2019;9:1–11.
Google Scholar 
LeBrun ES, King RS, Back JA, Kang S. Microbial community structure and function decoupling across a phosphorus gradient in streams. Microb Ecol. 2018;75:64–73.CAS 
Article 

Google Scholar 
McNichol J, Berube PM, Biller SJ, Fuhrman JA. Evaluating and improving small subunit rRNA PCR primer coverage for bacteria, archaea, and eukaryotes using metagenomes from global ocean surveys. Msystems. 2021;6:e00565–21.CAS 
Article 

Google Scholar 
De Bie T, De Meester L, Brendonck L, Martens K, Goddeeris B, Ercken D, et al. Body size and dispersal mode as key traits determining metacommunity structure of aquatic organisms. Ecol Lett. 2012;15:740–7.Article 

Google Scholar 
Soininen J, Korhonen JJ, Karhu J, Vetterli A. Disentangling the spatial patterns in community composition of prokaryotic and eukaryotic lake plankton. Limnol. Oceanogr. 2011;56:508–20.Article 

Google Scholar 
Wu W, Lu H-P, Sastri A, Yeh Y-C, Gong G-C, Chou W-C, et al. Contrasting the relative importance of species sorting and dispersal limitation in shaping marine bacterial versus protist communities. ISME J. 2018;12:485–94.Article 

Google Scholar 
Kraemer S, Ramachandran A, Colatriano D, Lovejoy C, Walsh DA. Diversity and biogeography of SAR11 bacteria from the Arctic Ocean. ISME J. 2020;14:79–90.Article 

Google Scholar 
Tsementzi D, Wu J, Deutsch S, Nath S, Rodriguez-R LM, Burns AS, et al. SAR11 bacteria linked to ocean anoxia and nitrogen loss. Nature. 2016;536:179–83.CAS 
Article 

Google Scholar 
Brown MV, Lauro FM, DeMaere MZ, Muir L, Wilkins D, Thomas T, et al. Global biogeography of SAR11 marine bacteria. Mol Syst Biol. 2012;8:595.Article 

Google Scholar 
Giovannoni SJ. SAR11 bacteria: the most abundant plankton in the oceans. Ann Rev Mar Sci. 2017;9:231–55.Article 

Google Scholar 
Thrash JC, Temperton B, Swan BK, Landry ZC, Woyke T, DeLong EF, et al. Single-cell enabled comparative genomics of a deep ocean SAR11 bathytype. ISME J. 2014;8:1440–51.Article 

Google Scholar 
Fernández-Gomez B, Richter M, Schüler M, Pinhassi J, Acinas SG, González JM, et al. Ecology of marine Bacteroidetes: a comparative genomics approach. ISME J. 2013;7:1026–37.Article 

Google Scholar 
Countway PD, Vigil PD, Schnetzer A, Moorthi SD, Caron DA. Seasonal analysis of protistan community structure and diversity at the USC Microbial Observatory (San Pedro Channel, North Pacific Ocean). Limnol Oceanogr. 2010;55:2381–96.Article 

Google Scholar 
Kim DY, Countway PD, Jones AC, Schnetzer A, Yamashita W, Tung C, et al. Monthly to interannual variability of microbial eukaryote assemblages at four depths in the eastern North Pacific. ISME J. 2014;8:515–30.Article 

Google Scholar 
Parris DJ, Ganesh S, Edgcomb VP, DeLong EF, Stewart FJ. Microbial eukaryote diversity in the marine oxygen minimum zone off northern Chile. Front Microbiol. 2014;5:543.Article 

Google Scholar 
Orsi W, Song YC, Hallam S, Edgcomb V. Effect of oxygen minimum zone formation on communities of marine protists. ISME J. 2012;6:1586–601.CAS 
Article 

Google Scholar 
Leibold MA, Holyoak M, Mouquet N, Amarasekare P, Chase JM, Hoopes MF, et al. The metacommunity concept: a framework for multi‐scale community ecology. Ecol Lett. 2004;7:601–13.Article 

Google Scholar 
Fuhrman JA, Hewson I, Schwalbach MS, Steele JA, Brown MV, Naeem S. Annually reoccurring bacterial communities are predictable from ocean conditions. Proc Natl Acad Sci USA. 2006;103:13104–9.CAS 
Article 

Google Scholar 
Chow C-ET, Sachdeva R, Cram JA, Steele JA, Needham DM, Patel A, et al. Temporal variability and coherence of euphotic zone bacterial communities over a decade in the Southern California Bight. ISME J. 2013;7:2259–73.CAS 
Article 

Google Scholar 
Parada AE, Fuhrman JA. Marine archaeal dynamics and interactions with the microbial community over 5 years from surface to seafloor. ISME J. 2017;11:2510–25.Article 

Google Scholar 
Mestre M, Ruiz-González C, Logares R, Duarte CM, Gasol JM, Sala MM. Sinking particles promote vertical connectivity in the ocean microbiome. Proc Natl Acad Sci USA. 2018;115:E6799–E807.CAS 
Article 

Google Scholar 
Mestre M, Ferrera I, Borrull E, Ortega‐Retuerta E, Mbedi S, Grossart HP, et al. Spatial variability of marine bacterial and archaeal communities along the particulate matter continuum. Mol Ecol. 2017;26:6827–40.CAS 
Article 

Google Scholar 
Wilson B, Müller O, Nordmann E-L, Seuthe L, Bratbak G, Øvreås L. Changes in marine prokaryote composition with season and depth over an Arctic polar year. Front Mar Sci. 2017;4:95.
Google Scholar 
Treusch AH, Vergin KL, Finlay LA, Donatz MG, Burton RM, Carlson CA, et al. Seasonality and vertical structure of microbial communities in an ocean gyre. ISME J. 2009;3:1148–63.Article 

Google Scholar 
Zinger L, Amaral-Zettler LA, Fuhrman JA, Horner-Devine MC, Huse SM, Welch DBM, et al. Global patterns of bacterial beta-diversity in seafloor and seawater ecosystems. PLoS ONE. 2011;6:e24570.CAS 
Article 

Google Scholar 
DeLong EF, Preston CM, Mincer T, Rich V, Hallam SJ, Frigaard N-U, et al. Community genomics among stratified microbial assemblages in the ocean’s interior. Science. 2006;311:496–503.CAS 
Article 

Google Scholar 
Agogué H, Lamy D, Neal PR, Sogin ML, Herndl GJ. Water mass‐specificity of bacterial communities in the North Atlantic revealed by massively parallel sequencing. Mol. Ecol. 2011;20:258–74.Article 

Google Scholar 
Walsh EA, Kirkpatrick JB, Rutherford SD, Smith DC, Sogin M, D’Hondt S. Bacterial diversity and community composition from seasurface to subseafloor. ISME J. 2016;10:979–89.Article 

Google Scholar 
Reji L, Tolar BB, Chavez FP, Francis CA. Depth-differentiation and seasonality of planktonic microbial assemblages in the Monterey Bay upwelling system. Front Microbiol. 2020;11:1075.Article 

Google Scholar 
Milici M, Vital M, Tomasch J, Badewien TH, Giebel HA, Plumeier I, et al. Diversity and community composition of particle‐associated and free‐living bacteria in mesopelagic and bathypelagic Southern Ocean water masses: Evidence of dispersal limitation in the Bransfield Strait. Limnol Oceanogr. 2017;62:1080–95.Article 

Google Scholar 
Crespo BG, Pommier T, Fernández‐Gómez B, Pedrós‐Alió C. Taxonomic composition of the particle‐attached and free‐living bacterial assemblages in the Northwest Mediterranean Sea analyzed by pyrosequencing of the 16S rRNA. Microbiologyopen. 2013;2:541–52.CAS 
Article 

Google Scholar 
Ganesh S, Parris DJ, DeLong EF, Stewart FJ. Metagenomic analysis of size-fractionated picoplankton in a marine oxygen minimum zone. ISME J. 2014;8:187–211.CAS 
Article 

Google Scholar 
Ghiglione J-F, Galand PE, Pommier T, Pedrós-Alió C, Maas EW, Bakker K, et al. Pole-to-pole biogeography of surface and deep marine bacterial communities. Proc Natl Acad Sci USA. 2012;109:17633–8.CAS 
Article 

Google Scholar 
Murillo AA, Ramírez-Flandes S, DeLong EF, Ulloa O. Enhanced metabolic versatility of planktonic sulfur-oxidizing γ-proteobacteria in an oxygen-deficient coastal ecosystem. Front Mar Sci. 2014;1:18.Article 

Google Scholar 
Hawley AK, Nobu MK, Wright JJ, Durno WE, Morgan-Lang C, Sage B, et al. Diverse Marinimicrobia bacteria may mediate coupled biogeochemical cycles along eco-thermodynamic gradients. Nat Commun. 2017;8:1–10.CAS 
Article 

Google Scholar 
Santoro AE, Buchwald C, McIlvin MR, Casciotti KL. Isotopic signature of N2O produced by marine ammonia-oxidizing archaea. Science. 2011;333:1282–5.CAS 
Article 

Google Scholar 
Aldunate M, De la Iglesia R, Bertagnolli AD, Ulloa O. Oxygen modulates bacterial community composition in the coastal upwelling waters off central Chile. Deep Sea Res Part II Top Stud Oceanogr. 2018;156:68–79.CAS 
Article 

Google Scholar 
Duret MT, Lampitt RS, Lam P. Prokaryotic niche partitioning between suspended and sinking marine particles. Environ Microbiol Rep. 2019;11:386–400.CAS 
Article 

Google Scholar 
Lindh MV, Sjöstedt J, Andersson AF, Baltar F, Hugerth LW, Lundin D, et al. Disentangling seasonal bacterioplankton population dynamics by high‐frequency sampling. Environ Microbiol. 2015;17:2459–76.Article 

Google Scholar 
Teeling H, Fuchs BM, Bennke CM, Krueger K, Chafee M, Kappelmann L, et al. Recurring patterns in bacterioplankton dynamics during coastal spring algae blooms. elife. 2016;5:e11888.Article 

Google Scholar 
Buchan A, LeCleir GR, Gulvik CA, González JM. Master recyclers: features and functions of bacteria associated with phytoplankton blooms. Nat Rev Microbiol. 2014;12:686–98.CAS 
Article 

Google Scholar 
Giovannoni SJ, Tripp HJ, Givan S, Podar M, Vergin KL, Baptista D, et al. Genome streamlining in a cosmopolitan oceanic bacterium. Science. 2005;309:1242–5.CAS 
Article 

Google Scholar 
Cram JA, Xia LC, Needham DM, Sachdeva R, Sun F, Fuhrman JA. Cross-depth analysis of marine bacterial networks suggests downward propagation of temporal changes. ISME J. 2015;9:2573–86.Article 

Google Scholar 
Salazar G, Cornejo‐Castillo FM, Borrull E, Díez‐Vives C, Lara E, Vaqué D, et al. Particle‐association lifestyle is a phylogenetically conserved trait in bathypelagic prokaryotes. Mol Ecol. 2015;24:5692–706.Article 

Google Scholar 
Mohit V, Archambault P, Toupoint N, Lovejoy C. Phylogenetic differences in attached and free-living bacterial communities in a temperate coastal lagoon during summer, revealed via high-throughput 16S rRNA gene sequencing. Appl Environ Microbiol. 2014;80:2071–83.Article 

Google Scholar 
Rieck A, Herlemann DP, Jürgens K, Grossart H-P. Particle-associated differ from free-living bacteria in surface waters of the Baltic Sea. Front Microbiol. 2015;6:1297.Article 

Google Scholar 
Pachiadaki MG, Brown JM, Brown J, Bezuidt O, Berube PM, Biller SJ, et al. Charting the complexity of the marine microbiome through single-cell genomics. Cell. 2019;179:1623–35. e11.CAS 
Article 

Google Scholar 
Giner CR, Pernice MC, Balagué V, Duarte CM, Gasol JM, Logares R, et al. Marked changes in diversity and relative activity of picoeukaryotes with depth in the world ocean. ISME J. 2020;14:437–49.Article 

Google Scholar 
Countway PD, Gast RJ, Dennett MR, Savai P, Rose JM, Caron DA. Distinct protistan assemblages characterize the euphotic zone and deep sea (2500 m) of the western North Atlantic (Sargasso Sea and Gulf Stream). Environ Microbiol. 2007;9:1219–32.CAS 
Article 

Google Scholar 
Ollison GA, Hu SK, Mesrop LY, DeLong EF, Caron DA. Come rain or shine: depth not season shapes the active protistan community at station ALOHA in the North Pacific Subtropical Gyre. Deep Sea Res Part I Oceanogr Res Pap. 2021;170:103494.Article 

Google Scholar 
Schnetzer A, Moorthi SD, Countway PD, Gast RJ, Gilg IC, Caron DA. Depth matters: microbial eukaryote diversity and community structure in the eastern North Pacific revealed through environmental gene libraries. Deep Sea Res Part I Oceanogr Res Pap. 2011;58:16–26.Article 

Google Scholar 
Martin P, Allen JT, Cooper MJ, Johns DG, Lampitt RS, Sanders R, et al. Sedimentation of acantharian cysts in the Iceland Basin: strontium as a ballast for deep ocean particle flux, and implications for acantharian reproductive strategies. Limnol Oceanogr. 2010;55:604–14.CAS 
Article 

Google Scholar 
Lampitt R, Salter I, Johns D. Radiolaria: Major exporters of organic carbon to the deep ocean. Glob Biogeochem Cycle. 2009;23:GB1010.Article 

Google Scholar 
Skovgaard A, Massana R, Balague V, Saiz E. Phylogenetic position of the copepod-infesting parasite Syndinium turbo (Dinoflagellata, Syndinea). Protist. 2005;156:413–23.CAS 
Article 

Google Scholar 
Bachvaroff TR, Kim S, Guillou L, Delwiche CF, Coats DW. Molecular diversity of the syndinean genus Euduboscquella based on single-cell PCR analysis. Appl Environ Microbiol. 2012;78:334–45.CAS 
Article 

Google Scholar 
Guillou L, Viprey M, Chambouvet A, Welsh R, Kirkham A, Massana R, et al. Widespread occurrence and genetic diversity of marine parasitoids belonging to Syndiniales (Alveolata). Environ Microbiol. 2008;10:3349–65.CAS 
Article 

Google Scholar 
Berdjeb L, Parada A, Needham DM, Fuhrman JA. Short-term dynamics and interactions of marine protist communities during the spring–summer transition. ISME J. 2018;12:1907.Article 

Google Scholar 
Hu SK, Connell PE, Mesrop LY, Caron DA. A hard day’s night: diel shifts in microbial eukaryotic activity in the north pacific subtropical gyre. Front Mar Sci. 2018;5:351.Article 

Google Scholar 
Clarke LJ, Bestley S, Bissett A, Deagle BE. A globally distributed Syndiniales parasite dominates the Southern Ocean micro-eukaryote community near the sea-ice edge. ISME J. 2019;13:734–7.CAS 
Article 

Google Scholar 
De Vargas C, Audic S, Henry N, Decelle J, Mahé F, Logares R, et al. Eukaryotic plankton diversity in the sunlit ocean. Science. 2015;348:1261605.Article 

Google Scholar 
Pernice MC, Giner CR, Logares R, Perera-Bel J, Acinas SG, Duarte CM, et al. Large variability of bathypelagic microbial eukaryotic communities across the world’s oceans. ISME J. 2016;10:945–58.Article 

Google Scholar 
Crutsinger GM, Collins MD, Fordyce JA, Gompert Z, Nice CC, Sanders NJ. Plant genotypic diversity predicts community structure and governs an ecosystem process. Science. 2006;313:966–8.CAS 
Article 

Google Scholar 
Hawkins BA, Porter EE. Does herbivore diversity depend on plant diversity? The case of California butterflies. Am Nat. 2003;161:40–9.Article 

Google Scholar 
Scherber C, Eisenhauer N, Weisser WW, Schmid B, Voigt W, Fischer M, et al. Bottom-up effects of plant diversity on multitrophic interactions in a biodiversity experiment. Nature. 2010;468:553–6.CAS 
Article 

Google Scholar 
Yang JW, Wu W, Chung C-C, Chiang K-P, Gong G-C, Hsieh C-h Predator and prey biodiversity relationship and its consequences on marine ecosystem functioning—interplay between nanoflagellates and bacterioplankton. ISME J. 2018;12:1532–42.Article 

Google Scholar 
Fuhrman JA, Comeau DE, Hagström Å, Chan AM. Extraction from natural planktonic microorganisms of DNA suitable for molecular biological studies. Appl Environ Microbiol. 1988;54:1426–9.CAS 
Article 

Google Scholar 
Lie AA, Kim DY, Schnetzer A, Caron DA. Small-scale temporal and spatial variations in protistan community composition at the San Pedro Ocean Time-series station off the coast of southern California. Aquat Microb Ecol. 2013;70:93–110.Article 

Google Scholar 
Yeh Y-C, Needham DM, Sieradzki ET, Fuhrman JA. Taxon disappearance from microbiome analysis reinforces the value of mock communities as a standard in every sequencing run. MSystems. 2018;3:e00023–18.Article 

Google Scholar 
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucl Acids Res. 2013;41:D590–D6.CAS 
Article 

Google Scholar 
Guillou L, Bachar D, Audic S, Bass D, Berney C, Bittner L, et al. The Protist Ribosomal Reference database (PR2): a catalog of unicellular eukaryote small sub-unit rRNA sequences with curated taxonomy. Nucl Acids Res. 2013;41:D579–D604.
Google Scholar 
Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.CAS 
Article 

Google Scholar 
Decelle J, Romac S, Stern RF, Bendif EM, Zingone A, Audic S, et al. Phyto REF: a reference database of the plastidial 16S rRNA gene of photosynthetic eukaryotes with curated taxonomy. Mol Ecol Resour. 2015;15:1435–45.CAS 
Article 

Google Scholar 
Oksanen J, Blanchet F, Friendly M, Kindt R, Legendre P, McGlinn D, et al. Vegan: community ecology package. R package version 2.5-7. 2020. https://CRAN.R-project.org/package=vegan.Wickham H. ggplot2-elegant graphics for data analysis. Cham, Switzerland: Springer International Publishing; 2016.Kolde R. Pheatmap: pretty heatmaps. R Package Version. 2012;1:726.
Google Scholar 
Schloerke B, Crowley J, Cook D. Package ‘GGally’. Extension to ‘ggplot2’See. 2018;713. More

China is the world’s most populous country, with a population of over 1.4 billion people, or 18% of the world human population1. However, China has only about 9% of the 1.4 billion hectares of total arable land in the world2. The question of “who will feed China?” raised by Dr. Lester R. Brown in 1995 is still worthy of consideration today, and ensuring food security remains a top priority for the Chinese government3.Rice is the staple food on dining-tables of over 65% of the population in China; thus, adequate rice production is critical to ensure food security in China4. In order to produce more rice on the limited amount of arable land available, double-season rice cropping systems, which involve successively growing early-season rice (ESR) and late-season rice (LSR) from March to November within a single calendar year, have been extensively developed in southern China5. The development of double-season rice cropping systems has made a considerable contribution toward achieving rice self-sufficiency in China6.Hunan and Jiangxi are the top two double-season rice producing provinces in China7. However, in recent years, the planting area devoted to double-season rice has sharply decreased in the Hunan-Jiangxi region as a result of the conversion from double- to single-season rice (SSR) cropping systems (referred as the rice “double-to-single” phenomenon) (Fig. 1A). A reduced rural labor supply and rising labor wages due to urbanization and economic growth are the key driving forces for the rice “double-to-single” phenomenon11. Fortunately, the rice “double-to-single” phenomenon has not resulted in a decrease in total rice production in the Hunan-Jiangxi region (Fig. 1B). During the most recent 10 years (2011–2020), the total rice production in the Hunan-Jiangxi region has been ranged from 45.3 to 48.7 million tons (Mt) with an average of 46.6 Mt, and the contribution of the Hunan-Jiangxi region to rice production in China has been maintained at ~ 22%.Figure 1(A) Planting areas (million hectares, Mha) for early-season rice (ESR), late-season rice (LSR), and single-season rice (SSR) in the Hunan-Jiangxi region and (B) total rice production (million tons, Mt) in the Hunan-Jiangxi region and the contribution of the Hunan-Jiangxi region to total rice production in China from 2011 to 2020. In (B), the dashed line represents the average rice production during 2011–2020. The rice planting area and total rice production in the Hunan-Jiangxi region were calculated based on data collected from the Hunan Provincial Bureau of Statistics8 and the Jiangxi Provincial Bureau of Statistics9. The contribution of the Hunan-Jiangxi region to rice production in China is the percentage of total rice production in the Hunan-Jiangxi region to the total rice production in China. Data for total rice production in China were collected from the National Bureau of Statistics of China10.Full size imageBecause China’s population is still growing12, China must continue to increase rice production. The domestic demand for rice grain in China is expected to reach 217 Mt by 2030, when the population of China is expected to stabilize6. To meet this demand, the Hunan-Jiangxi region will need to produce 47.7 Mt of rice grains, assuming that the contribution of the Hunan-Jiangxi region to rice production in China remains at the level of the most recent 10 years (~ 22%) (Fig. 1B). This expected rice production (ERP) is 1.1 Mt higher than the average total rice production during the most recent 10 years. In order to avoid the negative effect of the “double-to-single” phenomenon on achieving the ERP in the Hunan-Jiangxi region by 2030, it is necessary to estimate how much planting area of double-season rice will be needed in this region by this point in time.The ERP can be expressed by the following formula: ERP = EPAESR × EGYESR + EPALSR × EGYLSR + EPASSR × EGYSSR, where EGYESR, EGYLSR, and EGYSSR are the estimated grain yields of ESR, LSR, and SSR, respectively; and EPAESR, EPALSR, and EPASSR are the estimated planting areas for ESR, LSR, and SSR, respectively. We assume the following conditions in the Hunan-Jiangxi region by 2030: (1) the total paddy field area will be maintained in the range of 4.57–5.02 million hectares (Mha) that was planted during the years 2011–20208,9; (2) the ratio of EPALSR to EPAESR is the same as the average ratio of planting area of LSR to ESR during 2011–2020 (i.e., 1.07) (Fig. 1A); (3) EPASSR is the difference between the total paddy field area and the EPALSR; and (4) EGYESR, EGYLSR, and EGYSSR are projected under three scenarios: (1) constant yield scenario, (2) 5% yield increase scenario, and (3) 10% yield increase scenario (Fig. 2). The baseline yield for all three scenarios is the average grain yields during 2011–2020. The EPAESR, EPALSR, and EPASSR in the Hunan-Jiangxi region needed to achieve the expected rice production by 2030 were obtained by solving the above formula and are shown in Fig. 3.Figure 2(A) Grain yields from 2011 to 2020 and (B) estimated grain yields by 2030 under three scenarios for early-season rice (ESR), late-season rice (LSR), and single-season rice (SSR) in the Hunan-Jiangxi region. The grain yields from 2011 to 2020 were calculated based on data collected from the Hunan Provincial Bureau of Statistics8 and the Jiangxi Provincial Bureau of Statistics9. In (A), ns denotes non-significant trend at the 0.05 probability level (Statistix 8.0, Analytical software, Tallahassee, FL, USA). In (B), the baseline yield for all three scenarios is the average grain yields during 2011–2020.Full size imageFigure 3Estimated planting areas for early-season rice (ESR), late-season rice (LSR), and single-season rice (SSR) in the Hunan-Jiangxi region that will be required to achieve the expected rice production by 2030 under three scenarios: (A) constant yield scenario, (B) 5% yield increase scenario, and (C) 10% yield increase scenario. Mha is million hectares.Full size imageThe results presented in Fig. 3 provide guidance and models for the government’s decision-making process in the planning planting areas for ESR, LSR, and SSR in the Hunan-Jiangxi region. In brief, farmers will need to plant 2.55–3.18 Mha of ESR, 2.73–3.40 Mha of LSR, and 1.17–2.29 Mha of SSR under the constant yield scenario, 2.09–2.72 Mha of ESR, 2.24–2.91 Mha of LSR, and 1.66–2.78 Mha of SSR under the 5% yield increase scenario, and 1.67–2.31 Mha of ESR, 1.79–2.47 Mha of LSR, and 2.10–3.23 Mha of SSR under the 10% yield increase scenario in the Hunan-Jiangxi region by 2030 depending on the total available paddy field area.One thing to note here is that the actual planting areas for ESR (2.44 Mha) and LSR (2.57 Mha) in 2020 are below the estimated lower limits of planting areas for ESR (2.55 Mha) and LSR (2.73 Mha) that will be needed by 2030 under the constant yield scenario, while the actual planting area for SSR in 2020 (2.42 Mha) is above the estimated upper limit for SSR (2.29 Mha) that will be needed by 2030 under the constant yield scenario (Figs. 1A and 3A). This finding indicates that it is urgent to avoid a further aggravated “double-to-single” phenomenon while maintaining the total paddy field area in the Hunan-Jiangxi region. Because it is not an easy task to maintain the total paddy field area under the projected scenario for urban expansion13, the government should prepare an alternative to reverse the “double-to-single” phenomenon in the Hunan-Jiangxi region. Increasing the mechanized level of farming operation and improving economic returns to farmers are two key aspects for the government to take into account to promote the development of double-season rice.Although the current planting area of double-season rice can fully meet the requirement for achieving the ERP in the Hunan-Jiangxi region by 2030 under both the 5% and 10% yield increase scenarios, there is some difficulty in reaching the yield increase targets. In recent years, the planting area of high-quality rice varieties has been rapidly increased in China14. However, grain yield is generally not very high for high-quality rice varieties, although no genetic linkage has been identified between grain yield and quality in rice15. Hence, great efforts are required to develop rice varieties with both high quality and high yield. In addition, rice yields are determined not only by the variety but also by environments and crop management practices. Soil nutrient deficiencies, unfavorable climatic conditions (e.g., heat, cold, and drought), and pest infestations have always been major yield-limiting factors for rice production in China16. Therefore, great efforts are also required to: (1) improve soil fertility of low- and medium-yielding rice fields and optimize nutrient management practices; (2) develop climate-smart agriculture practices for alleviating climatic stresses; and (3) promote integrated pest management practices. More

Experimental setupThis research was approved by the Carnegie Mellon University Committee of the Use of Human Subjects. All methods were performed in accordance with the relevant guidelines and regulations. Informed consent was obtained from all participants. Our data includes no identifying information of human participants. We conducted experiments from February to August 2021 (except for the preliminary sessions of random information; we ran the condition from June to November 2020). We preregistered the main experiment settings using AsPredicted (https://aspredicted.org/sm4k5.pdf).A total of 2786 subjects participated in our incentivized decision-making game experiments. We recruited subjects using Amazon Mechanical Turk (MTurk)52,53. Supplementary Table 1 shows the subject demographics. Our participants interacted anonymously over the Internet using customized software playable in a browser window (available at http://breadboard.yale.edu). All participants provided explicit consent and passed a series of human verification checks and a screening test of understanding game rules and payoffs before playing the game (see SI). We prohibited subjects from participating in more than one session of the experiment by using unique identifications for each subject on MTurk.In each session, subjects were paid a $2.00 show-up fee and a bonus depending on whether they took the appropriate disaster decision in four rounds. Furthermore, subjects earned $1.00 when they completed all four rounds. In each round, when a disaster stroke before they evacuated, the subjects earned no bonus. Otherwise, they earned a bonus of $1.00 without disaster or $0.50 with disaster by spending $0.50 for evacuation, plus $0.05 per other players who took the correct action accordingly (Supplementary Table 2). We have confirmed with prior work that the amount of evacuation cost, if any, makes no significant difference in the game’s performance23.At the start, subjects were required to pass a series of human verification checks. They needed to pass Google’s reCAPTCHA using the “I’m not a robot” checkbox. They were also requested to answer whether they were human players. The exact question asked was: “Please select an applicable answer about you.” The options were: “I am not a bot. I am a real person.” “I am not a real person. I am a bot.” “I am anything but a human.” and “I am a computer program working for a person.” The option’s order was randomized. Only the participants who selected “I am not a bot. I am a real person.” moved to the step of informed consent.When subjects provided explicit consent, they were asked to take a tutorial before the actual game would begin. In the tutorial, each subject separately interacted with three dummy players in two rounds of a 45-s practice game. In the actual game, some subjects would be informed in advance whether a disaster would indeed strike or not. In the practice game, while all subjects were not informed of such information in the first round, they were informed of the information in the second round. Thus, they practiced both conditions in terms of prior information on the disaster (see SI).After the practice game, subjects were assessed for their comprehension of the game rules and payment structure using four multiple-choice questions with three options. If they failed to select the correct answer in one of the questions, they could reselect it only once through the entire test. If they failed to select the correct answer more than once, they were unable to join the actual game.At 720 s after the tutorial beginning, a “Ready” button became visible simultaneously to all the subjects who completed the tutorial and passed the comprehension tests. The actual games started 30 s after the “Ready” button showed up. If subjects did not click the button before the game started, they were dropped. The game required a certain number of subjects. When the subjects who successfully clicked the button were more than 16, surplus subjects, randomly selected, were dropped from the game. When the number of qualified subjects was less than 12, the game did not start. As a result, subjects started the game in a group with an average size of 15.5 (s.d. = 1.1).At the start of the actual game, we selected one subject (the “informant”) at random who was informed in advance whether a disaster would indeed strike or not. The other subjects were informed that some players had accurate information about the disaster, but they were not informed who the informant was. The exact sentence that the informants received in their game screen was “A disaster is going to strike!” when a disaster would strike or “There is no disaster.” when a disaster would not strike. The one that the other uninformed subjects received was “A disaster may or may not strike.” Then, the group had the same informant across the four rounds except for a supplement condition of random informants. In the random informant condition, an informant was randomly selected every round.To prevent an end-of-game effect, we randomly set the game time with a normal distribution of a mean of 75 s and a standard deviation of 10 s. Prior work has confirmed that the game time is sufficient for players to communicate and make an evacuation decision23. As a result, each round ended at 75.0 s on average (s.d. = 9.5) without prior notice. In half of the sessions, a disaster struck at the end of the game. We did not inform any subjects, including the informants, when their sessions would end, the global network structure they were embedded in, or how many informants were in the game. After making their evacuation choice, subjects were informed of their success and failure along with overall results in their group. Then, subjects played another round of the evacuation game until they completed four total rounds. They had the same local network environment across four rounds except for the dynamic network condition.Network structure and tie rewiringIn the network sessions, subjects played the game in a directed network with a random graph configuration. A certain number of ties were present at the game’s onset as the initial density was set to 0.25.In the dynamic network conditions, subjects also could change their neighbors by making or breaking ties between rounds. In the tie-rewiring step, 40% of all the possible subject pairs were chosen at random. Thus, subjects could choose every other player at least once throughout the entire session (i.e., a set of four rounds) with a probability of about 80%. When the chosen pairs were connected, the pairs (the ties) were dissolved if the predecessor subject of the directed ties chose to break the tie. When the chosen pairs were not connected, the pairs (the ties) were newly created when the predecessor of the potential tie chose to create the tie. Subjects were not informed of the rewiring rate.To equalize the game time, we made subjects in the independent and static network conditions wait for additional 10 s after each game round ended. Despite the adjustment, the game time was significantly longer in the dynamic network sessions than in the independent and static network sessions. The average game time is 429.5 s (s.d. = 20.2) for the independent condition; 428.8 s (s.d. = 19.0) for the static network condition; and 564.7 s (s.d. = 36.3) for the dynamic network condition.To clarify mechanisms for dynamic networks to facilitate collective intelligence, we added one supplementary condition. In the supplementary condition, subjects were assigned to one of the 40 isomorphic networks that other subjects had developed with tie-rewiring options through the three rounds in the dynamic network condition (567 subjects in 40 groups). Network structure and other game settings (i.e., whether a disaster stroke, how long the game was, and which node was the informant) were identical to where the others played the game at the final round. However, players were different, and they had no prior experience in the game. They played the game in a network with a topology created by others ostensibly to optimize the accurate flow of information. In contrast to other conditions, subjects played only one round in the isomorphic network condition.Signal buttonsDuring the game of network sessions, subjects were allowed to share information about the possibly impending “disaster” by using “Safe” and “Danger” buttons that indicated their assessment (see SI). The default node color was grey. Then, when they clicked the Safe button, their node turned blue and, after 5 s, automatically returned to grey. Likewise, the Danger button turned their node to red for 5 s. Subjects could see only the colors of neighbors to whom they were directly connected. Since the signal exchange occurred through directed connections, an individual could send, but not receive, information from another subject (and vice versa). Once subjects chose to evacuate, they could no longer send signals, and their node showed grey (the default color) for the rest of the game. The neighbors of evacuated subjects were not informed of their evacuation. We have confirmed with prior work that collective performance does not vary with the communication continuity and the evacuation visibility23. Subjects could use the Safe and Danger buttons any time unless they evacuated, or they did not have to.Players dropping during the gameAfter each game round, when a player was inactive for 10 s, they were warned about being dropped. When they remained inactive after 10 s, they were dropped. When the selected informant was dropped, the session stopped at the round, and we did not use the data. Furthermore, as too many dropped players could affect the network structure and the behavioral dynamics of remaining players, we did not use the sessions where more than 25% of initial players were dropped during the game. Overall, 4 players dropped in 15 sessions; 3 players dropped in 22 sessions; 2 players dropped in 41 sessions; 1 player dropped in 44 sessions; and no player dropped in 58 sessions. The dropped players were prohibited from joining another session of this experiment.As noted above, players took the additional tie-rewiring step every round in the dynamic network sessions. Thus, the total game time was longer in the dynamic network sessions than in the independent and static network sessions even with the adjustment. As a result, more players were dropped in the dynamic network sessions than in the independent and static network sessions. The average number of dropped players across the four rounds is 0.40 (s.d. = 0.60) for the independent condition; 1.15 (s.d. = 0.86) for the static network condition; 1.75 (s.d. = 1.19) for the dynamic network condition. Although group size could affect collective performance, we found the differences in group size small enough for our study. We have confirmed the dynamic network’s performance improvement with a comprehensive analysis controlling the effect of group size (Supplementary Table 3). Also, there was no statistically significant difference in the dropped players’ performance of the dynamic network condition, compared with the other two conditions. The rate of correct actions of dropped players is 0.456 (s.d. = 0.322) for the independent condition, 0.594 (s.d. = 0.387) for the static network condition, and 0.558 (s.d. = 0.411) for the dynamic network condition; P = 0.106 between the independent condition and the dynamic network condition; P = 0.599 between the static network condition and the dynamic network condition (Welch two-sample t test).Analysis of signal diffusionsTo examine the change in signal diffusion, we analyzed “diffusion chains” for each signal type in the network sessions. We first identified the subjects who sent a signal when their neighbors had never sent one as spontaneous “diffusion sources.” When a subject sent a signal after at least one neighbor had sent the same type of signal, we regarded the subject’s signaling (and evacuation with danger signals) as occurring in a chain of signal diffusion and the total number of the responded subjects (including the diffusion source) as the diffusion size.We analyzed the distribution of signal diffusion chains with complementary cumulative distribution functions, measuring the fraction of diffusion chains that exhibit a given number of diffusion sizes. We found that the number of diffusions of both signals did not change across rounds. Safe-signal diffusions were more likely to occur than danger-signal diffusions regardless of whether a “disaster” would strike and how many rounds subjects played. On the other hand, the diffusion size varied greatly across rounds in disaster situations. With “disaster,” false safe signals spread further than true danger signals at the first round, but after that, warnings outperformed safe signals in terms of diffusion size. Figure 2B and Supplementary Fig. 3 scrutinize the changes in diffusion chains with their distributions.Analysis of individual responsivenessWe analyzed how individual evacuation behavior varies with exposure to signals from neighbors54. Let$${a}_{i}^{evacuate},, (t)=left{begin{array}{ll}1&quad text{if subject } i text{ evacuates at time } t\ 0&quad text{otherwise}end{array}right.$$$${a}_{i}^{show, safe},, (t)=left{begin{array}{ll}1&quad text{if subject } i mathrm{ shows a safe signal at time } t\ 0&quad text{otherwise}end{array}right.$$$${a}_{i}^{show , danger} ,, (t)=left{begin{array}{ll}1&quad text{if subject } i text{ shows a danger signal at time } t\ 0&quad text{otherwise}end{array}right.$$The hazard function, or instantaneous rate of occurrence of subject (i)’s evacuation at time t, is defined as:$${lambda }_{i},, (t)=underset{mathit{dt}to 0}{{mathrm{lim}}}frac{{mathrm{Pr}}({a}_{i}^{evacuate}=1;,, tt)}{dt}$$To model the time to evacuation, We used a Cox proportional hazards model with time-varying covariates for the number of signals, incorporating an individual actor-specific random effect55:$${lambda }_{i} ,, left{t|{{P}_{i}, X}_{i}(t), {G}_{i},{Y}_{i}(t)right}={lambda }_{0}(t)mathrm{exp}left{{{beta }_{P}^{{prime}}{P}_{i}+beta }_{X}^{{prime}}{X}_{i}(t)+{beta }_{G}^{{prime}}{G}_{i}+{beta }_{Y}^{{prime}}{Y}_{i}(t)+{gamma }_{i}right}$$where λ0(t) is a baseline hazard at time t.In the model, the hazard λi(t) depends on the covariates Pi, Xi(t), Gi, and Yi(t). The covariate Pi is the vector of subject i’s experiences before the sessions; that is, the number of rounds, the number of disasters that she has experienced, and the number of disasters that she has been struck by.The covariate Xi(t) is the vector of the number of safe signals ({x}_{i}^{safe} (t)), the number of danger signals ({x}_{i}^{danger} (t)). When subject j is a neighbor of subject i (i.e., (jin {N}_{i})), subject i is exposed to the signal of subject j, so that:$${x}_{i}^{safe},, (t)=sum_{jin {N}_{i}}{a}_{j}^{show, safe}(t)$$$${x}_{i}^{danger},, (t)=sum_{jin {N}_{i}}{a}_{j}^{show, danger}(t)$$The covariate Gi is the vector of the properties of the network in which subject i is embedded, out-degree, in-degree, and a network plasticity indicator. The covariate Yi(t) is the vector of the number of the subject i’s actions of sending safe and danger signals before time t. The coefficients β are the fixed effects and γi is the random effect for individual i. We assumed that waiting times to evacuation in different actors are conditionally independent given the sequence of signals they receive from network neighbors. This model shows how the hazard of an individual’s evacuation depends on the signaling actions of others, their network position, and experience (Supplementary Table 4). We applied the same model to the first signaling behavior. More

Wang, C. & Yang, J. Rhizospheric and heterotrophic components of soil respiration in six Chinese temperate forests. Glob. Change Biol. 13, 123–131 (2007).ADS 
Article 

Google Scholar 
Zhao, X., Li, L., Xie, Z. & Li, P. Effects of nitrogen deposition and plant litter alteration on soil respiration in a semiarid grassland. Sci. Total Environ. 740, 1–10 (2020).Article 

Google Scholar 
Jia, X., Shao, M. & Wei, X. Responses of soil respiration to N addition, burning and clipping in temperate semiarid grassland in northern China. Agr. For. Meteorol. 166, 32–40 (2012).Article 

Google Scholar 
Meyer, N., Meyer, H. & Welp, G. Soil respiration and its temperature sensitivity (Q10): rapid acquisition using mid-infrared spectroscopy. Geoderma 323, 31–40 (2018).ADS 
Article 

Google Scholar 
Gao, Q. et al. Effects of litter manipulation on soil respiration under short-term nitrogen addition in a subtropical evergreen forest. For. Ecol. Manag. 429, 77–83 (2018).Article 

Google Scholar 
Wang, Z. et al. Soil respiration response to alterations in precipitation and nitrogen addition in a desert steppe in northern China. Sci. Total Environ. 688, 231–242 (2019).ADS 
CAS 
Article 

Google Scholar 
Luo, J. et al. Temporal-spatial variation and controls of soil respiration in different primary succession stages on glacier forehead in Gongga Mountain China. PLoS ONE 7, 1–9 (2012).Article 

Google Scholar 
Tong, X. et al. Ecosystem carbon exchange over a warm-temperate mixed plantation in the lithoid hilly area of the North China. Atmos Environ. 49, 257–267 (2012).ADS 
CAS 
Article 

Google Scholar 
Hursh, A. et al. The sensitivity of soil respiration to soil temperature, moisture, and carbon supply at the global scale. Glob. Change Biol. 23, 2090–2103 (2017).ADS 
Article 

Google Scholar 
Huang, S. D. et al. Autotrophic and heterotrophic soil respiration responds asymmetrically to drought in a subtropical forest in the southeast China. Soil Biol. Biochem. 123, 242–249 (2018).CAS 
Article 

Google Scholar 
Zeng, X., Song, Y., Zhang, W. & He, S. Spatio-temporal variation of soil respiration and its driving factors in semi-arid regions of north China. Chin. Geogr. Sci. 28, 12–24 (2018).Article 

Google Scholar 
Li, X. et al. Contribution of root respiration to total soil respiration in a semi-arid grassland on the Loess Plateau China. Sci Total Environ. 627, 1209–1217 (2018).ADS 
CAS 
Article 

Google Scholar 
Luo, Y. & Zhou, X. Soil Respiration and the Environment. 3–4, (Elsevier, 2006).Bhupinderpal, S. et al. Tree root and soil heterotrophic respiration as revealed by girdling of boreal Scots pine forest: extending observations beyond the first year. Plant Cell Environ. 26, 1287–1296 (2003).Article 

Google Scholar 
Lavigne, M. et al. Soil respiration responses to temperature are controlled more by roots than by decomposition in balsam fir ecosystems. Can J Forest Res. 33, 1744–1753 (2003).CAS 
Article 

Google Scholar 
Rey, A. et al. Annual variation in soil respiration and its components in a coppice oak forest in Central Italy. Glob. Change Biol. 8, 851–866 (2002).ADS 
Article 

Google Scholar 
Hartley, I., Heinemeyer, A., Evans, S. & Ineson, P. The effect of soil warming on bulk soil vs rhizosphere respiration. Glob. Change Biol. 13, 2654–2667 (2007).ADS 
Article 

Google Scholar 
Zheng, Y., Zhang, Z., Hu, Y., Yao, D. & Chen, X. Seasonal variation of soil respiration and its environmental effect factors on refactoring soil in coal mine reclamation area. J. China Coal Soc. 39, 2300–2306 (2014).
Google Scholar 
Ren, Z. et al. Effect of weathered coal on soil respiration of reconstructed soils on mining area’s earth disposal sites in Shanxi-Shaanxi-Inner Monglia adjacent area. Trans. CSAE 31, 230–237 (2015).
Google Scholar 
Wang, F. Effect of coversoil thickness on reconstruction soil respiration characteristics in coal mining areas-A case from Panji mining area in Huainan China. Huainan Anhui Univ. Sci. Technol. 1, 59–60 (2017).
Google Scholar 
Sun, Z. H., Han, J. C. & Wang, H. Y. Soft rock for improving crop yield in sandy soil of Mu Us sandy land China. Arid Land Res Manag. 33, 136–154 (2019).CAS 
Article 

Google Scholar 
Sun, Z. H. & Han, J. C. Effect of soft rock amendment on soil hydraulic parameters and crop performance in Mu Us sandy land China. Field Crop Res. 222, 85–93 (2018).Article 

Google Scholar 
Liu, Y. S., Yang, Y. Y., Li, Y. Y. & Li, J. T. Conversion from rural settlements and arable land under rapid urbanization in Beijing during 1985–2010. J. Rural Stud. 51, 141–150 (2017).Article 

Google Scholar 
Lei, N. & Han, J. C. Effect of precipitation on soil respiration of different reconstructed soils. SCI REP-UK 10, 7328 (2020).ADS 
CAS 
Article 

Google Scholar 
Jin, Z., Qi, Y., Yun, S. & Domroes, M. Seasonal patterns of soil respiration in three types of communities along grass-desert shrub transition in Inner Mongolia, China. Adv atmos Sci. 26, 503–512 (2009).CAS 
Article 

Google Scholar 
Wang, X. et al. Soil respiration under climate warming: differential response of heterotrophic and autotrophic respiration. Glob. Change Biol. 20, 3229–3237 (2014).ADS 
Article 

Google Scholar 
Zhao, C., Zhao, Z., Hong, Z. & Jun, L. Contribution of root and rhizosphere respiration of Haloxylon ammodendron to seasonal variation of soil respiration in the Central Asian desert. Quatern Int. 244, 304–309 (2011).Article 

Google Scholar 
Hanson, P., Edwards, N., Garten, C. & Andrews, J. Separating root and soil microbial contributions to soil respiration: a review of methods and observations. Biogeochemistry 48, 115–146 (2000).CAS 
Article 

Google Scholar 
Liu, H. & Li, F. Effects of shoot excision on in situ soil and root respiration of wheat and soybean under drought stress. Plant Growth Regul. 50, 1–9 (2006).ADS 
CAS 
Article 

Google Scholar 
Han, X., Zhou, G. & Xu, Z. Research and prospects for soil respiration of farmland ecosystems in China. J Plant Ecol. 32, 719–733 (2008).CAS 

Google Scholar 
Tong, D., Xiao, H., Li, Z., Nie, X. & Huang, J. Stand ages adjust fluctuating patterns of soil respiration and decrease temperature sensitivity after revegetation. Soil Sci. Soc. Am. J. 84, 760–774 (2020).ADS 
CAS 
Article 

Google Scholar 
Gromova, M., Matvienko, A., Makarov, M., Cheng, C. & Menyailo, O. Temperature Sensitivity (Q10) of soil basal respiration as a function of available carbon substrate, temperature, and moisture. Eurasian Soil ence. 53, 377–382 (2020).ADS 
CAS 
Article 

Google Scholar 
Meyer, N., Welp, G. & Amelung, W. The temperature sensitivity (Q10) of soil respiration: controlling factors and spatial prediction at regional scale based on environmental soil classes. Glob. Biogeochem. Cy. 32, 204–210 (2018).Article 

Google Scholar 
Tang, X., Shao, H. & Liang, H. Soil respiration and net ecosystem production in relation to intensive management in moso bamboo forests. CATENA 137, 219–228 (2016).Article 

Google Scholar 
Zhou, Y., Wang, F., Chen, X., Chen, M. & Liu, B. Effects of ecological restoration patterns on diurnal variation of CO2 flux from rehabilitated soil of coal mining areas in Huainan City. Bull. Soil Water Conserv. 36, 40–46 (2016).
Google Scholar 
Lellei, K. et al. Temperature dependence of soil respiration modulated by thresholds in soil water availability across European shrub land ecosystems. Ecosystems 19, 1460–1477 (2016).Article 

Google Scholar 
Zhan, X., Yu, G., Zheng, Z. & Wang, Q. Carbon emission andspatial pattern of soil respiration of terrestrial ecosystems in China: based on geostatistic estimation of flux measurement. Adv. Earth Sci. 31, 97–108 (2012).
Google Scholar  More

Karr, J.R., & Chu, E.W. Introduction: sustaining living rivers. In Assessing the Ecological Integrity of Running Waters, Developments in Hydrobiology, vol 149 (eds. Jungwirth, M., Muhar, S., & S. Schmutz, S.) 1–14. (Springer: Dordrecht, 2000).Lu, S., Dai, W., Tang, Y. & Guo, M. A review of the impact of hydropower reservoirs on global climate change. Sci. Total Environ. 711, 134996 (2020).ADS 
CAS 
Article 

Google Scholar 
Liu, C., Ahn, C. R., An, X. & Lee, S. H. Life-cycle assessment of concrete dam construction: comparison of environmental impact of rock-filled and conventional concrete. J. Constr. Eng. Manage. 20139(12), A4013009. https://doi.org/10.1061/(ASCE)CO.1943-7862.0000752 (2013).Article 

Google Scholar 
Maavara, T. et al. River dam impacts on biogeochemical cycling. Nat. Rev. Earth Environ. 1, 103–116 (2020).ADS 
Article 

Google Scholar 
Grigg, N. S. Global water infrastructure: state of the art review. Int. J. Water Resour. Dev. 35(2), 181–205. https://doi.org/10.1080/07900627.2017.1401919 (2019).Article 

Google Scholar 
European Environment Agency. European waters: Assessment of status and pressures 2018. https://www.eea.europa.eu/publications/state-of-water (Publications Office of the European Union (2018).Belletti, B. et al. More than one million barriers fragment Europe’s rivers. Nature 588, 436–441 (2020).ADS 
CAS 
Article 

Google Scholar 
Grill, G. et al. An index-based framework for assessing patterns and trends in river fragmentation and flow regulation by global dams at multiple scales. Environ. Res. Lett. 10(1), 015001 (2015).ADS 
Article 

Google Scholar 
Kim, J. & An, K. G. Integrated ecological river health assessments, based on water chemistry, physical habitat quality and biological integrity. Water 7(11), 6378–6403. https://doi.org/10.3390/w7116378 (2015).ADS 
CAS 
Article 

Google Scholar 
Vörösmarty, C. J. et al. Global threats to human water security and river biodiversity. Nature 467, 555–561. https://doi.org/10.1038/nature09440 (2010).ADS 
CAS 
Article 

Google Scholar 
McCartney, M. Living with dams: managing the environmental impacts. Water Policy 11(S1), 121–139 (2009).MathSciNet 
Article 

Google Scholar 
Van Cappellen, P. & Maavara, T. Rivers in the Anthropocene: global scale modifications of riverine nutrient fluxes by damming. Ecohydrol. Hydrobiol. 16(2), 106–111 (2016).Article 

Google Scholar 
Drouineau, H. et al. Freshwater eels: a symbol of the effects of global change. Fish Fish 19(5), 903–930 (2018).Article 

Google Scholar 
Jones, J. et al. A comprehensive assessment of stream fragmentation in Great Britain. Sci. Total Environ. 673, 756–762. https://doi.org/10.1016/j.scitotenv.2019.04.125 (2019).ADS 
CAS 
Article 

Google Scholar 
Reid, A. J. et al. Emerging threats and persistent conservation challenges for freshwater biodiversity. Biol. Rev. 94, 849–873 (2019).Article 

Google Scholar 
Hermoso, V., Clavero, M., Blanco-Garrido, F. & Prenda, J. Invasive species and habitat degradation in Iberian streams: an analysis of their role in freshwater fish diversity loss. Ecol. Appl. 21(1), 175–188 (2011).Article 

Google Scholar 
Maceda-Veiga, A. Towards the conservation of freshwater fish: Iberian Rivers as an example of threats and management practices. Rev. Fish Biol. Fish. 23(1), 1–22 (2013).Article 

Google Scholar 
Sánchez-Pérez, A. et al. Seasonal use of fish passes in a modified Mediterranean river: first insights of the LIFE+ Segura-Riverlink. FiSHMED 008, 3. https://doi.org/10.29094/FiSHMED.2016.008 (2016).Article 

Google Scholar 
Schiermeir, Q. Dam removal restores rivers. Nature 557, 290–291. https://doi.org/10.1038/d41586-018-05182-1 (2018).ADS 
CAS 
Article 

Google Scholar 
Benjankar, R. et al. Dam operations may improve aquatic habitat and offset negative effects of climate change. J. Environ. Manage. 213, 126–134. https://doi.org/10.1016/j.jenvman.2018.02.066 (2018).Article 

Google Scholar 
Tupiño Salinas, C. E., Pinto Vidal de Oliveira, V., Brito, L., Ferreira, A. V. & de Araújo, J. C. Social impacts of a large-dam construction: the case of Castanhão, Brazil. Water Int. 44(8), 871–885. https://doi.org/10.1080/02508060.2019.1677303 (2019).Article 

Google Scholar 
Opperman, J. J. et al. Valuing Rivers: How the diverse benefits of healthy rivers underpin economies. WWF Global Science (2018).Kellner, E. Social acceptance of a multi-purpose reservoir in a recently deglaciated landscape in the Swiss Alps. Sustainability 11, 3819. https://doi.org/10.3390/su11143819 (2019).Article 

Google Scholar 
Boyé, H., & de Vivo, M. The environmental and social acceptability of dams. Field Actions Sci. Rep. http://journals.openedition.org/factsreports/4055 (2016).Wiejaczka, Ł, Piróg, D. & Fidelus-Orzechowska, J. Cost-benefit analysis of dam projects: the perspectives of resettled and non-resettled communities. Water Resour. Manag. 34(1), 343–357 (2020).Article 

Google Scholar 
Rodeles, A. A., Galicia, D. & Miranda, R. Recommendations for monitoring freshwater fishes in river restoration plans: a wasted opportunity for assessing impact. Aquat. Conserv. 27(4), 880–885. https://doi.org/10.1002/aqc.2753 (2017).Article 

Google Scholar 
Birnie-Gauvin, K., Tummers, J. S., Lucas, M. C. & Aarestrup, K. Adaptive management in the context of barriers in European freshwater ecosystems. J. Environ. Manag. 204, 436–441. https://doi.org/10.1016/j.jenvman.2017.09.023 (2017).Article 

Google Scholar 
Yousefi-Sahzabi, A. et al. Turkish challenges for low-carbon society: current status, government policies and social acceptance. Renew. Sustain. Energy Rev. 68, 596–608. https://doi.org/10.1016/j.rser.2016.09.090 (2017).Article 

Google Scholar 
Jiang, H., Lin, P. & Qiang, M. Public-opinion sentiment analysis for large hydro projects. J. Construct. Eng. Manage. 142(2), 05015013. https://doi.org/10.1061/(ASCE)CO.1943-7862.0001039 (2016).Article 

Google Scholar 
Schulz, C., Martin-Ortega, J. & Glenk, K. Understanding public views on a dam construction boom: the role of values. Water Resour. Manage. 33, 4687–4700. https://doi.org/10.1007/s11269-019-02383-9 (2019).Article 

Google Scholar 
Kirchherr, J., Pohlner, H. & Charles, K. J. Cleaning up the big muddy: A meta-synthesis of the research on the social impact of dams. Environ. Impact Assess. Rev. 60, 115–125. https://doi.org/10.1016/j.eiar.2016.02.007 (2016).Article 

Google Scholar 
Piróg, D., Fidelus-Orzechowska, J., Wiejaczka, L. & Łajczak, A. Hierarchy of factors affecting the social perception of dam reservoirs. Environ. Impact Assess. Rev. 79, 106301. https://doi.org/10.1016/j.eiar.2019.106301 (2019).Article 

Google Scholar 
Arboleya, E., Fernandez, S., Clusa, L., Dopico, E. & Garcia-Vazquez, E. River connectivity is crucial for safeguarding biodiversity but may be socially overlooked. Insights from Spanish University students. Front. Environ. Sci. 9, 643820. https://doi.org/10.3389/fenvs.2021.643820 (2021).Article 

Google Scholar 
Gilg, A., & Barr, S. Behavioural attitudes towards water saving? Evidence from a study of environmental actions. Ecol. Econ. 57(3), 400–414. doi:https://doi.org/10.1016/j.ecolecon.2005.04.010 (2006)Schapper, A., Unrau, C., & Killoh, S. Social mobilization against large hydroelectric dams: a comparison of Ethiopia, Brazil, and Panama. Sustain. Develop. 28, 413–423. doi:https://doi.org/10.1002/sd.1995 (2020)Flaminio, S., Piégay, H., & Le Lay, Y-F. To dam or not to dam in an age of anthropocene: insights from a genealogy of media discourses. Anthropocene. 36, 100312, doi:https://doi.org/10.1016/j.ancene.2021.100312 (2021)Bellmore, J. R. et al. Conceptualizing ecological responses to dam removal: If you remove it, what’s to come?. Bioscience 69(1), 26–39. https://doi.org/10.1093/biosci/biy152 (2019).Article 

Google Scholar 
Heberlein, T. A. Navigating environmental attitudes. Conserv. Biol. 26(4), 583–585. https://doi.org/10.1111/j.1523-1739.2012.01892.x (2012).Article 

Google Scholar 
Lewandowsky, S., Gignac, G. E. & Vaughan, S. The pivotal role of perceived scientific consensus in acceptance of science. Nat. Clim. Change. 3, 399–404. https://doi.org/10.1038/NCLIMATE1720 (2013).ADS 
Article 

Google Scholar 
Schuldt, J. P., Roh, S. & Schwarz, N. Questionnaire design effects in climate change surveys: Implications for the partisan divide. Ann. Am. Acad. Pol. Soc. Sci. 658(1), 67–85. https://doi.org/10.1177/0002716214555066 (2015).Article 

Google Scholar 
Bowden, V., Nyberg, D. & Wright, C. Planning for the past: local temporality and the construction of denial in climate change adaptation. Glob. Environ. Change 57, 101939. https://doi.org/10.1016/j.gloenvcha.2019.101939 (2019).Article 

Google Scholar 
Venus, T. E., Hinzmann, M., Bakken, T. H., Gerdes, H., Nunes Godinho, F., Hansen, B., Pinheiro, A., & Sauer, J. The public’s perception of run-of-the-river hydropower across Europe. Energy Policy. 140, 111422. doi:https://doi.org/10.1016/j.enpol.2020.111422 (2020)Schober, M. F. The future of face-to-face interviewing. Qual. Assur. Educ. 26(2), 290–302. https://doi.org/10.1108/QAE-06-2017-0033 (2018).MathSciNet 
Article 

Google Scholar 
Couper, M. P. The future of modes of data collection. Public Opin. Q. 75, 889–908. https://doi.org/10.1093/poq/nfr046 (2011).Article 

Google Scholar 
Zhang, X., Kuchinke, L., Woud, M. L., Velten, J. & Margraf, J. Survey method matters: Online/offline questionnaires and face-to-face or telephone interviews differ. Comput. Hum. Behav. 71, 172–180. https://doi.org/10.1016/j.chb.2017.02.006 (2017).Article 

Google Scholar 
Garcia de Leaniz, C., Berkhuysen, A., & Belletti, B. Beware small dams, they can do damage, too. Nature 570, 164–164; doi:https://doi.org/10.1038/d41586-019-01826-y (2019).Belletti, B., et al. Small isn’t beautiful: the impact of small barriers on longitudinal connectivity of European rivers. Geophys. Res. Abst. 20: EGU2018-PREVIEW (2018).Hophmayer-Tokich, S. & Krozer, Y. Public participation in rural area water management: experiences from the North Sea countries in Europe. Water Int. 33(2), 243–257. https://doi.org/10.1080/02508060802027604 (2008).Article 

Google Scholar 
San-Martín, E., Larraz, B. & Gallego, M. S. When the river does not naturally flow: a case study of unsustainable management in the Tagus River (Spain). Water Int. 45(3), 189–221. https://doi.org/10.1080/02508060.2020.1753395 (2020).Article 

Google Scholar 
Dunlap, R. E. Environmental concern. The Wiley‐Blackwell Encyclopedia of Globalization. (Wiley, Amsterdam, 2012).European Commission Ethics for researchers. Facilitating Research Excellence in FP7. https://doi.org/10.2777/7491 (Publications Office of the European Union, 2013).Jenner, B. M. & Myers, K. C. Intimacy, rapport, and exceptional disclosure: a comparison of in-person and mediated interview contexts. Int. J. Soc. Res. Methodol. 22(2), 165–177. https://doi.org/10.1080/13645579.2018.1512694 (2019).Article 

Google Scholar 
Given, L. M. 100 questions (and answers) about qualitative research (Sage, 2015).
Google Scholar 
Saris, W. E. & Gallhofer, I. N. Design, evaluation, and analysis of questionnaires for survey research (Wiley, 2014).Book 

Google Scholar 
Avella, J. R. Delphi panels: research design, procedures, advantages, and challenges. IJDS 11(1), 305–321. https://doi.org/10.28945/3561 (2016).Article 

Google Scholar 
Vandenplas, C. & Loosveldt, G. Modeling the weekly data collection efficiency of face-to-face surveys: six rounds of the European social survey. J. Surv. Stat. Methodol. 5(2), 212–232. https://doi.org/10.1093/jssam/smw034 (2017).Article 

Google Scholar 
Barbero-García, M. I., Vila-Abad, E. & Holgado-Tello, F. P. Tests adaptation in cross-cultural comparative studies. Acción Psicol. 5, 7–16. https://doi.org/10.5944/ap.5.2.454 (2008).Article 

Google Scholar 
Flick, U. Triangulation in data collection. The SAGE Handbook of Qualitative Data Collection. (Sage, London, 2018).Heesen, R., Bright, L. K. & Zucker, A. Vindicating methodological triangulation. Synthese 196(8), 3067–3081. https://doi.org/10.1007/s11229-016-1294-7 (2019).MathSciNet 
Article 

Google Scholar 
DeVellis, R. F. Scale development: Theory and applications (Sage, 2012).
Google Scholar 
Hammer, Ø., Harper, D.A.T., & Ryan, P.D. PAST: paleontological statistics software package for education and data analysis. Palaeontol. Elect. 4(1), 9. http://palaeo-electronica.org/2001_1/past/issue1_01.htm (2001). More

Deterministic logistic modelThe following population dynamics model was applied to reconstruct the initial dugong population size in 1894 from fishery statistics between 1894 and 1914:$$N_{t + 1} = N_{t} left( {1 , + r{-}r , N_{t} /K} right) – C_{t} ,$$where r is the intrinsic rate of population increase, Nt is the population size in year t, K is the carrying capacity, and Ct is the number of individuals removed from the waters near the Ryukyu Islands in year t. The carrying capacity (K) in 1893 was sufficient to sustain the initial population of dugongs at that time (N1894). The intrinsic rate of population increase (r) was given between 1 and 5% within a range of natural one.Approximate Bayesian calculationWe conducted approximate Bayesian calculation (ABC)32 to estimate the number of individuals in 1979 based on bycatch data between 1979 and 2019, and the constraints of the numbers of individuals were 11 in 1997, three in 2007, and almost extinct in 2019. We denoted fecundity as f, the survival rate until 1 year old as s0, the annual survival rate after 1 year old as s, the age at maturity as am, and the physiological longevity as A. We assumed that the sex ratio at birth was 1:1 on average; the age at maturity am was eight years of age33, and the physiological longevity A was 73 years6. We ignored environmental stochasticity because no mass deaths caused by infectious diseases or changes in survival or mortality rates due to environmental fluctuations have not been recorded during this period. We also ignored density effects because the carrying capacity of the location was sufficiently greater than the initial population size, and our goal was to investigate the possibility of population recovery after a decrease in population using a population dynamics model and estimate the natural growth rate during this period. The detailed extinction risk depends on age structure.According to the life history parameters, except the physiological longevity compiled by (ref.33), the annual survival probability of an a year-old individual is s for a = 1, 2, …, 72; s0 for a = 0, and 0 for a = 73; the reproductive probability of an adult female  > 8 years old is 2f. As the number of years for a population to become extinct or recover depends on age composition, age-specific survival, and reproductive rates, we obtain the population growth rate by the maximum eigenvalue of the following Leslie matrix, L = {Lij} (i = 1,…73, j = 1,…,73) as:$$L_{i1} = s_{0} f/2quad {text{for}}quad i ge a_{m} ,L_{i+ 1,i} = squad {text{for}}quad i = 1, ldots ,72,quad {text{and}}quad L_{ij} = 0,{text{otherwise}}{.}$$We used the population growth rate λ, defined by the maximum eigenvalue of L, as an indicator of the population growth rate.We assumed that the sex of each individual in 1979 was randomly sampled by the 1:1 sex ratio, and its age was randomly sampled by the stable age structure that is given by the eigenvector of the Leslie matrix with the maximum eigenvalue. We assumed that the number of individuals at age 1 year in year t + 1, denoted by N1,t+1, is determined by the binomial distribution:$$Prleft[ {N_{1,t + 1} = x} right] = left( {begin{array}{*{20}c} {N_{f} } \ x \ end{array} } right)left( {s_{0} f} right)^{x} left[ {1 – left( {s_{0} f} right)} right]^{{N_{f} – x}} ,$$where Nf represents the number of adult females in year t. We assumed that no twins were born. We assumed that the probability that an individual with age x survived in the next year is s if x = 1 or s0 if x = 0. We also assumed that Ct individuals who died by bycatch were randomly chosen from any sex and age because the age of individuals caught by bycatch is rarely known. We do not know the sex of some individuals.We assumed the following prior distributions for N1997, f, and s: N1979 (in) U(11, 80), f (in) U(1/14, 1/6) if at least one adult male existed in the population, s0 (in) U(0.1, 0.85); and s (in) U(0.8, 0.97), where U(a, b) is the uniform random variable between a and b. These probabilities were constant for each simulation trial from 1997 to 2019. We selected the set of parameters with the population growth rate (λ) obtained when the maximum eigenvalue of the Leslie matrix was between 0.96 and 1.01.We rejected trials that did not satisfy the following summary statistics: N1997 ≥ 11 (intensive survey in 1997), Nt ≥ 3 during 2004–2017 (monitoring), and N2019 ≤ 1 (“local extinction”). We obtained the prior distributions of N1997, f, s0, s, and N2004, and of the  > 130,000 trials in the prior distribution with natural population growth rates λ of 96.1–98.8%, 99.3% were rejected. For 95% of the 1000 adopted trials, N1979 ranged from 14 to 58. If λ  > 98%, N1997 was ≤ 45 for the adopted trials (Extended Data Fig. 7. Even if all the stranding deaths were due to anthropogenic factors, such as the release of dugongs after bycatch or boat strike, the range of N1997 changed to  98%, with only a slight upward shift, but positive natural growth rate (or λ  > 1) was again very unlikely (0.3%) among the adopted trials.Population viability analysis to assess the impact of bycatch on the extinction riskWe re-evaluated the extinction risk with and without bycatch using the 1000 parameter sets of N1979, f, s0, and s that satisfied the summary statistics in the ABC and stochastic individual-based model, beginning from N1979 for the corresponding parameters. For each parameter set, 100 trials were conducted for each scenario to compare the extinction risks. More

Animal culture and mediumFive different clonal lineages of the water flea Daphnia magna were sampled from two ponds on agricultural land in Belgium (Vleteren: 50°55′06.7″ N, 2°43′27.0″ E and De Haan 51°13′53.8″ N, 3°01′49.2″). They were cultured separately in 210 ml glass jars under optimized laboratory conditions (20 ± 1 °C, 14:10 h light:dark cycle). Seed shrimp and rotifer resting eggs were obtained from a commercial supplier (MicroBioTests Inc., H. incongruens strain MBT/1999/10, product code TB36; B. calyciflorus, product code TK21, Belgium) and represent laboratory cultured, single clonal lineages. More details on animal culture are reported in the online supplementary methods (Appendix 3).Natural pond water was used as medium both in animal cultures and the experiment. It was extracted from a Belgian region (50°59′00.92″ N, 5°19′55.85″ E, Zonhoven) with soft, poorly buffered water (Alkalinity 3–8°d; pH 6.5–8.5) which is likely to be susceptible to acidification under elevated pCO2. More information on medium and mineral composition is reported in the online supplementary information (Appendix 3; Table S3, Appendix 1).Experimental set-upOrganisms were exposed to three pCO2 treatments, an ambient control (C; 1,520 ppm ± 702 SD), an elevated (T1; 25,609 ppm ± 4,541 SD) and an extreme pCO2 level (T2; 83,201 ppm ± 15,533 SD). The control pCO2 level represents the current global mean that is measured in lentic freshwaters considering most ponds and lakes are already supersaturated10,12. The T1 level is currently only observed in more extreme cases11. However, it reflects a pCO2 level that could be encountered more commonly in the field in the future. The T2 treatment represents an extreme test of the tolerance limits of extant species. These treatments are a necessary simplification of reality since pCO2 can experience strong fluctuations in ponds and lakes. An overview of freshwater pCO2 concentrations from literature can be found in Table S1 (Appendix 1).The elevated pCO2 concentrations were manipulated in the water by injecting pure CO2 (99.998% pure, ALPHAGAZ CO2 SFC * B50-N48, Airliquide, Belgium) from gas cylinders into the water (cf.49) at a constant flowrate, using a high-pressure regulator (HBS 200–10.2,5; AirLiquide, Belgium) and a flow controller (Sho-rate model 1350G, Brooks Instruments, USA). In the control treatment, ambient air was supplied at a similar rate as the CO2 to ensure equal perturbation levels across all containers. Water of all experimental containers (including control) were also injected with ambient air to keep the water oxygenated. A relatively constant pCO2 was ensured by continuously monitoring pH and kept between a range of ~ 20,000–30,000 ppm (pH 6.9–6.7) for T1 and ~ 70,000–120,000 ppm (pH 6.4–6.1) for T2 (Figure S2, Appendix 2).Each treatment included 13 replicate 210 mL glass jars per species, resulting in a total of 117 experimental units. Per replicate, one mature water flea (8–11 days old) was inoculated in a jar containing aerated pond water. The five clonal lineages were distributed evenly over the experimental conditions so that each condition had the same number of replicates per clone. Seed shrimp replicates each contained one newly hatched ( More

From the popular algorithms, we chose the Random forest model as the most suitable for our case. The data required for predictions can be divided into plant occurrence records and environmental features. Bioclimatic variables and soil properties were selected as the main environmental features. All of the data were obtained from open sources.Heracleum Sosnowskiy plant descriptionHeracleum sosnowskyi is a monocarpic perennial plant of the Apiaceae family. The height is up to 3–5 m with a straight stem up to 12 cm in diameter. HS compound steam leaves can reach 150 cm, both long and wide38. The blooming period starts in July and continues until the end of September. Plant reproduction is performed by seeds only. The seeds’ depth of germination is reported as mainly in the upper 5 cm down to 15 cm of soil. One plant can produce 10–20,000 seeds39,40. Seeds germinate in the early spring, while some have reported that a period of cold stratification for the dormancy break is obligatory for germination development. Suitable conditions for HS include a temperate climate with warm humid summers and cold winters, while it is probably not drought resistant. Plants of HS tend to neutral soils with a pH range from 6 to 7, rich in nutrients, and being reported as nitrophilous, so the eutrophication of the environment favours HS development. HS plants do not tolerate shade conditions in the first growing period.HS is mostly spread in artificial and semi-natural habitats, including grasslands, pastures, parks, roadsides, agricultural fields, riverbanks or canal sides, and other distributed habitats. Currently, the main pathways of spread include an involuntary entry with soil on vehicles, machinery, footwear or the use of soil as a commodity (as the growing medium rich in organic matter)39.Study areaThe area for modelling extends from approximately 41(^{circ }) to 70(^{circ }) N and from 27(^{circ }) to 60(^{circ }) E, and Kaliningrad region, it equals to approximately 4 mln km2 (Fig. 4).Figure 4Map of the study area: white colour represents the territory used for prediction, red points correspond to the dataset of HS occurrence, collected from the available sources.Full size imageThe European part of Russia is the most inhabited part of the country, and it is the home of approximately 80% of the total population of Russia. It includes the East European Plain, Caucasus mountains and Ural mountains, with the predominance of the East European Plain. Environmental characteristics across the territory of study vary significantly. The climate is changing from semi-arid in the south to subarctic in the north, including humid continental climate conditions. Natural vegetation is represented by almost all types of biomes with the prevalence of different types of forests: broadleaf and mixed forests, coniferous forests, and boreal forests (taiga), while the area of arable lands is reported to be approximately 650,000 km241,42. The territory is subjected to the constant land-use types and cover changes due to the urbanization and switch of the status of arable lands—i.e. reduction of croplands and development of fallows and forests, and, vice versa, returning of some of them into the cultivation process43. The soil cover is represented by the contrast by their physicochemical properties groups, in the northern part of Luvisols, Podzols, Histosols, while of the southern part—by Chernozems, Kastanozems, Solonetz44.Collection of the input dataPlant occurrence dataPlant occurrence coordinates were collected from several publicly available sources related to citizen science projects: the Global Biodiversity Information Facility database45, iNaturalist database46, and the database of the “Antiborschevik” community47. Records were documented by human observation and collected from 2000 to 2021. The overall number of initial occurrence points from combined sources is 7637.Environmental predictorsClimate data Modelling was performed for current and future climate conditions at its two scenarios, selected year ranges were 2000–2018 and 2040–2060 respectively.Climatic variables were collected from the Worldclim database48, containing the average seasonal information relevant to the physiological characteristics of species and available at different resolutions. We chose 10 arc-minutes spatial resolution taking into account the size of the studied area. Table 1 provides a short description of the used bioclimatic features, and we refer the reader to the Worldclim project for detailed information on the variables’ calculation.For the future climate scenarios, we used two Shared Socioeconomic Pathways (SSPs)49—1-2.6 and 5-8.5, corresponding to the lowest (keeping global mean temperature increase below 2 (^{circ })C) and the highest (at the increase of population without technological change) predicted future greenhouse gases emission scenarios. For these data, we took the same resolution (10 arc-minutes) as discussed above.We used the Equilibrium Climate Sensitivity to select the climate model to model future HS distribution. Equilibrium climate sensitivity (ECS) is defined as the global mean surface air temperature change due to a rapid doubling of carbon dioxide concentrations as soon as the associated ocean-atmosphere-sea ice system reaches equilibrium. As the ECS value increases, the model’s sensitivity to the CO(_2) concentration in the atmosphere increases. We have chosen CanESM5 model (ECS—5.6), CNRM-CM6-1 model (ECS—4.3) and BCC-CSM2-MR model (ECS—3.0)50.Table 1 Description of used bioclimatic variables.Full size tableFor the future climate scenarios we selected three climate models:

BCC-CSM2-MR Beijing Climate Center climate system model developed in Beijing Climate Center, China Meteorological Administration51. Model has horizontal resolution 1.125(^{circ }) by 1.125(^{circ }).

CanESM5 Canadian Earth System Model version 5 developed in Canadian Center for Climate Modelling and Analysis, Canada52. Horizontal resolution 2.81(^{circ }) by 2.81(^{circ }).

CNRM-CM6-1 Climate model developed in National Center of Meteorological Research, France53. Horizontal resolution 1.4(^{circ }) by 1.4(^{circ }).

Authors of the WorldClim project prepared historical and future climate data to a uniform spatial (10 arc-minutes) and temporal resolution.Soil data Soil data were downloaded from the SoilGrids database54—a system for global digital soil mapping. SoilGrids provides continuous data at several depths of the spatial distribution of soil properties across the globe with selected resolution. It uses a machine learning approach to reconstruct continuous data from 230,000 soil profile observations from the WoSIS (The World Soil Information Service) database and a series of environmental covariates.From the whole set of the data provided by SoilGrids several properties were chosen for the forecasting: relative percentage of silt (Silt, %), sand (Sand, %), a volumetric fraction of coarse fragments (CF, %), cation exchange capacity (CEC, ({text{cmol}}_{c}/{text{kg}})) and soil organic carbon (SOC, g/kg) at the depth 5–15 cm, where the HS seeds are assumed to be located. These variables are expected to be more stable over time than bioclimatic predictors; thus, chosen soil properties could be implemented for the future time the same as in the present.Data pre-processingAll the data were transformed to the ASCII format by R script and using software DIVA-GIS following the tutorial for the preparation of WorldClim files for use in SDM (http://www.lep-net.org/wp-content/uploads/2016/08/WorldClim_to_MaxEnt_Tutorial.pdf) with unified selected resolution 340 sq.km.Optimization of the occurrence points amountThe general problem in using the available data collected from the databases of the citizen science projects is that the points of observation are distributed non-uniformly. For instance, the frequency of the records depends on the density of the population directly. The spatial filtering of the data (reducing the number of points) can be performed to reduce the sampling bias55. We prepared three datasets with a distance between points of 4, 7 and 10 km with 2402, 1846 and 1504 occurrence points correspondingly filtering the initial dataset. For the thinning step thin() function was used within the R package spThin with 100 iterations for each of chosen thinning distances. To understand how much data we could lose, we used the analysis of feature distribution and evaluated the general fairness of the model performance.Pseudo-absence generationDue to the availability only of the presence points, it is important to generate the absence points for further implementation of the selected algorithm. Although the generation of pseudo-absence points in SDM research is a widespread solution, a closer look at the literature reveals several gaps and shortcomings. Since the raw dataset of the HS distribution demonstrates strong sampling bias, the generation of pseudo-absence points using the usual ‘random’ strategy can aggravate the sampling bias problem. Thus, the combination of the ‘disk’ and ‘random’ strategies was applied for the generation of the pseudo-absence points using the biomod R package17.

The ‘disk’ strategy is established on the geographic distance works as separation from truth presence and possible absence points. The optimal geographic distance for HS was chosen as 25 km. This distance was chosen empirically by trial-and-error. We started with 18 km (because the size of the cell is   9–18 km depending on location) and finished with 50 km. Using distances such as 30–50 km lead to a positive spatial autocorrelation. Thus, we decided to set 25 km which finally provided both optimal model performance and reduced spatial autocorrelation.

The second part of the generation was based on the ‘random’ strategy with filtration: according to the different range of climate conditions on the territory of Russia, there are several places where HS is not detected, thus not growing. The selection of unsuitable places for HS related to the north of Russia, where it is might be too cold for plant species. From all amount of randomly generated generated points we selected points with condition latitude ( > 64^{circ }), according to tundra board line.

Features selection procedureTo avoid over-fitting and to choose the most conscientious set of parameters for final modelling, two approaches were combined. We searched features that are not correlated with others by a selected threshold is equal to 0.8 in absolute values56 and estimated variable importance using the Mean Decrease Gini (MDG) and the Mean Decrease Accuracy (MDA) as the result of modelling on enumerated parameters’ combinations. MDG score is related to the homogeneity of the nodes and leaves coefficient. With the rise of the MDG score the importance of the corresponding feature is also increasing. MDA describes how much accuracy decrease by removing the feature. We selected the most important features according to the MDG and MDA scores by the highest values of both metrics using a sequential search from an initial set of variables.Modelling approachRandom forestChoosing the appropriate method for creating the tool for accurate SDM is crucial because the overall performance could vary dramatically, depending on the selected model and particular use case. There is a limited amount of acceptable machine learning methods that can be used in SDM. Several popular methods demonstrated high performance in modelling on large areas: GBM, RF, and GLM. In particular, for modelling and prediction of the potential distribution of invasive species, GLM and RF were used57. We decided to use RF because this model was successfully implemented for solving a variety of tasks such as predictions of animal and plant distributions, and also was used for making predictions on a large territory58. The other important advantage that should be noticed is the straightforward interpretability of RF, which means that it is possible to evaluate the impact of each environmental parameter on the occurrence of the invasive species.Approach to the cross-validation of the modelA unique approach for the model calibration is needed to reduce spatial autocorrelation caused by the absence of a strict sampling design. In our case, the data was split into training and testing folds using the spatial blocks technique in a scheme of 13-fold cross-validation. Random spatial splitting was performed 20 times to calibrate the model, with a distance between blocks set as 100 km. To calibrate the model we used a spatial blocks approach with random type from R package blockCV.Evaluation of the model performanceTo evaluate the performance of the model a classic approach for ecology was used—Area Under Curve (AUC) or Receiver operating characteristic (ROC), related to the independent threshold techniques16. The principle of methods lies in the standard confusion matrix, where rows and columns represent actual and predicted classes. The construction of ROC curves uses all possible thresholds to obtain different confusion matrices which leads to the reproduction of the curve with two-dimensional space: (1) on y-axis is True Positive Rate (sensitivity, recall); (2) on x-axis is False Positive Rate (equal to 1 − specificity). In our case true positive (TP, sensitivity) rate means that predicted places where HS grows correspond to actual. Similarly, true negative rate (TN, specificity) indicates correctly classified locations as absence points. In contrast, the missteps when the model predicted places as presence points for plants that are incorrect are False Positive, FP, and places where HS is absent, according to the model, while this is not true are recognised as False Negative, FN. More