Ecology

Test surfacesWe used smooth and hairy, homogeneous and striped (with different widths of black and white stripes) test surfaces. In order to mimic the curved zebra back, all test surfaces had a convex cylindrical shape (length: 15 cm, radius: 7 cm). In the schlieren measurements (see subsection “Schlieren imaging”), the cylinder’s horizontal long axis was perpendicular (Supplementary Fig. S1A,B) or parallel (Supplementary Fig. S1C) to the collimated horizontal light beam illuminating the target area, while the stripes were perpendicular (Supplementary Fig. S1A) or parallel (Supplementary Fig. S1B) to the cylinder’s long axis. Supplementary Tables S1 and S2 list the patterns, colours and names of the 8 smooth and 10 hairy test surfaces. The smooth surfaces were composed of cardboard squares (15 cm × 15 cm) (Supplementary Fig. S2). The hairy surfaces were composed of cattle, horse and zebra hides glued by dextrin to a cylindrical (length: 15 cm, radius: 7 cm) gypsum base (Supplementary Fig. S3). Surfaces hsc1(7b7w)perp, hsc1(8b7w)par, hsc3(3b2w)perp and hsgc3(3s2L)perp were used to model that the black stripes of zebras could be separately erected, while the white remained flat7. No animals were killed, the horse and cattle hides were provided by Hungarian horse and cattle keepers, while the zebra hide was obtained from a Hungarian zoological garden.Thermography of lamplit test surfacesUsing a thermocamera (VarioCAM, Jenoptik Laser Optik Systeme GmbH, Jena, Germany, nominal precision of ± 1.5 K, with relative pixel-to-pixel precision  1 s. Since in the upper rectangular window of the filtered schlieren images the lifespan t of local minima of I(x) was very short, we performed statistics only for the lower rectangular window.For the statistical analysis of the characteristics of I(x) (Nmin, ΔI, dave) and air stream behaviour (t, dcovered, v, dmax, dse), we applied Principal Component Analysis (PCA) and fitted ellipses to component scores with 95% confidence interval. We applied Wilcoxon rank sum test with Bonferroni correction for the datasets. Because of the large sample size of the dataset Nmin, ΔI, dave—for smooth test surfaces 4475 observations per surface and for hairy test surfaces 5369 observations per surface—, the Wilcoxon rank sum test would result in highly significant differences for almost all comparisons15. Therefore, we used a Monte Carlo approach, where we randomly selected 250 (≈ 5% of the full dataset) samples, ran the Wilcoxon rank sum test with Bonferroni correction, recorded the results, repeated this 499 times (to have 500 runs) and finally the average of the results of the 500 runs was calculated. The data evaluation was made by our custom-written scripts in Python programming language. For statistical analyses we used the R statistical package 3.6.316.Disturbance caused by an artificial wind and a butterflyTo demonstrate the influence of very weak winds on the air streams above lamplit test surfaces, we blew air by a compressor with a press of 1.6 bar from 1.5 m above the following test surfaces: ss1(8b7w), hsz(8b7w)perp, hsh1(8b8w)perp and hhbc. The compressor tube ended in an air-blow gun. Before the measurements, the air pressure was 4 bar in the 3-litre compressor tank, and the regulator was set to 1.6 bar. The compressor was turned off prior to recording a given sequence and during the measurement the air-blow gun was used to generate horizontal “wind” until the pressure decreased to 1 bar (atmospheric pressure). The artificial wind speed created by the compressor was measured (with accuracy  More

Rearing methodologyEristalis flower flies were reared in laboratory conditions from egg clutches laid by wild-caught females in the summer of 2019 (see Supplementary Materials for detailed rearing methodology). Only flies that emerged on the same day were used in the experiments. An artificial diapause protocol (see Supplementary Materials for detailed protocol) was used to prolong the lifespan of lab-reared flower flies, as adult Eristalis flower flies in lab colonies have shorter lifespans than adult Eristalis flies in the wild48. Once the adult flies emerged, all siblings were placed in artificial diapause in a refrigerator and fed 10% sucrose ad libitum until the experiment began. These Eristalis flower fly rearing and artificial diapause protocols are a modification of previously published protocols48,49.Osmia lignaria (n = 50; Crown Bees, Woodinville, WA, USA) and Megachile rotundata (n = 50; Watts Solitary Bees, Bothell, WA, USA) were purchased and allowed to emerge in an incubator kept at 23 °C and 65% humidity. Bumble bees (Bombus impatiens) used as C. bombi source colonies or as uninfected sources of bees for the dose–response trials were purchased from Biobest (Biobest, Leamington, Ontario, Canada) and maintained in the lab by feeding sucrose and pollen from a mixture of honey bee-collected poly-floral pollen (Bee Pollen Granules, CC Pollen High Desert, Phoenix AZ, USA). To ensure the commercial colonies were free of parasites, we pulled 20 workers and screened them for parasites via microscopy. No parasites were found in any of the colonies used for the dose–response trials.Evaluating whether the European drone fly, Eristalis tenax, is a host of Crithidia bombi
After breaking artificial diapause, the E. tenax flower flies were allowed to groom, but not feed, for one hour. Each fly was then placed abdomen-first into a 1.5 mL microcentrifuge tube harness to collect defecation events (Supplementary Figure S5). The size of these tubes allowed the flies to feed comfortably, but the tubes were also tapered at the bottom, which prevented the flies from stepping in their feces. Holes were placed along the side of the tubes so the fly could respire. One large hole was placed on the lid of the tube so the fly could be inoculated directly with a pipette.Flies were randomly divided into treatment and control groups. E. tenax flies in a roughly 1:2 F:M sex ratio were used in both the treatment (n = 30) and control groups (n = 30), for a total of 60 replicates. The flies that emerged from the same egg clutch with this 1:2 F:M sex ratio were the only siblings that could accommodate the replicates needed for this experiment, which is why this sex ratio was used.The C. bombi inoculum was made fresh from infected B. impatens individuals the morning of the experiment using established protocols. Briefly, we dissected the gut of infected B. impatiens workers from a laboratory source colony that sustained a strain collected from wild B. impatiens workers from Massachusetts, USA (GPS coordinates: 42.363911 N, – 72.567747 W). We homogenized the bee guts in distilled water and diluted the mixture to 1280 C. bombi cells μL−1, which we then combined 1:1 with 30% sucrose solution for an inoculum of 640 cells μL−1, a standard inoculum concentration for infecting bumble bees with C. bombi35,50. Control groups were fed 5 μL of a 30% sucrose and blue dye (Butler Extract Co., Lancaster, PA, USA) that in pilot experiments was not found to influence host or parasite survival. Treatment groups were inoculated with 5 μL (3200 cells total) of C. bombi, 30% sucrose and blue dye solution. The number of cells used in the inoculum is similar to levels of C. bombi found in the feces of bumblebees with recently established infections37. Blue dye was used to better visualize when fecal events occurred and flies that did not drink the entire 5 μL inoculum were not used in the experiment.After feeding, the flies were monitored continuously until defecation occurred. As these flies recently emerged from artificial diapause and were starved pre-experiment, every hour post-inoculation the flies were fed a 30% sucrose and blue dye solution ad libitum to encourage defecation. Once a fly defecated, the feces were collected via pipette and diluted to a 10 μL solution with deionized water to observe and count parasites using Kova Glasstic slides. The fly was then placed in an individual 60 mL plastic portion cup with filter paper (Sigma–Aldrich, St Louis, MO, USA) and a 1.5 mL microcentrifuge tube feeder containing 500 μL of a 30% sucrose and blue dye solution for 10 days. Feeders and filter papers were replaced every 3 days to prevent mold growth. As C. bombi typically replicates in high numbers after 10 days in the guts of bumble bees51, both control and treatment flies were dissected and C. bombi gut counts were performed 10 days post-inoculation. Since actively swimming, and thus live, C. bombi is infective to susceptible bumble bee hosts35, only actively swimming C. bombi were counted. The fecal volume, dilution factor and counts of C. bombi were quantified for each individual fly to calculate the exact amount of C. bombi in the individual’s first defecation event.Dose–response dataCrithidia bombi inoculum was made from infected B. impatiens individuals the morning of each trial using the protocols described above, with two exceptions. First, the C. bombi strain was collected from wild B. impatiens workers from New York, USA (GPS coordinates: 42.457350, − 76.426907). Second, a range of serially diluted doses were used to inoculate uninfected B. impatiens workers. The doses were: 25,000 cells, 12,500 cells, 6250 cells, 3125 cells, 1563 cells, 781 cells, 391 cells, 195 cells, 98 cells, 49 cells, 24 cells, and 12 cells. To obtain these doses, we homogenized bee guts in distilled water and diluted the mixture to 5000 C. bombi cells μL−1 with 30% sucrose solution. Serial dilutions were then conducted with a 10% sucrose solution to ensure the same osmolarity of each inoculum.We conducted four replicate dose–response trials over a period of four weeks. Each week, five uninfected workers per dose from each of two colonies were administered 5 μL of C. bombi inoculum. The ten highest doses were administered for the first 2 weeks, and two additional doses (24 cells, and 12 cells) were added for the final 2 weeks. Inoculated bees were kept individually in vials and fed 30% sucrose ad libitum for 7 days at 23 °C and 65% humidity. After 7 days, the bees were dissected and C. bombi loads were quantified using a hemocytometer as described above. In addition, the right forewing was removed from each bee and marginal cell length was measured as a proxy for size52. In total, 220 bees were inoculated (20 replicates for each of the ten highest doses, 10 replicates for the two lowest doses).Defecation patterns on a shared floral resourceAll pollinators (O. lignaria, M. rotundata, E. arbustorum and E. tenax) were placed in individual 60 mL plastic portion cups lined with filter paper. Each pollinator received a 1.5 mL microcentrifuge tube feeder containing 500 μL of fluorescent dye via 2.5 g of fluorescent powder (Stardust Micas) dissolved into 500 mL 30% sucrose feeders to visualize fecal deposition on flowers. After 24 hours, filter papers were collected (for analysis of fecal volume and defecation frequency, see below) and a total of five, randomly selected pollinators of the same species were placed in 12 × 12 × 12″ mesh cages (Bioquip Products, Rancho Dominguez, CA, USA) containing inflorescences of similar sized Solidago dansolitum ‘Little Lemon’ goldenrod each replicate trial. Goldenrod was used in this experiment because both bees and flower flies were observed foraging on this abundant floral resource. Only pollinators with filter papers containing fluorescing defecation events were released in the mesh cages.All E. arbustorum cages (n = 10) contained 2:3 F:M sex ratios, except one cage contained a 3:2 F:M sex ratio. All E. tenax cages (n = 20) contained 3:2 F:M sex ratios, except four cages contained 2:3 F:M sex ratios. All O. lignaria cages (n = 10) contained 4:1 F:M sex ratios, except one cage contained a 3:2 F:M sex ratio. For the two fly species, sample sizes and F:M sex ratios were determined by the greatest, same-day sibling emergence. For O. lignaria, sample sizes and F:M sex ratios were determined by emergence availability. M. rotundata floral deposition data was not collected, as the F:M emergence was heavily skewed to males that did not interact with, and therefore defecate on, the flowers.After 24-hours, the pollinators were removed and the defecation events on the goldenrod from all cages were counted under a blacklight. The location of the defecation events on the goldenrod was recorded. The plant parts were divided into six categories: ‘inside’ the flower (inside the corolla), ‘outside’ the flower (surface of the corolla), on the sepal, on the bract (the leaflike structure beneath the flower), on the stem or on a leaf.Defecation frequency and fecal volumesThe diameter of the smallest and largest defecation events per filter paper was measured by a digital caliper and an average diameter was calculated from these two values for all pollinators. The average diameter of the defecation events was converted to an average volume (in μL) using a standard curve (Supplementary Figure S6; R2 = 0.99 for the calibration data). The collected fecal volumes defecated by control flies from the E. tenax inoculation experiment (see above) were compared to the average fly fecal volumes calculated here. This was done to analyze whether flies in a confined environment, where they were inoculated with C. bombi, defecated similar volumes to flies allowed to move freely in an individual cup, which the average volumes were estimated from. In addition, the number of defecation events (frequency) over a 24-hour period on the collected E. arbustorum (n = 46) and E. tenax (n = 100) filter papers were counted for each fly.Statistical analysesFor the E. tenax inoculation experiment, we evaluated the amount of C. bombi cells in the first defecation event using a negative binomial generalized linear model (GLM), with fly sex as predictor. We chose negative binomial over Poisson to account for overdispersion, which we evaluated using Pearson residuals. Significance of sex was evaluated using a likelihood ratio test (LRT).Data from the B. impatiens inoculation experiment were used to fit two dose–response curves, the first for infection probability, and the second for infection intensity among infected bees. Infection intensity was defined using the loads estimated from the hemocytometer. A bee was considered infected if the counts were nonzero. We first tested whether the dose ingested, wing length (as a proxy for body size) and the colony the bee came from affected its response. For infection probability, this was done using a GLM with log10(dose), colony, wing length and their interactions as predictors, and infection status as the Bernoulli response. For infection intensity, this was done using a linear model (LM) with the same predictors, and log10(intensity) as response, using only infected bees. Doses were log-transformed in accordance to how the experimental doses were varied, while intensities were log-transformed to achieve normality of the residuals. Significance of predictors were tested in accordance with the principle of marginality.While we found that wing length and colony were significant predictors, in practice the colony-specific response of a wild bee is unknown (since it would not have come from any of the experimental colonies), while the dependence on wing length is only useful in a size-based epidemiological model. Hence, we generated dose–response curves by marginalizing across colony and wing length. Finally, we tested whether linear relationships between the link function and log10(dose), assumed in LMs and GLMs, were sufficient to capture the shape of the dose–response curves, by fitting the data to shape-constrained additive models and then comparing AIC values53. SCAMs are generalized additive models (GAMs) on which additional constraints such as monotonicity have been imposed; being more flexible, they can better capture the shapes of the dose–response curves should the linear relationships be inadequate.We evaluated whether fecal volume depended on pollinator species and sex with a linear model (LM), fitted using weighted least squares to account for unequal variances between group (detected using Levene’s test). Since the transformation from diameter (of feces on filter paper) to volume introduced a noticeable skew to the distribution, we transformed the volume back to diameter and further performed a Box-Cox transformation to achieve normality54, which we verified using the Shapiro–Wilk and D’Agostino’s K2 test. The transformed volume was used as the response in the abovementioned linear model.For E. tenax, fecal volume was also manually collected from the 1.5 microcentrifuge tubes during the inoculations experiment. We compared the fecal volume from the two methods using a LM with method and fly sex as well as their interaction as predictors. Volumes were log-transformed to achieve normality, while the linear model was fitted using ordinary least squares since Levene’s test indicated no significant deviation from the assumption of equal variance.We evaluated whether defecation frequency depended on pollinator species and sex with a LM, again fitted using weighted least squares to account for unequal variances between groups. While two of the groups showed deviation from normality using the Shapiro–Wilk test, the deviations were only marginally significant and hence not expected to qualitatively affect the results55.Finally, we evaluated defecation patterns on goldenrod using a negative binomial GLM, with feces counts as the response, and pollinator species, plant location and their interaction as predictors. We did not use a mixed model with cage number as a random effect since there was only one count value per cage per location, so pseudo-replication was not an issue. Significance of predictors were evaluated using LRT in accordance to the principle of marginality56 (i.e., main effects were tested only when their interactions were insignificant and hence dropped). Post-hoc tests of pairwise contrasts with Tukey corrections were performed for predictors that were significant. We recognize that the principle sex ratio and its interactions with other predictors could also be included among the predictors; however, since each species had cages with predominantly one sex ratio (E. tenax 3F:2M; E. arbustorum 2F:3M; O. lignaria mix of 4F:1M and 5F:0M), this meant that species and sex ratios were highly correlated, making it impossible to separate their effects. Nonetheless, since female Eristalis flies do not provision their brood, the differences between sexes (e.g., time spent foraging on plants) may be less pronounced than in bees. More

1.Somveille, M., Rodrigues, A. S. L. & Manica, A. Why do birds migrate? A macroecological perspective. Glob. Ecol. Biogeogr. 24, 664–674 (2015).Article 

Google Scholar 
2.Van Der Graaf, S., Stahl, J., Klimkowska, A., Bakker, J. P. & Drent, R. H. Surfing on a green wave—How plant growth drives spring migration in the Barnacle Goose Branta leucopsis. Ardea 94, 567–577 (2006).
Google Scholar 
3.Shariatinajafabadi, M. et al. Migratory herbivorous waterfowl track satellite-derived green wave index. PLoS ONE 9, 1–11 (2014).Article 
CAS 

Google Scholar 
4.Bennetts, R. E. & Kitchens, W. M. Factors influencing movement probabilities of a nomadic food specialist: Proximate foraging benefits or ultimate gains from exploration?. Oikos 91, 459–467 (2000).Article 

Google Scholar 
5.Trierweiler, C. et al. A Palaearctic migratory raptor species tracks shifting prey availability within its wintering range in the Sahel. J. Anim. Ecol. 82, 107–120 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
6.Ma, Z., Cai, Y., Li, B. & Chen, J. Managing wetland habitats for waterbirds: An international perspective. Wetlands 30, 15–27 (2010).CAS 
Article 

Google Scholar 
7.Smit, I. P. J. Resources driving landscape-scale distribution patterns of grazers in an African savanna. Ecography (Cop.) 34, 67–74 (2011).Article 

Google Scholar 
8.Donnelly, J. P. et al. Synchronizing conservation to seasonal wetland hydrology and waterbird migration in semi-arid landscapes. Ecosphere 10, 1–12 (2019).Article 

Google Scholar 
9.Bennitt, E., Bonyongo, M. C. & Harris, S. Habitat selection by African buffalo (Syncerus caffer) in response to landscape-level fluctuations in water availability on two temporal scales. PLoS ONE 9, 1–14 (2014).Article 

Google Scholar 
10.Kleyheeg, E. et al. Movement patterns of a keystone waterbird species are highly predictable from landscape configuration. Mov. Ecol. 5, 1–14 (2017).Article 

Google Scholar 
11.Roshier, D. A., Doerr, V. A. J. & Doerr, E. D. Animal movement in dynamic landscapes: Interaction between behavioural strategies and resource distributions. Oecologia 156, 465–477 (2008).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Henry, D. A. W., Ament, J. M. & Cumming, G. S. Exploring the environmental drivers of waterfowl movement in arid landscapes using first-passage time analysis. Mov. Ecol. 4, 1–18 (2016).Article 

Google Scholar 
13.Cook, M. I., Call, E. M., Kobza, R., Mac Hill, S. D. & Saunders, C. J. Seasonal movements of crayfish in a fluctuating wetland: Implications for restoring wading bird populations. Freshw. Biol. 59, 1608–1621 (2014).Article 

Google Scholar 
14.Weimerskirch, H. et al. Lifetime foraging patterns of the wandering albatross: Life on the move!. J. Exp. Mar. Bio. Ecol. 450, 68–78 (2014).Article 

Google Scholar 
15.Krüger, S., Reid, T. & Amar, A. Differential range use between age classes of southern African bearded vultures Gypaetus barbatus. PLoS ONE 9, e114920 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
16.Wolfson, D. W., Fieberg, J. R. & Andersen, D. E. Juvenile Sandhill Cranes exhibit wider ranging and more exploratory movements than adults during the breeding season. Ibis 162, 556–562 (2019).Article 

Google Scholar 
17.Péron, C. & Grémillet, D. Tracking through life stages: Adult, immature and juvenile Autumn migration in a long-lived seabird. PLoS ONE 8, e72713 (2013).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
18.Hake, M., Kjellén, N. & Alerstam, T. Age-dependent migration strategy in honey buzzards Pernis apivorus tracked by satellite. Oikos 103, 385–396 (2003).Article 

Google Scholar 
19.Gschweng, M., Kalko, E. K. V., Querner, U., Fiedler, W. & Berthold, P. All across Africa: Highly individual migration routes of Eleonora’s falcon. Proc. R. Soc. B Biol. Sci. 275, 2887–2896 (2008).Article 

Google Scholar 
20.Miller, T. A. et al. Limitations and mechanisms influencing the migratory performance of soaring birds. Ibis 158, 116–134 (2016).Article 

Google Scholar 
21.Rotics, S. et al. The challenges of the first migration: movement and behaviour of juvenile vs. adult white storks with insights regarding juvenile mortality. J. Anim. Ecol. 85, 938–947 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
22.Sergio, F. et al. Individual improvements and selective mortality shape lifelong migratory performance. Nature 515, 410–413 (2014).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
23.Thorup, K. et al. Resource tracking within and across continents in long-distance bird migrants. Sci. Adv. 3, e1601360 (2017).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
24.Howison, R. A., Piersma, T., Kentie, R., Hooijmeijer, J. C. E. W. & Olff, H. Quantifying landscape-level land-use intensity patterns through radar-based remote sensing. J. Appl. Ecol. 55, 1276–1287 (2018).Article 

Google Scholar 
25.Wang, X. et al. Stochastic simulations reveal few green wave surfing populations among spring migrating herbivorous waterfowl. Nat. Commun. 10, 2187 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.McFeeters, S. K. The use of the normalized difference water index (NDWI) in the delineation of open water features. Int. J. Remote Sens. 17, 1425–1432 (1996).ADS 
Article 

Google Scholar 
27.Mcfeeters, S. K. Using the normalized difference water index (NDWI) within a geographic information system to detect swimming pools for mosquito abatement: A practical approach. Remote Sens. 5, 3544–3561 (2013).ADS 
Article 

Google Scholar 
28.Yang, X., Zhao, S., Qin, X., Zhao, N. & Liang, L. Mapping of urban surface water bodies from Sentinel-2 MSI Imagery at 10m resolution via NDWI-based image sharpening. Remote Sens. 9, 1–19 (2017).ADS 

Google Scholar 
29.Choi, C. Y. et al. Where to draw the line? Using movement data to inform protected area design and conserve mobile species. Biol. Conserv. 234, 64–71 (2019).Article 

Google Scholar 
30.Guillet, A. Distribution and conservation of the shoebill (Balaeniceps rex) in the southern Sudan. Biol. Conserv. 13, 39–49 (1978).Article 

Google Scholar 
31.BirdLife International. Species factsheet: Balaeniceps rex. https://www.birdlife.org. Accessed on Apr 14, 2020 (2020).32.Dodman, T. International single species plan for the conservation of the Shoebill Balaeniceps rex. AEWA Technical Series 51 (2013).33.Guillet, A. Aspects of the foraging behaviour of the shoebill. Ostritch J. Afr. Ornithol. 50, 252–255 (1979).Article 

Google Scholar 
34.Mullers, R. H. E. & Amar, A. Shoebill Balaeniceps rex foraging behaviour in the Bangweulu Wetlands, Zambia. Ostritch J. Afr. Ornithol. 86, 113–118 (2015).Article 

Google Scholar 
35.Roxburgh, L. & Buchanan, G. M. Revising estimates of the Shoebill (Balaeniceps rex) population size in the Bangweulu Swamp, Zambia, through a combination of aerial surveys and habitat suitability modelling. Ostrich J. Afr. Ornithol. 81, 25–30 (2010).Article 

Google Scholar 
36.John, J. R. M., Nahonyo, C. L., Lee, W. S. & Msuya, C. A. Observations on nesting of Shoebill Balaeniceps rex and Wattled Crane Bugeranus carunculatus in Malagari wetlands, western Tanzania. Afr. J. Ecol. 51, 184–187 (2013).Article 

Google Scholar 
37.Mullers, R. H. E. & Amar, A. Parental nesting behavior, chick growth and breeding success of Shoebills (Balaeniceps rex) in the Bangweulu Wetlands, Zambia. Waterbirds 38, 1–9 (2015).Article 

Google Scholar 
38.Elliott, A., Garcia, E. F. J. & Boesman, P. Shoebill (Balaeniceps rex). in Handbook of the Birds of the World (eds. del Hoyo, J., Elliott, A., Sargatal, J., Christie, D. A. & de Juana, E.) (Lynx Edicions, 2020).39.African Parks. African Parks: Unlocking the value of protected areas. African Parks Annual Report 2018. (2018).40.Möller, W. Beobachtungen zum Nahrungserwerb des Schuhschnabels (Balaeniceps rex). J. Ornithol. 123, 19–28 (1982).Article 

Google Scholar 
41.Christensen, K. D., Falk, K., Jensen, F. P. & Petersen, B. S. Abdim’s Stork Ciconia abdimii in Niger: Population size, breeding ecology and home range. Ostritch J. Afr. Ornithol. 79, 177–185 (2008).Article 

Google Scholar 
42.McCann, K. I. & Benn, G. A. Land use patterns within Wattled Crane (Bugeranus carunculatus) home ranges in an agricultural landscape in KwaZulu-Natal, South Africa. Ostritch J. Afr. Ornithol. 77, 186–194 (2006).Article 

Google Scholar 
43.El-Hacen, E.-H.M., Overdijk, O., Lok, T., Olff, H. & Piersma, T. Home Range, habitat selection, and foraging rhythm in Mauritanian Spoonbills (Platalea leucorodia balsaci): A satellite tracking study. Waterbirds 36, 277–286 (2013).Article 

Google Scholar 
44.King, D. T. et al. Winter and summer home ranges of American White Pelicans (Pelecanus erythrorhynchos) captured at loafing sites in the Southeastern United States. Waterbirds 39, 287–294 (2016).Article 

Google Scholar 
45.Shaw, A. K. Drivers of animal migration and implications in changing environments. Evol. Ecol. 30, 991–1007 (2016).Article 

Google Scholar 
46.Folmer, E. O., Olff, H. & Piersma, T. The spatial distribution of flocking foragers: Disentangling the effects of food availability, interference and conspecific attraction by means of spatial autoregressive modeling. Oikos 121, 551–561 (2012).Article 

Google Scholar 
47.Folmer, E. O. & Piersma, T. The contributions of resource availability and social forces to foraging distributions: A spatial lag modelling approach. Anim. Behav. 84, 1371–1380 (2012).Article 

Google Scholar 
48.Mendez, L. & Weimerskirch, H. Ontogeny of foraging behaviour in juvenile red-footed boobies (Sula sula). Sci. Rep. 7, 13886 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
49.Patrick, S. C. & Weimerskirch, H. Consistency pays: Sex differences and fitness consequences of behavioural specialization in a wide-ranging seabird. Biol. Lett. 10, 20140630 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
50.Patrick, S. C. & Weimerskirch, H. Personality, foraging and fitness consequences in a long lived seabird. PLoS ONE 9, e87269 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
51.Doherty, T. S. & Driscoll, D. A. Coupling movement and landscape ecology for animal conservation in production landscapes. Proc. R. Soc. B Biol. Sci. 285, 20172272 (2018).Article 

Google Scholar 
52.Riotte-lambert, L. & Weimerskirch, H. Do naive juvenile seabirds forage differently from adults?. Proc. R. Soc. B Biol. Sci. 280, 20131434 (2013).Article 

Google Scholar 
53.Buxton, L., Slater, J. & Brown, L. The breeding behaviour of the shoebill or whale-headed stork Balaeniceps rex in the Bangweulu Swamps, Zambia. Afr. J. Ecol. 16, 201–220 (1978).Article 

Google Scholar 
54.Roshier, D. A., Robertson, A. I. & Kingsford, R. T. Responses of waterbirds to flooding in an arid region of Australia and implications for conservation. Biol. Conserv. 106, 399–411 (2002).Article 

Google Scholar 
55.Chevallier, D. et al. Human activity and the drying up of rivers determine abundance and spatial distribution of Black Storks Ciconia nigra on their wintering grounds determine abundance and spatial distribution of Black Storks Ciconia nigra on their wintering grounds. Bird Study 3657, 369–380 (2010).Article 

Google Scholar 
56.Ng’onga, M., Kalaba, F. K., Mwitwa, J. & Nyimbiri, B. The interactive effects of rainfall, temperature and water level on fish yield in Lake Bangweulu fishery, Zambia. J. Therm. Biol. 84, 45–52 (2019).57.Grissac, S. D., Bartumeus, F., Cox, S. L. & Weimerskirch, H. Early-life foraging: Behavioral responses of newly fledged albatrosses to environmental conditions. Ecol. Evol. 7, 6766–6778 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
58.Bolduc, F. & Afton, A. D. Relationships between wintering waterbirds and invertebrates, sediments and hydrology of coastal marsh ponds. Waterbirds 27, 333–341 (2004).Article 

Google Scholar 
59.Ratcliffe, C. The fishery of the Lower Shire River area. Malawy Fisheries Bulletin No. 3. Fisheries Department, Malawi (1972).60.Xu, H. Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. Int. J. Remote Sens. 27, 3025–3033 (2006).ADS 
Article 

Google Scholar 
61.Tian, S., Zhang, X., Tian, J. & Sun, Q. Random forest classification of wetland landcovers from multi-sensor data in the arid region of Xinjuang, China. Remote Sens. 8, 1–14 (2016).CAS 

Google Scholar 
62.Kamweneshe, B. M. Ecology, Conservation and Management of the Black Lechwe (Kobus leche smithemani) in the Bangweulu Basin, Zambia. University of Pretoria (2000).63.BirdLife International. Important Bird Areas Factsheet: Bangweulu Swamps. https://www.birdlife.org. Accessed on Oct 14, 2020 (2020).64.Thurlow, J., Zhu, T. & Diao, X. The impact of climate variability and change on economic growth and poverty in Zambia. International Food Policy Research Institute (2009).65.Evans, D. W. Lake Bangweulu: A study of the complex and fishery. Fisheries Service Reports, Zambia (1978).66.Kolding, J. & van Zwieten, P. A. M. Relative lake level fluctuations and their influence on productivity and resilience in tropical lakes and reservoirs. Fish. Res. 115–116, 99–109 (2012).Article 

Google Scholar 
67.Howard, G. W. & Aspinwall, D. R. Aerial censuses of Shoebills, Saddlebilled Storks and Wattled Cranes at the Bangweulu Swamps and Kafue Flats, Zambia. Ostrich J. Afr. Ornithol. 55, 207–212 (1984).Article 

Google Scholar 
68.Microwave Telemetry. Microwave Telemetry Solar Argos/GPS 70g PTT. https://www.microwavetelemetry.com/. Accessed on Oct 14, 2020 (2020).69.Calenge, C. The package “adehabitat” for the R software: A tool for the analysis of space and habitat use by animals. Ecol. Model. 197, 516–519 (2006).Article 

Google Scholar 
70.Hijmans, R. J. geosphere: Spherical trignometry. https://cran.r-project.org/package=geosphere (2019).71.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. http://www.r-project.org/ (2019).72.Pebesma, E. J. & Bivand, R. S. Classes and methods for spatial data in R. https://cran.r-project.org/package=sp (2005).73.Bivand, R. S., Pebesma, E. J. & Gomez-Rubio, V. Applied Spatial Data Analysis with R. (Springer, 2013).74.E. Vermote. MOD09A1 MODIS/Terra Surface Reflectance 8-Day L3 Global 500m SIN Grid V006. NASA EOSDIS Land Processes DAAC. https://doi.org/10.5067/MODIS/MOD09A1.006 (2015).75.Hijmans, R. J. raster: Geographic data analysis and modeling. https://cran.r-project.org/package=raster (2019).76.Bivand, R. S., Keitt, T. & Rowlingson, B. rgdal: Bindings for the ‘geospatial’ abstraction library. https://cran.r-project.org/package=rgdal (2019).77.Bivand, R. S. & Rundel, C. rgeos: Interface to geometry engine – Open source (GEOS). https://cran.r-project.org/package=rgeos (2019).78.Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
79.Barton, K. MuMIn: Multi-Model Inference. https://cran.r-project.org/package=MuMIn (2019). More

The CV represents the discrete degree of trait values, that is, the size of the trait space (Fig. 6b; Supplementary Note S1); S and K are generally used to describe the shape of trait distribution (Fig. 6c,d, Supplementary Note S1). Environment filtering can force a trait to deviate from the original distribution, with characteristics of smaller CV and larger S and K values16,17. Partly consistent with our hypothesis, MAT significantly exerted positive effects on the total concentration of CV, S, and K, but weaker negative effects on the three values of Chl a/b (Fig. 6a). That is, the distributions of Chl concentration and composition shifted in opposite directions under global warming: Chl concentration was distributed in a broader but more differential way (Fig. 7a), while Chl a/b was distributed in a narrower but more uniform way (Fig. 7b).Figure 7Theoretical sketches of distribution shifts for (a) chlorophyll concentration and (b) composition (Chl a/b) under global warming. Purple curves denote the current distributions, and pink ones represent the scenario under global warming. Dashed lines denote the normal distribution in the respective scenarios. It is supposed that the distribution of chlorophyll concentration will shift toward higher mean, CV, S and K values, while Chl a/b shifts toward higher mean but lower CV, S and K values under warming. Chl chlorophyll, CV coefficient of variation, S skewness, K kurtosis.Full size imageThe trait distributional shift under warming is possibly caused by the relative role of species turnover and intraspecific variation (due to plasticity and/or heritable differentiation)25. For Chl concentration and composition, very weak phylogenetic signals were found in three plateaus (Supplementary Table S2), indicating the phenotypic plasticity of Chl, which environments have influenced during the long-term evolution. However, plasticity and intraspecific variation are not the focus of the discussion. Because the species compositions were significantly different among the three plateaus: with only a few species overlapping (Supplementary Fig. S3), and the dominant species and co-existing species gradually varied along the 30 sites (Supplementary Table S3). Shifts in Chl distributions under warming may be interpreted mainly by the alternation of species composition.For Chl concentration, a broader trait space (higher CV) and a more skewed distribution (higher S and K) under warming conditions indicate several new species that differ in functions (here refer to rare species with higher Chl concentration) appeared or increased. This contributed to the long tail of the curve and raised the average Chl concentration. At the same time, most of the other species converged at lower Chl concentrations; that is, Chl concentration undergoes more substantial differentiation and functional contrasting species co-exist under warming. The concentration of Chl is representative of plant growth rate and production ability. Its distribution shift may imply a possible trend of polarisation in functions: both acquisitive and conservative species occur simultaneously. This alteration in species composition indicates changing biotic interactions26. The co-existence of functional contrasting species allows individuals to avoid competition and enhance the exploitation of resources and niche27,28, which is of great importance in optimising community functions28,29. In desert and alpine regions, functional contrasting species with large inter-specific trait variations improve community multi-functionality and enable better resistance to climate change17,30.However, despite the shift in species composition, the distribution of Chl a/b only changed slightly compared to the Chl concentration under warming. The ratio of Chl a to Chl b represents the plant allocation to RC and LHC in PS and the efficiency trade-offs between light capture and light conversion6,7. This ratio is characteristic of conservatism which is mainly manifested in the following aspects: (1) Chl a/b is independent of Chl concentration (orthogonal relationship of the two; Supplementary Fig. S2); (2) Chl a/b distributed more converged with higher K and lower CV (Supplementary Table S1); (3) relative fixed allometric relationships were found between Chl a and Chl b (beside TP; Fig. 8). Plants may adjust their RC and LHC allocation to a common ratio of 3:1 despite large variations in light availability or Chl concentration, which has also been confirmed by a study from forests14. Considering that RC is costlier than LHC, plants tend to sustain the Chl a/b as low as possible unless there is a functional imbalance caused by environmental stress such as warming9,31.Figure 8Standardised major axis regression of chlorophyll a to chlorophyll b in three plateaus. Slopes were given and compared among regions; different lowercase words denote significant (p  More

This study uses SincNet according to the instructions provided by the authors for its application in a different dataset32. This section provides an introduction to SincNet and NIPS4Bplus before detailing the experimental procedure.SincNetThe first convolutional layer of a standard CNN trained on the raw waveform learns filters from the data, where each filter has a number of parameters that matches the filter length (Eq. 1).$$yleft[ n right] = xleft[ n right] times fleft[ n right] = mathop sum limits_{i = 0}^{I – 1} xleft[ i right] cdot fleft[ {n – i} right],$$
(1)

where (xleft[ n right]) is the chunk of the sound, (fleft[ n right]) is the filter of length (I), and (yleft[ n right]) is the filtered output. All the elements of the filter ((i)) are learnable parameters. SincNet replaces (fleft[ n right]) with another function (g) that only depends on two parameters per filter: the lower and upper frequencies of a rectangular bandpass filter (Eq. 2).$$gleft[ {n,f_{l} ,f_{h} } right] = 2f_{h} sincleft( {2pi f_{h} n} right) – 2f_{l} sincleft( {2pi f_{l} n} right),$$
(2)

where (f_{l} text{ and } f_{h}) are the learnable parameters corresponding to the low and high frequencies of the filter and (sincleft( x right) = frac{sinleft( x right)}{x}). The function (g) is smoothed with a Hamming window and the learnable parameters are initialised with given cut-off frequencies in the interval (left[ {0,frac{{f_{s} }}{2}} right]), where (f_{s}) is the sampling frequency.This first layer of SincNet performs the sinc-based convolutions for a set number and length of filters, over chunks of the raw waveform of given window size and overlap. A conventional CNN architecture follows the first layer, that in this study maintains the architecture and uses both standard and enhanced settings. The standard settings used are those of the TIMIT speaker recognition experiment27,32. They include two convolutional layers after the first layer with 60 filters of length 5. All three convolutions use layer normalisation. Next, three fully-connected (leaky ReLU) layers with 2048 neurons each follow, normalised with batch normalisation. To obtain frame-level classification, a final softmax output layer, using LogSoftmax, provides a set of posterior probabilities over the target classes. The classification for each file derives from averaging the frame predictions and voting for the class that maximises the average posterior. Training uses the RMSprop optimiser with the learning rate set to 0.001 and minibatches of size 128. A sample of sinc-based filters generated during this study shows their response both in the time and the frequency domains (Fig. 4).Figure 4Examples of learned SincNet filters. The top row (a–c) shows the filters in the time domain, the bottom row (d–f) shows their respective frequency response.Full size imageThe SincNet repository32 provides an alternative set of settings used in the Librispeech speaker recognition experiment27. Tests of the alternative settings, which include changes in the hidden CNN layers, provided similar results to those of the TIMIT settings and are included as Supplementary Information 1.NIPS4BplusNIPS4Bplus includes two parts: sound files and rich labels. The sound files are the training files of the 2013 NIPS4B challenge for bird song classification23. They are a single channel with a 44.1 kHz sampling rate and 32 bit depth. They comprise field recordings collected from central and southern France and southern Spain15. There are 687 individual files with lengths from 1 to 5 s for a total length of 48 min. The tags in NIPS4Bplus are based on the labels released with the 2013 Bird Challenge but annotated in detail by an experienced bird watcher using dedicated software15. The rich labels include the name of the species, the class of sound, the starting time and the duration of each sound event for each file. The species include 51 birds, 1 amphibian and 9 insects. For birds there can be two types of vocalisations: call and song; and there is also the drumming of a woodpecker. Calls are generally short sounds with simple patterns, while songs are usually longer with greater complexity and can have modular structures or produced by one of the sexes8,13. In the dataset, only bird species have more than one type of sound, with a maximum of two types. The labels in NIPS4Bplus use the same 87 tags present in the 2013 Bird Challenge training dataset with the addition of two other tags: “human” and “unknown” (for human sounds and calls which could not be identified). Tagged sound events in the labels typically correspond to individual syllables although in some occasions the reviewer included multiple syllables into single larger events15. The tags cover only 569 files of the original training set of 687 files. Files without tags include 100 that, for the purpose of the challenge, had no bird sounds but only background noise. Other files were excluded for different reasons such as vocalisations hard to identify or containing no bird or only insect sounds15. The 2013 Bird Challenge also includes a testing dataset with no labels that we did not use15.The total number of individual animal sounds tagged in the NIPS4Bplus labels is 5478. These correspond to 61 species and 87 classes (Fig. 5). The mean length of each tagged sound ranges from ~ 30 ms for Sylcan_call (the call of Sylvia cantillans, subalpine warbler) to more than 4.5 s for Plasab_song (the song of Platycleis sabulosa, sand bush-cricket). The total recording length for a species ranges from 0.7 s for Turphi_call (the call of Turdus philomelos, song thrush) to 51.4 s for Plasab_song. The number of individual files for each call type varies greatly from 9 for Cicatr_song (the call of Cicadatra atra, black cicada) to 282 for Sylcan_call.Figure 5Distribution of sound types by number of calls (number of files) and total length in seconds. Sound types are sorted first by taxonomic group and then by alphabetical order.Full size imageProcessing NIPS4BplusThe recommended pre-processing of human speech files for speaker recognition using SincNet includes the elimination of silent leading and trailing sections and the normalisation of the amplitude27. This study attempts to replicate this by extracting each individual sound as a new file according to the tags provided in the NIPS4Bplus labels. A Python script42 uses the content of the labels to read each wavefile, apply normalisation, select the time of origin and length specified in each individual tag and save it as a new wavefile. The name of the new file includes the original file name and a sequential number suffix according to the order in which tags are listed in the label files (the start time of the sound) to match the corresponding call tags at the time of processing. Each wavefile in the new set fully contains a sound according to the NIPS4Bplus labels. A cropped file may contain sounds from more than one species15, with over 20% of the files in the new set overlapping, at least in part, with sound from another species. The machine learning task does not use files containing background noise or the other parts of the files that are not tagged in the NIPS4Bplus labels. A separate Python script42 generates the lists of files and tags that SincNet requires for processing. The script randomly generates a 75:25 split into lists of train and test files and a Python dictionary key that assigns each file to the corresponding tag according to the file name. The script selects only files confirmed as animal sounds (excluding the tags “unknown” and “human”) and generates three different combinations of tags, as follows: (1) “All classes”: includes all the 87 types of tags originally included in the 2013 Bird Challenge training dataset; (2) “Bird classes”: excludes tags for insects and one amphibian species for a total of 77 classes; and (3) “Bird species”: one class for each bird species independently of the sounds type (call, songs and drumming are merged for each species) for a total of 51 classes. The script also excludes three very short files (length shorter than 10 ms) which could not be processed without code modifications.To facilitate the repeatability of the results, this study attempts to maintain the default parameters of SincNet used in the TIMIT speaker identification task27,32. The number and length of filters in the first sinc-based convolutional layer was set to the same values as the TIMIT experiment (80 filters of length 251 samples) as was the architecture of the CNN. The filters were initialised following the Mel scale cut-off frequencies. We did change the following parameters: (1) reduced the window frame size (cw_len) from 200 to 10 ms to accommodate for the short duration of some of the sounds in the NIPS4Bplus tags (such as some bird vocalisations); (2) reduced the window shift (cw_shift) from 10 to 1 ms in proportion to the reduction in window size (a value a 0.5 could not be given without code modifications); (3) updated the sampling frequency setting (fs) from the TIMIT 16,000 to the 44,100 Hz of the present dataset; and (4) updated the number of output classes (class_lay) to match the number of classes in each training run.To evaluate performance, the training sequence was repeated with the same settings and different random train and test file splits. Five training runs took place for each of the selection of tags: “All classes”, “Bird classes” and “Bird species”.Enhancements and comparisonsChanges in the parameters of SincNet result in different levels of performance. To assess possible improvements and provide baselines to compare against other models we attempted to improve the performance by adjusting a series of parameters, but did not modify the number of layers or make functional changes to the code other than the two outlined below. The parameters tested include: the length of the window frame size, the number and length of the filters in the first layer, number of filters and lengths of the other convolutional and fully connected layers, the length and types of normalisation in the normalisation layers, alternative activation and classification functions, and the inclusion of dropouts (Supplementary Information 1). In addition the SincNet code includes a hard-coded random amplification of each sound sequence; we also tested changing the level and excluding this random amplification through changes in the code. In order to process window frames larger than some of the labelled calls in the NIPS4Bplus dataset, the procedure outlined earlier in which files are cut according to the labels was replaced by a purpose-built process. The original files were not cut, instead a custom python script42 generated train and test file lists that contain the start and length of each labelled call. A modification of the SincNet code42 uses these lists to read the original files and select the labelled call. When the call is shorter than the window frame the code randomly includes the surrounding part of the file to complete the length of the window frame. Grid searches for individual parameters or combinations of similar parameters, over a set number of epochs, selected the best performing values. We also tested the use of the Additive Margin Softmax (AM-softmax) as a cost function37. The best performing models reported in the results use combinations of the best parameter values (Supplementary Information 1). All enhancements and model comparisons use the same dataset selection, that is the same train and test dataset split, of the normalised files for each set of tagged classes.The comparison using waveform + CNN models trained directly on the raw waveform, replaces the initial sinc-based convolution of SincNet with a standard 1d convolutional layer27, thus retaining the same network architecture as SincNet. As with SincNet enhancements, a series of parameter searches provided the best parameter combinations to obtain the best performing models.The pre-trained models used for comparison are DenseNet121, ResNet50 and VGG16 with architectures and weights sourced from the Torchvision library of PyTorch33. We tested three types of spectrograms: Fast Fourier Transform (FFT), Mel spectrum (Mel) and Mel-frequency cepstral coefficient (MFCC) to fine-tune the pre-trained models. FFT calculations used a frame length of 1024 samples, 896 samples overlap and a Hamming window. Mel spectrogram calculations used 128 Mel bands. Once normalised and scaled to 255 pixel intensity three repeats of the same spectrogram represented each of the three input channels of the pre-trained models. The length of sound used to generate the spectrograms was 3 s, and similarly with routines above, for labelled calls shorter than 3 s the spectrogram would randomly include the surrounding sounds. That is, the extract would randomly start in the interval between the end of the labelled call minus 3 s and the start of the call plus 3 s. This wholly includes the labelled call but its position is random within the 3 s sample. A fully connected layer replaced the final classifying layer of the pre-trained models to output the number of labelled classes. In the fine-tuning process the number of trainable layers of the model was not limited to the final fully connected layer, but also included an adjustable number of final layers to improve the results. The learning rate set initially to 0.0001 was halved if the validation loss stopped decreasing for 10 epochs.MetricsMeasures of performance include accuracy, ROC AUC, precision, recall, F1 score, top 3 accuracy and top 5 accuracy. Accuracy, calculated as part of the testing routine, is the ratio between the number of correctly predicted files of the test set and the total number of test files. The calculation of the other metrics uses the Scikit-learn module43 relying on the predicted values provided by the model and performing weighted averages. The ROC AUC calculation uses the mean of the posterior probabilities provided by SincNet for each tagged call. In the pre-trained models the ROC AUC calculations used the probabilities obtained after normalising the output with a softmax function. More

Sample collection and filtrationWe used 40 water samples for eDNA metabarcoding from 27 sites in 9 rivers and 13 lakes in Japan from 2016 to 2018 (Fig. 4). Sampling ID and detailed information for each site are listed in Supplementary Table S1. In the river water sampling, 1-L water samples were collected from the surface of at the shore of each river using bleached plastic bottles. In the field, a 1-ml Benzalkonium chloride solution (BAC, Osvan S, Nihon Pharmaceutical, Tokyo, Japan)33 was added to each water sample to suppress eDNA degeneration before filtering the water samples. We did not include field negative control samples in the HTS library, considering the aim of the presents study. The lake samples were provided by Doi et al. (2020)34 as DNA extracted samples. In the lake samples, 1-L water samples were collected from the surface at shore sites at each lake. The samples were then transported to the laboratory in a cooler at 4 °C. Each of the 1-L water samples was filtered through GF/F glass fiber filter (normal pore size = 0.7 μm; diameter = 47 mm; GE Healthcare Japan Corporation, Tokyo, Japan) and divided into two parts (maximum 500-ml water per 1 GF/F filter). To prevent cross-contamination among the water samples, the filter funnels, and the measuring cups were bleached after filtration. All filtered samples were stored at -20 ℃ in the freezer until the DNA extraction step.Figure 4Sampling sites used in the present study. Blue circles and orange triangles show the locations of the river and lake samples, respectively. Detailed information on each site is listed in Supplementary Table S1. This map has been illustrated using QGIS ver.3.10 (http://www.qgis.org/en/site/) based on the Administrative Zones Data (http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_3.html) which were obtained from free download service of the National Land Numerical Information (http://nlftp.mlit.go.jp/ksj/index.html, edited by RN). There was no need of obtaining permissions for editing and publishing of map data.Full size imageDNA extraction and library preparationThe total eDNA was extracted from each filtered sample using the DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany). Extraction methods were according to Uchii et al.35, with a few modifications. A filtered sample was placed in the upper part of a Salivette tube and 440 μL of a solution containing 400 μL Buffer AL and 40 μL Proteinase K added. The tube with the filtered sample was incubated at 56 °C for 30 min. Afterward, the tube was centrifuged at 5000 × g for 3 min, and the solution at the bottom part of the tube was collected. To increase eDNA yield, 220-μL Tris–EDTA (TE) buffer was added to the filtered sample and the sample re-centrifuged at 5000 × g for 1 min. Subsequently, 400 μL of ethanol was added to the collected solution, and the mixture was transferred to a spin column. Afterward, the total eDNA was eluted in 100-μL buffer AE according to the manufacturer’s instructions. All eDNA samples were stored at -20 °C until the library preparation step.In the present study, we used a universal primer set “MiFish” for eDNA metabarcoding9. The amplicon library was prepared according to the following protocols. In the first PCR, the total reaction volume was 12 μL, containing 6.0μL 2 × KOD buffer, 2.4 μL dNTPs, 0.2 μL KOD FX Neo (TOYOBO, Osaka, Japan), 0.35 μL MiFish-U-F (5ʹ-ACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNGTCGGTAAAACTCGTGCCAGC-3ʹ), MiFish-U-R (5ʹ-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNNCATAGTGGGGTATCTAATCCCAGTTTG-3ʹ), MiFish-E-F (5ʹ-ACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNRGTTGGTAAATCTCGTGCCAGC-3ʹ) and MiFish-E-R (5ʹ-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNNGCATAGTGGGGTATCTAATCCTAGTTTG -3ʹ) primers with Illumina sequencing primer region and 6-mer Ns, and 2 μL template DNA. The thermocycling conditions were 94 ℃ for 2 min, 35 cycles of 98 ℃ for 10 s, 65 ℃ for 30 s, 68 ℃ for 30 s, and 68 ℃ for 5 min. The first PCR was repeated four times for each sample, and the replicated samples were pooled as a single first PCR product for use in the subsequent step. The pooled first PCR products were purified using the Solid Phase Reversible Immobilization select Kit (AMPure XP; BECKMAN COULTER Life Sciences, Indianapolis, IN, USA) according to the manufacturer’s instructions. The DNA concentrations of purified first PCR products were measured using a Qubit dsDNA HS assay kit and a Qubit 3.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). All purified first PCR products were diluted to 0.1 ng/μL with H2O, and the diluted samples were used as templates for the second PCR. In the first PCR step, the PCR negative controls (four replicates) were included in each experiment. A total of three PCR negative controls were included in the library (PCR Blank 1–3 samples in Supplementary Table S1, S2, S4, and S5).The second PCR was performed to add HTS adapter sequences with 8-bp dual indices. The total reaction volume was 12 μL, containing 6.0 μL 2 × KAPA HiFi HotStart ReadyMix, 1.4 μL forward and reverse primer (2.5 μM), 1 μL purified first PCR product, and 2.2 μL H2O. The thermocycling conditions were 95 ℃ for 3 min, 12 cycles of 98 ℃ for 20 s, 72 ℃ for 15 s, and 72 ℃ for 5 min.Each Indexed second PCR product was pooled in the equivalent volume, and 25 μL of the pooled libraries were loaded on a 2% E-Gel SizeSelect agarose gels (Thermo Fisher Scientific), and a target library size (ca. 370 bp) was collected. The quality of the amplicon library was checked using an Agilent 2100 Bioanalyzer and Agilent 2100 Expert (Agilent Technologies Inc., Santa Clara, CA, USA), and the DNA concentrations of the amplicon library were measured using Qubit dsDNA HS assay Kit using a Qubit 3.0 fluorometer.High-throughput sequencingAmplicon library was sequenced using iSeq and MiSeq platforms (Illumina, San Diego, CA, USA). To normalize the percentage of pass-filtered read numbers, the sequencing runs using the same libraries were performed using iSeq i1 Reagent and MiSeq Reagent Kit v2 Micro. Both sequencing was performed with 8 million pair-end reads and 2 × 150 bp read lengths. Each library was spiked with approximately 20% PhiX control (PhiX Control Kit v3, Illumina, San Diego, CA, USA) before sequencing runs according to the recommendation of Illumina. The wells of cartridges in the iSeq run were loaded with 20 μL of 50 pM library pool, and sequencing performed at Yamaguchi University, Yamaguchi, Japan. The wells of cartridges for MiSeq runs were loaded with 600 μL of 16 pM library pool, and sequencing performed at Illumina laboratories (Minato-ku, Tokyo, Japan). Subsequently, the sequencing dataset outputs from iSeq and MiSeq were subjected to pre-processing and taxonomic assignments. All sequence data are registered in the DNA Data Bank of Japan (DDBJ) Sequence Read Archive (DRA, Accession number: DRA10593).Pre-processing and taxonomic assignmentsWe used the USEARCH v11.066736 for all data pre-processing activities and taxonomic assignment of the HTS datasets obtained from the iSeq and MiSeq platforms16,37. First, pair-end reads (R1 and R2 reads) generated from iSeq and MiSeq platforms were assembled using the “fastq_mergepairs” command with a minimum overlap of 10 bp. In the process, the low-quality tail reads with a cut-off threshold at a Phred score of 2, and the paired reads with too many mismatches ( > 5 positions) in the aligned regions were discarded38. Secondly, the primer sequences were removed from the merged reads using the “fastx_truncate” command. Afterward, read quality filtering was performed using the “fastq_filter” command with thresholds of max expected error  > 1.0 and  > 50 bp read length. The pre-processed reads were dereplicated using the “fastx_uniques” command, and the chimeric reads and less than 10 reads were removed from all samples as the potential sequence errors. Finally, an error-correction of amplicon reads, which checks and discards the PCR errors and chimeric reads, was performed using the “unoise3” command in the unoise3 algorithm39. Before the taxonomic assignment, the processed reads from the above steps were subjected to sequence similarity search using the “usearch_global” command against reference databases of fish species that had been established previously (MiFish local database v34). The sequence similarity and cut off E-value were 99% and 10–5, respectively. If there was only one species with ≧ 99% similarity, the sequence was assigned to the top-hit species. Conversely, sequences assigned to two or more species in the ≧ 99% similarity were merged as species complex and listed in the synonym group. Generally, the species complexes were assigned to the genus level (e.g., Asian crucian carp Carassius spp.). Species that were unlikely to inhabit Japan were excluded from the candidate list of species complexes. For example, the sequence of one of bitterling Acheilignathus macropterus included other different two species, A. barbatus and A. chankaensis, as the species of the 2nd hit candidate; however, the two species are not currently found in Japan. Therefore, the sequence was assigned to A. macropterus in the present study. Because we used only freshwater fish species, we removed the operational taxonomic units (OTUs) assigned to marine and brackish fishes from each sample. Finally, sequence reads of each fish species were arranged into the matrix, with the rows and columns representing the number of sites and fish species (or genus), respectively.We evaluated sequence quality based on (1) the percentage of clustering passing filter (% PF) and (2) sequencing quality score ≧ % Q30 (Read1 and Read2) between iSeq and MiSeq platforms. The % PF value is an indicator of signal purity for each cluster40. The condition leads to poor template generation, which decreases the % PF value40. In the present study, a  > 80% PF value was set as the threshold of sequence quality in iSeq and MiSeq runs. Sequence quality scores (Q score) measure the probability that a base is called incorrectly. Higher Q scores indicate lower probability of sequencing error, and lower Q scores indicate probability of false-positive variant calls resulting in inaccurate conclusions41. In the present study, the % Q30 values (error rate = 0.001%) were used for the comparison of sequence quality between iSeq and MiSeq. The parameters were collected directly using Illumina BaseSpace Sequence Hub. We also evaluated changes in sequence reads in pre-processing steps between iSeq and MiSeq platforms. Sequence reads were assessed based (1) merge pairs, (2) quality filtering, and (3) denoising. In each step, the change in the number of reads before and after processing was calculated. The calculated numbers of sequence reads are listed in Supplementary Table S2 and S3 in series.Comparing sequence quality and fish fauna between iSeq and MiSeqTo test a relationship of remained sequence reads between iSeq and MiSeq in each pre-processing part, we performed spearman’s rank correlation test in each step. In the present study, however, the sequencing run by iSeq and MiSeq was performed only once each for the same sample. Therefore, we could not assess the variabilities of the sequence read in quality checks and taxonomic assignment in the same samples between iSeq and MiSeq.Before the comparison of fish fauna, rarefaction curves were illustrated for each sample in both iSeq and MiSeq to confirm that the sequencing depth adequately covered the species composition using the “rarecurve” function of the “vegan” package ver. 2.5–6 (https://github.com/vegandevs/vegan) in R ver. 3.6.242. In the present study, the differences in the numbers of sequence reads among samples were confirmed in the two sequencers, but rarefaction curves were saturated in all iSeq and MiSeq samples (Supplementary Fig. S6 and S7). We performed a rarefaction using the “rrarefy” function in “vegan” package to match up the iSeq sequence depths of each sample with that of MiSeq. However, the number of species in each sample on the iSeq have not changed before or after the rarefaction. Therefore, we have used the raw data set before the rarefaction for the subsequent analyses.We compared the species detection capacities of iSeq and MiSeq based on environmental DNA metabarcoding. Using fish faunal data obtained from iSeq and MiSeq, non-metric multidimensional scaling (NMDS) was performed in 1000 separate runs using the “metaMDS” functions in the “vegan” package ver. 2.5–6. For NMDS, the dissimilarity of the fish fauna was calculated based on the incidence-based Jaccard indices. To evaluate the differences in species composition and variance across sites between the two HTS, we performed a permutational multivariate analysis of variance (PERMANOVA) and the permutational analyses of multivariate dispersions (PERMDISP) with 10,000 permutations, respectively. For the PERMANOVA and PERMDISP, we used the “adonis”, and “betadisper” functions in the “vegan” package ver. 2.5-6.Comparison of fish species detectability between eDNA metabarcoding and conventional methodsWe evaluated species detectability between the two HTS by comparing the fish species lists of the two HTS with lists from conventional methods. Five sampling sites were selected from Kyushu and Chugoku districts (R23–27 in Fig. 4). The fish fauna data obtained by conventional methods were based on the results of a previous study43. The conventional surveys were conducted through hand-net sampling and visual observation by snorkeling (see a previous study43 for the detailed methods). The count data of each species were replaced with the incidence-based datasets (presence or absence) for comparing with the eDNA metabarcoding datasets. Fish sequence reads of each sampling site obtained by eDNA metabarcoding were also replaced with the incidence-based data.To test the detectability of species observed by conventional methods, the fish species compositions in five rivers were compared between the eDNA metabarcoding (iSeq and MiSeq) and the conventional methods. To visualize the differences in the species composition between HTSs and conventional methods, heat maps were illustrated for each sampling site. To assess differences in the number of species among methods at each river, the repeated measures analysis of variance (ANOVA) was performed among iSeq, MiSeq, and conventional methods. If a significant difference was found in repeated measures ANOVA, the Tukey–Kramer multiple comparison test was performed to analyze differences among methods.Using fish faunal data obtained from iSeq, MiSeq, and conventional methods, the NMDS was performed in 1000 separate runs with Jaccard indices. The PERMANOVA was performed with 1000 permutations to assess the differences in fish fauna among the methods and sites. Furthermore, to evaluate variance across sites among methods, the PERMDISP was also performed with 1000 permutations. To visualize the number of species in each method and the number of common species between methods, Venn diagrams were illustrated for each river using the “VennDiagram” package ver. 1.6.2 in R44. More

Site descriptionField experiments were conducted at the Changwu Experimental Station (35° 12′ N, 107° 40′ E, altitude 1200 m) located in Shaanxi Province, China. The experimental site was in the typical dryland farming area on the Loess Plateau. Annual precipitation in the area averaged 582 mm between 1957 and 2013, with a mean annual temperature of 9.7 °C over that period. Rainfall and temperature during the two study years are shown in Fig. S1. Soils were generally of the Calcaric Regosol group, according to the FAO/UNESCO soil classification system52, and were composed of 4% sand, 59% silt, and 37% clay53. The 0–20 cm soil properties were the following: pH, 8.4; organic matter content, 11.8 g kg−1; total N content, 0.87 g kg−1; and Olsen-P, 14.4 mg kg−1.Experimental design and field managementTwo-year experiment was arranged in a randomized complete block design with three replicate plots during 2012 and 2013 growing seasons25,54_ENREF_53. The study was conducted using the soybean cultivar (Glycine max L.) cv. Zhonghuang 24 and the maize cultivar (Zea mays L.) cv. Zhengdan 958 grown in cereal–legume agricultural systems. Zhonghuang 24 was bred from Jilin 21 and fendou 31 × Zhongdou 19 (deposition number 2008003); Zhengdan 958 was the offspring of inbred Zheng 58 and Chang 7-2 (deposition number 20000009), which are approved in China. The cropping system treatments were as follows:

1.

Sole-cropped soybean (S).

2.

Sole-cropped maize (M).

3.

Two rows of maize intercropped with two rows of soybean (M2S2).

4.

Two rows of maize intercropped with four rows of soybean (M2S4).

Each plot measured 6 m × 4 m, with row spacing of 50 cm for maize and soybean both in sole crops and intercrops. Individual plants were spaced at 22 cm and 19 cm for maize and soybean, respectively, with one plant per stand for maize and two plants per stand for soybean to attain densities of 90,000 and 210,000 plants ha−1, respectively. In 2012, seeds of maize and soybean were sown on 25 April and harvested on 28 September, and in 2013, seeds were sown on 20 April and harvested on 25 September. Before sowing, basal fertilizer was applied at a rate of 90 kg N ha−1 as urea (46% N) and 150 kg P2O5 ha−1 as superphosphate (12%, P2O5), and then additional fertilizers were uniformly spread in each plot, which were then ploughed into the 0–30 cm soil layer using a rotary tiller. All of the plots received 67.5 kg N ha−1 as urea at the bell and silking stages using a hole-seeding machine. No irrigation was applied, and weeds were removed by hand when sighted. The research on plants complied with relevant institutional, national, and international guidelines and legislation.Above- and below-ground measurementsThe Pn was measured with a LI-6400 portable photosynthesis system (LI-COR Inc., Lincoln, NE, USA) from 9:00 to 11:00 h at 120 days after sowing, which corresponds to the milk stage in maize and full seed stage in soybean7,13. We measured photosynthesis of ear leaves of maize, the first spreading leaves at the top of soybean in both the sole crops and intercrops. The Pn values were calculated as the sum of the mean readings for five leaves in each plot. The LAI values, DIFN were recorded using a Plant Canopy Analyzer (Li-2200, LiCor Inc., Lincoln, NE, USA) without direct sunlight at milk stage of maize. One above-canopy measurement and three below-canopy measurements at the soil surface were taken for four replicates in each plot. SPAD were collected using a hand-held dual wavelength meter (SPAD 502, Chlorophyll meter, Minolta Camera Co., Ltd., Japan) at milk stage of maize. Measurements were taken midway along the ear leaves of maize and the first spreading leaves at the top of soybean from five adjacent plants at the center of row in each plot.The SWS was measured gravimetrically using a soil auger at 10 cm intervals over a depth of 100 cm and at 20 cm intervals over a depth of 200 cm at milk stage of maize for three replicates in each plot. The SWS was calculated for each plot in the 0–200 cm soil profile for the soil moisture using the following formula: SWS = SWC × SD × SBD, where SWC represents soil water content, SD represents soil depth, and SBD represents soil bulk density. Apparent water use during crop growth season was expressed as evapotranspiration (ET), which was determined according to the following formula: ET = ΔSWS + P, where ΔSWS is the change in soil water storage in the top 200 cm and P is the rainfall (mm) between planting and at milk stage in maize. The six adjacent plant samples were collected at milk stage of maize in the middle two rows of each plots (Fig. S2). The sampling included shoots and roots of maize and soybean. At the cotyledonary node, above-ground parts were separated from below-ground parts. Soil core samples (9 cm diameter × 15 cm) at the intra-row of crop were collected to a depth of 100 cm using an auger and separated in 10-cm sections to determine the root growth in sole-cropping and intercropping systems. The samples were exposed to 105 °C for 30 min and then dried to a constant weight at 75 °C. The oven-dried samples were put in small plastic bags after grinding. The study of N and P uptake are the most common among mineral elements55,56. Concentrations of N and P in the plant dry matter were determined after digestion with H2SO4 and H2O2; N concentration was measured according to the Kjeldahl method20, whereas P concentration was measured by the molybdenum-antimony anti-spectrophotometric method16. Crop N and P uptake were calculated by the actual above-ground biomass multiplied by plant tissue N and P concentrations. Grain yield was estimated at harvest from 6 m2 for maize and soybean based on the average of three plot replicates.Data analysisThe LER for assessment of land use advantage. LER is sum of ratio of intercrop to sole crop for maize and soybean yield57:$$ LER = LER_{m} + LER_{s} ,;LER_{m} = frac{{Y_{im} }}{{Y_{sm} }}, ;LER_{s} = frac{{Y_{is} }}{{Y_{ss} }} $$where LERm and LERs are patial LER for maize and soybean, respectively. Yim and Yis are yields of maize and soybean under intercrops, respectively. Ysm and Yss are the yield of maize and soybean under sole crop, respectively.The water equivalent ratio (WER) was calculated to measure water use advantage of intercropping58:$$ WER = WER_{m} + WER_{s} ,;WER_{m} = frac{{Y_{im} /ET_{im} }}{{Y_{sm} /ET_{sm} }},;WER_{s} = frac{{Y_{is} /ET_{is} }}{{Y_{ss} /ET_{ss} }} $$where WERm and WERs are patial WER for maize and soybean, respectively. ETim and ETis are ET of maize and soybean under intercrops, respectively. ETsm and ETss are the ET of maize and soybean under sole crop, respectively.All analyses were conducted in SPSS Statistics 17.0 (SPSS Inc., Chicago, IL, USA). Treatment means showing significant differences among different cropping systems were separated using one-way ANOVA or least significant difference (LSD) at a threshold of 5% to compare the effect of yield, above- and below-ground related parameters (Pn, LAI, SPAD, DIFN, SWS, N and P uptake) in different maize–soybean intercropping. The variation in Pn, LAI, SPAD, DIFN, SWS, N, and P uptake of crop, and the effects of cropping system × year were made using Univariate General Linear Models. Pearson’s correlation test was used to analyze between LER and above-and below-ground biomass of maize and soybean. The effects of above- and below-ground factors on biological yield were quantified, by calculating the contribution value of some key factors to yield. The effects of between above- (LAI, SPAD, DIFN) and below-ground (SWS, N and P uptake) competition on the biological yield and contribution rate were conducted by the linear regression model59:$$ Y = beta_{0} LAI + beta_{1} SPAD + beta_{2} DIFN + beta_{3} SWS + beta_{4} {text{N}} + beta_{5} {text{P}} + beta_{6} X + beta_{7} $$
(1)
where Y represents biological yield, LAI represents leaf area index, SPAD represents chlorophyll, DIFN represents diffuse non interceptance, SWS represents soil water storage, N represents crop nitrogen uptake, P represents crop phosphorus uptake, X represents interaction for LAI, SPAD, DIFN, SWS, N, and P, and β0, β1, β2, β3, β4, β5, β6 and β7 represent the fitted parameters. The standard regression coefficients (Beta) of LAI, SPAD, DIFN, SWS, N, and P were determined on the basis of Eq. (1) to split their influence on the biological yield by the following equations:$$ beta_{0}^{prime } = beta_{0} times left( {LAI^{prime } /Y^{prime } } right) $$
(2)
$$ beta_{1}^{prime } = beta_{1} times left( {SPAD^{prime } /Y^{prime } } right) $$
(3)
$$ beta_{2}^{prime } = beta_{2} times left( {DIFN^{prime } /Y^{prime } } right) $$
(4)
$$ beta_{3}^{prime } = beta_{3} times left( {SWS^{prime } /Y^{prime } } right) $$
(5)
$$ beta_{4}^{prime } = beta_{4} times left( {{text{N}}^{prime } /Y^{prime } } right) $$
(6)
$$ beta_{5}^{prime } = beta_{5} times left( {{text{P}}^{prime } /Y^{prime } } right) $$
(7)
where β0′, β1′, β2′, β3′, β4′, and β5′ represent the standard regression coefficients for LAI, SPAD, DIFN, SWS, N, and P. LAI′, SPAD′, DIFN′, SWS′, N′, and P′ represent the standard deviations for LAI, SPAD, DIFN, SWS, N, and P. Y′ is the standard deviation for the modeled biological yield. More

1.Litchman, E., Ohman, M. D. & Kiørboe, T. Trait-based approaches to zooplankton communities. J. Plankton Res. 35, 473–484. https://doi.org/10.1093/plankt/fbt019 (2013).Article 

Google Scholar 
2.Kiørboe, T. & Hirst, A. G. Shifts in mass scaling of respiration, feeding, and growth rates across life-form transitions in marine pelagic organisms. Am. Nat. 183, E118–E130. https://doi.org/10.1086/675241 (2014).Article 
PubMed 

Google Scholar 
3.Andersen, K. H. et al. Characteristic sizes of life in the oceans, from bacteria to whales. Annu. Rev. Mar. Sci. 8, 1–25. https://doi.org/10.1146/annurev-marine-122414-034144 (2015).Article 

Google Scholar 
4.Bergmann, C. Über die Verhältnisse der Wärmeökonomie der Thiere zu ihrer Grösse. Göttinger Studien 3, 595–708 (1847).
Google Scholar 
5.Woodson, C., Schramski, J. R. & Joye, S. B. A unifying theory for top-heavy ecosystem structure in the ocean. Nat. Commun. 9, 1–8. https://doi.org/10.1038/s41467-017-02450-y (2018).CAS 
Article 

Google Scholar 
6.Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. & West, G. B. Toward a metabolic theory of ecology. Ecology 85, 1771–1789. https://doi.org/10.1890/03-9000 (2004).Article 

Google Scholar 
7.Gardner, J. L., Peters, A., Kearney, M. R., Joseph, L. & Heinsohn, R. Declining body size: A third universal response to warming?. Trends Ecol. Evol. 26, 285–291. https://doi.org/10.1016/j.tree.2011.03.005 (2011).Article 
PubMed 

Google Scholar 
8.Angilletta, M. J., Steury, T. D. & Sears, M. W. Temperature, growth rate, and body size in ectotherms: Fitting pieces of a life-history puzzle. Integr. Comp. Biol. 44, 498–509. https://doi.org/10.1093/icb/44.6.498 (2004).Article 
PubMed 

Google Scholar 
9.Atkinson, D. Temperature and organism size: A biological law for ectotherms?. Adv. Ecol. Res. 25, 1–58. https://doi.org/10.1016/S0065-2504(08)60212-3 (1994).Article 

Google Scholar 
10.Atkinson, D. & Sibly, R. M. Why are organisms usually bigger in colder environments? Making sense of a life history puzzle. Trends Ecol. Evol. 12, 235–239. https://doi.org/10.1016/S0169-5347(97)01058-6 (1997).CAS 
Article 
PubMed 

Google Scholar 
11.Sunagawa, S. et al. Structure and function of the global ocean microbiome. Science 348, 1261359. https://doi.org/10.1126/science.1261359 (2015).CAS 
Article 
PubMed 

Google Scholar 
12.Audzijonyte, A. et al. Is oxygen limitation in warming waters a valid mechanism to explain decreased body sizes in aquatic ectotherms?. Glob. Ecol. Biogeogr. 28, 64–77. https://doi.org/10.1111/geb.12847 (2018).Article 

Google Scholar 
13.Begon, M., Townsend, C. R. & Harper, J. L. Ecology: From Individuals to Ecosystems 4th edn. (Blackwell Publishing, New York, 2006).
Google Scholar 
14.Hirata, T., Aiken, J., Hardman-Mountford, N., Smyth, T. J. & Barlow, R. G. An absorption model to determine phytoplankton size classes from satellite ocean colour. Remote Sens. Environ. 112, 3153–3159. https://doi.org/10.1016/j.rse.2008.03.011 (2008).ADS 
Article 

Google Scholar 
15.Kostadinov, T., Siegel, D. & Maritorena, S. Global variability of phytoplankton functional types from space: Assessment via the particle size distribution. Biogeosciences 7, 3239–3257. https://doi.org/10.5194/bg-7-3239-2010 (2010).ADS 
Article 

Google Scholar 
16.Brun, P., Payne, M. R. & Kiørboe, T. Trait biogeography of marine copepods: An analysis across scales. Ecol. Lett. 19, 1403–1413. https://doi.org/10.1111/ele.12688 (2016).Article 
PubMed 

Google Scholar 
17.Horne, C. R., Hirst, A. G., Atkinson, D., Neves, A. & Kiørboe, T. A global synthesis of seasonal temperature–size responses in copepods. Glob. Ecol. Biogeogr. 25, 988–999. https://doi.org/10.1111/geb.12460 (2016).Article 

Google Scholar 
18.Garzke, J., Hansen, T., Ismar, S. M. H. & Sommer, U. Combined effects of ocean warming and acidification on copepod abundance, body size and fatty acid content. PLoS ONE 11, e0155952. https://doi.org/10.1371/journal.pone.0155952 (2006).CAS 
Article 

Google Scholar 
19.Stelzer, C. P. Phenotypic plasticity of body size at different temperatures in a planktonic rotifer: Mechanisms and adaptive significance. Funct. Ecol. 16, 835–841. https://doi.org/10.1046/j.1365-2435.2002.00693.x (2002).Article 

Google Scholar 
20.Riemer, K., Anderson-Teixeira, K. J., Smith, F. A., Harris, D. J. & Ernest, S. K. M. Body size shifts influence effects of increasing temperatures on ectotherm metabolism. Global Ecol. Biogeogr. 27, 958–967. https://doi.org/10.1111/geb.12757 (2018).Article 

Google Scholar 
21.Gorsky, G. et al. Digital zooplankton image analysis using the ZooScan integrated system. J. Plankton Res. 32, 285–303. https://doi.org/10.1093/plankt/fbp124 (2010).Article 

Google Scholar 
22.Ibarbalz, F. M. et al. Global trends in marine plankton diversity across kingdoms of life. Cell 179, 1084–1097. https://doi.org/10.1016/j.cell.2019.10.008 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
23.Hoefnagel, K. N. & Verberk, W. C. Is the temperature-size rule mediated by oxygen in aquatic ectotherms?. J. Therm. Biol. 54, 56–65. https://doi.org/10.1016/j.jtherbio.2014.12.003 (2015).Article 
PubMed 

Google Scholar 
24.Wojewodzic, M. W., Kyle, M., Elser, J. J., Hessen, D. O. & Andersen, T. Joint effect of phosphorus limitation and temperature on alkaline phosphatase activity and somatic growth in Daphnia magna. Oecologia 165, 837–846. https://doi.org/10.1007/s00442-010-1863-2 (2011).ADS 
Article 
PubMed 

Google Scholar 
25.Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251. https://doi.org/10.1126/science.1061967 (2001).ADS 
CAS 
Article 
PubMed 

Google Scholar 
26.Czarnoleski, M., Ejsmont-Karabin, J., Angilletta, M. K. & Kozlowski, J. Colder rotifers grow larger but only in oxygenated waters. Ecosphere 6, 1–5. https://doi.org/10.1890/ES15-00024.1 (2015).Article 

Google Scholar 
27.Kiørboe, T. How zooplankton feed: Mechanisms, traits and trade-offs. Biol. Rev. 86, 311–339. https://doi.org/10.1111/j.1469-185X.2010.00148.x (2011).Article 
PubMed 

Google Scholar 
28.Benedetti, F., Gasparini, S. & Ayata, S.-D. Identifying copepod functional groups from species functional traits. J. Plankton Res. 38, 159–166. https://doi.org/10.1093/plankt/fbv096 (2016).Article 
PubMed 

Google Scholar 
29.Brun, P., Payne, M. R. & Kiørboe, T. A trait database for marine copepods. Earth Syst. Sci. Data 9, 99–113. https://doi.org/10.5194/essd-9-99-2017 (2017).ADS 
Article 

Google Scholar 
30.Anderson, T. R. Plankton functional type modelling: Running before we can walk?. J. Plankton Res. 27, 1073–1081. https://doi.org/10.1093/plankt/fbi076 (2005).ADS 
Article 

Google Scholar 
31.Biard, T. et al. In situ observations unveil an unexpectedly large biomass of Radiolaria and Phaeodaria (Rhizaria) in the oceans. Nature 532, 504–507. https://doi.org/10.1038/nature17652 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
32.Takagi, H. et al. Characterizing photosymbiosis in modern planktonic foraminifera. Biogeosciences 16, 3377–3396. https://doi.org/10.5194/bg-16-3377-2019 (2019).ADS 
CAS 
Article 

Google Scholar 
33.Rink, S., Kühl, M., Bijma, J. & Spero, H. J. Microsensor studies of photosynthesis and respiration in the symbiotic foraminifer Orbulina universa. Mar. Biol. 131, 583–595. https://doi.org/10.1007/s002270050350 (1998).Article 

Google Scholar 
34.Lombard, F., Erez, J., Michel, E. & Labeyrie, L. Temperature effect on respiration and photosynthesis of the symbiont-bearing planktonic foraminifera Globigerinoides ruber, Orbulina universa, and Globigerinella siphonifera. Limnol. Oceanogr. 54, 210–218. https://doi.org/10.4319/lo.2009.54.1.0210 (2009).ADS 
CAS 
Article 

Google Scholar 
35.Lesser, M. P. Coral Bleaching: Causes and Mechanisms. In Coral Reefs: An Ecosystem in Transition (eds Dubinsky, Z. & Stambler, N.) (Springer, 2011).
Google Scholar 
36.Villar, E. et al. Symbiont chloroplasts remain active during bleaching-like response induced by thermal stress in Collozoum pelagicum (Collodaria, Retaria). Front. Mar. Sci. 5, 387. https://doi.org/10.3389/fmars.2018.00387 (2018).Article 

Google Scholar 
37.Hemleben, C., Spindler, M. & Anderson, O. R. Modern Planktonic Foraminifera (Springer-Verlag, 1989).Book 

Google Scholar 
38.Suzuki, N. & Not, F. Biology and ecology of radiolaria. In Marine Protists: Diversity and Dynamics (eds Ohtsuka, S. et al.) (Springer, 2015).
Google Scholar 
39.de Puelles, F. et al. Zooplankton abundance and diversity in the tropical and subtropical ocean. Diversity 11, 203. https://doi.org/10.3390/d11110203 (2019).CAS 
Article 

Google Scholar 
40.Beaugrand, G., Edwards, M. & Legendre, L. Marine biodiversity, ecosystem functioning, and carbon cycles. PNAS 107, 10120–10124. https://doi.org/10.1073/pnas.0913855107 (2010).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
41.Brun, P. et al. Climate change has altered zooplankton-fuelled carbon export in the North Atlantic. Nat. Ecol. Evol. 3, 416–423. https://doi.org/10.1038/s41559-018-0780-3 (2019).Article 
PubMed 

Google Scholar 
42.Buitenhuis, E. T., Le Quéré, C., Bednaršek, N. & Schiebel, R. Large contribution of pteropods to shallow CaCO3 export. Glob. Biogeochem. Cyc. 33, 458–468. https://doi.org/10.1029/2018GB006110 (2019).ADS 
CAS 
Article 

Google Scholar 
43.Follows, M. J., Dutkiewicz, J., Grant, S. & Chisholm, S. W. Emergent biogeography of microbial communities in a model ocean. Science 315, 1843–1846. https://doi.org/10.1126/science.1138544 (2007).ADS 
CAS 
Article 
PubMed 

Google Scholar 
44.Ward, B. A., Dutkiewicz, S., Jahn, O. & Follows, J. F. A size-structured food-web model for the global ocean. Limnol. Oceanogr. 57, 1877–1891. https://doi.org/10.4319/lo.2012.57.6.1877 (2012).ADS 
Article 

Google Scholar 
45.Sailley, S. F. et al. Comparing food web structures and dynamics across a suite of global marine ecosystem models. Ecol. Model. 261, 43–57. https://doi.org/10.1016/j.ecolmodel.2013.04.006 (2013).Article 

Google Scholar 
46.Le Quéré, C. et al. Role of zooplankton dynamics for Southern Ocean phytoplankton biomass and global biogeochemical cycles. Biogeosciences 13, 4111–4133. https://doi.org/10.5194/bg-13-4111-2016 (2016).ADS 
CAS 
Article 

Google Scholar 
47.Kwiatkowski, L. et al. Emergent constraints on projections of declining primary production in the tropical oceans. Nat. Clim. Change 7, 355–358. https://doi.org/10.1038/nclimate3265 (2017).ADS 
CAS 
Article 

Google Scholar 
48.Sunagawa, S. et al. Tara Oceans: Towards global ocean ecosystems biology. Nat. Rev. Microbiol. 18, 428–445. https://doi.org/10.1038/s41579-020-0364-5 (2020).CAS 
Article 
PubMed 

Google Scholar 
49.Pesant, S. et al. Open science resources for the discovery and analysis of Tara Oceans data. Sci. Data 2, 150023. https://doi.org/10.1038/sdata.2015.23 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
50.Picheral, M. et al. Vertical profiles of environmental parameters measured on discrete water samples collected with Niskin bottles at station TARA_147 during the Tara Oceans expedition 2009–2013. PANGAEA https://doi.org/10.1594/PANGAEA.839235 (2014).51.Guidi, L. et al. Environmental context of all samples from the Tara Oceans Expedition (2009–2013), about sensor data in the targeted environmental feature. PANGAEA https://doi.org/10.1594/PANGAEA.875576 (2017).52.Guidi, L. et al. Environmental context of all samples from the Tara Oceans Expedition (2009–2013), about pigment concentrations (HPLC) in the targeted environmental feature. PANGAEA https://doi.org/10.1594/PANGAEA.875569 (2017).53.Guidi, L. et al. Environmental context of all samples from the Tara Oceans Expedition (2009–2013), about nutrients in the targeted environmental feature. PANGAEA https://doi.org/10.1594/PANGAEA.875575 (2017).54.Speich, S. et al. Environmental context of all samples from the Tara Oceans Expedition (2009–2013), about the water column features at the sampling location. PANGAEA https://doi.org/10.1594/PANGAEA.875579 (2017).55.de Boyer-Montegut, C., Madec, G., Fischer, A. S., Lazar, A. & Iudicone, D. Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. J. Geophys. Res. 109, C12003. https://doi.org/10.1029/2004JC002378 (2004).ADS 
Article 

Google Scholar 
56.Aminot, A., Kérouel, R. & Coverly, S. C. Nutrients in seawater using segmented flow analysis. In Practical Guidelines for the Analysis of Seawater (ed. Wurl, O.) (CRC Press, 2009).
Google Scholar 
57.Uitz, J., Claustre, H., Morel, A. & Hooker, S. B. Vertical distribution of phytoplankton communities in Open Ocean: An assessment based on surface chlorophyll. J. Geophys. Res. 111, C08005. https://doi.org/10.1029/2005JC003207 (2006).ADS 
Article 

Google Scholar 
58.Pante, E. & Simon-Bouhet, B. marmap: A Package for importing, plotting and analyzing bathymetric and topographic data in R. PLoS ONE 8(9), e73051. https://doi.org/10.1371/journal.pone.0073051 (2013).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
59.Picheral, M., Colin, S. & Irisson J.-O. EcoTaxa, A Tool for the Taxonomic Classification of Images. http://ecotaxa.obs-vlfr.fr (2017).60.Wood, S. N. Generalized Additive Models: An Introduction with R 2nd edn. (Chapman and Hall/CRC, 2017).Book 

Google Scholar 
61.Dormann, C. F. et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46. https://doi.org/10.1111/j.1600-0587.2012.07348.x (2013).Article 

Google Scholar 
62.Giorgino, T. Computing and visualizing dynamic time warping alignments in R: The dtw package. J. Stat. Softw. 31, 1–24. https://doi.org/10.18637/jss.v031.i07 (2009).Article 

Google Scholar 
63.Park, H.-S. & Jun, C.-H. A simple and fast algorithm for K-medoids clustering. Expert Syst. Appl. 36, 3336–3341. https://doi.org/10.1016/j.eswa.2008.01.039 (2009).Article 

Google Scholar 
64.R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2018). https://www.R-project.org/.65.Wickham, H. et al. Welcome to the Tidyverse. J. Open Source Softw. 4, 1686. https://doi.org/10.21105/joss.01686 (2019).ADS 
Article 

Google Scholar 
66.Heiberger, R. M. HH: Statistical Analysis and Data Display: Heiberger and Holland. R package version 3.1–40, https://CRAN.R-project.org/package=HH (2020).67.Lê, S., Josse, J. & Husson, F. FactoMineR: An R package for multivariate analysis. J. Stat. Softw. 25, 1–18. https://doi.org/10.18637/jss.v025.i01 (2008).Article 

Google Scholar 
68.Sarda-Espinosa, A. dtwclust: Time Series Clustering Along with Optimizations for the Dynamic Time Warping Distance. R package version 5.5.6, https://CRAN.R-project.org/package=dtwclust (2019). More