Philippa Kaur

Seguin, R. et al. Nat. Sustain. https://doi.org/10.1038/s41893-022-00981-x (2022).Article 

Google Scholar 
Roberts, C. M. & Polunin, N. V. C. Rev. Fish Biol. Fish. 1, 65–91 (1991).Article 

Google Scholar 
Cinner, J. E. et al. Soc. Nat. Resour. 27, 994–1005 (2014).Article 

Google Scholar 
MacNeil, M. A. et al. Nature 520, 341–344 (2015).Article 
CAS 

Google Scholar 
Morais, R. A. & Bellwood, D. R. Coral Reefs 39, 1221–1231 (2020).Article 

Google Scholar 
Morais, R. A., Connolly, S. R. & Bellwood, D. R. Glob. Change Biol. 26, 1295–1305 (2020).Article 

Google Scholar 
Di Lorenzo, M. et al. Fish Fish. 21, 906–915 (2020).Article 

Google Scholar 
Ban, N. C. et al. Nat. Sustain. 2, 524–532 (2019).Article 

Google Scholar 
Rogers, A. et al. Ecology 99, 450–463 (2018).Article 

Google Scholar 
Robinson, J. P. W. et al. Nat. Ecol. Evol. 3, 183–190 (2019).Article 

Google Scholar  More

Dataset acquisitionPlant-frugivore network data were obtained through different online sources and publications (Supplementary Table 1). Only networks that met the following criteria were retrieved: (i) the network contains quantitative data (a measure of interaction frequency) from a location, pooling through time if necessary; (ii) the network includes avian frugivores. Importantly, we removed non-avian frugivores from our analyses because only 28 out of 196 raw networks (before data cleaning) sampled non-avian frugivores, and not removing non-avian frugivores would generate spurious apparent turnover between networks that did vs. did not sample those taxa. In addition, the removal of non-avian frugivores did not strongly decrease the number of frugivores in our dataset (Supplementary Fig. 20a) or the total number of links in the global network of frugivory (Supplementary Fig. 20b). Furthermore, non-avian frugivores, as well as their interactions, were not shared across ecoregions and biomes (Supplementary Fig. 21), so their inclusion would only strengthen the results we found (though as noted above, we believe that this would be spurious because they are not as well sampled); (iii) the network (after removal of non-avian frugivores) contains greater than two species in each trophic level. Because this size threshold was somewhat arbitrary, we used a sensitivity analysis to assess the effect of our network size threshold on the reported patterns (see Sensitivity analysis section in the Supplementary Methods and Supplementary Figs. 22–24); and (iv) network sampling was not taxonomically restricted, that is, sampling was not focused on a specific taxonomic group, such as a given plant or bird family. Note, however, that authors often select focal plants or frugivorous birds to be sampled, but this was not considered as a taxonomic restriction if plants and birds were not selected based on their taxonomy (e.g., focal plants were selected based on the availability of fruits at the time of sampling, or focal birds were selected based on previous studies of bird diet in the study site). The first source for network data was the Web of Life database42, which contains 33 georeferenced plant-frugivore networks from 28 published studies, of which 12 networks met our criteria.We also accessed the Scopus database on 04 May 2020 using the following keyword combination: (“plant-frugivore*” OR “plant-bird*” OR “frugivorous bird*” OR “avian frugivore*” OR “seed dispers*”) AND (“network*” OR “web*”) to search for papers that include data on avian frugivory networks. The search returned a total of 532 studies, from which 62 networks that met the above criteria were retrieved. We also contacted authors to obtain plant-frugivore networks that were not publicly available, which provided us a further 110 networks. The remaining networks (N = 12) were obtained by checking the database from a recently published study12. In total, 196 quantitative avian frugivory networks were used in our analyses.Generating the distance matrices to serve as predictor and response variablesEcoregion and biome distancesWe used the most up-to-date (2017) map of ecoregions and biomes3, which divides the globe into 846 terrestrial ecoregions nested within 14 biomes, to generate our ecoregion and biome distance matrices. Of these, 67 ecoregions and 11 biomes are represented in our dataset (Supplementary Figs. 1 and 2). We constructed two alternative versions of both the ecoregion and biome distance matrices. In the first, binary version, if two ecological networks were from localities within the same ecoregion/biome, a dissimilarity of zero was given to this pair of networks, whereas a dissimilarity of one was given to a pair of networks from distinct ecoregions/biomes (this is the same as calculating the Euclidean distance on a presence–absence matrix with networks in rows and ecoregion/biomes in columns).In the second, quantitative version, we estimated the pairwise environmental dissimilarity between our ecoregions and biomes using six environmental variables recently demonstrated to be relevant in predicting ecoregion distinctness, namely mean annual temperature, temperature seasonality, mean annual rainfall, rainfall seasonality, slope and human footprint38. We obtained climatic and elevation data from WorldClim 2.143 at a spatial resolution of 1-km2. We transformed the elevation raster into a slope raster using the terrain function from the raster package44 in R45. As a measure of human disturbance, we used human footprint—a metric that combines eight variables associated with human disturbances of the environment: the extent of built environments, crop land, pasture land, human population density, night-time lights, railways, roads and navigable waterways26. The human footprint raster was downloaded at a 1-km2 resolution26. Because human footprint data were not available for one of our ecoregions (Galápagos Islands xeric scrub), we estimated human footprint for this ecoregion by converting visually interpreted scores into the human footprint index. We did this by analyzing satellite images of the region and following a visual score criterion26. Given the previously demonstrated strong agreement between visual score and human footprint values26, we fitted a linear model using the visual score and human footprint data from 676 validation plots located within the Deserts and xeric shrublands biome – the biome in which the Galápagos Islands xeric scrub ecoregion is located – and estimated the human footprint values for our own visual scores using the predict function in R45.We used 1-km2 resolution rasters and the extract function from the raster package44 to calculate the mean value of each of our six environmental variables for each ecoregion in our dataset. Because biomes are considerably larger than ecoregions (which makes obtaining environmental data for biomes more computationally expensive) we used a coarser spatial resolution of 5-km2 for calculating the mean values of environmental variables for each biome. Since a 5-km2 resolution raster was not available for human footprint, we transformed the 1-km2 resolution raster into a 5-km2 raster using the resample function from the same package.To combine these six environmental variables into quantitative matrices of ecoregion and biome environmental dissimilarity, we ran a Principal Component Analysis (PCA) on our scaled multivariate data matrix (where rows are ecoregions or biomes and columns are environmental variables). From this PCA, we selected the scores of the four and three principal components, which represented 89.6% and 88.7% of the variance for ecoregions and biomes, respectively, and converted it into a distance matrix by calculating the Euclidean distance between pairs of ecoregions/biomes using the vegdist function from the vegan package46. Finally, we transformed the ecoregion or biome distance matrix into a N × N matrix where N is the number of local networks. In this matrix, cell values represent the pairwise environmental dissimilarity between the ecoregions/biomes where the networks are located. The main advantage of using this quantitative approach is that, instead of simply evaluating whether avian frugivory networks located in distinct ecoregions or biomes are different from each other in terms of network composition and structure (as in our binary approach), we were also able to determine whether the extent of network dissimilarity depended on how environmentally different the ecoregions or biomes are from one another.Local-scale human disturbance distanceTo generate our local human disturbance distance matrix, we extracted human footprint data at a 1-km2 spatial resolution26 and calculated the mean human footprint values within a 5-km buffer zone around each network site. For the networks located within the Galápagos Islands xeric scrub ecoregion (N = 4), we estimated the human footprint index using the same method described in the previous section for ecoregion- or biome-scale human footprint. We then calculated the pairwise Euclidean distance between human footprint values from our network sites. Thus, low cell values in the local human disturbance distance matrix indicate pairs of network sites with a similar level of human disturbance, while high values represent pairs of network sites with very different levels of human disturbance.Spatial distanceThe spatial distance matrix was generated using the Haversine (i.e., great circle) distance between all pairwise combinations of network coordinates. In this matrix, cell values represent the geographical distance between network sites.Elevational differenceWe calculated the Euclidean distance between pairwise elevation values (estimated as meters above sea level) of network sites to generate our elevational difference matrix. Elevation values were obtained from the original sources when available or using Google Earth47. In the elevational difference matrix, low cell values represent pairs of network sites within similar elevations, whereas high values represent pairs of network sites within very different elevations.Network sampling dissimilarityWe used the metadata retrieved from each of our 196 local networks to generate our network sampling dissimilarity matrices, which aim to control statistically for differences in network sampling. There are many ways in which sampling effort could be quantified, so we began by calculating a variety of metrics, then narrowed our options by assessing which of these was most related to network metrics. We divided the sampling metrics into two categories: time span-related metrics (i.e., sampling hours and months) and empirical metrics of sampling completeness (i.e., sampling completeness and sampling intensity), which aim to account for how complete network sampling was in terms of species interactions (Supplementary Table 2).We selected the quantitative sampling metrics to be included in our models based on (i) the fit of generalized linear models evaluating the relationship between number of sampling hours and sampling months of the study and network-level metrics (i.e., bird richness, plant richness and number of links), and (ii) how well time span-related metrics, sampling completeness and sampling intensity predicted the proportion of known interactions that were sampled in each local network (hereafter, ratio of interactions) for a subset of the data. This latter metric, defined as the ratio between the number of interactions in the local network and the number of known possible interactions in the region involving the species in the local network, captures raw sampling completeness. Therefore, ratio of interactions estimates, for a given set of species, the proportion of all their interactions known for a region that are found to occur among those same species in the local network. To calculate this metric, we needed high-resolution information on the possible interactions, so we used a subset of 14 networks sampled in Aotearoa New Zealand, since there is an extensive compilation of frugivory events recorded for this country48. After this process, we selected number of sampling hours, number of sampling months and sampling intensity for inclusion in our statistical models (Supplementary Figs. 7 and 8; Supplementary Table 2). We generated the corresponding distance matrices by calculating the Euclidean distance between metric values. Similarly, we generated a Euclidean distance matrix for differences in sampling year between pairs of networks, which aims to account for long-term changes in the environment, species composition and network sampling methods. We obtained the sampling year of our local networks from the original sources and calculated the mean sampling year value for those networks sampled across multiple years.Because sampling methods, such as sampling design, focus (i.e., focal taxa, which determines whether a zoocentric or phytocentric method was used), interaction frequency type (i.e., how interaction frequency was measured) and coverage (total or partial) might also affect the observed plant-frugivore interactions49, we combined these variables into a single distance matrix to estimate the overall differences in sampling methods between networks. Because most of these variables were categorical with multiple levels (Supplementary Table 3), we generated our method’s dissimilarity matrix by using a generalization of Gower’s distance method50, which allows the treatment of different types of variables when calculating distances. For this, we used the dist.ktab function from the ade4 package51. We ran a Principal Coordinates Analysis (PCoA) on this distance matrix, selected the first four axes, which explained 81.2% of the variation in method’s dissimilarity, and calculated the Euclidean distance between pairs of networks using the vegdist function from the vegan package46 in R45.Network dissimilarityWe generated three network dissimilarity matrices to be our response variables in the statistical models. In the first, cell values represent the pairwise dissimilarity in species composition between networks (beta diversity of species; βS)27. Second, we measured interaction dissimilarity (beta diversity of interactions; βWN), which represents the pairwise dissimilarity in the identity of interactions between networks27. Importantly, we did not include interaction rewiring (βOS) in our main analysis because this metric can only be calculated for networks that share interaction partners (i.e., it estimates whether shared species interact differently)27, which limited the number and the spatial distribution of networks available for analysis (but see the Rewiring analysis section for an analysis on the subset of our dataset for which this was possible). Metrics were calculated using the network_betadiversity function from the betalink package52 in R45.Finally, we calculated a third dissimilarity matrix to capture overall differences in network structure. We recognize that there are many potential metrics of network structure, and that many of these are strongly correlated with one another53,54,55,56. We therefore chose a range of metrics that captured the number of links, their relative weightings (including across trophic levels), and their arrangement among species, then combined these into a single distance matrix. Specifically, we quantified network structural dissimilarity using the following metrics: weighted connectance, weighted nestedness, interaction evenness, PDI and modularity.Weighted connectance represents the number of links relative to the number of possible links, weighted by the frequency of each interaction55, and is therefore a measure of network-level specialization (higher values of weighted connectance indicate lower specialization). Importantly, it has been suggested that connectance affects persistence in mutualistic systems54. We measured nestedness (i.e., the pattern in which specialist species interact with proper subsets of the species that generalist species interact with) using the weighted version of nestedness based on overlap and decreasing fill (wNODF)57. Notably, nested structures have been commonly reported in plant-frugivore networks33. Interaction evenness is Shannon’s evenness index applied for species interactions and represents how evenly distributed the interactions are in the network21,58. This metric has been previously demonstrated to decline with habitat modification as a consequence of some interactions being favored over others in high-disturbance environments21. PDI (Paired Difference Index) is a measure of species-level specialization on resources and a reliable indicator not only of specialization, but also of absolute generalism59. Thus, this metric contributes to understanding of the ecological processes that drive the prevalence of specialists or generalists in ecological networks59. In order to obtain a network-level PDI, we calculated the weighted mean PDI for each local network. Finally, we calculated modularity (i.e., the level of compartmentalization within networks) using the DIRTPLAwb+ algorithm60. Modularity estimates the extent to which species within modules interact more with each other than with species from other modules61, and it has been demonstrated to affect the persistence and resilience of mutualistic networks54. All the selected network metrics are based on weighted (quantitative) interaction data, as these have been suggested to be less biased by sampling incompleteness62 and to better reflect environmental changes21. All network metrics were calculated using the bipartite package63 in R45.We ran a Principal Component Analysis (PCA) on our scaled multivariate data matrix (N × M where N is the number of local networks in our dataset and M is the number of network metrics), selected the scores of the three principal components, which represented 89.9% of the variance in network metrics, and converted it into a network structural dissimilarity matrix by calculating the Euclidean distance between networks. In this distance matrix, cell values represent differences in the overall architecture of networks (over all the network metrics calculated), and therefore provide a complementary approach for evaluating how species interaction patterns vary across large-scale environmental gradients.Statistical analysisWe employed a two-tailed statistical test that combines Generalized Additive Models (GAM)29 and Multiple Regression on distance Matrices (MRM)30 to evaluate the effect of each of our predictor distance matrices on our response matrix. With this approach, we were able to fit GAMs where the predictor and responsible variables are distance matrices, while accounting for the non-independence of distances from each local network by permuting the response matrix30. The main advantage of using GAMs is their flexibility in modeling non-linear relationships through smooth functions, which are represented by a sum of simpler, fixed basis functions that determine their complexity29. Using GAM-based MRM models allowed us to obtain F values for each of the smooth terms (i.e., smooth functions of the predictor variables in our model), and test statistical significance at the level of individual variables. The binary versions of ecoregion and biome distance matrices (with two levels, “same” or “distinct”) were treated as categorical variables in the models, and t values were used for determining statistical significance. We fitted GAMs with thin plate regression splines64 using the gam function from the mgcv package29 in R45. Smoothing parameters were estimated using restricted maximum likelihood (REML)29. Our GAM-based MRM models were calculated using a modified version of the MRM function from the ecodist package65, which allowed us to combine GAMs with the permutation approach from the original MRM function (see Code availability). All the models were performed with 1000 permutations (i.e., shuffling) of the response matrix.We explored the unique and shared contributions of our predictor variables to network dissimilarity using deviance partitioning analyses. These were performed by fitting reduced models (i.e., GAMs where one or more predictor variables of interest were removed) using the same smoothing parameters as in the full model and comparing the explained deviance. We fixed smoothing parameters for comparisons in this way because these parameters tend to vary substantially (to compensate) if one of two correlated predictors is dropped from a GAM.Assessing the influence of individual studies on the reported patternsBecause our dataset comprises 196 local frugivory networks obtained from 93 different studies, and some of these studies contained multiple networks, we needed to evaluate whether our results were strongly biased by individual studies. To do this, we followed the approach from a previous study66 and tested whether F values of smooth terms and t values of categorical variables (binary version of ecoregion and biome distances) changed significantly when jackknifing across studies. We did this by dropping one study from the dataset and re-fitting the models, and then repeating this same process for all the studies in our dataset.We found a number of consistent patterns within different subsets of the data (Supplementary Figs. 15 and 16); however, some of the patterns we observed appear to be driven by individual studies with multiple networks, and hence are less representative. For instance, the study with the greatest number of networks in our dataset (study ID = 76), which contains 35 plant-frugivore networks sampled across an elevation gradient in Mt. Kilimanjaro, Tanzania67, had an overall high influence on the results when compared with the other studies. By re-running our GAM-based MRM models after removing this study from our dataset, we found that the effect of biome boundaries on interaction dissimilarity is no longer significant, whereas the effects of ecoregion boundaries, human disturbance distance, spatial distance and elevational differences remained consistent with those from the full dataset (Supplementary Table 33). Nevertheless, all the results were qualitatively similar to those obtained for the entire dataset when using network structural dissimilarity as the response variable (Supplementary Table 34).Rewiring analysisInteraction rewiring (βOS) estimates the extent to which shared species interact differently27. Because this metric can only be calculated for networks that share species from both trophic levels, we selected a subset of network pairs that shared plants and frugivorous birds (N = 1314) to test whether interaction rewiring increases across large-scale environmental gradients. Importantly, since not all possible combinations of network pairs contained values of interaction rewiring (i.e., not all pairs of networks shared species), a pairwise distance matrix could not be generated for this metric. Thus, we were not able to use the same statistical approach used in our main analysis, which is based on distance matrices (see Statistical analysis section). Instead, we performed a Generalized Additive Mixed-effects Model (GAMM) using ecoregion, biome, human disturbance, spatial, elevational, and sampling-related distance metrics as fixed effects and network IDs as random effects (to account for the non-independence of distances) (Supplementary Table 35). We also performed a reduced model with only ecoregion and biome distance metrics as predictor variables (Supplementary Table 36). The binary version of ecoregion and biome distance metrics (with two levels, “same” or “distinct”) were used as categorical variables in both models. Interaction rewiring (βOS) was calculated using the network_betadiversity function from the betalink package52 in R45. Although it has been recently argued that this metric may overestimate the importance of rewiring for network dissimilarity68, our main focus was not the partitioning of network dissimilarity into species turnover and rewiring components, but rather simply detecting whether the sub-web of shared species interacted differently. In this case, βOS (as developed by ref. 27) is an adequate and useful metric68. We fitted our models using the gamm4 function from the gamm4 package69 in R45. Smoothing parameters were estimated using restricted maximum likelihood (REML)29.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

All methods were performed in accordance with the local relevant guidelines, regulations and legislation.InstrumentsLI-6400 photosynthesis system (LI-6400 Inc., Lincoln, NE, USA) and PAM-2500 portable chlorophyll fluorescence apparatus (PAM-2500, Walz, Germany) were used in the study.MaterialsAbout 90 living E. brevicornu plants were collected from Taihang Mountains in October 2018. The E. brevicornu was not in endangered or protected. The collection of these E. brevicornu plants was permitted by local government. These plants were averagely planted in nine plots of 2 m2. The roots of E. pubescens were planted 6–8 cm below ground. These plots were placed on farmland near Taihang Mountains and covered with sunshade net (about 70% light transmittance). These plants were timely irrigated after planting to ensure that they grew well but not fertilized.Determination of photosynthetic characteristicsThe photosynthetic characteristics of mature leaves on the E. brevicornu plants were determined between June 6–8, 2019 with the Li-6400 photosynthesis system. The diurnal variation of photosynthesis in three leaves of three plants was determined. When the light response curve was determined, the temperature of the leaf chamber was set at 28 °C, and the concentration of CO2 in the leaf chamber was set at 400 µmol mol−1. When determining the CO2 response curve, the light intensity in the leaf chamber was set at 1000 µmol m−2 s−1, and the temperature of the leaf chamber was set at 28 °C. The light response curve and CO2 response curve were determined three times in three leaves of three different plants.Determination of chlorophyll fluorescence characteristicsThe fluorescence characteristics of chlorophyll in E. brevicornu leaves were determined with PAM-2500 portable chlorophyll fluorescence apparatus between June 8–9, 2019. The leaves underwent dark adaptation for 30 min before determining slow kinetics of chlorophyll fluorescence. Then the light curves of chlorophyll fluorescence were determined. All of these determinations were repeated three times on three mature leaves of three plants.The data was analysed with SPSS (Statistical Product and Service Solutions, International Business Machines Corporation, USA). The light response curves were fitted with following modified rectangular hyperbola model11,12.$${text{Photo}}, = ,{text{E}}cdotleft( {{1} – {text{M}}cdot{text{PAR}}} right)cdotleft( {{text{PAR}} – {text{LCP}}} right)/({1}, + ,{text{N}}cdot{text{PAR}})$$PAR is the value of light intensity in leaf chamber of Li-6400 photosynthesis system. Photo is net photosynthetic rate. LCP is the light compensation point. E is the apparent quantum yield. M and N are parameters. The dark respiration rate under the LCP is calculated according to E·LCP. The light saturation point (LSP) is calculated according to (((M + N) ·(1 + N·LCP)/M)½)/−1)/N.The net photosynthetic rate under the light saturation point (LSP) can be calculated according to the above model.The CO2 response curves were fitted with below modified rectangular hyperbola model11,12.$${text{Photo}}, = ,{text{E}}cdotleft( {{1} – {text{M}}cdot{text{PAR}}} right)cdotleft( {{text{PAR}} – {text{CCP}}} right)/({1}, + ,{text{N}}cdot{text{PAR}})$$PAR is the value of light intensity in leaf chamber of Li-6400 photosynthesis system. Photo is net photosynthetic rate. CCP is CO2 compensation point. E is also the apparent quantum yield. M and N are parameters. The dark respiration rate under the CO2 calculated according to E·CCP. The CO2 saturation point (CSP) is calculated according to (((M + N) ·(1 + N·CCP)/M)½)/−1)/N.The net photosynthetic rate under the CO2 saturation point (CSP) can be alternatively calculated according to the above model.The light curves of chlorophyll fluorescence were fitted according to the below model of Eilers and Peeters12,13.$${text{ETR}}, = ,{text{PAR}}/({text{a}}cdot{text{PAR}}^{{2}} , + ,{text{b}}cdot{text{PAR}}, + ,{text{c}})$$ETR is the electron transport rate of photosynthetic system II. PAR is fluorescence intensity. The letters a, b and c are parameters. More

Balvanera, P. et al. Bioscience 64, 49–57 (2014).Article 

Google Scholar 
Hooper, D. U. et al. Ecol. Monogr. 75, 3–35 (2005).Article 

Google Scholar 
Balvanera, P. et al. Ecol. Lett. 9, 1146–1156 (2006).Article 
PubMed 

Google Scholar 
Cardinale, B. J. et al. Am. J. Bot. 98, 572–592 (2011).Article 
PubMed 

Google Scholar 
Cardinale, B. J. B. J. et al. Nature 486, 59–67 (2012).Article 
CAS 
PubMed 

Google Scholar 
Manning, P. et al. in Advances in Ecological Research (eds Eisenhauer N. et al.) 323–356 (Academic, 2019).Le Provost, G. et al. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-022-01918-5 (2022).Felipe-Lucia, M. R. et al. Proc. Natl Acad. Sci. USA 117, 28140–28149 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Foley, J. A. et al. Science 309, 570–574 (2005).Article 
CAS 
PubMed 

Google Scholar 
Cardinale, B. J. et al. Ecology 94, 1697–1707 (2013).Article 
PubMed 

Google Scholar 
Teles da Mota, V. & Pickering, C. J. Outdoor Recreat. Tour. 30, 100295 (2020).Article 

Google Scholar 
Mitchell, M. G. E. et al. Trends Ecol. Evol. 30, 190–198 (2015).Article 
PubMed 

Google Scholar 
Raudsepp-Hearne, C. & Peterson, G. D. Ecol. Soc. 21, 16 (2016).Article 

Google Scholar 
Chaplin-Kramer, R. & Kremen, C. Ecol. Appl. 22, 1936–1948 (2012).Article 
PubMed 

Google Scholar  More

The plots of Igeo, PERI, and PLI of HMs in the topsoil of the tourist area of Sayram Lake (Fig. 5) reveal the degree of HM pollution and eco-risk in this study area on the one hand and, on the other hand, indicate the direction for the relevant agencies to target soil environmental protection and HM pollution prevention and control measures. In this study, the Igeo results showed that Cd was the most highly enriched HM, and Pb, Zn, Cd, and Ni were slightly enriched in a few sample sites. The unnatural accumulation of these elements is usually closely associated with human activities in the area34. Tourism is the main economic activity in the district, and published studies have reported that tourism infrastructure construction (e.g., roads, buildings, etc.) and tourism wastes (e.g., plastic bags, batteries, hotel wastewater) release Cd into the soil35. Additionally, the accumulation of Pb, Zn, Cu and Ni in soils is usually associated with traffic emissions36. The PERI showed that the study area was at low risk overall, with only point ss04 exhibiting medium risk; however, this result was caused by the abnormally high Cd concentration value (Fig. 4) at point ss04 (Cd (concentration): 1.08 mg/kg, Cd (background): 0.34 mg/kg). This anomalous concentration value has a large influence on the PERI calculated based on the measured concentration, the background value and the toxicity coefficient. Therefore, references to this point can be appropriately removed when considering eco-risk. The PLI of each sampling point was greater than 1 and less than 2, which means that the area was in a moderately contaminated state. In general, the degree of soil HM contamination in this area was low; however, due to HM toxicity, bioaccumulation, and persistence37, the HM contamination of this area still requires sustained attention.Figure 5Contamination and ecological risk indices: (a) geoaccumulation index (Igeo) of HMs; (b) ecological risk of individual HMs; (c) potential ecological risk index (PERI) of HMs; (d) pollution load index (PLI) of HMs.Full size imageCorrelation analysis is an efficient way to reveal correlations among HMs through Pearson correlation coefficients, and HMs with significant correlations may originate from the same source38. As shown in Table S5, the elemental pairs Cd-Cu (p  More

Chauhan, B. S. Grand challenges in weed management. Front. Agron. https://doi.org/10.3389/fagro.2019.00003 (2020).Article 

Google Scholar 
Oerke, E. C. Crop losses to pests. J. Agric. Sci. 144, 31–43 (2006).
Google Scholar 
Samejima, H. & Sugimoto, Y. Recent research progress in combatting root parasitic weeds. Biotechnol. Biotechnol. Equip. 32(2), 221–240 (2018).CAS 

Google Scholar 
Aly, R. Conventional and biotechnological approaches for control of parasitic weeds. In Vitro Cell. Dev. Biol. Plant 43(4), 304–317 (2007).
Google Scholar 
Fernández-Aparicio, M., Delavault, P. & Timko, M. P. Management of infection by parasitic weeds: A review. Plants 9(9), 1184 (2020).PubMed Central 

Google Scholar 
Rodenburg, J., Demont, M., Zwart, S. J. & Bastiaans, L. Parasitic weed incidence and related economic losses in rice in Africa. Agric. Ecosyst. Environ. 235, 306–317 (2016).
Google Scholar 
Weisberger, D., Nichols, V. & Liebman, M. Does diversifying crop rotations suppress weeds? A meta-analysis. PLoS One 14(7), e0219847 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ejeta, G. The Striga scourge in Africa: A growing pandemic. In Integrating New Technologies for Striga Control: Towards Ending the Witch-hunt 3–16 (World Scientific, 2007). https://doi.org/10.1142/9789812771506_0001.Chapter 

Google Scholar 
Netting, R. M. & Stone, M. P. Agro-diversity on a farming frontier: Kofyar smallholders on the Benue plains of central Nigeria. Africa 66(1), 52–70 (1996).
Google Scholar 
Pimentel, D. et al. Conserving biological diversity in agricultural and forestry systems. Bioscience 42, 354–362 (1992).
Google Scholar 
Khoshbakht, K. & Hammer, K. How many plant species are cultivated?. Genet. Resour. Crop Evol. 55(7), 925–928. https://doi.org/10.1007/s10722-008 (2008).Article 

Google Scholar 
Hajjar, R., Jarvis, D. I. & Gemmill-Herren, B. The utility of crop genetic diversity in maintaining ecosystem services. Agric. Ecosyst. Environ. 123(4), 261–270 (2008).
Google Scholar 
He, H. M. et al. Crop diversity and pest management in sustainable agriculture. J. Integr. Agric. 18(9), 1945–1952 (2019).
Google Scholar 
Ofori, F. & Stern, W. R. Cereal–legume intercropping systems. Adv. Agron. 41, 41–90 (1987).
Google Scholar 
Tanveer, M., Anjum, S. A., Hussain, S., Cerdà, A. & Ashraf, U. Relay cropping as a sustainable approach: Problems and opportunities for sustainable crop production. Environ. Sci. Pollut. Res. 24(8), 6973–6988 (2017).
Google Scholar 
Hartwig, N. L. & Ammon, H. U. Cover crops and living mulches. Weed Sci. 50(6), 688–699 (2002).CAS 

Google Scholar 
Raseduzzaman, M. D. & Jensen, E. S. Does intercropping enhance yield stability in arable crop production? A meta-analysis. Eur. J. Agron. 91, 25–33 (2017).
Google Scholar 
Davis, A. S., Hill, J. D., Chase, C. A., Johanns, A. M. & Liebman, M. Increasing cropping system diversity balances productivity, profitability and environmental health. PLoS One 7(10), e47149 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Himmelstein, J., Ares, A., Gallagher, D. & Myers, J. A meta-analysis of intercropping in Africa: Impacts on crop yield, farmer income, and integrated pest management effects. Int. J. Agric. Sustain. 15(1), 1–10 (2017).
Google Scholar 
Abson, D. J., Fraser, E. D. & Benton, T. G. Landscape diversity and the resilience of agricultural returns: A portfolio analysis of land-use patterns and economic returns from lowland agriculture. Agric. Food Secur. 2(1), 1–15 (2013).
Google Scholar 
Renard, D. & Tilman, D. National food production stabilized by crop diversity. Nature 571(7764), 257–260 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Gaudin, A. C. et al. Increasing crop diversity mitigates weather variations and improves yield stability. PLoS One 10(2), e0113261 (2015).PubMed 
PubMed Central 

Google Scholar 
Bowles, T. M. et al. Long-term evidence shows that crop-rotation diversification increases agricultural resilience to adverse growing conditions in North America. One Earth 2(3), 284–293 (2020).ADS 

Google Scholar 
Chauhan, B. S., Singh, R. G. & Mahajan, G. Ecology and management of weeds under conservation agriculture: A review. Crop Prot. 38, 57–65 (2012).
Google Scholar 
Nichols, V., Verhulst, N., Cox, R. & Govaerts, B. Weed dynamics and conservation agriculture principles: A review. Field Crop Res. 183, 56–68 (2015).
Google Scholar 
Banik, P., Midya, A., Sarkar, B. K. & Ghose, S. S. Wheat and chickpea intercropping systems in an additive series experiment: Advantages and weed smothering. Eur. J. Agron. 24(4), 325–332 (2006).
Google Scholar 
Workayehu, T. & Wortmann, C. S. Maize–bean intercrop weed suppression and profitability in Southern Ethiopia. Agron. J. 103(4), 1058–1063 (2011).
Google Scholar 
Haugaard-Nielsen, H., Ambus, P. & Jensen, E. S. Interspecific competition, N use and interference with weeds in pea barley intercropping. Field Crop Res. 70, 101–109 (2001).
Google Scholar 
Jensen, E. S. Intercropping of Cereals and Grain Legumes for Increased Production, Weed Control, Improved Product Quality and Prevention of N-losses in European Organic Farming Systems, Final Report on Intercrop Project (QLK5-CT-2002-02352) (Risø National Laboratory, 2006).Arlauskienė, A., Šarūnaitė, L., Kadžiulienė, Ž, Deveikytė, I. & Maikštėnienė, S. Suppression of annual weeds in pea and cereal intercrops. Agron. J. 106(5), 1765–1774 (2014).
Google Scholar 
Szumigalski, A. & van Acker, R. Weed suppression and crop production in annual intercrops. Weed Sci. 53(6), 813–825 (2005).CAS 

Google Scholar 
Stoltz, E. & Nadeau, E. Effects of intercropping on yield, weed incidence, forage quality and soil residual N in organically grown forage maize (Zea mays L.) and faba bean (Vicia faba L.). Field Crop Res. 169, 21–29 (2014).
Google Scholar 
Sauerborn, J., Müller-Stöver, D. & Hershenhorn, J. The role of biological control in managing parasitic weeds. Crop Prot. 26(3), 246–254 (2007).
Google Scholar 
Jamil, M., Rodenburg, J., Charnikhova, T. & Bouwmeester, H. J. Pre-attachment Striga hermonthica resistance of New Rice for Africa (NERICA) cultivars based on low strigolactone production. New Phytol. 192(4), 964–975. https://doi.org/10.1111/j.1469-8137.2011.03850.x (2011).Article 
CAS 
PubMed 

Google Scholar 
Yoneyama, K. et al. Nitrogen deficiency as well as phosphorus deficiency in sorghum promotes the production and exudation of 5-deoxystrigol, the host recognition signal for arbuscular mycorrhizal fungi and root parasites. Planta 227(1), 125–132. https://doi.org/10.1007/s00425-007-0600-5 (2007).Article 
CAS 
PubMed 

Google Scholar 
Sauerborn, J. Legumes used for weed control in agroecosystems in the tropics. Plant Res. Dev. 50, 74–82 (1999).
Google Scholar 
Ejeta, G. & Butler, L. G. Host-parasite interactions throughout the Striga life cycle, and their contributions to Striga resistance. Afr. Crop Sci. J. 1(2), 75–80. https://doi.org/10.4314/acsj.v1i2.69889 (1993).Article 

Google Scholar 
Carsky, R. J., Singh, L. & Ndikawa, R. Suppression of Striga hermonthica on sorghum using a cowpea intercrop. Exp. Agric. 30(3), 349–358. https://doi.org/10.1017/s0014479700024467 (1994).Article 

Google Scholar 
Hsiao, A. I., Worsham, A. D. & Moreland, D. E. Effects of temperature and dl-strigol on seed conditioning and germination of witchweed (Striga asiatica). Ann. Bot. 61(1), 65–72. https://doi.org/10.1093/oxfordjournals.aob.a087528 (1988).Article 
CAS 

Google Scholar 
Patterson, D. T. Effects of Environment on Growth and Reproduction of Witchweed and the Ecological Range of Witchweed (Monograph Series of the Weed Science Society of America, 1990).Stewart, G. R. & Press, M. C. The physiology and biochemistry of parasitic angiosperms. Annu. Rev. Plant Biol. 41(1), 127–151. https://doi.org/10.1146/annurev.pp.41.060190.001015 (1990).Article 
CAS 

Google Scholar 
Anil, L., Park, R. H. P. & Miller, F. A. Temperate intercropping of cereals for forage: A review of the potential for growth and utilization with particular reference to the UK. Grass Forage Sci. 53, 301–317 (1998).
Google Scholar 
Mamolos, A. & Kalburtji, K. Significance of allelopathy in crop rotation. J. Crop Prod. 4, 197–218 (2001).
Google Scholar 
Khan, T. D., Chung, M. I., Xuan, T. D. & Tawata, S. The exploitation of crop allelopathy in sustainable agricultural production. J. Agron. Crop Sci. 191(3), 172–184 (2005).
Google Scholar 
Cissoko, M., Boisnard, A., Rodenburg, J., Press, M. C. & Scholes, J. D. New Rice for Africa (NERICA) cultivars exhibit different levels of post-attachment resistance against the parasitic weeds Striga hermonthica and Striga asiatica. New Phytol. 192(4), 952–963 (2011).CAS 
PubMed 

Google Scholar 
Rodenburg, J. et al. Do NERICA rice cultivars express resistance to Striga hermonthica (Del.) Benth. and Striga asiatica (L.) Kuntze under field conditions?. Field Crop Res. 170, 83–94 (2015).
Google Scholar 
Randrianjafizanaka, M. T., Autfray, P., Andrianaivo, A. P., Ramonta, I. R. & Rodenburg, J. Combined effects of cover crops, mulch, zero-tillage and resistant varieties on Striga asiatica (L.) Kuntze in rice-maize rotation systems. Agric. Ecosyst. Environ. 256, 23–33 (2018).
Google Scholar 
Rodenburg, J. et al. Genetic variation and host–parasite specificity of Striga resistance and tolerance in rice: The need for predictive breeding. New Phytol. 214(3), 1267–1280. https://doi.org/10.1111/nph.14451 (2017).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Nickrent, D. L. & Musselman, L. J. Introduction to parasitic flowering plants. Plant Health Instr. 13(6), 300–315 (2004).
Google Scholar 
Parker, C. Parasitic weeds: A world challenge. Weed Sci. 60(2), 269–276 (2012).CAS 

Google Scholar 
Shai Vaingast 2014. im2graph. Retrieved from: https://www.im2graph.co.il/free-downloads/windows-3264bit/ (2014).Google Maps 2021. https://maps.google.com [Accessed February 2021–December 2022].Kambach, S. et al. Consequences of multiple imputation of missing standard deviations and sample sizes in meta-analysis. Ecol. Evol. 10(20), 11699–11712 (2020).PubMed 
PubMed Central 

Google Scholar 
Nakagawa, S. & Freckleton, R. P. Missing inaction: The dangers of ignoring missing data. Trends Ecol. Evol. 23(11), 592–596 (2008).PubMed 

Google Scholar 
Idris, N. R. N. & Robertson, C. The effects of imputing the missing standard deviations on the standard error of meta analysis estimates. Commun. Stat. Simul. Comput. 38(3), 513–526. https://doi.org/10.1080/03610910802556106 (2009).Article 
MathSciNet 
MATH 

Google Scholar 
van Buuren, S. & Groothuis-Oudshoorn, K. mice: Multivariate imputation by chained equations in R. J. Stat. Softw. 45, 1–67 (2011).
Google Scholar 
van Buuren, S. Flexible Imputation of Missing Data (CRC Press, 2018).MATH 

Google Scholar 
Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315. https://doi.org/10.1002/joc.5086 (2017).Article 

Google Scholar 
O’Donnell, M. S. & Ignizio, D. A. Bioclimatic predictors for supporting ecological applications in the conterminous United States. US Geol. Surv. Data Ser. 691(10), 4–9 (2012).
Google Scholar 
Reuter, H. I., Nelson, A. & Jarvis, A. An evaluation of void filling interpolation methods for SRTM data. Int. J. Geogr. Inf. Sci. 21(9), 983–1008 (2007).
Google Scholar 
CGIAR—Consortium for Spatial Information. http://srtm.csi.cgiar.org © 2004–2021. Accessed September 19, 2021, via: http://srtm.csi.cgiar.org/srtmdata/.QGIS Development Team. QGIS Geographic Information System http://qgis.osgeo.org (Open Source Geospatial Foundation Project, 2020).Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest Package: Tests in linear mixed effects models. J. Stat. Softw. 82(13), 26. https://doi.org/10.18637/jss.v082.i13 (2017).Article 

Google Scholar 
Song, C., Peacor, S. D., Osenberg, C. W. & Bence, J. R. An assessment of statistical methods for non-independent data in ecological meta-analyses. Ecology 101(12), e03184. https://doi.org/10.1002/ecy.3184 (2020).Article 
PubMed 

Google Scholar 
Del Rey, A. C. compute.es: Compute Effect Sizes. R package version 0.2-2. https://cran.r-project.org/package=compute.es (2013).R Core Team. R: A language and environment for statistical computing. http://www.R-project.org/ (R Foundation for Statistical Computing, 2020).Wickham, H., Francois, R., Henry, L. & Müller, K. dplyr: A grammar of data manipulation. R package version 0.4. 3 (2015)Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67(1), 48. https://doi.org/10.18637/jss.v067.i01 (2015).Article 

Google Scholar 
Liebman, M. & Dyck, E. Crop rotation and intercropping strategies for weed management. Ecol. Appl. 3(1), 92–122 (1993).PubMed 

Google Scholar 
Pumariño, L. et al. Effects of agroforestry on pest, disease and weed control: A meta-analysis. Basic Appl. Ecol. 16(7), 573–582 (2015).
Google Scholar 
Kuyah, S., Whitney, C. W., Jonsson, M., Sileshi, G. W., Öborn, I., Muthuri, C. W. & Luedeling, E. Agroforestry delivers a win-win solution for ecosystem services in sub-Saharan Africa. A meta-analysis (2019).Evidente, A., Fernandez-Aparicio, M., Andolfi, A., Rubiales, D. & Motta, A. Trigoxazonane, a mono-substituted trioxazonane from Trigonella foenum-graecum root exudates, inhibits Orobanche crenata seed germination. Phytochemistry 68, 2487–2492 (2007).CAS 
PubMed 

Google Scholar 
Khan, Z. R. et al. Control of witchweed Striga hermonthica by intercropping with Desmodium spp., and the mechanism defined as allelopathic. J. Chem. Ecol. 28(9), 1871–1885 (2002).CAS 
PubMed 

Google Scholar 
Nakagawa, S. et al. Methods for testing publication bias in ecological and evolutionary meta-analyses. Methods Ecol. Evol. 13(1), 4–21 (2022).
Google Scholar 
Bakker, A. et al. Beyond small, medium, or large: Points of consideration when interpreting effect sizes. Educ. Stud. Math. 102(1), 1–8 (2019).
Google Scholar 
Scott, D. et al. Mapping the drivers of parasitic weed abundance at a national scale: A new approach applied to Striga asiatica in the mid-west of Madagascar. Weed Res. 60(5), 323–333 (2020).
Google Scholar 
Scott, D. et al. Identifying existing management practices in the control of Striga asiatica within rice–maize systems in mid-west Madagascar. Ecol. Evol. 11(19), 13579–13592 (2021).PubMed 
PubMed Central 

Google Scholar 
Rubiales, D. & Fernández-Aparicio, M. Innovations in parasitic weeds management in legume crops. A review. Agron. Sustain. Dev. 32(2), 433–449 (2012).CAS 

Google Scholar 
Bir, M. S. H. et al. Weed population dynamics under climatic change. Weed Turfgrass Sci. 3(3), 174–182 (2014).
Google Scholar 
Mohamed, K. I., Bolin, J. F., Musselman, L. J. & Townsend, P. A. Genetic diversity of Striga and implications for control and modelling future distributions. In Integrating New Technologies for Striga Control—Towards Ending the Witch-Hunt (eds Ejeta, G. & Gressel, J.) 71–84 (World Scientific, 2007).
Google Scholar 
Mandumbu, R., Mutengwa, C. S., Mabasa, S. & Mwenje, E. Predictions of the Striga scourge under new climate in southern Africa. J. Biol. Sci. 17, 192–201. https://doi.org/10.3923/jbs.2017.194.201 (2017).Article 

Google Scholar 
Mudereri, B. T. et al. Multi-source spatial data-based invasion risk modelling of Striga (Striga asiatica) in Zimbabwe. GIScience Remote Sens. 57(4), 553–571. https://doi.org/10.1080/15481603.2020.1744250 (2020).Article 

Google Scholar  More

Réale, D., Reader, S. M., Sol, D., McDougall, P. T. & Dingemanse, N. J. Integrating animal temperament within ecology and evolution. Biol. Rev. 82, 291–318 (2007).PubMed 

Google Scholar 
Kaiser, M. I. & Müller, C. What is an animal personality?. Biol. Philos. 36, 1 (2021).
Google Scholar 
Sih, A., Bell, A. & Johnson, J. C. Behavioral syndromes: An ecological and evolutionary overview. Trends Ecol. Evol. 19, 372–378 (2004).PubMed 

Google Scholar 
Gosling, S. D. From mice to men: What can we learn about personality from animal research?. Psychol. Bull. 127, 45–86 (2001).PubMed 

Google Scholar 
Biro, P. A. & Stamps, J. A. Are animal personality traits linked to life-history productivity?. Trends Ecol. Evol. 23, 361–368 (2008).PubMed 

Google Scholar 
Koolhaas, J. M. Coping style and immunity in animals: Making sense of individual variation. Brain Behav. Immun. 22, 662–667 (2008).PubMed 

Google Scholar 
Réale, D. et al. Personality and the emergence of the pace-of-life syndrome concept at the population level. Philos. Trans. R. Soc. B Biol. Sci. 365, 4051–4063 (2010).
Google Scholar 
Stamps, J. A. Growth-mortality tradeoffs and ‘personality traits’ in animals. Ecol. Lett. 10, 355–363 (2007).PubMed 

Google Scholar 
Finkemeier, M. A., Langbein, J. & Puppe, B. Personality research in mammalian farm animals: Concepts, measures, and relationship to welfare. Front Vet. Sci. 10(5), 355–363 (2018).
Google Scholar 
Murphy, E., Nordquist, R. E. & van der Staay, F. J. A review of behavioural methods to study emotion and mood in pigs. Sus. Scrofa. Appl. Anim. Behav. Sci 159, 9–28 (2014).
Google Scholar 
Lauber, M. C. Y., Hemsworth, P. H. & Barnett, J. L. The effects of age and experience on behavioural development in dairy calves. Appl. Anim. Behav. Sci. 99, 41–52 (2006).
Google Scholar 
Neave, H. W., Costa, J. H. C., Weary, D. M. & von Keyserlingk, M. A. G. Personality is associated with feeding behavior and performance in dairy calves. J. Dairy Sci. 101, 7437–7449 (2018).PubMed 

Google Scholar 
Foris, B., Zebunke, M., Langbein, J. & Melzer, N. Evaluating the temporal and situational consistency of personality traits in adult dairy cattle. Plos One 13, e0204619 (2018).PubMed 
PubMed Central 

Google Scholar 
Dingemanse, N. J. & Dochtermann, N. A. Quantifying individual variation in behaviour: Mixed-effect modelling approaches. J. Anim. Ecol. 82, 39–54 (2013).PubMed 

Google Scholar 
Dingemanse, N. J., Kazem, A. J. N., Réale, D. & Wright, J. Behavioural reaction norms: Animal personality meets individual plasticity. Trends Ecol. Evol. 25, 81–89 (2010).PubMed 

Google Scholar 
Nakagawa, S. & Schielzeth, H. Repeatability for Gaussian and non-Gaussian data: A practical guide for biologists. Biol. Rev. https://doi.org/10.1111/j.1469-185X.2010.00141.x (2010).Article 
PubMed 

Google Scholar 
Bell, A. M., Hankison, S. J. & Laskowski, K. L. The repeatability of behaviour: A meta-analysis. Anim. Behav. 77, 771–783 (2009).PubMed 
PubMed Central 

Google Scholar 
Neave, H. W., Costa, J. H. C., Benetton, J. B., Weary, D. M. & von Keyserlingk, M. A. G. Individual characteristics in early life relate to variability in weaning age, feeding behavior, and weight gain of dairy calves automatically weaned based on solid feed intake. J. Dairy Sci. 102, 10250–10265 (2019).PubMed 

Google Scholar 
Berckmans, D. Precision livestock farming technologies for welfare management in intensive livestock systems. Rev. Sci. Tech. OIE 33, 189–196 (2014).
Google Scholar 
Carslake, C., Vázquez-Diosdado, J. A. & Kaler, J. Machine learning algorithms to classify and quantify multiple behaviours in dairy calves using a sensor: Moving beyond classification in precision livestock. Sensors 21, 88 (2020).ADS 
PubMed Central 

Google Scholar 
Hertel, A. G., Niemelä, P. T., Dingemanse, N. J. & Mueller, T. A guide for studying among-individual behavioral variation from movement data in the wild. Mov. Ecol. 8(1), 1–18 (2020).
Google Scholar 
Occhiuto, F., Vázquez-Diosdado, J. A., Carslake, C. & Kaler, J. Personality and predictability in farmed calves using movement and space-use behaviours quantified by ultra-wideband sensors. R. Soc. Open Sci. 9, 212019 (2022).ADS 
PubMed 
PubMed Central 

Google Scholar 
Carslake, C., Occhiuto, F., Vázquez-Diosdado, J. A. & Kaler, J. Repeatability and predictability of calf feeding behaviors—quantifying between- and within-individual variation for precision livestock farming. Front. Vet. Sci. https://doi.org/10.3389/fvets.2022.827124 (2022).Article 
PubMed 
PubMed Central 

Google Scholar 
Tolkamp, B. J. & Kyriazakis, I. To split behaviour into bouts, log-transform the intervals. Anim. Behav. 57, 807–817 (1999).PubMed 

Google Scholar 
Houslay, T. M. & Wilson, A. J. Avoiding the misuse of BLUP in behavioural ecology. Behav. Ecol. 28, 948 (2017).PubMed 
PubMed Central 

Google Scholar 
R Core Team. R. Preprint at (2021).Bürkner, P.-C. Advanced bayesian multilevel modeling with the R package brms. R. J. 10, 395 (2018).
Google Scholar 
Dancey, C. P. & Reidy, J. Statistics without maths for psychology (Pearson education, 2007).
Google Scholar 
von Keyserlingk, M. A. G., Brusius, L. & Weary, D. M. Competition for teats and feeding behavior by group-housed dairy calves. J. Dairy Sci. 87, 4190–4194 (2004).
Google Scholar 
Fraley, R. C. & Roberts, B. W. Patterns of continuity: A dynamic model for conceptualizing the stability of individual differences in psychological constructs across the life course. Psychol. Rev. 112, 60–74 (2005).PubMed 

Google Scholar 
Ashcroft, J., Semmler, C., Carnell, S., van Jaarsveld, C. H. M. & Wardle, J. Continuity and stability of eating behaviour traits in children. Eur. J. Clin. Nutr. 62, 985–990 (2008).PubMed 

Google Scholar 
Neave, H. W., Costa, J. H. C., Weary, D. M. & von Keyserlingk, M. A. G. Long-term consistency of personality traits of cattle. R. Soc. Open Sci. 7, 191849 (2020).ADS 
PubMed 
PubMed Central 

Google Scholar 
Müller, R. & von Keyserlingk, M. A. G. Consistency of flight speed and its correlation to productivity and to personality in Bos taurus beef cattle. Appl. Anim. Behav. Sci. 99, 193–204 (2006).
Google Scholar 
Neja, W., Sawa, A., Jankowska, M., Bogucki, M. & Krężel-Czopek, S. Effect of the temperament of dairy cows on lifetime production efficiency. Arch. Anim. Breed 58, 193–197 (2015).
Google Scholar 
Haskell, M. J., Simm, G. & Turner, S. P. Genetic selection for temperament traits in dairy and beef cattle. Front Genet. https://doi.org/10.3389/fgene.2014.00368 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Whalin, L., Neave, H. W., Føske Johnsen, J., Mejdell, C. M. & Ellingsen-Dalskau, K. The influence of personality and weaning method on early feeding behavior and growth of Norwegian red calves. J. Dairy Sci. 105, 1369–1386 (2022).PubMed 

Google Scholar 
Dammhahn, M., Dingemanse, N. J., Niemelä, P. T. & Réale, D. Pace-of-life syndromes: A framework for the adaptive integration of behaviour, physiology and life history. Behav. Ecol. Sociobiol. 72(3), 1–8 (2018).
Google Scholar 
Kelly, D. N. et al. Large variability in feeding behavior among crossbred growing cattle. J. Anim. Sci. 98, 1–10 (2020).
Google Scholar 
Neave, H. W., Weary, D. M. & von Keyserlingk, M. A. G. Review: Individual variability in feeding behaviour of domesticated ruminants. Animal 12, S419–S430 (2018).PubMed 

Google Scholar 
DeVries, T. J., von Keyserlingk, M. A. G., Weary, D. M. & Beauchemin, K. A. Measuring the feeding behavior of lactating dairy cows in early to peak lactation. J. Dairy Sci. 86, 3354–3361 (2003).PubMed 

Google Scholar 
Kelly, D. N., Sleator, R. D., Murphy, C. P., Conroy, S. B. & Berry, D. P. Phenotypic and genetic associations between feeding behavior and carcass merit in crossbred growing cattle. J. Anim. Sci. 99, skab285 (2021).PubMed 

Google Scholar 
Weary, D. M., Huzzey, J. M. & von Keyserlingk, M. A. G. Board-invited review: Using behavior to predict and identify ill health in animals. J. Anim. Sci. 87, 770–777 (2009).PubMed 

Google Scholar 
Carter, A. J., Feeney, W. E., Marshall, H. H., Cowlishaw, G. & Heinsohn, R. Animal personality: What are behavioural ecologists measuring?. Biol. Rev. 88, 465–475 (2013).PubMed 

Google Scholar 
Biro, P. A. Do rapid assays predict repeatability in labile (behavioural) traits?. Anim Behav 83, 1295–1300 (2012).
Google Scholar 
Percie du Sert, N. et al. Reporting animal research: Explanation and elaboration for the ARRIVE guidelines 20. Plos Biol. 18, e3000411 (2020).PubMed 
PubMed Central 

Google Scholar  More

Study siteWe conducted laboratory studies at the field site in Finca Irlanda, which is a 300-hectare organic shaded coffee farm located at 1100-m altitude, in the municipality of Tapachula, the state of Chiapas in Southern Mexico (92° 20′ 29″ W and 15° 10′ 65″ N). For the laboratory experiments, all organisms were freshly collected from Finca Irlanda or reared in the lab from insects collected from the field close by. The lab and field work was performed with a permit from the farm owner the Peters family.Ant aggression testTo examine the effect of phorid flies (P. lascinosus) on the aggressivity of ants (A. sericeasur) towards the parasitoids of the beetle larvae (A. orbigera), we conducted an ant aggression test with two treatments: with and without phorids. In the first treatment, a small coffee branch containing two leaves with scale insects (C. viridis) and 20 ant workers were both introduced into a one-liter plastic container. This was done to mimic as much as possible field conditions where the ants are tending scale insects. After waiting for at least 15 min for the ants to calm down and start tending the scale insects, one third- or fourth-instar larva of the beetle was introduced. In the second treatment, all settings were the same except for the addition of 3–4 phorid flies. Once the two treatments were established, one female parasitoid wasp (H. shuvakhinae) was released into each container. During a forty-minute trial, each time that a parasitoid wasp encountered an ant worker, the response of the ant individual was recorded. Ant responses to parasitoids were classified into two categories: (1) the ant ignores the wasp; (2) the ant attacks the wasp. All insects were used for a single replicate and then discarded. A total of four replicates were completed for both the presence and absence of phorids. For each trial, we calculated the proportion of actions (either aggressive or none) by ants when encountering the parasitoid wasp in the treatments with and without phorid flies. We used R36 to conduct a two-sample Mann–Whitney U test on the proportion of ant actions.Parasitism experiments and analysesTo examine the parasitoid wasp’s host preference and the effect of the 1st degree and the 2nd degree HOIs on the beetle’s parasitism and sex ratio, we conducted a laboratory experiment in insect tents (60 cm × 60 cm × 60 cm) with three treatments: (1) no ants (no HOIs but only the wasp and the beetle larvae), (2) ants (1st degree HOI), and (3) ants and phorids (1st and 2nd degree HOIs) (Fig. 1-B). We randomly assigned insect tents to each treatment in each trial, and the tents for each treatment were also shuffled in each trial. All beetle larvae used for these experiments were reared in the laboratory for at least two generations from freshly collected beetle adults. In each tent we placed a coffee branch with 4–6 leaves infested with approximately 100 scale insects inside a plastic container at the center of an insect tent. The set up for the three treatments of species combinations were as follows: (1) 4–5 third or fourth instar beetle larvae and a parasitoid wasp; (2) 4–5 third or fourth instar beetle larvae, a parasitoid wasp, and about 60–80 ant workers; (3) 4–5 third or fourth instar beetle larvae, a parasitoid wasp, about 60–80 ant workers and 3–4 phorid flies. Organism densities in these treatments were close to those observed in the field. To allow for acclimation, we introduced organisms into the tents in the following order: first, we introduced the coffee branch containing scales, immediately followed by the ants (in treatments 2 and 3). After the ants settled down and started tending the scale insects, we introduced the beetle larvae. Once the larvae began moving on the coffee leaves, we introduced the phorids (in treatment 3). When the three treatments were established, and the organisms exhibit normal behavior, we released one lab-reared female parasitoid wasp (H. shuvakhinae) in each tent (treatments 1, 2, and 3). We allowed the organisms to interact for 24 h. After 24 h, we collected all beetle larvae in each treatment and reared them with sufficient scale insects as food, until beetle adults emerged or parasitism symptoms appeared (parasitized larvae turned into hardened black mummies). The treatments of no HOI and 1st + 2nd degree HOI were repeated for 10 consecutive times, and the treatment of 1st degree HOI was repeated for 11 consecutive times, with new individuals of each organism. We recorded parasitism instances and beetle sexes upon emergence. To estimate the sex ratio without parasitoid influence, 78 randomly selected beetle individuals were reared on coffee leaves with scale insects without any interaction with other organisms.To analyze the effect of the parasitoid, the ant and the phorid fly on the parasitism rate and the sex ratio of the beetle, we developed a nested model, starting from$$logitleft(widehat{P}(S)right)=a+bA$$where (widehat{P}(S)) is the probability of an individual being parasitized, A is a binary variable, standing for the absence (0) and presence (1) of ants, a is the baseline probability of parasitism, and b is the magnitude of parasitism altered by ants in the logistic function. We further hypothesized that phorid attacks modify the strength of the interaction modification that ants exert upon the host-parasitoid interaction. Therefore,$$b=g+hP$$where P is another binary variable, standing for the presence (1) and absence (0) of phorids. Substituting b, we obtain the following function,$$logitleft(widehat{P}(S)right)=a+gA+hAP$$where g represents the effect of ants on the parasitism rate of A. orbigera larvae, and h represents the effect of the fly’s facilitation, via interfering with the ant’s interference on the parasitism rate of A. orbigera larvae. We used binary responses (1: survival; 0: parasitized) of all available beetle individuals across the three treatments. We performed model selection based on the Akaike Information Criterion (AIC) and likelihood ratio tests. For the latter, we started model selection by fitting the full model and preceding each step by eliminating the term that had the least significance (the greatest p-value) on the explanation of the dependent variable. The analysis was performed with the application of the bbmle package in R. By doing this, we determined the maximum likelihood estimates of survival probability of the beetle, (widehat{P}(S)), in the three treatments: (1) A = 0, AP = 0 (no HOI); (2) A = 1, AP = 0 (one HOI: ant interference) and (3) A = 1, AP = 1 (interacting HOIs: phorid interference with ant interference), and errors associated with these estimates.The same idea applies to the sex ratio of the beetle under the influence of various organisms. We developed the following equation,$$logitleft(widehat{P}(F|S)right)= r+mA+nAP$$where (widehat{P}(F|S)) is the probability of a parasitism survivor being female. A and P are both binary variables. Respectively, they represent the ant and the phorid fly, and the numeric attributes, 0 and 1, denote their absence and presence. As before, model selection and parameter estimates were conducted with AIC. By doing this, we determined (widehat{P}(F|S)), the estimate of being a female beetle given survival, for the three treatments: (1) A = 0, AP = 0 (no HOI); (2) A = 1, AP = 0 (one HOI: ant interference) and (3) A = 1, AP = 1 (interacting HOIs: phorid interference with ant interference), and errors associated with these estimates. We employed the mle2 function in the bbmle package in R to estimate the female probability (1) in the absence of HOI (the beetle and the parasitoid alone), (2) in the presence of the 1st degree HOI (the beetle, the parasitoid and the ant), and (3) in the presence of the 1st and the 2nd degree HOIs (the beetle, the parasitoid, the ant and the phorid fly).Probabilities of per capita female and per capita male survival from parasitism under the influence of ant and the phorid flyTo test whether the sex ratio of beetle survivors’ population is due to sex-differential survival probability, Bayes’ theorem was employed. Per capita female survival probability from parasitism in each treatment of the parasitism experiment was derived based on (widehat{P}(F)), (widehat{P}left(F|Sright),) and (widehat{P}(S)), and per capita male survival probability was derived based on (widehat{P}(M)), (widehat{P}left(M|Sright),) and (widehat{P}(S)). According to the Central Limit Theorem, the estimates of proportions, (widehat{P}left(S|Fright)) and (widehat{P}left(S|Mright)), are approximately normally distributed,$$widehat{P}left(S|Fright)sim Nleft(widehat{P}left(S|Fright), sqrt{frac{widehat{P}(S|F)times left(1-widehat{P}left(S|Fright)right)}{{n}^{*}}}right)$$$$widehat{P}left(S|Mright)sim Nleft(widehat{P}left(S|Mright), sqrt{frac{widehat{P}(S|M)times left(1-widehat{P}left(S|Mright)right)}{{n}^{*}}}right)$$with means (widehat{P}left(S|Fright)) and (widehat{P}(S|M)), and standard deviations (sqrt{frac{widehat{P}left(S|Fright)times (1-widehat{P}left(S|Fright))}{{n}^{*}}}) and (sqrt{frac{widehat{P}left(S|Mright)times (1-widehat{P}left(S|Mright))}{{n}^{*}}}), where (widehat{P}(S|F)) and (widehat{P}(S|M)), respectively, are the population proportions of females and males. Here we employ n*, the smallest sample size among those of the three variables in the Bayesian formulas for males and females. Since the three variables have different sample sizes, n* guarantees a conservative estimate of standard error, and thus confidence interval, of each derived probability. More