Rubber agroforestry in Thailand provides some biodiversity benefits without reducing yields

1School of Environmental Sciences, University of East Anglia, Norwich, UK; 2Department of Biology, University of York, York, UK; 3Department of Animal and Plant Sciences, University of Sheffield, Sheffield, UK; 4Department of Science, Pitchalai Preparatory School, Songkhla, Thailand; 5Department of Biology, Faculty of Science, Prince of Songkla University, Songkhla, Thailand; 6Department of Zoology, University of Cambridge, Cambridge, UK; 7CIRAD, UPR Systèmes de pérennes, Hevea Research Platform in Partnership (HRPP), Kasetsart University, Bangkok, Thailand; 8CIRAD, UPR Systèmes de pérennes, Univ Montpellier, Montpellier, France; 9CIRAD, UMR Innovation, Montpellier, France; 10CIRAD, INRA, Montpellier SupAgro, Montpellier, France and 11Faculty of Technology and Community Development, Thaksin University, Phatthalung, Thailand

There is a trade-off between rubber yields, and biodiversity and ecosystem function Drescher et al., 2016). At the two extremes, monocultural rubber yields are approximately double or triple those in 'jungle' rubber agroforests (Villamor et al., 2014;Warren-Thomas et al., 2015). However, yields in simpler agroforests containing only a few additional plant species, particularly those using clonal rubber varieties with even-aged trees, may not suffer this penalty. Evidence for nonlinear relationships between livelihoods and biodiversity in other systems Teuscher et al., 2015) opens the question: could rubber production become more biodiversity friendly without reducing yields?
In Thailand, the world's largest rubber producer, much lowland forest has been converted to rubber, 90% of which is grown by smallholders (Somboonsuke & Wettayaprasit, 2013). Intensification of rubber production to monocultures was incentivized via the Office of Rubber Replanting Aid Fund (now the Rubber Authority of Thailand, RAOT, Romyen, Sausue, & Charenjiratragul, 2018). Over 85% of Thai rubber is now grown monoculturally using clonal planting material at standard planting densities (Simien & Penot, 2011). However, <15% is grown in agroforests: most commonly in relatively intensive systems combining modern rubber cultivation methods with additional crops (Simien & Penot, 2011;Stroesser, Penot, Michel, Tongkaemkaew, & Chambon, 2018), rather than jungle-type systems. In response to price fluctuations, incentives have recently shifted towards rubber conversion to fruit or oil palm, and a policy promoting diversification through rubber agroforestry was approved in 2014 (ORRAF, 2006;Stroesser et al., 2018). However, rubber remains the dominant choice for smallholders in southern Thailand due to the well-established supply chain and local processing facilities, while further north, conversion to oil palm will be restricted by its intolerance of long dry seasons and cold (Warren-Thomas et al., 2018). Although we are not aware of any biodiversity-specific policies relating to rubber agroforestry, farmer motivations to undertake agroforestry include: receiving intercropping planting material from RAOT, anticipation of increased incomes and awareness of environmental benefits for soils, water and climate (Romyen et al., 2018). Rubber remains the main cash crop in such systems (Stroesser et al., 2018). To our knowledge, no study has assessed the biodiversity value of intensive rubber agroforests relative to monocultures, while research assessing measures to improve biodiversity in rubber monocultures is also scarce.
Structurally diverse jungle rubber in Sumatra supports frugivorous birds absent from monocultures, seven to 13 additional threatened species and eight additional forest specialists (Beukema et al., 2007;Prabowo et al., 2016). Non-rubber trees also provide additional canopy cover that could influence species composition, particularly of ectotherms (Koh, 2007;Wanger et al., 2010). Increasing height and density of the plantation understorey could also increase species richness, as shown for birds and butterflies (Aratrakorn, Thunhikorn, habitat extent in the landscape. Conservation priority and forest-dependent birds were not supported within rubber.

Synthesis and applications.
Rubber agroforestry using clonal varieties provides modest biodiversity benefits relative to monocultures, without compromising yields. Agroforests may also generate ecosystem service and livelihood benefits.
Management of monocultural rubber production to increase inter-row vegetation height and complexity may further benefit biodiversity. However, biodiversity losses from encroachment of rubber onto forests will not be offset by rubber agroforestry or rubber plot management. This evidence is important for developing guidelines around biodiversity-friendly rubber and sustainable supply chains, and for farmers interested in diversifying rubber production. & Donald, 2006;Azhar et al., 2011;Barbosa Cambui, Nogueira de Vasconcelos, Mariano-Neto, Felipe Viana, & Zikán Cardoso, 2017).
In Thailand, monocultures with vegetated understories had greater bird species richness than those without (Aratrakorn et al., 2006).
In Brazil, rubber plantations with 10-to 20-year-old understorey supported four additional butterfly species, and were more similar in composition to forest fragments than intensive rubber (Barbosa Cambui et al., 2017). Finally, biodiversity responds to land use at multiple spatial scales, and on-farm biodiversity is influenced by the wider landscape (Tscharntke, Klein, Kruess, Steffan-Dewenter, & Thies, 2005). In rubber monoculture-dominated landscapes of Southwest China, bird species richness was greater in landscapes with more forest cover (Zhang, Chang, & Quan, 2017), and richness may also increase in response to greater land-use diversity.
We investigated whether Thai agroforests can offer benefits for both rubber yields and biodiversity, using three contrasting taxa (birds, fruit-feeding butterflies and reptiles), selected because they are taxonomically resolved, relatively well studied in other agroforestry systems and likely to respond to different aspects of management.
First, we examined whether yield, species richness and species composition differed between agroforests and monocultures. Second, whether richness and composition varied in response to understorey vegetation, and the types and densities of non-rubber crops and trees. Third, how richness and composition were influenced by landscape composition, and potential interactions with plot management.

| Study region
The study was conducted in southern Thailand ( Figure S2), where lowland landscapes are smallholder rubber dominated (farm sizes <1 ha to several ha), with smaller areas of oil palm, orchards, rice and forest. The largest forest fragments were ~320 ha of karst hilltop forest and 400 ha of secondary forest; others were small (~4 ha) and degraded. Three protected forests cover mostly upland areas (100-1,350 m asl). Rain is frequent May-December, with a dry season January-March (Phommexay, Satasook, Bates, Pearch, & Bumrungsri, 2011). Biodiversity data were collected March-June 2016, during unusually low rainfall and high temperatures linked to a strong El Niño-Southern Oscillation ( Figure S3).

| Sampling plots
Biodiversity data were collected from 64 rubber 'plots' in Songkhla and Phattalung provinces. Plots were defined as even-aged units ≥1 ha (mean 1.92 ha, range 1.0-7.6 ha) at least 100 × 100 m in dimensions (measured using GPS and laser rangefinder), termed the 'biodiversity dataset'. Plots were categorized as monocultures (MO, n = 25) or agroforests (AF, n = 39). Plots were chosen by first locating an appropriately sized agroforestry plot with the help of local agroforestry groups, because agroforests are relatively rare in the landscape. Almost all agroforestry plots above the minimum size requirements were surveyed, and we then co-located monoculture plots in the vicinity of each agroforest. Clonal rubber trees were usually planted at 3 m intervals in rows 7 m apart in both systems (476 stems per ha, Phommexay et al., 2011). Agroforests contained cultivated non-rubber tree, shrub or herbaceous plants, or naturally regenerated wild trees, mostly in the inter-row (mean 162 ± 207 SD stems per ha). Complexity ranged from one or two cultivated species (Figure 1b), to multiple native tree species (Figure 1d). Three monocultures contained non-rubber crops at densities too low to be considered agroforests (e.g. a single bamboo clump, 17 ± 51 stems per ha). Latitude of biodiversity dataset plots ranged from 6.641374 Agroforests and monocultures were represented in each block wherever possible, but five blocks contained only agroforests.
Blocks were clumped within five 'districts' (not administrative districts, up to 127 km apart, Figure S2), with plots <9.00 km apart within a district (Figure 2b).

| Biodiversity data
Biodiversity data were collected from up to nine agroforestry and monocultural plots each day, to control for potential weather and seasonal effects. The number of rainy days during sampling was recorded based on field or overnight observations (within 10 km of plot). Strong El Niño conditions during sampling could have influenced sampling by decreasing the suitability of relatively open habitats due to increased heat and decreased humidity (such as monocultures with lower canopy cover), or by forcing some species into cooler, wetter, refugia (e.g. stream beds).
Birds were surveyed by one experienced ornithologist (L.N.) on three consecutive mornings (06:00-09:30 hr) in 10-min point counts at each plot centre, alternating the order in which plots were visited (Gilroy, Woodcock, et al., 2014b). Birds within a 50-m radius (excluding fly-throughs) were identified to species using sight or sound, and number of detections recorded. A 2-week pilot phase was used to familiarize L.N. with local species and the 50 m sampling radius. In the smallest plots (100 m × 100 m), the sampling radius reached the plot boundaries, and during the pilot phase, these plots were used to check the location of calling birds: this experience was used to estimate the 50-m radius in larger plots. Sound recordings were made of each count (Olympus LS-11 Linear Recorder), and sounds checked against existing recordings (Xeno-Canto Foundation, 2017). Species' habitat associations (del Hoyo, Elliott, Sargatal, Christie, & Juana, 2017) and conservation status (IUCN, 2016) were defined.
Five traps were set per plot, one at the centre and four 50 m away in cardinal directions, baited with two tablespoons of fermented banana mixture (750 ml ripe mashed bananas, one teaspoon quick action yeast, two tablespoons sugar, 1 tablespoon rum, fermented for 48 hr). Traps were set with the base 1.5 m from the ground. In plantation systems with simple canopy structures, this effectively samples butterflies from all strata (Barlow, Overal, Araujo, Gardner, & Peres, 2007). Traps were set on the first day and checked on four subsequent afternoons (13:00-18:00 hr), replacing bait each day and discarding old bait away from the plot. Without trap loss or damage, this gave 20 trap-days per plot. Trapped butterflies were photographed (Canon 700D D-SLR, 105 mm macro lens) and marked before release. No individuals were re-trapped across sampling plots. All analyses use presence-absence data, except where comparing total catch (as a proxy for total abundance across all species) between AF and MO, as the relative abundance of species is unlikely to be reflected in relative catches (Hughes, Daily, & Ehrlich, 1998). Individuals were identified to subspecies following (Corbet & Pendlebury, 1992;Ek-Amnuay, 2012)

| Rubber plot management
Rubber plot management (and resulting vegetation structure) data were collected from each plot in the biodiversity dataset. Cultivated plant species names were recorded during a farmer interview.
Vegetation structure was quantified through field measurements ( Figure S4 and text) giving plot-level values for herb height (cm), canopy cover (%), small stem density (ha −1 ), non-rubber tree stem density (ha −1 ), fruit tree stem density (ha −1 ) and the number of agroforestry species.

| Landscape composition
Landscape composition across each block was quantified by recording land covers with a GPS during semi-exhaustive surveys on foot (access and terrain permitting) resulting in a mean 139 ± 43 SD GPS points per block. The points were used to define land use at 100-m intervals along the block perimeter, once within each sampled plot and once in the adjacent management units in each cardinal direction, giving 39 points per block ( Figure 2). Where sampled plots were adjacent (as in Figure 2), land cover of the next management unit was recorded. Where only two plots were sampled within a block, land use was recorded in one additional management unit and its neighbours.
Land cover was recorded as: rubber agroforestry (AF), rubber monoculture (MO), immature rubber, bare ground, scrub, urban (village, road or town), natural forest, fruit orchard, home garden, cassava, oil palm, rice paddy, timber plantation or coconut grove. Using Google Earth, we mapped streams and rivers, and calculated riparian feature length per block. These data were summarized into six variables: percentage of the points (hereafter termed extents) in rubber (total AF and MO), open habitats (total of immature rubber, bare ground, cassava or rice paddy) and natural forest, the ratio of AF to MO, Shannon-Wiener diversity index of land uses (using point-frequency data) and riparian length.
This point-based sampling approach was validated by comparison with area-based measures of landscape composition, extracted from manual mapping of management units using Google Earth imagery for a subset of 10 blocks, informed by all available GPS points for each block. Automated land-use classification was impossible due to the identical canopy characteristics of monocultures and agroforests. The two datasets correlated strongly (Pearson correlation, r ≥ .77, p ≤ .05, tested separately for each land cover type; Figure S5) supporting the validity of the point-sampling approach to quantify land use.

| Rubber yields
Rubber yield estimates were collected in 2016, via questionnaires with farmers, for a separate set of 47 agroforest and 34 monocultural rubber plots in Phattalung province with considerable spatial overlap to the 'biodiversity dataset' ( Figure S2). These data were collected as part of a wider study of farmer livelihoods and rubber agroforestry, which also included data on agroforestry species composition, and is termed the 'yield dataset' (full methods in Stroesser et al., 2018). Annual rubber yields were either reported for latex (coagulated, reported as dry weight) or as 'cup-lump' values per plot (converted to dry weight assuming dry rubber content of 65%) to give yields in t ha −1 year −1 .
Plot characteristics with the potential to affect rubber yield were compared between the yield and biodiversity datasets using generalized linear models and Mann-Whitney U tests. Rubber yields in Southeast Asia are affected by elevation, cold and drought (Ahrends et al., 2015). Mean latitude of yield plots was 7.473321°N (SD 0.199305°), and mean elevation (proxy for temperature and drought There were no differences in cultivated plant species richness per plot, or the stem densities of rubber, fruit and timber trees respectively (stems per ha; Figures S6 and S7). The agroforestry systems are thus similar in both datasets, and any effects of structural complexity or intercropping density on yields and biodiversity are likely to be highly comparable.
Soil type may affect rubber yield. Soils for the region are defined as Acrisols (Food & Agriculture Organization of the United Nations, 2012), but according to a national-level soil survey, the plots fell within six soil types, dominated by Udults (freely drained, humus poor, light soils; Department of Land Development, 2002;USDA, 1999). A greater proportion of yield plots had marginal 'skeletal' soil types than biodiversity plots (53% vs. 22%; Table S1). It is possible that soils in yield plots are generally poorer quality than biodiversity plots. However, within the yield dataset, 47% of plots were skeletal, of which 64% were agroforest, suggesting yield variation due to soil type is likely to be present within both agroforests and monocultures. Support from RAOT during rubber planting/re-planting means farmers tend to follow standard practices for immature rubber, but experienced rubber farmers in southern Thailand have diverse ongoing management practices (e.g. fertilization rates, tapping intensity) once trees are mature (Besson, 2002;Chambon, Promkhambut, Duangta, Lesturgez, & Sainte-Beuve, 2017). However, we cannot conceive any reason why management practices would differ systematically between the yield and biodiversity datasets.

| Analysis
All analyses were conducted in r separately for each taxon (R Core Team, 2017). Estimated rubber yields (t ha −1 year −1 ) were compared between AF and MO yield plots using a general linear model including plot type, soil type (loamy/clayey, skeletal or slope complex, Table   S1) and their interactions. The AIC c of this model was compared to a null model; where ∆AIC c of the alternative model was ≥2 lower than the null, it was considered to be informative. Power analysis of the result was conducted using package pwr (Champely, 2018), to estimate the effect size (difference in yields) detectable with our sample.
Sampling completeness within each plot type (AF or MO) was calculated as the percentage of the estimated asymptotic species richness that was observed, based on four estimators (Jack1, Jack2, Bootstrap and Mmean). These were compared between AF and MO using a Mann-Whitney U test. Cumulative species richness was compared between AF and MO using sample-based rarefaction extrapolated to the largest sample size (n = 39 for AF) using the iNEXT package (Chao & Colwell, 2014).
Species richness and number of detections or catches per plot (proxies for abundance) were compared between AF and MO plots using generalized linear mixed models (GLMMs), fitted using maximum likelihood, with a Poisson distribution and log link function, using the lmE4 package (Bates, Maechler, Bolker, & Walker, 2015).
Block was included as an intercept-only random effect to account for the nested sampling design, with district and rainfall index additional intercept-only random effects for butterfly models, as these influenced butterfly species richness (Figures S8-S10). The AIC c of the full models was compared to a null model containing only random effects (Burnham & Anderson, 2002). Systematic spatial autocorrelation of model residuals, examined using a Monte-Carlo permutation test for Moran's I using package spdEp with 1,000 iterations (Bivand & Wong, 2018), was not found for any model. Residuals were tested for overdispersion, but theta (Pearson residuals/residual degrees of freedom) was <1 in all cases (Burnham & Anderson, 2002). Species composition response to agroforestry was investigated with a partial redundancy analysis (RDA) with Block as a conditional effect, using detection data for birds and reptiles (maximum individuals, or detections, recorded on one sampling day, divided by the variance of each species; Oksanen et al., 2017) and presence-absence data for butterflies (Hughes et al., 1998).
Species richness response to plot management was then investigated using multimodel inference (Burnham & Anderson, 2002). A global GLMM (containing all six vegetation structure variables, centred and standardized to zero mean and 0.5 SD so that effect sizes were on comparable scales; Grueber, Nakagawa, Laws, & Jamieson, 2011), a null (intercept-only) model containing the random effects and a candidate model set (comprising all possible model subsets with four or fewer variables, ensuring at least 15 observations for each candidate variable) were generated using the mumiN package (Bartoń, 2016;Grueber et al., 2011). The 57 candidate models were ranked according to AIC c and weights using the bbmlE package (Bolker & R Development Core Team, 2017). Models with a cumulative weight ≥95% were averaged by the full (zero) averaging method (Burnham & Anderson, 2002;Grueber et al., 2011), using the mumiN package (Bartoń, 2016). If the 95% confidence intervals of the averaged parameter estimate excluded zero, it was considered influential (Grueber et al., 2011). Influential rubber plot management variables were included in the next step of analysis.
The same approach was used to generate a final set of candidate models for species richness responses. Explanatory variables were plot type (AF vs. MO), influential rubber plot management variables, landscape composition variables and two interaction terms (plot type and AF:MO ratio, plot type and natural forest). Model averaging was applied to the final model set. To estimate effect sizes and investigate interactions, species richness was predicted from the final averaged models for each taxon using the mumiN package, holding all other continuous explanatory variables at the mean and including mean levels of random effects (Bartoń, 2016 (Oksanen et al., 2017;Wagner, 2004), was not found for any taxon. Automatic backward model selection was performed on the global model (9,999 permutations) with a threshold of p > .05 (a 'pseudo-F' test statistic: ratio of constrained and unconstrained total inertia in the RDA, divided by their respective ranks) to drop terms from the model (Legendre, Oksanen, & Braak, 2011;Oksanen et al., 2017). The significance of each term in the final model was examined using the same method.

| Rubber yields
Rubber yields did not appear to differ between AF (mean 1.34 t ha −1 year −1 ± 0.61 SD) and MO plots (1.51 ± 0.54 t ha −1 year −1 ; Figure 3), and soil type also had no effect ( Figure S11). With the sample size available (n = 77), power analysis showed that we would have been able to detect a difference in yield of 0.24 t ha −1 year −1 (17% change relative to mean yields of all plots) given a power of 0.80. We would have required a sample of 117 plots to detect a difference of 0.14 t ha −1 year −1 (i.e. a 10% change). It is therefore possible that we have failed to detect a real, but small, difference between the two systems. Yield variability of AF plots was slightly higher (SD of 0.61, range 0.05-2.78 t ha −1 year −1 ) than in MO plots (SD of 0.54, range 0.12-2.42 t ha −1 year −1 ), but both systems show considerable yield variability, suggesting that neither option necessarily reduces farmer risk in terms of yield variation.

| Rubber plot management and structure
AF plots in the biodiversity dataset had greater fruit and timber tree species richness (mean 6.6 ± 0.6 SE non-rubber tree species, 0.2 ± 0.1 in MO), greater density of timber, fruit and wild trees (162 ± 33 stems per ha, 17 ± 10 in MO), smaller rubber basal area (16.8 ± 1.0 m 2 /ha, 21.5 ± 1.7 in MO) and greater canopy density (76% ± 0.8, 71% ± 1.4 in MO), small stem density (2,552 ± 522, 601 ± 224 in MO) and understorey density (index 6.6 ± 0.7, 4.3 ± 0.7 in MO; Figure S12). However, mean herb height was similar between AF and MO, because understorey management varied among plots of both types. Three tree species recorded in agroforests are globally threatened, and a total of 37 plant species were recorded across all plots grown for a variety of uses, including timber, fruit and resin (Table S2).

| Species richness and composition of birds, butterflies and reptiles, in agroforests compared to monocultures
In total, 1,204 registrations of 69 bird species, 544 detections of 17 reptile species, and 809 individuals of 44 fruit-feeding butterfly species (excluding female Mycalesis) were recorded. Plot-level detections of each species are given in Table S3. Species richness estimators showed that ≥74% of species were detected, and sampling completeness did not differ between AF and MO ( Figure S13).
Agroforestry supported greater butterfly species richness, both at habitat level (cumulative species richness across all plots; Figure 4c) and plot level (Figure 4f). AF plots contained 6.1 ± 3.4 (mean ± SD) species, while MO contained only 3.6 ± 3.5. However, bird and reptile species richness did not differ between AF and MO, whether at habitat (Figure 4a, b) or plot level (Figure 4d, e). MO plots contained 11.9 ± 3.0 bird and 4.6 ± 1.7 reptile species, and AF 12.8 ± 3.1 bird and 4.3 ± 1.4 reptile species. Detections of birds and reptiles did not differ between AF and MO, but total butterfly catches were greater in AF ( Figure S14). Species composition of all three taxa was unaffected by agroforestry (Table S4).

| Species richness and composition of birds, butterflies and reptiles, in response to agroforestry, rubber plot management and landscape composition
Bird, butterfly and reptile species richness each responded to different combinations of plot type (AF or MO), plot management and landscape composition variables. Plot-level bird richness increased with greater herb height (Figure 5a), but no other variables affected richness (Figure 5d). Predictions from the final averaged model showed that increasing herb height from 37 cm (25% quantile) to 98 cm (maximum), increased mean bird richness per plot from 11.8 (weighted across AF and MO) to 14.8 (Figure 6a).
Butterfly richness was greater in AF than MO (Figure 4c, f), and was affected by an interaction between plot type (AF or MO) and forest extent (Figure 5f), although we note that few plots were in landscapes with relatively high forest extent. Within AF, butterfly richness increased with natural forest extent, but not in MO F I G U R E 3 Rubber yield of agroforest (AF) and monoculture (MO) plots. Boxes = 25% and 75% quartiles; thick lines = median; notches = 95% CI; diamonds = mean; whiskers = 1.5× interquartile range. ∆AIC c is for the general linear model incorporating plot type (AF or MO), soil type ( Figure S11) and their interactions, relative to the null model (i.e. the null model has lower AIC c ). Yield is unaffected by plot type  (Figure 6b). An agroforest in a block containing no forest (25% quantile) was predicted to contain 4.8 species, which increased to 11.2 species in a block containing 51% natural forest (maximum). In Bird species composition showed differing patterns to those of richness: herb height affected composition (as for richness), but non-rubber tree density and the surrounding extent of rubber, forest and open habitat in the block were also influential ( Figure 7a; Table S5). The RDA model explained 16% of total inertia (pseu- Butterfly species composition altered in response to forest extent, the interaction between plot type and ratio of AF:MO in block and canopy cover (Figure 7c). This RDA model explained 15% of total inertia (pseudo-F = 1.95, p < .001; Table S5). To explore the interaction between plot type and the ratio of AF:MO, the RDA was rerun separately for AF and MO ( Figure S16), showing that composition only responded to the ratio of AF:MO within monocultures, with some species strongly associated with landscapes containing a greater proportion of agroforestry. Species composition in AF was also affected by forest extent, but not in monocultures, mirroring patterns in species richness.
In contrast to both birds and butterflies, reptile species composition only changed in response to canopy cover (positive response to canopy openness by some species) and open habitat extent F I G U R E 4 Species richness of three taxa in rubber agroforest (AF, circles) and monoculture (MO, triangles), showing rarefaction and plotlevel richness. Upper panels (a-c) show sample-based rarefaction and extrapolation of detection data for birds and reptiles, and presenceabsence data for butterflies (excluding Mycalesis females). Dashed lines (a-c) show extrapolation of MO sample (n = 25) to AF sample size (n = 39), rescaled to number of individuals/detections (birds/reptiles) or plots (butterflies); grey shading shows 95% CI. Lower panels (d-f) show mean and 95% CI (whiskers) of species richness per plot; ∆AIC c is for the alternative model incorporating plot type, relative to the null. The alternative model was no better than the null for birds and reptiles, but was best for butterflies (asterisks indicate support for difference between plot types)  Table S5). This model explained 14% of inertia (pseudo-F = 1.79, p = .005).
Only two bird species of conservation concern were recorded, both as singletons within AF plots (IUCN NT), while one forest-dependent species was recorded in both AF and MO (Table S3). Fifteen bird species were open-habitat specialists, but the majority were small generalist insectivores. The five reptile species with IUCN assessments were all LC (nine unassessed), and no forest specialists (Table S3) were associated with agroforestry, plot management variables or landscape composition.

F I G U R E 5
Species richness response of three taxa to plot management and landscape context. Upper panels (a-c) show response to management, lower panels (d-f) to plot type (MO vs. AF reference), influential structural/management variables (*in a-c) and landscape composition. Parameter estimates (predicted change in species richness with a one-unit change of the standardized predictor variable) are from full model averaging across the 95% confidence set of models, each containing a maximum of four predictors. The central bar line shows the mean of parameter estimate and the bar encompasses 95% CI; those excluding zero are considered influential, marked "*". Relative variable importance (proportion of models in set containing each predictor) shown above each bar. Structural/management variables: Can_Cov = canopy cover; Fru_stha = fruit trees stem density; Hrb_h = herb height; n_AF_spp = n agroforestry species; Non_rub_stha = non-rubber tree stem density; Sml_stha = small stem density. Landscape composition variables: AF_ratio = agroforestmonoculture ratio; Lduse_Shannon = Shannon diversity index of land covers; NF_prop = natural forest extent; Rub_prop = rubber extent; Stream = riparian feature length Future rubber demand is likely to drive expansion and/or intensification of plantations, driving biodiversity and ecosystem loss (Warren- Thomas et al., 2015). This generates a need for biodiversity-friendly but high-yielding cultivation methods. Clonal rubber agroforests in southern Thailand retained yields and benefitted butterfly species richness, while plot management (understorey plant height, nonrubber tree density) and landscape composition (proportion of rubber cultivation in agroforestry, natural forest extent) also increased species richness and altered composition of butterflies and birds. Our findings relating to forest extent must be interpreted in the context of a relatively small number of high-forest-cover landscapes, and further evidence for the importance of forest fragment cover on biodiversity in rubber is needed. We found no evidence that agroforests supported forest-dependent or conservation-priority species, but our findings suggest there is room for improvement of biodiversity within rubber plantations, without negatively impacting rubber yields.
Our finding that butterfly richness was greater in agroforests, and increased further with increasing forest extent (reflected in compositional change), suggests that agroforests may act as transitional habitat for fruit-feeding butterflies that can move several kilometres in a lifetime (Marchant et al., 2015). Contrastingly, monocultures are relatively impermeable for butterflies (Scriven, Beale, Benedick, & Hill, 2017). This corroborates findings of similar butterfly composition between forest fragments and structurally complex F I G U R E 6 Predicted species richness response to variables influential in averaged models of rubber plot management and landscape composition. Panels show: (a) bird response to herb height (no interaction with plot type; line shows predicted effect in both plot types) and b) butterfly response to natural forest extent, and interaction with plot type (dark grey line = AF, light grey line = MO). Original data points shown (black circle = AF plot, grey triangle = MO plot). Lines fitted to predicted species richness values (points not shown) with a linear model; CI not plotted as SE cannot be reliably computed for mixed effects models rubber in Brazil (Barbosa Cambui et al., 2017). These patterns could be driven by increased food plant availability, or microclimatic similarities between forest and agroforests (affecting flight ability or thermal tolerance; Koh, 2007 Unlike butterflies, bird richness was unaffected by agroforestry, but altered in agroforests with greater non-rubber tree density, and was similarly, but independently, influenced by greater natural forest extent around plots of either type. A similar conclusion was drawn for cacao, where distance to forest and shade tree density had independent but similar effects on bird richness and composition (Clough et al., 2009), and corroborates findings of altered bird diversity and composition in response to increased tree diversity and/or forest area in monocultural rubber (Zhang et al., 2017) and oil palm (Gilroy, Prescott, et al., 2014a;Teuscher et al., 2015).
Birds and butterflies using agroforests and forest fragments may differ functionally from forest-dependent or agriculture-tolerant species (Koh, 2007;Sekercioglu, 2012). Agroforests and forest fragments may therefore be refugia for species that cannot utilize monocultures (Bhagwat, Willis, Birks, & Whittaker, 2008), and their presence may increase landscape beta diversity (Faria, Paciencia, Dixo, Laps, & Baumgarten, 2007) and connectivity (Bhagwat et al., 2008). However, this does not mean that agroforests support forest-dependent or threatened species. The lack of such species in rubber agroforests in our study echoes findings in polycultural oil palm (Azhar et al., 2011(Azhar et al., , 2015, but contrasts with findings from industrial tree plantations of Albizia and Acacia (Sheldon, Styring, & Hosner, 2010), highlighting the importance of gathering evidence on biodiversity responses to different plantation systems. Together, this suggests that conservation of contiguous forest (Barbosa Cambui et al., 2017;Edwards et al., 2010) or complex jungle rubber (Prabowo et al., 2016) remains essential for forest-dependent and conservation priority birds and butterflies.
In addition to modest biodiversity benefits, rubber agroforests could provide additional ecosystem functions and services.
Integration of native trees into rubber improves water infiltration, improving and stabilizing soil (Langenberger, 2017), while reducing herbicide decreases runoff, soil erosion and loss of total organic soil carbon (Liu et al., 2016). Maintenance of herbaceous vegetation in rubber could be a simple management tool to both modestly increase bird diversity (by up to three species) and protect soils. This could also improve soil macrofauna diversity and abundance, as found in oil palm (Ashraf et al., 2018;Ashton-Butt et al., 2018). Evidence from oil palm suggests that reduced herbicide use does not reduce soil fertility (Ashton-Butt et al., 2018), but the potential effects of understorey plant competition for water and nutrients on rubber yields (Langenberger, Cadisch, Martin, Min, & Waibel, 2016), and the effect of understorey plant diversity on yields and birds, warrant further investigation. Further ecosystem services provided by more complex rubber agroforests could include enhanced pest control, pollination and decomposition (Maas et al., 2016;Tscharntke et al., 2005), climatic stability and carbon storage   (Wanger et al., 2010). Other studies assessing reptile response to agroforest complexity (Deheuvels et al., 2014), and forest cover around plantations (Faria et al., 2007;Gilroy, Prescott, et al., 2014a), have similarly failed to find effects. This may be because relevant variables were not measured, such as log piles, temperature or leaf litter (Wanger et al., 2010). Alternatively, reptile detectability could have been lower in agroforests due to increased habitat complexity. Similarly, arboreal reptiles (e.g. Draco. spp), that may be most likely to respond to agroforestry, are more difficult to detect than ground-dwelling species, meaning a real effect may have been missed.
Our finding of similar reported rubber yields in agroforests and monocultures is supported by existing empirical evidence that rubber tree growth is unaffected by inter-planting (Wibawa, Joshi, Noordwijk, & Penot, 2006). Our data were collected under strong El Niño conditions, resulting in higher temperatures and reduced rainfall than an average year, which may have reduced yields through drought stress. We cannot determine whether yields in agroforests or monocultures were relatively more resilient under these conditions, but this question is likely to become increasingly important as El Niño frequencies increase. Relationships between yields, agroforest structure and biodiversity warrant further direct investigation, including in relation to climatic changes. As rubber is a canopy tree, and secondary plants grow either below the canopy or in shared canopy space, the relationship between yields and agroforest complexity may differ from other agroforestry systems.
Many factors affect the sustainability of agroforestry for farmers aside from rubber yields: shade-grown crops yield less than unshaded, and labour constraints can make additional crop cultivation unfeasible Langenberger et al., 2016). However, in southern Thailand, simply structured high-yielding agroforestry not only improves incomes relative to monocultures, but also provides a social function by generating fruit crops to share within communities, and can provide as good or better return for labour (Romyen et al., 2018;Stroesser et al., 2018). Moreover, farmers with smaller plots were more likely to practice agroforestry, suggesting it provides additional benefits for poorer farmers (Romyen et al., 2018). Despite concerns about economic viability of agroforestry in some regions of Asia (Langenberger et al., 2016), appropriately designed rubber agroforestry seems to provide sustainable livelihood benefits while maintaining rubber yields in southern Thailand.

| Synthesis and applications
Overall, we found modest benefits for birds and butterflies in intensive-clonal rubber agroforests in southern Thailand, while yields were maintained. Increases in non-rubber tree density and understorey vegetation height enhanced bird diversity, while forest fragments were important at the landscape scale for birds and butterflies. Our findings echo those of studies in oil palm and cacao, and contribute further evidence to support sustainable intensification of tropical agriculture.
Small forest fragments may work synergistically with agroforestry for some taxa, such as butterflies, and should be retained, but conserving contiguous forest tracts is preferable for forest-dependent and threatened species (Edwards et al., 2010). Biodiversity losses from continued encroachment of rubber onto protected forests in Thailand (Aratrakorn et al., 2006)