Ecology

Ant hostWe used workers of the invasive Argentine ant, Linepithema humile, as host species. As typical for invasive ants, populations of this species lack territorial structuring and instead consist of interconnected nests forming a single supercolony with constant exchange of individuals between nests40. We collected L. humile queens, workers and brood in 2011, 2016 and 2022 from its main supercolony in Europe that extends more than 6,000 km along the coasts of Portugal, Spain and France40,41,42, from a field population close to Sant Feliu de Guíxols, Spain (41° 49’ N, 3° 03’ E). Field-collected ants were reared in large stock colonies in the laboratory. For the experiments, we sampled worker ants from outside the brood chambers and placed them into petri dishes with plastered ground (Alabastergips, Boesner, BAG), subjected to their respective treatments as detailed below. Experiments were carried out in a temperature- and humidity-controlled room at 23 °C, 65% relative humidity and a 12 h day/night light cycle. During experiments, ants were provided with ad libitum access to a sucrose-water solution (100 g l−1) and plaster was watered every 2–3 d to keep humidity high.Collection of this unprotected species from the field was in compliance with international regulations, such as the Convention on Biological Diversity and the Nagoya Protocol on Access and Benefit-Sharing (ABS, permit numbers ABSCH-IRCC-ES-260624-1 ESNC126 and SF0171/22). All experimental work followed European and Austrian law and institutional ethical guidelines.Fungal pathogensAs pathogen, we used the obligate-killing entomopathogenic fungus Metarhizium, whose infectious conidiospores naturally infect ants43,44,45 by penetrating their cuticles, killing them and growing out to produce highly infectious sporulating carcasses23,46. We used a total of six strains of the two species M. robertsii and M. brunneum, all isolated from the soil of the same natural population—an agricultural field at the Research Centre Årslev, Denmark27,47, which makes host co-infections with these sympatric strains in the field likely. As in ref. 24, we used three strains of M. robertsii (R1: KVL 12-36, R2: KVL 12-38, R3: KVL 12-35) and three of M. brunneum (B1: KVL 13-13, B2: KVL 12-37, B3: KVL 13-14), all obtained from the University of Copenhagen, Denmark (B. M. Steinwender, J. Eilenberg and N. V. Meyling).We started our selection experiment by exposing the ants to a mix of the six strains in equal proportions. To this end, each strain was grown separately from monospore cultivates from its long-term storage (43% glycerol (Sigma-Aldrich, G2025) in skimmed milk, −80 °C) on SDA plates (Sabouraud-4% dextrose agar, Sigma-Aldrich, 84088-500G) at 23 °C until sporulation. Conidiospores (abbreviated to ‘spores’) were collected by suspending them in sterile 0.05% Triton X-100 (Sigma-Aldrich, X-100; in milliQ water, autoclaved) and mixed in equal amounts to a total concentration of 1 × 106 spores ml−1. Before mixing, we confirmed that all strains had ≥98% germination.We exposed worker ants individually to the fungal pathogen by dipping them into the spore suspension using clean forceps (Gebrüder Martin; bioform, B32d). Afterwards, each ant was brieftly placed on filter paper (Whatman; VWR, 512-1027) to remove excess liquid before being placed into its experimental Petri dish.Serial passage experimentWe tested for the long-term effect of social immunity on pathogen selection, in which the pathogen was serially cycled through the host in the absence or presence of social immunity while the host population remained constant.Experimental design and procedureAfter exposure to the fungal spore mix, worker ants were either kept alone (individual host treatment, n = 10 replicate lines) or together with two untreated nestmates (social host treatment, n = 10 replicate lines; Fig. 1a). Individual ants could only protect themselves by individual immunity (selfgrooming behaviour and their physiological immune system), while the attended ants experienced both individual and social immunity due to the additional allogrooming by their caregiving nestmates. Thus, comparing the two host conditions revealed the effect of social immunity.As sanitary care by the nestmates reduces the pathogens’ success to kill the exposed individuals, we had to set up more experimental dishes of the social host treatment to obtain equal numbers of sporulating carcasses under both selection treatments, from which we then collected the spores for the next host infection cycle. For the individual treatment, we exposed an average of 23 workers per cycle, while an average of 40 workers per cycle were exposed in the social host treatment. The experiment was run for 10 host passages, that is, 27 weeks. In total, 6,312 workers (2,299 in the individual and 4,013 in the social host treatment) were exposed during the course of the experiment, and 8,026 nestmates were used. To obtain the spore suspensions for the next steps, we then collected and pooled the outgrowing spores of the first 8 carcasses produced per replicate line and cycle (that is, a total of n = 800 carcasses from the individual and n = 800 carcasses from the social host treatment, over the 10 host passages). Dead nestmates were not considered (see below).In detail, at each host cycle, the freshly exposed ants were placed into Petri dishes with plastered, humidified ground (Ø 3.5 cm for the individual and Ø 6 cm for the social host condition; both Bioswisstec AG, 10035 and 10060) in the absence (individual host treatment) or presence (social host treatment) of two untreated nestmates. We checked survival daily for 8 d. Ants that died within 24 h after exposure were excluded from the experiment as their mortality could not yet have resulted from infection, but rather from treatment procedures. Ants dying from days 2 to 8 were checked for internal Metarhizium infections by surface-sterilization (washing the carcass in 70% ethanol (Honeywell; Bartelt, 24194-2.5l; diluted with water) for a few seconds, rinsing it in distilled water, incubating in 3% bleach (Sigma-Aldrich, 1056142500) in sterile 0.05% Triton X-100 for 3 min and rinsing it again three times in water48), followed by incubation in a Petri dish on humidified filter paper at 23 °C until day 13, when they were checked for Metarhizium spore outgrowth. This timeline was chosen as preliminary work showed that the exposed ants die mostly on days 4 to 8 (median day 5, for both individual and social host treatments) after exposure and that sporulation required no longer than 5 d in our experimental conditions, so that a duration of 13 d per cycle also allowed for the later dying ants to complete sporulation. Preliminary work further revealed that in cases where nestmates contracted the disease, they died at a delayed timepoint and with spore outgrowth mostly around the mouthparts. These characteristics were used to distinguish between the directly exposed ants and infected nestmates in the experiment where ants were not colour-marked. The carcasses of sporulating nestmates were excluded from further procedures. An additional control experiment using 120 sham-treated ants showed no Metarhizium outgrowth, so that all Metarhizium outgrowth in our experiment could be attributed to our experimental infections. Carcasses with saprophytic outgrowth were not considered. For each host passage and each replicate line, we collected the spores of the first 8 ants dying after day 1 from their Metarhizium-sporulating carcasses at day 13 in 0.05% Triton X-100, pooled and counted them using an automated cell counter (Cellometer Auto M10, Nexcelom Bioscience). The concentration of each pool was then adjusted to 1 × 106 spores ml−1, and was used directly (that is, in the absence of any intermediate fungal growth step on agar plates) for exposing the ants in the next host infection cycle. The ants of each host passage were thus dipped in the same spore concentration. The remaining spore suspension was frozen at −80 °C in a long-term storage for further analysis.Pathogen diversity and strain compositionWe analysed which strains were present and in which proportion after 5 and 10 passages in each of the 10 individual and 10 social replicate lines. To this end, we first extracted total DNA from the respective spore pools (n = 40), which we analysed (1) quantitatively for the respective representation of M. robertsii vs M. brunneum (using species-specific real-time PCR targeting the PR1-gene sequence; detailed below) and (2) qualitatively for which of the 6 original strains were still present in the pool (using strain-specific microsatellite analysis; detailed below). We used this first estimate of remaining strain diversity and composition of each pool to determine how many spores we had to analyse separately for their strain identity after individualization by FACS sorting and growing them individually as colony forming units (c.f.u.s). This clone-level strain identification was again performed using microsatellite analysis (n = 1,347 individualized clones from the 40 spore mixes, in addition to n = 27 spores from the 6 ancestral strains; detailed below). Such clonal separation was needed since expansion of the spore mix by growth on SDA plates was not representative of the genetic composition of the strains in the pool, due to strong strain–strain growth inhibition when growing in a mix.In detail, we extracted the DNA of the 6 ancestral strains and the 40 spore mixes (10 each for individual and social lines at passages 5 and 10), as well as of 27 individualized clones of the ancestral strains and 1,374 clones from the 40 pools of passages 5 and 10, by centrifuging 100 µl of the spore suspensions in 1.5 ml tubes (Eppendorf, 0030120086) at full speed for 1 min and discarding the supernatant. Nuclease-free water (50 µl) was added and the spores were crushed in a bead mill (Qiagen TissueLyser II, 85300) at 30 Hz for 10 min using acid-washed glass beads (425–600 µm; Sigma-Aldrich, G8772). DNA was extracted using a DNeasy blood and tissue kit (Qiagen, 69506) following the manufacturer’s instructions, using a final elution volume of 50 µl buffer AE.For the quantitative species-level analysis of the pools, we performed quantitative real-time PCR (qPCR) using primers and differently labelled probes24 that we had developed on the basis of the sequence of the PR1 gene49 (forward: 5′ TCGATATTTTCGCTCCTG, reverse 5′-TTGTTAGAGCTGGTTCTGAAG, PR1 probe M. brunneum: 5′-(6-carboxyfluorescein (6FAM))TATTGTACCTACCTCGATAAGCTTAGAGAC(BHQ1), PR1 probe M. robertsii: 5′-(hexachloro-fluorescein (HEX))AGTATTGTACCTCGATAAGCTCGGAGAC(BHQ1)). Reactions were performed in 20 μl volumes using 10 μl iQ Multiplex Powermix (Bio-Rad, 1725849), with 600 nM of each primer (Sigma-Aldrich), 200 nM of each probe (Sigma-Aldrich) and 2 μl of extracted DNA. The amplification programme was initiated with a first step at 95 °C for 3 min, followed by 40 cycles of 10 s at 95 °C and 45 s at 60 °C. Primer efficiency was above 92% for both primer/probe combinations using standard curves of 10-fold dilutions of known input amounts. Data were analysed using Bio-Rad CFX Manager software.For the strain-specific analysis of both the pools and the individualized clones, we used two microsatellite loci, Ma30750 and Ma205451. Microsatellite locus Ma307 (forward: 5′-(6FAM)CATGCTCCGCCTTATTCCTC-3′, reverse: 5′-GGGTGGCGAAGAAGTAGACG-3′) allowed distinction of all strains except two of the M. brunneum strains (B1 and B3), which were distinguished by microsatellite locus Ma2054 (forward: 5′-(6FAM)GCCTGATCCAGACTCCCTCAGT-3′, reverse: 5′-GCTTTCGTACCGAGGGCG-3′). We analysed the microsatellites by E-Gel high-resolution 4% agarose gels (ILife Technologies, G501804) and fragment length analysis (done by Eurofins MWG) using Peak Scanner software 2.For clone individualization, we used flow cytometry to sort single spores out of the 40 spore pools (and the 6 ancestral strains for comparison) on 96-well plates (TPP; Biomedica, TP-92696) containing SDA (100 µl per well). The unstained spore population was detected using the FSC (forward scatter)/SSC (side scatter) in linear mode (70 μm nozzle, FACS ARIA III, BD Biosciences, as exemplified in Supplementary Fig. 1). Purity mode was set to ‘single cell’ and spore clones were obtained by sorting 1 particle event into each well. Sorting and data analysis were performed using Diva 6.2 software. The number of spores that we obtained for microsatellite analysis varied for each replicate, as it was adjusted to the remaining strain diversity estimate that we obtained from the quantitative and qualitative analysis of the pools. In total, we analysed 4–5 clones per ancestral strain (total n = 27) and a median of 5, but up to 101 different clones for the pools (total n = 1,347), as we intensified analysis for the strains that were revealed to be present at low frequency on the basis of previous analysis.Common garden experimentExperimental design and procedureWe then tested whether the successful lines at the end of the experiment (that is, after 10 host passages) differed in their virulence (induced host mortality) and investment into transmission stages (produced spore number) depending on their selection history (individual vs social), when current host social context either reflected the selection history or not. This common garden experiment thus led to 20 matched combinations of selection history and current condition (10 each of the individual lines in current individual host conditions (individual–individual) and the social lines in current social host conditions (social–social)) and 20 non-matched conditions (10 each of the individual lines in current social host conditions (individual–social) and the social lines in current individual host conditions (social–individual)).We obtained the lines for performance of the common garden experiment by the following procedure: (1) for the 16 out of the 20 replicate lines, where a single strain was the sole remaining representative at the end of the experiment (Fig. 1b), we expanded one of the c.f.u.s grown after FACS sorting (see above) by plating on SDA; (2) for the 4 remaining replicates in which two strains had remained (two individual and two social replicate lines), we expanded one c.f.u. of each of the remaining strains and mixed the spores in their representative proportion, as determined above.Virulence and transmissionFor the 10 individual and 10 social lines, we determined the induced host mortality as a measure of virulence and the outgrowing spore number as transmission stage production under their matched and non-matched current host conditions. We exposed the workers as in the selection treatment, kept them either alone or with two untreated nestmates, and monitored their mortality daily for 8 d. Again, ants dying in the first 24 h after treatment and dying nestmates were excluded from the analysis. In total, we obtained survival data of 797 ants (19–20 ants exposed for each of the 10 replicates from each of 4 combinations of selection history and current host condition). Dead ants were treated as above and their outgrowing spores collected by a needle dipped in sterile 0.05% Triton X-100 directly from the carcass, and resuspended in 100 µl of sterile 0.05% Triton X-100. The number of spores per carcass was counted individually using the automated cell counter, as described above (n = 215; median of 5 per replicate). We excluded one outlier carcass(from replicate I5) where we expected a counting error as this single carcass showed approx. 100-fold higher spore count than the other carcasses of this replicate. Exclusion of this outlier did not affect the statistical outcome. The proportion of ants dying per replicate line for each combination of selection history and current host condition and the number of spores produced by all carcasses per replicate were respectively used as measures of virulence and transmission (mean carcass spore load per replicate plotted in Fig. 2).Allogrooming elicitation by the fungal linesWe determined the allogrooming elicited by the individual and the social lines. To this end, we exposed workers as above and observed the allogrooming performed by two untreated nestmates towards the exposed ant. In detail, we performed 3 biological replicates for each of the 20 replicate lines (n = 10 individual and 10 social lines, resulting in a total of 60 videos), where the exposed ant was placed with two untreated nestmates within 10 min after exposure, and filmed with Ueye cameras for 30 min (whereby 4 cameras were used in parallel, each filming 3 replicates simultaneously, and using StreamPix 5 software (NorPix 2009-2001) for analysis). Videos were obtained in a randomized manner and labels did not contain treatment information so that the observer was blind to both the selection history and individual treatment during the behavioural annotations. For each ant, we observed both self- and allogrooming. Start and end times for each grooming event were determined, supported by use of the software BioLogic (Dimitri Missoh, 2010 (https://sourceforge.net/projects/biologic/)).As the ants in our serial passage and common garden experiments were not colour-marked, we also used unmarked ants for this behavioural experiment to keep conditions the same. This was possible as preliminary data with colour-coded nestmates (n = 18 videos) had shown that exposure alters the ant’s behaviour and that of its untreated nestmates in a predictable way that allows reliable classification of the pathogen-exposed individuals from the untreated nestmates; we used the following rules to classify an ant as the exposed individual: (1) the individual spent >5% more time (of the 30 min observation period) selfgrooming than the other individuals; (2) if the difference in selfgrooming time between the individuals was More

No plants were collected or harmed during this study, and all research involving plants followed relevant national, and international guidelines and legislation.Study areaThe study site encloses a wetland area bordering Ramalhete Channel, in the western part of the Ria Formosa lagoon, a mesotidal system located in southern Portugal (Fig. 1). Lunar tides are semi-diurnal, with a mean tidal range of about 2 m that can reach up to 3.5 m during spring tides. Offshore waves have no major propagation inside the lagoon33,34. Water circulation inside the lagoon is mostly driven by tides. The lagoon extends over 55 km along the coast and is connected to the ocean through six tidal inlets35. The three westmost inlets of the system (Ancão, Faro-Olhão, and Armona), which together capture ca. 90% of the total prism, are highly interconnected, with a strong residual circulation from Faro-Olhão Inlet directed towards Ancão and Armona inlets (located in Fig. 1), during both spring and neap tides36. The tidal currents in Ramalhete Channel, connecting the Faro-Olhão and Ancão Inlet, have high tidal asymmetry and shifts in tidal dominance, from flood to ebb. There are no significant fluvial inputs into the lagoon, with a yearly average terrestrial sediment influx of around 2 × 105 m3/yr37, reaching the system through small streams. The main sediment delivery to the system is through the inlets, though there are few studies assessing related fluxes. The net sediment entry through the stabilized Faro-Olhão Inlet is estimated at 1.4 × 105 m3/year38. Recent sedimentation rates in the marsh of the westmost edge of the lagoon were estimated at 1.1 ± 0.1 mm/yr39.The lagoon system is composed of large salt marsh patches, tidal flats and a complex net of natural, and partially dredged tidal channels. The tidal flats (vegetated and non-vegetated) and salt marshes represent more than 2/3 of the total lagoon area. The salt marshes comprise silt and fine sand40, while coarser (sand to shingle) shell-rich sediment, of marine provenance, is found on tidal channels and the lower domain of intertidal flats41. The dominant intertidal species are Spartina maritima and the seagrass Zostera noltei, the latter occupying an estimated area of 1304 ha, which represent 45% of the total intertidal area42.Figure 1Location of the field site in the Ria Formosa lagoon western sector over a satellite image collected in 2019 (South Portugal; upper panel); zoom to monitoring stations S1 to S4 (left lower panel); and field view of the studied site (right lower panel). Map generated with ArcGIS 10.8 (http://www.esri.com) and Adobe Illustrator 2022. Map data: Google Earth 7.3, image Landsat / Copernicus.Full size imageExperimental setup and data analysisAn experimental setup was deployed in the study area to assess dominant local topography, hydrodynamics (water levels and current velocities), Suspended Sediment Concentrations (SSCs), Deposition Rates (DRs), vegetation characteristics, and bed sediment grain size and organic matter content. Measurements were made during a full tide cycle, on a spring tide (tidal range = 3.2 m), and on a neap tide (tidal range = 1.8 m). Sampling was conducted in four wetland stations: S1 and S2 in a vegetated tidal flat comprising Zostera noltei; S3 in the low marsh comprising Spartina maritima; and S4 in the mid-upper marsh with the most abundant species of Sarcocornia perennis and Atriplex portucaloides (see S1 to S4, Fig. 1); the tidal flat is interrupted by a small oblique secondary tidal creek that flows near S2 station.Stations of sediment sampling and equipment deployment along the transect are illustrated in Fig. 2. During neap tide there was no data collection in S4, since the inundation time of the station was very short. The profile elevation was measured using Real Time Kinematic Differential Global Positioning System (RTK-DGPS, Trimble R6; vertical error in the order of few centimetres), and the slope of each habitat within a transect was calculated and expressed in percentage (%). Vegetation at each point was characterized by the canopy height, calculated as the average shoot length.Suspended Sediment Samplers (SSSs) were installed during low tide in the monitored stations using siphon samplers (Fig. 2) and recovered in the next low tide. These samplers consist of 0.5 L bottles with two holes on the cap, one for water intake and the other for air exhaust, according to the method described in13. Each intake tube is adjusted to form a siphon (i.e., inverse U), allowing to control the water level at which intake starts. Siphons were aligned at the same elevation along the transect for spring and neap tides, which means that all SSSs were collecting at the same time within the tidal cycle. During spring tide, in S1 and S2 at the tidal flat, SSSs were sampling at 0.1, 0.9, and 1.2 m from the bed, while at S3 SSSs were sampling at 0.7 and 1.0 m from the bed, and at S4 the SSS was sampling at 0.1 m from the bed (Fig. 2). During neap tide, in S1 and S2, SSSs were sampling at 0.1 and 0.9 m from the bed, while at S3 the SSS was sampling at 0.7 m from the bed.Surficial sediment samples were collected in each habitat to characterize the sediment grain size (d50) and content of organic matter (% OM). Sediment traps were installed in 3 replicates, during low tide, at each sampling point to measure the short-term sediment deposition rate (i.e., deposition over a tidal cycle, following procedures of43). Traps consisted of 3 cm diameter pre-labeled cylindrical tubes (Falcon® tubes, 50 ml). Traps and sediment samples were transported to the laboratory and maintained in a fridge. The sediment content was washed, and both the inorganic and organic weights were determined.The measured inorganic DR (g/m2/hr) was calculated as:$${text{DR}} = {raise0.7exhbox{${{text{W}}_{{{text{DS}}}} }$} !mathord{left/ {vphantom {{{text{W}}_{{{text{DS}}}} } {{text{A}} cdot {text{T}}}}}right.kern-0pt} !lower0.7exhbox{${{text{A}} cdot {text{T}}}$}}$$
(1)
where WDS is the weight of deposited sediment (in grams), A is the area of the sediment trap opening (m2), and T is in hours. Two different tide durations were considered to compute DRs, one assuming T equal to the hydroperiod in each station, and one assuming T equal to the entire tide duration (~ 12.4 h). These measured DRs are hereon mentioned as flood and tide DRs (DRflood and DRtide, respectively). The former is an expression of the actual deposition rate within the flood phase, during the period in which each station is inundated (and therefore active deposition can take place). The latter is the value used to compare with DRs in literature, which typically corresponds to values averaged over multiple tidal cycles (thus accounting for the entire tide duration).Tide levels were measured in the field using pressure sensors (PT, InSitu Inc. Level TROLL; ~ 2 cm from the bed), deployed from S2 towards S4 (Fig. 2). Velocity currents were measured at 20 cm from the bed, using an electromagnetic current meter (EMCM; Infinity Series JFE Advantech Co., Ltd; in S2 to S4; Fig. 2), and raw data (recording interval: 30 s) were filtered using a 10 min moving average for cross-shore and longshore components. To identify tidal asymmetry and assess the related phase dominance, tidal current skewness was calculated through the formula described in44 by which:
$$Sk_{U} = frac{{frac{1}{N – 1}mathop sum nolimits_{t = 1}^{N} left( {U_{t} – overline{U}} right)^{3} }}{{left( {frac{1}{N – 1}mathop sum nolimits_{t = 1}^{N} left( {U_{t} – overline{U}} right)^{2} } right)^{{{raise0.7exhbox{$3$} !mathord{left/ {vphantom {3 2}}right.kern-0pt} !lower0.7exhbox{$2$}}}} }}$$
(2)
where N is the number of recordings, Ut is the input velocity signal and (overline{U}) is the mean velocity. Positive/negative skewness indicates flood/ebb dominance (assuming that flood currents are positive).Figure 2Deployment of the sediment traps, SSSs and devices (electromagnetic current meter EMCM; pressure transducer PT) in the stations (S1 to S4) during spring tide (sketch is exaggerated in the vertical).Full size imageComplementary to the measured DRs, theoretical DRs were also determined from the data, allowing us to link the sediment and flow data collected, and validate the deposition patterns from the traps. The theoretical deposition rate was determined based on45 formula:$${text{DR}} = left{ {begin{array}{*{20}c} {{text{C}}_{{text{b}}} cdot {text{w}}_{{text{s}}} cdot left( {1 – frac{{{uptau }_{{text{b}}} }}{{{uptau }_{{{text{cd}}}} }}} right)} & {{uptau }_{{text{b}}} < {uptau }_{{{text{cd}}}} } \ 0 & {{uptau }_{{text{b}}} ge {uptau }_{{{text{cd}}}} } \ end{array} } right.$$ (3) where Cb is the SSC at the bed, ws is the flock settling velocity, τb is the bed shear stress and τcd is the corresponding critical value for deposition.To determine the settling rate of the flocculates, the modified Stokes’ velocity for cohesive sediment was used, taking shape factors α and β (α = β = 1 for perfectly spherical particles):$${text{w}}_{{text{s}}} = frac{{upalpha }}{{upbeta }} cdot frac{{left( {{uprho }_{{text{s}}} - {uprho }_{{text{w}}} } right) cdot {text{g}} cdot {text{D}}_{50}^{2} }}{{{uprho }_{{text{w}}} cdot 18 cdot {upnu }}}$$ (4) where ρw and ρs are the densities of the water and sediment, respectively and ν is the kinematic viscosity of water (~ 106 m2/s).The bed shear stress τb was calculated from the measured current magnitude, |U| using the law of the wall:$$begin{array}{*{20}c} \ {{uptau }_{{text{b}}} = {uprho }_{{text{w}}} cdot {text{u}}_{*}^{2} , {text{u}}_{*} = frac{left| U right| cdot kappa }{{ln left( {{raise0.7exhbox{$z$} !mathord{left/ {vphantom {z {z_{0} }}}right.kern-0pt} !lower0.7exhbox{${z_{0} }$}}} right)}} } \ end{array} { }$$ (5) where κ is the von Kármán constant (~ 0.4) and z0 is the roughness length. For Zostera noltei, the roughness length was estimated at 5 mm46, value that was also used in the other stations, in lack of related estimate for marsh plants.The critical shear for deposition, τcd, was calculated using the formula47:$$sqrt {frac{{{uptau }_{{{text{cd}}}} }}{{{uprho }_{{text{w}}} }}} = left{ {begin{array}{*{20}c} {0.008} & {{text{w}}_{{text{s}}} le 5 cdot 10^{ - 5} {text{m}}/{text{s}}} \ {0.094 + 0.02 cdot {text{log}}_{10} left( {{text{w}}_{{text{s}}} } right)} & {3 cdot 10^{ - 4} le {text{w}}_{{text{s}}} le 5 cdot 10^{ - 5} {text{m}}/{text{s}}} \ {0.023} & {{text{w}}_{{text{s}}} ge 3 cdot 10^{ - 4} {text{m}}/{text{s}}} \ end{array} } right.$$ (6) Theoretical values of minimum SSCs needed for these DRs were also calculated, assuming that there is constant deposition (i.e., setting τb = 0), and compared with the field results. More

Himalayan marmot habitat analysisThe environmental factors including temperature, vegetation and elevation are the key drivers for the wildlife in alpine ecosystems32. Specific landform attributes such as slope and elevation and vegetation cover affect the population and burrowing of rodents33. For example, rodent burrows in the Western Usambara Mountains in Tanzania were only found at an elevation of above 1600 m33. However, the Himalayan marmot seems to prefer to inhabit areas with low elevation and high land surface temperature34. In this study, the data showed that 76.25% of the Himalayan marmots were found in areas with elevation values of 3400–4600 m. The majority of marmots were found in areas with slopes of 5–20° and vegetation cover higher than 60%. Most marmots were found in alpine meadows, a few were found in temperate grasslands and alpine grasslands, and none were found in other grassland types.Preliminary statistical analysis of vegetation cover, grass type, vegetation type, and Himalayan marmot distribution sample sites obtained using spatial geographic information technology revealed that the meadow grassland areas with lush grass growth, more dominant plants, and abundant food had more marmots. When the vegetation cover reached 0.60–1.00, the number of marmot distribution sample sites was the highest. Dense grass is an ideal habitat and provides concealment for Himalayan marmots, and the abundant plant types provide sufficient food for marmots. In contrast, no marmots were distributed in the alpine scrub, coniferous forest, and alpine snow/ice covered areas where vegetation growth was poor, vegetation cover was low, and food was relatively scarce. Moreover, 70.24% of Himalayan marmots were found in alpine meadows with a wide variety of plant species, including Poaceae, Cyperaceae, and grasses. This finding indicated that alpine meadows are more suitable for Himalayan marmots and have more advantageous habitat conditions compared with other grassland types. The elevation of alpine meadows is 3236–5126 m, and the vegetation is mainly meadows with simple vegetation structure, substantial vegetation cover and dense vegetation growth, and a wide variety of plants, rich food, soft grass, and good palatability. Therefore, alpine meadows provide good natural habitats and foraging sites for marmots.Habitat selection of large rodents is influenced by a combination of vegetation cover availability, food availability, and population density35. Vegetation cover is an important parameter that describes vegetation communities and ecosystems and is closely related to vegetation quantity and productivity. The quality of habitat vegetation is an important factor that affects the spatial distribution of plateau rodents. Both feeding and concealment depend on vegetation, and the height and cover of edible plants and vegetation suitable for concealment determine the choice of vegetation type by marmots. Thus, vegetation cover becomes an important factor for habitat selection by marmots. Different grassland types determine different plant conditions, and selection of different vegetation conditions can increase the chances of survival and improve the reproductive success of marmots; therefore, grassland type is an important ecological factor in habitat selection by marmots. A study showed that the ecological factors affecting habitat selection of Himalayan marmots are mainly topography, anthropogenic disturbance, and vegetation8. Another study concluded that habitat selection by Himalayan marmots is closely related to elements such as topography, landform, temperature, precipitation, and vegetation24.The functions of burrows’ physical parameters is to protect the Himalayan marmots from natural enemies and bad weather36. There is clearly influence of slope on habitat selection by marmots. When the slope is large, wind is strong, and burrows are not well hidden; this makes them difficult to defend against enemies, unsafe for survival, and not conducive to hibernation during winter. In addition, Himalayan marmots prefer to burrow on sunny aspect, because the temperature is suitable and the vegetation is lush, which is suitable for marmots to breed. Therefore, the number of marmot burrows gradually decreases with increasing slope and ubac. Although flat and low-lying areas with small slopes are good for marmots to create dens, rainwater will easily flow into the dens during summer rainfall, which will kill marmots. Therefore, a suitable slope and sunny aspect are also very important for habitat selection by marmots.Application of the predictive spatial distribution map of Himalayan marmots in Qinghai provincePlague surveillance is the main measure used for plague prevention and control in China. Although we have made many improvements in plague surveillance, the traditional method of dragnet surveillance still consumes a lot of human and material resources, is inefficient. The pasture area of Qinghai province is approximately 380,000 km2, and the identified natural plague focus is approximately 180,000 km2; therefore, there is still 200,000 km2 of pasture where the distribution of Himalayan marmots and plague have not been identified. Currently, RS technology is widely used in the fields of mapping and ecological surveillance18,19,21,22,37.Applications of RS technology in areas such as malaria, dengue, schistosomiasis and plague have been previously reported27,37. Using GIS combined with remotely sensed data, Proches Hieronimo et al. found that the presence of small mammals was positively influenced by elevation, whereas the presence of fleas was clearly influenced by land management features, and thus these observations have positive implications for plague surveillance27. In this study, RS technology combined with field validations were used to determine the distribution and areas of different types of grasslands in Qinghai province, and the average density of Himalayan marmot distribution in different types of grasslands. The high-, low-, and very low-density areas of Himalayan marmot distribution were identified. The soil map, vegetation map, administrative map, and marmot density statistics were merged to form the spatial data and attribute data basis for the information system to map the distribution of Himalayan marmot and determine the area of Himalayan marmot distribution. Generally speaking, the occurrence of human plague epidemic is closely related to the local animal plague epidemic2. However, a large part of the high-density distribution of Himalayan marmots is located in uninhabited areas and the areas are generally sparsely populated, which also indicates that we should reasonably allocate plague prevention and control resources to areas where human plague is most likely to occur to prevent the occurrence of human plague epidemics.Field validation for verificationThrough field validation and information from local farmers and herdsmen, we confirmed that Himalayan marmots inhabited 68 sample sites in Tongde, Zeku, Guinan, Xunhua, Haiyan, Ulan, Qilian, Hualong, and Huzhu counties. Among them, Tongde, Zeku, Guinan, Xunhua, Haiyan, Ulan, and Qilian counties have all historically experienced marmot plague outbreaks and can be considered as reliable natural plague foci38. The data from this field validation are consistent with the previous survey data and the epidemic history of the counties in Qinghai province39.MAE can better reflect the actual number of errors in prediction values; the smaller the MAE value, the higher the prediction accuracy. The MAE derived from the field validation data was 0.1331 and the prediction accuracy was 0.8669. The accuracy of the predicted Himalayan marmot spatial distribution reached 87%, which indicated that the predicted probability map of the Himalayan marmot spatial distribution can better predict the potential marmot distribution.The predicted spatial distribution map of Himalayan marmot in Qinghai province was then compared with environmental information such as elevation, vegetation, grass type, slope, and aspect of 352 field survey sites. The obtained RS data showed that the prediction results were excellent, and the predicted spatial distribution map of Himalayan marmot in Qinghai province was drawn with high accuracy. The prediction map visually reflects the different density distribution of Himalayan marmots; this allows us to optimize the settings and reasonable spatial layout of animal plague surveillance sites and improve surveillance efficiency.Application of marmot information collection system V3.0Marmot information collection system V3.0 was developed based on the “3S” technology standardizing the collection of surveillance data, and makes the management and analysis of information more convenient and faster. This study revolutionized the traditional method of considering plague-stricken counties as the plague foci, and effectively reduces the work intensity of operators and improves the data collection efficiency. In 2016 and 2017, we applied this system to the animal plague surveillance tasks in the plague-stricken counties of Haidong, Hainan, and Haibei in Qinghai province, and standardized the collection of provincial geographic location data of animal plague surveillance (data not shown). In 2018, we also applied this system in Wulan County, which frequently experiences plague, and achieved a good application effect (data not shown).In the next step, we will expand the pilot areas (mainly national and provincial plague surveillance sites), collect surveillance data from each surveillance site, continuously optimize and update the system, improve the efficiency of data analysis and utilization, detect the plague epidemic in marmot in a timely and accurate manner, correctly determine the epidemic trend of plague in marmots, and attempt to strictly prevent the plague from spreading to humans. We plan to use a new model of drone surveillance to create a multidimensional, three-dimensional, real-time big data plague surveillance information reporting system to enhance early plague warnings and prediction in Qinghai province and even in the country, which will be of positive practical significance to serve and guarantee the Belt and Road Initiative. These approaches are expected to provide new technical means for plague investigation and research, and to provide references for setting up plague surveillance programs and prediction for the natural Himalayan marmot plague focus in Qinghai province and the QTP. More

Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251 (2001).Article 
CAS 

Google Scholar 
Davidson, E. A. & Janssens, I. A. Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. Nature 440, 165–173 (2006).Article 
CAS 

Google Scholar 
Lopez-Urrutia, A., San Martin, E., Harris, R. P. & Irigoien, X. Scaling the metabolic balance of the oceans. Proc. Natl Acad. Sci. USA 103, 8739–8744 (2006).Article 
CAS 

Google Scholar 
Yvon-Durocher, G. et al. Reconciling the temperature dependence of respiration across timescales and ecosystem types. Nature 487, 472–476 (2012).Article 
CAS 

Google Scholar 
Crowther, T. W. et al. Quantifying global soil carbon losses in response to warming. Nature 540, 104–108 (2016).Article 
CAS 

Google Scholar 
Bar-On, Y. M., Phillips, R. & Milo, R. The biomass distribution on Earth. Proc. Natl Acad. Sci. USA 115, 6506–6511 (2018).Article 
CAS 

Google Scholar 
Rivkin, R. B. & Legendre, L. Biogenic carbon cycling in the upper ocean: effects of microbial respiration. Science 291, 2398–2400 (2001).Article 
CAS 

Google Scholar 
Friedlingstein, P. et al. Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks. J. Clim. 27, 511–526 (2014).Article 

Google Scholar 
Smith, T. P. et al. Community-level respiration of prokaryotic microbes may rise with global warming. Nat. Commun. 10, 5124 (2019).Article 

Google Scholar 
Antwis, R. E. et al. Fifty important research questions in microbial ecology. FEMS Microbiol. Ecol. 93, fix044 (2017).Bardgett, R. D., Freeman, C. & Ostle, N. J. Microbial contributions to climate change through carbon cycle feedbacks. ISME J. 2, 805–814 (2008).Article 
CAS 

Google Scholar 
Enquist, B. J. et al. Scaling from traits to ecosystems: developing a general trait driver theory via integrating trait-based and metabolic scaling theories. Adv. Ecol. Res. 52, 249–318 (2015).Article 

Google Scholar 
Allen, A. P., Gillooly, J. F. & Brown, J. H. Linking the global carbon cycle to individual metabolism. Funct. Ecol. 19, 202–213 (2005).Article 

Google Scholar 
Schramski, J. R., Dell, A. I., Grady, J. M., Sibly, R. M. & Brown, J. H. Metabolic theory predicts whole-ecosystem properties. Proc. Natl Acad. Sci. USA 112, 2617–2622 (2015).Article 
CAS 

Google Scholar 
Alster, C. J., Koyama, A., Johnson, N. G., Wallenstein, M. D. & von Fischer, J. C. Temperature sensitivity of soil microbial communities: an application of macromolecular rate theory to microbial respiration. J. Geophys. Res. Biogeosci. 121, 1420–1433 (2016).Article 

Google Scholar 
Yvon-Durocher, G. et al. Five years of experimental warming increases the biodiversity and productivity of phytoplankton. PLoS Biol. 13, e1002324 (2015).Article 

Google Scholar 
Garzke, J., Connor, S. J., Sommer, U. & O’Connor, M. I. Trophic interactions modify the temperature dependence of community biomass and ecosystem function. PLoS Biol. 17, e2006806 (2019).Foster, K. R. & Bell, T. Competition, not cooperation, dominates interactions among culturable microbial species. Curr. Biol. 22, 1845–1850 (2012).Article 
CAS 

Google Scholar 
Coyte, K. Z., Schluter, J. & Foster, K. R. The ecology of the microbiome: networks, competition, and stability. Science 350, 663–666 (2015).Article 
CAS 

Google Scholar 
Machado, D. et al. Polarization of microbial communities between competitive and cooperative metabolism. Nat. Ecol. Evol. 5, 195–203 (2021).Article 

Google Scholar 
Bradford, M. A. et al. Cross-biome patterns in soil microbial respiration predictable from evolutionary theory on thermal adaptation. Nat. Ecol. Evol. 3, 223–231 (2019).Article 

Google Scholar 
Garcia-Martin, E. E., McNeill, S., Serret, P. & Leakey, R. J. G. Plankton metabolism and bacterial growth efficiency in offshore waters along a latitudinal transect between the UK and Svalbard. Deep Sea Res. I 92, 141–151 (2014).Article 
CAS 

Google Scholar 
Davidson, E. A., Richardson, A. D., Savage, K. E. & Hollinger, D. Y. A distinct seasonal pattern of the ratio of soil respiration to total ecosystem respiration in a spruce-dominated forest. Glob. Change Biol. 12, 230–239 (2006).Article 

Google Scholar 
Dutkiewicz, S., Follows, M. J. & Bragg, J. G. Modeling the coupling of ocean ecology and biogeochemistry. Glob. Biogeochem. Cycles 23, GB4017 (2009).Article 

Google Scholar 
Follows, M. J., Dutkiewicz, S., Ward, B. & Follett, C. in Microbial Ecology of the Oceans 3rd edn (eds Gasol, J. & Kirchman, D.) Ch. 12 (John Wiley, 2018).Letten, A. D. & Stouffer, D. B. The mechanistic basis for higher-order interactions and non-additivity in competitive communities. Ecol. Lett. 22, 423–436 (2019).Article 

Google Scholar 
Grilli, J., Barabás, G., Michalska-Smith, M. J. & Allesina, S. Higher-order interactions stabilize dynamics in competitive network models. Nature 548, 210–213 (2017).Article 
CAS 

Google Scholar 
Maynard, D. S., Crowther, T. W. & Bradford, M. A. Competitive network determines the direction of the diversity–function relationship. Proc. Natl Acad. Sci. USA 114, 11464–11469 (2017).Article 
CAS 

Google Scholar 
Fiegna, F., Moreno-Letelier, A., Bell, T. & Barraclough, T. G. Evolution of species interactions determines microbial community productivity in new environments. ISME J. 9, 1235–1245 (2015).Article 

Google Scholar 
Lawrence, D. et al. Species interactions alter evolutionary responses to a novel environment. PLoS Biol. 10, e1001330 (2012).Article 
CAS 

Google Scholar 
Harcombe, W. R., Chacón, J. M., Adamowicz, E. M., Chubiz, L. M. & Marx, C. J. Evolution of bidirectional costly mutualism from byproduct consumption. Proc. Natl Acad. Sci. USA 115, 12000–12004 (2018).Article 
CAS 

Google Scholar 
Goldford, J. E. et al. Emergent simplicity in microbial community assembly. Science 361, 469–474 (2018).Article 
CAS 

Google Scholar 
Yvon-Durocher, G. et al. Methane fluxes show consistent temperature dependence across microbial to ecosystem scales. Nature 507, 488–491 (2014).Article 
CAS 

Google Scholar 
Fox, J. W. & Harpole, W. S. Revealing how species loss affects ecosystem function: the trait-based price equation partition. Ecology 89, 269–279 (2008).Article 

Google Scholar 
Kontopoulos, D., Smith, T. P., Barraclough, T. G. & Pawar, S. Adaptive evolution shapes the present-day distribution of the thermal sensitivity of population growth rate. PLoS Biol. 18, e3000894 (2020).Article 
CAS 

Google Scholar 
Wilson, W. G. & Lundberg, P. Biodiversity and the Lotka–Volterra theory of species interactions: open systems and the distribution of logarithmic densities. Proc. R. Soc. Lond. B 271, 1977–1984 (2004).Article 

Google Scholar 
Rossberg, A. G. in Food Webs and Biodiversity 181–191 (John Wiley & Sons, 2013).Schloss, P. D. et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541 (2009).Article 
CAS 

Google Scholar 
Garcia, F. C., Bestion, E., Warfield, R. & Yvon-Durocher, G. Changes in temperature alter the relationship between biodiversity and ecosystem functioning. Proc. Natl Acad. Sci. USA 115, 10989–10999 (2018).Article 
CAS 

Google Scholar 
Padfield, D., O’Sullivan, H. & Pawar, S. rTPC and nls.multstart: a new pipeline to fit thermal performance curves in R. Methods Ecol. Evol. 12, 1138–1143 (2021).Article 

Google Scholar  More

Study areaThe Qinghai–Tibet Plateau (26°00′-39°47′N, 73°19′-104°47′E), one of the most important pastoral areas in the world, straddles the southwest regions of China, and it includes 244 counties, which belong to six provinces: Tibet, Qinghai, Xinjiang, Gansu, Sichuan, and Yunnan. It is characterized by rich natural grassland resources, including desert steppes, alpine steppes, and alpine meadows (Fig. 1a). The grassland areas account for over 56% of this region34. The grassland plays a vital role in providing regional and national animal husbandry products and fodder35, which enables the local herders to obtain almost all of the resources required for survival36. The grazing density distribution is extremely unbalanced (Fig. 1a) owing to the high spatial heterogeneity of economic development (Fig. 1b-1) and grassland production (Fig. 1b-2), resulting from the differences in resources and environmental factors37. Over the past few decades, there has been a significant change in the number of livestock animals, and the number of sheep exceeded 160 million by 2020. Therefore, it is urgent to obtain a high-resolution gridded grazing dataset for its evaluating spatiotemporal changes and coordinating the relationship between human beings and the grassland ecosystem.Fig. 1Location of the Qinghai–Tibet Plateau: (a) grassland type and distribution, and grazing density (GD) in 244 counties; (b) spatial heterogeneity of economic development (ED) and grassland production (GP) in 244 counties. GD, ED, and GP are represented by sheep unit per grassland area per county (SU/hm2), human footprint index per pixel (HF/pixel) per county, and net primary production per grassland area per county (gC/m2), respectively.Full size imageFig. 2Methodological framework for grazing spatialization.Full size imageMethodological frameworkWe developed a methodological framework for high-resolution gridded grazing dataset mapping. The framework mainly includes four parts: (i) identifying features affecting grazing, (ii) extracting theoretical suitable grazing areas, (iii) building grazing spatialization model, and (iv) correcting the grazing spatialization dataset. Each step is explained in more detail below (Fig. 2).Step 1: Identifying features affecting grazingGrazing activities are affected by the spatial heterogeneity of resources and environmental factors, regulated by the grazing behavior of herders and the foraging behavior of herds, and restricted by ecological protection policies. Therefore, the specific implications of the 14 influencing factors from the above four aspects are presented in Table 1. These factors are necessary for spatializing the county-level grazing data.Table 1 The identified features affecting grazing.Full size tableStep 2: Extracting theoretical suitable grazing areasThe decision tree approach38 was adopted to extract the theoretical suitable grazing areas for further grazing spatialization (step 2 in Fig. 2). First, the potential grazing area was identified according to the boundary of the grassland ecosystem, because grazing behavior only occurs in the grassland. Then, the unsuitable areas for grazing, i.e., extremely-high-altitude areas and areas adjacent to towns, were removed from the potential grazing area stepwise. The areas strictly prohibited for grazing, i.e., the core areas of national nature reserves39 within grassland areas, were also deemed unsuitable for grazing. Finally, the extracted areas were the theoretically suitable grazing areas.Step 3: Building grazing spatialization model(i) Extracting cross-scale feature (CSFs)In the traditional method, the spatial resolution of the training data (i.e., the average value at the administrative level) differs from that of the predicting data (i.e., the value at the pixel level), and the trained model can only capture the characteristics within the training data. However, the extreme value of the predicting data inevitably exceeds the range of the training data, which can result in underestimation in these parts40. To reduce these mismatches, we built an improved method for CSFs extraction (Fig. 2, first part of step 3).First, the census grazing data are simply distributed from county level to pixel level using the weight of the absolute disturbance (AD) index as Eq. (1). The AD index is measured by Mahalanobis distance using Eq. (2), which is calculated according to the deviation between the potential and observed normalized difference vegetation index (NDVI) values22. Second, the distributed grazing data are graded via the hierarchical clustering method, and the optimal number of the group can be determined using the Davies–Bouldin index (DBI)41 as Eq. (3), an index for evaluating the quality of clustering algorithm. The smaller the DBI, the smaller the distance within each group. Therefore, the DBI can be used to select the best similar values to minimize the deviation within each group. Finally, we can obtain the scope of the groups within each county using the above two steps and obtain the average value of all independent variables and the dependent variable accordingly. As expected, we can decompose the average value at the county level (traditional features in Fig. 2) into the average value at the group level (improved features in Fig. 2).$$S{U}_{i}=S{U}_{j}^{C}frac{{w}_{A{D}_{i}}}{{w}_{A{D}_{j}}}$$
(1)
where SUi and (S{U}_{j}^{C}) are the grazing value for pixel i and the census grazing value for county j; ({w}_{A{D}_{i}}) is the weight of the AD index for pixel i and ({w}_{A{D}_{j}}) represents the summed weight of the AD index values for all pixels in county j.$$begin{array}{cll}A{D}_{i} & = & sqrt{{({D}_{i}-u)}^{T}co{v}^{-1}({D}_{i}-u)}\ {D}_{i} & = & NDV{I}_{i}^{T}-NDV{I}_{i}^{P}end{array}$$
(2)
where ADi is the AD index for pixel i; the vector composed of its observed NDVI (left(NDV{I}_{i}^{T}right)) and potential NDVI (left(NDV{I}_{i}^{P}right)) time-series data could be considered as two points in the feature space for pixel i, and Di and u are the difference and the mean value of the vector, respectively; cov is the covariance matrix.$$DB{I}_{k}=frac{1}{k}{sum }_{x=1}^{k}ma{x}_{yne x}left(frac{overline{{a}_{x}}+overline{{a}_{y}}}{left|{delta }_{x}-{delta }_{y}right|}right)$$
(3)
where DBIk is the DBI coefficient when the cluster number is k; (overline{{a}_{x}}) and (overline{{a}_{y}}) are the average distances of the group xth and the group yth, respectively; δx and δy are the center distance of the group xth and the group yth, respectively.Different from the traditional method, our method can decompose features into multiple features using the grading AD index. The differences among counties will not be easily averaged out. Moreover, our method is less affected by scale mismatch and can be transferred to cross-scale modeling26.(ii) Building RF model with partitioningA single model cannot accurately obtain the variation information of the Qinghai–Tibet Plateau with high spatial heterogeneity. The partition model, a widely used method for estimating population distribution and others42,43, can be incorporated into the proposed model to improve its performance. The thresholds (0.43, 0.35 and 0.21 SU/hm2), determined according to the theoretical livestock carrying capacity (equation S1), were calculated and used to separate independent variables and dependent variable for each grassland types: alpine meadow, alpine steppe and alpine desert steppe (see Section 6.1 for details). Then, the RF models were established, and the training and testing samples were randomly divided in the proportion of 3:1. It is notable that transforming the response variable using natural log prior to RF model fitting is necessary to achieve higher prediction accuracies44. Finally, the independent variables at the pixel level were inputted into the two trained RF models, and the corresponding grid grazing dataset was output by combining the two results (Fig. 2, second part of step 3).(iii) Validating the accuracy of the methodsThe performance of the grazing spatialization model was evaluated through a comparison of the predicted value with census value26. Accuracy validation indexes, including coefficients of determination (R2), root mean square error (RMSE), and mean absolute error (MAE), were used to evaluate the performances of the proposed RF-based models (Table 2), as presented in Eq. (4).$$begin{array}{ccc}{R}^{2} & = & 1-frac{{sum }_{j=1}^{N}{left(S{U}_{j}^{C}-S{U}_{j}^{P}right)}^{2}}{{sum }_{j=1}^{N}{left(S{U}_{j}^{C}-overline{S{U}^{C}}right)}^{2}}\ RMSE & = & sqrt{frac{{sum }_{j=1}^{N}{left(S{U}_{j}^{C}-S{U}_{j}^{P}right)}^{2}}{N}}\ MAE & = & frac{{sum }_{j=1}^{N}| S{U}_{j}^{C}-S{U}_{j}^{P}| }{N}end{array}$$
(4)
where (S{U}_{j}^{C}) and (S{U}_{j}^{P}) are the census grazing value and the predicted grazing value for county j, respectively; (overline{S{U}^{C}}) is the average census data for all counties; and N is the number of all counties.Table 2 The proposed methods and their descriptions.Full size tableStep 4: Correcting grazing spatialization dataset(i) Correcting residuals of datasetCorrecting residuals is necessary to obtain datasets with higher accuracy45,46, because propagating the cross-scale relationship in the RF models will inevitably generate errors47. The residuals, calculated by the difference between the average census grazing and predicted grazing values at the administrative level, were used to calibrate the errors related to all pixels within this county. The revised dataset after residual correction is the final product provided in this study. The residual correction method is expressed by Eq. (5), and the process is shown in the fourth step in Fig. 2.$$S{U}_{i}^{RP}=S{U}_{i}^{P}+{R}_{j}$$
(5)
where (S{U}_{i}^{RP}) denotes the predicted grazing value revised by the residuals for pixel i, (S{U}_{i}^{P}) denotes the predicted grazing for pixel i, and Rj denotes the residuals calculated from the difference between census grazing and predicted grazing data for county j.(ii) Validating the accuracy of datasetTwo goodness-of-fit indexes were used to validate the consistency of spatial distribution and the temporal trend between predicted grazing data and census grazing data. Generally, the coefficient of determination (R2), defined in Eq. (4), is used to verify the consistency of spatial distribution, and the Nash–Sutcliffe efficiency (NSE, Eq. (6)) is used to verify the consistency of temporal trend. An index value closer to 1 corresponds to a more accurate dataset. Meanwhile, we also collected field grazing data from 56 sites to further validate the spatial accuracy of the dataset, and it measured using the R2 in Eq. (4).$$NSE=1-frac{{sum }_{t=1}^{T}{left(S{U}_{t}^{RP}-S{U}_{t}^{C}right)}^{2}}{{sum }_{t=1}^{T}{left(S{U}_{t}^{C}-overline{S{U}^{{C}^{{prime} }}}right)}^{2}}$$
(6)
where (S{U}_{t}^{RP}) and (S{U}_{t}^{C}) are the predicted grazing value revised by residuals and the census grazing value of all counties in year t, respectively; (overline{S{U}^{{C}^{{prime} }}}) is the average census grazing value of all years; and T is the number of time steps.(iii) Identifying uncertainties associated with datasetThe uncertainties associated with the dataset originate from the following two aspects: First, the unreasonableness of our method, owing to the errors related to cross-scale modeling or the inappropriate selection of influencing factors, is an important source of uncertainties. Second, the incompleteness of auxiliary variables also introduces uncertainties. In this instance, grassland-free areas are not accurately identified in some counties, but livestock animals are raised in these counties. These counties have no effective value for livestock density prediction. Overall, the uncertainties can be identified in terms of the mean relative error (MRE) in Eq. (7).$$MRE=frac{{sum }_{j=1}^{N}left|frac{S{U}_{j}^{C}-S{U}_{j}^{RP}}{S{U}_{j}^{C}}right|}{N}ast 100 % $$
(7)
where (S{U}_{j}^{C}) is the census grazing value for county j, (S{U}_{j}^{RP}) is the predicted grazing value revised by residuals for county j, and N is the number of counties.Data sourceCensus grazing data at county levelEight types of livestock, namely cattle, yaks, horses, donkeys, mules, camels, goats, and sheep, were considered according to the regional characteristics, and livestock stocking quantity at the end of year for each county can be determined from statistical yearbooks. However, the numbers of livestock at the county level for some years between 1982 and 2015 were not recorded. The missing data were indirectly approximated from city- or provincial-level data (e.g., interpolation using their temporal trends). Each type of livestock stocking quantity was converted into standard sheep unit (SU) according to the national standards using Eq. (8)48, namely the calculation of rangeland carrying capacity (NY/T 635-2015). Of the 244 counties of the Qinghai–Tibet Plateau, only 242 counties were considered, as the census grazing data for the other 2 counties were unavailable. The unit of grazing statistics data at the county level is defined as SU per county per year (SU·county−1·year−1).$$begin{array}{l}SU={N}_{sheep}+0.8times {N}_{goats}+5times {N}_{cattle}+5times {N}_{yaks+}+\ 6times {N}_{horses}+3times {N}_{donkeys}+6times {N}_{mules}+7times {N}_{camels}end{array}$$
(8)
where Nsheep, Ngoats, Ncattle, Nyaks, Nhorses, Ndonkeys, Nmules, Ncamels are the number of sheep, goats, cattle, yaks, horses, donkeys, mules, and camels at the year-end, respectively. SU denotes the standard sheep unit (SU·county−1·year−1).Data of grazing influencing factors at pixel levelThe types of features affecting grazing were obtained from the first step described in Methods, and the detailed information, such as original spatiotemporal resolution, format, and source, is shown in Table 3. The format (i.e., GeoTIFF), spatial resolution (i.e., 0.083°), and the number of rows and columns of the gridded features were leveraged to further produce a high-resolution grazing dataset.Table 3 Data source of grazing influence factors.Full size table More

Now consider a generalized model in which the species interactions are heterogeneous. A natural way of introducing heterogeneity in the system is by having a species diversify into subpopulations with different interaction strengths12,13,14,15. This way of modeling heterogeneity is useful as it can describe different kinds of heterogeneity. For example, the subpopulations could represent polymorphic traits that are genetically determined or result from plastic response to heterogeneous environments. A population could also be divided into local subpopulations in different spatial patches, which can migrate between patches and may face different local predators. We can also model different behavioral modes as subpopulations that, for instance, spend more time foraging for food or hiding from predators. We study several kinds of heterogeneity after we introduce a common mathematical framework. By studying these different scenarios using variants of the model, we show that our main results are not sensitive to the details of the model.We focus on the simple case where only the prey species splits into two types, (C_1) and (C_2), as illustrated in Fig. 1b. The situation is interesting when predator A consumes (C_1) more readily than predator B and B consumes (C_2) more readily than A (i.e., (a_1 / a_0 > b_1 / b_0) and (b_2 / b_0 > a_2 / a_0), which is equivalent to the condition that the nullclines of A and B cross, see section “Resources competition and nullcline analysis”). The arrows between (C_1) and (C_2) in Fig. 1b represent the exchange of individuals between the two subpopulations, which can happen by various mechanisms considered below. Such exchange as well as intraspecific competition between (C_1) and (C_2) result from the fact that the two prey types remain a single species.The system is now described by an enlarged Lotka-Volterra system with four variables, A, B, (C_1), and (C_2): $$begin{aligned} dot{A}&= varepsilon _A ,alpha _{A1} , A , C_1 + alpha _{A2} , A , C_2 – beta _A , A end{aligned}$$
(3a)
$$begin{aligned} dot{B}&= varepsilon _B , alpha _{B1} , B , C_1 + alpha _{B2} , B , C_2 – beta _B , B end{aligned}$$
(3b)
$$begin{aligned} dot{C_1}&= C_1 , (beta _C – alpha _{CC} , C)-alpha _{A1} , C_1 A-alpha _{B1} , C_1 B – sigma _1 , C_1 + sigma _2 , C_2 end{aligned}$$
(3c)
$$begin{aligned} dot{C_2}&= C_2 , (beta _C – alpha _{CC} , C) -alpha _{A2} , C_2 A -alpha _{B2} , C_2 B + sigma _1 , C_1 – sigma _2 , C_2 end{aligned}$$
(3d)
The parameters in these equations and their meanings are listed in Table 1. Here we assume that the prey types (C_1) and (C_2) have the same birth rate and intraspecific competition strength, but different interaction strengths with A and B. Note that (C_1) and (C_2) are connected by the (sigma _i) terms, which represent the exchange of individuals between these subpopulations through mechanisms studied below; these terms indicate a major difference between our model of a prey with intraspecific heterogeneity and other models of two prey species. For the convenience of analysis, we transform the variables (C_1) and (C_2) to another pair of variables C and (lambda), where (C equiv C_1 + C_2) is the total population of C as before, and (lambda equiv C_2 / (C_1 + C_2)) represents the composition of the population (Fig. 1c). After this transformation and rescaling of variables (described in “Methods”), the new dynamical system can be written as: $$begin{aligned} dot{A}&= A , big ( C , (a_1 (1-lambda ) + a_2 lambda ) – a_0 big ) end{aligned}$$
(4a)
$$begin{aligned} dot{B}&= B , big ( C , (b_1 (1-lambda ) + b_2 lambda ) – b_0 big ) end{aligned}$$
(4b)
$$begin{aligned} dot{C}&= C , big ( 1 – C – A (a_1 (1-lambda ) + a_2 lambda ) – B (b_1 (1-lambda ) + b_2 lambda ) big ) end{aligned}$$
(4c)
$$begin{aligned} dot{lambda }&= lambda (1-lambda ) , big ( A (a_1 – a_2) + B (b_1 – b_2) big ) + eta _1 (1-lambda ) – eta _2 lambda end{aligned}$$
(4d)
Here, (a_i) and (b_i) are the (rescaled) feeding rates of the predators on the prey type (C_i); (a_0) and (b_0) are the death rates of the predators as before; (eta _1) and (eta _2) are the exchange rates of the prey types (Table 1). The latter can be functions of other variables, representing different kinds of heterogeneous interactions that we study below. Notice that Eqs. (4a–4c) are equivalent to the homogeneous Eqs. (2a–2c) but with effective interaction strengths (a_text {eff} = (1-lambda ) , a_1 + lambda , a_2) and (b_text {eff} = (1-lambda ) , b_1 + lambda , b_2) that both depend on the prey composition (lambda) (Fig. 1c).Table 1 Model parameters (before/after rescaling) and their meanings.Full size tableThe variable (lambda) can be considered an internal degree of freedom within the C population. In all of the models we study below, (lambda) dynamically stabilizes to a special value (lambda ^*) (a bifurcation point), as shown in Fig. 3a. Accordingly, a new equilibrium point (P_N) appears (on the line (mathscr {L}) in Fig. 2), at which all three species coexist. For comparison, Fig. 3b shows the equilibrium points if (lambda) is held fixed at any other values, which all result in the exclusion of one of the predators. Thus, heterogeneous interactions give rise to a new coexistence phase (see Fig. 4 below) by bringing the prey composition (lambda) to the value (lambda ^*), instead of having to fine-tune the interaction strengths. The exact conditions for this new equilibrium to be stable are detailed in “Methods”.Figure 3(a) Time series of (lambda) for systems with each kind of heterogeneity. All three systems stabilize at the same (lambda ^*) value, which is the bifurcation point in panel (b). (b) Equilibrium population of each species (X = A), B, or C, with (lambda) held fixed at different values. Solid curves represent stable equilibria and dashed curves represent unstable equilibria (see Eq. (9) in “Methods”). The vertical dashed line is where (lambda = lambda ^*), which is also the bifurcation point. Notice that the equilibrium population of C is maximized at this point (for (a_1 > a_2) and (b_2 > b_1)). Parameters used here are ((a_0, a_1, a_2, b_0, b_1, b_2, rho , eta _1, eta _2, kappa ) = (0.25, 0.5, 0.2, 0.4, 0.2, 0.6, 0.5, 0.05, 0.05, 50)).Full size imageInherent heterogeneityWe first consider a scenario where individuals of the prey species are born as one of two types with a fixed ratio, such that a fraction (rho) of the newborns are (C_2) and ((1-rho )) are (C_1). This could describe dimorphic traits, such as the winged and wingless morphs in aphids12 or the horned and hornless morphs in beetles13. We call this “inherent” heterogeneity, because individuals are born with a certain type and cannot change in later stages of life. The prey type given at birth determines the individual’s interaction strength with the predators. This kind of heterogeneity can be described by Eq. (4d) with (eta _1 = rho (1-C)) and (eta _2 = (1-rho ) (1-C)) (see “Methods”).Figure 4Phase diagrams showing regions of parameter space identified by the stable equilibrium points. Yellow region represents (P_C) (predators A, B both extinct), red represents (P_A) (A excludes B), blue represents (P_B) (B excludes A), and green represents (P_N) (A, B coexist). The middle point (black dot) is where the preferences of the two predators are identical, (a_2/a_0=b_2/b_0) and (b_1/b_0=a_1/a_0). The coexistence phase appears in all three kinds of heterogeneity modeled here. (a–d) Inherent heterogeneity: Individuals of the prey population are born in two types with a fixed composition (rho). In the extreme cases of (rho = 0) and 1, the prey is homogeneous and there is no coexistence of the predators. (e–h) Reversible heterogeneity: Individual prey can switch types with fixed switching rates (eta _1) and (eta _2). As the switching rates increase, the coexistence region shrinks because the prey population becomes effectively homogeneous (the occasional green spots are numerical artifacts because the time to reach the equilibrium becomes long in this limit). (i–l) Adaptive heterogeneity: The switching rates (eta _i) dynamically adapt to the predator densities, so as to maximize the growth rate of the prey. As the sharpness (kappa) of the sigmoidal decision function is increased, the prey adapts more optimally and the region of coexistence expands. Parameters used here are ((a_0, a_1, b_0, b_2) = (0.3, 0.5, 0.4, 0.6)).Full size imageThe stable equilibrium of the system can be represented by phase diagrams that show the identities of the species at equilibrium. We plot these phase diagrams by varying the parameters (a_2) and (b_1) while keeping (a_1) and (b_2) constant. As shown in Fig. 4a–d, the equilibrium state depends on the parameter (rho). In the limit (rho = 0) or 1, we recover the homogeneous case because only one type of C is produced. The corresponding phase diagrams (Fig. 4a, d) contain only two phases where either of the predators is excluded, illustrating the competitive exclusion principle. For intermediate values of (rho), however, there is a new phase of coexistence that separates the two exclusion phases (Fig. 4b, c). There are two such regions of coexistence, which touch at a middle point and open toward the bottom left and upper right, respectively. The middle point is at ((a_2/a_0 = b_2/b_0, b_1/b_0 = a_1/a_0)), where the feeding preferences of the two predators are identical (hence their niches fully overlap). Towards the origin and the far upper right, the predators consume one type of C each (hence their niches separate). The coexistence region in the bottom left is where the feeding rates of the predators are the lowest overall. There can be a region (yellow) where both predators go extinct, if their primary prey type alone is not enough to sustain each predator. Increasing the productivity of the system by increasing the birth rate ((beta _C)) of the prey eliminates this extinction region, whereas lowering productivity causes the extinction region to take over the lower coexistence region. Because the existence and identity of the phases is determined by the configuration of the equilibrium points (Fig. 2, see also section “Mathematical methods”), the qualitative shape of the phase diagram is not sensitive to changes of parameter values.The new equilibrium is not only where the predators A and B can coexist, but also where the prey species C grows to a larger density than what is possible for a homogeneous population. This is illustrated in Fig. 3b, which shows the equilibrium population of C if we hold (lambda) fixed at different values. The point (lambda = lambda ^*) is where the system with a dynamic (lambda) is stable, and also where the population of C is maximized (when A and B prefer different prey types). That means the population automatically stabilizes at the optimal composition of prey types. Moreover, the value of (C^*) at this coexistence point can even be larger than the equilibrium population of C when there is only one predator A or B. This is discussed further in section “Multiple-predator effects and emergent promotion of prey”. These results suggest that heterogeneity in interaction strengths can potentially be a strategy for the prey population to leverage the effects of multiple predators against each other to improve survival.Reversible heterogeneityWe next consider a scenario where individual prey can switch their types. This kind of heterogeneity can model reversible changes of phenotypes, i.e., trait changes that affect the prey’s interaction with predators but are not permanent. For example, changes in coat color or camouflage14,16,17, physiological changes such as defense15, and biomass allocation among tissues18,19. One could also think of the prey types as subpopulations within different spatial patches, if each predator hunts at a preferred patch and the prey migrate between the patches20,21. With some generalization, one could even consider heterogeneity in resources, such as nutrients located in different places, that can be reached by primary consumers, such as swimming phytoplankton22. We can model this “reversible” kind of heterogeneity by introducing switching rates from one prey type to the other. In Eq. (4d), (eta _1) and (eta _2) now represent the switching rates per capita from (C_1) to (C_2) and from (C_2) to (C_1), respectively. Here we study the simplest case where both rates are fixed.In the absence of the predators, the natural composition of the prey species given by the switching rates would be (rho equiv eta _1 / (eta _1 + eta _2)), and the rate at which (lambda) relaxes to this natural composition is (gamma equiv eta _1 + eta _2). Compared to the previous scenario where we had only one parameter (rho), here we have an additional parameter (gamma) that modifies the behavior of the system. Fig. 4e–h shows phase diagrams for the system as (rho) is fixed and (gamma) varies. We again find the new equilibrium (P_N) where all three species coexist. When (gamma) is small, the system has a large region of coexistence. As (gamma) is increased, this region is squeezed into a border between the two regions of exclusion, where the slope of the border is given by (eta _1/eta _2) as determined by the parameter (rho). However, this is different from the exclusion we see in the case of inherent heterogeneity, which happens only for (rho rightarrow 0) or 1, where the borders are horizontal or vertical (Fig. 4a,d). Here the predators exclude each other despite having a mixture of prey types in the population.This special limit can be understood as follows. For a large (gamma), (lambda) is effectively set to a constant value equal to (rho), because it has a very fast relaxation rate. In other words, the prey types exchange so often that the population always maintains a constant composition. In this limit, the system behaves as if it were a homogeneous system with effective interaction strengths (a_text {eff} = (1-rho ) , a_1 + rho , a_2) and (b_text {eff} = (1-rho ) , b_1 + rho , b_2). As in a homogeneous system, there is competitive exclusion between the predators instead of coexistence. This demonstrates that having a constant level of heterogeneity is not sufficient to cause coexistence. The overall composition of the population must be able to change dynamically as a result of population growth and consumption by predators.An interesting behavior is seen when we examine a point inside the shrinking coexistence region as (gamma) is increased. Typical trajectories of the system for such parameter values are shown in Fig. 5. As (gamma) increases, the system relaxes to the line (mathscr {L}) quickly, then slowly crawls along it towards the final equilibrium point (P_N). This is because increasing (gamma) increases the speed that (lambda) relaxes to (lambda ^*), and when (lambda rightarrow lambda ^*), (mathscr {L}) becomes marginally stable. Therefore, the attraction to (mathscr {L}) in the perpendicular direction is strong, but the attraction towards the equilibrium point along (mathscr {L}) is weak. This leads to a long transient behavior that makes the system appear to reach no equilibrium in a limited time23,24. It is especially true when there is noise in the dynamics, which causes the system to diffuse along (mathscr {L}) with only a weak drift towards the final equilibrium (Fig. 5). Thus, the introduction of a fast timescale (quick relaxation of (lambda) due to a large (gamma)) actually results in a long transient.Figure 5Trajectories of the system projected in the A-B plane, with parameters inside the coexistence region (by holding the position of (P_N) fixed). As (gamma) increases, the system tends to approach the line (mathscr {L}) quickly and then crawl along it. The grey trajectory is with independent Gaussian white noise ((sim mathscr {N}(0,0.5))) added to each variable’s dynamics. Noise causes the system to diffuse along (mathscr {L}) for a long transient period before coming to the equilibrium point (P_N). Parameters used here are ((a_0, a_1, a_2, b_0, b_1, b_2) = (0.2, 0.8, 0.5, 0.2, 0.6, 0.9)), chosen to place (P_N) away from the middle of (mathscr {L}) to show the trajectory drifting toward the equilibrium.Full size imageAdaptive heterogeneityA third kind of heterogeneity we consider is the change of interactions in time. By this we mean an individual can actively change its interaction strength with others in response to certain conditions. This kind of response is often invoked in models of adaptive foraging behavior, where individuals choose appropriate actions to maximize some form of fitness25,26. For example, we may consider two behaviors, resting and foraging, as our prey types. Different predators may prefer to strike when the prey is doing different things. In response, the prey may choose to do one thing or the other depending on the current abundances of different predators. Such behavioral modulation is seen, for example, in systems of predatory spiders and grasshoppers27. Phenotypic plasticity is also seen in plant tissues in response to consumers28,29,30.This kind of “adaptive” heterogeneity can be modeled by having switching rates (eta _1) and (eta _2) that are time-dependent. Let us assume that the prey species tries to maximize its population growth rate by switching to the more favorable type. From Eq. (4c), we see that the growth rate of C depends linearly on the composition (lambda) with a coefficient (u(A,B) equiv (a_1 – a_2) A + (b_1 – b_2) B). Therefore, when this coefficient is positive, it is favorable for C to increase (lambda) by switching to type (C_2). This can be achieved by having a positive switching rate (eta _2) whenever (u(A,B) > 0). Similarly, whenever (u(A,B) < 0), it is favorable for C to switch to type (C_1) by having a positive (eta _1). In this way, the heterogeneity of the prey population constantly adapts to the predator densities. We model such adaptive switching by making (eta _1) and (eta _2) functions of the coefficient u(A, B), e.g., (eta _1(u) = 1/(1+mathrm {e}^{kappa u})) and (eta _2(u) = 1/(1+mathrm {e}^{-kappa u})). The sigmoidal form of the functions means that the switching rate in the favorable direction for C is turned on quickly, while the other direction is turned off. The parameter (kappa) controls the sharpness of this transition.Phase diagrams for the system with different values of (kappa) are shown in Fig. 4i–l. A larger (kappa) means the prey adapts its composition faster and more optimally, which causes the coexistence region to expand. In the extreme limit, the system changes its dynamics instantaneously whenever it crosses the boundary where (u(A,B) = 0), like in a hybrid system31. Such a system can still reach a stable equilibrium that lies on the boundary, if the flow on each side of the boundary points towards the other side32. This is what happens in our system and, interestingly, the equilibrium is the same three-species coexistence point (P_N) as in the previous scenarios. The region of coexistence turns out to be largest in this limit (Fig. 4l).Our results suggest that the coexistence of the predators can be viewed as a by-product of the prey’s strategy to maximize its own benefit. The time-dependent case studied here represents a strategy that involves the prey evaluating the risk posed by different predators. This is in contrast to the scenarios studied above, where the prey population passively creates phenotypic heterogeneity regardless of the presence of the predators. These two types of behavior are analogous to the two strategies studied for adaptation in varying environments, i.e., sensing and bet-hedging33,34. The former requires accessing information about the current environment to make optimal decisions, whereas the latter relies on maintaining a diverse population to reduce detrimental effects caused by environmental changes. Here the varying abundances of the predators play a similar role as the varying environment. From this point of view, the heterogeneous interactions studied here can be a strategy of the prey species that is evolutionarily favorable. More

Ren, W. et al. Changes of periphyton abundance and biomass driven by factors specific to flooding inflow in a river inlet area in Erhai Lake, China. Front. Environ. Sci. 9, 680718. https://doi.org/10.3389/fenvs.2021.680718 (2021).Article 

Google Scholar 
Woodruff, S. L. et al. The effects of a developing biofilm on chemical changes across the sediment-water interface in a freshwater environment. Int. Rev. Hydrobiol. 84(5), 509–532 (1999).CAS 

Google Scholar 
Muñoz, I., Real, M., Guasch, H., Navarro, E. & Sabater, S. Effects of atrazine on periphyton under grazing pressure. Aquat. Toxicol. 55(3–4), 239–249 (2001).
Google Scholar 
Hoagland, K. D., Roemer, S. C. & Rosowski, J. R. Colonization and community structure of two periphyton assemblages, with emphasis on the diatoms (Bacillariophyceae). Am. J. Bot. 69, 188–213. https://doi.org/10.2307/2443006 (1982).Article 

Google Scholar 
Steinman, A. D. & McIntire, C. D. Effects of current velocity and light energy on the structure of periphyton assemblages in laboratory streams. J. Phycol. 22, 352–361. https://doi.org/10.1111/J.1529-8817.1986.TB00035.X (1986).Article 

Google Scholar 
Tonkin, J. D., Death, R. G. & Barquín, J. Periphyton control on stream invertebrate diversity: Is periphyton architecture more important than biomass?. Mar. Freshw. Res. 65(9), 818–829 (2014).
Google Scholar 
Beck, W. S., Markman, D. W., Oleksy, I. A., Lafferty, M. H. & Poff, N. L. Seasonal shifts in the importance of bottom-up and top-down factors on stream periphyton community structure. Oikos 128, 680–691. https://doi.org/10.1111/oik.05844 (2018).Article 
CAS 

Google Scholar 
Hogsden, K. L. & Harding, J. S. Consequences of acid mine drainage for the structure and function of benthic stream communities: A review. Freshw. Sci. 31, 108–120. https://doi.org/10.1899/11-091.1 (2012).Article 

Google Scholar 
Sofi, M. S., Bhat, S. U., Rashid, I. & Kuniyal, J. C. The natural flow regime: A master variable for maintaining river ecosystem health. Ecohydrology 13(8), e2247. https://doi.org/10.1002/eco.2247 (2020).Article 

Google Scholar 
Biggs, B. J. F. Eutrophication of streams and rivers: Dissolved nutrient-chlorophyllrelationship for benthic algae. J. N. Am. Benthol. Soc. 19, 17–31 (2000).
Google Scholar 
Ormerod, S. J., Dobson, M., Hildrew, A. G. & Townsend, C. Multiple stressors in freshwater ecosystems. Freshw. Biol. 55, 1–4 (2010).
Google Scholar 
Poff, et al. The natural flow regime: A paradigm for river conservation and restoration. Bioscience 47, 769–784 (1997).
Google Scholar 
Naiman, R. J., Décamps, H., & McClain, M. E. Riparia: Ecology, Conservation and Management of Streamside Communities, (Elsevier/Academic Press, 2005).Gleick, P. H. Water use. Annu. Rev. Environ. Resour. 28, 275–314 (2003).
Google Scholar 
Jenkins, K. M. & Boulton, A. J. Connectivity in a dryland river: Short-term aquatic macroinvertebrate recruitment following floodplain inundation. Ecology 84(10), 2708–2723 (2003).
Google Scholar 
Biggs, B. J. F. Patterns in benthic algae of streams. In Algal Ecology in Freshwater Benthic Ecosystems (eds. Stevenson, R. J., Bothwell, M. L., & Lowe, R. L.) 31–56 (Academic Press, 1996).Smolar-Žvanut, N. & Mikoš, M. The impact of flow regulation by hydropower dams on the periphyton community in the Soča River, Slovenia. Hydrol. Sci. J. 59(5), 1032–1045. https://doi.org/10.1080/02626667.2013.834339 (2014).Article 
CAS 

Google Scholar 
Curry, C. J. & Baird, D. J. Habitat type and dispersal ability influence spatial structuring of larval Odonata and Trichoptera assemblages. Freshw. Biol. 60, 2142–2152 (2015).
Google Scholar 
Wu, N., Cai, Q. & Fohrer, N. Contribution of microspatial factors to benthic diatom communities. Hydrobiologia 732, 49–60. https://doi.org/10.1007/s10750-014-1843-3 (2014).Article 
CAS 

Google Scholar 
Mueller, M., Pander, J. & Geist, J. The effects of weirs on structural stream habitat and biological communities. J. Appl. Ecol 48(6), 1450–1461. https://doi.org/10.1111/j.1365-2664.2011.02035.x (2011).Article 

Google Scholar 
Davies, P. M. et al. Flow–ecology relationships: closing the loop on effective environmental flows. Mar. Freshw. Res. 65(2), 133–141 (2013).
Google Scholar 
Jun, Y. C. et al. Spatial distribution of benthic macroinvertebrate assemblages in relation to environmental variables in Korean nationwide streams. Water 8(1), 27. https://doi.org/10.3390/w8010027 (2016).Article 

Google Scholar 
Biggs, B. J. F. & Close, M. E. Periphyton biomass dynamics in gravel bed rivers: The relative effects of flows and nutrients. Freshw. Biol. 22, 209–231 (1989).CAS 

Google Scholar 
Jowett, I. & Biggs, B. J. F. Flood and velocity effects on periphyton and silt accumulation in two New Zealand rivers. N. Zeal. J. Mar. Freshw. Res. 31, 287–300 (1997).
Google Scholar 
Biggs, B. J. F., Goring, D. G. & Nikora, V. I. Subsidy and stress responses of stream periphyton to gradients in water velocity as a function of community growth form. J. Phycol. 34, 598–607 (1998).
Google Scholar 
Malmqvist, B. & Englund, G. Effects of hydropower-induced flow perturbations on mayfly (Ephemeroptera) richness and abundance in north Swedish river rapids. Hydrobiologia 341(2), 145–158 (1996).
Google Scholar 
Poff, N. L. & Ward, J. V. Herbivory under different flow regimes: A field experiment and test of a model with a benthic stream insect. Oikos 72, 179–188 (1995).
Google Scholar 
Poff, L. N., Wellnitz, T. & Monroe, J. B. Redundancy among three herbivorous insects across an experimental current velocity gradient. Oecologia 134, 262–269. https://doi.org/10.1007/s00442-002-1086-2 (2003).Article 

Google Scholar 
Vaughn, C. C. The role of periphyton abundance and quality in the microdistribution of a stream grazer, Helicopsyche borealis (Trichoptera: Helicopsychidae). Freshw. Biol. 16, 485–493 (1986).
Google Scholar 
Francoeur, S. N. Meta-analysis of lotic nutrient amendment experiments: Detecting and quantifying subtle responses. J. N. Am. Benthol. Soc. 20, 358–368 (2001).
Google Scholar 
Elser, J. J. et al. Global analysis of nitrogen and phosphorus limitation of primary producers in freshwater, marine and terrestrial ecosystems. Ecol. Lett. 10, 1135–1142 (2007).
Google Scholar 
Hillebrand, H. Meta-analysis of grazer control of periphyton biomass across aquatic ecosystems. J. Phycol. 45, 798–806 (2009).
Google Scholar 
Lamberti, G. A. The role of periphyton in benthic food webs. In Algal Ecology—Freshwater Benthic Ecosystems, 533–572 (eds. Stevenson, R. J., Bothwell, M. L. & Lowe, R. L.) (Academic Press, 1996).Lamberti, G. A. et al. Influence of grazer type and abundance on plant–herbivore interactions in streams. Hydrobiologia 306, 179–188 (1995).
Google Scholar 
Gregory, S. V. Plant–herbivore interactions in stream systems. In Stream Ecology (eds. Barnes, J. R. & Minshall, G. W.) 157–189 (Plenum, 1983).Lamberti, G. A. & Moore, J. W. Aquatic insects as primary consumers. In The Ecology of Aquatic Insects (eds Resh, V. H. & Rosenberg, D. M.) 164–195 (Praeger, 1984).
Google Scholar 
Sterner, R. W., Elser, J. J. & Hessen, D. O. Stoichiometric relationships among producers, consumers and nutrient cycling in pelagic ecosystems. Biogeochemistry 17, 49–67 (1992).CAS 

Google Scholar 
Kahlert, M. & Baunsgaard, M. T. Nutrient recycling—A strategy of a grazer community to overcome nutrient limitation. J. N. Am. Benthol. Soc. 18, 363–369 (1999).
Google Scholar 
Burkholder, J. M., Wetzel, R. G. & Klomparens, K. L. Direct comparison of phosphate uptake by adnate and loosely attached microalgae within and intact biofilm matrix. Appl. Environ. Microbiol. 56, 2882–2890 (1990).CAS 

Google Scholar 
Steinman, A. D. Effects of grazers on freshwater benthic algae. In Algal Ecology: Freshwater Benthic Ecosystems (eds. Stevenson, R. J., Bothwell & Lowe, R. L.) 341–366 (Academic Press, 1996).Smucker, N. J. & Vis, M. L. Spatial factors contribute to benthic diatom structure in streams across spatial scales: Considerations for biomonitoring. Ecol. Indic. 11, 1191–1203 (2011).
Google Scholar 
Myers, et al. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).Article 
CAS 

Google Scholar 
Wang, J., Pan, F., Soininen, J., Heino, J. & Shen, J. Nutrient enrichment modifies temperature-biodiversity relationship in large scale field experiments. Nat. Commun. 7, 13 (2016).
Google Scholar 
Wu, et al. Flow regimes filter species traits of benthic diatom communities and modify the functional features of lowland streams-a nationwide scale study. Sci. Total Environ. 651, 357–366 (2019).CAS 

Google Scholar 
Nisar, M. A. Geospatial approach to study environmental characterization of a Kashmir wetland (Anchar) catchment with special reference to land use/land cover and changing climate. Ph.D Thesis, Sher-e-Kashmir University of Agricultural Sciences and Technology, Kashmir. Weblink. http://krishikosh.egranth.ac.in/handle/1/91309 (2012).Bhat, S. U., Sofi, A. H., Yaseen, T., Pandit, A. K. & Yousuf, A. R. Macro invertebrate community from Sonamarg streams of Kashmir Himalaya. Pak. J. Biol. Sci. 14(3), 182–194. https://doi.org/10.3923/pjbs.2011.182.194 (2011).Article 
CAS 

Google Scholar 
Baba, A. I., Sofi, A. H., Bhat, S. U., & Pandit, A. K. Periphytic algae of river Sindh in the Sonamarg area of Kashmir valley. J. Phytol. 3(6) (2011).Sofi, M. S., Rautela, K. S., Bhat, S. U., Rashid, I. & Kuniyal, J. C. Application of geomorphometric approach for the estimation of hydro-sedimentological flows and cation weathering rate: Towards understanding the sustainable land use policy for the Sindh Basin, Kashmir Himalaya. Water Air Soil Pollut. 232(7), 1–11. https://doi.org/10.1007/s11270-021-05217-w (2021).Article 
CAS 

Google Scholar 
Romshoo, S. A., & Fayaz, M. Use of high resolution remote sensing for improving environmental Friendly tourism master planning in the Alpine Himalaya: A case study of Sonamarg tourist resort, Kashmir. J. Himalayan Ecol. Sustain. Dev. 14 (2019).Biggs, B. J. F. & Kilroy, C. Stream periphyton monitoring manual. Published by NIWA for Ministry for the Environment, 226 Christchurch, New Zealand: NIWA (2000).APHA. Standard Methods for Examination of Water and Wastewater, 22nd edn. (American Public Health Association, 2012).Cox, E. J. Identification of Freshwater Diatoms from Live Material. (Chapman and Hall, 1996). https://doi.org/10.1017/S0025315400041023.Krammer, K., & Lange-Bertalot, H. Bacillariophyceae, Part 5. English and French Translation of the Keys. (VEB Gustav Fisher Verlag, 2000).Reichardt, E. A remarkable association of diatoms in a spring habitat in the Grazer Bergland, Austria. In Iconographia Diatomologica (ed. Lange-Bertalot, H.) 419–479 (2004).Żelazna-Wieczorek, J. Diatom flora in springs of Lódz Hills (Central Poland). Biodiversity, taxonomy and temporal changes of epipsammic diatom assemblages in springs affected by human impact, 419. Volume 13 of Diatom monographs. Gantner. https://books.google.co.in/books?id=bdxeewAACAAJ (2011).Stark, J. D., Boothroyd, I. K. G., Harding, J. S., Maxted, J. R. & Scarsbrook, M. R. Protocols for sampling macroinvertebrates in wadeable streams. In New Zealand Macroinvertebrate Working Group Report no. 1. Prepared for the Ministry for the Environment. Sustainable Management Fund Project, 5103 (2001).Winterbourn, M. J. Sampling stream invertebrates. In Biological Monitoring of Freshwaters. Proceedings of the Seminar. Water and Soil Miscellaneous Publication No. 83 (eds. Pridmore, R. D., Cooper, A. B.) 241–258. (National Water and Soil Conservation Authority, 1985).Barbour, M. T., Gerritsen, J., Snyder, B. D., Stribling, J. B. Rapid bioassessment protocols for use in streams and wadeable rivers: periphyton, benthic macroinvertebrates and fish, 339. (United States Environmental Protection Agency, Office of Water, 1999).Malmqvist, B. & Hoffsten, P. O. Macroinvertebrate taxonomic richness, community structure and nestedness in Swedish streams. Fundam. Appl. Limnol. 150(1), 29–54. https://doi.org/10.1127/archiv-hydrobiol/150/2000/29 (2000).Article 

Google Scholar 
Ilmonen, J. & Paasivirta, L. Benthic macrocrustacean and insect assemblages in relation to spring habitat characteristics: Patterns in abundance and diversity. Hydrobiologia 533(1–3), 99–113. https://doi.org/10.1007/s10750-004-2399-4 (2005).Article 

Google Scholar 
Munasinghe, D. S. N., Najim, M. M. M., Quadroni, S. & Musthafa, M. M. Impacts of streamflow alteration on benthic macroinvertebrates by mini-hydro diversion in Sri Lanka. Sci. Rep. 11(1), 546. https://doi.org/10.1038/s41598-020-79576-5 (2021).Article 
CAS 

Google Scholar 
Edmondson, W. T. Fresh-Water Biology, 2nd ed. 1050–1056 (Wiley, 1959).Pennak, R. W. Freshwater Invertebrates of United States. (Wiley, 1978).McCafferty, W. P., Provonsha, A. V. Aquatic entomology: The fishermen’s and ecologists’ Illustrated Guide to Insects and their Relatives. (Jones and Bartlett Publishers, 1983).Borror, D., Triplehorn, C., Johnson, N. An Introduction to the Study of Insects, 6th ed. (Saunders College Publishing, 1989).Ward, J. V. Aquatic Insect Ecology, Biology and Habitat. (Wiley, 1992).Engblom, E. & Lingdell, P.E. Analyses of Benthic Invertebrates (ed. Nyman, L.) (1999).Bouchard, R. W. Guide to Aquatic Invertebrates of the Upper Midwest: Identification Manual for Students (University of Minnesota, 2004).
Google Scholar 
Subramanian, K. A. & Sivaramakrishnan, K. G. Aquatic Insects for Biomonitoring Freshwater Ecosystems—A Methodology Manual. (Ashoka Trust for Ecology and Environment (ATREE), 2007).Thorp, J. H., & Covich, A. P. (eds.) Ecology and Classification of North American Freshwater Invertebrates. (Academic Press, 2009).Allan, J. D. & Castillo, M.M. An introduction to fluvial ecosystems. In Stream Ecology: Structure and Function of Running Waters, 1–12 (2007).Oksanen, et al. Vegan: Community ecology package. In: R package version 2.4-3.McCune, B. & Grace, B. Analysis of Ecological Communities (MjM Software Design, 2016).Clarke, K. R. & Gorley, R. N. Primer v6 Permanova+ (Primer-E Ltd., 2006).
Google Scholar 
Salazar, G. EcolUtils: Utilities for Community Ecology Analysis. R package version 0.1 software (2018).Anderson, M. J., Ellingsen, K. E. & McArdle, B. H. Multivariate dispersion as a measure of beta diversity. Ecol. Lett. 9(6), 683–693 (2006).
Google Scholar 
Gardener, M. Community Ecology: Analytical Methods in Using R and Excel. (Pelagic Publishing, 2014).Chao, A. & Bunge, J. Estimating the number of species in a stochastic abundance model. Biometrics 58, 531–539. https://doi.org/10.1111/j.0006-341X.2002.00531.x (2002).Article 
MATH 

Google Scholar 
Peres-Neto, P. R., Legendre, P., Dray, S. & Borcard, D. Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology 87, 2614–2625 (2006).
Google Scholar 
Meng, X. L. et al. Responses of macroinvertebrates and local environment to short-term commercial sand dredging practices in a flood-plain lake. Sci. Total Environ. 631, 1350–1359 (2018).
Google Scholar 
Core Team, R. R: A Language and Environmental for Statistical Computing. (R Foundation for Statistical Computing, 2017).Wood, P. J. & Armitage, P. D. Biological effects of fine sediment in the lotic environment. Environ. Manag. 21, 203–217 (1997).CAS 

Google Scholar 
Marchant, R. Changes in the benthic invertebrate communities of the Thomson River, southeastern Australia, after dam construction. Regul. Rivers Res. Manag. 4, 71–89 (1989).
Google Scholar 
Gray, L. J. & Ward, J. V. Effects of sediment releases from a reservoir on stream macroinvertebrates. Hydrobiologia 96, 177–184 (1982).
Google Scholar 
Sand-Jensen, K., Moller, J. & Olesen, B. H. Biomass regulation of microbenthic algae in Danish lowland streams. Oikos 53, 332–340 (1988).
Google Scholar 
Lewis, M. A., Weber, D. E., Stanley, R. S. & Moore, J. C. Dredging impact on an urbanized Florida bayou: Effects on benthos and algal-periphyton. Environ. Pollut. 115(2), 161–171 (2001).CAS 

Google Scholar 
Biggs, B. J. Algal ecology in freshwater benthic ecosystems geology and landuse to the habitat template of periphyton in stream ecosystems. Freshw. Biol. 33, 419–438 (1995).
Google Scholar 
Taylor, et al. Can diatom-based pollution indices be used for biomonitoring in South Africa? A case study of the Crocodile West and Marico water management area. Hydrobiologia 592, 455–464 (2007).
Google Scholar 
Porter, et al. Efficacy of algal metrics for assessing nutrient and organic enrichment in flowing waters. Freshw. Biol. 53, 1036–1054 (2008).
Google Scholar 
Wetzel, R. G. & Likens, G. E. Limnological analyses, 3rd ed. In Nitrogen, Phosphorus, and Other Nutrients, 85–113. (Springer, 2000). https://doi.org/10.1007/978-1-4757-3250-4.Wetzel, R. G. Attached algal-substrata interactions: Fact or myth, and when and how? vol. 17. In Periphyton of Freshwater Ecosystems (ed. Wetzel, R.) 207–215 (Springer, 1983). https://doi.org/10.1007/978-94-009-7293-3_28.Krajenbrink, H. J. et al. Diatoms as indicators of the effects of river impoundment at multiple spatial scales. PeerJ 7, e8092. https://doi.org/10.7717/peerj.8092 (2019).Article 

Google Scholar 
Poff, N. L., Voelz, N. J., Ward, J. V. & Lee, R. E. Algal colonization under four experimentally-controlled current regimes in a high mountain stream. J. N. Am. Benthol. Soc. 9, 303–318 (1990).
Google Scholar 
Dodds, W. K. & Marra, J. L. Behaviors of the midge, Cricotopus (Diptera; Chironomidae) related to mutualism with Nostoc parmeloides (Cyanobacteria). Aquat. Insects 11, 201–208 (1989).
Google Scholar 
Tang, T., Niu, S. Q. & Dudgeon, D. Responses of epibenthic algal assemblages to water abstraction in Hong Kong streams. Hydrobiologia 703(1), 225–237. https://doi.org/10.1007/s10750-012-1362-z (2013).Article 
CAS 

Google Scholar 
Maheshwari, K., Vashistha, J., Paulose, P. V. & Agarwal, T. Seasonal changes in phytoplankton community of lake Ramgarh, India. Int. J. Curr. Microbiol. Appl. Sci. 4(11), 318–330 (2015).CAS 

Google Scholar 
Luttenton, M. R., & Baisden, C. The relationships among disturbance, substratum size and periphyton community structure. In Advances in Algal Biology: A Commemoration of the Work of Rex Lowe 111–117. (Springer, 2006).Uehlinger, U. Spatial and temporal variability of periphyton biomass in a prealpine river (Necker, Switzerland). Arch. Fur. Hydrobiol. 123, 219–237 (1991).
Google Scholar 
Hill, W. R. Effects of light. In Algal Ecology in Freshwater Benthic Ecosystems. 121–148 (eds. Stevenson, R. J., Bothwell, M. L., Lowe, R. L.) (Academic Press, 1996).DeNichola, D. M. Periphyton responses to temperature at different ecological levels. In Algal Ecology in Freshwater Benthic Ecosystems. (eds. Stevenson, R. J., Bothwell, M. L., Lowe, R. L.) 149–181 (Academic Press, 1996).O’Reilly, C. M. Seasonal dynamics of periphyton in a large tropical lake. Hydrobiologia 553, 293–301. https://doi.org/10.1007/s10750-005-0878-x (2006).Article 

Google Scholar 
Borduqui, M. & Ferragut, C. Factors determining periphytic algae succession in a tropical hypereutrophic reservoir. Hydrobiologia 683, 109–122. https://doi.org/10.1007/s10750-011-0943-6 (2012).Article 
CAS 

Google Scholar 
De Souza, M. L., Pellegrini, B. G. & Ferragut, C. Periphytic algal community structure in relation to seasonal variation and macrophyte richness in a shallow tropical reservoir. Hydrobiologia 755, 183–196. https://doi.org/10.1007/s10750-015-2232-2 (2015).Article 

Google Scholar 
Prowse, T. D. River-ice hydrology. In Encyclopedia of Hydrological Sciences, vol. 4 (ed. Anderson, M. G.). (Wiley, 2005).Rusanov, A. G., Stanislavskaya, E. V. & Ács, É. Periphytic algal assemblages along environmental gradients in the rivers of the Lake Ladoga basin, Northwestern Russia: Implication for the water quality assessment. Hydrobiologia 695(1), 305–327 (2012).CAS 

Google Scholar 
Sofi, M. S., Hamid, A., Bhat, S. U., Rashid, I. & Kuniyal, J. C. Impact evaluation of the run-of-river hydropower projects on the water quality dynamics of the Sindh River in the Northwestern Himalayas. Environ. Monit. Assess. 194(9), 1–6 (2022).
Google Scholar 
MCCormick, P. V. Resource competition and species coexistence in freshwater algal assemblages. In Algal ecology—Freshwater Benthic Ecosystems (eds. Stevenson, R. J., Bothwell, M. L., Lowe, R. L.) 229–252 (Academic, 1996).Hillebrand, H., Worm, B. & Lotze, H. K. Marine microbenthic community structure regulated by nitrogen loading and grazing pressure. Mar. Ecol. Prog. Ser. 204, 27–38 (2000).CAS 

Google Scholar  More

During the COVID-19 pandemic, behavioral restrictions were imposed, after which various health problems were reported in many countries45,46. The pandemic has also increased food insecurity worldwide; consequently, panic buying has been observed in many countries, including Japan47. However, even in such situations, we found that diversity in local food access, ranging from self-cultivation to direct-to-consumer sales, was significantly associated with health and food security variables. Specifically, our results revealed the following five key discussion points.Urban agriculture in walkable neighborhoods bore fruit for health and food system resilience. However, the magnitude of its contribution differed depending on the type of urban agricultureThe results of this study showed that those who grew food by themselves at allotment farms and home gardens had significantly better subjective well-being and physical activity levels than those who did not. This result is in line with previous studies conducted during times free from the impact of infectious disease pandemics38,39,40. The use of direct sales was not related to subjective well-being but was significantly associated with physical activity. The reason might be that farm stand users tend to live in areas with farmland and travel to purchase fruits and vegetables at farm stands on foot or by bicycle. This result is consistent with that of a previous study demonstrating that the food environment in neighborhoods is an important component in promoting physical activity17.Our results also showed that those who grew food by themselves at allotment farms and those who purchased local foods at farm stands were significantly less anxious about the availability of fresh food both during the state of emergency and in the future than their counterparts. In contrast, home garden users showed significant differences only for the state of emergency. This result might be due to the differences in the size and yield of cultivation at allotment farms and home gardens. One lot in allotment farms in Tokyo can produce as much as or more than the average annual vegetable consumption per household in Japan48. However, home gardens are generally smaller and produce limited fresh foods for consumption, which may have influenced food security concerns.As in other countries, Japan imports much food from overseas and is deeply integrated into the large-scale global food system. However, as shown in this study, urban agriculture in Japanese suburbs forms small-scale, decentralized, and community-based local food systems. This multilayered food system can complement the disruptions and shortages of the global system when various problems occur for climatic, sociopolitical, or other reasons, such as pandemics. In fact, our empirical evidence suggests that urban agriculture in walkable neighborhoods, particularly allotment farms and direct-to-consumer sales at farm stands, contributed to the mitigation of food security concerns in neighborhood communities. This means that urban agriculture could enhance the resilience of the urban food system at a time when the global food system has been disrupted due to a pandemic. This validates recent discussions about the potential of urban agriculture to facilitate food system resilience10. Furthermore, our findings imply that the types of urban agriculture employed matter in determining the degree of contribution to food system resilience.To summarize the overall results, urban agriculture in walkable neighborhoods bore fruit for health and food system resilience during the COVID-19 pandemic. However, different types of urban agriculture exhibited varying associations with health and resilience. Allotment farms were positively related to all of the following: subjective well-being, physical activity, and food security concerns, both during the state of emergency and in the future. Home gardens were positively related to subjective well-being, physical activity, and food security concerns only during the state of emergency. Farm stands were positively related to physical activity and food security concerns both during the state of emergency and in the future.These differences may be due to the characteristics of the respective spaces. It is suggested that this diversity of urban agriculture has led to different types of people benefiting from various kinds of urban agriculture. Allotment farms were found to be associated with high subjective well-being, physical activity, and food security, but they may not be feasible for those who do not have enough physical strength because users are responsible for cultivating their lots, which measure 10–30 square meters40. In contrast, home gardens can be created even by those who are not confident in their physical strength. In fact, our study showed that women and older people engaged in home gardening more than men and younger people. In addition, direct-to-consumer sales at farm stands are the easiest way to obtain local fresh foods for those who do not have the time and space for allotment farms and home gardens. The need for urban agriculture has been argued in many countries2,3. However, little attention has been paid to its scale, accessibility, and diversity. Our study suggests that it is worthwhile to create diverse food production spaces within walkable neighborhoods while considering the diversity of people who access these spaces.Compared to other urban greenery and food retailers, the benefits of urban agriculture on subjective well-being and food security could be greaterCompared to the use of other urban green spaces, including urban parks, our results indicated that self-cultivation at allotment farms and home gardens was more strongly associated with subjective well-being. Previous studies have offered limited perspectives on the differences among various types of urban green spaces33. Our study further suggests that urban parks, allotment farms, and home gardens are differently associated with human health. However, as the reason was not determined, further research is needed.Furthermore, compared to other food retailers, such as supermarkets, convenience stores, and co-op deliveries, allotment farms and farm stands were more strongly associated with less anxiety about fresh food availability in the future. The availability of local fresh foods within walkable neighborhoods might have mitigated food security concerns because residents could grow food by themselves or directly observe farmers’ production processes, which may have made the difference from purchasing at places where the food systems were not visible.Flexibility in work style might promote urban agriculture in walkable neighborhoodsThere was an association between work style—working from home—and access to local food. According to the Ministry of Health, Labor and Welfare (https://www.mhlw.go.jp/english), 52% of Tokyo office workers worked from home during the first emergency declaration. Long commute times and high train congestion rates have been a problem in Tokyo suburbs, but remote workers have gained more time at and around their homes by reducing their commute times, increasing their opportunities to access local food in their walkable neighborhoods. Those who worked from home sought outdoor activities for refreshment and exercise and used a variety of urban green spaces during the pandemic49. Allotment farms and home gardens might be used as such urban green spaces. This result is consistent with previous studies assessing the characteristics of Canadian gardeners during the COVID-19 pandemic28,30.Until now, urban planners and policymakers have rarely taken work style into account. However, the flexibility of work styles and work hours may bring new insights; for example, those who work from home may become important players in urban agriculture. It has been pointed out that cities have a large hidden potential for urban agriculture by cultivating underused lands50. Our study suggests that such underused lands could be converted into productive urban landscapes for remote workers to engage in farming or gardening in between jobs as a hobby or as a side business.Food equity might be improved by urban agriculture in walkable neighborhoodsLocal fresh food is generally considered more expensive than junk food in high-income countries, creating social issues of food inequity. Therefore, past discussions on urban agriculture and food security have focused primarily on low-income households in socioeconomically disadvantaged areas24,25,26.In contrast, our study covered people from all income groups and found no statistically significant relationship between access to local food and income. This finding might be due to two urban cultural backgrounds regarding local food in Tokyo, that is, accessibility and affordability. First, residential segregation by income levels is not noteworthy in Tokyo and people from various income brackets live mixed in the same neighborhoods51. Therefore, most urban residents living in the suburbs have geographically equitable opportunities to access local foods. Second, local foods sold at farm stands are affordable. Prices are almost the same or cheaper than buying food at food retailers. While prices increase because of middleman margins related to shipping in the wholesale market, such increases are unnecessary when selling directly to consumers at farm stands. In addition, the allotment farm lots are not expensive to rent, particularly those operated by local municipalities (Supplementary Note 1).These two backgrounds make local fresh food physically and economically accessible to consumers of all income levels, resulting in food equity. This is particularly important because the concept of food system resilience includes the equitability perspective27.The integration of urban agriculture into walkable neighborhoods is a fruitful wayWhile the current discussion on walkable neighborhoods does not emphasize urban agriculture, our evidence indicated its effectiveness. The concept of walkable neighborhoods (e.g., the 15-min city model) stresses the decarbonization benefit of limiting vehicle travel, as well as the health benefits of promoting walking and cycling13,14,15,16. In addition, our research indicated that urban agriculture in walkable neighborhoods benefited health and well-being by increasing recreational outdoor opportunities to neighborhood communities, including remote workers. It also contributed to food system resilience by providing local foods to all people, including low-income households, when the global food system was disrupted due to the pandemic. Furthermore, recent studies on urban agriculture reported the decarbonization benefit of reducing carbon footprints in food production and distribution7,8. Small-scale and community-based urban agriculture in walkable neighborhoods might especially bring this benefit because neighborhood communities travel to farms on foot or by bicycle, which means almost no emission by distribution. While urban green spaces have various health benefits32,33,34,35, urban agriculture also contributes to food system resilience as well as carbon emission reduction, which makes it unique.Urban agriculture was once considered a failure of urban planning in Japan because it symbolized uncontrolled sprawl. This is analogous to the Western view, as urban agriculture was once considered the ultimate oxymoron1. However, our empirical evidence suggests that the urban‒rural mixture at neighborhood scales is a reasonable urban form that contributes to the resilience of the urban food system and to the health and well-being of neighborhood communities. It is no longer a failure of urban planning but a legacy of urban sprawl in the current urban context.Our study showed that integrating urban agriculture into walkable neighborhoods is a fruitful way of creating healthier cities and developing more resilient urban food systems during times of uncertainty. In cities where there is no farmland in intraurban areas, it would be considered effective to utilize underused spaces such as vacant lots and rooftops as productive urban landscapes. In growing cities where urban areas are still expanding, it would be advantageous to conserve agricultural landscapes within their urban fabrics. Our study could provide referential insights and robust evidence for urban policy to integrate urban agriculture into walkable neighborhoods.This study has potential limitations, including the timing of the survey and the measurement method that was utilized. We conducted the survey between June 4 and 8, 2020, just after the end of the first declaration of a state of emergency by the Japanese government. During this period, the main cultivation activities were planting and growing, and the harvest was just beginning. This seasonal constraint may have influenced the results. Because the survey was conducted during the pandemic, we used subjective methods to measure health and well-being status. However, the results might be different using objective methods52, thus further research is necessary. In addition, a longitudinal study is needed to determine whether the trends observed in this study were specific to the emergency period or whether they will persist after the COVID-19 pandemic. More