Philippa Kaur

Data were analyzed from samples collected, processed, and published previously [21, 25, 29] and have been summarized here. The present analysis, which consisted of sequence data processing, the calculation of taxon-specific isotopic signatures, and subsequent analyses, reflects original work.Sample collection and isotope incubationTo generate experimental data, three replicate soil samples were collected from the top 10 cm of plant-free patches in four ecosystems along the C. Hart Merriam elevation gradient in Northern Arizona. From low to high elevation, these sites are located in the following environments: desert grassland (GL; 1760 m), piñon-pine juniper woodland (PJ; 2020 m), ponderosa pine forest (PP; 2344 m), and mixed conifer forest (MC; 2620 m). Soil samples were air-dried for 24 h at room temperature, homogenized, and passed through a 2 mm sieve before being stored at 4 °C for another 24 h. This produced three distinct but homogenous soil samples from each of the four ecosystems that were subject to experimental treatments. Three treatments were applied to bring soils to 70% water-holding capacity: water alone (control), water with glucose (C treatment; 1000 µg C g−1 dry soil), or water with glucose and a nitrogen source (CN treatment; [NH4]2SO4 at 100 µg N g−1 dry soil). To track growth through isotope assimilation, both 18O-enriched water (97 atom %) and 13C-enriched glucose (99 atom %) were used. In all treatments isotopically heavy samples were paired with matching “light” samples that received water with a natural abundance isotope signatures. For 18O incubations, this design resulted in three soil samples per ecosystem per treatment (across four ecosystems and three treatments, n = 36) while 13C incubations were limited to only C and CN treatments (n = 24). Previous analyses suggest that three replicates is sufficient to detect growth of 10 atom % 18O in microbial DNA with a power of 0.6 and a growth of 5 atom % 18O with a power of 0.3 (12 and 6 atom % respectively for 13C) [30]. All soils were incubated in the dark for one week. Following incubation, soils were frozen at −80 °C for one week prior to DNA extraction.Quantitative stable isotope probingThe procedure of qSIP (quantitative stable isotope probing) is described here but has been applied to these samples as previously published [17, 21, 25]. DNA extraction was performed on soils using a DNeasy PowerSoil HTP 96 Kit (MoBio Laboratories, Carlsbad, CA, USA) and following manufacturer’s protocol. Briefly, 0.25 g of soils from each sample were carefully added to deep, 96-well plates containing zirconium dioxide beads and a cell lysis solution with sodium dodecyl sulfate (SDS) and shaken for 20 min. Following cell lysis, supernatant was collected and centrifuged three times in fresh 96-well plates with reagents separating DNA from non-DNA organic and inorganic materials. Lastly, DNA samples were collected on silica filter plates, rinsed with ethanol and eluted into 100 µL of a 10 mM Tris buffer in clean 96-well plates. To quantify the degree of 18O or 13C isotope incorporation into bacterial DNA (excess atom fraction or EAF), the qSIP protocol [31] was used, though modified slightly as reported previously [21, 24, 32]. Briefly, microbial growth was quantified as the change in DNA buoyant density due to incorporation of the 18O or 13C isotopes through the method of density fractionation by adding 1 µg of DNA to 2.6 mL of saturated CsCl solution in combination with a gradient buffer (200 mM Tris, 200 mM KCL, 2 mM EDTA) in a 3.3 mL OptiSeal ultracentrifuge tube (Beckman Coulter, Fullerton, CA, USA). The solution was centrifuged to produce a gradient of increasingly labeled (heavier) DNA in an Optima Max bench top ultracentrifuge (Beckman Coulter, Brea, CA, USA) with a Beckman TLN-100 rotor (127,000 × g for 72 h) at 18 °C. Each post-incubation sample was thus converted from a continuous gradient into approximately 20 fractions (150 µL) using a modified fraction recovery system (Beckman Coulter). The density of each fraction was measured with a Reichart AR200 digital refractometer (Reichert Analytical Instruments, Depew, NY, USA). Fractions with densities between 1.640 and 1.735 g cm−3 were retained as densities outside this range generally did not contain DNA. In all retained fractions, DNA was cleaned and purified using isopropanol precipitation and the abundance of bacterial 16S rRNA gene copies was quantified with qPCR using primers specific to bacterial 16S rRNA genes (Eub 515F: AAT GAT ACG GCG ACC ACC GAG TGC CAG CMG CCG CGG TAA, 806R: CAA GCA GAA GAC GGC ATA CGA GGA CTA CVS GGG TAT CTA AT). Triplicate reactions were 8 µL consisting of 0.2 mM of each primer, 0.01 U µL−1 Phusion HotStart II Polymerase (Thermo Fisher Scientific, Waltham, MA), 1× Phusion HF buffer (Thermo Fisher Scientific), 3.0 mM MgCl2, 6% glycerol, and 200 µL of dNTPs. Reactions were performed on a CFX384 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) under the following cycling conditions: 95 °C at 1 min and 44 cycles at 95 °C (30 s), 64.5 °C (30 s), and 72 °C (1 min). Separate from qPCR, retained sample-fractions were subject to a similar amplification step of the 16S rRNA gene V4 region (515F: GTG YCA GCM GCC GCG GTA A, 806R: GGA CTA CNV GGG TWT CTA AT) in preparation for sequencing with the same reaction mix but differing cycle conditions – 95 °C for 2 min followed by 15 cycles at 95 °C (30 s), 55 °C (30 s), and 60 °C (4 min). The resulting 16S rRNA gene V4 amplicons were sequenced on a MiSeq sequencing platform (Illumina, Inc., San Diego, CA, USA). DNA sequence data and sample metadata have been deposited in the NCBI Sequence Read Archive under the project ID PRJNA521534.Sequence processing and qSIP analysisIndependently from previous publications, we processed raw sequence data of forward and reverse reads (FASTQ) within the QIIME2 environment [33] (release 2018.6) and denoised sequences within QIIME2 using the DADA2 pipeline [34]. We clustered the remaining sequences into amplicon sequence variants (ASVs, at 100% sequence identity) against the SILVA 138 database [35] using a pre-trained open-reference Naïve Bayes feature classifier [36]. We removed samples with less than 3000 sequence reads, non-bacterial lineages, and global singletons and doubletons. We converted ASV sequencing abundances in each fraction to the number of 16S rRNA gene copies per gram dry soil based on qPCR abundances and the known amount of dry soil equivalent added to the initial extraction. This allowed us to express absolute population densities, rather than relative abundances. Across all replicates, we identified 114 543 unique bacterial ASVs.We calculated the 18O and 13C excess atom fraction (EAF) for each bacterial ASV using R version 4.0.3 [37] and data.table [38] with custom scripts available at https://www.github.com/bramstone/. Negative enrichment values were corrected using previously published methods [17]. ASVs that appeared in less than two of the three replicates of an ecosystem-treatment combination (n = 3) and less than three density fractions within those two replicates were removed to avoid assigning spurious estimates of isotope enrichment to infrequent taxa. Any ASVs filtered out of one ecosystem-treatment group were allowed to be present in another if they met the frequency threshold. Applying these filtering criteria, we limited our analysis towards 3759 unique bacterial ASVs which accounted for a small proportion of the total diversity but represented 68.0% of all sequence reads, and encompassed most major bacterial groups (Supplementary Fig. 1).Analysis of life history strategies and nutrient responseAll statistical tests were conducted in R version 4.0.3 [37]. We assessed the ability of phylum-level assignment of life history strategy to predict growth in response to C and N addition, as proxied by the incorporation of heavy isotope during DNA replication [39, 40]. Phylum-level assignments (Table 1) were based on the most frequently observed behavior of lineages with a representative phylum (or subphylum) as compiled previously [23]. We averaged 18O EAF values of bacterial taxa for each treatment and ecosystem and then subtracted the values in control soils from values in C-amended soils to determine C response (∆18O EAFC) and from the 18O EAF of bacteria in CN-amended soils to determine C and N response (Δ18O EAFCN). Because an ASV must have a measurable EAF in both the control and treatment for a valid Δ18O EAF to be calculated, we were only able to resolve the nutrient response for 2044 bacterial ASVs – 1906 in response to C addition and 1427 in response to CN addition.We used Gaussian finite mixture modeling, as implemented by the mclust R package [41], to demarcate plausible multi-isotopic signatures for oligotrophs and copiotrophs. For each treatment, we calculated average per-taxon 13C and 18O EAF values. To compare both isotopes directly, we divided 18O EAF values by 0.6 based on the estimate that this value (designated as µ) represents the fraction of oxygen atoms in DNA derived from the 18O-water, rather than from 16O within available C sources [42]. Two mixture components, corresponding to oligotrophic and copiotrophic growth modes, were defined using the Mclust function using ellipsoids of equal volume and shape. We observed several microorganisms with high 18O enrichment but comparatively low 13C enrichment, potentially indicating growth following the depletion of the added glucose, and that were reasonably clustered as oligotrophs in our mixture model.We tested how frequently mixture model clustering of each microorganism’s growth (based on average 18O–13C EAF in a treatment) could predict its growth across replicates (n = 12 in each treatment—although individual). We applied the treatment-level mixture models defined above to the per-taxon isotope values in each replicate, recording when a microorganism’s life history strategy in a replicate agreed with the treatment-level cluster, and when it didn’t. We used exact binomial tests to test whether the number of “successes” (defined as a microorganism being grouped in the same life history category as its treatment-level cluster) was statistically significant. To account for type I error across all individual tests (one per ASV per treatment), we adjusted P values in each treatment using the false-discovery rate (FDR) method [43].To determine the extent that life history categorizations may be appropriately applied at finer levels of taxonomic resolution, we constructed several hierarchical linear models using the lmer function in the nlme package version 3.1-149 [44]. To condense growth information from both isotopes into a single analysis, 18O and 13C EAF values were combined into a single variable using principal components analysis separately for each treatment. Across the C and CN treatments, the first principal component (PC1) was able to explain – respectively – 86% and 91% of joint variation of 18O and 13C EAF values. In all cases, we applied PC1 as the response variable and treated taxonomy and ecosystem as random model terms to limit the potential of pseudo-replication to bias significance values. We used likelihood ratio analysis and Akaike information criterion (AIC) values to compare models where life history strategy was determined based on observed nutrient responses at different taxonomic levels (Eq. 1) against a model with the same random terms but without any life history strategy data (Eq. 2). Separate models were applied to each treatment. To reduce model overfitting, we removed families represented by fewer than three bacterial ASVs as well as phyla represented by only one order. In addition, we removed bacterial ASVs with unknown taxonomic assignments (following Morrissey et al. [21]). This limited our analysis to 1 049 ASVs in the C amendment and 984 in the CN amendment.$${{{{{rm{PC}}}}}}{1}_{{18{{{{{rm{O}}}}}} – 13{{{{{rm{C}}}}}}}}sim {{{{{rm{strategy}}}}}} + 1|{{{{{rm{phylum}}}}}}/{{{{{rm{class}}}}}}/{{{{{rm{order}}}}}}/{{{{{rm{family}}}}}}/{{{{{rm{genus}}}}}}/{{{{{rm{eco}}}}}}$$
(1)
$${{{{{rm{PC}}}}}}{1}_{{18{{{{{rm{O}}}}}} – 13{{{{{rm{C}}}}}}}}sim 1 + 1|{{{{{rm{phylum}}}}}}/{{{{{rm{class}}}}}}/{{{{{rm{order}}}}}}/{{{{{rm{family}}}}}}/{{{{{rm{genus}}}}}}/{{{{{rm{eco}}}}}}$$
(2)
Here, life history strategy was defined at each taxonomic level using the mixture models above and based on the mean 18O and 13C EAF values of each bacterial lineage (Supplemental Fig. 2). We compared these models with the no-strategy model (Eq. 2) directly using likelihood ratio testing. More

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

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

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

Roemmich, D. et al. Unabated planetary warming and its ocean structure since 2006. Nat. Clim. Change 5, 240–245 (2015).Article 
ADS 

Google Scholar 
Frölicher, T. L. et al. Dominance of the Southern Ocean in anthropogenic carbon and heat uptake in CMIP5 models. J. Clim. 28, 862–886 (2015).Article 
ADS 

Google Scholar 
Buesseler, K. O. & Boyd, P. W. Shedding light on processes that control particle export and flux attenuation in the twilight zone of the open ocean. Limnol. Oceanogr. 54, 1210–1232 (2009).Article 
ADS 
CAS 

Google Scholar 
Arteaga, L., Haentjens, N., Boss, E., Johnson, K. S. & Sarmiento, J. L. Assessment of export efficiency equations in the Southern Ocean applied to satellite-based net primary production. J. Geophys. Res.-Oceans 123, 2945–2964 (2018).Article 
ADS 

Google Scholar 
Siegel, D. A. et al. Prediction of the export and fate of global ocean net primary production: the EXPORTS science plan. Front. Marine Sc. 3, 22 (2016).Perissinotto, R. & Pakhomov, E. A. The trophic role of the tunicate Salpa thompsoni in the Antarctic marine ecosystem. J. Mar. Syst. 17, 361–374 (1998).Article 

Google Scholar 
Atkinson, A., Siegel, V., Pakhomov, E. & Rothery, P. Long-term decline in krill stock and increase in salps within the Southern Ocean. Nature 432, 100–103 (2004).Article 
ADS 
CAS 

Google Scholar 
Perissinotto, R. & Pakhomov, E. A. Contribution of salps to carbon flux of marginal ice zone of the Lazarev sea, southern ocean. Mar. Biol. 131, 25–32 (1998).Article 
CAS 

Google Scholar 
Phillips, B., Kremer, P. & Madin, L. P. Defecation by Salpa thompsoni and its contribution to vertical flux in the Southern Ocean. Mar. Biol. 156, 455–467 (2009).Article 

Google Scholar 
Stone, J. P. & Steinberg, D. K. Salp contributions to vertical carbon flux in the Sargasso Sea. Deep-Sea Res. Part I 113, 90–100 (2016).Article 
CAS 

Google Scholar 
Ramaswamya, V., Sarin, M. M. & Rengarajan, R. Enhanced export of carbon by salps during the northeast monsoon period in the northern Arabian Sea. Deep-Sea Res. Part II 52, 1922–1929 (2005).Article 
ADS 

Google Scholar 
Smith, K. L. et al. Large salp bloom export from the upper ocean and benthic community response in the abyssal northeast Pacific: day to week resolution. Limnol. Oceanogr. 59, 745–757 (2014).Article 
ADS 
CAS 

Google Scholar 
Madin, L. P. & Kremer, P. Determination of the filter feeding rates of salps (Tunicata, Thaliacea). ICES J. Mar. Sci. 52, 583–595 (1995).Article 

Google Scholar 
Wiebe, P. H., Madin, L. P., Haury, L. R., Harbison, G. R. & Philbin, L. M. Diel vertical migration by Salpa aspera and its potential for large-scale particulate organic matter transport to the deep-sea. Mar. Biol. 53, 249–255 (1979).Article 

Google Scholar 
Dadon-Pilosof, A., Lombard, F., Genin, A., Sutherland, K. R. & Yahel, G. Prey taxonomy rather than size determines salp diets. Limnol. Oceanogr. 64, 1996–2010 (2019).Article 
ADS 

Google Scholar 
Stukel, M. R., Décima, M., Selph, K. E. & Gutiérrez-Rodríguez, A. Size-specific grazing and competitive interactions between large salps and protistan grazers. Limnol. Oceanogr. 66, 2521–2534 (2021).Madin, L. P. Production, composition and sedimentation of salp fecal pellets in oceanic waters. Mar. Biol. 67, 39–45 (1982).Article 

Google Scholar 
Michaels, A. F. & Silver, M. W. Primary production, sinking fluxes and the microbial food web. Deep-Sea Res. Part I 35, 473–490 (1988).Article 
ADS 

Google Scholar 
Luo, J. Y. et al. Gelatinous zooplankton‐mediated carbon flows in the global oceans: a data‐driven modeling study. Glob. Biogeochem. Cycle 34, e2020GB006704 (2020).Kremer, P. & Madin, L. P. Particle retention efficiency of salps. J. Plankton Res. 14, 1009–1015 (1992).Article 

Google Scholar 
Harbison, G. R. & Gilmer, R. W. Feeding rates of pelagic tunicate Pegea confederata and 2 other salps. Limnol. Oceanogr. 21, 517–528 (1976).Article 
ADS 
CAS 

Google Scholar 
Harbison, G. R. & McAlister, V. L. The filter-feeding rates and particle retention efficiencies of 3 species of Cyclosalpa (Tunicata, Thaliacea). Limnol. Oceanogr. 24, 875–892 (1979).Article 
ADS 

Google Scholar 
Mullin, M. M. In situ measurement of filtering rates of the salp Thalia democratica, on phytoplankton and bacteria. J. Plankton Res. 5, 279–288 (1983).Article 

Google Scholar 
Deibel, D. Clearance rates of the salp Thalia democratica fed naturally occurring particles. Mar. Biol. 86, 47–54 (1985).Article 

Google Scholar 
Fender, C. K. et al. Prey size spectra and predator:prey ratios of 7 species of New Zealand salps. Mar. Biol. (in press).Chiswell, S. M., Bostock, H. C., Sutton, P. J. H. & Williams, M. J. M. Physical oceanography of the deep seas around New Zealand: a review. N.Z. J. Mar. Freshw. Res. 49, 286–317 (2015).Article 

Google Scholar 
Henschke, N. et al. Salp-falls in the Tasman Sea: a major food input to deep-sea benthos. Mar. Ecol. Prog. Ser. 491, 165–175 (2013).Article 
ADS 

Google Scholar 
Childerhouse, S., Dix, B. & Gales, N. Diet of new Zealand sea lions (Phocarctos hookeri) at the Auckland islands. Wildl. Res. 28, 291–298 (2001).Article 

Google Scholar 
Horn, P. L., Burrell, T., Connell, A. & Dunn, M. R. A comparison of the diets of silver (Seriolella punctata) and white (Seriolella caerulea) warehou. Mar. Biol. Res. 7, 576–591 (2011).Article 

Google Scholar 
Horn, P. L., Forman, J. S. & Dunn, M. R. Dietary partitioning by two sympatric fish species, red cod (Pseudophycis bachus) and sea perch (Helicolenus percoides), on Chatham Rise, New Zealand. Mar. Biol. Res. 8, 624–634 (2012).Article 

Google Scholar 
Forman, J. S., Horn, P. L. & Stevens, D. W. Diets of deepwater Oreos (Oreosomatidae) and orange roughy Hoplostethus atlanticus. J. Fish. Biol. 88, 2275–2302 (2016).Article 
CAS 

Google Scholar 
Carroll, E. L. et al. Multi-locus DNA metabarcoding of zooplankton communities and scat reveal trophic interactions of a generalist predator. Sci. Rep. 9, 1–14 (2019).Savoye, N. et al. 234Th sorption and export models in the water column: a review. Mar. Chem. 100, 234–249 (2006).Article 
CAS 

Google Scholar 
Sutton, P. J. H. The Southland Current: a subantarctic current. N.Z. J. Mar. Freshw. Res. 37, 645–652 (2003).Article 

Google Scholar 
Foxton, P. The distribution and life history of Salpa thompsoni Foxton with observations on a related species, Salpa gerlachei Foxton. Discovery Rep. 34, 1–116 (1966).Loeb, V. J. & Santora, J. A. Population dynamics of Salpa thompsoni near the Antarctic Peninsula: growth rates and interannual variations in reproductive activity (1993–2009). Prog. Oceanogr. 96, 93–107 (2012).Article 
ADS 

Google Scholar 
Pakhomov, E. A. & Hunt, B. P. V. Trans-Atlantic variability in ecology of the pelagic tunicate Salpa thompsoni near the Antarctic Polar Front. Deep-Sea Res. Part II 138, 126–140 (2017).Article 

Google Scholar 
Lüskow, F., Pakhomov, E. A., Stukel, M. R. & Décima, M. Biology of Salpa thompsoni at the Chatham Rise, New Zealand: demography, growth, and diel vertical migration. Mar. Biol. 167, 1–18 (2020).Pakhomov, E. A. & Froneman, P. W. Zooplankton dynamics in the eastern Atlantic sector of the Southern Ocean during the austral summer 1997/1998—Part 2: grazing impact. Deep-Sea Res. Part II 51, 2617–2631 (2004).Article 
ADS 

Google Scholar 
Iversen, M. H. et al. Sinkers or floaters? Contribution from salp pellets to the export flux during a large bloom event in the Southern. Ocean. Deep-Sea Res. Part II 138, 116–125 (2017).Article 
CAS 

Google Scholar 
Buesseler, K. O., Boyd, P. W., Black, E. E. & Siegel, D. A. Metrics that matter for assessing the ocean biological carbon pump. Proc. Natl Acad. Sci. USA 117, 9679 (2020).Article 
ADS 
CAS 

Google Scholar 
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).Article 

Google Scholar 
Hall, J., Safi, K. & Cumming, A. Role of microzooplankton grazers in the subtropical and subantarctic waters to the east of New Zealand. N.Z. J. Mar. Freshw. Res. 38, 91–101 (2004).Article 

Google Scholar 
Zeldis, J. R. & Décima, M. Mesozooplankton connect the microbial food web to higher trophic levels and vertical export in the New Zealand Subtropical Convergence Zone. Deep-Sea Res. Part I 155, 103146 (2020).Article 
CAS 

Google Scholar 
Hall, J. A., James, M. R. & Bradford-Grieve, J. M. Structure and dynamics of the pelagic microbial food web of the Subtropical Convergence region east of New Zealand. Aquat. Micro. Ecol. 20, 95–105 (1999).Article 

Google Scholar 
Bradford-Grieve, J. M. et al. Pelagic ecosystem structure and functioning in the subtropical front region east of New Zealand in austral winter and spring 1993. J. Plankton Res. 21, 405–428 (1999).Article 

Google Scholar 
Nodder, S. & Gall, M. Pigment fluxes from the Subtropical Convergence region, east of New Zealand: relationships to planktonic community structure. N.Z. J. Mar. Freshw. Res. 32, 441–465 (1998).Article 
CAS 

Google Scholar 
Nodder, S. D., Chiswell, S. M. & Northcote, L. C. Annual cycles of deep-ocean biogeochemical export fluxes in subtropical and subantarctic waters, southwest Pacific Ocean. J. Geophys. Res.: Oceans 121, 2405–2424 (2016).Article 
ADS 
CAS 

Google Scholar 
Kiko, R. et al. Biological and physical influences on marine snowfall at the equator. Nat. Geosci. 10, 852-+ (2017).Article 
ADS 
CAS 

Google Scholar 
Kelly, R. P., Shelton, A. O. & Gallego, R. Understanding PCR processes to draw meaningful conclusions from environmental DNA studies. Sci. Rep. 9, 12133 (2019).Article 
ADS 

Google Scholar 
Caron, D. A., Madin, L. P. & Cole, J. J. Composition and degradation of salp fecal pellets: implications for vertical flux in oceanic environments. J. Mar. Res. 47, 829–850 (1989).Article 
CAS 

Google Scholar 
Sempere, R., Yoro, S., Wambeke, F. V. & Charriere, B. Microbial decomposition of large organic particles in the northwestern Mediterranean Sea: an experimental approach. Mar. Ecol. Prog. Ser. 198, 61–72 (2000).Article 
ADS 

Google Scholar 
Dell’Anno, A. & Corinaldesi, C. Degradation and turnover of extracellular DNA in marine sediments: ecological and methodological considerations. Appl. Environ. Microbiol. 70, 4384–4386 (2004).Article 
ADS 

Google Scholar 
Torti, A., Lever, M. A. & Jørgensen, B. B. Origin, dynamics, and implications of extracellular DNA pools in marine sediments. Mar. Genomics 24, 185–196 (2015).Article 

Google Scholar 
Norris, R. Sediments of the Chatham Rise. N.Z. Dep. Sci. Ind. Res. Res. Bull. 159, 38 (1964).
Google Scholar 
Waite, A. M., Safi, K. A., Hall, J. A. & Nodder, S. D. Mass sedimentation of picoplankton embedded in organic aggregates. Limnol. Oceanogr. 45, 87–97 (2000).Article 
ADS 

Google Scholar 
Gafar, N. A. & Schulz, K. G. A three-dimensional niche comparison of Emiliania huxleyi and Gephyrocapsa oceanica: reconciling observations with projections. Biogeosciences 15, 3541–3560 (2018).Article 
ADS 
CAS 

Google Scholar 
Eynaud, F., Giraudeau, J., Pichon, J. J. & Pudsey, C. J. Sea-surface distribution of coccolithophores, diatoms, silicoflagellates and dinoflagellates in the South Atlantic Ocean during the late austral summer 1995. Deep-Sea Res. Part I 46, 451–482 (1999).Article 

Google Scholar 
Hagino, K., Okada, H. & Matsuoka, H. Coccolithophore assemblages and morphotypes of Emiliania huxleyi in the boundary zone between the cold Oyashio and warm Kuroshio currents off the coast of Japan. Mar. Micropaleontol. 55, 19–47 (2005).Article 
ADS 

Google Scholar 
Rhodes, L. L., Peake, B. M., Mackenzie, A. L. & Marwick, S. Coccolithophores Gephyrocapsa oceanica and Emiliana Huxleyi (Prymnesiophyceae = Haptophyceae) in New Zealand’s coastal waters: characteristics of blooms and growth in laboratory culture. N.Z. J. Mar. Freshw. Res. 29, 345–357 (1995).Article 

Google Scholar 
Ziveri, P., de Bernardi, B., Baumann, K.-H., Stoll, H. M. & Mortyn, P. G. Sinking of coccolith carbonate and potential contribution to organic carbon ballasting in the deep ocean. Deep-Sea Res. Part II 54, 659–675 (2007).Article 
ADS 

Google Scholar 
Buesseler, K. O. et al. VERTIGO (VERtical Transport In the Global Ocean): a study of particle sources and flux attenuation in the North Pacific. Deep Sea Res. II 55, 1522–1539 (2008).Article 
ADS 

Google Scholar 
Billett, D. S. M., Lampitt, R. S., Rice, A. L. & Mantoura, R. F. C. Seasonal sedimentation of phytoplankton to the deep-sea benthos. Nature 302, 520–522 (1983).Article 
ADS 
CAS 

Google Scholar 
Martin, J. H., Fitzwater, S. E., Gordon, R. M., Hunter, C. N. & Tanner, S. J. Iron, primary production and carbon nitrogen fluxes during the JGOFS North Atlantic Bloom Experiment. Deep-Sea Res. Part II 40, 115–134 (1993).Article 
ADS 
CAS 

Google Scholar 
Buesseler, K. O. et al. The effect of marginal ice-edge dynamics on production and export in the Southern Ocean along 170 degrees W. Deep-Sea Res. Part II 50, 579–603 (2003).Article 
ADS 
CAS 

Google Scholar 
Kiko, R. et al. Zooplankton-mediated fluxes in the eastern tropical north. Atl. Front. Mar. Sci. 7, 21 (2020).
Google Scholar 
Kelly, T. B. et al. The importance of mesozooplankton diel vertical migration for sustaining a mesopelagic food web. Front. Mar. Sci. 6, 508 (2019).Maiti, K., Charette, M. A., Buesseler, K. O. & Kahru, M. An inverse relationship between production and export efficiency in the Southern Ocean. Geophys. Res. Lett. 40, 1557–1561 (2013).Article 
ADS 

Google Scholar 
Loeb, V. et al. Effects of sea-ice extent and krill or salp dominance on the Antarctic food web. Nature 387, 897–900 (1997).Article 
ADS 
CAS 

Google Scholar 
Steinberg, D. K. et al. Long-term (1993–2013) changes in macrozooplankton off the Western Antarctic Peninsula. Deep-Sea Res. Part I 101, 54–70 (2015).Article 

Google Scholar 
Cabanes, D. J. E. et al. First evaluation of the role of salp fecal pellets on iron biogeochemistry. Front. Mar. Sci. 3, 10 (2017).Article 

Google Scholar 
Belcher, A. et al. Krill faecal pellets drive hidden pulses of particulate organic carbon in the marginal ice zone. Nat. Commun. 10, 1–8 (2019).Manno, C. et al. Continuous moulting by Antarctic krill drives major pulses of carbon export in the north Scotia Sea, Southern Ocean. Nat. Commun. 11, 6051 (2020).Article 
ADS 
CAS 

Google Scholar 
Law, C. S. et al. Did dilution limit the phytoplankton response to iron addition in HNLCLSi sub-Antarctic waters during the SAGE experiment? Deep-Sea Res. Part II 58, 786–799 (2011).Article 
ADS 
CAS 

Google Scholar 
Gutiérrez‐Rodríguez, A. et al. Decoupling between phytoplankton growth and microzooplankton grazing enhances productivity in Subantarctic waters on Campbell Plateau, southeast of New Zealand. J. Geophys. Res.: Oceans 125, e2019JC015550 (2020).Sherman, J., Gorbunov, M. Y., Schofield, O. & Falkowski, P. G. Photosynthetic energy conversion efficiency in the West Antarctic Peninsula. Limnol. Oceanogr. 65, 1–14 (2020).Peterson, B. J. Aquatic primary productivity and the 14C-CO2 method: a history of the productivity problem. Annu Rev. Ecol. Syst. 11, 359–385 (1980).Article 

Google Scholar 
Landry, M. R. & Hassett, R. P. Estimating the grazing impact of marine microzooplankton. Mar. Biol. 67, 283–288 (1982).Article 

Google Scholar 
Landry, M. R., Haas, L. W. & Fagerness, V. L. Dynamics of microbial plankton communities—experiments in Kaneohe Bay, Hawaii. Mar. Ecol. Prog. Ser. 16, 127–133 (1984).Article 
ADS 
CAS 

Google Scholar 
Gutierrez-Rodriguez, A., Latasa, M., Estrada, M., Vidal, M. & Marrase, C. Carbon fluxes through major phytoplankton groups during the spring bloom and post-bloom in the Northwestern Mediterranean Sea. Deep Sea Res. Part I 57, 486–500 (2010).Article 
CAS 

Google Scholar 
Gutierrez-Rodriguez, A. & Latasa, M. Pigment-based measurements of phytoplankton rates. in Phytoplankton Pigments Characterization, Chemotaxonomy and Applications in Oceanography (eds Roy, S. et al.) (Cambridge University Press, 2011) 472–495.Lorenzen, C. J. Determination of chlorophyll and pheo-pigments: spectrophotometric equations. Limnol. Oceanogr. 12, 343–346 (1967).Article 
ADS 
CAS 

Google Scholar 
Conover, R. J., Durvasula, R., Roy, S. & Wang, R. Probable loss of chlorophyll-derived pigments during passage through the gut of zooplankton, and some of the consequences. Limnol. Oceanogr. 31, 878–887 (1986).Article 
ADS 
CAS 

Google Scholar 
Latasa, M. A simple method to increase sensitivity for RP-HPLC phytoplankton pigment analysis. Limnol. Oceanogr. Meth 12, 46–53 (2014).Article 

Google Scholar 
Caporaso, J. G., Paszkiewicz, K., Field, D., Knight, R. & Gilbert, J. A. The Western English Channel contains a persistent microbial seed bank. ISME J. 6, 1089–1093 (2012).Article 
CAS 

Google Scholar 
Callahan, B. J. et al. DADA2: high-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581-+ (2016).Article 
CAS 

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

Google Scholar 
McMurdie, P. J. & Holmes, S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8, e61217 (2013).Oksanen, J. et al. vegan: community ecology package. R package version 2.5-6. (2019).Piredda, R. et al. Diversity and temporal patterns of planktonic protist assemblages at a Mediterranean Long Term Ecological Research site. FEMS Microbiol. Ecol. 93, fiw200 (2017).Guillou, L. et al. The Protist Ribosomal Reference database (PR2): a catalog of unicellular eukaryote Small Sub-Unit rRNA sequences with curated taxonomy. Nucleic Acids Res. 41, D597–D604 (2013).Article 
CAS 

Google Scholar 
Pike, S. M., Buesseler, K. O., Andrews, J. & Savoye, N. Quantification of Th-234 recovery in small volume sea water samples by inductively coupled plasma-mass spectrometry. J. Radioanal. Nucl. Chem. 263, 355–360 (2005).Article 
CAS 

Google Scholar 
Benitez-Nelson, C. R. et al. Testing a new small-volume technique for determining Th-234 in seawater. J. Radioanal. Nucl. Chem. 248, 795–799 (2001).Article 
CAS 

Google Scholar 
Bone, Q. The Biology of Pelagic Tunicates (Oxford University Press, 1998).Foxton, P. An aid to the detailed examination of salps [Tunicata: Salpidae]. J. Mar. Biol. Assoc. UK 45, 679–681 (1965).Article 

Google Scholar 
Thompson, H. Pelagic Tunicates of Australia (Commonwealth Council for Scientific and Industrial Research, 1948).Iguchi, N. & Ikeda, T. Metabolism and elemental composition of aggregate and solitary forms of Salpa thompsoni (Tunicata: Thaliacea) in waters off the Antarctic Peninsula during austral summer 1999. J. Plankton Res. 26, 1025–1037 (2004).Article 
CAS 

Google Scholar 
von Harbou, L. et al. Salps in the Lazarev Sea, Southern Ocean: I. feeding dynamics. Mar. Biol. 158, 2009–2026 (2011).Article 

Google Scholar 
Pakhomov, E. A., Dubischar, C. D., Strass, V., Brichta, M. & Bathmann, U. V. The tunicate Salpa thompsoni ecology in the Southern Ocean. I. Distribution, biomass, demography and feeding ecophysiology. Mar. Biol. 149, 609–623 (2006).Article 

Google Scholar 
Huntley, M. E., Sykes, P. F. & Marin, V. Biometry and trophodynamics of Salp thompsoni Foxton (Tunicata, Thaliacea) near the Antarctic peninsula in austral summer 1983–1984. Polar Biol. 10, 59–70 (1989).Article 

Google Scholar 
Knauer, G. A., Martin, J. H. & Bruland, K. W. Fluxes of particulate carbon, nitrogen, and phosphorus in the upper water column of the northeast Pacific. Deep Sea Res. Part I 26, 97–108 (1979).Article 
ADS 
CAS 

Google Scholar 
Karl, D. M. et al. Seasonal and interannual variability in primary production and particle flux at Station ALOHA. Deep-Sea Res. Part II 43, 539–568 (1996).Article 
ADS 
CAS 

Google Scholar  More

Anna Lena Bercht interviewed fishers in Lofoten, Norway, to assess how climate change was affecting their livelihoods.Credit: Anna Lena Bercht

I remember February 2011, when, in the Chinese megacity of Guangzhou, an older man finally overcame his scepticism about being interviewed and invited me to sit down next to him on a stone bench under a shady tree. I held my notebook on my lap, and we sat on either side of a translator and talked about his life and world for more than two hours. It was one of the most informative and revealing interviews that I had done during my fieldwork in the city.
Making it in the megacity
One of the most fundamental challenges in qualitative fieldwork is gaining access to research participants. This is often time-consuming and labour-intensive, particularly when the topic requires in-depth methods and addresses a sensitive subject.Advice that goes beyond the usual recommendations of establishing relationships with gatekeepers, ensuring anonymity for interviewees and relying on the snowball sampling technique (in which one research participant suggests further ones) is rare. In this light, I’m happy to share some simple, but often neglected, examples from my qualitative fieldwork in the lively Guangzhou (where I worked for 12 months)1 and on the remote, Arctic island chain of Lofoten, Norway (done over 4 months)2, that might offer some inspiration and encouragement.I have a background in human geography, and did my PhD on experiences of stress, coping and resilience among the Chinese population of Guangzhou in the face of the city’s rapid urbanization. I travelled there five times to help to establish research cooperation with Chinese scholars, make field observations, select a case-study site and interview locals. I, together with other PhD students, stayed in a typical Chinese high-rise apartment in a neighbourhood that wasn’t a common choice for expatriates. Living side-by-side with the locals gave us a perfect opportunity to experience genuine everyday life and Chinese culture.My first postdoctoral project after my PhD brought me to Lofoten, where I looked at psychological barriers to climate adaptation in small-scale coastal fisheries. I went to Lofoten twice. On my first visit, I travelled across the whole archipelago by bus for one month to get a profound overview of the fishing villages and local living conditions, and to conduct first interviews. During my second visit, I stayed for a total of three months in rental locations near fishing harbours, and conducted more extensive interviews.In both China and Norway, I used in-depth interviews to learn about the challenges that people face. I asked people about unemployment, about the possibility of being forced to move elsewhere and about how climate change might affect their livelihoods. This required a sensitive and thoughtful approach to ‘getting invited’ into people’s lives. In Guangzhou, German- and English-speaking Chinese students assisted me as translators (and interpreters, when needed). On Lofoten, I conducted the interviews myself in English.There are two ways to access research participants: physical access, which refers to the ability of the researcher to get in direct face-to-face contact with people, and mental access. Successful mental access means that interlocutors open up about why they think, feel and behave as they do. Physical access is a necessary condition for mental access; however, in my experience, both are equally valuable.

Chinese interviewees in Guangzou shared their feelings about the rapid urbanization of their city.Credit: Anna Lena Bercht

Compared with Lofoten, it took longer to get physical access to local inhabitants in China. Presumably, this was because of the language barrier and reliance on translators, as well as cultural differences. Trust is considered a central tenet in Chinese relationships, and time and effort are needed to let it grow. During my time in Guangzhou, I occasionally benefited from being a foreigner: people were touched that someone from abroad showed genuine interest in their well-being. In Lofoten, fishers appreciated talking to a social scientist instead of a natural scientist who would have mainly asked questions about fishing quotas and catch volume.My advice for other social scientists hoping to gain access to research participants falls into those two categories.How to get good physical accessUse local public transport. Using local public transport creates many unexpected opportunities to bump into people, get into conversations and gain relevant information. For example, while waiting at a bus stop in Lofoten, I came across an art-gallery owner from a fishing village. He wondered why I was travelling out of the peak tourism season. I ended up with an invitation to his gallery, where he introduced me to two retired fishers whom he had also invited. Without the gallerist and his proactive networking, I probably would not have been given the chance to interview these two very informative and engaging fishers.In a metro station in Guangzhou, a toddler kept staring at me and tried to touch my light hair. This small interaction led me to chat to the toddler’s father, who recommended that I talk to a local teacher to learn more about the area’s history. His advice opened up important insights into urban-restructuring processes that I would have missed otherwise.
Nine ‘brain food’ tips for researchers
Use local media. In Norway, a journalist was at the harbour to get first-hand information on the year’s cod catch, when he saw me interviewing fishers. He became curious and eager to learn more about my work. In the end, he wrote an article about my research, which was published a few days later across Lofoten. His article was a door-opener for me.People recognized me from my photo in the article and contacted me to tell me about their lives and the cod fisheries. They also invited me on their vessels and put me in touch with other key informants.Change your workplace. During fieldwork, a workplace is often needed for interview transcription, literature research and interim data analysis. Moving the workplace outside wherever you are staying during a field trip allows you to immerse yourself in the daily lives of local people and interact with them more easily. For me, such agile ‘mini-office’ locations were cafes, public libraries and picnic tables. In this way, I was able to recruit interview partners on the spot.How to create deeper mental accessWear appropriate outfits. First impressions count, always. Researchers are judged not only on what they say and how they say it, but also on how they look. Certain clothes, such as those with a political slogan or religious symbol, have certain meanings and connotations. Depending on the context and whom you talk to, your appearance could promote or impede making connections and building rapport. For instance, whereas my practical ‘outdoorsy’ get-dirty outfit was appropriate for interviews on fishing vessels, a modest appearance (non-branded clothes and a simple style) was useful in rural areas of Guangzhou.Show respect. Just like in any other relationship, respect and humility play a crucial part in building a trustworthy interviewer–interviewee relationship. Showing respect can be subtly embedded in conversations in many ways, including in the content of questions and the manner in which they are asked. When interviewees started to close down when asked about painful issues, such as underemployment or loss of identity, I upheld their privacy, comfort and security by not probing when given an evasive answer. Instead, I changed the interview focus and, when appropriate, cautiously reapproached the sensitive issue by using interview techniques such as roleplaying. Interviewees were asked to put themselves in the position of someone else, such as a spatial planner or politician, and assess the issue at hand from this perspective. Taking such an imaginary role can help to make the interviewees feel more secure and face pain more openly.Be humble. Having a modest view of yourself is essential to communicate at eye level with people. As a scientist, you can easily fall into the trap of thinking that your thoughts and concepts are somehow more valuable because you are well-educated and established. However, you are the one asking questions — and the interviewees, whether they are fishers, farmers or homeless people, often know more about many things than you do. Being aware of this is an expression of humility. I let the interviewees know that they were the local experts and I was the foreign learner.Use small talk. Small talk — including non-verbal communication, such as smiling, or connective gestures, for example handing out a handkerchief or offering some tea — has an essential bonding function. Talking about ‘safe’ topics can help the interviewee to overcome the feelings of otherness, newness and discomfort that can emerge in an interview, and fosters social cohesiveness. This can help to counteract the asymmetrical power relationship between the researcher (who asks) and the researched (who answers). For example, before substantive questioning, I created shared experiences by talking about last night’s storm or the world cod-fishing championship, which takes place every year in Lofoten. This took the relationship to a greater level of intimacy and togetherness — which small talk after finishing the interview can strengthen. I remember joking about my stamina for eating properly with chopsticks to one interviewee.Use self-disclosure. Revealing selected information about yourself and sharing your own thoughts with interlocutors can help to create and reaffirm a sphere of confidentiality and trust. Fishers in Norway would, for instance, often ask “What interested you in Lofoten coastal fisheries?” or “Why do you ask me and not the scientists from Tromsø University?” I answered such questions honestly, which assisted in creating a more balanced relationship, encouraging the interviewees to address sensitive subjects more openly and readily.Change interview sites. In several interviews, I found that the answers given tended to depend on where the interview was held and which identity that site evoked for the interviewee. For example, a fisher did not talk about climate-change concerns on his fishing vessel (any concern was masked by his existential fear of losing his livelihood as a coastal fisher), but he later that day freely discussed his worries in his home. Changing the interview site can be a helpful technique to access hidden thoughts and feelings.Above all, be realistic. You will probably make mistakes; I regretted not dressing warmly enough on a fishing vessel in Arctic weather. Locals will find you amusing, weird or impolite. They will keep out of your way, and you will never know why. And they will terminate interviews prematurely with no excuse. And that’s all right. In the end, fieldwork is a combination of planning, resources, time, skills, hard work, commitment, headache, joy — and luck. Learn from your mistakes, and accept the things you cannot change. More

Arunrat, N., Sereenonchai, S., Chaowiwat, W. & Wang, C. Climate change impact on major crop yield and water footprint under CMIP6 climate projections in repeated drought and flood areas in Thailand. Sci. Total Environ. 807, 150741 (2022).ADS 
CAS 

Google Scholar 
Chandio, A. A., Shah, M. I., Sethi, N. & Mushtaq, Z. Assessing the effect of climate change and financial development on agricultural production in ASEAN-4: the role of renewable energy, institutional quality, and human capital as moderators. Environ. Sci. Pollut. Res. 29, 13211–13225 (2022).
Google Scholar 
Masood, N., Akram, R., Fatima, M., Mubeen, M., Hussain, S., Shakeel, M., Khan, N., Adnan, M., Wahid, A., Shah, A. N. and Ihsan, M. Z. (2022) Insect pest management under climate change. In Building climate resilience in agriculture. Springer, ChamOzdemir, D. The impact of climate change on agricultural productivity in Asian countries: A heterogeneous panel data approach. Environ. Sci. Pollut. Res. 29, 8205–8217 (2022).
Google Scholar 
Aidoo, O. F. et al. Climate-induced range shifts of invasive species (Diaphorina citri Kuwayama). Pest Manag. Sci. 78, 2534–2549 (2022).CAS 

Google Scholar 
Hebbar, K. B. et al. Predicting the Potential Suitable Climate for Coconut (Cocos nucifera L.) Cultivation in India under Climate Change Scenarios Using the MaxEnt Model. Plants. 11, 731 (2022).
Google Scholar 
Martín-Vélez, V. & Abellán, P. Effects of climate change on the distribution of threatened invertebrates in a Mediterranean hotspot. Insect Conserv. Divers. 15, 370–379 (2022).
Google Scholar 
Williams, J. J., Freeman, R., Spooner, F. & Newbold, T. Vertebrate population trends are influenced by interactions between land use, climatic position, habitat loss and climate change. Glob. Chang. Biol. 28, 797–815 (2022).CAS 

Google Scholar 
Aidoo, O. F. et al. Lethal yellowing disease: insights from predicting potential distribution under different climate change scenarios. J. Plant Dis. Prot. 128, 1313–1325 (2021).
Google Scholar 
Sofaer, H. R. et al. Development and delivery of species distribution models to inform decision-making. Bioscience 69, 544–557 (2019).
Google Scholar 
Mead FW, The Asiatic citrus psyllid, Diaphorina citri Kuwayama (Homoptera: Psyllidae). Florida Department of Agriculture Conservation Service, Division of Plant Industry Entomological Circular No. 180.Bové, J. M. Huanglongbing: A destructive, newly-emerging, century-old disease of citrus. Plant Pathol. J. 1, 7–37 (2006).
Google Scholar 
Li, S., Wu, F., Duan, Y., Singerman, A. & Guan, Z. Citrus greening: Management strategies and their economic impact. HortScience 55, 604–612 (2020).
Google Scholar 
Jia, H. et al. Genome editing of the disease susceptibility gene Cs LOB 1 in citrus confers resistance to citrus canker. Plant Biotechnol. J. 15, 817–823 (2017).CAS 

Google Scholar 
Ehsani, R., Dewdney, M. & Johnson, E. Controlling HLB with thermotherapy: What have we learned so far?. Citrus Ind. News 9, 26–28 (2016).
Google Scholar 
Spreen, T. H., Baldwin, J. P. & Futch, S. H. An economic assessment of the impact of Huanglongbing on citrus tree plantings in Florida. J. Hortic. Sci. 49, 1052–1055 (2014).
Google Scholar 
Djeddour, D., Pratt, C., Constantine, K., Rwomushana, I. and Day, R., (2021) The Asian citrus greening disease (Huanglongbing). Evidence note on invasiveness and potential economic impacts for East Africa. CABI Working Paper, 24, 94Hu, J., Jiang, J. & Wang, N. Control of citrus Huanglongbing via trunk injection of plant defense activators and antibiotics. Phytopathology 108, 186–195 (2018).CAS 

Google Scholar 
Fan, G. C. et al. Evaluation of thermotherapy against Huanglongbing (citrus greening) in the greenhouse. J. Integr. Agric. 15, 111–119 (2016).
Google Scholar 
Nguyen, V. A., Bartels, D. & Gilligan, C. Modelling the spread and mitigation of an emerging vector-borne pathogen: citrus greening in the US. Biorxiv https://doi.org/10.1101/2022.05.04.490566 (2022).Article 

Google Scholar 
Milosavljević, I. et al. Post-release evaluation of Diaphorencyrtus aligarhensis (Hymenoptera: Encyrtidae) and Tamarixia radiata (Hymenoptera: Eulophidae) for biological control of Diaphorina citri (Hemiptera: Liviidae) in Urban California, USA. Agronomy 12, 583 (2022).
Google Scholar 
Maluta, N., Castro, T. & Lopes, J. R. S. Entomopathogenic fungus disrupts the phloem-probing behavior of Diaphorina citri and may be an important biological control tool in citrus. Sci. Rep. 12, 1–10 (2022).
Google Scholar 
Hall, D. G., Richardson, M. L., Ammar, E. D. & Halbert, S. E. Asian citrus psyllid, Diaphorina citri, vector of citrus huanglongbing disease. Entomol. Exp. Appl. 146, 207–223 (2013).
Google Scholar 
Vázquez-García, M. et al. Insecticide resistance in adult Diaphorina citri Kuwayama1 from lime orchards in central west Mexico. Southwest. Entomol. 38, 579–596 (2013).
Google Scholar 
Naeem, A., Freed, S., Jin, F. L., Akmal, M. & Mehmood, M. Monitoring of insecticide resistance in Diaphorina citri Kuwayama (Hemiptera: Psyllidae) from citrus groves of Punjab Pakistan. Crop Prot. 86, 62–68 (2016).CAS 

Google Scholar 
Hulme, P. E. et al. Grasping at the routes of biological invasions: A framework for integrating pathways into policy. J. Appl. Ecol. 45, 403–414 (2008).
Google Scholar 
Oke, A. O., Oladigbolu, A. A., Kunta, M., Alabi, O. J. & Sétamou, M. First report of the occurrence of Asian citrus psyllid Diaphorina citri (Hemiptera: Liviidae), an invasive species in Nigeria. West Africa. Sci. Rep. 10, 1–8 (2020).
Google Scholar 
Tang, Y.Q. (1990) On the parasite complex of Diaphorina citri Kuwayama (Homoptera: Psyllidae) in Asian-Pacific and other areas. In proceedings 4th international conference on citrus rehabilitation, Chiang Mai, Thailand. 4: 240 245Chien, C. C., Chiu, S. C. & Ku, S. C. Biological control of Diaphorina citri in Taiwan. Fruits 44, 401–407 (1989).
Google Scholar 
Hoddle, M. S. Foreign exploration for natural enemies of Asian citrus psyllid, Diaphorina citri (Hemiptera: Psyllidae), in the Punjab of Pakistan for use in a classical biological control program in California USA. Pakistan Entomol. 34, 1–5 (2012).
Google Scholar 
Étienne, J., Quilici, S., Marival, D., Franck, A. & Gonzalez Fernandez, C. Biological control of Diaphorina citri (Hemiptera: Psyllidae) in Guadeloupe by imported Tamarixia radiata (Hymenoptera: Eulophidae). Fruits 56, 307–315 (2001).
Google Scholar 
Qureshi, J. A., Rogers, M. E., Hall, D. G. & Stansly, P. A. Incidence of invasive Diaphorina citri (Hemiptera: Psyllidae) and its introduced parasitoid Tamarixia radiata (Hymenoptera: Eulophidae) in Florida citrus. J. Econ. Entomol. 102, 247–256 (2009).
Google Scholar 
Chen, X., Triana, M. & Stansly, P. A. Optimizing production of Tamarixia radiata (Hymenoptera: Eulophidae), a parasitoid of the citrus greening disease vector Diaphorina citri (Hemiptera: Psylloidea). Biol. Control. 105, 13–18. https://doi.org/10.1016/j.biocontrol.2016.10.010 (2017).Article 

Google Scholar 
Kistner, E. J., Amrich, R., Castillo, M., Strode, V. & Hoddle, M. S. Phenology of Asian citrus psyllid (Hemiptera: Liviidae), with special reference to biological control by Tamarixia radiata, in the residential landscape of southern California. J. Econ. Entomol. 109, 1047–1057. https://doi.org/10.1093/jee/tow021 (2016).Article 

Google Scholar 
Ramos Aguila, L. C. et al. Temperature-dependent biological control effectiveness of Tamarixia radiata (Hymenoptera: Eulophidea) under laboratory conditions. J. Econ. Entomol. 114, 2009–2017 (2021).
Google Scholar 
Ramos Aguila, L. C. et al. Temperature-dependent demography and population projection of Tamarixia radiata (Hymenoptera: Eulophidea) reared on Diaphorina citri (Hemiptera: Liviidae). J. Econ. Entomol. 113, 55–63 (2020).
Google Scholar 
Ashraf, H. J. et al. Comparative microbiome analysis of Diaphorina citri and its associated parasitoids Tamarixia radiata and Diaphorencyrtus aligarhensis reveals Wolbachia as a dominant endosymbiont. Environ. Microbiol. 24, 1638–1652 (2022).CAS 

Google Scholar 
Chow, A. & Sétamou, M. Parasitism of Diaphorina citri (Hemiptera: Liviidae) by Tamarixia radiata (Hymenoptera: Eulophidae) on residential citrus in Texas: Importance of colony size and instar composition. Biol. Control 165, 104796 (2022).
Google Scholar 
Ajene, I. J. et al. Habitat suitability and distribution potential of Liberibacter species (“Candidatus Liberibacter asiaticus” and “Candidatus Liberibacter africanus”) associated with citrus greening disease. Environ. Microbiol. 26, 575–588 (2020).
Google Scholar 
Shabani, F., Kumar, L. & Ahmadi, M. A comparison of absolute performance of different correlative and mechanistic species distribution models in an independent area. Ecol. Evol. 6, 5973–5986 (2016).
Google Scholar 
Kearney, M. & Porter, W. Mechanistic niche modelling: Combining physiological and spatial data to predict species’ ranges. Ecol 12, 334–350 (2009).
Google Scholar 
Byeon, D. H., Jung, S. & Lee, W. H. Review of CLIMEX and MaxEnt for studying species distribution in South Korea. J. Asia-Pac. Biodivers. 1, 325–333 (2018).
Google Scholar 
Kriticos, D. J., Yonow, T. & McFadyen, R. E. The potential distribution of Chromolaena odorata (Siam weed) in relation to climate. Weed Res 45, 246–254 (2005).
Google Scholar 
Wharton, T. N. & Kriticos, D. J. The fundamental and realized niche of the Monterey pine aphid, Essigella californica (Essig) (Hemiptera: Aphididae): implications for managing softwood plantations in Australia. Divers. Distrib. 10, 253–262 (2004).
Google Scholar 
Sutherst, R., Maywald, G. and Kriticos, D., CLIMEX version 3: user’s guide. (2007).Ramirez-Cabral, N. Y., Kumar, L. & Shabani, F. Global alterations in areas of suitability for maize production from climate change and using a mechanistic species distribution model (CLIMEX). Sci. Rep. 7, 1–3 (2017).CAS 

Google Scholar 
McCalla, K. A., Keçeci, M., Milosavljević, I., Ratkowsky, D. A. & Hoddle, M. S. The influence of temperature variation on life history parameters and thermal performance curves of Tamarixia radiata (Hymenoptera: Eulophidae), a parasitoid of the Asian citrus psyllid (Hemiptera: Liviidae). J. Econ. Entomol. 112, 1560–1574 (2019).
Google Scholar 
Gonzalez-Cabrera, J., Moreno-Carrillo, G., Sanchez-Gonzalez, J. A. & Bernal, H. C. Natural and augmented parasitism of tamarixia radiata (Hymenoptera Eulophidae) in Urban Areas of western Mexico. Entomol. Sci. 53, 486–492. https://doi.org/10.18474/JES17-112.1 (2018).Article 

Google Scholar 
Chavez, Y. et al. Tamarixia radiata (Waterston) and Cheilomenes sexmaculata (Fabricius) as biological control agents of Diaphorina citri Kuwayama in Ecuador. Chil. J. Agric. Res. 77, 180–184. https://doi.org/10.4067/S0718-58392017000200180 (2017).Article 

Google Scholar 
Flores, D. & Ciomperlik, M. Biological control using the ectoparasitoid, Tamarixia radiata, against the Asian citrus psyllid, Diaphorina citri, in the lower Rio Grande valley of Texas. Southwest. Entomol. 42, 49–59. https://doi.org/10.3958/059.042.0105 (2017).Article 

Google Scholar 
Parra, J. R., Alves, G. R., Diniz, A. J. & Vieira, J. M. Tamarixia radiata (Hymenoptera: Eulophidae) × Diaphorina citri (Hemiptera: Liviidae): Mass rearing and potential use of the parasitoid in Brazil. J. Integr. Pest. Manag. https://doi.org/10.1093/jipm/pmw003 (2016).Article 

Google Scholar 
Diniz, A. J. F., Otimização da criação de Diaphorina citri Kuwayama, 1908 (Hemiptera: Liviidae) e de Tamarixia radiata (Waterston, 1922) (Hymenoptera: Eulophidae), visando a produção em larga escala do parasitoide e avalliação do seu estabelecimento em campo. Tese (Doutorado em Entomologia)—Escola Superior de Agricultura “Luiz de Queiroz”, Universidade de São Paulo, São Paulo. (2013)Hoddle, M. S. & Pandey, R. Host range testing of Tamarixia radiata (Hymenoptera: Eulophidae) sourced from the Punjab of Pakistan for classical biological control of Diaphorina citri (Hemiptera: Liviidae: Euphyllurinae: Diaphorinini) in California. J. Econ. Entomol. 107, 125–136. https://doi.org/10.1603/EC13318 (2014).Article 

Google Scholar 
Gómez-Torres, M. L., Nava, D. E. & Parra, J. R. Thermal hygrometric requirements for the rearing and release of Tamarixia radiata (Waterston) (Hymenoptera, Eulophidae). Rev. Bras. Entomol. 58, 291–295. https://doi.org/10.1590/S0085-56262014000300011 (2014).Article 

Google Scholar 
Gómez-Torres, M. L., Nava, D. E. & Parra, J. R. Life table of Tamarixia radiata (Hymenoptera: Eulophidae) on Diaphorina citri (Hemiptera: Psyllidae) at different temperatures. J. Econ. Entomol. 105, 338–343 (2012).
Google Scholar 
Chong, J. H., Roda, A. L. & Mannion, C. M. Density and natural enemies of the Asian Citrus Psyllid, Diaphorina citri (Hemiptera: Psyllidae), in the residential landscape of Southern Florida. J. Agric. Urban Entomol. 27, 33–49. https://doi.org/10.3954/11-05.1 (2010).Article 

Google Scholar 
Pluke, R. W., Qureshi, J. A. & Stansly, P. A. Citrus flushing patterns, Diaphorina citri (Hemiptera: Psyllidae) populations and parasitism by Tamarixia radiata (Hymenoptera: Eulophidae) in Puerto Rico. Florida Entomol. 91, 36–42 (2008).
Google Scholar 
Ashraf, H. J. et al. Genetic diversity of Tamarixia radiata populations and their associated endosymbiont Wolbachia species from China. Agronomy 11, 2018 (2021).CAS 

Google Scholar 
Jung, J. M., Lee, W. H. & Jung, S. Insect distribution in response to climate change based on a model: Review of function and use of CLIMEX. Entomol. Res. 46, 223–235 (2016).
Google Scholar 
Kriticos, D. J. et al. CLIMEX Version 4, 184p (2015).
Google Scholar 
Gomez-Marco, F., Gebiola, M., Baker, B. G., Stouthamer, R. & Simmons, G. S. Impact of the temperature on the phenology of Diaphorina citri (Hemiptera: Liviidae) and on the establishment of Tamarixia radiata (Hymenoptera: Eulophidae) in urban areas in the lower Colorado Desert in Arizona. Environ. Entomol. 48, 514–523 (2019).
Google Scholar 
Vieira, J. M. Biologia em temperaturas alternantes e exigências térmicas de Diaphorina citri Kuwayama, 1908 (Hemiptera: Liviidae) e Tamarixia radiata (Waterston, 1922) (Hymenoptera: Eulophidae) visando ao seu zoneamento em regiões citrícolas do estado (Doctoral dissertation, Universidade de São Paulo).Castillo, J., Jacas, J. A., Peña, J. E., Ulmer, B. J. & Hall, D. G. Effect of temperature on life history of Quadrastichus haitiensis (Hymenoptera: Eulophidae), an endoparasitoid of Diaprepes abbreviatus (Coleoptera: Curculionidae). Biol. Control. 36, 189–196 (2006).
Google Scholar 
McFarland, C. D. & Hoy, M. A. Survival of Diaphorina citri (Homoptera: Psyllidae), and its two parasitoids, Tamarixia radiata (Hymenoptera: Eulophidae) and Diaphorencyrtus aligarhensis (Hymenoptera: Encyrtidae), under different relative humidities and temperature regimes. Fla. Entomol. 84, 227–233 (2001).
Google Scholar 
Fauvergue, X. & Quilici, S. Etude de certains parametres de la biologie de Tamarixia radiata (Waterston, 1992)(Hymenoptera: Eulophidae), ectoparasitoide primaire de Diaphorina citri Kuwayama (Hemiptera: Psyllidae) vecteur du greening des agrumes. Paris Fruits 46, 179–179 (1991).
Google Scholar 
Araújo, F. H. et al. Modelling climate suitability for Striga asiatica, a potential invasive weed of cereal crops. Crop Prot. 1(160), 106050 (2022).
Google Scholar 
Silva, D. A. & RS, Kumar L, Shabani F and Picanço MC,. Potential risk levels of invasive Neoleucinodes elegantalis (small tomato borer) in areas optimal for open-field Solanum lycopersicum (tomato) cultivation in the present and under predicted climate change. Pest Manag. Sci 73, 616–627 (2017).
Google Scholar 
Kumar, S., Neven, L. G. & Yee, W. L. Evaluating correlative and mechanistic niche models for assessing the risk of pest establishment. Ecosphere 5, 1–23. https://doi.org/10.1890/ES14-00050.1 (2014).Article 
CAS 

Google Scholar 
Kriticos, D. J. et al. CliMond: global high-resolution historical and future scenario climate surfaces for bioclimatic modelling. Methods Ecol. Evol. 1, 53–64 (2012).
Google Scholar 
Santana Júnior PA, Worldwide spatial distribution of Tuta absoluta (Lepidoptera: Gelechiidae) and its natural enemies under current and future climatic change conditions through modelling. 136 f 2019 (Tese (Doutorado em Fitotecnia) – Universidade Federal de Viçosa, 2019).
Google Scholar 
Kriticos, D. J., Maywald, G. F., Yonow, T., Zurcher, E. J., Herrmann, N. I. and Sutherst, R. W., CLIMEX Version 4: Exploring the effects of climate on plants, animals and diseases. CSIRO, Canberra.156, (2015)Ramos Aguila, L. C. et al. Temperature-dependent demography and population projection of Tamarixia radiata (Hymenoptera: Eulophidea) reared on Diaphorina citri (Hemiptera: Liviidae). J. Econ. Entomol. 113, 55–63 (2019).
Google Scholar 
Oliveira, R. C., Modelagem de nicho ecológico para Helicoverpa punctigera (Wallengren, 1860) (Lepidoptera: Noctuidae) no mundo: Potencial invasão e riscos diante das mudanças climáticas. (2021). http://www.repositorio.ufc.br/handle/riufc/61961Bazzocchi, G. G., Lanzoni, A., Burgio, G. & Fiacconi, M. R. Effects of temperature and host on the pre-imaginal development of the parasitoid Diglyphus isaea (Hymenoptera: Eulophidae). Biol. Control 26, 74–82 (2003).
Google Scholar 
Hondo, T., Koike, A. & Sugimoto, T. Comparison of thermal tolerance of seven native species of parasitoids (Hymenoptera: Eulophidae) as biological control agents against Liriomyza trifolii (Diptera: Agromyzidae) in Japan. Appl. Entomol. Zool. 41, 73–82 (2006).
Google Scholar 
Duale, A. Effect of temperature and relative humidity on the biology of the stem borer parasitoid Pediobius furvus (Gahan) (Hymenoptera: Eulophidae) for the management of stem borers. Environ. Entomol. 34, 1–5 (2005).
Google Scholar 
Ashraf, H. J. et al. Comparative transcriptome analysis of Tamarixia radiata (Hymenoptera: Eulophidae) reveals differentially expressed genes upon heat shock. Comp. Biochem. Physiol. D: Genom. Proteom. 41, 100940 (2022).CAS 

Google Scholar 
van Doan, C. et al. Natural enemies of herbivores maintain their biological control potential under short-term exposure to future CO2, temperature, and precipitation patterns. Ecol. Evol. 11, 4182–4192 (2021).
Google Scholar 
Thomson, L. J., Macfadyen, S. & Hoffmann, A. A. Predicting the effects of climate change on natural enemies of agricultural pests. Biol. Control. 52, 296–306 (2010).
Google Scholar 
Rosenblatt, A. E. & Schmitz, O. J. Climate change, nutrition, and bottom-up and top-down food web processes. Trends Ecol. Evol. 31, 965–975 (2016).
Google Scholar 
Aidoo, O. F. et al. A machine learning algorithm-based approach (MaxEnt) for predicting invasive potential of Trioza erytreae on a global scale. Ecol. Inform. 71, 101792 (2022).
Google Scholar 
Aidoo, O. F. et al. The Impact of Climate Change on Potential Invasion Risk of Oryctes monoceros Worldwide. Front. Ecol. Evol. https://doi.org/10.3389/fevo.2022.895906 (2022).Article 

Google Scholar 
Hao, M. et al. Global potential distribution of Oryctes rhinoceros, as predicted by Boosted Regression Tree model. Glob. Ecol. Conserv. 1(37), e02175 (2022).
Google Scholar 
Aidoo, O. F. et al. Model-based prediction of the potential geographical distribution of the invasive coconut mite, Aceria guerreronis Keifer (Acari: Eriophyidae) based on MaxEnt. Agric. For. Entomol. 24, 390–404 (2022).
Google Scholar  More