Ecology

NEWS AND VIEWS
04 August 2021

A cocktail of pesticides, parasites and hunger leaves bees down and out

Pollinators are under threat. A meta-analysis reveals that the combination of agrochemicals, parasites and malnutrition has a cumulative negative effect on bees, and that pesticide–pesticide interactions increase bee mortality.

Adam J. Vanbergen

 ORCID: http://orcid.org/0000-0001-8320-5535

0

Adam J. Vanbergen

Adam J. Vanbergen is in the Department of Plant Health and Environment, INRAE (the National Research Institute for Agriculture, Food and Environment), Dijon 21000, France.

View author publications

You can also search for this author in PubMed
 Google Scholar

Share on Twitter
Share on Twitter

Share on Facebook
Share on Facebook

Share via E-Mail
Share via E-Mail

Download PDF

Worldwide, an estimated 20,000 species of wild and managed insects pollinate flowers, aiding plant reproduction1. In doing so, they form a key link in the tangled web of species interactions that support biodiverse and healthy ecosystems1,2. Moreover, humans enjoy a variety of sociocultural and economic benefits from pollinator biodiversity2,3, and pollination secures crop yields that supply essential nutrients and healthy, diverse diets1,4. Writing in Nature, Siviter et al.5 report a pollinator threat that jeopardizes these benefits.Pollinators and pollination are threatened by environmental pressures, including many that are a consequence of human activity (Fig. 1). These pressures include land-use and climate change2,6, intensive agriculture7, the spread of invasive alien species and problems with pests and disease-causing agents (pathogens)2,8. The individual effects of these pressures on pollinators are well established1,2, raising the question of whether an interplay between these various pressures exacerbates the overall risk that they pose to pollinators and pollination9–11. This issue has been recognized by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, which stated2 in 2016 that “many drivers that directly impact the health, diversity and abundance of pollinators … can combine in their effects and thereby increase the overall pressure on pollinators”.

Figure 1 | The effect of multiple stressors on bees. Agricultural intensification has put pollinators under pressure. a, Practices associated with intensive farming reduce food availability for pollinators1,2,7,12. b, Managed high-density bee colonies for crop pollination are associated with these pollinators being at risk of disease and parasite infection2. This poses a risk of illness spreading to wild bees. c, Exposure to a variety of pesticides poses another risk to pollinators. d, Siviter et al.5 present a meta-analysis that indicates the consequences for bees of combinations (some examples shown) of these challenges.

Intensive agriculture is a multifactorial source of stress on pollinator populations1,7,10,11. Pollinating insects, such as bees, face the physiological challenge of acute or chronic harm from exposure to various agrochemicals, including fungicides and pesticides, that are used to protect crop plants. They also face nutritional stress arising from the lack of pollen- and nectar-providing wild flowers in large-scale, intensive crop monocultures1,2,7,12. Moreover, the industrial transport and use of managed high-density colonies of honey bees (Apis mellifera) for crop pollination can increase pollinator exposure to parasites or pathogens2, and might result in disease spillover to wild pollinators13. Over the past decade, the lethal or sublethal effect of combinations of agrochemical, pathogenic or nutritional stressors on bees has been tested in many individual experiments2,9,10.
Read the paper: Agrochemicals interact synergistically to increase bee mortality
Siviter et al. advance this knowledge through a quantitative meta-analysis of the effect of interactions between agrochemical, pathogenic and nutritional stressors on multiple aspects of bee health and fitness. Their analysis is notable because of the breadth of bee responses considered (for example, foraging behaviour, memory, mortality and colony reproduction), and for comparisons of the interactions of multiple classes of stressor (for example, agrochemical–parasite, parasite–nutrition, agrochemical–agrochemical and parasite–parasite interactions).The authors conducted a monumental literature search that yielded almost 15,000 relevant individual studies. Siviter and colleagues combed through these publications to focus on the experiments that investigated the combined effect of parasites (microorganisms and invertebrates), agrochemicals and nutritional stressors on bee health. The authors selected studies that used a balanced and replicated experimental design, and that provided accessible data (means, standard deviations and sample sizes) for each treatment. This rigorous focus and quality control resulted in a final set of 90 studies being selected for further analysis.These studies provide a total of 356 effect sizes (measurements indicating the magnitude of a relationship between factors of interest and a particular outcome) for different stressor and bee-response combinations. The authors accounted for data issues that might have confounded their accurate detection of bee responses. Such challenges included those arising from statistical non-independence of multiple effects reported from a single study, publication biases (for example, the lack of negative results), species skews (honey bee data sets predominated), and how experimental treatments such as pesticide dose compare with what might be realistically encountered in the field (termed field realism).
Robust evidence of declines in insect abundance and biodiversity
Siviter and colleagues tested whether the stressor interactions were synergistic, meaning that their combined effect was greater than the sum of their individual effects, as would be the case if the effect of one stressor on a bee elevates the effect of another stressor. The authors also examined alternative scenarios in which the effects of multiple stressors were antagonistic (the effect of one stressor lessens the effect of another) or additive (the combined effect is equivalent to the sum of the individual effects).A consistent message from their analysis is that bee mortality is increased by a synergistic interaction between multiple stressors — the worst-case scenario, indicating a disproportionate effect of multiple stressors on bee survival. Interactions between different agrochemicals, rather than other stressors, drove this overall effect, and this finding held true when accounting for the field realism of the agrochemical doses. This result confirms that the cocktail of agrochemicals that bees encounter in an intensively farmed environment can create a risk to bee populations1,2,9,14. Multi-stressor interactions involving parasites and nutritional stress (including in combination with agrochemicals) produced additive effects on bee mortality overall.The authors’ deeper analysis of the biological complexity, however, revealed large differences between particular parasite groups in terms of the full range of additive, antagonistic and synergistic effects on bee mortality, when considering interactions between different parasites or between different parasites and nutritional stress. This variability in response, together with the lower sample sizes for the interactions involving stressors other than agrochemicals, indicate a caveat to consider and also suggest a need for more research on the combined effects of biological sources of stress.It is intriguing that Siviter and colleagues found that additive, not synergistic, effects predominated for the non-lethal effects of stressors on fitness proxies (such as modifications of bee behaviour or reproduction, changes in parasite load or immune function). Such non-lethal changes could ultimately affect bee mortality rates. Consequently, how the observed synergistic effects of agrochemicals on bee mortality arise remains to be established. More work is therefore needed to identify the mechanism that links exposure to behavioural or physiological changes and mortality.
An alternative to controversial pesticides still harms bumblebees
The majority of the studies in the data set were of managed populations of A. mellifera, so the authors also separately analysed responses at the level of bee genus (Apis, Bombus, Megachile and Osmia). Apis mortality was affected by a synergistic multi-stressor interaction qualitatively similar to the full analysis of all bee genera. Other bee genera exhibited additive or antagonistic mortality responses from many fewer studies. This raises an important point. There is a need for research efforts and regulators to widen their focus from A. mellifera — a single, mostly managed bee species — to other pollinator model organisms, whose different ecology and evolutionary history might result in different responses to stressors10.Siviter and colleagues’ findings of the cumulative negative effect of multi-stressor interactions on bees reinforces the call to evaluate such interactions to avoid unforeseen risks to biodiversity and healthy ecosystems1,9,10. In some regions of the world, regulatory risk-assessment frameworks for plant-protection products are being developed to deal with sublethal, long-term and potentially synergistic effects among stressors15 (see go.nature.com/3f4ax5r), but their biological and geographical scope must be extended. The authors acknowledge that the high levels of variability between the studies and parameters investigated demand an appropriately cautious interpretation. However, this highlights the need for worldwide reconsideration of risk-assessment approaches for pesticide regulation.Given the widespread loss of habitat resources — such as pollen and nectar sources — from intensively managed agricultural landscapes7,12, nutritional deficits occurred surprisingly infrequently as a mechanism underlying bees’ physiological stress (they accounted for only 58 out of the 365 measurements of effect sizes). A greater consideration of how nutritional stress interacts with exposure to pathogens and agrochemicals is therefore an obvious research gap to fill. Moreover, ensuring that experimental treatments are calibrated to simulate realistic environmental conditions would greatly aid risk assessments. This might include three-way combinations of field-realistic chemical doses and parasite levels, and a spatio-temporal dietary diversity similar to that found in semi-natural or highly human-modified landscapes.The next challenge is to look beyond these parasite–nutrition–agrochemical interactions to consider other risks to pollination. Future studies must ultimately consider, through a combination of correlative and experimental approaches, the interplay of nutrition–pathogen–agrochemical interactions alongside the effects of other human-driven changes (such as climate change, pollution, land-use changes and the spread of invasive species)1,2,11. Although such assessments would be non-trivial to carry out, they will be vital for understanding and ranking the relative risks to pollinators and pollination that are coming from multiple combinations of pressures resulting from global changes.

doi: https://doi.org/10.1038/d41586-021-02079-4

References1.Potts, S. G. et al. Nature 540, 220–229 (2016).PubMed 
Article 

Google Scholar 
2.IPBES. The Assessment Report of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services on Pollinators, Pollination and Food Production (IPBES Secretariat, 2016).
Google Scholar 
3.Hill, R. et al. Nature Sustain. 2, 214–222 (2019).Article 

Google Scholar 
4.Chaplin-Kramer, R. et al. Science 366, 255–258 (2019).PubMed 
Article 

Google Scholar 
5.Siviter, H. et al. Nature https://doi.org/10.1038/s41586-021-03787-7 (2021).Article 

Google Scholar 
6.Soroye, P., Newbold, T. & Kerr, J. Science 367, 685–688 (2020).PubMed 
Article 

Google Scholar 
7.Kovács-Hostyánszki, A. et al. Ecol. Lett. 20, 673–689 (2017).PubMed 
Article 

Google Scholar 
8.Vanbergen, A. J., Espíndola, A. & Aizen, M. A. Nature Ecol. Evol. 2, 16–25 (2018).PubMed 
Article 

Google Scholar 
9.Goulson, D., Nicholls, E., Botías, C. & Rotheray, E. L. Science 347, 1255957 (2015).PubMed 
Article 

Google Scholar 
10.Vanbergen, A. J. & the Insect Pollinators Initiative. Front. Ecol. Environ. 11, 251–259 (2013).Article 

Google Scholar 
11.González-Varo, J. P. et al. Trends Ecol. Evol. 28, 524–530 (2013).PubMed 
Article 

Google Scholar 
12.Baude, M. et al. Nature 530, 85–88 (2016).PubMed 
Article 

Google Scholar 
13.Manley, R. et al. Ecol. Lett. 22, 1306–1315 (2019).PubMed 

Google Scholar 
14.Woodcock, B. A. et al. Nature Commun. 7, 12459 (2016).PubMed 
Article 

Google Scholar 
15.European Food Safety Authority. EFSA J. 11, 3295 (2013).Article 

Google Scholar 
Download references

Competing Interests
The author declares no competing interests.

Related Articles

Read the paper: Agrochemicals interact synergistically to increase bee mortality

An alternative to controversial pesticides still harms bumblebees

Robust evidence of declines in insect abundance and biodiversity

See all News & Views

Subjects

Ecology

Agriculture

Latest on:

Ecology

The world’s species are playing musical chairs: how will it end?
News Feature 04 AUG 21

Agrochemicals interact synergistically to increase bee mortality
Article 04 AUG 21

Biodiversity needs every tool in the box: use OECMs
Comment 26 JUL 21

Agriculture

Europe’s reformed agricultural policy disappoints
Correspondence 27 JUL 21

Italy: scientists petition against biodynamic farming law
Correspondence 13 JUL 21

Does the fight against hunger need its own IPCC?
Editorial 13 JUL 21

Jobs

Associate Editor or Senior Editor, BMC Series

Springer Nature
New York, NY, United States

Faculty Positions in Computational Oncology, Memorial Sloan Kettering Cancer Center

Memorial Sloan Kettering Cancer Center (MSKCC)
New York, NY, United States

Journals Operations Coordinator – BioMed Central

Springer Nature
London, United Kingdom

Molecular Scientist – Team Leader

Meridian Bioscience
London, United Kingdom

Nature Briefing
An essential round-up of science news, opinion and analysis, delivered to your inbox every weekday.

Email address

Yes! Sign me up to receive the daily Nature Briefing email. I agree my information will be processed in accordance with the Nature and Springer Nature Limited Privacy Policy.

Sign up More

Temperature profiles, sampling methods, and physical descriptions were compiled for 153 globally-distributed lakes. We selected vertical temperature profiles measured in the same location for multiple years, typically the deepest part of the lake, using a single mid-summer profile for each lake and year. Data are presented for all years for which profiles were available, with both raw vertical temperature profiles and interpolated vertical temperature profiles presented. Depth increments for raw profiles vary by lake based on sampling methodology, as described for each lake below. For interpolation, the selected temperature profiles were linearly interpolated or binned to 0.5 m increments throughout the water column. Mid-summer profiles for each lake were selected as described in (7). Briefly, we calculated relative thermal resistance to mixing28,29 (RTR) for all profiles for each lake for each year. The day of year with the maximum RTR value for each year was selected for each lake, and the median day of year across all years per lake was considered the target day for single mid-summer profile selection. For each lake and year, the temperature profile nearest to this median day of year ± 21 days was selected. If no profile spanning surface to near the maximum lake depth was available, then no profile was selected. Details on sampling methods from each contributing research group follows, organized by major geographic region and alphabetized within geographic region. Additional details for water chemistry sampling methods can be found in the Supplementary Information (Supplementary Table S1).Western North AmericaCastle Lake in California, USA has been sampled annually since 1959. Temperature data were measured between 12:00 and 14:00 at the deepest part of the lake in 1 m depth intervals. Measurements were taken using a reversing thermometer (1959–1975), Hydrolab (1975–1982; 2011–2013), and Yellow Springs Instruments (YSI) model 85 (1999–2011) (sampling instrument is unknown for most of the 1980s and 90s).Crater Lake is located in the northwestern USA at the crest of the Cascade Mountains in Oregon. Temperature data were collected between 8:00 and 17:00 from the middle of the lake in the deepest basin. From 1988 to present, temperature profiles were measured by lowering a Seabird Instruments SBE19 CTD through the water column at a rate of approximately 0.5 m per second, collecting two readings per second. Data were binned to 1 m increments. Prior to 1988, temperature profiles were measured manually with a Montedoro-Whitney thermistor with a 250 m cable.Emerald Lake, California, USA has temperature profiles measured from the deepest part of the lake beginning in 1984. Temperature was measured in 1 m increments throughout the water column from 1984 to 2006 using a YSI 58, after which a thermistor chain was deployed with Onset Water Temp Pro spaced at 0.5 to 1 m intervals.Flathead Lake in northwest Montana, USA is the largest freshwater lake west of the Mississippi River by surface area. Temperature data were collected fortnightly in June, July, and August 1978–2013 with various instruments. From 1978–1984, data were collected with a Hydrolab I; from 1985–1993, with a Hydrolab Surveyor II; from 1993–2003, with a Hydrolab Surveyor III; and from 2003–2001, with a Hydrolab Surveyor IV. Since 2011, data have been collected with a Hydrolab DS5 unit. Measurements were taken between 10:00 and 14:00 at the deepest point of the lake in 1 to 10 m depth intervals throughout the water column, with greater intervals between measurements at deeper depths. The intervals between depths on a given sampling date also changed periodically early in the data set before consistent intervals were established. Instrumentation was calibrated before each sampling event to account for the elevation of the sampling site and highly variable barometric pressure.Washington Lake is located in Washington, USA. Temperature profiles were measured between 9:00 and 16:00 in the main trench of the lake. Temperature was measured using various bathythermographs between 1933–1986, a Kahl digital temperature meter 202WA510 beginning in 1974, and a YSI 6600 V2 sonde beginning in 2012. For measurements with the bathythermographs, temperature was measured every 1 m close to the thermocline and every 5 m elsewhere. When the Kahl temperature meter was used, temperature was recorded for every meter through 20 m and in 5 m increments after 20 m depth. Data from the YSI sonde were recorded continuously, so temperature readings were used at similar intervals to those from the previous instruments.Central North AmericaActon Lake is a eutrophic reservoir in southeastern Ohio, USA. Temperature was measured in the deepest part of the lake at 0.5 m depth increments from 0 to 5 m and 1 m increments for the remainder of the water column. Temperature data were typically collected in the mid-morning and occasionally in the early afternoon in mid-July of each year from 1992–2012 using a handheld YSI temperature/dissolved oxygen sensor or a YSI sonde.Lakes Bighorn, Harrison, Pipit, and Snowflake are located along the eastern front range of the Rocky Mountains in the Cascade Valley of Banff National Park (Alberta, Canada). Temperature data were collected in the deepest point of the lake in 1 m depth intervals between 10:00 and 15:00 using a MK II Thermistor (Flett Research Ltd, Winnipeg, Canada) for measurements taken beginning in the 1990s. Earlier measurements (1960s and 1970s) were taken using a YSI Model 425 C thermistor thermometer, which was calibrated against a mercury thermometer.Douglas Lake is a mildly eutrophic glacially formed multiple ice-block kettle lake in northwestern Cheboygan County, Michigan, USA. Temperature data were collected from 1913–2014 at 1 m intervals from 0 to 24 m in the center of South Fishtail Bay Kettle with a reversing thermometer at approximately 1 m increments (1913–1970), a YSI model 54 A oxygen electrode-thermistor thermometer (1971–1982), a Hydrolab MS-5 multiprobe (1983–2009), and an 8-node MHL thermistor string (2010–2014).Lakes Eucha, Grand Lake O’ the Cherokees, Spavinaw, Texoma, and Thunderbird are reservoirs located in Oklahoma, USA. Temperature profiles were measured at the deepest point in these reservoirs between morning and mid-day with various sensors. YSI temperature probes were used for all samples collected in Lake Thunderbird, through 1991 in Grand Lake O’ the Cherokees, through 1995 in Lakes Eucha and Spavinaw, and through 2000 in Lake Texoma. Beginning in 2000 in Lake Texoma, Hydrolabs were used for temperature measurements, with a Hydrolab H2O through 2008, and Hydrolab DSX5 thereafter. A Hydrolab was used for measurements in Grand Lake O’ the Cherokees from 2011–2013. Lakes Eucha and Spavinaw were sampled with a Hydrolab H2O through 2005, a YSI 6930 V2 for samples from 2006–2012, and a YSI EXO1 for samples after 2012.Katepwa Lake is a eutrophic, riverine site located in the Qu’Appelle River drainage basin in southern Saskatchewan, Canada. The lake has been sampled since 1994 as part of the Qu’Appelle Long-term Ecological Research network (QU-LTER). Temperature data were collected from a standard deep site in 1 m intervals between 10:30 and 13:00 using a YSI 85 or similar multi-parameter probe.Northeastern North AmericaTemperature profiles were collected beginning in 1969 for six lakes (Lakes 222, 224, 239, 240, 373, and 442) at the IISD Experimental Lakes Area (International Institute for Sustainable Development, Northwestern Ontario, Canada). Profiles were measured in the deepest part of the lake in 1 m intervals, except in the thermocline where temperature was measured every 0.25 m. Montedoro-Whitney thermistors (models TC-5A and TC-5C) were used through 1983, a Flett Research Mark II digital telethermometer was used for temperature measurements from 1984–2009, and a RBR XRX620 multifunction probe with integrated temperature sensor for measurements beginning in 2010.The Dorset “A lakes”, Blue Chalk, Chub, Crosson, Dickie, Harp, Heney, Plastic, and Red Chalk Lakes are located in the Muskoka-Haliburton region of south-central Ontario, Canada. The study sites are primarily small headwater lakes, with the exception of Red Chalk Lake which is located downstream of Blue Chalk Lake. Temperature data were collected from the deepest point in each lake using a YSI 58 temperature/dissolved oxygen meter (or occasionally a digital YSI 95 meter) beginning in the late 1970s. Measurements were collected between 9:00 and 16:00, with readings taken every 1 m from the lake surface (0.1 m depth) to within approximately 1 m of lake sediments.Bubble Pond, Eagle Lake, and Jordan Pond are located on Mount Desert Island off the coast of Maine, USA. Temperature data were collected from the location of the maximum depth at 1 m increments using a YSI 600XL multiparameter water-quality monitor (sonde) from 2006 to present and a YSI 54ARC before this time.Lake Champlain is located in the northeastern USA, on the border of Vermont and New York state and partly extending into Quebéc, Canada. Temperature data were collected from a sampling station in Mallet’s Bay beginning in 1992. Measurements of temperature were taken at 1 m intervals between 8:00 and 17:00 using a Hydrolab MS-5 multi-probe sonde.Clearwater, Sans Chambre, and Whitepine Lakes are located in northeastern Ontario, Canada, and Hawley Lake is in the Hudson Bay Lowlands area of subarctic Ontario, Canada. Temperature profiles were measured from near the area of maximum depth, beginning in the 1970s in Clearwater and Hawley, and beginning in the 1980s in Sans Chambre and Whitepine. Measurements were taken in 1 m intervals through the water column between 12:00 and 17:00. Clearwater, Sans Chambre, and Whitepine Lakes’ temperature profiles were measured using various YSI temperature/dissolved oxygen meters (models 50B, 51B, 52, 54, or 58) beginning in 1998, with a YSI model 54 temperature/dissolved oxygen meter used for Sans Chambre and Whitepine for most earlier years’ measurements. In Clearwater Lake, a YSI model 432D telethermometer was used through 1975, a Montedoro-Whitney TC-5C thermistor from 1976–1981, and a Mark II Telethermometer from 1982–1998. In Hawley Lake, various YSI temperature probes were used in earlier years, and since 2009, a YSI Pro ODO meter was used to measure water temperature.Lakes Giles and Lacawac are located in the Pocono Mountains region of Pennsylvania, USA. Temperature data were collected from the deepest point in the lake in late July or early August each year beginning in 1988, between 9:00 and 16:00. Temperature measurements were measured in 1 m increments using a YSI 58 temperature/dissolved oxygen meter (1988–1992), or with a rapid recording Biospherical Instruments PUV 500 (1993–2003) or BIC 2104 P (2004-present) recording at 4 Hz while being lowered through the water column. Profiles taken with the YSI were linearly interpolated to 0.5 m depth increments, and profiles taken with the PUV and BIC were binned to 0.5 m depth increments.Lake Lillinonah is a reservoir on the Housatonic River located in western Connecticut, USA. Temperature data were collected at the deepest point of the lake between 9:00 and 17:00 beginning in 1996 by First Light Power. Measurements were collected every 5 m using a YSI 58 temperature/dissolved oxygen meter.Mohonk Lake is a small glacial lake with a single deep basin located at the Shawangunk Ridge, New York, USA. Temperature data were collected from the northern end of the lake from 1984–2013 during daytime. Temperature measurements were measured manually in 1 m increments using Digi-sense Economical Thermistor 400 series (Model #93210-00). Additional weekly temperature profiles are publicly available30.Temperature profiles were compiled from 11 lakes from the North Temperate Long Term Ecological Research Network in Wisconsin, USA. Fish, Mendota, Monona, and Wingra Lakes are located in southern Wisconsin, and Allequash, Big Muskellunge, Crystal Bog, Crystal Lake, Sparkling, Trout Bog, and Trout Lake are in northern Wisconsin. Temperature profiles were measured in 1 m increments from the surface to lake bottom at the deepest location in each lake since 1981 (northern lakes) and 1995–1996 (southern lakes), and since 1894 in Mendota31. Various temperature/dissolved oxygen probes were used to collect these data, and were calibrated in the field prior to data collection. For the northern lakes, a Montedoro Whitney CTU-3B sensor was used for some data collected between 1981–1986, a Whitney TC-5C for 1982–1983, Whitney DOR-2A in 1984, YSI-57 in 1983 and 1985, and a YSI-58 for 1985 onward. Temperature data in the southern lakes were collected with a YSI-58 temperature/dissolved oxygen sensor.Lake Opeongo is the largest oligotrophic deepwater lake in Algonquin Provincial Park, Ontario, Canada. Temperature data were collected by Ontario Ministry of Natural Resources and Forestry staff at Harkness Fisheries Research Station in the west basin of Opeongo’s South Arm in midsummer in 1958–1965, and between 9:00 and 16:00 in 1998–2014. Temperature loggers (Hobo Tidbits, 1998–2004; Onset Water Temperature Pro V1, 2004–2008; and Onset Pro V2, 2009–2014) were installed after ice out in the same location at approximately 1.5 m intervals to a 15 m depth, with a deepwater thermistor placed at approximately 20 m. Temperatures were recorded at high temporal resolution (10 to 15 minute intervals) throughout the summer using these temperature strings. Earlier (1958–1965) temperature profiles were measured using handheld thermistors.Lake Sunapee is the fifth largest lake located within New Hampshire, USA. Temperature data were collected in the morning from the deepest point of the lake in the central basin beginning in 1986. Temperature measurements were measured manually in 1 m increments using a YSI 52 temperature meter.Lake Wallenpaupack is a reservoir in the Pocono Mountains region of Pennsylvania, USA. Temperature profiles were measured in the center of the lake during daytime hours, with measurements taken every 0.5 to 1 m. Temperature was measured with various YSI instruments including a YSI 610-DM/600XL (2002–2005), YSI 85 (2008–2010), YSI 600XL sonde with 600D datalogger (2011–2012), and various temperature sensors prior to 2002.Southeastern North AmericaLake Annie is located in the Lake Wales Ridge region, Florida, USA. Temperature data were collected from the deepest point in the lake on a monthly basis between 9:00 and 16:00 beginning in 1984. Temperature measurements were measured manually in 1 m increments using a Montedoro Corporation Thermistor Model TC-5c (1984–2008) and a YSI Pro Plus (2009–2014) with values measured while lowering the meter to the bottom of the lake and again when raising it to the surface. The data reported are the means of the two depth-specific values.Temperature data were measured in Lake Okeechobee in Florida, USA beginning in 1973. Temperature data were measured near the surface at approximately 0.2 m between 8:00 and 12:00 using a Hydrolab through 1995 and a YSI 58 temperature sensor from 1996–2014.SubarcticAleknagik, Beverley, Chignik, Hidden, Kulik, Little Togiak, Lynx, and Nerka Lakes are located in Alaska, USA. Temperature profiles have been measured on these lakes since the 1960s (Aleknagik, Beverley, Kulik, Little Togiak, and Nerka Lakes), and since the early 2000s (Chignik, Hidden, and Lynx Lakes). Measurements were taken during the day, typically between 10:00 and 19:00. A bathythermograph was used for temperature measurements through 1967, a digital thermistor from 1968–1998, a YSI 660 sonde for samples from 1999–2012, and a YSI Castaway beginning in 2013, with the exception of samples from Chignik Lake, for which a handheld thermometer was used to measure the water temperature from Van Dorn casts.Toolik Lake is a kettle lake located in Alaska, USA. Temperature data were collected between June and August in the south basin of the lake between 9:00 and 11:00. Temperature was recorded in 1 m increments throughout the water column with a Hydrolab profiler sampling Surveyor 4a datalogger and a datasonde 4a multiprobe32,33,34,35.Vulture Lake is located in the sub-Arctic region of the Northwest Territories, Canada. Temperature data were collected from one of the deeper parts of the lake in late July or early August during the open-water season from 1997–2014. Temperature measurements were recorded using a multi-probe sonde (e.g., YSI 6820) at 0.5 m or 1 m increments (every 1 m for most years, except 1999 [every 0.5 m] and 2009 [every 0.2 m]). The corresponding depth was measured using the readings from the depth sensor. Measurements were collected continuously from just below the lake’s surface to approximately 0.5 m above the water-sediment interface. Sensors equipped on profiling instruments were calibrated following the manufacturer’s recommended frequency and methods to ensure accurate and reliable operation of the sensors in the field.Central and South AmericaLake Atitlán is a deep tropical mountain lake in the Guatemalan highlands, sampled in 1968–1969 and 2010–2011. Temperature profiles at the lake center were measured manually between 7:00 and 13:00 in variable increments (1 to 5 m) in the first 30 m using a YSI 51 or YSI 95 temperature/dissolved oxygen meter. For depths below 30 m, samples were collected using a Van Dorn bottle and temperature was measured immediately upon the sample reaching the surface.Lake Mascardi is located in the North Patagonian Andes in Argentina. Temperature data were collected near the deepest point in the Catedral arm of the lake in mid-summer (January-February) between 12:00 and 14:00 beginning in 1994. Temperature measurements were taken using rapid recording Biospherical Instruments PUV 500 (1994) or PUV 500B (1996–2014) recording at 4 Hz while being lowered through the water column.AfricaLake Kivu is located on the border between Rwanda and the Democratic Republic of Congo, and is one of the seven African Great Lakes. Temperature data were collected in the Ishungu basin beginning in 2002 between 9:00 and 16:00. Temperature was measured in 5 m increments using a YSI 55 temperature/dissolved oxygen meter (2002–2005), or with a suite of instruments (YSI 6600 V2, Hydrolab DS4a, DataSonde 4a 42071, Sea and Sun 725 and 257) recording at high frequency while being lowered through the water column. All temperature profiles were vertically interpolated to a regular vertical grid with 1 m increments using piecewise cubic Hermite interpolation18,26.Lake Nkugute is located in Uganda. Temperature data for Lake Nkugute were collected in the deepest part of the lake at a depth interval of 1 to 5 m using a liquid-in-glass thermometer set in a Ruttner sampler in 1964 with the contemporary (2002-present) data being measured with a YSI sonde.Lake Nkuruba is located in western Uganda, in the vicinity of Kibale National Park (northern sector). Temperature data were collected beginning in 1992 from the deepest part of the lake in 1 m increments through a depth of approximately 30 m. Temperature was measured manually using a YSI 50 or 51B temperature/dissolved oxygen meter.Lake Tanganyika has temperature profile data dating back to 1912. Temperature profiles in this lake were typically measured in the morning, between 9:00 and 12:00 from the north basin of the lake near Kigoma, Tanzania. Temperature profiles over this century of data collection have been taken using a multitude of instruments. Since 1993, temperature profiles have been measured using a YSI 6600 V2 sonde, titanium RBRduo TD, Seacat Profiler V3.1b, Onset HOBO U22 temperature loggers, CTD Seabird 19, STD-12 Plus CTP profiler, and a YSI 58 temperature/dissolved oxygen meter. Prior to 1993 various methods were used, including measurement of water temperature from Van Dorn casts with a mercury thermometer, various data loggers, reversing thermometers, and Niskin bottles and a bathythermograph36.Lake Victoria is shared between Tanzania, Uganda, and Kenya. Temperature data were collected from stations across the lake during acoustic surveys each year in 2000–2001 and in 2008 using a submersible Conductivity Temperature-Depth profiling system (CTD, Sea-bird Electronics, Sea Cat SBE 19).Scandinavia and Northern EuropeLakes Allgjuttern, Brunnsjön, Fiolen, Fracksjön, Övre Skärjön, Remmarsjön, Rotehogstjärnen, St. Skärsjön, Stensjön, and Stora Envättern are relatively small, boreal lakes in Sweden. In contrast, Lake Vänern is the largest lake in Sweden. Temperature data have been collected since 1988 (since 1973 for Lake Vänern) between morning and mid-afternoon. Temperature measurements were taken from the deepest point in each lake from the surface through 1 m above the lake bottom at depth intervals varying between 1 m and up to 10 m (in Lake Vänern).Lakes Byglandsfjorden, Hornindalsvatnet, Mjøsa, Øyeren, Selbusjøen, and Strynevatnet are all large and deep lakes located in the central region of western Norway. Temperature profiles were measured in the deepest part of the lake between 9:00 and 16:00 beginning in the mid-1990s. Temperature was measured using an Aanderaa 4060 every meter in the upper part of the water column, with greater than 1 m intervals between measurements in depths below 20 m. For the past 2 to 3 years of data collection, a Castaway CTD recorded temperatures while lowered. Data collection was made by The Norwegian Water Resources and Energy Directorate (NVE).Temperature data were collected from Sweden’s Lake Erken from the deepest point in the lake in 1 m intervals between 7:30 and 9:30 beginning in 1940. Temperature was measured at 1 m intervals using a variety of instruments. In early years, a thermometer inside a transparent Ruttner sampler was read to obtain the temperature of water collected from different depths. Later, underwater thermistors were used from a variety of manufacturers. In recent years, combined temperature and dissolved oxygen sensors have been used to collect water temperature measurements: a YSI model 52 (1996–2006), WTW Oxi 340i (2006–2012), and Hach HQ40d sensor system (2012-present).Lakes Inarijärvi, Kallavesi, Konnevesi, Näsijärvi, Päijänne, Pielinen, and Pyhäjärvi are generally large lakes located throughout Finland, from southern Finland to the northern-most part of the country (Lapland). Lake Pesiöjärvi in the same region has a significantly smaller surface area than others. Lakes Konnevesi and Päijänne have two different temperature profile sites (Konnevesi: Näreselkä and Pynnölänniemi, Päijänne: Linnasaari and Päijätsalo). Temperature data were typically collected in the deepest part of the lake, or in case of large fragmented lakes, the deepest part of the respective basin. Reversing mercury thermometers (with or without Ruttner water samplers) were used until the digital temperature measurements were introduced in the early 1970s. HL Hydrolab Ab PT77A (approximately 1975–1995), DeltaOhm HD8601P (1995–2005), and HT Hydrotechnik Type 110 (2005–2014) have been used for measuring water temperature profiles. However, all device types have been used at each station as long as they have worked. Unfortunately, site-specific documentation of devices used during different years are not available. Before the 1980s, measurements were collected every 5 m before and every 10 m after a depth of 20 m. Since 1980–1981, measurements have been made in 1 m intervals from the surface through 20 m, every 2 m from 20 to 50 m, and every 5 m past 50 m.Lake Pyhäselkä (Pyhaselka) is a large, humic lake located in North Karelia, Finland. It is the northernmost basin of the Saimaa lake system. Temperature data were collected from the deepest point in the lake (Kokonluoto) in 5 m depth intervals beginning in 1962 between 8:00 and 16:00 using a thermometer in a Ruttner water sampler.Central EuropeLakes Annecy, Bourget, and Geneva are located in eastern France37. Temperature data were collected from the deepest point in the lake between 9:30 and 11:00, beginning in 1991 for Lake Annecy, in 1984 for Lake Bourget, and in 1974 for Lake Geneva. In Lake Annecy, various multiparameter probes, including Meerestechnik Elektronik (1991–2001), CTD 90 (2003–2005), CTD 90 M (2008–2011, 2013), and RBR (2012) were used for temperature measurements at depth intervals between 0.01 and 1.8 m. Temperature profiles were measured in Lake Bourget at depth intervals between 0.01 and 10 m using various multiparameter probes, including ISMA probe DNTC (1984–1985), Meerestechnik Elektronik ECO 236 (1986–1998), CTD SBE 19SeaCAT Profiler (1992–2002), and CTD SBE 19plus V2 SeaCAT (2003–2013). In Lake Geneva, water temperature was measured from water samples with a thermometer until 1990 from discrete depths of 5 to 10 m intervals through 50 m depth and approximately 50 m intervals for the rest of the water column. After 1990 various multiparameter probes were used for temperature measurements, including Meerestechnik Elektronik (1991–2001), CTD 90 (2002–2007), CTD 90 M (2008–2011, 2013), and RBR (2012). Data for Annecy, Bourget, and Geneva are available38 at https://data.inrae.fr/dataset.xhtml?persistentId=doi:10.15454/YOLA0Y.The Cumbrian lakes in the English Lake District, Bassenthwaite Lake, Blelham Tarn, Derwent Water, Esthwaite Water, Grasmere, and Windermere North Basin, are located in northwest England. Temperature data were collected in midsummer at the deepest point of each between 9:00 and 14:00 beginning in 1991. Temperature profiles were measured at depth intervals judged in the field depending on the stratification pattern and depth of the lake. Temperature was measured with various combined temperature/oxygen sensors, including a YSI 58 (1991–2002) and WTW Oxi 340 (2002–2013) in Windermere North Basin, and a YSI 58 (1991–2002), WTW Oxi 340i (2002–2010), and Hach HQ 30d and LDO probe (2010–2013) in all other lakes.Lake Constance is located on the border of Germany, Switzerland, and Austria. Temperature profiles were measured at the deepest part of the central basin (Upper Lake Constance) of the lake beginning in 1964, with measurements taken between 9:00 and 10:00 using a thermometer. Depth intervals between samples increased with depth, with measurements taken every 2.5 to 5 m through 20 m, every 10 m through 50 m, and every 50 m down to a depth of 250 m39.Lake Mondsee is located in the Lake District “Salzkammergut” of Austria. Temperature data were collected from the deepest point of the lake approximately monthly near noon beginning in 1968. Temperature measurements were usually measured at depth intervals of 1 to 2 m in the epilimnion and 5 to 10 m intervals in the hypolimnion using a thermometer housed in a Schindler sampler prior to 1998. Beginning in 1998, data were extracted from continuous YSI 6920 profiler readings, with a YSI 6600 used beginning in 2008 and a thermistor chain from 2010–2013.Lake Müggelsee, located in Berlin, Germany was sampled weekly beginning in 1978 between 8:00 and 9:00. Temperature was measured every 0.5 to 1 m at the deepest part of the lake using a Hydrolab H2O sensor beginning in 1992, and a thermistor probe for years prior to 1992.Lake Piburgersee is located in Tyrol, Austria. Temperature was measured using a calibrated thermometer every 3 m throughout the water column in the deepest part of the lake, beginning in 1970. Additional information and analyses are available elsewhere40.Plusssee (Plußsee) is located in northern Germany (Schleswig-Holstein). Temperature data were collected from the deepest point of the lake between 9:00 and 15:00 beginning in 1971. Temperature was measured manually in 1 m increments from 0 to 15 m and in 5 m increments from 15 to 25 m using a thermometer mounted into a Ruttner sampler prior to 1976, and after 1976 with a WTW temperature/dissolved oxygen probe.Traunsee is a large and deep oligotrophic lake in the Salzkammergut lake district of Austria. Temperature data were collected at the deepest point of the lake between 9:00 and 12:00 beginning in 1965. Temperature was measured at 2 to 5 m intervals through 20 m, and at 20 m intervals through the rest of the water column using a mercury thermometer mounted in a 5-liter water sampler.Lower Lake Zurich is located in Switzerland. Temperature profiles were measured at the deepest part of this lake between 8:30 and 12:00 beginning in 1936. Measurements were taken at the deepest point in the lake using a range of sensors, mainly NTC thermistors (1936–2000). Beginning in 2001 various sondes were used, including FLP-10 multisonde (2001–2008), multisonde Hydrolab DS5 (2008–2015) at 0.5 to 1 m intervals through 30 m, 5 m intervals through 50 m, and 10 m intervals throughout the rest of the water column.Southern EuropeLake Garda is one of the largest lakes in Europe, and the largest Italian lake. Owing to its deep depth, Lake Garda is characterized by long periods of incomplete vertical winter water circulation, which are interrupted by full mixing of the water column after the occurrence of harsh winters. Limnological investigations have been carried out since 1991 in a pelagic station located at the point of maximum depth of the northwest basin. Profiles of water temperature were recorded during the summer months using Idronaut Ocean Seven 401 (1991–1997), Seacat SBE 19–03 (1998–2008), and Idronaut Ocean Seven 316Plus since 200925.Lake Iseo is located in northern Italy (Lombardy Region). It is a deep mesotrophic lake, characterized by long periods of incomplete vertical winter water circulation41. Temperature was measured at the deepest point in the lake using an automatic thermistor probe coupled with an oxygen sensor from 1993 to 2011 with Microprocessor Oximeter WTW OXI 320 and from 2012 to 2016 with Microprocessor WTW multi 3410. Temperature was measured for at least ten discrete depths and all measurements were regularly checked with a mercury-filled Celsius reversing thermometer.Lake Lugano is located in the foothills of the Central Alps, on the border between Switzerland and Italy. The lake is divided into a northern and southern basin, which are separated by a causeway (built on a natural moraine). Due to reduced connectivity42 (flow of approximately 0.38 km3 year−1 from north to south) and different morphometric characteristics, the two basins were considered separately in the data set. Temperature profiles were collected at sites near the deepest point of each basin. Temperature was measured using reversing thermometers from 1974–1979 and multiparameter probes thereafter (Hydropolyester HTP 77 during 1980–1985, Ocean Seven 401 during 1986–1993, Ocean Seven 316 during 1994–2015). As an exception, temperatures at depths greater than 100 m were also measured using a reversing thermometer between 1980–1985. Temperature was measured for at least nine or seven discrete depths for the northern and southern basins, respectively, from 1974–1986, whereas full temperature profiles with vertical resolution of 0.5 to 1 m were measured between 1987–2015.Lake Maggiore is a deep lake located in northwestern Italy and Switzerland, south of the Alps. Lake Maggiore can be classified as holo-oligomictic, with complete overturns only occurring at the end of particularly cold and windy winters. Temperature data have been collected at the deepest point of the lake since 1981, usually between 10:00 and 12:00 at the Ghiffa station. Temperature has been measured at discrete depths of 0, 5, 10, 20, 30, 50 m, and every 50 m through 360 m using mercury-filled thermometers connected to the bottle used for water sampling43.Eastern EuropeLakes Batorino, Myastro, and Naroch are located in the northwest part of Belarus, in the glacial landscape. Temperature data were collected monthly in the center of the lake during the vegetative season of May to October beginning in the 1950s and 60s. Measurements of water temperature were taken every 2 to 4 m through the water column between 9:00 and 14:00 with a mercury deepwater thermometer with a scale resolution of 0.1°C mounted in a metal frame, or in a Ruttner sampler.RussiaLake Baikal is located in Siberia, Russia. Temperature data were collected from a station situated 2.8 km from the shoreline near the Bolshie Koty settlement at depths of 0 and 50 m between 9:00 and 12:00 from 1948–2009. Temperature was measured with a mercury thermometer inside a Van Dorn bottle.Lake Glubokoe is located in Central European Russia, Moscow Province. Temperature profile data were collected from the deepest point of the lake from the surface through 10 m. Since 1982, water was taken from the required depth with a 10 L bathometer and its temperature was measured with a mercury thermometer; instrumentation prior to 1982 is unknown.Kurilskoye Lake in Kamchatka, Russia was sampled beginning in 1942. Temperature profiles were measured in the deepest part of the lake between 8:00 and 15:00. Various temperature sensors were used over time, including a reversing thermometer (1942–1965), bathythermograph (1980–2003), Hydrolab (2004–2008), and RINKO profiler (2009–2014).Lake Shira is located in the south of Siberia, Russia. Temperature profiles were measured in the deepest part of the lake between 11:00 and 15:00 from the surface to the depth of 20 to 24 m with various temperature sensors including multisonde Hydrolab 4 A (2000–2008) and a YSI 6600 (2009–2014).Middle EastLake Kinneret is located in Israel (Jordan Valley). Temperature data were collected near the deepest point of the lake (Station A) between 7:00 and 16:00 from 1969–2013, measured every centimeter with an error of ± 0.005 °C, and averaged to every 1 m. Temperature measurements were taken from 1969 to 1986 using an underwater thermometer (Whitney-Montedoro), from 1987 to 2003 using a STD-12 Plus (Applied Microsystems), and from 2003 to 2013 using AML Oceanographic Minos•X.AsiaLake Biwa is located in the central part of the Japanese archipelago (Shiga Prefecture, Japan). Temperature data were collected from a station near the deepest part of the lake from 1958–2010. Temperature was mostly measured between 9:00 and 12:00 at 5 m intervals from 1959–2005, and at 1 m intervals beginning in 2006. Measurements were made using an electronic thermometer (Murayama Denki Ltd.) from 1958–1970, a thermistor thermometer (Shibaura electronics, HCB III) from 1970–1994, a CTD profiler (Alec electronics, ABT-1) from 1994–2006, and a CTD profiler (JFE Advantech, compact-CTD) from 2007–2010.AustraliaLake Burley Griffin is a reservoir constructed in 1963 by damming the Molonglo River. It is located in the geographic center of Canberra, the capital of Australia. Temperature data were collected from the deepest point of the reservoir, near the dam wall. Profiles were measured in 1 m depth intervals (reduced to 3 m intervals in 1992) between 8:45 and 16:15 from 1981–2010 by the National Capital Authority.Lakes Samsonvale (North Pine) and Somerset are located on the east coast of Australia, in southeast Queensland. Temperature in each lake was measured at a site approximately 100 m from the dam wall. Samsonvale’s (North Pine’s) temperature was measured using a YSI 6560 sensor on a YSI 6600 V2 sonde beginning in 2009 continuously over a 24-hour period at 1 m intervals. Prior to 2009, temperature was measured via a thermistor string with various unknown instruments. Somerset’s temperature profiles from 2000–2002 were measured using temperature sensors on a thermistor string, spaced at 0.5 m intervals through 3 m, 1 m intervals through 7 m, 2 m intervals through 17 m, and 3 m intervals for the rest of the water column.New ZealandLakes Brunner and Taupo are located in New Zealand, in the West Coast and Waikato regions, respectively. Temperature profiles were measured at the deepest part of these lakes between 9:00 and 16:00, beginning in the early 1990s. Temperature was measured by lowering CTD profilers through the water column, using an YSI EXO sonde (Brunner) and RBR profiler (Taupo).Lakes Okareka, Okaro, Okataina, Rerewhakaaitu, Rotoehu, Rotoiti, Rotoma, Rotorua, Tarawera, and Tikitapu are located in Rotorua, Bay of Plenty, New Zealand. Temperature profiles were measured between 10:00 and 14:00 in the central basin using Seatech CTD casts with a Seabird 19Plus or 19PlusV2 beginning in 2003. Temperature was measured at a frequency of 4 Hz during each cast and data were binned to 1 m depth intervals. Profiles before 2003 were measured in 1 m depth intervals with either a YSI Water Quality Logger 3800 or YSI Sonde model 3815.AntarcticaLakes Heywood (1962–1995), Moss (1972–2003), and Sombre (1973–2003) are located in the South Orkney Islands, Antarctica. Temperature data from these lakes were collected using a Mackareth-type probe at 1 to 2 m intervals in the deepest part of the lake. More

1.Orams, M. B. Feeding wildlife as a tourism attraction: A review of issues and impacts. Tour. Manag. 23, 281–293 (2002).Article 

Google Scholar 
2.Balmford, A. et al. Walk on the wild side: Estimating the global magnitude of visits to protected areas. PLoS Biol. 13, e1002074 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
3.Knight, J. Making wildlife viewable: Habituation and attraction. Soc. Anim. 17, 167–184 (2009).Article 

Google Scholar 
4.Carter, N. H. et al. Coupled human and natural systems approach to wildlife research and conservation. Ecol. Soc. 19, 43 (2014).Article 

Google Scholar 
5.Balasubramaniam, K. N. et al. Addressing the challenges of research on human–wildlife interactions using the concept of coupled natural & human systems. Biol. Conserv. 257, 109095 (2021).Article 

Google Scholar 
6.Knight, J. The ready-to-view wild monkey: The convenience principle in Japanese wildlife tourism. Ann. Tour. Res. 37, 744–762 (2010).Article 

Google Scholar 
7.Okello, M. M., Manka, S. G. & D’Amour, D. E. The relative importance of large mammal species for tourism in Amboseli National Park, Kenya. Tour. Manag. 29, 751–760 (2008).Article 

Google Scholar 
8.Penteriani, V. et al. Consequences of brown bear viewing tourism: A review. Biol. Conserv. 206, 169–180 (2017).Article 

Google Scholar 
9.Ewen, J. G., Walker, L., Canessa, S. & Groombridge, J. J. Improving supplementary feeding in species conservation. Conserv. Biol. 29, 341–349 (2015).PubMed 
Article 

Google Scholar 
10.Jones, C. G. et al. The restoration of the Mauritius Kestrel Falco punctatus population. Ibis 137, S173–S180 (1995).Article 

Google Scholar 
11.Murray, M. H., Becker, D. J., Hall, R. J. & Hernandez, S. M. Wildlife health and supplemental feeding: A review and management recommendations. Biol. Conserv. 204, 163–174 (2016).Article 

Google Scholar 
12.Oro, D., Genovart, M., Tavecchia, G., Fowler, M. S. & Martínez-Abraín, A. Ecological and evolutionary implications of food subsidies from humans. Ecol. Lett. 16, 1501–1514 (2013).PubMed 
Article 

Google Scholar 
13.Civitello, D. J., Allman, B. E., Morozumi, C. & Rohr, J. R. Assessing the direct and indirect effects of food provisioning and nutrient enrichment on wildlife infectious disease dynamics. Philos. Trans. R. Soc. B Biol. Sci. 373, 20170101 (2018).Article 

Google Scholar 
14.Lappan, S., Malaivijitnond, S., Radhakrishna, S., Riley, E. P. & Ruppert, N. The human–primate interface in the new normal: Challenges and opportunities for primatologists in the COVID-19 era and beyond. Am. J. Primatol. 82, e23176 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Gibb, R. et al. Zoonotic host diversity increases in human-dominated ecosystems. Nature 584, 398–402 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
16.Dunay, E., Apakupakul, K., Leard, S., Palmer, J. L. & Deem, S. L. Pathogen transmission from humans to great apes is a growing threat to primate conservation. EcoHealth 15, 148–162 (2018).PubMed 
Article 

Google Scholar 
17.Fuentes, A., Shaw, E. & Cortes, J. Qualitative assessment of macaque tourist sites in Padangtegal, Bali, Indonesia, and the Upper Rock Nature Reserve, Gibraltar. Int. J. Primatol. 28, 1143–1158 (2007).Article 

Google Scholar 
18.Dellatore, D. F., Waitt, C. D. & Foitovà, I. The impact of tourism on the behavior of rehabilitated orangutans (Pongo abelii) in Bukit Lawang, North Sumatra, Indonesia. In Primate Tourism: A Tool for Conservation (eds Russon, A. E. & Wallis, J.) 98–120 (Cambridge University Press, 2014).Chapter 

Google Scholar 
19.Berman, C. M., Matheson, M. D., Ogawa, H. & Ionica, C. S. Tourism, infant mortality and stress indicators among Tibetan macaques at Huangshan, China. In Primate Tourism: A Tool for Conservation (eds Russon, A. E. & Wallis, J.) 21–43 (Cambridge University Press, 2014).Chapter 

Google Scholar 
20.Kurita, C. M. Provisioning and tourism infree-ranging Japanese macaques. In Primate Tourism: A Tool for Conservation (eds Russon, A. E. & Wallis, J.) 44–55 (Cambridge University Press, 2014).Chapter 

Google Scholar 
21.Long, Y., Bleisch, W. & Richardson, M. Rhinopithecus bieti. The IUCN Red List of Threatened Species 2020: e.T19597A8986243. https://doi.org/10.2305/IUCN.UK.2008.RLTS.T19597A8986243.en. IUCN Red List of Threatened Species https://www.iucnredlist.org/en (2020).22.Li, B., Pan, R. & Oxnard, C. E. Extinction of snub-nosed monkeys in China during the past 400 years. Int. J. Primatol. 23, 1227–1244 (2002).Article 

Google Scholar 
23.Wong, M. H. G., Li, R., Xu, M. & Long, Y. An integrative approach to assessing the potential impacts of climate change on the Yunnan snub-nosed monkey. Biol. Conserv. 158, 401–409 (2013).Article 

Google Scholar 
24.Li, L., Xue, Y., Wu, G., Li, D. & Giraudoux, P. Potential habitat corridors and restoration areas for the black-and-white snub-nosed monkey Rhinopithecus bieti in Yunnan, China. Oryx 49, 719–726 (2015).Article 

Google Scholar 
25.Long, Y., Kirkpatrick, C. R., Zhongtai, & Xiaolin,. Report on the distribution, population, and ecology of the Yunnan snub-nosed monkey (Rhinopithecus bieti). Primates 35, 241–250 (1994).Article 

Google Scholar 
26.Afonso, E. et al. Creating small food-habituated groups might alter genetic diversity in the endangered Yunnan snub-nosed monkey. Glob. Ecol. Conserv. 26, e01422 (2021).Article 

Google Scholar 
27.Cui, Z., Li, J., Chen, Y. & Zhang, L. Molecular epidemiology, evolution, and phylogeny of Entamoeba spp. Infect. Genet. Evol. 75, 104018 (2019).PubMed 
Article 

Google Scholar 
28.Jacob, A. S., Busby, E. J., Levy, A. D., Komm, N. & Clark, C. G. Expanding the Entamoeba universe: New hosts yield novel ribosomal lineages. J. Eukaryot. Microbiol. 63, 69–78 (2016).CAS 
PubMed 
Article 

Google Scholar 
29.Verweij, J. J. et al. Entamoeba histolytica infections in captive primates. Parasitol. Res. 90, 100–103 (2003).PubMed 
Article 

Google Scholar 
30.Tachibana, H. et al. Isolation and characterization of a potentially virulent species Entamoeba nuttalli from captive Japanese macaques. Parasitology 136, 1169–1177 (2009).CAS 
PubMed 
Article 

Google Scholar 
31.Levecke, B. et al. Molecular identification of Entamoeba spp. in captive nonhuman primates. J. Clin. Microbiol. 48, 2988–2990 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Levecke, B. et al. Transmission of Entamoeba nuttalli and Trichuris trichiura from nonhuman primates to humans. Emerg. Infect. Dis. 21, 1871–1872 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
33.Rivera, W. L., Yason, J. A. D. L. & Adao, D. E. V. Entamoeba histolytica and E. dispar infections in captive macaques (Macaca fascicularis) in the Philippines. Primates 51, 69 (2009).PubMed 
Article 

Google Scholar 
34.Regan, C. S., Yon, L., Hossain, M. & Elsheikha, H. M. Prevalence of Entamoeba species in captive primates in zoological gardens in the UK. PeerJ 2, e492 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
35.Elsheikha, H. M., Regan, C. S. & Clark, C. G. Novel Entamoeba findings in nonhuman primates. Trends Parasitol. 34, 283–294 (2018).PubMed 
Article 

Google Scholar 
36.Tuda, J. et al. Identification of Entamoeba polecki with unique 18S rRNA gene sequences from celebes crested macaques and pigs in Tangkoko Nature Reserve, North Sulawesi, Indonesia. J. Eukaryot. Microbiol. 63, 572–577 (2016).CAS 
PubMed 
Article 

Google Scholar 
37.Nolan, M. J. et al. Molecular characterisation of protist parasites in human-habituated mountain gorillas (Gorilla beringei beringei), humans and livestock, from Bwindi Impenetrable National Park, Uganda. Parasit. Vectors 10, 340 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
38.Ruiz-López, M. J., Monello, R. J., Gompper, M. E. & Eggert, L. S. The effect and relative importance of neutral genetic diversity for predicting parasitism varies across parasite taxa. PLoS ONE 7, e45404 (2012).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
39.Acevedo-Whitehouse, K. et al. Contrasting effects of heterozygosity on survival and hookworm resistance in California sea lion pups. Mol. Ecol. 15, 1973–1982 (2006).CAS 
PubMed 
Article 

Google Scholar 
40.Grueter, C. C. et al. Ranging of Rhinopithecus bieti in the Samage Forest, China. I. Characteristics of range use. Int. J. Primatol. 29, 1121–1145 (2008).Article 

Google Scholar 
41.Li, D. et al. Ranging of Rhinopithecus bieti in the Samage Forest, China. II. Use of land cover types and altitudes. Int. J. Primatol. 29, 1147 (2008).Article 

Google Scholar 
42.Xue, Y. et al. Analysis of habitat connectivity of the Yunnan snub-nosed monkeys (Rhinopithecus bieti) using landscape genetics. Shengtai Xuebao Acta Ecol. Sin. 31, 5886–5893 (2011).ADS 

Google Scholar 
43.Fu, R., Li, L., Yu, Z., Afonso, E. & Giraudoux, P. Spatial and temporal distribution of Yunnan snub-nosed monkey, Rhinopithecus bieti, indices. Mammalia 83, 103 (2018).Article 

Google Scholar 
44.Vlčková, K. et al. Diversity of Entamoeba spp. in African great apes and humans: An insight from Illumina MiSeq high-throughput sequencing. Int. J. Parasitol. 48, 519–530 (2018).PubMed 
Article 
CAS 

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

Google Scholar 
46.Paradis, E. & Schliep, K. ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528 (2019).CAS 
PubMed 
Article 

Google Scholar 
47.Pagès, H., Aboyoun, P., Gentleman, R. & DebRoy, S. Biostrings: Efficient manipulation of biological strings. (2017).48.Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Wright, E. S., Yilmaz, L. S. & Noguera, D. R. DECIPHER, a search-based approach to chimera identification for 16S rRNA sequences. Appl. Environ. Microbiol. 78, 717–725 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Schliep, K. P. phangorn: Phylogenetic analysis in R. Bioinformatics 27, 592–593 (2011).CAS 
PubMed 
Article 

Google Scholar 
51.McMurdie, P. J. & Holmes, S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8, e61217 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Morgan, M. et al. ShortRead: A bioconductor package for input, quality assessment and exploration of high-throughput sequence data. Bioinformatics 25, 2607–2608 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Galan, M. et al. 16S rRNA Amplicon sequencing for epidemiological surveys of bacteria in wildlife. mSystems 1, e00032 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
54.Smith, D. P. & Peay, K. G. Sequence depth, not PCR replication, improves ecological inference from next generation DNA sequencing. PLoS ONE 9, e90234 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
55.Stensvold, C. R. et al. Increased sampling reveals novel lineages of Entamoeba: Consequences of genetic diversity and host specificity for taxonomy and molecular detection. Protist 162, 525–541 (2011).PubMed 
Article 

Google Scholar 
56.Burnham, K. P. & Anderson, D. R. Data-based selection of an appropriate biological model: The key to modern data analysis. In Wildlife 2001: Populations (eds McCullough, D. R. & Barrett, R. H.) 16–30 (Springer, 2001). https://doi.org/10.1007/978-94-011-2868-1_3.Chapter 

Google Scholar 
57.Oksanen, J. et al. vegan: Community Ecology Package. R package version 2.5-6. (2019).58.Matsubayashi, M. et al. First detection and molecular identification of Entamoeba bovis from Japanese cattle. Parasitol. Res. 117, 339–342 (2018).PubMed 
Article 

Google Scholar 
59.Balloux, F., Amos, W. & Coulson, T. Does heterozygosity estimate inbreeding in real populations?. Mol. Ecol. 13, 3021–3031 (2004).CAS 
PubMed 
Article 

Google Scholar 
60.Szulkin, M., Bierne, N. & David, P. Heterozygosity-fitness correlations: A time for reappraisal. Evolution 64, 1202–1217 (2010).PubMed 

Google Scholar 
61.Feng, M. et al. Prevalence and genetic diversity of Entamoeba species infecting macaques in southwest China. Parasitol. Res. 112, 1529–1536 (2013).PubMed 
Article 

Google Scholar 
62.Guan, Y. et al. Comparative analysis of genotypic diversity in Entamoeba nuttalli isolates from Tibetan macaques and rhesus macaques in China. Infect. Genet. Evol. 38, 126–131 (2016).CAS 
PubMed 
Article 

Google Scholar 
63.Ponce Gordo, F., Martı́nez Dı́az, R. A. & Herrera, S. Entamoeba struthionis n.sp. (Sarcomastigophora: Endamoebidae) from ostriches (Struthio camelus). Vet. Parasitol. 119, 327–335 (2004).CAS 
PubMed 
Article 

Google Scholar 
64.Ai, S. et al. The first survey and molecular identification of Entamoeba spp. in farm animals on Qinghai-Tibetan Plateau of China. Comp. Immunol. Microbiol. Infect. Dis. 75, 101607 (2021).PubMed 
Article 

Google Scholar 
65.Stensvold, C. R., Lebbad, M. & Clark, C. G. Genetic characterisation of uninucleated cyst-producing Entamoeba spp. from ruminants. Int. J. Parasitol. 40, 775–778 (2010).CAS 
PubMed 
Article 

Google Scholar 
66.Fadrosh, D. W. et al. An improved dual-indexing approach for multiplexed 16S rRNA gene sequencing on the Illumina MiSeq platform. Microbiome 2, 6 (2014).PubMed 
PubMed Central 
Article 

Google Scholar  More

1.Hutchinson, G. E. Cold Spring Harbor Symposia on Quantitative Biology. Concluding Remarks 22 415–427 (1957).2.Smith, E. P. Niche breadth, resource availability, and inference. Ecology 63, 1675–1681. https://doi.org/10.2307/1940109 (1982).Article 

Google Scholar 
3.Leibold, M. A. The niche concept revisited: Mechanistic models and community context. Ecology 76, 1371–1382. https://doi.org/10.2307/1938141 (1995).Article 

Google Scholar 
4.Sexton, J. P., Montiel, J., Shay, J. E., Stephens, M. R. & Slatyer, R. A. Evolution of ecological niche breadth. Annu. Rev. Ecol. Evol. Syst. 48, 183–206. https://doi.org/10.1146/annurev-ecolsys-110316-023003 (2017).Article 

Google Scholar 
5.Jaenike, J. Host specialization in phytophagous insects. Annu. Rev. Ecol. Syst. 21, 243–273. https://doi.org/10.1146/annurev.es.21.110190.001331 (1990).Article 

Google Scholar 
6.Thompson, J. N. The Coevolutionary Process (University of Chicago Press, 1994).Book 

Google Scholar 
7.Krasnov, B. R., Mouillot, D., Shenbrot, G. I., Khokhlova, I. S. & Poulin, R. Geographical variation in host specificity of fleas (Siphonaptera) parasitic on small mammals: The influence of phylogeny and local environmental conditions. Ecography 27, 787–797. https://doi.org/10.1111/j.0906-7590.2004.04015.x (2004).Article 

Google Scholar 
8.Poullain, V., Gandon, S., Brockhurst, M. A., Buckling, A. & Hochberg, M. E. The evolution of specificity in evolving and coevolving antagonistic interactions between bacteria and its phage. Evolution 62, 1–11. https://doi.org/10.1111/j.1558-5646.2007.00260.x (2008).Article 
PubMed 

Google Scholar 
9.Whitlock, M. C. The Red Queen beats the Jack-Of-All-Trades: The limitations on the evolution of phenotypic plasticity and niche breadth. Am. Nat. 148, S65–S77. https://doi.org/10.1086/285902 (1996).Article 

Google Scholar 
10.Gandon, S. Local adaptation and the geometry of host–parasite coevolution. Ecol. Lett. 5, 246–256. https://doi.org/10.1046/j.1461-0248.2002.00305.x (2002).Article 

Google Scholar 
11.Alizon, S. & Michalakis, Y. Adaptive virulence evolution: The good old fitness-based approach. Trends Ecol. Evol. 30, 248–254. https://doi.org/10.1016/j.tree.2015.02.009 (2015).Article 
PubMed 

Google Scholar 
12.Frank, S. A. & Schmid-Hempel, P. Mechanisms of pathogenesis and the evolution of parasite virulence. J. Evol. Biol. 21, 396–404. https://doi.org/10.1111/j.1420-9101.2007.01480.x (2008).CAS 
Article 
PubMed 

Google Scholar 
13.Beadell, J. S. et al. Global phylogeographic limits of Hawaii’s avian malaria. Proc. R. Soc. B: Biol. Sci. 273, 2935–2944. https://doi.org/10.1098/rspb.2006.3671 (2006).Article 

Google Scholar 
14.Krasnov, B. R. Functional and Evolutionary Ecology of Fleas: A Model for Ecological Parasitology (Cambridge University Press, 2008).Book 

Google Scholar 
15.Poulin, R. Evolutionary Ecology of Parasites (Princeton University Press, 2011).Book 

Google Scholar 
16.Välimäki, P. et al. Geographical variation in host use of a blood-feeding ectoparasitic fly: Implications for population invasiveness. Oecologia 166, 985–995. https://doi.org/10.1007/s00442-011-1951-y (2011).ADS 
Article 
PubMed 

Google Scholar 
17.Theodosopoulos, A. N., Hund, A. K. & Taylor, S. A. Parasites and host species barriers in animal hybrid zones. Trends Ecol. Evol. 34, 19–30. https://doi.org/10.1016/j.tree.2018.09.011 (2019).Article 
PubMed 

Google Scholar 
18.Mackenzie, A. A trade-off for host plant utilization in the black bean aphid, Aphis fabae. Evolution 50, 155–162. https://doi.org/10.1111/j.1558-5646.1996.tb04482.x (1996).Article 
PubMed 

Google Scholar 
19.Harrington, L. C., Edman, J. D. & Scott, T. W. Why do female Aedes aegypti (Diptera: Culicidae) feed preferentially and frequently on human blood?. J. Med. Entomol. 38, 411–422. https://doi.org/10.1603/0022-2585-38.3.411 (2001).CAS 
Article 
PubMed 

Google Scholar 
20.Dick, C. W. & Patterson, B. D. Against all odds: Explaining high host specificity in dispersal-prone parasites. Int. J. Parasitol. 37, 871–876. https://doi.org/10.1016/j.ijpara.2007.02.004 (2007).Article 
PubMed 

Google Scholar 
21.Torchin, M. E. & Mitchell, C. E. Parasites, pathogens, and invasions by plants and animals. Front. Ecol. Environ. 2, 183–190. https://doi.org/10.1890/1540-9295(2004)002[0183:PPAIBP]2.0.CO;2 (2004).Article 

Google Scholar 
22.Clark, N. J. & Clegg, S. M. The influence of vagrant hosts and weather patterns on the colonization and persistence of blood parasites in an island bird. J. Biogeogr. 42, 641–651. https://doi.org/10.1111/jbi.12454 (2015).Article 

Google Scholar 
23.Kawecki, T. J. Red Queen meets Santa Rosalia: Arms races and the evolution of host specialization in organisms with parasitic lifestyles. Am. Nat. 152, 635–651. https://doi.org/10.1086/286195 (1998).CAS 
Article 
PubMed 

Google Scholar 
24.Egas, M., Dieckmann, U. & Sabelis, M. W. Evolution restricts the coexistence of specialists and generalists: The role of trade-off structure. Am. Nat. 163, 518–531. https://doi.org/10.1086/382599 (2004).Article 
PubMed 

Google Scholar 
25.Poulin, R. & Keeney, D. B. Host specificity under molecular and experimental scrutiny. Trends Parasitol. 24, 24–28. https://doi.org/10.1016/j.pt.2007.10.002 (2008).CAS 
Article 
PubMed 

Google Scholar 
26.Lyimo, I. N. & Ferguson, H. M. Ecological and evolutionary determinants of host species choice in mosquito vectors. Trends Parasitol. 25, 189–196. https://doi.org/10.1016/j.pt.2009.01.005 (2009).Article 
PubMed 

Google Scholar 
27.Visher, E. & Boots, M. The problem of mediocre generalists: Population genetics and eco-evolutionary perspectives on host breadth evolution in pathogens. Proc. R. Soc. B: Biol. Sci. 287, 20201230. https://doi.org/10.1098/rspb.2020.1230 (2020).Article 

Google Scholar 
28.Sarfati, M. et al. Energy costs of blood digestion in a host-specific haematophagous parasite. J. Exp. Biol. 208, 2489. https://doi.org/10.1242/jeb.01676 (2005).Article 
PubMed 

Google Scholar 
29.Fry, J. D. The evolution of host specialization: Are trade-offs overrated?. Am. Nat. 148, S84–S107. https://doi.org/10.1086/285904 (1996).Article 

Google Scholar 
30.Fessl, B. et al. Galápagos landbirds (passerines, cuckoos, and doves): Status, threats, and knowledge gaps. Galápagos Rep. 2016, 149 (2015).
Google Scholar 
31.Fessl, B., Heimpel, G. E. & Causton, C. E. Invasion of an avian nest parasite, Philornis downsi, to the Galapagos Islands: colonization history, adaptations to novel ecosystems, and conservation challenges. In Disease Ecology: Galapagos Birds and their Parasites (ed Patricia G. Parker) 213–266 (Springer International Publishing, 2018).32.Frankham, R. Do island populations have less genetic variation than mainland populations?. Heredity 78, 311–327. https://doi.org/10.1038/hdy.1997.46 (1997).Article 
PubMed 

Google Scholar 
33.Reichard, M. et al. The bitterling–mussel coevolutionary relationship in areas of recent and ancient sympatry. Evolution 64, 3047–3056. https://doi.org/10.1111/j.1558-5646.2010.01032.x (2010).Article 
PubMed 

Google Scholar 
34.Wiedenfeld, D. A., Jiménez, G. U., Fessl, B., Kleindorfer, S. & Carlos Valarezo, J. Distribution of the introduced parasitic fly Philornis downsi (Diptera, Muscidae) in the Galápagos Islands. Pacific Conserv. Biol. 13, 14–19. https://doi.org/10.1071/PC070014 (2007).Article 

Google Scholar 
35.Fessl, B., Sinclair, B. J. & Kleindorfer, S. The life-cycle of Philornis downsi (Diptera: Muscidae) parasitizing Darwin’s finches and its impacts on nestling survival. Parasitology 133, 739–747. https://doi.org/10.1017/S0031182006001089 (2006).CAS 
Article 
PubMed 

Google Scholar 
36.Kleindorfer, S. & Dudaniec, R. Y. Host-parasite ecology, behavior and genetics: A review of the introduced fly parasite Philornis downsi and its Darwin’s finch hosts. BMC Zool. 1, 1. https://doi.org/10.1186/s40850-016-0003-9 (2016).Article 

Google Scholar 
37.Galligan, T. H. & Kleindorfer, S. Naris and beak malformation caused by the parasitic fly, Philornis downsi (Diptera: Muscidae), in Darwin’s small ground finch, Geospiza fuliginosa (Passeriformes: Emberizidae). Biol. J. Lin. Soc. 98, 577–585. https://doi.org/10.1111/j.1095-8312.2009.01309.x (2009).Article 

Google Scholar 
38.Kleindorfer, S., Custance, G., Peters Katharina, J. & Sulloway Frank, J. Introduced parasite changes host phenotype, mating signal and hybridization risk: Philornis downsi effects on Darwin’s finch song. Proc. R. Soc. B: Biol. Sci. 286, 20190461. https://doi.org/10.1098/rspb.2019.0461 (2019).Article 

Google Scholar 
39.Kleindorfer, S., Peters, K. J., Custance, G., Dudaniec, R. Y. & O’Connor, J. A. Changes in Philornis infestation behavior threaten Darwin’s finch survival. Curr. Zool. 60, 542–550. https://doi.org/10.1093/czoolo/60.4.542 (2014).Article 

Google Scholar 
40.O’Connor, J. A., Sulloway, F. J., Robertson, J. & Kleindorfer, S. Philornis downsi parasitism is the primary cause of nestling mortality in the critically endangered Darwin’s medium tree finch (Camarhynchus pauper). Biodivers. Conserv. 19, 853–866. https://doi.org/10.1007/s10531-009-9740-1 (2010).Article 

Google Scholar 
41.Knutie, S. A. et al. Galápagos mockingbirds tolerate introduced parasites that affect Darwin’s finches. Ecology https://doi.org/10.1890/15-0119 (2016).Article 
PubMed 

Google Scholar 
42.Peters, K. J., Evans, C., Aguirre, J. D. & Kleindorfer, S. Genetic admixture predicts parasite intensity: Evidence for increased hybrid performance in Darwin’s tree finches. R. Soc. Open Sci. 6, 181616. https://doi.org/10.1098/rsos.181616 (2019).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
43.Kleindorfer, S. The ecology of clutch size variation in Darwin’s Small Ground Finch Geospiza fuliginosa: Comparison between lowland and highland habitats. Ibis 149, 730–741. https://doi.org/10.1111/j.1474-919X.2007.00694.x (2007).Article 

Google Scholar 
44.Fessl, B. & Tebbich, S. Philornis downsi– a recently discovered parasite on the Galápagos archipelago: A threat for Darwin’s finches?. Ibis 144, 445–451. https://doi.org/10.1046/j.1474-919X.2002.00076.x (2002).Article 

Google Scholar 
45.Dudaniec, R. Y., Fessl, B. & Kleindorfer, S. Interannual and interspecific variation in intensity of the parasitic fly, Philornis downsi, Darwin’s finches. Biol. Cons. 139, 325–332. https://doi.org/10.1016/j.biocon.2007.07.006 (2007).Article 

Google Scholar 
46.Cimadom, A. et al. Invasive parasites, habitat change and heavy rainfall reduce breeding success in Darwin’s Finches. PLoS ONE 9, e107518. https://doi.org/10.1371/journal.pone.0107518 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
47.Cimadom, A. et al. Weed management increases the detrimental effect of an invasive parasite on arboreal Darwin’s finches. Biol. Cons. 233, 93–101. https://doi.org/10.1016/j.biocon.2019.02.025 (2019).Article 

Google Scholar 
48.Kleindorfer, S. & Dudaniec, R. Y. Love thy neighbour? Social nesting pattern, host mass and nest size affect ectoparasite intensity in Darwin’s tree finches. Behav. Ecol. Sociobiol. 63, 731–739. https://doi.org/10.1007/s00265-008-0706-1 (2009).Article 

Google Scholar 
49.Common, L. K., Dudaniec, R. Y., Colombelli-Négrel, D. & Kleindorfer, S. Taxonomic shifts in Philornis larval behaviour and rapid changes in Philornis downsi Dodge & Aitken (Diptera: Muscidae): An invasive avian parasite on the Galápagos Islands. in Life Cycle and Development of Diptera (ed Muhammad Sarwar) (IntechOpen, 2019).50.McNew, S. M. et al. Annual environmental variation influences host tolerance to parasites. Proc. R. Soc. B: Biol. Sci. 286, 20190049. https://doi.org/10.1098/rspb.2019.0049 (2019).CAS 
Article 

Google Scholar 
51.McNew, S. M. & Clayton, D. H. Alien invasion: Biology of Philornis flies highlighting Philornis downsi, an introduced parasite of Galápagos birds. Annu. Rev. Entomol. 63, 369–387. https://doi.org/10.1146/annurev-ento-020117-043103 (2018).CAS 
Article 
PubMed 

Google Scholar 
52.Kleindorfer, S. & Dudaniec, R. Y. Hybridization fluctuates with rainfall in Darwin’s tree finches. Biol. J. Lin. Soc. 130, 79–88. https://doi.org/10.1093/biolinnean/blaa029 (2020).Article 

Google Scholar 
53.Peters, K. J., Myers, S. A., Dudaniec, R. Y., O’Connor, J. A. & Kleindorfer, S. Females drive asymmetrical introgression from rare to common species in Darwin’s tree finches. J. Evol. Biol. 30, 1940–1952. https://doi.org/10.1111/jeb.13167 (2017).CAS 
Article 
PubMed 

Google Scholar 
54.Kleindorfer, S. et al. Species collapse via hybridization in Darwin’s Tree Finches. Am. Nat. 183, 325–341. https://doi.org/10.1086/674899 (2014).Article 
PubMed 

Google Scholar 
55.Loo, W. T., Dudaniec, R. Y., Kleindorfer, S. & Cavanaugh, C. M. An inter-island comparison of Darwin’s finches reveals the impact of habitat, host phylogeny, and island on the gut microbiome. PLoS ONE 14, e0226432. https://doi.org/10.1371/journal.pone.0226432 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
56.Galapagos Conservancy. Galapagos Vital Signs: A satellite-based environmental monitoring system for the Galapagos Archipelago, https://galapagosvitalsigns.org (2021).57.Couri, M. Considerações sobre as relações ecológicas das larvas de Philornis Meinert, 1890 (Diptera, Muscidae) com aves. Revista Brasileira de Entomologia 29, 17–20. https://doi.org/10.1017/S0031182006001089 (1985).Article 

Google Scholar 
58.Skidmore, P. The Biology of the Muscidae of the World Vol. 29 (Springer, 1985).
Google Scholar 
59.O’Connor, J. A., Robertson, J. & Kleindorfer, S. Video analysis of host–parasite interactions in nests of Darwin’s finches. Oryx 44, 588–594. https://doi.org/10.1017/S0030605310000086 (2010).Article 

Google Scholar 
60.O’Connor, J. A., Robertson, J. & Kleindorfer, S. Darwin’s finch begging intensity does not honestly signal need in parasitised nests. Ethology 120, 228–237. https://doi.org/10.1111/eth.12196 (2014).Article 

Google Scholar 
61.Kleindorfer, S. & Sulloway, F. J. Naris deformation in Darwin’s finches: Experimental and historical evidence for a post-1960s arrival of the parasite Philornis downsi. Glob. Ecol. Conserv. 7, 122–131. https://doi.org/10.1016/j.gecco.2016.05.006 (2016).Article 

Google Scholar 
62.Lahuatte, P. F., Lincango, M. P., Heimpel, G. E. & Causton, C. E. Rearing larvae of the avian nest parasite, Philornis downsi (Diptera: Muscidae), on chicken blood-based diets. J. Insect Sci. https://doi.org/10.1093/jisesa/iew064 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
63.Kleindorfer, S. Nesting success in Darwin’s small tree finch, Camarhynchus parvulus: Evidence of female preference for older males and more concealed nests. Anim. Behav. 74, 795–804. https://doi.org/10.1016/j.anbehav.2007.01.020 (2007).Article 

Google Scholar 
64.Nijhout, H. F. & Callier, V. Developmental mechanisms of body size and wing-body scaling in insects. Annu. Rev. Entomol. 60, 141–156. https://doi.org/10.1146/annurev-ento-010814-020841 (2015).CAS 
Article 
PubMed 

Google Scholar 
65.Singh, D. & Bala, M. The effect of starvation on the larval behavior of two forensically important species of blow flies (Diptera: Calliphoridae). For. Sci. Int. 193, 118–121. https://doi.org/10.1016/j.forsciint.2009.09.022 (2009).Article 

Google Scholar 
66.Coulson, S. J. & Bale, J. S. Characterisation and limitations of the rapid cold-hardening response in the housefly Musca domestica (Diptera: Muscidae). J. Insect Physiol. 36, 207–211. https://doi.org/10.1016/0022-1910(90)90124-X (1990).Article 

Google Scholar 
67.R Core Team. R: A language and environment for statistical computing. R version 4.0.3 (R Foundation for Statistical Computing, Vienna, Austria, 2020).68.Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. https://doi.org/10.18637/jss.v067.i01 (2015).Article 

Google Scholar 
69.Venables, B. & Ripley, B. Modern Applied Statistics with S-PLUS (Springer Science & Business Media, 2002).70.Fox, J. & Weisberg, S. An R Companion to Applied Regression (Sage publications, 2011).
Google Scholar 
71.Sarkar, D. Lattice: Multivariate Data Visualization with R (Springer, 2008).Book 

Google Scholar 
72.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2009).Book 

Google Scholar 
73.Fox, J. Effect displays in R for generalised linear models. J. Stat. Softw. https://doi.org/10.18637/jss.v008.i15 (2003).Article 

Google Scholar 
74.Burnham, K. P. & Anderson, D. R. Model Selection and Multimodal Inference: A Practical Information-Theoretic Approach (eds Kenneth P. Burnham & David R. Anderson) 75–117 (Springer New York, 1998).75.Grueber, C. E., Nakagawa, S., Laws, R. J. & Jamieson, I. G. Multimodel inference in ecology and evolution: Challenges and solutions. J. Evol. Biol. 24, 699–711. https://doi.org/10.1111/j.1420-9101.2010.02210.x (2011).CAS 
Article 
PubMed 

Google Scholar 
76.Haaland, T. R., Wright, J. & Ratikainen, I. I. Generalists versus specialists in fluctuating environments: A bet-hedging perspective. Oikos 129, 879–890. https://doi.org/10.1111/oik.07109 (2020).Article 

Google Scholar 
77.Davies, N. Cuckoos, Cowbirds and Other Cheats (Bloomsbury Publishing, 2010).
Google Scholar 
78.Dudaniec, R. Y., Gardner, M. G. & Kleindorfer, S. Offspring genetic structure reveals mating and nest infestation behaviour of an invasive parasitic fly (Philornis downsi) of Galápagos birds. Biol. Invas. 12, 581–592. https://doi.org/10.1007/s10530-009-9464-x (2010).Article 

Google Scholar 
79.Fredensborg, B. L. & Poulin, R. Larval helminths in intermediate hosts: Does competition early in life determine the fitness of adult parasites?. Int. J. Parasitol. 35, 1061–1070. https://doi.org/10.1016/j.ijpara.2005.05.005 (2005).CAS 
Article 
PubMed 

Google Scholar 
80.Begon, M., Harper, J. L. & Townsend, C. R. Ecology: Individuals, Populations and Communities (Blackwell Scientific Publications, 1986).
Google Scholar 
81.Fraik, A. K. et al. Disease swamps molecular signatures of genetic-environmental associations to abiotic factors in Tasmanian devil (Sarcophilus harrisii) populations. Evolution 74, 1392–1408. https://doi.org/10.1111/evo.14023 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
82.Dvorak, M. et al. Conservation status of landbirds on Floreana: The smallest inhabited Galápagos Island. J. Field Ornithol. 88, 132–145. https://doi.org/10.1111/jofo.12197 (2017).Article 

Google Scholar 
83.Hedrick, P. W., Kim, T. J. & Parker, K. M. Parasite resistance and genetic variation in the endangered Gila topminnow. Anim. Conserv. 4, 103–109. https://doi.org/10.1017/S1367943001001135 (2001).Article 

Google Scholar 
84.Lewontin, R. C. & Birch, L. C. Hybridization as a source of variation for adaptation to new environments. Evolution 20, 315–336. https://doi.org/10.2307/2406633 (1966).CAS 
Article 
PubMed 

Google Scholar 
85.Wolinska, J., Lively, C. M. & Spaak, P. Parasites in hybridizing communities: The Red Queen again?. Trends Parasitol. 24, 121–126. https://doi.org/10.1016/j.pt.2007.11.010 (2008).Article 
PubMed 

Google Scholar 
86.Floate, K. D. & Whitham, T. G. The, “Hybrid Bridge” Hypothesis: Host shifting via plant hybrid swarms. Am. Nat. 141, 651–662. https://doi.org/10.1086/285497 (1993).CAS 
Article 
PubMed 

Google Scholar 
87.Le Brun, N., Renaud, F., Berrebi, P. & Lambert, A. Hybrid zones and host-parasite relationships: Effect on the evolution of parasitic specificity. Evolution 46, 56–61. https://doi.org/10.1111/j.1558-5646.1992.tb01984.x (1992).Article 
PubMed 

Google Scholar 
88.Fritz, R. S., Moulia, C. & Newcombe, G. Resistance of hybrid plants and animals to herbivores, pathogens, and parasites. Annu. Rev. Ecol. Syst. 30, 565–591. https://doi.org/10.1146/annurev.ecolsys.30.1.565 (1999).Article 

Google Scholar 
89.Moulia, C., Brun, N. L., Loubes, C., Marin, R. & Renaud, F. Hybrid vigour against parasites in interspecific crosses between two mice species. Heredity 74, 48–52. https://doi.org/10.1038/hdy.1995.6 (1995).Article 
PubMed 

Google Scholar 
90.Gibson, A. K., Refrégier, G., Hood, M. E. & Giraud, T. Performance of a hybrid fungal pathogen on pure-species and hybrid host plants. Int. J. Plant Sci. 175, 724–730. https://doi.org/10.1086/676621 (2014).Article 

Google Scholar 
91.Arnold, M. L. & Martin, N. H. Hybrid fitness across time and habitats. Trends Ecol. Evol. 25, 530–536. https://doi.org/10.1016/j.tree.2010.06.005 (2010).Article 
PubMed 

Google Scholar 
92.Ben-Yosef, M. et al. Host-specific associations affect the microbiome of Philornis downsi, an introduced parasite to the Galápagos Islands. Mol. Ecol. 26, 4644–4656. https://doi.org/10.1111/mec.14219 (2017).Article 
PubMed 

Google Scholar 
93.Loo, W. T., García-Loor, J., Dudaniec, R. Y., Kleindorfer, S. & Cavanaugh, C. M. Host phylogeny, diet, and habitat differentiate the gut microbiomes of Darwin’s finches on Santa Cruz Island. Sci. Rep. 9, 18781. https://doi.org/10.1038/s41598-019-54869-6 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
94.Knutie, S. A. Relationships among introduced parasites, host defenses, and gut microbiota of Galapagos birds. Ecosphere 9, e02286. https://doi.org/10.1002/ecs2.2286 (2018).Article 

Google Scholar 
95.Knutie, S. A., Chaves, J. A. & Gotanda, K. M. Human activity can influence the gut microbiota of Darwin’s finches in the Galapagos Islands. Mol. Ecol. 28, 2441–2450. https://doi.org/10.1111/mec.15088 (2019).Article 
PubMed 

Google Scholar  More

1.Barringer, L. E., Donovall, L. R., Spichiger, S. E., Lynch, D. & Henry, D. The first New World record of Lycorma delicatula (Insecta: Hemiptera: Fulgoridae). Entomol. News 125, 20–23 (2015).Article 

Google Scholar 
2.Dara, S. K., Barringer, L. & Arthurs, S. P. Lycorma delicatula (Hemiptera: Fulgoridae): A new invasive pest in the United States. J. Integr. Pest Manag. 6, 20 (2015).Article 

Google Scholar 
3.Jung, J. M., Jung, S., Byeon, D. & Lee, W. Model-based prediction of potential distribution of the invasive insect pest, spotted lanternfly Lycorma delicatula (Hemiptera: Fulgoridae), by using CLIMEX. J. Asia-Pac. Biodivers. 10, 532–538 (2017).Article 

Google Scholar 
4.Lee, D.-H., Park, Y.-L. & Leskey, T. C. A review of biology and management of Lycorma delicatula (Hemiptera: Fulgoridae), an emerging global invasive species. J. Asia-Pac. Entomol. 22, 589–596 (2019).Article 

Google Scholar 
5.(NYSIPM) New York State Integrated Pest Management. 2020. Spotted lanternfly. https://nysipm.cornell.edu/environment/invasive-species-exotic-pests/spotted-lanternfly/. Accessed 18 January 2021.6.Park, M., Kim, S. M. & Lee, J. H. Genetic structure of Lycorma delicatula (Hemiptera: Fulgoridae) populations in Korea: implication for invasion processes in heterogeneous landscapes. Bull. Entomol. Res. 103, 414–424 (2013).CAS 
Article 

Google Scholar 
7.Keller, J. A. et al. Dispersal of Lycorma delicatula (Hemiptera: Fulgoridae) nymphs through contiguous, deciduous forest. Environ. Entomol. 49, 1012–1018 (2020).Article 

Google Scholar 
8.Smyers, E. C. et al. Spatio-temporal model for predicting spring hatch of the spotted lanternfly (Hemiptera: Fulgoridae). Environ. Entomol. 50, 126–137 (2020).Article 

Google Scholar 
9.Wakie, T. T., Neven, L. G., Yee, W. L. & Lu, Z. The establishment risk of Lycorma delicatula (Hemiptera: Fulgoridae) in the United States and globally. J. Econ. Entomol. 113, 306–314 (2019).
Google Scholar 
10.Urban, J. M. Perspective: shedding light on spotted lanternfly impacts in the USA. Pest Manag. Sci. 76, 10–17 (2020).CAS 
Article 

Google Scholar 
11.Harper, J. K., Stone, W., Kelsey, T. W. & Kime, L. F. Potential Economic Impact of the Spotted Lanternfly on Agriculture and Forestry in Pennsylvania (The Center for Rural Pennsylvania, 2019).
Google Scholar 
12.Song, M. K. Damage by Lycorma delicatula and chemical control in vineyards. M.S. thesis. Chunbuk National University, Korea (2010).13.Tedders, W. L. & Smith, J. S. Shading effect on pecan by sooty mold growth. J. Econ. Entomol. 69, 551–553 (1976).Article 

Google Scholar 
14.Lemos-Filho, J. P. D. & Paiva, É. A. S. The effects of sooty mold on photosynthesis and mesophyll structure of mahogany (Swietenia macrophylla King., Meliaceae). Bragantia 65, 11–17 (2006).Article 

Google Scholar 
15.Han, J. M. et al. Lycorma delicatula (Hemiptera: Auchenorrhyncha: Fulgoridae: Aphaeninae) finally, but suddenly arrived in Korea. Entomol. Res. 38, 281–286 (2008).Article 

Google Scholar 
16.Park, J. D. et al. Biological characteristics of Lycorma delicatula and the control effects of some insecticides. Korean J. Appl. Entomol. 48, 53–57 (2009).Article 

Google Scholar 
17.Liu, H. Oviposition substrate selection, egg mass characteristics, host preference, and life history of the spotted lanternfly (Hemiptera: Fulgoridae) in North America. Environ. Entomol. 48, 1452–1468 (2019).Article 

Google Scholar 
18.Barringer, L. E. & Ciafré, C. M. Worldwide feeding host plants of spotted lanternfly, with significant additions from North America. Environ. Entomol. 49, 999–1011 (2020).Article 

Google Scholar 
19.Uyi, O. et al. Spotted lanternfly (Hemiptera: Fulgoridae) can complete development and reproduce without access to the preferred host, Ailanthus altissima. Environ. Entomol. 49, 1185–1190 (2020).Article 

Google Scholar 
20.Murman, K. Distribution, survival, and development of spotted lanternfly on host plants found in north America. Environ. Entomol. 49, 1270–1281 (2020).Article 

Google Scholar 
21.Magnusson, A. glmmTMB: Generalized linear mixed models using template model builder. R package v. 0.1.3. https://github.com/glmmTMB (2017).22.R Development Core Team. R: A Language and Environment for Statistical Computing Computer Program, Version 3.6.3 (R Development Core Team, 2020).
Google Scholar 
23.Kariyat, R. R. & Portman, S. L. Plant–herbivore interactions: Thinking beyond larval growth and mortality. Am. J. Bot. 103, 1–3 (2016).Article 

Google Scholar 
24.Fordyce, J. A. & Shapiro, A. M. Another perspective on the slow-growth/high-mortality hypothesis: chilling effects on swallowtail larvae. Ecology 84, 263–268 (2003).Article 

Google Scholar 
25.Uesugi, A. The slow-growth high-mortality hypothesis: direct experimental support in a leaf mining fly. Ecol. Entomol. 40, 221–228 (2015).Article 

Google Scholar 
26.Song, S., Kim, S., Kwon, S. W., Lee, S.-I. & Jablonski, P. G. Defense sequestration associated with narrowing of diet and ontogenetic change to aposematic colours in the spotted lanternfly. Sci. Rep. 8, 16831 (2018).ADS 
Article 

Google Scholar 
27.Domingue, M. J. & Baker, T. C. Orientation of flight for physically disturbed spotted lanternflies, Lycorma delicatula, (Hemiptera, Fulgoridae). J. Asia Pac. Entomol. 22, 117–120 (2019).Article 

Google Scholar 
28.Lee, J. E. et al. Feeding behavior of Lycorma delicatula (Hemiptera: Fulgoridae) and response on feeding stimulants of some plants. Korean J. Appl. Entomol. 48, 467–477 (2009).Article 

Google Scholar 
29.Liu, H. Seasonal development, cumulative growing degree-days, and population density of spotted lanternfly (Hemiptera: Fulgoridae) on selected hosts and substrates. Environ. Entomol. 49, 1171–1184 (2020).Article 

Google Scholar 
30.Jaenike, J. On optimal oviposition behavior in phytophagous insects. Theor. Popul. Biol. 14, 350–356 (1978).CAS 
Article 

Google Scholar 
31.Gripenberg, S., Mayhew, P. J., Parnell, M. & Roslin, T. A meta-analysis of preference–performance relationships in phytophagous insects. Ecol. Lett. 13, 383–393 (2010).Article 

Google Scholar 
32.Fujiyama, N., Torii, C., Akabane, M. & Katakura, H. Oviposition site selection by herbivorous beetles: a comparison of two thistle feeders: Cassida rubiginosa and Henosepilachnniponica. Entomol. Exp. Appl. 128, 41–48 (2008).Article 

Google Scholar 
33.Wolfin, M. S., Myrick, A. J. & Baker, T. C. Flight duration capabilities of dispersing adult spotted lanternflies, Lycorma delicatula. J. Insect Behav. 33, 125–137 (2020).Article 

Google Scholar 
34.Faraji, F., Janssen, A. & Sabelis, M. W. Oviposition patterns in a predatory mite reduce the risk of egg predation caused by prey. Ecol. Entomol. 27, 660–664 (2002).Article 

Google Scholar 
35.Behmer, S. T. & Joern, A. Coexisting generalist herbivores occupy unique nutritional feeding niches. Proc. Natl. Acad. Sci. USA 105, 1977–1982 (2008).ADS 
CAS 
Article 

Google Scholar 
36.Behmer, S. T. Insect herbivore nutrient regulation. Annu. Rev. Entomol. 54, 165–187 (2009).CAS 
Article 

Google Scholar  More

Acoustic measurementsA set of acoustic echo sounder data was used to analyze fish density. This was collected within the Mycto-3D-MAP program using split-beam echo sounders at 38 and 120 kHz. The Mycto-3D-MAP program included multiple large-scale oceanographic surveys during 2 years and a dedicated cruise in the Kerguelen area. The dataset was collected during 4 large-scale surveys in 2013 and 2014, both in summer (including both northward and southward transects) and in winter, corresponding to 6 acoustic transects of 2860 linear kilometers (see Table 1 for more details). Note that all legs except summer 2014 (MYCTO-3D-Map cruise) were logistic operations, during which environmental in situ data (such as temperature or salinity profiles) could not be collected. The data were then treated with a bi-frequency algorithm, applied to the 38 and 120 kHz frequencies (details of data collection and processing are provided in37). This provides a quantitative estimation of the concentration of gas-bearing organisms, mostly attributed to fish containing a gas-filled swimbladder in the water column38. Most mesopelagic fish present swimbladders and several works pointed out that myctophids are the dominant mesopelagic fish in the region39. Therefore, we considered the acoustic signal as mainly representative of myctophids concentration. Data were organized in acoustic units, averaging acoustic data over 1.1 km along the ship trajectory on average. Myctophid school length is in the order of tens of meters40. For this reason, acoustic units were considered as not autocorrelated. Every acoustic unit is composed of 30 layers, ranging from 10 to 300 meters (30 layers in total).The dataset was used to infer the Acoustic Fish Concentration (AFC) in the water column. We considered as AFC of the point ((x_i), (y_i)) of the ship trajectory the average of the bifrequency acoustic backscattering of 29 out of 30 layers (the first layer, 0-10 m, was excluded due to surface noise) AFC quantity is dimensionless.As the previous measurements were performed through acoustic measurements, a non-invasive technique, fishes were not handled for this study.Table 1 Details of the acoustic transects analyzed.Full size tableRegional data selectionThe geographic area of interest of the present study is the Southern Ocean. To select the ship transects belonging to this region, we used the ecopartition of41. Only points falling in the Antarctic Southern Ocean region were considered. We highlight that this choice is consistent with the ecopartition of42 (group 5), which is specifically designed for myctophids, the reference fish family (Myctophidae) of this study. Furthermore, this choice allowed us to exclude large scale fronts (i.e., fronts that are visible on monthly or yearly averaged maps) which have been the subject of past research works43,44. In this way our analysis is focused specifically on fine-scale fronts.Day-night data separationSeveral species of myctophids present a diel vertical migration. They live at great depths during the day (between 500 and 1000 m), ascending at dusk in the upper euphotic layer, where they spend the night. Since the maximal depth reached by the echo sounder we used is 300 m, in the analysis reported in Figs. 2 and 3 we excluded data collected during the day. However, their analysis is reported in SI.1. A restriction of our acoustic analyses to the ocean upper layer is also consistent with the fact that the fine-scale features we computed are derived in this work by satellite altimetry, thus representative of the upper part of the water column ((sim 50) m). Finally, we note that the choice of considering the echo sounder data in the first 300 m of the water columns is coherent with the fact that LCS may extend almost vertically in depth even at 600 m depth45,46 and with the fact that altimetry-derived velocity fields are consistent with subsurface currents in our region of interest down to 500 m20.Satellite dataVelocity current data and Finite-Size Lyapunov Exponent (FSLE) processing. Velocity currents are obtained from Sea Surface Height (SSH), which is measured by satellite altimetry, through geostrophic approximation. Data, which were downloaded from E.U. Copernicus Marine Environment Monitoring Service (CMEMS, http://marine.copernicus.eu/), were obtained from DUACS (Data Unification and Altimeter Combination System) delayed-time multi-mission altimeter, and displaced over a regular grid with spatial resolution of (frac{1}{4}times frac{1}{4}^circ) and a temporal resolution of 1 day.Trajectories were computed with a Runge-Kutta scheme of the 4th order with an integration time of 6 hours. Finite-size Lyapunov Exponents (FSLE) were computed following the methods in14, with initial and final separation of (0.04^circ) and (0.4^circ) respectively. This Lagrangian diagnostic is commonly used to identify Lagrangian Coherent Structures. This method determines the location of barriers to transport, and it is usually associated with oceanic fronts9. Details on the Lagrangian techniques applied to altimetry data, including a discussion of its limitation, can be found in10.Temperature field and gradient computation The Sea Surface Temperature (SST) field was produced from the OSTIA global foundation Sea Surface Temperature (product id: SST_GLO_ SST_L4_NRT_OBSERVATIONS_010_001) from both infrared and microwave radiometers, and downloaded from CMEMS website. The data are represented over a regular grid with spatial resolution of (0.05times 0.05^circ) and daily-mean maps. The SST gradient was obtained from:$$begin{aligned} Grad SST=sqrt{g_x^2+g_y^2} end{aligned}$$where (g_x) and (g_y) are the gradients along the west-east and the north-south direction, respectively. To compute (g_x), the following expression was used:$$begin{aligned} g_x=frac{1}{2 dx}cdot (SST_{i+1}-SST_{i-1}) end{aligned}$$where the SST values of the adjacent grid cells (along the west-east direction: cells (i+1) and (i-1)) were employed. dx identifies the kilometric distance between two grid points along the longitude (which varies with latitude). The definition is analog for (g_y), considering this time the north-south direction and (dysimeq 5) km (0.05(^circ)).Chlorophyll field Chlorophyll estimations were obtained from the Global Ocean Color product (OCEANCOLOUR_ GLO_CHL_L4_REP_OBSERVATIONS_009_082-TDS) produced by ACRI-ST. This was downloaded from CMEMS website. Daily observations were used, in order to match the temporal resolution of the acoustic and satellite observations. The spatial resolution of the product is 1/24(^{circ }).Estimation of satellite data along ship trajectory For each point ((x_i), (y_i)) of the ship trajectory, we extracted a local value of FSLE, SST gradient, and chlorophyll concentration. These were obtained by considering the respective average value in a circular around of radius (sigma) of each point ((x_i), (y_i)) . Different (sigma) were tested (ranging from 0.1(^circ) to 0.6(^circ)), and the best results were obtained with (sigma =0.2^circ), reference value reported in the present work. This value is consistent with the resolution of the altimetry data.Statistical processingFront identification FSLE and SST gradient were used as diagnostics to detect frontal features, since they are typically associated with front intensity and Lagrangian Coherent Structures10. Note that the two diagnostics provide similar but not identical information. FSLEs are used to analyze the horizontal transport and to identify material lines along which a confluence of waters with different origins occur. These lines typically display an enhanced SST gradient because water masses of different origin have often contrasted SST signature. However, this is not a necessary condition. FSLE ridges may be invisible in SST maps if transport occurs in a region of homogeneous SST. Conversely, SST gradient unveils structures separating waters of different temperatures, whose contrast is often – but not always – associated with horizontal transport. Therefore, even if they usually detect the same structures, these two metrics are complementary. Frontal features were identified by considering a local FSLE (or SST gradient, respectively) value larger than a given threshold. The threshold values have been chosen heuristically but consistently with the ones found in previous works. For the FSLEs, we used 0.08 days(^{-1}), a threshold value in the range of the ones chosen in18 and47. For the SST gradient, we considered representative of thermal front values greater than 0.009({^circ })C/km, which is about half the value chosen in47. However, in that work, the SST gradient was obtained from the advection of the SST field with satellite-derived currents for the previous 3 days, a procedure which is known to enhance tracer gradients.Bootstrap method The threshold value splits the AFC into two groups: AFC co-located with FSLE or SST gradient values over the threshold are considered as measured in proximity of a front (i.e., statistically associated with a front), while AFC values below the threshold are considered as not associated with a frontal structure. The statistical independence of the two groups was tested using a Mann-Whitney or U test. Finally, bootstrap analysis is applied following the methodologies used in47. This allows estimating the probability that the difference in the mean AFC values, over and under the threshold, is significant, and not the result of statistical fluctuations. Other diagnostics tested are reported in SI.1.Linear quantile regression Linear quantile regression method48 was employed as no significant correlation was found between the explanatory and the response variables. This can be due to the fact that multiple factors (such as prey or predator distributions) influence the fish distribution other than the frontal activity considered. The presence of these other factors can shadow the relationship of the explanatory variables (in this case, the FSLE and the SST gradient) with the mean value of the response variable (the AFC). A common method to address this problem is the use of the quantile regression48,49, that explores the influence of the explanatory variables on other parts of the response variable distribution. Previous studies, adopting this methodology, revealed the limiting role played by the explanatory variables in the processes considered50. The percentiles values used are 75th, 90th, 95th, and 99th. The analysis is performed using the statistical package QUANTREG in R (https://CRAN.R-project.org/package=quantreg, v.5.3848,51), using the default settings.Chlorophyll-rich waters selection The AFC observations were considered in chlorophyll-rich waters if the corresponding chlorophyll concentration was higher than a given threshold. All the other AFC measurements were excluded, and a linear regression performed between the remaining AFC and FSLE (or SST gradient) values. The corresponding thresholds (one for FSLE and one for SST gradient case) were selected though a sensitivity test reported in SI.1. These resulted in 0.22 and 0.17 mg/m(^3) for FSLE and for SST gradient, respectively. These values are consistent among them and, in addition, they are in coherence with previous estimates of chlorophyll concentration used to characterise productive waters in the Southern Ocean (0.26mg/m(^3)52).Gradient climbing modelAn individual-based mechanistic model is built to describe how fish could move along frontal features. We assume that the direction of fish movement along a frontal gradient is influenced by the noise of the prey field (SI. 2). Specifically, we consider a Markovian process along the (one dimensional) cross-front direction. Time is discretized in timesteps of length (varDelta tau). We presuppose that, at each timestep, the fish chooses between swimming in one of the two opposite cross-front directions (“left” and “right”). This decision depends on the comparison between the quantity of a tracer (a cue) present at its actual position and the one perceived at a distance (p_R) from it, where (p_R) is the perceptual range of the fish. We consider the fish immersed in a tracer whose concentration is described by the function T(x).An expression for the average velocity of the fish, (U_F(x)), can now be derived. We assume that the fish is able to observe simultaneously the tracer to its right and its left. Fish sensorial capacities are discussed in SI.2. The tracer observed is affected by noise. Noise distribution is considered uniform, with (-xi _{MAX}{tilde{T}}(x_0-varDelta x)), the fish will move to the right, and, vice versa, to the left. We hypothesize that the observational range is small enough to consider the tracer variation as linear, as illustrated in Fig. S7 (SI. 3). In this way:$$begin{aligned}&{tilde{T}}(x_0+varDelta x)=T(x_0)+ p_R,frac{partial T}{partial x}+xi _1 \&{tilde{T}}(x_0-varDelta x)=T(x_0)- p_R,frac{partial T}{partial x}+xi _2 ;. end{aligned}$$In case of absence of noise, or with (xi _{MAX}p_R,frac{partial T}{partial x}). If (T(x_0+varDelta x) >T(x_0-varDelta x)) (as in Fig. S7), and the fish will swim leftward if$$begin{aligned} xi _1-xi _2 >2p_R,frac{partial T}{partial x}; . end{aligned}$$Since (xi _1) and (xi _2) range both between (-xi _{MAX}) and (xi _{MAX}), we can obtain the probability of leftward moving P(L). This will be the probability that the difference between (xi _1) and (xi _2) is greater than (2p_R,frac{partial T}{partial x})$$begin{aligned} P(L)&=frac{1}{8xi _{MAX}^2} bigg (2 xi _{MAX} – 2 p_R,frac{partial T}{partial x}bigg )^2\&=frac{1}{2} bigg (1-frac{p_R}{xi _{MAX}},frac{partial T}{partial x}bigg )^2 end{aligned}$$.The probability of moving right will be$$begin{aligned} P(R)&=1-P(L) end{aligned}$$and their difference gives the frequency of rightward moving$$begin{aligned} P(R)-P(L)&=1-2P(L)=1-bigg (1-frac{p_R}{xi _{MAX}},frac{partial T}{partial x}bigg )^2\&=frac{p_R}{xi _{MAX}}frac{partial T}{partial x}bigg (2-frac{p_R}{xi _{MAX}}bigg |frac{partial T}{partial x}bigg |bigg ); , end{aligned}$$where the absolute value of (frac{partial T}{partial x}) has been added to preserve the correct climbing direction in case of negative gradient. The above expression leads to:$$begin{aligned} U_F(x)=frac{V p_R}{xi _{MAX}}frac{partial T}{partial x}bigg (2-frac{p_R}{xi _{MAX}}bigg |frac{partial T}{partial x}bigg |bigg );. end{aligned}$$
(1)
We then assume that, over a certain value of tracer gradient (frac{partial T}{partial x}_{MAX}), the fish are able to climb the gradient without being affected by the noise. This assumption, from a biological perspective, means that the fish does not live in a completely noisy environment, but that under favorable circumstances it is able to correctly identify the swimming direction that leads to higher prey availability. This means that$$begin{aligned} p_R*frac{partial T}{partial x}_{MAX}=xi _{MAX},. end{aligned}$$
(2)
Substituting then (2) into (1) gives:$$begin{aligned} U_F(x)=V frac{frac{partial T}{partial x}}{frac{partial T}{partial x}_{MAX}}bigg (2-frac{big |frac{partial T}{partial x}big |}{frac{partial T}{partial x}_{MAX}}bigg );. end{aligned}$$
(3)
Finally, we can include an eventual effect of transport by the ocean currents, considering that the tracer is advected passively by them, simply adding the current speed (U_C) to Expr. (3).The representations provided are individual based, with each individual representing a single fish. However, we note that all the considerations done are also valid if we consider an individual representing a fish school. (U_F) will then simply represent the average velocity of the fish schools. This aspect should be stressed since many fish species live and feed in groups, especially myctophids (see SI.2 for further details).Continuity equation in one dimension The solution of this model will now be discussed. The continuity equation is first considered in one dimension. The starting scenario is simply an initially homogeneous distribution of fish, that are moving in a one dimensional space with a velocity given by (U_{F}(x)).We assume that in the time scales considered (few days to some weeks), the fish biomass is conserved, so for instance fishing mortality or growing rates are neglected. In that case, we can express the evolution of the concentration of the fish (rho) by the continuity equation$$begin{aligned} frac{partial rho }{partial t}+nabla cdot mathbf{j },=,0 end{aligned}$$
(4)
in which (mathbf{j }=rho ;U_{F}(x)), so that Eq. (4) becomes$$begin{aligned} frac{partial rho }{partial t}+nabla cdot big (rho ;U_{F}(x)big ),=,0;. end{aligned}$$
(5)
In one dimension, the divergence is simply the partial derivate along the x-axis. Eq. (5) becomes$$begin{aligned} frac{partial rho }{partial t}=-frac{partial }{partial x} bigg (rho ;U_{F}bigg ) end{aligned}$$
(6)
Now, we decompose the fish concentration (rho) in two parts, a constant one and a variable one (rho ,=,rho _0+{tilde{rho }}). Eq. (6) will then become$$begin{aligned} frac{partial rho }{partial t}=-U_Ffrac{partial {tilde{rho }}}{partial x}-rho frac{partial U_F}{partial x};. end{aligned}$$
(7)
Using Expr. (3), Eq. (7) is numerically simulated with the Lax method. In Expr. (3) we impose that (U_F(x)=V) when (U_F(x) >V). This biological assumption means that fish maximal velocity is limited by a physiological constraint rather than by gradient steepness. Details of the physical and biological parameters are provided in SI.6. More

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