A meta-analysis on the impact of trees on yield of intercrops in alley-cropping systems of temperate climates

Review A meta-analysis on the impact of trees on yield of intercrops in alley-cropping systems of temperate climates Anna Panozzo a,* , Paul Quataert b , Tom De Swaef b , Paul Pardon b , Teofilo Vamerali a , Kris Verheyen c , Bert Reubens b a Department of Agronomy, Food, Natural Resources, Animals and the Environment, University of Padova, Legnaro, Italy b Flanders research institute for agriculture, fisheries and food (ILVO), Merelbeke, Belgium c Forest & Nature Lab, Ghent University, Gontrode, Belgium HIGHLIGHTS GRAPHICAL ABSTRACT •Meta-analysis on the grain yield of arable intercrops in temperate alleycropping. •30 % average yield reduction in alleycropping compared to full sun from 18 trials. •The negative effect of trees on yield increases with tree vicinity and tree aging. •Trees have lower detrimental impact on winter, especially barley, vs. summer crops. •Less detrimental or positive impact under low rainfall and high evapotranspiration. ARTICLE INFO Editor: Mark van Wijk Keywords: Agroforestry Distance Maize Silvoarable systems Soybean Tree age Winter cereals ABSTRACT CONTEXT: The environmental benefits of agroforestry have been highlighted worldwide, although improved intercrop productivity has been clearly demonstrated only in the tropics. OBJECTIVE AND METHODS: This meta-analysis aimed at summarizing knowledge from 18 trials on grain yield of arable intercrops in alley-cropping systems of temperate climates, within a mixed-effect model framework. RESULTS AND CONCLUSIONS: A general negative impact of trees on crop grain yield was documented, with an average reduction by 30 % compared to full sun, across the whole inter-row of wheat, barley, soybean and maize. Key findings included: (i) distance from trees is the major driver of the grain yield response, with increasing impact in the vicinity of the tree row, (ii) significance of crop phenology and species choice, with lower impact on winter vs. summer crops; (iii) tree age is the only relevant variable of the woody component, with increasing impact with aging; and (iv) available rainfall and potential evapotranspiration are key moderators, with a less detrimental or positive impact of trees under low rainfall and high evapotranspiration. This study also describes the implications of some tree design and management practices: (i) branchless tree rows allows for halving the alley width where crop yield is negatively affected by the tree line compared to a hedgerow, (ii) crop yield * Corresponding author. E-mail addresses: [email protected] (A. Panozzo), [email protected] (P. Quataert), [email protected] (T. De Swaef), paul. [email protected] (P. Pardon), [email protected] (T. Vamerali), [email protected] (K. Verheyen), [email protected] (B. Reubens). Contents lists available at ScienceDirect Agricultural Systems journal homepage: www.elsevier.com/locate/agsy https://doi.org/10.1016/j.agsy.2025.104578 Received 9 July 2025; Received in revised form 12 November 2025; Accepted 18 November 2025 Agricultural Systems 232 (2026) 104578 Available online 28 November 2025 0308-521X/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). recovers at distances (D) from the tree row at which D/tree heigh approaches 1, this suggesting the interrow should be at least twice the maximum height of trees, (iii) intercrop should be varied across the tree cycle, by cultivating high shade-tolerant species/varieties when the tree age is >8 years. SIGNIFICANCE: This meta-analysis underscores the need for further empirical studies on other intercrops and within climatic zones where limited data is currently available. 1. Introduction The role of silvoarable practices in mitigating the impact of climate change on farming systems, while providing a considerable variety of ecosystem services, has been largely documented both in tropical and temperate environments. Numerous review studies and meta-analyses, particularly from tropics and subtropics, have summarized the benefits of integrating trees with crops, with a focus on carbon sequestration (Kumara et al., 2023; Cerda et al., 2019; Santos et al., 2019) and soilmediated ecosystem services (Muchane et al., 2020; Kuyah et al., 2019). At global scale, the transition from conventional agriculture to agroforestry might significantly increase soil organic carbon (SOC) stocks by 26 %, 40 %, and 34 % in the 0–15, 0–30, and 0–100 cm depth interval respectively (De Stefano and Jacobson, 2018), with averagely 126 t ha −1 additional carbon within 1 m depth soil layer (Shi et al., 2018) and 46.1 t ha −1 within the tree biomass compared with sole croplandor pasture-based land uses (Ma et al., 2020). In tropical areas, SOC can peak within five years, reaching equilibrium rapidly (Ma et al., 2020), whereas in temperate regions, SOC increases at a slower rate, requiring approximately a decade to observe measurable benefits and likely 35 years to reach equilibrium (Pardon et al., 2017; Ivezi´ c et al., 2022; Mayer et al., 2022). While capturing CO 2 from the atmosphere into stable carbon pools, agroforestry tree systems provide numerous ecosystem services at both below and aboveground level that sustain stable crop production. In humid and sub-humid tropics, agroforestry can reduce soil erosion by up to 50 % and improve water balance compared to monoculture (Muchane et al., 2020). Higher soil moisture content in tropical agroforestry systems is linked to improved soil infiltration rates and reduced evapotranspiration, facilitated by better soil hydraulic conductivity, soil porosity, and microclimate modifications associated to windbreak effect (Nyamadzawo et al., 2008; Kuyah et al., 2019; Siriri et al., 2013). Improvements in the nutrients cycle have also been extensively documented in the tropics, as a result of higher overall above and belowground C input by pruning residues, litterfall, root exudates and turnover, leading to 46 % and 11 % increases of N and available P, respectively, under agroforestry compared to monoculture (Muchane et al., 2020). In temperate climates, similar patterns have been highlighted by various empirical studies, followed by some attempts to summarize the key findings across the whole temperate zone (Dmuchowski et al., 2024; Sollen-Norrlin et al., 2020) or Europe (Torralba et al., 2016). The major findings indicate that trees in agroforestry can mitigate microclimate extremes, reducing evapotranspiration, soil erosion, while improving soil nutrients availability (Cardinael et al., 2021; Dmuchowski et al., 2024) and enhancing biodiversity (Torralba et al., 2016; Edo et al., 2024). Despite extensive recognition of agroforestry’s environmental benefits, evidence for enhanced intercrop productivity is mostly limited to tropical climates, where under-canopy crops like coffee and cacao perform well in shaded conditions (Nesper et al., 2017; Mokondoko et al., 2022; Mattalia et al., 2022), and in poor soils and arid climates (Bayala and Prieto, 2020; Kuyah et al., 2019; Durand-Bessart et al., 2020). In temperate zones, crop yield response is variable, due to interspecific competition for natural resources, land use, and diverse pedo-climatic conditions. Evidence suggests that silvoarable systems often outperform monocultures regarding land use efficiency, as supported by favorable land equivalent ratios (LER) (Mead and Willey, 1980) ranging from 0.95 to 2, depending on pedo-climatic zones and agricultural management (Sereke et al., 2015; García de Jal´ on et al., 2018; Xu et al., 2019; Lehmann et al., 2020; Pent, 2020). While early research focused primarily on environmental benefits, recent studies have demonstrated that higher LERs in silvoarable systems result from improved use efficiency of both agricultural land and natural resources, such as solar radiation and water. These resources limited availability concerns an increasing number of farmers in temperate climates (Graves et al., 2007; Wolz and DeLucia, 2018). Although farmers recognize these benefits, they often perceive silvoarable practices as economically unviable, as revealed by several surveys assessing farmers’ opinions on barriers to the integration of trees into agricultural land (Smith et al., 2012; Trozzo et al., 2014; Rois-Díaz et al., 2018; Sereke et al., 2015). This is among the main obstacles to a broader adoption of silvoarable practices across temperate climate zones (Eichhorn et al., 2006; García de Jal´ on et al., 2018; Bettles et al., 2021; Tavernier et al., 2024; Tranchina et al., 2024), depending on scattered and contrasting information on yield of intercrops available in the literature. Indeed, across various crop and tree species, design characteristics, management choices and pedo-climatic conditions, crop yield in silvoarable systems was either reduced (Burgess et al., 2005; Tsonkova et al., 2012; Swieter et al., 2019; Inurreta-Aguirre et al., 2018; Pardon et al., 2018) or enhanced (Moreno, 2008; Piotto et al., 2024; Arenas-Corraliza et al., 2021) as compared to full sun conditions without trees. Comprehensive information on agroforestry systems of the temperate climates derived from case studies that are summarized by meta-analyses, as attempts to provide synthetic key drivers of crop yield response. Van Vooren et al. (2017), by focusing on silvoarable systems with hedgerows (HR) highlighted a significant effect of the distance from trees on crop yield, with reduced productivity close to the HR and gradually restoring when the distance-totree height ratio (D/H) increases. Over 60 studies from temperate climates, these authors have demonstrated a negative impact on crop yield within a D/H 0–2.1 interval, estimating a relative crop yield within this area of 79 %, although neither crop choice nor HR characteristics (e.g., species, age, orientation) were found significant explanatory variables. The adverse effect of reduced distance from trees on relative crop yield was also confirmed by Ivezi´ c et al. (2021), who investigated modern alley cropping systems and traditional systems with scattered trees such as Dehesas, making comparisons between Northern and Southern Europe. Over 13 European studies, these researchers documented that the relative crop yield tend to be higher in Dehesa than in alley cropping, while the latter provided similar impact in northern and southern Europe, and relative yield of fodder tend to be lower than those of cereals in alley cropping (Ivezi´ c et al., 2021). More recently, Scordia et al. (2023) investigated the impact of the tree species, tree cover, crop species, and agronomic management on crop productivity, with a focus on agroforestry systems from semiarid to warm temperate Mediterranean areas. In such a study, the distance from trees was not investigated, but crop and tree species choice, and soil coverage by trees resulted as significant explanatory variables, across 21 Mediterranean studies (Scordia et al., 2023). Within this framework, a comprehensive summary of crop yield responses in temperate alley-cropping systems is still lacking. There currently is an urgent need to summarize results from a wide range of pedo-climatic conditions in the temperate area, for identifying the best system designs, tree management and tree-crop combinations for maximizing crop productivity in alley silvoarable systems. This study aims to integrate existing knowledge on crop yield in silvoarable alley A. Panozzo et al. Agricultural Systems 232 (2026) 104578 2 cropping systems from diverse temperate regions, with a focus on arable grain crops and alley-cropping systems. Grain crops play a major role in the diet of a vast part of the global population, but increasingly affected by climate extremes (Challinor et al., 2014; Rezaei et al., 2023). Alleycropping is the most promising design for the implementation of highproductive silvoarable systems in temperate regions, as a result of the facilitation of mechanical operation as compared to scattered designs. Specifically, adressed research questions (RQ) and hypotheses (H) were as follows: RQ.1. Does the impact of trees on the grain yield of intercrops vary according to the distance from the tree row? H.1. The most relevant yield reductions are expected close to the tree line. RQ.2. Does the impact of trees vary according to crop phenology and species? H.2. Summer crops are expected to be more negatively affected by trees as compared to winter crops. RQ.3. Are there specific tree characteristics and system design that attenuate the impact on grain yield of intercrops? H.3. Tree age and row orientation are expected to be relevant moderators. RQ.4. Does the impact of trees vary according to climatic variables? H.4. Yield reduction of intercrops in the vicinity of the tree row is expected to be attenuated under high temperature and low rainfall. 2. Materials and Methods 2.1. Literature search A systematic literature search was conducted using Scopus and Web of Science, for peer-reviewed studies focused on grain crop yield response in silvoarable agroforestry systems (AF) compared to adjacent controls under full sun without trees (C). This search covered studies from January 1980 to December 2022 (see detailed search string in SI.1). Silvoarable agroforestry is defined as intercropping with widely spaced tree rows (Dupraz and Newman, 1997), where trees can be distributed following various designs, such as alley cropping, isolated/ scattered trees or line belts (Mosquera-Losada et al., 2009). The main purpose here was to summarize the effect of row-arranged trees on the yield of the neighboring arable crops, thereby the definition of silvoarable was restricted to alley cropping and line belts, excluding studies with isolated/scattered trees. Case studies on hedgerow were included only when trees were pruned regularly and shrubs were absent under the trees; the other cases of hedgerow, as well as pollarded trees were excluded from the analysis. The systematic literature search was performed accordingly to the PRISMA guidelines (Moher et al., 2009) and the process is described in SI.2. After the removal of duplicates, the literature search yielded 1109 references. A primary screening of the candidate studies for further analysis was carried out based on title and abstract, with the following criteria: (i) study site is located within temperate regions defined by K¨ oppen-Geiger climate classification, regardless Northern and Southern hemisphere (Beck et al., 2018; SI.3), (ii) empirical data on intercrop grain yield are available, but reviews and modelling studies were excluded, (iii) studies are located in natural open field with real tree rows: pot and greenhouse studies, as well as artificial shading trials were excluded, (iv) true controls are present nearby the silvoarable field/plot, allowing yield comparison with and without tree rows. With the use of these criteria, 15 studies were retained. Additionally, to retrieve as many relevant studies as possible, a further search through the reference lists of these papers was carried out. In total, 18 original research studies that met the inclusion criteria were selected for further analysis of their data collection (see SI.4, SI.5). Of these, 10 studies were from Europe, 1 from North America, 3 from South America, 1 from Africa and 3 from Asia (SI.6) (Fig. 1). 2.2. Data extraction from selected primary studies For each selected study, the mean, sample size (when available) and relative estimates of uncertainty such as standard error or standard deviation of intercrop grain yield under the presence of the tree row (AF) and without (C), were extracted. Data were retrieved directly from tables or the main text, and the software WebPlotDigitizer (Rohatgi, 2015) was used to retrieve data available from graphs only. When available, data across different sites (i.e., coordinate locations), and years within the same study were extracted. Following the first hypothesis (H.1) regarding yield variation across the inter-row, for each site and year, all the available data of crop yield at different distances and both sides from the tree rows were extracted. Each observation (equal to one row in the dataset) represented a comparison between the grain yield of a specific crop in AF at a defined distance and side from the tree row, and the yield of the same crop in an adjacent control setting without trees (C). When only an average yield for the entire area between two adjacent tree rows was reported, the midpoint of the inter-row was used as the distance reference, with the following number of observations for each crop: 102 for wheat, 33 for barley, 5 for triticale, 26 for soybean, and 51 for maize (SI.6). Since the Fig. 1. Geographical distribution of the 18 original studies identified through systematic literature search. For each study, point size depict the number of observations available and the colors the temperate climate subgroups, i.e., mediterranean (MC; including both warm and hot summers – mediterranean climates sub groups), oceanic (OC), and warm temperature with hot summer (WT). A. Panozzo et al. Agricultural Systems 232 (2026) 104578 3 number of observations for triticale was excessively low to allow any robust analysis, this crop was excluded from the study. Wheat and barley, being always cultivated during winter-spring months, were also referred to as winter crops, while soybean and maize as summer crops. Additionally, a set of moderator variables was collected to potentially explain heterogeneity in effect sizes. These moderators included (i) study site characteristics (latitude, longitude, altitude, country, and city), (ii) climatic variables, (iii) crop and tree attributes, and (iv) agroforestry design parameters. Since the climatic variables were often absent or reported for incomplete measurement period, they were retrieved from the Comprehensive Environmental Data Archive (CEDA), using the Climatic Research Unit (CRU) Time-Series (TS) version 4.04 of high-resolution (0.5 ×0.5 degree) gridded data of month-by-month variation in climate (Jan. 1901Dec. 2019) (Harris et al., 2020). TS are produced by CRU at the University of East Anglia funded by the UK National Centre for Atmospheric Science (NCAS), as a NERC collaborative Centre. Given the coordinates of each study site, the monthly average air temperature (temp; degree Celsius) and potential evapotranspiration (pet; mm), and the monthly rainfall (rain; mm) were extracted from CEDA archive. The aridity index (arid) was calculated as the ratio between rain and pet (Zomer et al., 2008). For each study site and year the above four climatic variables were expressed either as growing season means (tempG, rainG, petG, aridG) and as yearly mean (tempY, rainY, petY, aridY). The G variables were calculated by averaging the monthly climatic data between sowing and harvesting. For winter crops, which grow across two subsequent calendar years, the Y variables considered the year of the vegetative and reproductive growth, i.e., the second calendar year only. Each study site was classified also by the temperate climate subgroup of belonging, according to the K¨ oppen-Geiger classification (Beck et al., 2018). Accordingly, three climate subgroups were identified: mediterranean (MC; including both the warm and the hot summers – mediterranean climates subgroups), oceanic (OC), and warm temperature with hot summer (WT). Lastly, each study site was defined by the categorical ‘contrast_zone_TR4’ variable, that was created to categorize rainfall and temperature combinations at each site, as described in Table 1, and further detailed in SI.7 and SI.8. Crop-specific data included crop species, sowing and harvesting times, and side and distance from the tree row; the cultivar name was also extracted as potential moderator, although available only in about half of the observations. Regarding tree characteristics, the following moderators were collected, when available: tree species, cultivar/clone, age, height, diameter at breast height (DBH) and canopy width. Design characteristics of the agroforestry systems, that were collected from primary studies, were: tree row orientation, inter-row distance, intrarow distance and population density. When only part of the above information was reported in a study, the authors were contacted to gaining the missing information. Different formats of potential explanatory variables were tested, either numerical or categorical, as described in Table 1, and further detailed in SI.7. 2.3. Effect size calculation The log response ratio (RR) was calculated as a measure of the effect size (Hedges et al., 1999): RR =ln(ResponseAF ResponseC) For each observation, the response ratio was calculated between the mean response, i.e., crop grain yield at harvest in agroforestry (Response AF ) at a specific distance from the tree row, and the mean response of that crop in the adjacent open field control (Response C ) without trees. In this way, a positive response ratio (RR >0) means that the presence of the tree row has a positive effect on crop yield, while a negative value (RR <0) a negative effect, and close to zero (RR ≈0) Table 1 List of explanatory variables tested along model screening, both in numerical (N) and categorical (C) formats, with reference to the model group (for more details see: Table 2, SI.7, SI.8 and SI.12). For climate variables: G =data for growing season, Y =yearly data. Model group Variable extended Variable code Variable format Description BASE Distance from tree-row ln_DH N ln(D/H) D/H C D/H ≤0.25; D/H >0.25 & ≤0.50; D/H >0.50 & ≤1.00; D/H >1.00 CORE Crop phenology crop_pheno C winter crops (barley, common wheat); summer crops (soybean, maize) Crop species crop_species C barley; wheat; maize; soybean AGE Tree age tree_age_orig N original variable tree_age_log2 N log2(original variable) tree_age_C2 C young (≤8 years old); old (>9 years old) tree_age_C3 C young (≤8 years old); middle (≤9 & >16 years old); old (>17 years old) TREEROW Tree height tree_height N original variable Tree species tree_species C 13 species (for a complete list, see Table SI.7) Tree growth cycle tree_growth_C2 C Fast versus slow growing (for a complete list, see Table SI.7) Intra-row distance tree_dinrow_C3 C low (≤4 m); intermediate (≤4 & >7 m); high (>7 m) Tree-row orientation tree_orient_C3 C N-S (North-South); E-W (East-West), NW-SE (North West – South East) Side of the tree-row tree_side_C3 C E (East); W (West); Others (North, South, Centre of the alley) CLIMATE Temperate climate subgroup (K¨ oppen-Geiger categories) clim_zone_KG C MC (warm summer and hot summer mediterranean climates); OC (oceanic climate); WT (warm temperature with hot summer climate) Contrasting climatic conditions clim_zone_TR4 C high ◦C & high mm (14–26 ◦C & 500–900 mm); low ◦C & high mm (7–14 ◦C & 500–900 mm); high ◦C & low mm (14–26 ◦C & 100–500 mm); low ◦C and low mm (7–14 ◦C & 100–500 mm) Air temperature (X =G/Y) tempX_orig N original variable tempX_log2 N log2 (original variable) tempX_pred C low, middle, high (for the class boundaries, see SI.8) Rainfall (X =G/Y) rainX_orig N original variable rainX_log2 N log2 (original variable) rainX_pred C low, middle, high (for the class boundaries, see SI.8) Potential evapotranspiration (X = G/Y) petX_orig N original variable petX_log2 N log2 (original variable) petX_pred C low, middle, high (for the class boundaries, see SI.8) Aridity index (X =G/Y) aridX_orig N original variable aridX_log2 N log2 (original variable) aridX_pred C low, middle, high (for the class boundaries, see SI.8) A. Panozzo et al. Agricultural Systems 232 (2026) 104578 4 little or no effect. Since most of the studies provided yield data for multiple distances from the tree rows, RR was calculated for each distance relative to the control fields. Since this calculation violates independence assumptions of the effect sizes (Olkin and Gleser, 2009; i.e., multiple Response AF values compared to one Response C ), this was accounted for the correlation among effect sizes by computing the variance–covariance (VCV) matrix proposed by Lajeunesse (2011), using the ‘v.calc’ function of the R package metafor. The VCV matrix models the dependencies that arise when using the same control group and estimating multiple effect sizes (Lajeunesse, 2016). The inverse of the square root of the sampling variance of the VCV matrix was used to weight the precision of the effect size (Lajeunesse, 2011). To compute the VCV matrix, the sampling variance of each pairwise comparison was calculated based on Hedges et al. (1999), as ~follows: σ (RR) = (SDC NCResponseC)+(SDAF NAFResponseAF) where SD is the standard deviation of the mean and N denotes the number of replicates (plots) for AF and C settings. As the estimate of uncertainty of grain yield was not always available in the studies, when this information was not obtained by the authors, missing SD values were imputed using the coefficient of variation from all complete cases using the ‘impute_SD’ function of the R package metagear (Lajeunesse, 2016). In this way, it was possible to apply a conventional meta-analysis with weighting of the different observations by their variances on the full dataset (Gurevitch et al., 2018). To test the overall effect of agroforestry on grain yield, the summary effect size was calculated, as well as the total heterogeneity ( τ 2 ) across different studies, considering 95 % credible (i.e., Bayesian confidence) intervals (CI). The existence of publication bias was assessed using funnel plots (Peters et al., 2008), and specifically the ‘funnel’ function of the R package metafor (Viechtbauer, 2010); the Egger’s test modification proposed by Nakagawa and Santos (2012) were applied to assess funnel plots’ asymmetry of the null models’ residuals. 2.4. A mixed-effect regression approach A mixed-effect regression model framework was set up to investigate and quantify the four research hypotheses (RQ1 – RQ4). Regression modelling was possible as the full datasets of the studies were available with yield expressed along with its distance to the tree row and many other background characteristics. The overarching format of the models was as follows: RR∼f1(Distance)*(f2(Crop)+f3(Treerow)+f4(Climate))+(1|Study:Site) The fixed effect terms refer to the groups of explanatory variables expressing the research hypotheses. Functions fi()() express that different formats of the model terms were explored: for instance, Climate was available as both a categorical variable (climate zone or climate group) and with the analytical variables rain, temperature, potential evapotranspiration and the aridity index. Moreover, to allow for nonlinear relationships, the numerical variables were expressed in the log-scale or were categorized in two (C2) or three classes (C3). The interaction terms with Distance and the other explanatory variables (as expressed with the interaction symbol *) allowed to modify the slope of the Distance function. The random effect term assigned to each study site modeled the correlation structure within the study sites because of the repeated measures. As expressed by Study:Site, different sites within the same study received a different random effect. Site indicates each experimental site where a crop was investigated for each year of study. Therefore, for each Study there were as many Sites as the number of crop ×year ×site. It was not possible to fit a hierarchical effect (Study/Site) as the dataset was highly unbalanced. Some studies contained many sites with few observations, while others contained few sites or even one site. 2.5. The model building strategy Because of the many variables to explore and in order to guarantee interpretable models, a forward regression strategy was used to build models, and variables were step by step included representing the research hypotheses as summarized in Table 2. At each step, alternative variables were screened, and the best model was selected based on regression diagnostics and visualized for its ecological meaning before continuing. In the first step, the yield expressed as RR was modeled as a function of the distance from the tree row; it was a major feature of the data as it is not meaningful to study the impact of other variables without taking into account this major factor. To allow comparisons among different studies and trials, the distance was expressed in relative terms to the height of the tree row as ratio of the distance from the tree row (D) to the height of the trees (H) (D/H). Tree height was reported for the majority of the studies; when not available, it was estimated based on tree species and age. To linearize the distance trend, D/H was log-transformed (ln (D/H)) and the relation RR ~ ln(D/H) was tested as the base model for step 1 (RQ1). Table 2 Details of the steps of the forward regression, in accordance with the research questions. Model step Research question Type of moderators Competing models List of moderators Output model Step 1 RQ1 Crop distance from the tree row Null model Crop distance Base Step 2 RQ2 Crop type Crop phenology Base model & Crop phenology Base model & Crop species Core Step 3 RQ3 Tree row characteristics Core model & Tree age Core model & Tree height Core model & Tree species Core model & Tree growth cycle Many competing models Core model & Tree row orientation Core model & Intra-row distance Core model & Side of the tree row Step 4 RQ4 Climatic characteristics Core model & Climatic zone Core model & Temperature Core model & Rainfall Core model & Potential evapotranspiration Core model & Aridity index A. Panozzo et al. Agricultural Systems 232 (2026) 104578 5 This first step resulted in the base model expressing the fundamental relation of yield and distance (RQ1). Next, the specific impact of the crop species and its phenology was investigated (RQ2). The special point of interest was whether crop phenology (winter or summer) could explain an itself important part of the total variability without referring to a specific crop species. Both an additive and interaction relation with ln (D/H) were explored and the two competing phenology and crop models were compared both numerically and graphically. The result of the second step was called the core model as it was the point of departure for further refinements. In the next two steps, many explanatory variables were screened groupwise to improve the core models. In step 3 the impact of the tree characteristics (including tree age) and the tree row design (RQ3) were explored extending the core model, as well as in step 4 with the climate variables (RQ4). As the relationship is not necessarily linear, log-transformed and categorized variables were considered as detailed in SI.7. This resulted in several competing models investigated globally with both the information indices AIC and BIC and more in detail graphically. The IC gives an overall appreciation of the goodness of fit considering the model complexity, while the graphical approach allowed to appreciate whether model quality was uniform over the range of the explanatory variables. The global screening relied on an information-theoretic perspective (Burnham and Anderson, 2014). Two complementary information criteria (IC) were used: Akaike’s Information Criterion (AIC) and Bayesian Information Criterion (BIC). Both criteria make a trade-off between goodness of fit as expressed by the deviance of the model (lower values are better) and model complexity as expressed by the number of parameters (lower number is better). However, for BIC, the penalty for model complexity is higher, favoring more simple models. Considering both criteria simultaneously provided insights on whether there was no overfitting. To oversee many competing models, a graphical tool was developed to visualize the trade-off between AIC and BIC. In this graph, by group of variables, the competing models were ordered as function of BIC, to facilitate the choice of the optimal model(s) and to contrast the evolution of BIC with AIC. Model screening followed some inclusion criteria: •In case of a severe conflict between AIC and BIC outputs, BIC was considered as a warning for model complexity and overfitting, and the model was rejected; •With similar AIC and BIC outputs, the choice of the optimal model considered (i) the regression diagnostics, (ii) the level of data coverage, and (iii) ecological meaning of the candidate models. Data coverage of a model is related to data balance: for each range or combination of variables, the model should contain sufficient data to verify the model. This is a criterion often overlooked and is linked with the concept of model extrapolation. Regression models are very good in interpolation, but function bad when extended beyond the observation range. Statistical analyses were performed in R version 3.4.4 (The R Foundation for Statistical Computing Platform, 2024) and RStudio version 2023.12.1 (Posit Software, PBC, 2009–2024). The model building with mixed-effects regression models was carried out within the lme4 environment, i.e., with the ‘lmer’ function of the R package lme4 (Bates et al., 2015) in combination with the ggeffects package (Lüdecke, 2018) to visualize the models. 3. Results 3.1. Overall effect, study heterogeneity and publication bias Response ratio ranged from −2.77 (94 % lower yield compared to the reference) for soybean at +2 m on the Northern side from an E-W oriented row of 9-year old Eucalyptus trees in Brazil (study ID#13) up to +0.57 (yield 77 % higher than the reference) for maize at +10 m from a N-S oriented row of 6-year old poplar trees in Italy (study ID#17), both cases falling in the warm temperature with a hot summer climate (WT). The meta-analysis of the entire dataset revealed a significant negative effect of tree row on crop grain yield, with an estimate of the summary effect size equal to −0.32 (relative yield loss by 30 %) [95 % confidence interval from −0.39 to −0.26] (p ≤0.001). Heterogeneity among studies also showed to be significant, indicating high diversity in environmental and weather conditions, and particularities of the site/ study, with an estimated amount of total heterogeneity τ 2 equal to 0.21 [95 % confidence interval =0.20, 0.33] (p ≤0.001), as detailed in SI.9. The model building process, from the null to the full model, showed that the successive inclusion of relevant moderators reduced funnel plots asymmetry (SI.10). The base model with crop distance revealed to be slightly unbalanced. However, fitting the best full models with relevant explanatory variables showed no evidence for publication bias based on both visual estimation of funnel plots and Egger’s test. 3.2. The model building process Fig. 2 and SI.11 show how the model performance changes by inclusion of predictors reflecting the research hypotheses (Table 2) in a forward regression. The top panel represents the starting point, while the panels below build further on the optimal model in the top model. Within each panel, the models are ranked by BIC, and it is checked whether AIC follows the same trend to control for overfitting. If a major discrepancy occurs, then the model is ruled out, as seen in case for instance for the model with tree species. Although AIC favored this model overall, it has very high BIC value indicating overfitting. It is the worst model within the panel investigating the impact of tree row variables. As there are many different tree species, the number of parameters in the model is ultimately very high. This high complexity results in a good fit, but the BIC criterion indicates overfitting due to the high number of parameters involved. As regards crop type variables, i.e., crop species and phenology, they were confirmed to be a significant explanatory variable of the response ratio, as reflected by the significant drop in both AIC and BIC as compared to the base model (SI.11). With crop type variables included in the model, the choice of the linear relation ‘RR ~ ln(D/H)’ over a polynomial one was also verified. Although AIC favored the polynomial fit, the linear relation was kept, due to (i) higher BIC value for the polynomial, warning for model complexity, and (ii) limited range of ln (D/H) covered by the observations, specifically the lack of data >ln(D/ H) =~0. Crop species and phenology allowed for a very similar improvement of the base model, with AIC favoring the first, and BIC the latter, although slightly (SI.11). As the gain in AIC and BIC achieved with these two moderators was the largest among all the competing variables, ‘RR ~ ln(D/H) * crop species’ (crop_DH) and ‘RR ~ ln(D/H) * crop phenology’ (pheno_DH) were both defined as the core models of the dataset and both tracks were investigated further. The inclusion of additional relevant moderators in the subsequent model building was pursued with both core models. Fig. 2 and SI.11 show the AIC and BIC results of model building with pheno_DH, since it scored consistently better on the BIC criterion, as in line with the inclusion criteria for model screening. In the supplementary material, SI.12 shows the results of the alternative track with crop_DH as a starting point. In the model building with pheno_DH, regardless of the variable included, the number of alternative models yielding better than the pheno_DH model was generally limited, highlighting the strength of the pattern fit by this model to the dataset (Fig. 2). Both AIC and BIC indicated that a relevant improvement of the core model could be achieved by including tree age. Across tree age-models, the trajectories of AIC and BIC outputs showed that the relation with age is not linear (as represented by age_orig: original tree age values), while the log2-scale of the original variable gives a better relationship, but does not succeed in modelling the nonlinearity properly, as the two categorizations of age A. Panozzo et al. Agricultural Systems 232 (2026) 104578 6 (three class ‘age_C3’ and two class ‘age_C2’ variables) scored better. More generally, although by definition AIC-values are lower than BICvalues, the ranking of tree age-models remained very similar as seen from both perspectives. In particular, the class variable with two binary categories of age was favored more by BIC than by AIC due to its simplicity, showing the best performance as compared to all the other age variable formats. Conversely, the addition of none of the other tree row characteristics performed better than the pheno_DH model, although similar BIC outputs to the core models were achieved with the distance between trees along the row as three class variable (dinrow_C3) (Fig. 2, SI.11). The addition of age_C2 to the pheno_DH model performed better than any other alternative models, as reflected by both the AIC and BIC indexes, resulting in the best model across all those investigated (Fig. 2). Nevertheless, climatic variables improved model quality (lower BIC values) as compared to the core model, particularly when including cumulative rainfall from sowing-to-harvesting (rainG) and year evapotranspiration (petY), both in the original and log2-scales, while the aridity index (aridG) scored similarly to the core model. In contrast, the two categorical variables related to the different climatic area —the temperate climate subgroup (K¨ oppen–Geiger) (clim_zone_KG) and an alternative classification based on the contrast between low/high for air temperature and rainfall (clim_zone_TR4) (Table 1, SI.7)— did not improve the core model. However, the latter scored much better than the former in terms of both AIC and BIC and gave similar output to the core model. When comparing the outputs of model building with pheno_DH (Fig. 2) and crop_DH (SI.12) as core models, the results appeared remarkably similar, although some variations can be observed. A first general feature is that with crop_DH as the core, the contrast between AIC and BIC was much larger than in Fig. 2, as expected due to the larger number of parameters to be estimated in the crop model. Moreover, the Fig. 2. Screening of competitor models starting from the base model with distance (dist_DH) and the core model with crop phenology. The top panel compares the base model with distance only and the two core models extending the base model with phenology or crop. The black squares indicate AIC/BIC for the base model dist_DH. These values are adjusted downwards to not distort the readability of the graph. The other panels below start from the core model chosen and include crop, tree age, tree row and climate characteristics to test the research hypotheses (see Table 2) and to find out the optimal model format. If none of the models in a certain group scores better than the core model, the group is considered as not relevant. The models (y-axis) are ordered by BIC to reveal discrepancies with AIC easily. Additional information about the characteristics and performances of these competitor models are available in SI.11. A. Panozzo et al. Agricultural Systems 232 (2026) 104578 7 BIC scores for the crop model were in most cases worse than with the phenology core model. Yet, the BIC scores of the tree age-models attained similar values, confirming the impact of tree age. In contrast, the tree row parameters became less important, confirming the analysis in Fig. 2. Also, including the continuous climate variables to the crop model did not imply an improvement from the BIC perspective, however, from the AIC perspective there seemed to be a large gain (SI.12), supporting the evidence of a relevant impact of climate. However, the dataset is not sufficiently large to rule out which climate format performed best, as depending on the perspective other models were classified as the best. 3.3. Model visualization: distance from the tree row (RQ1) and crop type (RQ2) According with the initial hypotheses (H2), a contrasting yield response between winter and summer crops was confirmed, with the first group being less negatively affected by tree presence than summer crops, as revealed in Fig. 3. Within winter crops, a different slope of RR ~ ln(D/H) regression lines was further revealed. In particular, barley was the least impaired by tree row proximity across all the investigated crops, revealing a nearly null effect of trees on grain yield, as reflected by the flat regression line. When examining at different D/H intervals, a significant barley RR decrease was observed only when D/H was lower than 0.25 (RR = − 0.16, [−0.28, −0.01]), with 15 % lower yield vs. C (SI.14). In wheat, the slope of the regression line indicated a greater yield decrease as the tree row approaches, with a significant negative effect until a D/H value between 0.5 and 1 (SI.13), i.e., ~ −33 % of grain yield from the tree row to D/H =0.5 and −16 % at the inter-row width interval between D/H =0.5 and D/H =1 (SI.14). Across the two summer crops, regression lines of the core model ‘RR ~ ln(D/H) * crop species’ were parallel, reflecting a similar pattern of yield decrease for maize and soybean when getting closer to the tree line. However, the largest yield loss, compared to full sun, was observed in soybean. This is shown by the lower intercept of the regression line, and significant yield reductions along the entire investigated inter-row, always exceeding 40 % yield loss vs. C and reaching −75 % at the closest D/H interval (D/H <0.25; SI.13 and SI.14). In maize, instead, significant yield losses compared to full sun were limited to the inter-row width interval between D/H =0 and D/H =0.5 (SI.13), with an average yield loss of −64 % and −37 % vs. C, respectively at D/H <0.25 and at 0.25–0.5 interrow width intervals. In general, ln(D/H) =0, i.e., when the distance from the tree row is equal to tree height, emerged as a relevant threshold. In fact, in the four crop species this is the distance at which RR approaches 0, meaning that grain yield is recovered and reaches a similar value as in the full sun control (Fig. 3). 3.4. Model visualization: tree and design characteristics (RQ3) Among tree and design characteristics, only tree age contributed to a relevant improvement of the core models, with the two-class variable ‘age_C2’ showing the greatest drop in both AIC and BIC (Fig. 2, SI.11). The graphical visualization of the age_C2 model highlights that the negative impact of the tree row on grain yield increases with tree age, with the ranking from young (≤8 years old: lower slope) to old trees (> 8 years old: higher slope) maintained across all crop species (Fig. 4). In winter crops, different patterns were revealed depending on tree age: in wheat, only the older trees caused a decreasing RR with tree proximity, while young trees had no effect; in barley, in contrast, old trees had a nearly null effect, whereas young trees had beneficial (RR >1) effects, increasing productivity as they approached the tree row. In summer crops, the grain yield was always negatively affected, regardless of the tree age, and was much more severely impacted than in winter crops with old trees. 3.5. Model visualization: climatic characteristics (RQ4) Given the interest in the impact of all quantitative climatic variables, all models with climate moderators (original scale) as growing season means (G format) were plotted for the crop_DH core model. As expected, pet and temp revealed patterns very similar to one another within the same crop species, and were opposite to rain and arid (Fig. 5). When considering pet and temp, the trend of RR decreases when getting close to the tree row was the lowest with high values of both parameters, and Fig. 3. Graphical visualization of the core model of the dataset: grain yield variation (RR yield) in response to the distance to the tree row, expressed as the ln of the ratio of the distance from the tree row (D) to the height of the trees (H) (ln(D/H), according to crop phenology and species. Black lines represent the slope estimate of the core model ‘RR ~ ln(D/H) * crop phenology’; colored lines represent the slope estimate of the core model ‘RR ~ ln(D/H) * crop species’, respectively for barley (red), winter wheat (green), maize (blue) and soybean (purple). Points represent individual observations, and grey bands represent the confidence bands (95 % confidence intervals). A positive response ratio (RR >0) means that the presence of the tree row has a positive effect on the crop yield, i.e., higher yield vs. C, while a negative response ratio (RR <0) means a negative effect of trees on crop yield, and effect sizes close to zero (RR ≈0) indicate little or no effect. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) A. Panozzo et al. Agricultural Systems 232 (2026) 104578 8 Fig. 4. Graphical visualization of the full model including tree age: grain yield variation (RR yield) in response to ln(D/H), according to tree age, expressed as categorical variable, and crop species. Lines represent the slope estimate of the full model ‘RR ~ ln(D/H) * (crop species +tree age_C)’, respectively for trees ≤8 years old (red) and >8 years old (blue). Points represent individual observations, and red or blue bands represent the confidence bands of the regression lines (95 % confidence intervals). A positive response ratio (RR >0) means that the presence of the tree row has a positive effect on the crop yield, i.e., higher yield vs. C, while a negative response ratio (RR <0) means a negative effect of trees on crop yield, and effect sizes close to zero (RR ≈0) indicate little or no effect. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 5. Graphical visualization of four of the full models including climatic characteristics: grain yield variation (RR yield) in response to ln(D/H), according to four climate variables and crop species. The four models included potential evapotranspiration (pet), rainfall (rain) and air temperature (air) as a categoric variable, considering the growing-to-harvesting (G) mean. Lines represent the slope estimate of the full model ‘RR ~ ln(D/H) * (crop species +climate)’, respectively for low (red), middle (green) and high (blue) categories (defined in SI.7 and SI.8) of the climatic variables considered. A positive response ratio (RR >0) means that the presence of the tree row has a positive effect on the crop yield, i.e., higher yield vs. C, while a negative response ratio (RR <0) means a negative effect of trees on crop yield, and effect sizes close to zero (RR ≈0) indicate little or no effect. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) A. Panozzo et al. Agricultural Systems 232 (2026) 104578 9 Van Vooren, L., Bert, R., Steven, B., Victoria, N., Paul, P., Kris, V., 2017. Ecosystem service delivery of Agri-environment measures: a synthesis for hedgerows and grass strips on arable land. Agric. Ecosyst. Environ. 244, 32–51. https://doi.org/10.1016/ j.agee.2017.04.015. Viechtbauer, W., 2010. Conducting meta-analyses in R with the metafor package. J. Stat. Softw. 36, 1–48. https://doi.org/10.18637/jss.v036.i03. Visscher, A.M., Vanek, S., Meza, K., Wellstein, C., Zerbe, S., Ccanto, R., Fonte, S.J., 2023. Tree-based land uses enhance the provision of ecosystem services in agricultural landscapes of the Peruvian highlands. Agric. Ecosyst. Environ. 342, 108213. https:// doi.org/10.1016/j.agee.2022.108213. Weninger, T., Scheper, S., Lack´ oov´ a, L., Kitzler, B., Gartner, K., King, N.W., Michel, K., 2021. Ecosystem services of tree windbreaks in rural landscapes—a systematic review. Environ. Res. Lett. 16, 103002. https://doi.org/10.1088/1748-9326/ ac1d0d. Wolz, K.J., DeLucia, E.H., 2018. Alley cropping: global patterns of species composition and function. Agric. Ecosyst. Environ. 252, 61–68. https://doi.org/10.1016/j. agee.2017.10.005. Wu, Y.S., Feng, Y.A.N.G., Gong, W.Z., Ahmed, S., Fan, Y.F., Wu, X.L., Yang, W.Y., 2017. Shade adaptive response and yield analysis of different soybean genotypes in relay intercropping systems. J. Integr. Agric. 16, 1331–1340. https://doi.org/10.1016/ S2095-3119(16)61525-3. Xu, Y., Lehmann, L.M., García de Jal´ on, S., Ghaley, B.B., 2019. Assessment of productivity and economic viability of combined food and energy (CFE) production system in Denmark. Energies 12, 166. https://doi.org/10.3390/en12010166. Zomer, R.J., Trabucco, A., Bossio, D.A., van Straaten, O., Verchot, L.V., 2008. Climate change mitigation: a spatial analysis of global land suitability for clean development mechanism afforestation and reforestation. Agric. Ecosyst. Environ. 126, 67–80. https://doi.org/10.1016/j.agee.2008.01.014. A. Panozzo et al. Agricultural Systems 232 (2026) 104578 16