Herb assemblage dynamics over seven years in different cocoa production systems

Both agronomic practices and spatial position can determine the assemblage of herbaceous species. We assess the dynamics and the contribution of these two aspects over time to the herb assemblages of different cocoa production systems. Braun-Blanquet surveys were performed over seven years in a long-term trial in Bolivia to compare different cocoa production systems: successional agroforestry (SA) with no external inputs, organic agroforestry (OA) and organic monoculture (OM), both including a leguminous perennial cover crop, conventional agroforestry (CA) and conventional monoculture (CM), where agrochemicals were applied. Using general linear models and multivariate analysis we found that assemblages were mainly driven by spatial position only at the beginning. After this, a very dynamic selection process related to the different management practices took place, which became more stable over the years. We observed a decline in species in both the CA and OA systems, due to the loss of heliophilous species and the low number of new species established in them. The OM presented the most conservative pattern, with the lowest number of new species and species lost, due to the presence of the cover crop. Both the most intensively managed system (CM) and the most diverse and least intensive one (SA) had the highest number of new species recorded over time, which led to highly specialized assemblages, with worldwide distributed and herbicide resistant species in the first case and secondary forest species, in the second. We conclude that the promotion of organic management and agroforestry systems, especially highly divers and successional agroforestry sytems would favour herb assemblages with high conservation value and prevent the establishment of globally distributed species.


Introduction
Herbs that inhabit agroecosystems are relevant for both agronomic production and biodiversity conservation. Farmers' experience indicates that greater effort in weed management relates to higher cocoa harvests (Aneani and Ofori-Frimpong 2013). In this line, there are species that can have negative impacts on production or diversity (Duguma et al. 2001;Afreen et al. 2017), and species that, managed in a certain way, can be useful for the cocoa plantation Abstract Both agronomic practices and spatial position can determine the assemblage of herbaceous species. We assess the dynamics and the contribution of these two aspects over time to the herb assemblages of different cocoa production systems. Braun-Blanquet surveys were performed over seven years in a long-term trial in Bolivia to compare different cocoa production systems: successional agroforestry (SA) with no external inputs, organic agroforestry (OA) and organic monoculture (OM), both including a leguminous perennial cover crop, conventional agroforestry (CA) and conventional monoculture (CM), where agrochemicals were applied. Using general linear models and multivariate analysis we found that assemblages were mainly driven by spatial position only at the beginning. After this, a very dynamic selection process related to the different management practices took place, which became more stable over the years. We observed a decline in species in both the CA and OA systems, due to the loss of heliophilous species and the low number of new species established in them. The OM presented 1 3 Vol:. (1234567890) (Chacon and Gliessman 1982). Herbs may host pathogenic organisms that are harmful to cocoa (Resende et al. 2000;Rich et al. 2009), as well as recruit their controllers-natural enemies- (Stenchly et al. 2012) or pollinators (Klein et al. 2003), or simply become part of the necessary resources for native animal species conservation (Gordon et al. 2007).
The herb richness and herb composition present in a crop field depend, on the one hand, on the ability of the species to "reach" the community or assemble through dispersion or regeneration processes, and, on the other, on the adaptations that allow the species to "thrive" in those specific ecological conditions, as well as on the availability of certain resources and the internal dynamics created by the crop (Booth and Swanton 2002). Agronomic management practices modify all these conditions and are a main selection factor of herbaceous species. For instance, in the case of cocoa and coffee, there are differences in the structure of the herb communities according to the levels of shade in the plantation, the structure of trees in the crop field (Soto-Pinto et al. 2002;Bobo et al. 2006;Cicuzza et al. 2011;Ramos et al. 2015), the use of herbicides (Konlan et al. 2019), and the application of fertilizers (Cicuzza et al. 2012).
Cocoa can be cultivated under a wide range of management intensity conditions, from full-sun monocultures with high external inputs, to very diverse and complex agroforestry systems with a large shade cover and no external inputs, or modified forests (Sambuichi et al. 2012;Somarriba and Lachenaud 2013;Zegada et al. 2020). In a previous study that compared the structure of herbal communities in different cacao production systems, including monocultures and agroforestry systems, as well as conventional and organic management, we found that the richness of strictly herbaceous species was not related to management intensity (Marconi and Armengot 2020). More importantly, we found huge differences in composition, indicating that the more intensive production systems, i.e., monocultures, particularly those under conventional management, welcome species that are worldwide distributed, while the species present in complex agroforestry systems tend to have a more restricted distribution.
The present study assesses the same communities but includes seven years of data and focuses in two different aspects. First, we aim to determine whether the previous pattern, observed with one year's data, varies over time. Herbaceous communities in crops are usually highly dynamic, adapted to the disturbances of management practices but also affected by spatial patterns (Colbach et al. 2000). For example, some invasive species are more common in cocoa plots near to roads, train tracks, and other exposed areas (Bakar 2004). Dispersion and regeneration processes might thus have an even stronger influence on the herbaceous community than agricultural practices, and the relative impact of all these factors can vary over time and depend on the type of production system.
Secondly, we aim to a better understanding of the dynamics of the herbaceous community, i.e., to know which and how many species persist or are new (immigrants) or lost in a certain area.
Considering, on the one hand, that herbal species in more intensively managed systems (e.g., monocultures under conventional management) have a wider geographical distribution range than those in less intensively managed systems (e.g., agroforestry under organic farming) (Marconi and Armengot 2020), and, on the other hand, that the invasiveness of a community increases with the number of disturbance events-such as weeding (Herben et al. 2004) and the fluctuation of resources (Davis et al. 2000)-, we hypothesise that more intensive production systems will be more dynamic, i.e., will experience more changes in the herbal community over time, with more local immigrations and extinctions than less intensively managed systems.

Study area
The study was carried out in Alto Beni, Bolivia (15° 27′ 30″ S, 67° 28′ 18″ W), at a long-term trial comparing different cacao production systems established in 2009 within the framework of the Syscom Programme (https:// syste ms-compa rison. fibl. org/ index. html). The trial site is located at an altitude of 365 m and receives an annual rainfall of 1519 ± 220 mm, 71 ± 6% of which accumulates between the months of October and March. The average annual temperature is 26 °C and the relative humidity is 78.9%.

Trial description
A secondary forest was cleared by cutting down of trees and chaqueo (controlled burning). Corn and Canavalia ensiformis were cultivated to test the homogeneity of the soil. At the end of 2008, 24 plots of 48 × 48 m were established and cultivated under five different production systems: conventional monoculture (CM), conventional agroforestry (CA), organic monoculture (OM), organic agroforestry (OA), and successional agroforestry (SA). In CM herbicides such as glyphosate or paraquat were applied four to five times per year and as granular fertilizers were applied around the cocoa trees twice per year; eventually, brush cutters were also used to control weeds. In CA the same agrochemicals were appliedd but at lower doses and frequency. Shade trees were planted, mainly legume trees (8 × 8 m spacing) such as "pacay" (Inga spp) and ceibo (Erythrina poeppigiana and E. fusca); also, in lower densities: fruit species such as avocado (Persea americana) and peach palm (Bactris gassipaes), and timber species such as mahogany (Swetenia macrophylla) and paquÍo (Hymenaea courbaril). In OM, glycine (Neonotonia wightii), a leguminous perennial cover crop was sown. During the first years, it covered the soil almost completely, then (due to the less lighting level taht reach the soil). It diminished over time due to the less light reaching the soil, i.e., in the last years the cover was around 25%. Compost was applied on cacao trees annually. In OA the same legume cover were planted, compost was applied annually at a lower dose until 2016 too and the same shade tree species described in the CA systems were planted. The reduction of N. wightii was more pronounced in the organic agroforests in comparison with the organic monocultures. In the SA, no external inputs were applied. The arrengement of trees was the same as the other agroforestry systems. But additional species, both cultivated and spontaneous, were included. Among the cultivated species, the most important were corn and rice produced in2009, cassava in 2010 and pineapple and Bixa orellanain 2011. In 2016, coffee was introduced in all agroforestry systems (CA, OA and SA).

Data collection
In the middle of each plot, a phytosociological Braun-Blanquet (Ellenberg and Mueller-Dombois 1974) survey was carried out in an area of 256 m 2 (16 × 16 m). The abundance scores used were: 1, for solitary individuals; 2, for groups of individuals with a coverage below 1%; 3, for a coverage > 1%; 4, for a coverage > 5%; 5, for a coverage > 25%; and 6, for a coverage > 50%. A total of seven sampling events were performed, on June 2010, December 2010, January 2013, January 2015, January 2016, December 2016 and October 2017. The periods between sampling events slightly varied during the study due to the accessibility to the trial and the adjustment with weeding events, i.e., surveys were postponed as much as possible after weeding in such a way that the herbaceous species were more present and easier to identify.
All plant species present in the herbaceous stratum were inventoried (0-1 m height). However, for all analyses, only herbal, sub-shrubs and herbaceous liana species were taken into account, i.e., the grouping criterion was that the species should be able (based on direct observations or knowledge on the species) to reproduce at heights less than 1 m. The plant material collected for identification was deposited at the National Herbarium of Bolivia.

Richness
In order to compare the total number of species, we applied a generalized linear mixed effect model, with time and type of production as independent variables and the plot as a random factor. Additionally, we tested the effect of the production system at each sampling event. A Poisson distribution and an identity function were used for this purpose. The statistical significance was evaluated using a likelihood ratio test; the post hoc tests were performed using multiple comparisons for parametric models with the Tukey method.
For the purpose of understanding the differences in species dynamics, we calculated the number of total immigrations in each plot, i.e., the number of species recorded for the first time in a plot during any of the six successive samplings carried out after the first sampling. The number of total species lost (local extinctions) was estimated from the number of species that were once recorded but no longer found in either of the last two samplings. Finally, the number of "permanent" species refers to the species recorded in the first sampling and still present in the last two samplings. The number of immigrant species, extinctions and permanent species was compared between the different systems using generalized linear models.

Composition
To analyse the community composition, a non-metric ordination, multivariate analysis based on permutations ("adonis" in R) was used to determine the contribution of time, plot and type of production to the variation in composition (Anderson 2001). ANOSIM and Mantel tests were carried out; in both cases the Bray Curtis ecological dissimilarity index was used (Kindt and Coe 2005). For the ANOSIM test, the independent variable was the type of production (5 levels), and for the Mantel test, the independent variable was the spatial distances between plots (measured in meters). The adjustment values (R) of both types of tests were correlated with the time measured in months (Considering Jannuary 2010 as month number one).
To evaluate the differences in community composition, we constructed a matrix for the 20 plots where the following scores were assigned to each group of species: 3, to successful immigrants, which immigrated to a plot at some point and persisted until the last two samplings; 2, to permanent species, which remained after the first sampling; 1, to species that immigrated to the plot but subsequently disappeared; and − 1, to extinctions, i.e., species that disappeared and were not recorded in the last two samplings. Based on this matrix, a dissimilarity matrix was also constructed, and, subsequently, a dendrogram was created with the average linkage algorithm.
At the species level, we present the species in groups based on a Q-mode hierarchical cluster analysis, conducted for the 25 most abundant species. Additionally, we performed an indicator species analysis for groups of sites, i.e., in our case, the different production systems (De Cáceres et al. 2010). For this purpose, we used the Spearman coefficient as statistical indicator value, and a matrix where the abundance of species corresponds to the sum of the Braun-Blanquet scores of all the samplings within each plot.
All analyses were performed using R statistical software, with vegan packages for ordinations and Bray Curtis distance tests (Oksanen et al. 2018); lmerTest, lmtest and multcomp for generalized linear and linear analyses (Hothorn et al. 2008); and indicspecies to associate the species with a specific type or combination of production systems (Cáceres and Legendre 2009).

Results
Throughout the different surveys, 99 species of spontaneous herbs were registered, belonging to 31 families of angiosperms and 2 of ferns. The most diverse families were grasses (Poaceae) with 16 species, followed by Euphorbiaceae, with 8 herbaceous species (5 corresponding to the genus Euphorbia); Asteraceae and Cyperacea with 7 species each; and Acanthaceae, with 6 species.

Species richness
The average number of species per sampling event was 62.4. It varied very little from one sampling to another (SD = 3.5) and had no clear relationship with time (Pearson's R = -0.16). At plot level, the species richness was estimated to be 16.7 ± 4.8 (mean ± SD) species.
Overall, we found that the number of species was affected by the type of production (X 2 = 41.1, P < 0.001), time (X 2 = 23, P < 0.001), and their interaction (X 2 = 6.2, P = 0.013). We observed that it significantly decreased over time in both the organic and conventional agroforestry systems (CA: − 0.067, − 1 species every 14 months, P = 0.0001; OA: − 0.046, − 1 species every 22 months, P = 0.0172), while it was constant over time in the other systems (Fig. 1).
In terms of dynamics, throughout all sampling events, we registered an average of 19.6 ± 5.5 (± SD) immigrant species per plot. Also, 17.9 ± 4.8 species were lost and about 6.2 ± 2.8 remained since the first sampling event. The total number of both immigrant species and species lost differed between production systems (X 2 = 12.9, P = 0.0119; and X 2 = 9.8, P = 0.0447, respectively), but that was not the case for the number of permanent species (X 2 = 12.9, P = 0.1288). Numerically, OM had the lowest number of immigrant species and also of extinct species (Fig. 2). However, the pairwise comparisons showed only marginally statistical differences, i.e., OM had less losses than CM (7.3 losses more in CM than in Fig. 1 Average species number over time in the different production systems. SA: Successional agroforestry system (black squares); OA: Organic agroforestry system (white squares); OM: Organic monoculture (white circles); CA: Conventional agroforestry system (grey squares); CM: Conventional monoculture (grey circles). Letters indicate subgroups by pairwise comparison in each survey. Complete lines indicate that the systems had a significant relation with time; dotted lines indicate no significant relation with time Fig. 2 Mean and standard deviation of total immigrant species and total species lost in the different production systems. SA: Successional agroforestry system; OA: Organic agroforestry system; OM: Organic monoculture; CA: Conventional agroforestry system; CM: Conventional monoculture. Homogenous subsets given by pairwise comparisons of total immigrants are indicated in capital letters, while total losses appear in small letters OM, P = 0.099), and less immigrants than both SA (7.8 immigrants more in SA, P = 0.07468) and CM (10.3 immigrants more in CM, P = 0.00883). The combination of local immigration and extinction events resulted in a stable number of species along the sampling events for the CM, OM and SA systems, but in a decline for both the OA and CA systems ( Figs. 1 and 2).

Composition
In the first sampling events (June and December 2010), the similarity in composition between plots was higher when the plots were closer (R (MANTEL) = 0.27, P = 0.003; R (MANTEl) = 0.19, P = 0.022, respectively, Fig. 3). This relationship did not exist in subsequent samplings and the adjustment value decreased over time. The only system that maintained spatial autocorrelation over time was the SA system (minimum and maximum R (MANTEL) = 0.48-0.98).
In contrast, the impact of the type of production system on community composition was significant as from December 2010 (R (ANOSIM) = 0.25, P = 0.015), and its importance increased with time, reaching a value of R (ANOSIM) = 0.87, P = 0.001, in October 2017 (Fig. 2).
The non-metric scale ordination, stress = 0.27 (Fig. 4), revealed significant differences in species composition in relation to time (accounting for 9.8% of the variation), type of production system (25.7%), and plot (23.5%); in all cases P = 0.001. When the initial samplings were excluded, i.e., when only the period from 2012 onwards was considered, the plot (spatial position) explained 24.7% of the variation, the influence of the type of production rose to 40.2%, and that of the sampling event (time) dropped to only 2.6%; again, for all three factors, P = 0.001. Considering periods of approximately 2 years (i.e., December 2010 to December 2012, December 2012 to January 2015, and January 2015 to December 2016), we saw that the average value of dissimilarity in composition was higher in the first period compared to the following ones (the Bray Curtis average values for the same plot in each of these periods were 0.51, 0.43 and 0.38, respectively).
The cluster analysis of the communities dynamics revealed a differentiated immigrant/lost/permanent species composition and allowed to identify three main groups: one containing the CM plots, another including the SA plots, and a larger one with the rest of the systems (Fig. 5). This pattern indicates a specialization process in both SA and CM plots, with species only or mainly present in Fig. 3 Importance of spatial position and type of production for species composition over time. Adjusted R value of the Mantel test with spatial positions (circles), and of the ANOSIM test with types of production (squares). Black represents significant differences (P < 0.05) these systems (Table 1). In addition, the indicator species analysis showed a group of species with low densities that appeared exclusively in the CM system: Boerhavia diffusa, Portulaca oleracea, Paspalum paniculatum, Torenia crustacea, Eleusine indica, Pilea microphylla, Kyllinga odorata and Pityrogramma calomelanos. It is interesting to note that the last three were recorded for the first time in January 2015, December 2010 and October 2017, respectively, while the rest had been identified already in 2010. In the case of the SA systems, there were also some species with low densities that only appeared in this production system:

Species number decrease in agroforestry systems
We found that species richness tends to decrease over time in OA and CA systems, i.e., in those systems the total mean value of immigrant species was lower than the total number of species lost. This pattern may have two explanations. First, it may be the result of increased shading over time. With the growing of the cocoa and agroforestry trees, the canopy cover was denser and the light reaching the soil diminished. We clearly observed, for instance, a decrease in heliophilous species such as Momordica charanthia and Bidens pilosa, which were found in a lower number of plots over time (Table 1). Different studies have shown a negative relation between levels of shade and richness of herbaceous species (Soto-Pinto et al. 2002;Bobo et al. 2006;Ramos et al. 2015), mainly non-forest species (Cicuzza et al. 2011). In addition to the reduction of existing heliophilous species, the higher levels of shade can Another explanation for the observed pattern is that the herb assemblages in the CA and OA systems were poor in terms of quantity, i.e., there were few individuals of each spontaneous herbal species. This made these systems more vulnerable to local extinctions (Lundberg et al. 2000).
Contrary to the situation in the CA and OA systems, no decline of species was observed in the SA system, even though the shade increased similarly over time. This could be related to the different composition of the SA system, with species that were found exclusively in it (see below). The presence of these species could be hampered by herbicide applications in the CA systems, and by competition with the leguminous cover crops in the OA system.

Resistance to invasion and cover crop
The OM system had fewer immigrants, i.e., less number of newly established species, than the CM system. This is explained by the presence of the perennial cover crop, which densely covered the soil. Our results confirm previous observations of temperate herb assemblages, which showed that perennial dominant species can generate a community with low susceptibility to invasions (Emery and Gross 2007). This proves the potential of this technique to prevent the establishment of worldwide distributed species, including exotic and invasive ones. In this line, we observed that the OA system was, numerically but not significant, more dynamic (more immigrants and losses) than the OM, which could be related to the faster reduction of the cover crop coverage in the former system compared to the latter.

Spatial position and change rate
In composition analisys, spatial autocorrelation was only detected at the beginning of the study, and was not relevant at all after 3 or 4 years. This finding coincides with those of a previous study conducted three years after the establishment of the trial, where the importance of the position of the plots for the species composition was very low (Marconi and Armengot 2020). Spatial autocorrelation indicates that the herb assemblages may strongly be influenced by the surrounding conditions (Legendre 1993) or by repopulation processes (Caners et al. 2009). Other studies performed in the same trial indicate that weeding was one of the most laborious activities, mainly during the first years (Armengot et al. 2016). For instance, after seven years of planting, the time spent in weeding tasks was only 45% of the time invested in the same tasks during the first year of the trial (data not published). Later on, a good management of the Fig. 5 Dendrogram representing the differences in the species composition dynamics. SA: Successional agroforestry system; OA: Organic agroforestry system; OM: Organic monoculture; CA: Conventional agroforestry system; CM: Conventional monoculture. Point size represents time. The length of the bar represents, on the left side, the total number of species lost, and, on the right side, the total number of immigrants. The colour represents the tendency in the plot based on Pearson correlation with time, where white stands for assemblages in decline (R < − 0.5), grey for stable assemblages (− 0.5 < R < 0.5), and black represents a growing pattern (R < 0.5) in species richness pre-existing or surrounding areas could have helped reduce the dispersion of the propagules.
Evidently, there is a filtering or selection process derived from the conditions of the different production systems, which began very early, already at the time of the second survey (one year after the cocoa planting). The percentage of compositional variation explained by time, the values of the Bray Curtis tests in different periods, and the shape of the curve of the adjusted R values in the ANOSIM test showed that the magnitude of this process decreased over time. It was expected that, at the beginning of the implementation of the new management conditions (e.g., agrochemicals or cover crops), the vegetation change rate would be higher than later on, when the new conditions were fully established. Specialization in the species composition of SA and CM Contrary to our expectations, both the SA system, the most complex and less intensively managed system in the trial, and the CM, the most intensively managed system, harboured more immigrant species than the other systems. However, the identity and composition of these new species was completely different in the two systems.
The new species found in the CM were identified as either resistant or tolerant to glyphosate. It is known that some glyphosate-resistant species can travel long distances, as is the case of Conyza canadensis (Dauer et al. 2007). Therefore, CM may interact with longer distance seed sources. Invasive species tend to be more fit (for instance, they have a higher number of seeds) and more resistant to water stress than non-invasive species (Van Kleunen et al. 2010), which enables them to travel along roads (water-stressed spaces) or open spaces. This phenomenon of immigration from distant sources seems to be one of the reasons why conventional systems promote biotic homogenization (Marconi and Armengot 2020).
The SA system contained a group of species that were not present in any other system. Furthermore, those species were the only ones in which spatial autocorrelation was maintained over time. We had previously found that SA systems tend to harbour species in geographically restricted areas (Marconi and Armengot 2020), and we could prove in this study that this remained true over time. These species also thrive in secondary forests (personal observations), which implies a close relationship with the pre-existing vegetation, basically consisting of secondary forests, therefore demonstrating that these systems promote species conservation. In this sense, it is important to highlight the importance of SA systems for biodiversity conservation. Other taxonomic groups such as birds (Naoki et al. 2017) and Thysanoptera (Zegada et al. 2020), a potential cocoa pollinator, were also more abundant in this production system. In addition, higher cellulase activity was observed in the SA systems (Alfaro-Flores et al. 2015). These studies show that SA agroforestry could contain complex ecological interactions (i.e., plant -fungus, plant-animal, plant -plant). Future research should address these interactions and their importance for biodiversity conservation in cocoa production.
Generalist species on OM, OA and SA systems The organic (OA and OM) systems and the SA system were close to each other in the NMDS ordination (Fig. 3). These types of production systems have in common the non-use of agrochemicals and soil protection. The compositional similarity of these systems is due, in part, to the presence of a group of species that grew in extension within the plots (see Table 1, first block). Trying to understand which morphological or physiological characteristics give them this ability to expand, we found that, except for Euphorbia poepiggi, they are all low-bearing species that reproduce by stolons. These possibly allow them to better survive mechanical weeding events. It is also possible that these species do not develop well in conventional plots because of the application of glyphosate. On the other hand, the fact that these species expanded into organic plots was due to the gradual disappearance of the coverage created by the legume Neonotonia wigthii.

Conclusions
Agronomic management affects the composition and richness of the herb communities in cocoa production systems, and leads to a selection process that begins very early and slows down over time. Each production system presents very different dynamics. Organic monocultures show resistance to species immigration and species losses; both types of agroforestry systems (organic and conventional) reduce the number of species over time due to shade increase; and successional agroforestry systems and conventional monocultures present similar levels of species immigration, but with very different and specialized compositions: in the first case, shade-tolerant species probably coming from nearby sources such as secondary forests, and, in the second, worldwide distributed heliophytes presumably arriving from roads and other perturbed habitats. Organic management and agroforestry especially complex successional ones, should be therefore promoted in cacao production since they both favour herb assemblages with a high conservation value and prevent the establishment of harmful and globally distributed species.