Data and R codes of "The effects of human population density on trophic interactions are contingent upon latitude"
Authors: Juan A. Hernández-Agüero, Ildefonso Ruiz-Tapiador3, Lucas A. Garibaldi4,5, Mikhail V. Kozlov6, Elina Mäntylä6,7,8, Marcos E. Nacif4,5, Norma Salinas9 & Luis Cayuela1

Methodology:
Herbivory
Insect herbivory data was obtained from the literature. Insect herbivory was quantified in the papers used as the proportion of leaf area consumed by defoliating insects. It was estimated with either image software or visually, as proposed by Alliende (1989). Because human population density is changing rapidly, we restricted our analysis to data collected since 2000. We extracted the greater part of data from a systematic global review (Kozlov et al., 2015), and we searched for data published between 2015 and December 2021 in the ISI Web of Science on 19 April 2022 with the words “Defoliation” and “Insect herbivory [i.e. Defoliation (All Fields) AND insect herbivory (All Fields) and 2021 or 2020 or 2019 or 2018 or 2017 or 2016 or 2015 (Publication Years)]”. "Defoliation” was used to avoid including other types of herbivory such as galler herbivory or miner herbivory, while “insect herbivory” was included to avoid papers referring to large herbivorous. Galler and miner herbivory were excluded from the search criteria because the data paper of Mendes et al., (2021) did not include them. We only included data of herbivory on woody plants following the searching criteria of Kozlov et al. (2015). From the 131 articles obtained only eight used the selected methodology to quantify herbivory, and we also included the submitted manuscript by Hernández-Agüero, Ruiz-Tapiador, Garibaldi, Kozlov, Nacif et al. (2023) (Figure 1a; Supplementary Material Table S1). A list of the data sources (10 sources) is found in Supplementary Material Table S1 and a complete list of studies reviewed (130 sources) and the justification for exclusion can be found in Supplementary material 3.
Predation
Bird predation on herbivorous insects was quantified by attack rates on artificial caterpillars made of odorless plasticine, which were placed on a woody plant branch or leaves. Following an exposure (from 2 to 64 days), the caterpillars were revisited and bird marks were counted following the methodology proposed in Low, Sam, McArthur, Posa & Hochuli (2014). We extracted the data from two global studies (Roslin et al., 2017; Zvereva et al., 2019) and searched for additional studies (published between January 2000 and 23 December of 2021) in the ISI Web of Science on 19 April 2022, using the keywords “bird”, “predation” and “larvae” [i.e. Results for bird (All Fields) AND predation (All Fields) AND ( larvae (All Fields) OR caterpillar (All Fields)) and 2021 or 2020 or 2019 or 2018 or 2017 or 2016 or 2015 or 2014 or 2013 or 2012 or 2011 or 2010 or 2009 or 2008 or 2007 or 2006 or 2005 or 2004 or 2003 or 2002 or 2001 or 2000 (Publication Years)]. Most of the studies identified included caterpillars of different colors, except for Roslin et al. (2017) which was based only on green-colored caterpillars. Because prey colors can influence attack rates (Hernández-Agüero, Polo, García, Simón, Ruiz-Tapiador et al., 2020) and colour preferences by birds can change latitudinally, we selected only green caterpillar data, which does not show changes in preference by predators along latitudinal gradients (Zvereva et al., 2019). Our search identified 217 publications, among which only 11 studies shared similar methods to those of Roslin et al. (2017) and Zvereva et al. (2019) and provide raw data on green plasticine models; we also included the submitted manuscript article by Hernández-Agüero et al. (2023) (Figure 1b; Supplementary Material Table S2). A list of the data sources (16 sources) is found in Supplementary Material Table S2 and a complete list of studies reviewed (320 sources) and the justification for exclusion can be found in Supplementary material 3.
To account for differences across studies in the length of the study period, we estimated the probability of a caterpillar to be attacked over one day (henceforth probability of bird predation) as:
P (X=1)=1-P (X=0)=〖1-[1-(N/T)]〗^((1/t) )
where P (X≥1) is the probability of having one or more caterpillars attacked by birds, N is the number of attacked caterpillars, T is the total number of caterpillars used per period, and t is the period length in days.
Human population
We used the Gridded Population of the World (GPWv4 2020) dataset, at 1 km spatial resolution, to estimate human population density at each study site. One site was consider when has unique coordinates. To do so, we created a 10 km radius buffer at each site with the ‘st_buffer’ function of the ‘sf’ package (Pebesma, 2018). This buffer size has been used in other studies to investigate the effects of human impact on vertebrate and plant species (Kim, Mizuno & Kobayashi, 2003; Pautasso, 2007) and is highly correlated (r ≥ 0.7) with values of population density obtained at smaller (e.g. 5 km radius buffer) and larger (e.g. 25 or 40 km radius buffer) scales (Supplementary Material Figure S1 and Figure S2). We then cropped the rasterized GPWv4 dataset with our site buffer of 10 km using the ‘crop’ function of the R package ‘raster’ (Hijmans, 2020) and created a new raster object using our buffered coordinates and the GPWv4 cropped dataset with the ‘mask’ function. Finally, we extracted the total human population at every site with the ‘extract’ function and then human population density per km2 was calculated by dividing by the area in km2. Correlation tests were made between the data obtained for HPD and other common used index to measure urbanization (i.e. proportion of “built area”). Built area data was obtained from Dynamic World V1 (Brown et al., 2022)
Data analyses
We used generalized linear mixed models (GLMM) with a beta error distribution and a logit link function to investigate the effects of human population and latitude on bird predation and insect herbivory. Beta is a family of continuous probability distributions defined on the interval [0,1], and therefore appropriate for the type of response variables we are modelling in this study. We used the natural logarithm of human population, absolute (i.e. unsigned) latitude and their interaction as predictors. In the GLMMs for herbivory, we included the following random factors: i) site, which accounted for potential spatial autocorrelation (there were from 1 to 49 data per site); and ii) plant species nested within genus, which accounted for differences in palatability and plant defenses against herbivory among plant taxa. In the GLMMs for bird predation, we included site as random factor (there were from 1 to 12 data per location). All GLMMs were fitted using maximum likelihood with the function ‘mixed_model’ of the R package ‘glmmTMB’ (Brooks, Kristensen, van Benthem, Magnusson, Berg et al., 2017).
For both herbivory and predation, alternative models were compared using the Akaike information criterion (AIC) to select the explanatory variables (i.e. fixed effects). Models with a difference in AIC > 2 indicated that the worst model could be omitted. Following Nakagawa & Schielzeth (2013), we estimated the R2 of all plausible linear or mixed models. This allowed two components of R2 to be calculated: (1) a marginal R2 (R2m) that only considers the variability explained by fixed effects; and (2) a conditional R2 (R2c) that accounts for the variability supported by both fixed and random effects. Model residuals were explored using a simulation-based approach to create readily interpretable scaled (quantile) residuals for the fitted GLMMs (Hartig, 2019). Moran’s index was used to estimate spatial autocorrelation in model residuals both for proportion of herbivory and probability of predation. To test for significance of spatial autocorrelation, this index was compared with a null model random distribution using the ‘spdep’ package (Bivand & Wong, 2018). When spatial autocorrelation in model residuals was significant, we re-fitted the model including a spatial autocorrelation function with exponential correlation structure. This uses an Euclidean distance matrix based on site coordinates (Brooks et al., 2017). We repeated all the analyses with human population density calculated at different buffer sizes (i.e. 1, 5, 10, 15, 20, 25, 30, 35 and 40 km radius buffer) to ensure that the spatial scales were not affecting our results (Supplementary Material Table S3 and Table S4).
Finally, mean annual temperature between 1970-2000 in each study site was estimated using WorldClim database obtained with the function ‘getData’ from ‘raster’ package for every coordinate. Correlation tests were made using the function ‘corrgram’ from ‘corrgram’ package (Wright, 2021) between the temperature and the absolute latitude.
## Description of the data and file structure
03-07-2023.R This is the main R script of the paper.
calculacion.R This is the secondary R script of the paper. It is a function to calculate human population density. It is needed to run 03-07-2023.R
herbivory2022.csv this is the herbivory dataset to be used in 03-07-2023.R
"site": is the study site named by the authors of the original paper
"lat": is the latitude (epsg 6933)
"long": is the longitude (epsg 6933)
"Country": is the Country (in english) where the study was made
"Species": is the plant species in which herbivory was evaluated
"Genus": is the Plant Genus in which herbivory was evaluated
"pop": is the human population 10 km arround the study site
"source": is the reference of the paper from which the data was obtained
"logpop": is the logarithm 10 of the human population 10 km arround each study site
"abslat": is the absolute latitude (epsg 6933)
"HPopulation": is the human population density per squarekm
"percentage": is the mean percentage of defoliation
Green.png This image created by Inés M. Alonso-Crespo is needed to create one of the figures of the paper in 03-07-2023.R
predation2022.csv this is the predation dataset to be used in 03-07-2023.R
"Site": is the study site named by the authors of the original paper
"SiteYear": is the study site named by the authors of the original paper and the year in which the experiment was made
"lat": is the latitude (epsg 6933)
"long": is the longitude (epsg 6933)
"Population": is the human population 10 km arround the study site
"Source": is the reference of the paper from which the data was obtained
"Hpopulation": is the human population 10 km arround the study site
"dat": is the number of days since the placement of the plasticine caterpillars and the evaluation of damage
"Damaged": is the number of caterpillars damaged
"Placed": is the number of caterpillars placed
"Rate": is the rate of caterpillars damaged
"Probability": is the probability of a caterpillar to be damaged one day
"Hpopkm": is the human population density per km square
"logHpopkm": is the logartihm of the human population density per km square
"abslat": is the absolute latitude (epsg 6933)
Defol.png This image created by Inés M. Alonso-Crespo is needed to create one of the figures of the paper in 03-07-2023.R
## Sharing/Access information
Data was derived from the following sources:
References:
- Allman, B. P., Kielland, K., & Wagner, D. (2018). Leaf herbivory by insects during summer reduces overwinter browsing by moose. BMC ecology, 18(1), 1-7.
- Guyot, V., Castagneyrol, B., Vialatte, A., Deconchat, M., & Jactel, H. (2016). Tree diversity reduces pest damage in mature forests across Europe. Biology Letters, 12(4), 20151037.
- Hernández-Agüero, J.A., Ruiz-Tapiador, I., Cosio, E., Garibaldi, L.A., Kozlov, M.V., Nacif, M.E., Salinas, N., Zverev, V., Zvereva, E.,. and Cayuela, L., Forest Management Affects Ecosystem Functioning (Predation and Herbivory) But Not Ecosystem Stability: A Comparative Study Across Four Forest Ecosystems Around the World. Available at SSRN: https://ssrn.com/abstract=4480117 or http://dx.doi.org/10.2139/ssrn.4480117
- Kozlov, M. V., Lanta, V., Zverev, V. & Zvereva, E. L. (2015). Global patterns in background losses of woody plant foliage to insects. Global Ecology and Biogeography, 24(10): 1126-1135.
- Mendes, G.M., Silveira, F.A., Oliveira, C., Dáttilo, W., Guevara, R., Ruiz‐Guerra, B., Boaventura, M.G., Sershen, , Ramdhani, S., Phartyal, S.S., Ribeiro, S.P., Pinto, V.D., Vasconcelos, H., Tito, R., Pereira, C.C., Carvalho, B., Carvalho, G.M., del Val, E., Buisson, E., Arruda, A., Toth, J.‐B., Roque, F., Souza, A.H., Bolzan, F., Neves, F., Kuchenbecker, J., Demetrio, G.R., Seixas, L., Romero, G., de Omena, P.M., Silva, J.O., Paolucci, L., Queiroz, E., Ooi, M.K., Mills, C.H., Gerhold, P., Merzin, A., Massante, J.C., Aguilar, R., Carbone, L.M., Campos, R., Gomes, I., Zorzal, G., Solar, R., Ramos, L., Sobrinho, T., Sanders, P. & Cornelissen, T. (2021), How much leaf area do insects eat? A data set of insect herbivory sampled globally with a standardized protocol. Ecology. E03301
- Muehleisen, A. J., Engelbrecht, B. M., Jones, F. A., Manzané‐Pinzón, E., & Comita, L. S. (2020), Local adaptation to herbivory within tropical tree species along a rainfall gradient. Ecology, 101(11), e03151.
- Nakajima, H. (2015), Defoliation by gypsy moths negatively affects the production of acorns by two Japanese oak species. Trees, 29, 1559-1566.
- Nakajima, H. (2018), Refoliation of deciduous canopy trees following severe insect defoliation: comparison of Fagus crenata and Quercus crispula. Plant Ecology, 219(6), 665-675.
- Valoy, M., Ordano, M., Bernacki, F., Palacio, F. X., López-Acosta, J. C., & Varela, O. (2018), Patrones de herbivoría en Vassobia breviflora (Solanaceae): variación en el daño foliar y selección natural mediada por herbívoros. Revista de Biología Tropical, 66(4), 1683-1700.
- Wang, B., Tian, C., & Liang, Y. (2021), Mixed effects of landscape structure, tree diversity and stand’s relative position on insect and pathogen damage in riparian poplar forests. Forest Ecology and Management, 479, 118555.
- Gomes, D.G.E., Toth, C.A., Cole, H.J., Francis, C.D. & Barber, J.R. (2021). Phantom rivers filter birds and bats by acoustic niche. Nature communications, 12(1): 1-8.
- Hernández-Agüero, J.A., Ruiz-Tapiador, I., Cosio, E., Garibaldi, L.A., Kozlov, M.V., Nacif, M.E., Salinas, N., Zverev, V., Zvereva, E.,. and Cayuela, L., Forest Management Affects Ecosystem Functioning (Predation and Herbivory) But Not Ecosystem Stability: A Comparative Study Across Four Forest Ecosystems Around the World. Available at SSRN: https://ssrn.com/abstract=4480117 or http://dx.doi.org/10.2139/ssrn.4480117
- Kozlov, M.V., Lanta, V., Zverev, V., Rainio, K., Kunavin, M.A. & Zvereva, E.L. (2017). Decreased losses of woody plant foliage to insects in large urban areas are explained by bird predation. Global Change Biology, 23(10): 4354-4364.
- Leuenberger, W., Cohen, J. B., Rustad, L., Wallin, K. F. & Parry, D. (2021). Short-Term Increase in Abundance of Foliage-Gleaning Insectivorous Birds Following Experimental Ice Storms in a Northern Hardwood Forest. Frontiers in Forests and Global Change, 3: 566376.
- Maas, B., Tscharntke, T., Saleh, S., Putra, D.D. & Clough, Y. (2015). Avian species identity drives predation success in tropical cacao agroforestry. Journal of Applied Ecology, 52(3): 735-743.
- Mäntylä, E., Alessio, G.A., Blande, J.D., Heijari, J., Holopainen, J.K., Laaksonen, T., Piirtola, P. & Klemola, T. (2008). From plants to birds: higher avian predation rates in trees responding to insect herbivory. PloS One, 3(7), e2832.
- Mäntylä, E., Blande, J. D. & Klemola, T. (2014). Does application of methyl jasmonate to birch mimic herbivory and attract insectivorous birds in nature?. Arthropod-Plant Interactions, 8(2): 143-153.
- Mrazova, A., & Sam, K. (2018), Application of methyl jasmonate to grey willow (Salix cinerea) attracts insectivorous birds in nature. Arthropod-Plant Interactions, 12(1), 1-8.
- Murray, K., Soh, M.C., Puan, C.L., Abu Hashim, A.K.B., Abdul Patah, P.B. & Peh, K.S.H. (2020). Predation on multiple prey types across a disturbance gradient in tropical montane forests of Peninsular Malaysia. Frontiers in Forests and Global Change, 3: 80.
- Roels, S.M., Porter, J.L. & Lindell, C.A. (2018). Predation pressure by birds and arthropods on herbivorous insects affected by tropical forest restoration strategy. Restoration ecology, 26(6): 1203-1211.
- Roslin, T., Hardwick, B., Novotny, V., Petry, W.K., Andrew, N.R., Asmus, A., Barrio, I.C., Basset, Y., Boesing, A.L., Bonebrake, T.C., Cameron, E.K., Dáttilo, W., Donoso, D.A., Drozd, P., Gray, C.L., Hik, D. S., Hill, S.J., Hopkins, T., Huang, S., Koane, B., Laird-Hopkins, B., Laukkanen, L., Lewis, O.T., Milne, S., Mwesige, I., Nakamura, A., Nell, C.S., Nichols, E., Prokurat, A., Sam, K., Schmidt, N.M., Slade, A., Slade, V., Suchanková, A., Teder, T., van Nouhuys, S., Vandvik, V., Weissflog, A., Zhukovich, V. & Slade, E.M. (2017). Higher predation risk for insect prey at low latitudes and elevations. Science, 356: 742– 744.
- Schwab, D., Wurz, A., Grass, I., Rakotomalala, A.A., Osen, K., Soazafy, M.R., Martín, D.A & Tscharntke, T. (2020). Decreasing predation rates and shifting predator compositions along a land‐use gradient in Madagascar’s vanilla landscapes. Journal of Applied Ecology, 58(2): 360-371.
- Valdés-Correcher, E., Moreira, X., Augusto, L., Barbaro, L., Bouget, C., Bouriaud, O., Branco, M., Centenaro, G., Csóka, G., Damestoy, T., Dobrosavljević, J., Duduman, M.L., Dulaurent, A.M., Eötvös, C.B., Faticov, M., Ferrante, M., Fürjes-Mikó, Á., Galmán, A., Gossner, M.M., Hampe, A., Harvey, D., Howe, A.G., Kadiri, Y., Kaennel-Dobbertin, M., Koricheva, J., Kozel, A., Kozlov, M.V., Lövei, G.L., Lupaștean, D., Milanović, S., Mrazova, A., Opgennoorth, L., Pitkänen, J.M., Popova, A., Popović, M., Prinzing, A., Queloz, V., Roslin, T., Sallé, A., Sam, K., Scherer-Lorenzen, M., Schuldt, A., Selikhovkin, A., Suominen, L., Tack, A.J.M., Tahadlova, M., Thomas, R. & Castagneyrol, B. (2021). Search for top‐down and bottom‐up drivers of latitudinal trends in insect herbivory in oak trees in Europe. Global Ecology and Biogeography, 30(3): 651-665.
- Zvereva, E.L., Castagneyrol, B., Cornelissen, T., Forsman, A., Hernández-Agüero, J.A., Klemola, T., Paolucci, L., Polo, V., Salinas., N., Theron, K.J., Xu, G., Zverev, V. & Kozlov, M. (2019). Opposite latitudinal patterns for bird and arthropod predation revealed in experiments with differently coloured artificial prey. Ecology and Evolution, 9(24): 14273e14285.
- Zvereva, E. L., & Kozlov, M. V. (2021), Seasonal variations in bird selection pressure on prey colouration. Oecologia, 196(4), 1017-1026.
- Zvereva, E.L., Zverev, V. & Kozlov, M.V. (2020). Predation and parasitism on herbivorous insects change in opposite directions in a latitudinal gradient crossing a boreal forest zone. Journal of Animal Ecology, 89(12): 2946-2957.