Philippa Kaur

Bond-Lamberty, B. & Thomson, A. Temperature-associated increases in the global soil respiration record. Nature 464, 579–582 (2010).Article 
CAS 
PubMed 

Google Scholar 
Crippa, M. et al. Food systems are responsible for a third of global anthropogenic GHG emissions. Nat. Food 2, 198–209 (2021).Article 
CAS 

Google Scholar 
Shakoor, A. et al. Effect of animal manure, crop type, climate zone, and soil attributes on greenhouse gas emissions from agricultural soils—A global meta-analysis. J. Clean Prod. 278, 124019. https://doi.org/10.1016/j.jclepro.2020.124019 (2021).Article 
CAS 

Google Scholar 
Wu, L. et al. Soil organic matter priming and carbon balance after straw addition is regulated by long-term fertilization. Soil Biol. Biochem. 135, 383–391 (2019).Article 
CAS 

Google Scholar 
Chen, F. et al. Effects of N addition and precipitation reduction on soil respiration and its components in a temperate forest. Agr. Forest. Meteorol. 271, 336–345 (2019).Article 

Google Scholar 
Lei, J. et al. Temporal changes in global soil respiration since 1987. Nat. Commun. 12, 1–9 (2021).
Google Scholar 
Wang, R. et al. Nitrogen application increases soil respiration but decreases temperature sensitivity: Combined effects of crop and soil properties in a semiarid agroecosystem. Geoderma 353, 320–330 (2019).Article 
CAS 

Google Scholar 
Du, K. et al. Influence of no-tillage and precipitation pulse on continuous soil respiration of summer maize affected by soil water in the North China Plain. Sci. Total Environ. 766, 144384. https://doi.org/10.1016/j.scitotenv.2020.144384 (2021).Article 
CAS 
PubMed 

Google Scholar 
Chen, X. & Chen, H. Y. Plant diversity loss reduces soil respiration across terrestrial ecosystems. Global Change Biol. 25, 1482–1492 (2019).Article 

Google Scholar 
Lang, A. K., Jevon, F. V., Ayres, M. P. & Matthes, J. H. Higher soil respiration rate beneath arbuscular mycorrhizal trees in a northern hardwood forest is driven by associated soil properties. Ecosystems 23, 1243–1253 (2020).Article 
CAS 

Google Scholar 
Huang, N. et al. Spatial and temporal variations in global soil respiration and their relationships with climate and land cover. Sci. Adv. 6, eabb8508 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Xiao, H. et al. The regulatory effects of biotic and abiotic factors on soil respiration under different land-use types. Ecol. Indic. 127, 107787. https://doi.org/10.1016/j.ecolind.2021.107787 (2021).Article 
CAS 

Google Scholar 
Liu, Y.-R. et al. New insights into the role of microbial community composition in driving soil respiration rates. Soil Biol. Biochem. 118, 35–41 (2018).Article 
CAS 

Google Scholar 
Wagg, C., Schlaeppi, K., Banerjee, S., Kuramae, E. E. & van der Heijden, M. G. Fungal-bacterial diversity and microbiome complexity predict ecosystem functioning. Nat. Commun. 10, 1–10 (2019).Article 
CAS 

Google Scholar 
Chen, L.-F. et al. Linkages between soil respiration and microbial communities following afforestation of alpine grasslands in the northeastern Tibetan Plateau. Appl. Soil Ecol. 161, 103882. https://doi.org/10.1016/j.apsoil.2021.103882 (2021).Article 

Google Scholar 
Choudhary, M. et al. Long-term effects of organic manure and inorganic fertilization on sustainability and chemical soil quality indicators of soybean-wheat cropping system in the Indian mid-Himalayas. Agr. Ecosyst. Environ. 257, 38–46 (2018).Article 

Google Scholar 
Zhang, M. et al. Increasing yield and N use efficiency with organic fertilizer in Chinese intensive rice cropping systems. Field Crop. Res. 227, 102–109 (2018).Article 

Google Scholar 
Bonanomi, G. et al. Repeated applications of organic amendments promote beneficial microbiota, improve soil fertility and increase crop yield. Appl. Soil Ecol. 156, 103714. https://doi.org/10.1016/j.apsoil.2020.103714 (2020).Article 

Google Scholar 
Gai, X. et al. Long-term benefits of combining chemical fertilizer and manure applications on crop yields and soil carbon and nitrogen stocks in North China Plain. Agr. Water Manage. 208, 384–392 (2018).Article 

Google Scholar 
Lai, R. et al. Manure fertilization increases soil respiration and creates a negative carbon budget in a Mediterranean maize (Zea mays L.)-based cropping system. Catena 151, 202–212 (2017).Article 
CAS 

Google Scholar 
Yan, T. et al. Negative effect of nitrogen addition on soil respiration dependent on stand age: Evidence from a 7-year field study of larch plantations in northern China. Agr. Forest Meteorol. 262, 24–33 (2018).Article 

Google Scholar 
Peng, Q. et al. Effects of nitrogen fertilization on soil respiration in temperate grassland in Inner Mongolia. China. Environ. Earth Sci. 62, 1163–1171 (2011).Article 
CAS 

Google Scholar 
Zeng, J. et al. Nitrogen fertilization directly affects soil bacterial diversity and indirectly affects bacterial community composition. Soil Biol. Biochem. 92, 41–49 (2016).Article 
CAS 

Google Scholar 
Levine, U. Y., Teal, T. K., Robertson, G. P. & Schmidt, T. M. Agriculture’s impact on microbial diversity and associated fluxes of carbon dioxide and methane. ISME J. 5, 1683–1691 (2011).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Chen, Q., Liu, Z., Zhou, J., Xu, X. & Zhu, Y. Long-term straw mulching with nitrogen fertilization increases nutrient and microbial determinants of soil quality in a maize–wheat rotation on China’s Loess Plateau. Sci. Total. Environ. 775, 145930. https://doi.org/10.1016/j.scitotenv.2021.145930 (2021).Article 
CAS 

Google Scholar 
Wang, J. et al. The impact of fertilizer amendments on soil autotrophic bacteria and carbon emissions in maize field on the semiarid Loess Plateau. Front. Microbiol. https://doi.org/10.3389/fmicb.2021.664120 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Subke, J. A., Inglima, I. & Francesca Cotrufo, M. Trends and methodological impacts in soil CO2 efflux partitioning: a metaanalytical review. Global Change Biol. 12, 921–943 (2006).Article 

Google Scholar 
Yan, W., Zhong, Y., Liu, J. & Shangguan, Z. Response of soil respiration to nitrogen fertilization: Evidence from a 6-year field study of croplands. Geoderma 384, 114829. https://doi.org/10.1016/j.geoderma.2020.114829 (2021).Article 
CAS 

Google Scholar 
Lamptey, S., Xie, J., Li, L., Coulter, J. A. & Jagadabhi, P. S. Influence of organic amendment on soil respiration and maize productivity in a semi-arid environment. Agronomy 9, 611. https://doi.org/10.3390/agronomy9100611 (2019).Article 
CAS 

Google Scholar 
Chen, Z. et al. Nitrogen fertilization stimulated soil heterotrophic but not autotrophic respiration in cropland soils: A greater role of organic over inorganic fertilizer. Soil Biol. Biochem. 116, 253–264 (2018).Article 
CAS 

Google Scholar 
Zheng, J., Zhang, X., Li, L., Zhang, P. & Pan, G. Effect of long-term fertilization on C mineralization and production of CH4 and CO2 under anaerobic incubation from bulk samples and particle size fractions of a typical paddy soil. Agr. Ecosyst. Environ. 120, 129–138 (2007).Article 
CAS 

Google Scholar 
Shen, J., Zhang, L., Guo, J., Ray, J. & He, J. Impact of long-term fertilization practices on the abundance and composition of soil bacterial communities in Northeast China. Appl. Soil Ecol. 46, 119–124 (2010).Article 

Google Scholar 
Chen, Q., An, X., Zheng, B., Ma, Y. & Su, J. Long-term organic fertilization increased antibiotic resistome in phyllosphere of maize. Sci. Total. Environ. 645, 1230–1237 (2018).Article 
CAS 
PubMed 

Google Scholar 
Zhang, W., Yu, C., Wang, X. & Hai, L. Increased abundance of nitrogen transforming bacteria by higher C/N ratio reduces the total losses of N and C in chicken manure and corn stover mix composting. Bioresource Technol. 297, 122410. https://doi.org/10.1016/j.biortech.2019.122410 (2020).Article 
CAS 

Google Scholar 
Chen, X. et al. Microbial carbon use efficiency, biomass turnover, and necromass accumulation in paddy soil depending on fertilization. Agr. Ecosyst. Environ. 292, 106816. https://doi.org/10.1016/j.agee.2020.106816 (2020).Article 
CAS 

Google Scholar 
Wang, J. et al. Nitrogen application increases soil microbial carbon fixation and maize productivity on the semiarid Loess Plateau. Plant Soil https://doi.org/10.1007/s11104-022-05457-7 (2022).Article 

Google Scholar 
Li, J. et al. The more straw we deep-bury, the more soil TOC will be accumulated: When soil bacteria abundance keeps growing. J. Soil Sediment 22, 162–171 (2022).Article 

Google Scholar 
Siczek, A., Frąc, M., Wielbo, J. & Kidaj, D. Benefits of flavonoids and straw mulch application on soil microbial activity in pea rhizosphere. Int. J. Environ. Sci. Te. 15, 755–764 (2018).Article 
CAS 

Google Scholar 
Zhao, S. et al. Change in straw decomposition rate and soil microbial community composition after straw addition in different long-term fertilization soils. Appl. Soil Ecol. 138, 123–133 (2019).Article 

Google Scholar 
Zhang, S. et al. Cow manure application effectively regulates the soil bacterial community in tea plantation. BMC Microbiol. 20, 1–11 (2020).Article 

Google Scholar 
Jiang, Y. et al. Crop rotations alter bacterial and fungal diversity in paddy soils across East Asia. Soil Biol. Biochem. 95, 250–261 (2016).Article 
CAS 

Google Scholar 
Drenovsky, R., Vo, D., Graham, K. & Scow, K. Soil water content and organic carbon availability are major determinants of soil microbial community composition. Microb. Ecol. 48, 424–430 (2004).Article 
CAS 
PubMed 

Google Scholar 
Rath, K. M., Fierer, N., Murphy, D. V. & Rousk, J. Linking bacterial community composition to soil salinity along environmental gradients. ISME J. 13, 836–846 (2019).Article 
CAS 
PubMed 

Google Scholar 
Zhao, F. et al. Changes of the organic carbon content and stability of soil aggregates affected by soil bacterial community after afforestation. CATENA 171, 622–631 (2018).Article 
CAS 

Google Scholar 
Goldfarb, K. C. et al. Differential growth responses of soil bacterial taxa to carbon substrates of varying chemical recalcitrance. Front. Microbiol. 2, 94. https://doi.org/10.3389/fmicb.2011.00094 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Zhao, J. et al. Response of soil microbial community to vegetation reconstruction modes in mining areas of the Loess Plateau, China. Front. Microbiol. https://doi.org/10.3389/fmicb.2021.714967 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Zhang, Y. et al. Fertilization shapes bacterial community structure by alteration of soil pH. Front. Microbiol. 8, 1325. https://doi.org/10.3389/fmicb.2017.01325 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Wang, X. et al. Organic amendments drive shifts in microbial community structure and keystone taxa which increase C mineralization across aggregate size classes. Soil Biol. Biochem. 153, 108062. https://doi.org/10.1016/j.soilbio.2020.108062 (2021).Article 
CAS 

Google Scholar 
Lin, Y. et al. Long-term manure application increases soil organic matter and aggregation, and alters microbial community structure and keystone taxa. Soil Biol. Biochem. 134, 187–196 (2019).Article 
CAS 

Google Scholar 
Woyke, T. et al. Symbiosis insights through metagenomic analysis of a microbial consortium. Nature 443, 950–955 (2006).Article 
CAS 
PubMed 

Google Scholar 
Zhang, B., Zhang, J., Liu, Y., Shi, P. & Wei, G. Co-occurrence patterns of soybean rhizosphere microbiome at a continental scale. Soil Biol. Biochem. 118, 178–186 (2018).Article 
CAS 

Google Scholar 
Wiens, J. J. et al. Niche conservatism as an emerging principle in ecology and conservation biology. Ecol. Lett. 13, 1310–1324 (2010).Article 
PubMed 

Google Scholar 
Deng, Y. et al. Molecular ecological network analyses. BMC Bioinf. 13, 1–20 (2012).Article 

Google Scholar 
Liao, H. et al. Complexity of bacterial and fungal network increases with soil aggregate size in an agricultural Inceptisol. Appl. Soil Ecol. 154, 103640. https://doi.org/10.1016/j.apsoil.2020.103640 (2020).Article 

Google Scholar 
Herren, C. M. & McMahon, K. D. Keystone taxa predict compositional change in microbial communities. Environ. Microbiol. 20, 2207–2217 (2018).Article 
PubMed 

Google Scholar 
Zhang, C., Jiao, S., Shu, D. & Wei, G. Inter-phylum negative interactions affect soil bacterial community dynamics and functions during soybean development under long-term nitrogen fertilization. Stress Biol. 1, 1–13 (2021).Article 
CAS 

Google Scholar 
Su, Y. G., Huang, G., Lin, Y. J. & Zhang, Y. M. No synergistic effects of water and nitrogen addition on soil microbial communities and soil respiration in a temperate desert. CATENA 142, 126–133 (2016).Article 
CAS 

Google Scholar 
Yang, C. et al. Assessing the effect of soil salinization on soil microbial respiration and diversities under incubation conditions. Appl. Soil Ecol. 155, 103671. https://doi.org/10.1016/j.apsoil.2020.103671 (2020).Article 

Google Scholar 
Banerjee, S. et al. Network analysis reveals functional redundancy and keystone taxa amongst bacterial and fungal communities during organic matter decomposition in an arable soil. Soil Biol. Biochem. 97, 188–198 (2016).Article 
CAS 

Google Scholar 
Chen, L.-F., He, Z.-B., Zhao, W.-Z., Kong, J.-Q. & Gao, Y. Empirical evidence for microbial regulation of soil respiration in alpine forests. Ecol. Indic. 126, 107710. https://doi.org/10.1016/j.ecolind.2021.107710 (2021).Article 
CAS 

Google Scholar 
Liu, S. et al. Decoupled diversity patterns in bacteria and fungi across continental forest ecosystems. Soil Biol. Biochem. 144, 107763. https://doi.org/10.1016/j.soilbio.2020.107763 (2020).Article 
CAS 

Google Scholar 
Lynch, M. D. & Neufeld, J. D. Ecology and exploration of the rare biosphere. Nat. Rev. Microbiol. 13, 217–229 (2015).Article 
CAS 
PubMed 

Google Scholar 
Chen, L. et al. Competitive interaction with keystone taxa induced negative priming under biochar amendments. Microbiome 7, 1–18 (2019).
Google Scholar 
Chiba, A. et al. Soil bacterial diversity is positively correlated with decomposition rates during early phases of maize litter decomposition. Microorganisms 9, 357 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Li, S., Wang, S., Fan, M., Wu, Y. & Shangguan, Z. Interactions between biochar and nitrogen impact soil carbon mineralization and the microbial community. Soil Till. Res. 196, 104437. https://doi.org/10.1016/j.still.2019.104437 (2020).Article 

Google Scholar 
Bao, S. Soil agrochemical analysis 30 (China Agricultural Press, Beijing, Chinese, 2000).
Google Scholar 
Zhai, L., Liu, H., Zhang, J., Huang, J. & Wang, B. Long-term application of organic manure and mineral fertilizer on N2O and CO2 emissions in a red soil from cultivated maize-wheat rotation in China. Agr. Sci. China 10, 1748–1757 (2011).Article 

Google Scholar 
Xia, W. et al. Autotrophic growth of nitrifying community in an agricultural soil. ISME J. 5, 1226–1236 (2011).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Pruesse, E. et al. SILVA: A comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids. Res. 35, 7188–7196 (2007).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B. 57, 289–300 (1995).MathSciNet 
MATH 

Google Scholar 
Layeghifard, M., Hwang, D. M. & Guttman, D. S. Disentangling interactions in the microbiome: A network perspective. Trends Microbiol. 25, 217–228 (2017).Article 
CAS 
PubMed 

Google Scholar 
Liaw, A. & Wiener, M. Classification and regression by randomForest. R news 2, 18–22 (2002).
Google Scholar 
Archer, E. rfPermute: Estimate permutation p-values for random Forest importance metrics. R package version 2(1), 81 (2020).MathSciNet 

Google Scholar 
Hooper, D., Coughlan, J. & Mullen, M. Structural equation modelling: Guidelines for determining model fit. Electron. J. Bus. Res. Methods 6(1), 53–60 (2008).
Google Scholar  More

Data sourcesWe analyzed records of confirmed and suspected livestock deaths attributed to anthrax occurring from 1 January 2006 to 31 December 2020 across Kenya (available online along with full code for the analysis in this paper https://github.com/spatialmodels/Kenyan_anthrax_model). The case records covering the entire country were reported from the Kenya Directorate of Veterinary Services (KDVS) located in Nairobi and the five Regional Veterinary Investigation Laboratories located in Nakuru, Eldoret, Karatina, Kericho, and Mariakani. The anthrax outbreaks were considered as any livestock (cattle, goats, sheep, pigs, camels) or wildlife deaths confirmed through clinical and laboratory diagnosis. Clinical diagnosis was defined as an acute disease accompanied by sudden death, bleeding from body orifices, swelling, lack of rigor mortis, and oedema of the neck and face in pigs. Laboratory confirmation was done through methylene blue staining to identify the characteristic bacterial capsule and the rod-shaped bacilli in clinical specimens collected from the infected carcasses.We extracted data from old paper records of livestock anthrax cases into Microsoft Excel. These records comprised the location of the livestock outbreaks, name of the farmer, number of animals dead and herd size, species affected, date, method of diagnosis, and the details of the reporting veterinary doctor. Since the locations of livestock anthrax outbreaks were reported at sub-county/district levels (districts refer to the old naming given to current sub-counties before the rollout of the current constitution), we recorded the geographic coordinates of livestock cases at the district level. During data cleaning, we removed duplicate coordinates, outliers, and entries with missing variables. In the end, we had 540 livestock cases that we used for analysis. The spatial granularity and prolonged surveillance period of these data allow for a more detailed perspective on the major drivers of anthrax across Kenya. We also collected wildlife data from the Kenya Wildlife Service (KWS). Most of the data from KWS was lacking information on the geographic coordinates of the outbreaks, so we visited the actual locations and collected the coordinates. We recorded 20 wildlife cases that we used to validate the performance of the spatial model.Processing socio-economic and ecological covariatesWe gathered geospatial data on ecological and socio-economic correlates of B. anthracis ecology and distribution. For the spatial model, we obtained the following variables: rainfall, vegetation, elevation, distance to permanent water bodies, and soil patterns. For the spatiotemporal models, we used human population estimates (total population, population density, and male and female population per sub-county), host population (livestock producing households, total number of indigenous, exotic dairy, and exotic beef cattle per sub-county), and agricultural practices that lead to soil disturbance (agricultural area under cultivation, number of farming households, and crop-producing households).We chose seven environmental covariates for the spatial model based on known correlates of B. anthracis suitability identified from previous peer-reviewed studies9,10,13,15,21,22,23. These comprised three soil variables, including soil pH (× 10) in H2O at a depth of 0 cm, exchangeable calcium at a depth of 0–20 cm, and soil water availability (volume of water per unit volume of soil) retrieved at a resolution of 250 m from the International Soil Reference and Information Centre (ISRIC) data hub (https://data.isric.org/geonetwork/srv/eng/catalog.search#/home). We used the shallowest depth available because although the bacterial spores can persist in the surface soil for up to five years and indefinitely in much deeper soils24, the spores in the surface soils are more likely to trigger host infection25. We retrieved monthly Enhanced Vegetation Index (EVI) data from 1 January 2006 to 31 December 2020 (180 tiles in total) from The Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MYD13A3 v.6) at a resolution of 1 km2 (https://lpdaac.usgs.gov/products/myd13a3v006/). The mean EVI was then calculated using QGIS by averaging all 180 tiles. EVI reduces variations in the canopy background and retains precision over dense vegetation conditions. Monthly Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) rainfall data from rain gauge and satellite observations was retrieved from the United States Geological Service (USGS) at a resolution of 0.05 degrees (https://climateserv.servirglobal.net/map). Since the rainfall data also comprised 180 tiles, the mean rainfall was calculated by averaging all 180 tiles using QGIS. We also collected data on the distance to permanent water bodies from a global hydrology map obtained from ArcGIS version 10.6.1.26 and elevation data at 1 km2 resolution from the Global Multi-resolution Terrain Elevation Data (GMTED2010) dataset available from USGS (Table 1).Table 1 Summary of the environmental variables used in the spatial model including variable name, unit, and spatial resolution.Full size tableFor the spatiotemporal sub-county-based models, we accessed the population data per sub-county (total population, male population, female population, and population density) from the 2019 Kenyan census report provided via the Humanitarian Data Exchange platform (https://data.humdata.org/dataset/kenya-population-per-county-from-census-report-2019). We also obtained data on livestock population (numbers of exotic dairy and beef cattle, and indigenous cattle), area of agricultural land in hectares, number of farming households, and the number of households actively practicing agriculture (crop production and livestock production) aggregated to the sub-county level from the 2019 Kenya Population and Housing Census volume IV provided by the OpenAfrica platform (https://open.africa/dataset/2019-kenya-population-and-housing-census).We conducted data exploration to check for outliers, collinearity, and the relationships between the covariates and the response variables. We used Cleveland dot plots to check for outliers. We measured collinearity using variance inflation factors (VIF), Pearson correlation coefficients, and pairs plots. For VIF scores, the covariates with scores higher than 3 were eliminated one-by-one until all the scores were equal to or less than 3. All the covariates included in the study had correlation coefficient values of less than 0.6 (Figs. 1, 2).Figure 1Results of correlation between covariates using Pearson’s correlation coefficient test for the spatial model. Correlation between covariates is shown by red numbers (negative correlation) and blue numbers (positive correlation). Correlations with a p-value  > 0.01 are regarded as insignificant and the correlation coefficient values are left blank. The figure was generated using R software v. 4.1.028.Full size imageFigure 2Results of correlation between covariates using Pearson’s correlation coefficient test for the spatiotemporal model. Correlation between covariates is shown by red numbers (negative correlation) and blue numbers (positive correlation). Correlations with a p-value  > 0.01 are regarded as insignificant and the correlation coefficient values are left blank. The figure was generated using R software v. 4.1.028.Full size imageSpatial model analysisWe used R version 4.1.0 together with the packages raster version 4.1.127, and R-INLA version 4.1.128 to conduct the data processing and statistical modelling. The R-INLA package applies the INLA framework in designing models. We used Quantum Global Information System (QGIS) version 3.16 (https://qgis.org) to create a 50 km buffer polygon around all the observed livestock outbreak points. We then created a 20 km2 grid within this buffer and counted the number of points within each grid cell to create a regular lattice with a given number of counts per cell. We then extracted the coordinates of the centroids of each cell to create marked locations with a given number of livestock cases per location. We essentially converted the data into a count process (number of livestock outbreaks per location). We had 95 cells with one or more counts which formed our new presence locations. We then randomly selected 95 pseudoabsences within the 50 km buffer polygon but at a distance of 10 km from the presence locations as shown in Fig. 3.Figure 3Spatial distribution of thinned livestock anthrax case locations across Kenya from 2006 to 2020. The map shows livestock anthrax case locations (n = 540) thinned to pixels of 20 km2 to form 95 new marked locations. The orange dots show the new presence locations which are marked points with colour intensity representing the number of livestock cases per location. The white triangles show the random pseudo-absence locations. The yellow squares are the wildlife cases obtained from the Kenya Wildlife Service. The green polygon is the background calibration buffer used to derive the random pseudo-absence locations. This map was generated using Quantum Geographical Information Systems (QGIS) v. 3.16.11 (https://www.qgis.org/en/site/forusers/download.html).Full size imageWe defined a Zero-inflated Poisson (ZIP) regression model with spatially correlated random effects, implemented as a generalized additive model (GAM) with anthrax incidence as the response variable. The model is defined as shown in Eqs. (1), (2), and (3)$${C}_{i} sim zero-inflated, Poisson left({mu }_{i},{p}_{i}right),$$
(1)
$$expectedleft({C}_{i}right)=left(1- {p}_{i}right)times {mu }_{i},$$
(2)
$$mathrm{log}left({mu }_{i}right)= alpha + sum_{j}{beta }_{j}{X}_{j,i}+ sum_{k}{delta }_{k,i}+{u}_{i},$$
(3)
where (Ci) denotes the observed number of anthrax livestock cases at location i, ({mu }_{i}) and ({p}_{i}) are parameters of the ZIP distribution. (expectedleft({C}_{i}right)) refers to the expected number of outbreaks at location i, (alpha) is the intercept, (beta) are the beta coefficients for the covariates, X is the matrix with all the covariates, (delta k) are the non-linear effects (cubic regression splines), and ({u}_{i}) is the spatial random effect at location i.To test whether the addition of the GAM smoothers and the spatially correlated random effects improved the fit of the model, we also considered candidate models without smoothers and spatial random effects. We tested three versions of the spatial model: the first used distance to water, elevation, and EVI as linear covariates without spatial random effects, the second applied non-linear terms to elevation and EVI also without spatial random effects, and the final model was similar to the second model but with the addition of spatial random effects. We then measured the DIC values of the candidate models to select the final spatial model.We conducted model validation by assessing the posterior distributions of the parameters and the residuals for adherence to the distributional assumptions. We checked whether the residuals were independent and normally distributed. We also plotted a sample variogram to check for any residual spatial auto-correlation using a well-defined method29. We then ran 1000 simulations to check whether the model was capable of handling zeros.The estimated model was used to map posterior predicted distributions for the incidence of anthrax disease (plotted as mean and 95% credible intervals). We validated the model using independent evaluation data withheld from the model calibration. This evaluation dataset comprises the wildlife cases collected from KWS. We then calculated the sensitivity by estimating the proportion of wildlife case locations correctly identified by the model, using the minimum presence training threshold (minimum value of the fitted presence training points).Spatiotemporal model analysisOur second objective was to investigate the socio-economic, population-based drivers of livestock anthrax risk at the sub-county level. These socioeconomic variables are usually collected at the sub-county level. Therefore, we developed a second areal model with the number of observations per sub-county as the new response variable. The occurrence data, gathered by the Kenya Directorate for Veterinary Services (KDVS), consisted of monthly case reports of livestock anthrax cases collected by all 290 sub-counties across Kenya between January 2006 to December 2020. We analyzed the whole monthly case time series from the year 2006 to 2020 and mapped out the annual counts of confirmed and suspected livestock anthrax cases across Kenya at the sub-county level to analyse the spatial and temporal trends throughout the surveillance period. The sub-county shapefiles that were used for mapping and modelling were derived from Humanitarian Data Exchange version 1.57.16 under a Creative Commons Attribution for Intergovernmental Organisations license (https://data.humdata.org/dataset/ken-administrative-boundaries).Due to the sparsity of data, we aggregated the monthly case counts and modelled the quarterly occurrence and incidence of anthrax at the sub-county-level scale, including spatial and temporal effects, to determine the spatial socio-economic drivers of livestock anthrax disease risk across Kenya. We used R-INLA version 4.1.1 (26) to conduct the data processing and statistical modelling. We used quarterly case counts that were confirmed per sub-county across the 15 years of surveillance (2006–2020) as a measure of anthrax incidence. Due to the zero-inflated and over-dispersed nature of the distribution, which is difficult to fit incidence counts, we employed a two-stage modelling approach using the hurdle model distribution to separately model anthrax occurrence (presence or absence) across all sub-counties via logistic regression, and incidence counts using a zero-inflated Poisson distribution. We were then able separately to estimate the contributions of the various socio-ecological factors that drive disease occurrence (the presence or absence of anthrax) and total incidence counts.We model the quarterly anthrax occurrence (n = 290 sub-counties over 60 quarters; 17,400 observations) where ({Y}_{i,t}) refers to the binary presence (denoted as 1) or absence (denoted as 0) of anthrax in sub-county i during year t, and ({P}_{i,t}) is the probability of anthrax occurrence, thus:$${Y}_{i,t} sim Bernoullileft({P}_{i,t}right).$$
(4)
We model quarterly anthrax incidence counts ({C}_{i,t}) using a zero-inflated Poisson process with parameters ({mu }_{i,t}) and ({p}_{i,t}) (see Eq. (5)). Equation (6) denotes the expected values for the ZIP distribution at sub-county i during year t.$${C}_{i,t} sim Zero-inflated, Poisson left({mu }_{i,t},{p}_{i,t}right),$$
(5)
$$expectedleft({C}_{i,t}right)=left(1- {p}_{i,t}right)times {mu }_{i,t}.$$
(6)
Both the Bernoulli and the ZIP distributions are modelled separately as functions of the covariates and the spatial and temporal random effects using a general linear predictor as shown in Eqs. (7) and (8):$$logit left({P}_{i,t}right)= alpha + sum_{j}{beta }_{j}{X}_{j,i}+{u}_{i,t}+{v}_{i,t}+{y}_{i,t},$$
(7)
$$mathrm{log}left({mu }_{i,t}right)= alpha + sum_{j}{beta }_{j}{X}_{j,i}+{u}_{i,t}+{v}_{i,t}+{y}_{i,t},$$
(8)
$${y}_{i,t}= {y}_{i,t-1}+ {w}_{i,t},$$
(9)
where α denotes the intercept; (X) signifies a matrix made up of the socio-economic covariates accompanied by their linear coefficients denoted as (beta); spatiotemporal reporting trends at the sub-county level were accounted for in the models using spatially structured (({u}_{i,t}); conditional autoregressive) and unstructured noise (({v}_{i,t}); i.i.d—independent and identically distributed) random-effects specified jointly as a Besag–York–Mollie model30,31, as well as temporally structured (({y}_{i,t})) random effects of the first order where ({w}_{i,t}) is a pure noise term that is normally distribute with a mean of zero and a variance of σ2. We used uninformative priors with a Gaussian distribution for the fixed effects and penalized complexity priors for the hyperparameters of all the random effects.For the two spatiotemporal models, we applied linear effects for all the variables: population density, total population, number of exotic dairy cattle, agricultural land area, and number of livestock producing households. We scaled the continuous covariates by standardizing them (to a mean of 0 and standard deviation of 1) before fitting the linear fixed effects.We used R-INLA to conduct model inference and selection and used DIC to evaluate the model fit for both the occurrence and incidence models. For both models (occurrence and incidence), we created 4 candidate models, compared them, and selected the model with the lowest DIC as the final model. The candidate models included: a baseline intercept only model; a second model with the intercept and covariates; a third model with the intercept, covariates, and the spatial random effects; and a fourth model with the intercept, covariates, spatial random effects, and a temporal trend.We evaluated the posterior distributions of the parameters and the residuals for adherence to the distributional assumptions. We assessed the residuals to check whether they were independent and normally distributed. We also plotted the residuals against the covariates to check for any non-linear patterns using a well-defined method29. We then ran 1000 simulations to check whether the model was capable of handling zeros.Ethics statementLicence to conduct the research was granted by the National Council for Science, Technology, and Innovation (NACOSTI) under reference number 651983, and the Kenya Wildlife Service under reference number KWS-0003-01-21. More

We used a hierarchical Bayesian modeling approach to evaluate the relationship between the spring phenology of tree and wildflower species and various climate drivers (see Methods). Following model selection, our final model structure included fixed effects of average spring (March–April) temperature and elevation, as well as species-level random effects. We show continental distributions of spring temperature values in Fig. 1b (means and standard deviations are listed in Table S2). We report estimates for spring temperature sensitivities from the final model structure in the main text. Parameter estimates for elevation sensitivities as well as the model performance of other potential drivers and combinations of drivers are reported in Tables S3 and S4. An extended discussion of model assumptions and limitations is included in the Supplementary Information.Sensitivity differences by strataTree leaf out phenology (LOD) was substantially more sensitive to average spring temperature in North America (mean = −3.62 days °C−1; 95% credible interval (CI) = [−3.76, −3.49]) than in Europe (mean = −2.79; CI = [−3.27, −2.30]) and Asia (mean = −2.62; CI = [−2.97, −2.26]; Fig. 2). These values are consistent with previously reported phenological sensitivities in North America7 (−5.5 to −3.3 days °C−1) and Europe8 (−4.1 to −3.0 days °C−1), as the credible intervals from our results overlap with the reported credible intervals of prior studies. However, the Asian LOD sensitivity was less sensitive than previously reported27 (−3.50 to −3.03 days °C−1), potentially owing to differences in species selection28 or model structure. Previously reported sensitivities were determined in separate studies using either observational data7,8 or long-term observation-based weather station data27. The general consistency between our findings suggests that phenology data from herbarium collections are good indicators of patterns in natural systems29,30,31, a point supported by a recent study of phenological sensitivity derived from herbaria and from observed citizen science data32. These herbarium-based results provide evidence that phenological sensitivity differs across the temperate forest biome (but see ref. 33 for evidence of differences in response to warming and chilling accumulation). To our knowledge, our study is the first to contrast overstory and understory phenology across multiple continents and, therefore, to find differences in phenological sensitivity between trees and forest wildflowers across continents. We recommend future studies explore these differences using alternative approaches and methodologies that focus on the physiological basis for and mechanisms that underlie these patterns.Fig. 2: Posterior estimated means and 95% credible intervals for spring temperature sensitivity.Shapes represent parameter estimates for wildflower First Flower Date (FFD, blue circles; n = 1418, 618, and 1060 for Asia, Europe, and North America, respectively) and canopy tree Leaf Out Date (LOD, yellow triangles; n = 899, 532, and 995, for Asia, Europe, and North America, respectively). Estimates are considered different from 0 if credible intervals do not overlap the dashed 0 line and are considered different from each other if credible intervals do not overlap.Full size imageIn contrast to trees, wildflower sensitivity to spring temperature was similar across all three continents and exhibited no strong differences (i.e., overlap in 95% Bayesian credible intervals) among continents (means and 95% credible intervals in brackets: North America = −3.14, [−3.28, −3.00]; Europe = −3.02, [−3.48, −2.56]; Asia = −3.12, [−3.36, −2.86]; Fig. 2). These values are also generally consistent with those reported elsewhere in the literature (i.e., 95% credible intervals overlap with those reported in other studies; −2.2, [−3.7, −0.76] days °C−1 in North America7 and −3.6, [−4.04, −3.18] days °C−1 in Europe9), although we are unaware of any studies that have estimated phenological sensitivity for Asian forest wildflowers in days °C−1. Ge et al.3 report herbaceous plant sensitivity of −5.71 days per decade in Asia (±7.90 standard deviation; based primarily on long-term observational data), which appears to be roughly consistent with our model results, but the difference in units makes this more speculative than the other comparisons. Discrepancies in mean responses between this study and others may be due in part to different types of data (herbarium specimens versus field observations) and to choice in focal taxa, as temperature sensitivity has been shown to vary widely across taxa28.Particularly noticeable in our results was that r2 coefficients of predicted versus observed phenology were much higher in North America (0.70 and 0.76 for wildflower and tree models, respectively) compared to Asian (0.40 and 0.44, respectively) and European models (0.41 and 0.25, respectively). This difference in model performance could be due to the higher interannual variability of spring temperatures in North America33, leading to greater selective pressure for strong sensitivity to spring temperatures in North American plants. This difference could explain why North American species exhibit higher correlation of phenology with average spring temperatures (Table S4). Alternatively, European and Asian species may have stronger phenological responses to alternative spring forcing windows, winter chilling temperatures, or photoperiod, relative to the March–April temperature period used in this study (see Methods). We think the latter explanation is unlikely, given the strong correlations of phenology with spring temperature across all continents (see Supplementary Information – Justification for March–April Temperature Window).Herbarium-based phenological models may be improved by accounting for spatial autocorrelation within the dataset. For example, Willems et al.9 found that including spatial autocorrelation significantly improved predictability of European herbaceous flowering phenology, even when accounting for multiple drivers of spring phenology. We followed a similar approach as their study and found similar improvements in model performance with the addition of spatial autocorrelation (Tables S3–S4) that had substantial positive effects on r2 values of Asian and European models. However, spatial distributions of specimens differed substantially among continents (see Figs. S2–S4), and these differences could lead to artifacts that make results unreliable to interpret (see Supplementary Information). Therefore, we focus here on results for models without spatial autocorrelation while acknowledging that spatial aggregation of herbarium specimens in Europe and Asia may be partially responsible for the relatively lower r2 values. We encourage other researchers to explore this question further both with our data set and other datasets.Climate change and spring light windowsThe relative difference between wildflower and tree sensitivity varied substantially among continents, with wildflowers being approximately equally as sensitive to spring temperature as trees in Asia and Europe but substantially less sensitive (i.e., 95% BCI do not overlap) than trees in North America (Fig. 2). Importantly, these differences were driven by changes in tree phenological sensitivities among continents and resulted in different expectations for spring light window duration (i.e., the difference in time between estimated wildflower flowering date and canopy tree leaf out date) on different continents under current climate conditions (Fig. 3), based on modeled leaf out and flowering under a climate scenario derived from average climate conditions from 2009–2018 (Fig. S5).Fig. 3: Current estimated phenological escape duration in northern temperate deciduous forests.Estimated mean difference between wildflower First Flower Date (FFD) and canopy tree Leaf Out Date (LOD) (in days) under current climate conditions (averaged from 2009–2018, see methods) in a Asia, b Europe, and c North America. Negative values indicate tree LOD is estimated to occur before wildflower FFD. Estimations were cropped by the estimated area of broadleaf and mixed-broadleaf forest (see methods). Dark gray regions indicate areas where the consensus land classification is More

Root morphological traitsThe wild-type parent ISR42-8 produced longer root length (RL) than the modern cultivar parent Golf and tested lines (Table S2, Fig. 1A.h,A.f [h = hydroponic; f = field]). The tested lines of the two evolving barley populations displayed significant variations under hydroponic conditions. Barley lines evolved under OCS had on average 3484 mm longer roots compared to CCS under hydroponic treatment (Fig. 1A.h, Table 2). Complementary results under field conditions show as well higher RL for the OCS lines, even though the variance was significantly less pronounced (Fig. 1A.f). In addition, a less evident variance was observed in the field within both groups compared to the hydroponic (Fig. 1A.d). Across both experimental setups, the observed range of RL was higher in the OCS lines [Standard deviation (SD)OCS = 883, SDCCS = 597] (Table 2).Figure 1Significantly variant root morphological phenotypes. Boxplots illustrate the overall distribution of observed data points for the parents Golf and ISR 42-8 as well as for the conventional (CCS) and organic (OCS) lines. Density plots highlight the overall distribution of organic and conventional adapted lines. (A)—Root length (RL)—the sum of all roots harvested in millimeters (mm), illustrated for all four groups. (A.h)—root length measured in the hydroponic experiment; (A.f)—field experiment; A.d—distribution histogram for root length in both field and hydroponic experiment for CCS and OCS adapted lines. (B)—the ratio of root length to volume (L/V). Data available for hydroponics (B.h), field (B.f), and distribution of the ratio of root length to volume illustrated in B.d. (C.f)—Root mass density (RMD) from the field; (D.f)—Root angle (RA) from the field, distribution of the root angle illustrated in (D.d); (E.f)—root tip per plant count from the field, corresponding histogram visualized in (E.d). (F.f)—root fork per plant count from the field.Full size imageTable 2 Comparison of organic and conventional population root phenotypes under field and hydroponic evaluation.Full size tableThe root length to volume (L/V) is an important indicator of the soil volume that can be explored by the roots. Under hydroponics conditions, variations were found for L/V between the parental genotypes as well as between the OCS and CCS populations (Tables 2 & S3). The organic lines were characterized by a significantly higher L/V, indicating a much more distinct exploration of the soil by these lines (Fig. 1B.h,B.f). In comparison to field, highest diversity in L/V was found under hydroponic experiments within both OCS and CCS populations (Fig. 1B.d).The root mass density (RMD) is the ratio of root volume for a given root mass and is a key indicator of root thickness. Although significant variations existed between ISR42-8 and Golf under hydroponics conditions, such significant variations were not found between the OCS and CCS groups (P = 0.09) (Fig. 1C.f and Tables 2 and S3).The root angle (RA) measurements were only performed under field conditions since plants grown under hydroponics conditions were placed in uniform growing vessels and the direction of root growth is restricted by tubes. Significant variation was observed for the RA between the two parental lines, which was also reflected in the CCS and OCS lines (Fig. 1D.f). ISR42-8 was characterized by an 11.5° average narrower RA than Golf (Table S3). The RA was 4.1° bigger in the OCS compared to the CCS population (P = 0.005) (Table 2). However, a higher diversity in RA was observed in the OCS compared to the CCS lines (Fig. 1D.d, Table 2).In addition to the RA, the number of root tips and forks was measured under field conditions only. Both tips and forks indicate a similar pattern, where the OCS lines produced on average more for both PForks = 0.014, PTips = 0.0041 (Fig. 1E.f,F.f). After applying a P-adjustment, the number of forks count was no longer significantly different between OCS and CCS (PForks = 0.07, Table 2). Complementary, ISR 42-8 was observed to produce more tips and forks than Golf, which remained highly significant even after probability adjustment (Fig. 1E.f,F.f, Table S3). The distribution and the standard deviation of observed phenotypes highlight once more the fact that the OCS lines tend to have a higher variation (Fig. 1E.d, Table 2). Similarly, a significant increasing trend was recorded in root surface area (RSA) and root average diameter (RAD) by ISR42-8 as compared to Golf under hydroponics (Table S3). Contrasting to the parental genotypes, no variation was observed between OCS and CCS lines for RSA (Table 2).Root anatomical traitsWithin the observed anatomical traits, four were considered due to their relevance and variation between the systems. In both hydroponic and field experiments, significant variations were observed for the late metaxylem number (LMXN) between the parental lines as well as OCS and CCS lines (Tables 2 & S3, Fig. 2A.h,A.f). An increased LMXN for ISR 42-8 compared to Golf was observed (Fig. 2A.h). Regarding the CCS and OCS lines, a heterogenic scenario was presented over both experimental setups. While the median LMXN under CCS was identical with ISR 42-8 in the seedling stages of plant development (Fig. 2A.h), it was much lower in flowering stages under field conditions (Fig. 1A.f). Additionally, the LMXN was significantly higher in the CCS lines in the seedling stage compared to OCS lines, vice-versa LMXN was observed at the flowering time point (Fig. 2A.d).Figure 2Significantly variant root anatomical traits. Boxplots illustrate the overall distribution of observed data points for the parents Golf and ISR 42-8 as well as for the conventional (CCS) and organic (OCS) lines. (A) –Late metaxylem number (LMN)—the sum of all roots harvested and expressed by plant−1, illustrated for all four groups. A.h—Late metaxylem number measured in the hydroponic experiment; A.f—field experiment; A.d—distribution histogram for late metaxylem number in both field and hydroponic experiment for CCS and OCS adapted lines. (B)—Aerenchyma area (AA). Data available for hydroponics (B.h), field (B.f), and distribution of the aerenchyma area illustrated in (B.d). (C.f)—Total cortical area (TCA) from the field; (D.f)—Root cross-section area (RA) from the field, distribution of the total cortical area and root cross-section area illustrated in C.d and D.d, respectively.Full size imageThe intercellular space, represented by the aerenchyma area (AA), was observed to be significantly more pronounced in the tested CCS compared to OCS lines in both environments (Fig. 2B.h,B.f). Furthermore, the OCS population did not show significant differences to both parents under hydroponics conditions, however, when grown under field conditions, it was noted that Golf had a significantly higher AA mean value as compared to the OCS population (Table S3). As illustrated by the values, the AA expended from early to late stages by a magnitude of 10-folds (Fig. 2B.d). In general, although the two parents did not indicate phenotypic variations, OCS and CCS lines showed significant variations (Table S4).A 0.12 mm2 decreased average total cortical area (TCA) was recorded in the OCS compared to the CCS population under field conditions (P = 0.003, Fig. 2C.f), although substantial variations for TCA was observed within OCS and CCS populations (Fig. 2C.d). The root cross-section area (RXA) is a two-dimensional axis of the root which is an important indicator of root thickness. In the hydroponic examination of the seedling stage, significant variations existed between the CCS and ISR42-8 as well as OCS population (Tables 2 and S3). The complementary study under field conditions observed a noticeable variation for OCS from both parental genotypes and the CCS (Table S3). About 0.13 mm2 increased average value for RXA was identified for CCS (Fig. 2D.f), while consistent significant variations were also observed between the populations in the under field experiment, where 0.13 mm2 increased average value for RXA was identified for CCS (Fig. 2D.f). Analog (Fig. 2D.d). Analogue to the AA, the RXA indicates a lower root extension in the OCS compared to the CCS population. For the stele area (SA), significant variations were only observed at the flowering stage, where ISR42-8 generally had the highest SA and varies significantly between Golf and its progeny lines (Tables 2 and S3).Shoot-related traitsBeyond the root-related phenotypic observations, above-ground characteristics were also recorded to assess the root-borne shoot dynamics (Figs. 3 and S2). Among the OCS and CCS populations and the parents, ISR42-8 had the longest duration of emergence. While CCS-adapted lines took on average 5.8 days of emergence (DE), OCS-adapted lines emerged 1.8 days later (7.6 days) (Fig. S2). No variation was observed for the tiller number (TN) throughout all tested groups, while ISR 42-8 tends to produce much more leaf number (LN), accompanied by a lower plant height (PH) and higher shoot dry weight (SDW) (Table S4, Fig. S2). The OCS and CCS plants significantly differed in PH as well as SDW (Fig. 3B,C). The LN was marginally above the probability threshold of 0.05 (p = 0.058, Fig. 3A), with a clear tendency of increased variability in phenotypic variation (Fig. 3D). Similar trend was recorded for the SDW (Fig. 3F).Figure 3Above-ground plant characteristics. Boxplots illustrate the overall distribution of observed data points for the parents Golf and ISR 42-8 as well as for the conventional (CCS) and organic (OCS) lines under the hydroponic experiment. (A)—Leaf number (LN) expressed by; (B)—Plant height (PH) and C-Shoot dry weight (SDW). The data distribution of the leaf number, plant height and shoot dry weight is illustrated in (D,E,F), respectively.Full size imageInterconnection of root-shoot traitsWe performed inter-trait correlation analysis to unravel association among root traits and in between root and shoot phenotypes (Fig. 4). Pearson correlation coefficient revealed significant correlations among root-shoot traits. LN, PH and SDW had strong positive associations with all root architectural traits under hydroponic conditions (P  0.30) in both CCS and OCS, while DE has negative association with all shoot traits (r = −0.17 to −0.48) (Fig. 4A.h,B.h). A consistent negative relationship was observed for L/V with shoot traits such as LN, PH and SDW and root morphological traits such as RL, RSA and RAD in both CCS and OSC populations (Fig. 4A.h,B.h). A strong negative association existed between RL and all shoot morphological, root architectural and anatomical traits in both populations, except for L/V where a weak negative (r = −0.09) association was displayed only in the OCS. Likewise, all above-ground traits and all root architectural traits exhibited significant positive associations with all root anatomical features in both groups with an exception for the AA (Fig. 4A.h,B.h). Moreover, correlation analysis revealed strong positive relationships in both groups of SDW and root dry weight (RDW) to all above-ground traits, below-ground traits including, RL, SA, and RAD, as well as in all root anatomical traits (Fig. 4A.h,B.h). This means that the growth of tissue and organ is proportional to the increase in total dry biomass. More importantly, we observed a significant positive correlation among most of the root morphological, architectural, and anatomical traits in both OCS and CCS adapted populations, with few exceptions such as L/V (Fig. 4A.h,B.h).Figure 4Correlation matrix for shoot morphological (only in hydroponic conditions; A.h and B.h), root architectural and anatomical traits in two groups of barley populations and their parental lines grown across two growing conditions. (A)—conventional and (B)—organic cropping systems. (A.h)—conventional under hydroponic, (B.h)—organic under hydroponic, (A.f)—conventional under field, and B.f—organic under field conditions. The color scale represents Spearman’s ranked correlation coefficient. A larger circle size indicates a smaller p-value; blank cells represent that correlation was non-significant at P  −0.90) and RDW (r =   > −0.80) (Fig. 4A.f,B.f) for CCS and OCS populations respectively, which means that narrower the angle of the nodal roots, the longer was the root system. The two root branching traits, the number of tips and number of roots forks which were known to be associated and dependent on the RL have a strong positive correlation reflected by r = 0.81 and 0.90 in CCS and r = 0.74 and 0.84 in OCS developed lines, respectively, while they have a significant negative correlation with RA (r = −0.72 to −0.77) in both barley groups. In addition, RA had also strong negative relationship to RDW contributing architectural traits including, RMD (r = −0.81 to −0.84) and L/V (r = −0.34 to −0.44). However, no positive associations were observed for RA and all root anatomical traits in both OCS and CCS populations (Fig. 4A.f,B.f).Allometry analysisThe correlation analysis identified interconnection among root and shoot-related traits. Therefore, we checked if these correlations can be explained by allometric relations (Tables 3 and 4).Table 3 Summary of allometric analysis of root-shoot system traits under hydroponic condition.Full size tableTable 4 Summary of allometric analysis of root-shoot system traits under field condition.Full size tableIn the hydroponic environment, we observed a total of ten allometric relations, from which six were annotated to the PH. The PH was allometrically related to the SDW, the RSA, the RV, the RDW, the SRL, and the RMD (Table 3). Besides, the SDW was allometrically associated with the RSD. Furthermore, the TCA was related to the RXA. Finally, an allometry relationship was detected between the LMXN and AA (Table 3).In the field experimental setup, we detected in total ten allometric relations (P  More

Data sourcesData for the trophic parameters and data categories listed in Tables 1 and 2 were gathered from peer-reviewed scientific publications, grey literature (e.g., agency reports, theses, and dissertations) and unpublished data by the authors of this paper. Data compilation on stomach contents, stable isotopes, FATM, and trophic positions, focussed on mesopelagic organisms, their potential prey and predators. For major and trace elements, energy density and estimates of diet proportions, our search concentrated on mesopelagic taxa. Nevertheless, we also gathered information from small or intermediate-sized epi-, bathy- or benthopelagic species found in the compiled data sources. These species were included because they play key roles in most marine ecosystems, both as important consumers of phytoplankton and zooplankton, and prey for many top predators, and can represent alternative energy pathways to mesopelagic organisms. However, we stress that the data coverage for these species in the current version of the database is very incomplete. Our main interest was on data from the central and eastern North Atlantic, and the Mediterranean, corresponding to the study regions of the SUMMER project. When we could not find suitable data within this region, we extended the geographic scope of our literature search to the western North Atlantic. We did not search for datasets in open access repositories since those data can be easily accessed and extracted. However, some of the data provided by the authors of this paper have been previously deposited in PANGAEA.DNA sequencing-based methods, such as metabarcoding and direct shotgun sequencing, are emerging as promising tools in dietary analyses due to the high resolution in taxonomic identification of many prey simultaneously, and the potential to provide quantitative diet estimates from relative read abundance29. However, recent studies have shown that various methodological and biological factors can break the correlation between the number and abundance of ingested prey and the prey DNA present in the sample, and lead to biased estimates of taxonomic diversity and composition of diet29,30. Given the uncertainties remaining in the interpretation of DNA sequencing-based diet data, we decided not to include these data in MesopTroph until additional research demonstrates that these techniques can be confidently applied for quantitative diet assessment.We identified available data sources in the literature through systematic searches on Web of Science, Google Scholar, ResearchGate, and the Google search engine. We used multiple combinations of terms related to specific data categories (Table 3), in conjunction with the common or scientific taxon names (from genus to order), and the ocean basin. For example, the search for stomach content data of fishes belonging to the family Myctophidae was undertaken using the following terms: “stomach content” OR “gut content” OR “prey composition” OR “diet composition”, AND “mesopelagic fish” OR “myctophid” OR “Myctophiformes” OR “Myctophidae”, AND “Atlantic” or “Mediterranean”. For the mesopelagic and predator species known to be numerically abundant in the SUMMER study regions, we performed a second literature search using the common or scientific name of the species, along with the terms “diet”, “feeding habits”, “trophic ecology”, “trophic markers”, or “food web”. We also examined the literature cited within each collected publication to locate additional data sources.Table 3 Terms used in the literature search for each data category.Full size tableWe next screened the full text of the compiled studies and retained data sources that: (1) were collected within the region of interest, (2) reported quantitative data for the trophic parameters of interest, (3) reported the number of samples for pooled or aggregated data, and (4) provided sufficient details on the methodology to enable a quality check. In the case of stable isotope data, we only included data sources reporting both δ13C and δ15N measurements.Data extraction, cleaning, and formattingWe created a template table for each data category in Microsoft Excel to assemble all datasets into a single file, and to facilitate cleaning and standardization of data records. We added a large number of metadata fields to the tables to annotate details about the sampling (e.g., location, date, methods), sampled specimen(s) (e.g., taxonomy, number and size of individuals, number of replicates, tissue analysed), and data source (e.g., full reference, DOI) for every record.Data contributors formatted and incorporated their datasets directly into the tables. For published sources, the data and associated metadata were extracted manually or digitized from the article text, tables, or supplementary material into the tables. Extraneous or hidden characters, and values such as “NA” (Not Available) or “ND” (Not Determined), were deleted from the parameter and metadata fields. Measurements of trophic parameters were standardized to the same units (see Tables 1 and 2). Parameter values that were clearly incorrect (e.g., δ15N  > 20, or the frequency of occurrence of a prey higher than the number of stomachs sampled) were corrected by searching for the value within the data source. When values could not be corrected, we deleted that data record.When available, we extracted information at the individual level. However, most studies reported data obtained from pooled samples of the same species. In some cases (e.g., small specimens such as planktonic organisms), a minimum and maximum number of individuals in the sample was provided instead of the actual number of individuals sampled. We added two columns to the tables presenting the minimum and maximum number of individuals in the sample. By filtering the column “Ind No (maximum per sample)” for values >1, users can easily identify records with aggregated data and differentiate them from records where information was drawn from a single individual (i.e., where “Ind No (maximum per sample)” =1). In addition, the tables Stomach contents and Estimates of diet proportions include a field “Sample ID” with a unique identifier of the sample. If data are reported at the individual level (i.e., “Ind No (maximum per sample)” =1) then Sample ID is the individual animal ID. If the data are from a group of individuals (i.e., “Ind No (maximum per sample)” >1), then Sample ID identifies that group.We standardized the taxonomic classification and nomenclature of fishes and elasmobranchs following the Eschmeyer’s Catalog of Fishes (http://researcharchive.calacademy.org/research/ichthyology/catalog/fishcatmain.asp)31,32. For the remaining taxa, we used the World Register of Marine Species (http://www.marinespecies.org/)33. Unaccepted or alternate taxon names were replaced by the most up-to-date valid name. When the identification of a taxon was uncertain, the taxonomic level of identification was decreased to a satisfactory level. For example, prey reported as “Cephalopods” were changed to “Cephalopoda”, “Sepiolids” to “Sepiolidae”, and “Myctophum punctatum?” to the genus “Myctophum”.Stomach contentsStomach contents analysis is a standard dietary assessment method that potentially enables quantifying diet components with high taxonomic resolution34. Three parameters are typically used to describe diet composition from stomach contents: the number of individuals of a prey type as a proportion of the total number of prey items (%N), the proportion of a prey item by weight or volume (%W), and the proportion of stomachs containing a particular prey item (i.e., percent frequency of occurrence, %F)35. When available, we collected data on the three parameters, as well as on the absolute number, weight, and frequency of occurrence of each prey type in the stomachs of each sampled individual or group of individuals. If stated in the data source, we indicate if prey weights were directly measured or reconstructed from hard remains (fish otoliths and vertebrae, cephalopod beaks), and if they represent dry or wet weight. Some datasets contained records of prey items without corresponding weights or numbers. As a result, the cumulative percent of all prey items did not sum to 100%. This occurred in 11 data records for the cumulative %W, and nine for the cumulative %N. While we checked the accuracy of percentage values and adjusted rounding errors, we did not attempt to fill in missing values nor did we remove records with missing values. When prey values were reported by an upper bound (e.g., “ More

MaterialsTrunk wash was collected from one male (Tembo, born 1985) and five female (Tonga, 1984; Numbi, 1992; Mongu, 2003; Iqhwa, 2013; Kibali, 2019) African elephants at the Vienna Zoo during routine procedures. Briefly, 100 mL of a sterile 0.9% saline solution is injected in each nostril of the trunk, which is kept in a lifted position, so that the solution is running up to the base of the trunk. The mixture of the solution and trunk mucus is collected in sterile plastic bags by active blowing of the elephant. Chemicals were all from Merck, Austria, unless otherwise stated. Restriction enzymes and kits for DNA extraction and purification were from New England Biolabs, USA. Oligonucleotides and synthetic genes were custom synthesised at Eurofins Genomics, Germany.Ethics declarationWe confirm that the trunk wash performed to provide a sample of the mucus was carried out as a routine procedure to monitor the health of elephants at the Vienna Zoo and in accordance with relevant guidelines and regulations.Trunk wash fractionationTrunk wash was centrifuged for 1 h at 10,000 g, the supernatant was dialyzed against 50 mM Tris–HCl buffer, pH 7.4 and concentrated by ultrafiltration in the Amicon stirred cell, then fractionated by anion-exchange chromatography on HiPrep-Q 16/10 column, 20 mL (Bio-Rad), along with standard protocols.Protein alkylation and digestion, and mass spectrometry analysisSDS-PAGE gel portions of proteins from whole elephant trunk wash (for component identification), chromatographic fractions of the elephant trunk wash (for PTMs analysis) or SDS-PAGE gel bands of LafrOBP1 expressed in P. pastoris were in parallel triturated, washed with water, in gel-reduced, S-alkylated, and digested with trypsin (Sigma, sequencing grade). Resulting peptide mixtures were desalted by μZip-TipC18 (Millipore) using 50% (v/v) acetonitrile, 5% (v/v) formic acid as eluent, vacuum-dried by SpeedVac (Thermo Fisher Scientific, USA), and then dissolved in 20 μL of aqueous 0.1% (v/v) formic acid for subsequent MS analyses by means of a nanoLC-ESI-Q-Orbitrap-MS/MS system, comprising an UltiMate 3000 HPLC RSLC nano-chromatographer (Thermo Fisher Scientific) interfaced with a Q-ExactivePlus mass spectrometer (Thermo Fisher Scientific) mounting a nano-Spray ion source (Thermo Fisher Scientific). Chromatographic separations were obtained on an Acclaim PepMap RSLC C18 column (150 mm × 75 μm ID; 2 μm particle size; 100 Å pore size, Thermo Fisher Scientific), eluting the peptide mixtures with a gradient of solvent B (19.92/80/0.08 v/v/v water/acetonitrile/formic acid) in solvent A (99.9/0.1 v/v water/formic acid), at a flow rate of 300 nL/min. In particular, solvent B started at 3%, increased linearly to 40% in 45 min, then achieved 80% in 5 min, remaining at this percentage for 4 min, and finally returned to 3% in 1 min. The mass spectrometer operated in data-dependent mode in positive polarity, carrying out a full MS1 scan in the range m/z 345–1350, at a nominal resolution of 70,000, followed by MS/MS scans of the 10 most abundant ions in high energy collisional dissociation (HCD) mode. Tandem mass spectra were acquired in a dynamic m/z range, with a nominal resolution of 17,500, a normalized collision energy of 28%, an automatic gain control target of 50,000, a maximum ion injection time of 110 ms, and an isolation window of 1.2 m/z. Dynamic exclusion was set to 20 s36.Bioinformatics for peptide identification and post-translational modification assignmentRaw mass data files were searched by Proteome Discoverer v. 2.4 package (Thermo Fisher Scientific), running the search engine Mascot v. 2.6.1 (Matrix Science, UK), Byonic™ v. 2.6.46 (Protein Metrics, USA) and Peaks Studio 8.0 (BSI, Waterloo, Ontario, Canada) software, both for peptide assignment/protein identification and for post-translational modification analysis.In the first case, analyses were carried out against a customized database containing protein sequences downloaded from NCBI (https://www.ncbi.nlm.nih.gov/) for superorder Afrotheria (consisting of 192,838 protein sequences, December 2021) plus the most common protein contaminants and trypsin. Parameters for database searching were fixed carbamidomethylation at Cys, and variable oxidation at Met, deamidation at Asn/Gln, and pyroglutamate formation at Gln. Mass tolerance was set to ± 10 ppm for precursors and to ± 0.05 Da for MS/MS fragments. Proteolytic enzyme and maximum number of missed cleavages were set to trypsin and 3, respectively. All other parameters were kept at default values. In the latter case, raw mass data were analyzed against a customized database containing LafrOBP1 (XP_023395442.1) protein sequence plus the most common protein contaminants and trypsin, allowing to search Lys-acetylation (Δm =  + 42.01), Ser/Thr/Tyr-phosphorylation (Δm =  + 79.97), and the most common mammals N-linked glycans at Asn and O-linked glycans at Ser/Thr/Tyr, using the same parameters previously set. The max PTM sites per peptide was set to 2.Proteome Discoverer peptide candidates were considered confidently identified only when the following criteria were satisfied: (i) protein and peptide false discovery rate (FDR) confidence: high; (ii) peptide Mascot score:  > 30; (iii) peptide spectrum matches (PSMs): unambiguous; (iv) peptide rank (rank of the peptide match): 1; (v) Delta CN (normalized score difference between the selected PSM and the highest-scoring PSM for that spectrum): 0. Byonic peptide candidates were considered confidently identified only when the following criteria were satisfied: (i) PEP 2D and PEP 1D:  More

DNA sample collectionTo analyze the population structure of wild medaka populations, we selected samples from the DNA collection of Takehana et al.29, deposited in University of Shizuoka. The original DNA collection had been made throughout 1980s and 2000s. The selected samples covered the major mitotypes and contained more than three individuals of each population (Table S11, Fig. 3), which were collected from three collection sites for O. sakaizumii and 12 collection sites for O. latipes. We also examined several artificial strains: HNI and Hd-rR, which are inbred strains derived from O. sakaizumii and O. latipes, respectively, and four Himedaka individuals from commercial stock (Uruma city, Okinawa Prefecture, Japan).In addition, samples were newly collected at Kunigami Village, Okinawa Prefecture. Live fish were anesthetized with MS-222 (aminobenzene methanesulfonate, FUJIFILM Wako Pure Chemical Corporation, Osaka, Japan) and then fixed in 99% ethanol. Genomic DNA was extracted using a DNeasy kit (Qiagen Inc., Hilden, Germany) from ethanol-fixed pectoral fin samples according to the manufacturer’s protocol. The DNA concentration was measured using a spectrophotometer (Nanodrop 1000, Thermo Fisher Scientific, Waltham, Massachusetts, USA), and the DNA was diluted with PCR-grade water to a concentration of c.a. 10 ng/µl (UltraPure™ DNase/RNase-Free Distilled Water, Thermo Fisher Scientific).Ethic statementAll methods were carried out in accordance with the Regulation for Animal Experiments at University of the Ryukyus for handling live fish. All experiments were approved by the Animal Care Ethics Committee of University of the Ryukyus (R2019035). All experimental methods are reported in accordance with ARRIVE guidelines.PCR primer designThe following steps were used to select primers for MAAS (Fig. 1). (1) All possible 10-mer sequence combinations (i.e., 410 = 1,048,576 sequences) were generated in silico. (2) The sequences containing simple sequence repeats, some of which had been used in the MIG-seq method17, were excluded. (3) Sequences containing a functional motif, such as a transcription factor-binding site, were also excluded because they may not be suitable for examining neutral genetic markers. We obtained a catalog of motifs from the JASPAR CORE40 (http://jaspar.genereg.net). (4) To avoid taxon-dependency in primer performance, we used information about the k-mer (k = 10) frequency of reference genomes from multiple phyla. Sequences that showed marked differences in frequency among taxa were excluded. The frequencies of each 10-mer sequence in the reference genomes of 17 species belonging to 12 phyla of metazoa were counted (Table S12) using the “oligonucleotideFrequency” function in the “Biostrings” package ver. 2.441. In each of these taxa, the frequencies of sequences were stratified into three grades ( 103). We then selected the sequences that showed the same grade in more than 80% (14/17) of the species. (5) To avoid synthesizing primer dimers, self-complementary sequences were excluded, taking Illumina adapter sequences (5′-CGCTCTTCCGATCT-3′ and 5′-TGCTCTTCCGATCT-3′) into account. Self-complementation of two bases at the 3′-end or every three continuous bases in primer sequences was then evaluated using a custom script in R ver. 3.5.0 (R Development Core Team, http://cran.r-project.org). Based on the selected 10-mer sequences (i.e., 129 sequences, Fig. 1), 7-mer primer sequences were designed by removing the 3 bases at the 3′ end. Finally, we selected 24 candidate sequences for both 10-mer and 7-mer primers for the subsequent step (Table S1).The primer sequence consisted of three parts17: partial sequence of the Illimina adapter, 7 N bases, and a short priming sequence, e.g., 5′-CGCTCTTCCGATCTNNNNNNNGTCGCCC-3′. PCR amplification was performed using the candidate primers using the first PCR protocol described below (Table S1). Banding patterns were observed by electrophoresis on 1% agarose gels (agarose S; TaKaRa, Japan). Of the candidate primers, we selected four 7-mer primers and four 10-mer primers that each gave a smeared banding pattern with amplification products ranging from 500 to 2000 bp, indicating uniform amplification of multiple target sequences (Table S1).Library construction and sequencingThe library was constructed by a two-step PCR approach using a modification of a MIG-seq protocol14. In the first PCR step, multiple regions of genomic DNA were amplified using a cocktail of primers with a Multiplex PCR Assay Kit Ver.2 (TaKaRa) (Table 1). The volume of the PCR reaction mixture was 7 μl, containing 1 μl of template DNA, 2 μM of each PCR primer, 3.5 μl of 2 × Multiplex PCR Buffer, and 0.035 μl of Multiplex PCR Enzyme Mix. PCR was performed under the following conditions: denaturation at 94 °C for 1 min; 25 cycles of 94 °C for 30 s, 38 °C for 1 min, and 72 °C for 1 min, followed by a final extension step at 72 °C for 10 min.The primers in the second PCR step contained the Illumina sequencing adapter and an index sequence to identify each sample. Following the Truseq indexes, we used the combinations of eight forward indexes (i5) and 12 reverse indexes (i7), which resulted in a total of 96 combinations. To be used as a template for the second PCR, the first PCR product from each sample was diluted 50 times with PCR-grade water. The second PCR was performed in a 15-μl reaction mixture containing, 3 μl of diluted first PCR product, 3 μl of 5 × PrimeSTAR GXL Buffer, 200 μM of each dNTP, 0.2 μM of forward index primer and reverse index primer, 0.375 U of PrimeSTAR GXL DNA Polymerase (TaKaRa). The PCR conditions were as follows: 12 cycles at 98 °C for 10 s, 54 °C for 15 s, and 68 °C for 30 s.The second PCR product of each sample was pooled by equal volume and size-selected from 600 to 1000 bp using solid phase reversible immobilization (SPRI) select beads (Beckman Coulter Inc, Brea, California, USA) according to the manufacturer’s protocol. The DNA concentration of the pooled library was measured using a Qubit fluorometer (Thermo Fisher Scientific). We sequenced the libraries using two NGS platforms, MiSeq (Illumina, MiSeq Reagent Kit v2 Micro, Paired-End (PE), 150 bp) and HiSeq X (Illumina, PE, 150 bp). Sequencing using the HiSeq X platform was performed by Macrogen Japan (Tokyo, Japan).To compare primer performance, the DNA libraries constructed using the 7-mer and 10-mer primers for one individual were sequenced using MiSeq. Then, a 7-mer primer cocktail containing four sets of mixed primers was used for the subsequent analyses (Table 1). We also constructed DNA libraries using 7-mer and MIG-seq primer cocktails for three individuals and sequenced them using the HiSeq X platform. Finally, we constructed DNA libraries using 7-mer primer cocktails for 67 wild individuals and six artificial strain individuals for population genetics analyses (Table S11, Fig. 3).Mapping and SNV callingGenotyping was conducted using the following BWA-GATK best-practices pipeline for each sample42. Primer sequences were removed using cutadapt with the –b option selected43. The Illumina adapter sequences were also removed and quality filtering was performed using fastp ver. 0.20.0 with the “–detect_adapter_for_pe, –cut_front” option selected44. The remaining reads were mapped on the reference genome of medaka, Hd-rR strain, GCA_002234675.1; ASM223467v127 using Burrows-Wheeler Alignment tool, BWA mem ver. 0.7.1745. After mapping, output files were converted to Binary Alignment/Map (BAM) format using SAMtools ver. 1.746. SNVs and InDels in the sample were determined following the best practice guidelines set out in the Genome Analysis Tool Kit (GATK ver. 3.8.0)42. We then filtered out SNVs and InDels based on the following criteria: “QD  60.0 || MQ  More

Howe, H. F. & Smallwood, J. Ecology of seed dispersal. Annu. Rev. Ecol. Syst. 13, 201–228 (1982).Article 

Google Scholar 
Wang, B. C. & Smith, T. B. Closing the seed dispersal loop. Trends Ecol. Evol. 17, 379–385 (2002).Article 

Google Scholar 
Fleming, T. H., Breitwish, R. & Whitesides, G. H. Patterns of tropical vertebrate frugivore diversity. Annu. Rev. Ecol. Syst. 18, 91–109 (1987).Article 

Google Scholar 
Schupp, E. W. Quantity, quality and the effectiveness of seed dispersal by animals. Vegetatio 107(108), 15–29 (1993).Article 

Google Scholar 
Garber, P. A. & Lambert, J. E. Primates as seed dispersers: Ecological processes and directions for future research. Am. J. Primatol. 45, 3–8 (1998).Article 
CAS 
PubMed 

Google Scholar 
Godínez-Alvarez, H. & Jordano, P. An empirical approach to analysing the demographic consequences of seed dispersal by frugivores. In Seed Dispersal: Theory and its Application in a Changing World (eds Dennis, A. J. et al.) 391–406 (CAB International, 2007).Chapter 

Google Scholar 
Schupp, E. W., Jordano, P. & Gomez, J. M. Seed dispersal effectiveness revisited: A conceptual review. New Phytol. 188, 333–353 (2010).Article 
PubMed 

Google Scholar 
McConkey, K. R., Brockelman, W. Y. & Saralamba, C. Mammalian frugivores with different foraging behavior can show similar seed dispersal effectiveness. Biotropica 46, 647–651 (2014).Article 

Google Scholar 
Culot, L., Huynen, M. C. & Heymann, E. W. Partitioning the relative contribution of one-phase and two-phase seed dispersal when evaluating seed dispersal effectiveness. Methods Ecol. Evol. 6, 178–186 (2015).Article 

Google Scholar 
McConkey, K. R., Brockelman, W. Y., Saralamba, C. & Nathalang, A. Effectiveness of primate seed dispersers for an “oversized’’ fruit, Garcinia benthamii. Ecology 96, 2737–2747 (2015).Article 
PubMed 

Google Scholar 
Camargo, P., Martins, M. M., Feitosa, R. M. & Christianini, A. V. Bird and ant synergy increases the seed dispersal effectiveness of an ornithochoric shrub. Oecologia 181, 507–518 (2016).Article 
ADS 
PubMed 

Google Scholar 
McConkey, K. R. et al. Different megafauna vary in their seed dispersal effectiveness of the megafaunal fruit Platymitra macrocarpa (Annonaceae). PLoS One 13, e0198960 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Balcomb, S. R. & Chapman, C. A. Bridging the gap: Influence of seed deposition on seedling recruitment in a primate-tree interaction. Ecol. Monogr. 73, 625–642 (2003).Article 

Google Scholar 
Russo, S. E. Responses of dispersal agents to tree and fruit traits in Virola calophylla (Myristicaceae): Implications for selection. Oecologia 136, 80–87 (2003).Article 
ADS 
PubMed 

Google Scholar 
Gross-Camp, N. D., Masozera, M. & Kaplin, B. A. Chimpanzee seed dispersal quantity in a tropical montane forest of Rwanda. Am. J. Primatol. 71, 901–911 (2009).Article 
PubMed 

Google Scholar 
Jordano, P. & Schupp, E. W. Seed disperser effectiveness: The quantity component and patterns of seed rain for Prunus mahaleb. Ecol. Monogr. 70, 591–615 (2000).Article 

Google Scholar 
Leighton, M. Modeling dietary selectivity by Bornean orangutans: evidence for integration of multiple criteria in fruit selection. Int. J. Primatol. 14, 257–313 (1993).Article 

Google Scholar 
Stevenson, P. R. Fruit choice by woolly monkeys in Tinigua National Park, Colombia. Int. J. Primatol. 25, 367–381 (2004).Article 

Google Scholar 
Palacio, F. X. & Ordano, M. The strength and drivers of bird-mediated selection on fruit crop size: A meta-analysis. Front. Ecol. Evol. 6, 18 (2018).Article 

Google Scholar 
Flörchinger, M., Braun, J., Böhning-Gaese, K. & Schaefer, H. M. Fruit size, crop mass, and plant height explain differential fruit choice of primates and birds. Oecologia 164, 151–161 (2010).Article 
ADS 
PubMed 

Google Scholar 
Wenny, D. G. Seed dispersal, seed predation, and seedling recruitment of a neotropical montane tree. Ecol. Monogr. 70, 331–351 (2000).Article 

Google Scholar 
Calviño-Cancela, M. Spatial patterns of seed dispersal and seedling recruitment in Corema album (Empetraceae): The importance of unspecialized dispersers for regeneration. J. Ecol. 90, 775–784 (2002).Article 

Google Scholar 
Masaki, T. & Nakashizuka, A. Seedling demography of Swida controversa: Effect of light and distance to conspecifics. Ecology 83, 3497–3507 (2002).Article 

Google Scholar 
Vulinec, K. Dung beetle communities and seed dispersal in primary forest and disturbed land in Amazonia. Biotropica 34, 297–309 (2002).Article 

Google Scholar 
Janzen, D. H. Herbivores and the number of tree species in tropical forests. Am. Nat. 104, 501–528 (1970).Article 

Google Scholar 
Connell, J. H. On the role of natural enemies in preventing competitive exclusion in some marine animals and in rain forest trees. In Dynamics of Populations (eds Den Boer, P. J. & Gradwell, G.) 298–312 (PUDOC, 1971).
Google Scholar 
Howe, H. F., Schupp, E. W. & Westley, L. C. Early consequences of seed dispersal for a Neotropical tree (Virola surinamensis). Ecology 66, 781–791 (1985).Article 

Google Scholar 
Valenta, K. & Fedigan, L. M. Spatial patterns of seed dispersal by white-faced capuchins in Costa Rica: Evaluating distant-dependent seed mortality. Biotropica 42, 223–228 (2010).Article 

Google Scholar 
Andresen, E. Seed dispersal by monkeys and the fate of dispersed seeds in a Peruvian rain forest. Biotropica 31, 145–158 (1999).
Google Scholar 
Dausmann, K. H., Glos, J., Linsenmair, K. E. & Ganzhorn, J. U. Improved recruitment of a lemur-dispersed tree in Malagasy dry forests after the demise of vertebrates in forest fragments. Oecologia 157, 307–316 (2008).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Jordano, P. & Herrera, C. M. Shuffling the offspring: uncoupling and spatial discordance of multiple stages in vertebrate seed dispersal. Ecoscience 2, 230–237 (1995).Article 

Google Scholar 
Rey, P. J. & Alcantara, J. M. Recruitment dynamics of a fleshy-fruited plant (Olea europaea): Connecting patterns of seed dispersal to seedling establishment. J. Ecol. 88, 622–633 (2000).Article 

Google Scholar 
Traveset, A., Gulias, J., Riera, N. & Mus, M. Transition probabilities from pollination to establishment in a rare dioecious shrub species (Rhamnus ludovici salvatoris) in two habitats. J. Ecol. 91, 427–437 (2003).Article 

Google Scholar 
Cordeiro, N. J. & Howe, H. F. Forest fragmentation severs mutualism between seed dispersers and an endemic African tree. Proc. Natl. Acad. Sci. USA 100, 14052–14056 (2003).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Blendinger, P. G., Blake, J. G. & Loiselle, B. A. Connecting fruit production to seedling establishment in two co-occurring Miconia species: consequences of seed dispersal by birds in upper Amazonia. Oecologia 167, 61–73 (2011).Article 
ADS 
PubMed 

Google Scholar 
Corlett, R. T. The impact of hunting on the mammalian fauna of tropical Asian forests. Biotropica 39, 292–303 (2007).Article 

Google Scholar 
Nuñez-Iturri, G., Olsson, O. & Howe, H. F. Hunting reduces recruitment of primate-dispersed trees in Amazonian Peru. Biol. Conserv. 141, 1536–1546 (2008).Article 

Google Scholar 
Effiom, E. O., Nuñez-Iturri, G., Smith, H. G., Ottosson, U. & Olsson, O. Bushmeat hunting changes regeneration of African rainforests. Proc. R. Soc. B-Biol. Sci. 280, 20130246 (2013).Article 

Google Scholar 
Harrison, R. D. et al. Consequences of defaunation for a tropical tree community. Ecol. Lett. 16, 687–694 (2013).Article 
PubMed 

Google Scholar 
Fuentes, E. R., Hoffmann, A. J., Poiani, A. & Alliende, M. C. Vegetation change in large clearings: patterns in the Chilean matorral. Oecologia 68, 358–366 (1986).Article 
ADS 
PubMed 

Google Scholar 
Holl, K. D. Do bird perching structures elevate seed rain and seedling establishment in abandoned tropical pasture?. Restor. Ecol. 6, 253–261 (1998).Article 

Google Scholar 
Beltran, L. C. & Howe, H. F. The frailty of tropical restoration plantings. Restor. Ecol. 28, 16–21 (2020).Article 

Google Scholar 
Bollen, A., Van Elsacker, L. & Ganzhorn, J. U. Relations between fruits and disperser assemblages in a Malagasy littoral forest: a community-level approach. J. Trop. Ecol. 20, 599–612 (2004).Article 

Google Scholar 
Ganzhorn, J. U. et al. Possible fruit protein effects of primate communities in Madagascar and the Neotropics. PLoS One 4, e8253 (2009).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Albert-Daviaud, A. et al. The ghost fruits of Madagascar: identifying dysfunctional seed dispersal in Madagascar’s endemic flora. Biol. Conserv. 242, 108438 (2020).Article 

Google Scholar 
Dew, J. L. & Wright, P. Frugivory and seed dispersal by four species of primates in Madagascar’s eastern rain forest. Biotropica 30, 425–437 (1998).Article 

Google Scholar 
Moses, K. L. & Semple, S. Primary seed dispersal by the black-and-white ruffed lemur (Varecia variegata) in the Manombo forest, south-east Madagascar. J. Trop. Ecol. 27, 529–538 (2011).Article 

Google Scholar 
Sato, H. Frugivory and seed dispersal by brown lemurs in a Malagasy tropical dry forest. Biotropica 44, 479–488 (2012).Article 

Google Scholar 
Razafindratsima, O. H., Jones, T. A. & Dunham, A. E. Patterns of movement and seed dispersal by three lemur species. Am. J. Primatol. 76, 84–96 (2014).Article 
PubMed 

Google Scholar 
Steffens, K. J., Sanamo, J. & Razafitsalama, J. The role of lemur seed dispersal in restoring degraded forest ecosystems in Madagascar. Folia Primatol. 93, 1–19 (2022).Article 

Google Scholar 
Ganzhorn, J. U., Fietz, J., Rakotovao, E., Schwab, D. & Zinner, D. Lemurs and the regeneration of dry deciduous forest in Madagascar. Conserv. Biol. 13, 794–804 (1999).Article 

Google Scholar 
Razafindratsima, O. H. et al. Consequences of lemur loss for above-ground carbon stocks in a Malagasy rainforest. Int. J. Primatol. 39, 415–426 (2018).Article 

Google Scholar 
Razafindratsima, O. H. Seed dispersal by vertebrates in Madagascar’s forests: Review and future directions. Madag. Conserv. Dev. 9, 90–97 (2014).Article 

Google Scholar 
Nathan, R. & Muller-Landau, H. C. Spatial patterns of seed dispersal, their determinants and consequences for recruitment. Trends Ecol. Evol. 15, 278–285 (2000).Article 
CAS 
PubMed 

Google Scholar 
Pérez-Ramos, I. M., Urbieta, I. R., Maranon, T., Zavala, M. A. & Kobe, R. K. Seed removal in two coexisting oak species: Ecological consequences of seed size, plant cover and seed-drop timing. Oikos 117, 1386–1396 (2008).Article 

Google Scholar 
Sato, H. Seasonal fruiting and seed dispersal by the brown lemur in a tropical dry forest, north-western Madagascar. J. Trop. Ecol. 29, 61–69 (2013).Article 

Google Scholar 
Chapman, C. A. & Chapman, L. J. Determinants of group size in primates: the importance of travel costs. In On the Move: How and Why Animals Travel in Groups (eds Boinski, S. & Garber, P. A.) 24–42 (The University of Chicago Press, 2000).
Google Scholar 
Sato, H. Diurnal resting in brown lemurs in a dry deciduous forest, northwestern Madagascar: Implications for seasonal thermoregulation. Primates 53, 255–263 (2012).Article 
PubMed 

Google Scholar 
Razanaparany, P. T. & Sato, H. Abiotic factors affecting the cathemeral activity of Eulemur fulvus in the dry deciduous forest of north-western Madagascar. Folia Primatol. 91, 463–480 (2020).Article 

Google Scholar 
Sato, H. Predictions of seed shadows generated by common brown lemurs (Eulemur fulvus) and their relationship to seasonal behavioral strategies. Int. J. Primatol. 39, 377–396 (2018).Article 

Google Scholar 
Agetsuma, N. & Nakagawa, N. Effects of habitat differences on feeding behaviors of Japanese monkeys: Comparison between Yakushima and Kinkazan. Primates 39, 275–289 (1998).Article 

Google Scholar 
Stevenson, P. R. Seed dispersal by woolly monkeys (Lagothrix lagothricha) at Tinigua National Park, Colombia: Dispersal distance, germination rates, and dispersal quantity. Am. J. Primatol. 50, 275–289 (2000).Article 
CAS 
PubMed 

Google Scholar 
Hladik, A. & Miquel, S. Seedling types and plant establishment in an African rain forest. In Reproductive Ecology of Tropical Forest Plants (eds Bawa, K. S. & Hadley, M.) 261–282 (UNESCO and The Parthenon Publishing Group, 1990).
Google Scholar 
Ibarra-Manriquez, G., Ramos, M. M. & Oyama, K. Seedling functional types in a lowland rain forest in Mexico. Am. J. Bot. 88, 1801–1812 (2001).Article 
CAS 
PubMed 

Google Scholar 
Blate, G. M., Peart, D. R. & Leighton, M. Post-dispersal predation on isolated seeds: A comparative study of 40 tree species in a Southeast Asian rainforest. Oikos 82, 522–538 (1998).Article 

Google Scholar 
Hosaka, T. et al. Responses of pre-dispersal seed predators to sequential flowering of Dipterocarps in Malaysia. Biotropica 49, 177–185 (2017).Article 

Google Scholar 
Iku, A., Itioka, T., Shimizu-Kaya, U., Kishimoto-Yamada, K. & Meleng, P. Differences in the fruit maturation stages at which oviposition occurs among insect seed predators feeding on the fruits of five dipterocarp tree species. Entomol. Sci. 21, 412–422 (2018).Article 

Google Scholar 
Kitajima, K. Impact of cotyledon and leaf removal on seedling survival in three tree species with contrasting cotyledon functions. Biotropica 35, 429–434 (2003).Article 

Google Scholar 
Marchand, P. et al. Seed-to-seedling transitions exhibit distance-dependent mortality but no strong spacing effects in a Neotropical forest. Ecology 101, e02926 (2020).Article 
PubMed 

Google Scholar 
Moles, A. T. & Westoby, M. Seedling survival and seed size: a synthesis of the literature. J. Ecol. 92, 372–383 (2004).Article 

Google Scholar 
Sonesson, L. K. Growth and survival after cotyledon removal in Quercus robur seedlings, grown in different natural soil types. Oikos 69, 65–70 (1994).Article 

Google Scholar 
Hubbell, S. P., Condit, R. & Foster, R. B. Presence and absence of density dependence in a neotropical tree community. Philos. Trans. R. Soc. Lond. Ser. B-Biol. Sci. 330, 269–281 (1990).Article 
ADS 

Google Scholar 
Terborgh, J. Enemies maintain hyperdiverse tropical forests. Am. Nat. 179, 303–314 (2012).Article 
PubMed 

Google Scholar 
Godínez-Alvarez, H., Valiente-Banuet, A. & Rojas-Martínez, A. The role of seed dispersers in the population dynamics of the columnar cactus Neobuxbaumia tetetzo. Ecology 83, 2617–2629 (2002).Article 

Google Scholar 
Carson, W. P., Anderson, J. T., Leigh, E. G. Jr. & Schnitzer, S. A. Challenges associated with testing and falsofying the Janzen-Connell hypothesis: a review and critique. In Tropical Forest Community Ecology (eds Carson, W. P. & Schnitzer, S. A.) 210–241 (Wiley-Blackwell, 2008).
Google Scholar 
Swamy, V. et al. Are all seeds equal? Spatially explicit comparisons of seed fall and sapling recruitment in a tropical forest. Ecol. Lett. 14, 195–201 (2011).Article 
PubMed 

Google Scholar 
Reid, J. L. et al. Multi-scale habitat selection of key frugivores predicts large-seeded tree recruitment in tropical forest restoration. Ecosphere 12, e03868 (2021).Article 

Google Scholar 
Steffens, K. J. E. Lemur food plants as options for forest restoration in Madagascar. Restor. Ecol. 28, 1517–1527 (2020).Article 

Google Scholar 
Chapman, C. A. & Dunham, A. E. Primate seed dispersal and forest restoration: an African perspective for a brighter future. Int. J. Primatol. 39, 427–442 (2018).Article 

Google Scholar 
Morris, P. & Hawkins, F. Birds of Madagascar: A Photographic Guide (Yale University Press, 1998).
Google Scholar 
Garbutt, N. Mammals of Madagascar: A Complete Guide (Yale University Press, 2007).
Google Scholar 
Mittermeier, R. A. et al. Lemurs of Madagascar 3rd edn. (Conservation International, 2010).
Google Scholar 
Sato, H., Ichino, S. & Hanya, G. Dietary modification by common brown lemurs (Eulemur fulvus) during seasonal drought conditions in western Madagascar. Primates 55, 219–230 (2014).Article 
PubMed 

Google Scholar 
Laman, T. G. Ficus seed shadows in a Bornean rain forest. Oecologia 107, 347–355 (1996).Article 
ADS 
PubMed 

Google Scholar 
Clark, C. J., Poulsen, J. R., Connor, E. F. & Parker, V. T. Fruiting trees as dispersal foci in a semi-deciduous tropical forest. Oecologia 139, 66–75 (2004).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Sato, H. Gut passage time and size of swallowed seeds in the common brown lemur and the mongoose lemur. Primate Res. 25, 45–54 (2009) (In Japanese with English summary).Article 

Google Scholar 
Menard, S. W. Applied Logistic Regression Analysis 2nd edn. (Sage Publication, 2002).Book 

Google Scholar 
Burnham, K. P. & Anderson, D. R. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach 2nd edn. (Springer, 2002).MATH 

Google Scholar 
Richards, S. A. Testing ecological theory using the information-theoretic approach: Examples and cautionary results. Ecology 86, 2805–2814 (2005).Article 

Google Scholar 
Symonds, M. R. E. & Moussalli, A. A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike’s information criterion. Behav. Ecol. Sociobiol. 65, 13–21 (2011).Article 

Google Scholar  More