Win-win opportunities combining high yields with high multi-taxa biodiversity in tropical agroforestry

Ethical statement

Ethics approval was obtained for this study from the ethics committee of the University of Goettingen (Chair: Prof. Dr. Peter-Tobias Stoll) under the reference number 17./04.22Wurz.

Study area

All plots were situated in northeastern Madagascar in the SAVA region (Supplementary Fig. 1). The natural vegetation is tropical lowland rainforest, but deforestation rates are high^30,67.

The region is globally and nationally one of the most biodiverse places with high levels of endemism^17,68. Forest loss is mainly driven by slash-and-burn shifting hill rice cultivation⁵⁸. The region is characterized by a warm and humid climate with an annual rainfall of 2255 mm and a mean annual temperature of 23,9 °C (mean value of 60 plots extracted from CHELSA climatology⁶⁹). Vanilla is the main cash crop in the SAVA region, making Madagascar the main vanilla producer globally^21,22. Vanilla prices have shown strong fluctuations over the past years, with a price boom between 2014 and 2019 triggering an expansion of vanilla agroforestry in the region^22,23.

Study design

We selected 10 villages based on the 60 villages selected within the Diversity Turn in Land Use Science project²² (Supplementary Fig. 1). We selected the villages based on the list of villages for our study region from official election lists which listed all villages within a fokontany individually²². Village boundaries, demographics, infrastructure were defined based on a rapid survey with the village chief. Among the 60 villages, we considered all villages without coconut plantations, with less than 40% water (river, sea, and lakes) to avoid a strong influence of water elements and with forest fragments and shifting cultivation present within a 2 km radius around the village. Two of these 17 villages overlapped within a 2 km radius of the villages, thus we randomly selected one of them, resulting in 14 villages. We visited these 14 villages in a randomized order and stopped after we found 10 villages which fulfilled the necessary criteria (all land-use types present, willing to participate). In each of the 10 villages, we selected three vanilla agroforests, one forest fragment, and two fallows. Overall, we studied 60 plots across 10 villages and 10 plots in one protected old-growth forest (Marojejy National Park). All plots had a minimum distance of 260 m and a mean minimum distance of 794 m (SD = 468 m) to each other. Plot elevation ranged between 10 and 819 m.a.s.l. (mean  = 205 m, SD = 213 m; Supplementary Table 20).

Plot selection

In each of the 10 villages, we selected three vanilla agroforests with low, medium, and high canopy closure, respectively, covering a within village canopy cover gradient. To refine our vanilla agroforest classification, we used interviews with the plot owners to categorize all vanilla agroforests based on land-use history into fallow- and forest-derived agroforests¹⁵. Forest-derived vanilla agroforests are established within forest fragments, which have been manually thinned of dense understory vegetation. Fallow-derived vanilla agroforests are established on formerly slashed and burned plots, where vegetation has been cleared for hill rice production (shifting cultivation system locally called tavy). Out of our 30 vanilla agroforests, 20 vanilla agroforests were fallow-derived and 10 vanilla agroforests were forest-derived, roughly matching the proportion of fallow- and forest-derived vanilla agroforests across the study region (70% are fallow-derived vanilla agroforests, 27% are forest-derived vanilla agroforests and 3% of unknown origin²².

In addition to vanilla agroforests, we selected one forest fragment in each village. Forest fragments were located inside the agricultural landscape and were remnants of the once continuous forest; these fragments are frequently used for natural product extraction. Forest fragments have not been burned or clear cut in living memory, yet the ongoing resource extraction results in a much simplified stand structure and fewer large trees compared to old-growth forest¹². Furthermore, we chose one herbaceous and one woody fallow in each of the 10 study villages. Both fallow types form part of the shifting hill rice production cycle and represent the fallow period at different stages after the crop production. Herbaceous fallows have been slashed and burned multiple times with the last cultivation cycle at the end of 2016, one year prior to the first species data collection in 2017, and thereafter left fallow¹¹. The continuous succession of herbaceous fallows turns them into woody fallows with the domination of woody plants including shrubs, trees, and sometimes bamboo. Our 10 woody fallows have last burned 4–16 years before data collection. In this study, we combine both herbaceous and woody fallows into the category “fallow”. Generally, fallows occur in different forms in the study region. The characteristics of fallows depend on the frequency of past fires and the length of fallow periods in between crop cultivation¹¹. Frequent burning results in a loss of native and woody species and a dominance of exotic species and grasses¹¹. In later fallow cycles, fern species increasingly appear¹¹.

Due to the commonly repeated slashing and burning, secondary forests are very rare in the study region. Shifting cultivation prevails in Madagascar⁷⁰, because it is an important option for people to grow food because means for agricultural intensification are scarce. According to our baseline survey (performed in 60 villages in our study region), 90% of the interviewed farmers grow rice for subsistence in addition to growing vanilla²². Out of this sample, 64% of farmers grow rice in irrigated paddies and 26% of farmers use shifting cultivation.

We also studied 10 plots at two sites in Marojejy National Park, the only remaining, continuous old-growth forest at a low altitude in our study area⁷¹. We chose accessible old-growth forest plots with a minimum distance of 250 m from the forest edge. Five of the 10 old-growth forest plots were located in Manantenina Valley, the other five old-growth forest plots were situated in the eastern part of Marojejy National Park, called Bangoabe area. Illegal selective logging has occurred in some parts of the park. During our plot selection, we avoided sites with traces of selective logging.

Land-use history classification

To collect information on the land-use history or farm history, interviews with farmers are common^72,73. We did interviews with the plot owner. Questions on land-use history were binary (forest-derived or fallow-derived) and did not include information on the detailed land-use history (e.g. frequency of burning, past crop systems). Thus, we consider this selfreported data very reliable. The land-use categorization derived by farmers was confirmed by our visual plot inspections (forest-derived vanilla agroforests do have a quite distinctive vegetation structure compared to fallow-derived vanilla agroforests). Additionally, data on tree species composition and soil characteristics show evident differences between the categories and back up the binary land-use history categorization. Analysis of tree species composition showed that fallow- and forest-derived vanilla agroforests differ significantly in tree species composition¹². Soil analysis (see Fig. S9) showed that our fallow-derived vanilla agroforests are associated with fertility-related variables such as an increase in calcium, pH, nitrogen, and phosphorus, which is common after slas-and-burn agriculture^74,75.

Plot design

We collected species data on plots with a radius of 25 m (1964 m², 0.1964 ha). We established our circular plots in a homogeneous area of the land-use type or forest. Adjacent land uses were usually different because farmers generally own small-scale land with a mean size of 0.66 ha (mean size of agroforests). We assessed vanilla plant data (yield, vine length, vine age, planting density) on 36 vanilla pieds on each of 30 circular vanilla plots (Supplementary Fig. 8). We defined one vanilla pied (foot in French) as the combination of a vanilla vine and a minimum of one support tree. The 36 vanilla pieds were evenly selected in each of the circular plots based on a sampling protocol to ensure comprehensive and unbiased sampling. We chose vanilla pieds independent of age, length or health condition. We marked the 36 selected vanilla pieds per plot with a unique barcode to assess vanilla yield (April 2018) and other plant health variables on the same plant (not used in this study). However, for 37 vanilla pieds (out of a total of 1080 marked vanilla pieds), the barcodes were lost or unreadable and we selected a new plant closest to the original position (independent of age, length, or condition) and marked it with a new unique barcode. We measured the size of the vanilla agroforest by walking with the agroforest owner and a hand-held GPS device at the perimeter of the plot.

Vanilla planting density

We counted each vanilla pied on each 25 m circular plot by dividing the plot in four-quarter segments. We calculated the area of each 25 m radius plot including slope correction and calculated vanilla planting density (vanilla pieds per hectare) by dividing the number of vanilla pieds by the slope-corrected plot area.

Vanilla yield

We measured yield on 30 vanilla plantations (10 forest-derived vanilla plantations and 20 fallow-derived vanilla plantations); three in each of our 10 study villages. We measured vanilla yield on a total of 36 vanilla pieds between March and April 2018. We assessed the vanilla yield before harvest to ensure an accurate yield assessment due to two reasons. Firstly, vanilla pods are commonly harvested successively due to their differing pollination date and maturity requiring multiple visits over several weeks. Secondly, theft of vanilla pods is commonplace around harvest time. We, therefore, estimated the weight of the on-plant-hanging vanilla pods by measuring pod volume and relating this to a prior established volume–weight correlation. This is possible because vanilla pods only grow in length and width in the first 8 weeks of their development⁷⁶. Our yield assessment consisted of one interview part with the plot owner and one measurement part. The interview part included questions about the occurrence of theft and early harvest on the plantation. During the measurement part, we assessed the number, diameter, and length of all vanilla pods. We measured vanilla pod length with a ruler starting at the junction of stem and pod until the tip of the pod without considering the bending of the pod. We measured the diameter at the widest part of the pod using a caliper. We firstly calculated pod volume based on the standard volume cylinder formula using the measured diameter (cm) and length (cm): V = πr²h.

Secondly, we calculated the weight (g) of each pod by using the linear regression equation (y = bx + a) of a weight–volume correlation of 114 vanilla pods from 114 different agroforests (weight, length, and diameter of these 114 green vanilla was assessed post-harvest in 2017). We calculated the weight of all measured pods of the harvest in 2018 based on the formula:

Here, we divided the pod diameter (mm) by 20 to obtain the radius and to transform millimeters to centimeters. Weight was defined as volume*0.5662 + 0.9699. No vanilla pods were stolen or already harvested on our 36 vanilla pieds and hence we did not need to account for it in our vanilla yield calculation.

Vanilla vine length

We assessed vanilla vine length for all 36 vanilla pieds (same vanilla pieds as used for the yield assessment) on each plot by measuring the total length of the vine from the lowest to the highest part with a measuring stick. If the vanilla vine was looped on the support tree (= vanilla vine is hanging in multiple loops on the support tree), we measured from the top height of the looping of the vanilla vine until the lowest height of the vine. At the medium height of the vanilla vine, we counted the number of times the vanilla vine passed through. We calculated the total length of the liana by multiplying the maximum height of the vanilla vine by the number of times the vine passed through the middle. In some cases, the vanilla vine looped at two different heights, we thus considered the middle between the two looping heights as the top height. If vanilla vines grew on two different support trees, we considered them as one vanilla pieds if support trees were <30 cm apart. If the distance between both support trees exceeded 30 cm, we considered only the support tree with the most vines for the measurement.

Pollination labor input

We performed a longitudinal survey with the plot owners of our 30 vanilla agroforests from October 2017 to October 2018. The questionnaire was pre-tested with 30 farmers in September 2017. All participants were trained by the four research assistants on how to use pictogram-supported questionnaires. Subsequently, a feedback workshop was held to adapt the pictograms and optimize the entire questionnaire. The pictograms had to be filled every day as a diary. Besides pollination labor input, we assessed the time spent on plantation establishment, planting, weeding, pruning, plantation safeguarding, harvesting, preparing (fermenting, drying, sorting), and selling of vanilla (not considered in this analysis). Every fortnight, trained assistants visited farmers to collect the diary questionnaires. Data entries that appeared unusual were verified with the farmers by the assistants. The diary questionnaires included questions on family labour input for pollination as well as other agricultural activities, such as weeding, harvesting, curing of vanilla, and others.

We decided to use pollination labour input as the only variable of labour input in our analysis for the following two reasons. Firstly, pollination is the most important and labour intensive part of the production⁷⁷. Secondly, pollination has a defined time frame because vanilla only flowers between October and December²². Thus, the hours of pollination labour input are easy to disentangle and define for the farmers. In contrast, other tasks such as weeding, pruning, and planting often happen continuously and in parallel and are thus harder for the farmer to depict in working hours. Pollination is a labor-intensive activity as every vanilla flower is pollinated manually and requires agricultural know-how. Thus, the working hours can be related to the number of flowers or/and the worker’s know-how. We calculated pollination labour input per hectare by summing all working hours (Oct 2017–Oct 2018) and dividing it by the slope-corrected agroforest size. For three out of 30 vanilla agroforests pollination labour input was missing; we thus used the mean value of all 30 agroforests for the three missing values.

The household head who filled the pictograms received 10,000 Ariary (roughly 2.50€) per month. Pictograms are drawings made by a local artist which visually describe each of the working steps of vanilla cultivation (e.g. planting vanilla vine, weeding plantation, pollination). The amount of compensation was recommended as a reasonable compensation by locally experienced Malagasy project members. The sum was handed out by the local research assistants at the end of each month. All participants of the surveys were informed that participation is voluntary, that they can leave the survey anytime, and that all data is anonymized, i.e., no personal data will be published or shared with third parties. Guidelines of “Good Scientific Practice” by the University of Göttingen were adopted (adapted based on the recommendations of the codex for good scientific practice from the DFG, German Research Foundation; https://www.unigoettingen.de/en/good+scientific+practice/567647.html).

The interview methodology (i.e. informed consent by the test persons as well as the questionnaires) was evaluated by the ethics committee of the University of Goettingen (https://www.unigoettingen.de/en/534983.html) and complied with their principles of the Higher Education Act of Lower Saxony (NHG) and the constitutionally protected right of academic freedom (Reference number: 17./04.22-Wurz). Co-author H.H. designed the survey and trained the local assistants with B.F. until the assistants were able to enter data in a consistent and standardized manner. F.A., F. S. B., and the trained assistants conducted the interviews. The research assistants collected the data bi-weekly but visited the households weekly. F.A. and F. S. B. checked entry data bi-weekly and clarified false entries. Our research assistants were recruited mainly from the student population of the regional CURSA university center in Antalaha, which is located in the research area. The students speak the local Malagasy dialects. Questionnaire in the original language, used pictograms, and the reporting sheet for farmers are available on Open Science Framework: https://osf.io/z5uxs/?view_only=1bd699c5cda64023963e058254a33eec.

Vanilla plant age

We assessed vanilla plant age by asking the farmer for each of the 36 vanilla pieds per plot. The Malagasy field researchers Evrard Benasoavina, Thorien Rabemanantsoa, and Gatien Rasolofonirina walked with the farmer to each vanilla plant (see the question “Ask farmer: How many years ago was the liana at this pied planted?” in uploaded original interview (yield assessment) on Open Science Framework: https://osf.io/wa2xn/?view_only=1bd699c5cda64023963e058254a33eec). Here, the age referred to the vanilla vine but not the support tree. In preparation for the farmer interviews, we prepared handouts in Malagasy language to inform the farmers about the scientific goals and the content of our data collection (see handout on Open Science Framework: https://osf.io/ndrxg/?view_only=1bd699c5cda64023963e058254a33eec).

Canopy closure

We measured mean canopy closure at five subplots of our circular plots by taking hemispherical images with a Nikon D5100 camera, equipped with a Sigma Circular Fisheye (180°) 4.5 mm 1:2.8 lens. The camera was fixed on a tripod at 2.4 m height above vanilla support trees and understory vegetation. We selected the images with the best contrast of sky and vegetation using the histogram-exposure protocol and calculated canopy closure using a minimum thresholding algorithm^78,79.

Slope and elevation

We used the 30 m-resolution digital surface model “ALOS World 3D” by Japan Aerospace Exploration Agency (JAXA) to assess the mean slope and the mean elevation of each plot.

For all values, we applied slope correction⁸⁰.

Landscape forest cover

We calculated forest cover in a 250 m radius around each plot center based on binary forest cover data from 2017 with a 30 m resolution³⁰ and the R-package raster⁸¹. To reliably identify forest, tree cover maps and satellite imagery with a tree cover threshold of 75% were combined³⁰. Here, agricultural lands such as tree plantations are excluded by combining historical forest maps with up-to-date forest cover change maps⁸². We chose a 250 m radius as a compromise between mobile and immobile taxa.

Understory vegetation cover

We estimated the vegetation cover (percentage woody and herbaceous cover) visually for the 0–2 m layers in % of five subplots on each plot (located in the plot center and at 16.6 m from the center in each cardinal direction) and calculated the mean understory vegetation cover per plot. We did not consider vanilla pieds in the estimation of the understory vegetation cover.

Soil characteristics (PC1)

We took soil samples with a MacFadyen soil corer (5 cm diameter, 295 ml, 0–15 cm depth). We divided the plots into eight subplots, four subplots at 8.3 m distance to the plot center (inner area) and four subplots at 16.6 m distance to the plot center (outer area). In total, we collected four cores in the inner and outer area each, resulting in two mixed soil samples per plot. We stored each soil sample in a zip-lock bag until laboratory analysis. In the laboratory, we measured pH (H₂O) with the fresh soil samples using 1:10 humus/water suspension after 24 h of equilibration. We measured pH (KCl) by adding 1.86 g KCl. Mean pH values by plot were calculated in logarithmic and back-transformed by using exponential. The remaining soil was dried at 70 °C and ground. We measured organic carbon (C_org) (mmol/g dry soil), total carbon (C) (mmol/g dry soil) and nitrogen (N) (mg/g dry soil) concentrations, and organic C-to-N ratio (mol/mol) by using the C/N elemental analyzer (Vario EL III, elementar, Hanau, Germany). Additionally, we determined effective cation exchange capacity in µmol/g dry soil of potassium (K), magnesium (Mg), calcium (Ca), iron (Fe), manganese (Mn), hydrogen (H), and aluminum (Al) by digesting oven-drier soil material in 65% HNO₃ at 195 °C for 8 h. Samples were analyzed with inductively coupled plasma optical emission spectrometry (ICP-OES) (Optima 3000 XL, Perkin Elmer, USA). We calculated the total effective cation exchange capacity (µmol/gTB) by the sum of H, P, Mg, Ca, Fe, Mn, and Al and total base saturation (%) as the sum of K, Mg, and Ca. Extractable phosphorus (P) Resin (µmol/gTB) was measured using resin bags, which were placed in a soil–water suspension. Then, P was re-exchanged with NaCl and NaOH solutions and quantified colorimetrically after blue-dyeing.

Due to the high number of soil characteristics measured and possible multicollinearity, we calculated Spearman correlations using the cor and corrplot function from the corrplot R-package⁸³. Based on collinearity, we excluded total effective cation exchange capacity, total base saturation, total C, organic C, effective cation exchange capacity of Mg, effective cation exchange capacity of pH percolate, and the effective cation exchange capacity of H percolate (Correlation matrix, Supplementary Fig. 9). We used the remaining variables, i.e. effective cation exchange capacity of Ca, K, Al and Mn, pH(KCl), total N, resin P, and organic C-to-N ratio to perform a principal component analysis (PCA) using the R-packages ggbiplot⁸⁴ and factoextra⁸⁵. The soil PC axis 1 (explained 45%) was mostly related to effective cation exchange capacity (µmol/gTB) of Ca, and K as well as pH(KCl), nitrogen (mmol/gTB), P(resin) (µmol/gTB) and the organic carbon–nitrogen ratio (mol/mol) while the soil PC axis 2 (explained 21%) was related to exchange capacity of Mg and Al and the Corg-to-N ratio (Supplementary Fig. 10). The coordinates of PC axis 1 were used as a proxy of soil characteristics for further analysis.

Trees

We sampled trees on all land-use types except herbaceous fallows between September 2018 and January 2019¹². Access was denied to two fallow-derived vanilla plantations, resulting in 58 plots assessed overall (including 28 vanilla agroforests). We did a full inventory of all trees with freestanding stems of ≥8 cm diameter at breast height in each plot. This included trees, arborescent palms, herbs, and tree ferns but excluded lianas. We identified tree species with the help of a local tree expert (Chrysostome Bevao) and a taxonomic expert (Patrice Antilahimena) from Missouri Botanical Garden (Antananarivo, Madagascar). We derived information on origin and endemism for each species from the Tropicos Madagascar Catalog⁸⁶. Voucher specimens are kept at the National Herbarium Tsimbazaza, Antananarivo (TAN) and the herbarium of the University of Mahajanga. Out of the 454 assessed species in this inventory, 276 (51%) were endemic to Madagascar.

Herbaceous plants

We sampled herbaceous plants in eight subplots of 4 m² each (32 m² overall) between September 2018 and December 2019. In each subplot, we assessed vascular plant species without apparent wood at maturity⁴⁰. We accounted for the possible seasonality variation of the plant phenology by sampling one village after another. In each village, we collected data on each land-use type except for the old-growth forests. Hence, the observations for each land-use type cover the study period along with the possible phenology variations. We determined each species’ endemism status from the Tropicos Madagascar Catalogue⁸⁶. We stored all herbarium specimens at the Plant Biology and Ecology Department at the University of Antananarivo in Madagascar. From the 299 species assessed in this study, 59 species (20%) were endemic to Madagascar.

Birds

We sampled birds during two 40 min point counts per plot with two observers per point count⁸⁷ following a commonly used standardized method⁸⁸. In all villages, we conducted one point count between September and December 2017 with co-author D.M. as the main observer and co-author R.A. as a second observer. The second point count was done between August and December 2018 with Eric Rakotomalala as the main observer and the co-authors D.M. or S.D. as the second observer. We exchanged the order of plot visits in the second year to minimize seasonal bias. For old-growth forest, we did all point counts in 2018; one in August 2018 with E.R. as the main observer and D.M. as the second observer, and a second count in December 2018 with E.R. as the main observer and S.D. as a second observer. E.R. and D.M. are both experienced birders and familiar with the encountered bird species due to previous fieldwork in Madagascar. S.D. and R.A. were trained prior to fieldwork to recognize calls and morphology of Malagasy bird species. The second observer was responsible for entering of data and supporting the identification. Main observers were responsible for spotting, hearing, and identifying the birds. On 63 plots we started one-point count around sunrise and one at least one hour after sunrise; on 7 plots this alteration was not possible due to logistical constraints. After arriving at the plot center, we waited for a minimum of three minutes to allow the birds to settle. We noted the conditions including rain (no rain, drizzle, light rain, heavy rain) and wind (Beaufort 1–12⁸⁹) before each point count. We only started point counts under good weather conditions, meaning no rain and wind equal to or less than Beaufort 4. If weather conditions deteriorated for more than 10 min from the start of the point count, we aborted the point count and started again later or the next day under better conditions. For calculating bird species richness per plot, we disregarded observations only in flight and outside the 25 m radius of plots. We defined species as endemic if only occurring in Madagascar according to BirdLife species fact sheets⁹⁰. Overall, we assessed 51 bird species of which 31 species (61%) were endemic to Madagascar.

Amphibians and reptiles

We sampled amphibians and reptiles using repeated time-standardized search walks for 45 min by two observers⁹¹. We visited each plot both during the day and at night both during the driest (one nocturnal and one diurnal search between October and December 2017; one nocturnal and one diurnal search between August and December 2018) and the wettest period (one nocturnal and one diurnal search between January and April 2018 or in February 2019). We did so during the driest period and the wettest period. We systematically walked the circular plot in a zigzag pattern always starting from the West part toward North, East, South, and end in West to avoid counting twice the same individual during observation in one of the plots. We actively checked microhabitats to detect individuals hiding therein (e.g., individuals hiding under rocks, in leaf axils, tree barks, tree holes, leaf litter, or deadwood). When encountering an individual, we stopped the standardized search time and identified the individual⁹². We identified individuals based on morphological characteristics with the help of field guides^93,94,95. We took DNA samples to determine the species for those individuals that proved difficult to identify using morphological characteristics only. To retrieve a DNA sample, we collected muscle or toe clips as tissue samples, conserved in 90% of alcohol. We stored DNA samples at the Evolutionary Biology laboratory at TU Braunschweig. We also took photos of specimens that we did not identify to species level (ventral, back, and flank view). Until release, we kept them in a ventilated bag to retain moisture. We released all specimens after completing the full-time-standardized search. We categorized endemic reptiles and endemic amphibians as species/morphospecies only occurring in the country of Madagascar^96,97. We found in total 58 amphibian species of which 57 species (98%) were endemic. In terms of reptiles, we found 61 species of which 74% (45 species) were endemic to Madagascar (and 7 species without defined endemism level).

Butterflies

We sampled butterflies with fruit traps and time-standardized netting between August and December 2018⁵³. We baited fruit traps with fermented bananas and deployed the cylindrical nets for 24 h. Before deployment, we fermented bananas for 48 h in an air-tight container.

On each plot, we installed a total of 8 fruit traps. We deployed four fruit traps at 16.6 m distance from the plot center in the four main cardinal directions and the other four fruit traps at 20 m distance from the center in the four intercardinal directions. We used fish lines covered with vaseline creme to hang the fruit traps. The Vaseline prevented ants to intrude on the fruit traps. In addition, we avoided any contact with branches on the fruit traps. We caught butterflies with a fruit trap with a 20 cm Cone Opening (90 cm long hanging 1.5 m above the ground. On plots without trees, we installed fruit traps on a support stick (in rice paddy and herbaceous fallow). During the 30 min time-standardized netting, we caught butterflies within an imaginary 2 m wide box to each side of the net while walking at a slow and steady speed in a zig-zag to equally cover the plot area. The net had a circular frame with a nylon mesh on a 1.5 m telescopic handle. We performed the timestandardized netting in dry and low-wind conditions only, either in the morning (8 a.m. to 12 p.m.) or in the afternoon (1 p.m. to 5 p.m.). We then collected and dried all captured butterflies and identified them to species level in the laboratory (moths excluded). Identification was done by Annemarie Wurz, David Lees and Sáfián Szabolcs. We categorized endemic butterflies as species/morphospecies only occurring in the country of Madagascar and updated with expert consultation by David Lees⁹⁸. All identified specimens remain at the insect collection in the Department of Crop Sciences, section Agroecology, University of Göttingen, Germany. Overall, we found 84 butterfly species of which 49 species (58%) were endemic to Madagascar.

Ants

We sampled ground-foraging ants using bait and pitfall traps⁵⁴. We conducted the sampling in all villages between October and December 2017, and in the old-growth forest in August and December 2018. We established five sampling stations per plot: one at the plot center, and four at 16 m distance from the plot center; one in each cardinal direction. We then set bait and pitfall traps 10 m apart at each sampling station. We baited the bait traps using sardine and sugar on two white flat plastic plates with a diameter of 13 cm and placed the two plates about 5 cm apart. We left the baited traps for 30 min before collecting ants for 30 s. We buried the pitfall traps (plastic cups of 9 cm top diameter, 11 cm depth, and 6 cm bottom diameter) in the soil and filled them one-third with 70% alcohol and a few drops of soapy water. We emptied the pitfall traps after 48 h and identified ants to species/morphospecies level in the laboratory. Identification was done by Anjaharinony Rakotomalala, using available identification keys^99,100,101. Cross-checking of the identification of the species was done with expert consultation by Jean Claude Rakotonirina (species of Leptogenys), Nicole Rasoamanana (species of Camponotus), and Manoa Ramamonjisoa (species of Tetramorium). We defined endemic ant species as those species only present in the country of Madagascar¹⁰². We stored voucher specimens at Madagascar Biodiversity Center, Antananarivo, Madagascar (MBC). We recorded in total 123 ant species of which 55 species (45%) were endemic to Madagascar.

Statistical analysis

We performed all statistical analyses in R version 3.6.3¹⁰³. To assess the overall species richness across all taxa, we calculated their mean normalized (endemic) richness. To ensure the equal weight of all taxa despite differences in the range of species richness per taxa, we normalized species richness measures with Min−Max Scaling, which is a common method in multi-taxa studies to compare species richness across different taxa and to calculate their overall richness^104,105,106. Using this method, the data are linearly transformed to y = (x−min(x))/(max(x)−min(x)), where x is the set of observed values of species richness. As a result, the normalized species richness (y) then ranges between 0 and 1, with 0 referring to the lowest observed value and 1 to the highest observed value, respectively. We then calculated the average of the normalized values to obtain the mean normalized species richness, i.e., the overall richness of all taxa at the plot level. Our results were robust (no relationship between biodiversity and yield) independently of the metric used, i.e., normalized richness or multidiversity.

We investigated the relationship of species richness or normalized (endemic) richness with vanilla yield in 30 vanilla agroforests using glmmTMB models¹⁰⁷ or a linear mixed-effects model¹⁰⁸, respectively. We used glmmTMB due to its bigger flexibility (especially in the case of zero inflation) and its higher speed when using multiple fixed effects as well as random effects¹⁰⁷. We treated vanilla yield (sqrttransformed) in interaction with land-use history (fallow vs. forest-derived) as an explanatory variable and site (village or old-growth forest site) as a random effect. We scaled and sqrt-transformed vanilla yield due to a few high-yielding plots inflating the data distribution.

We assessed the environmental and management-related covariates of species richness with a glmmTMB model¹⁰⁷ and vanilla yield using a linear mixed effect model¹⁰⁸ with yield sqrt-transformed. We tested for correlation among the covariates using the Spearman coefficient (Supplementary Fig. 11). For both models (species richness and vanilla yield) we used canopy closure, soil characteristics, slope, landscape forest cover, understory vegetation cover, elevation, planting density, pollination labour input, vanilla vine length, vanilla plant age as explanatory variables. We added site as a random effect and scaled all explanatory variables with the scale function. The function divides the centered columns (value per land-use type−mean value across land-use types) by their standard deviation across all land-use types. For the variable pollination labour input, data from three plots was missing. To avoid omitting these plots from the analysis, we performed a linear regression of available household labour with pollination labour/ha and found a good correlation between the two (p-value: 0.005, estimate = 52.82, SE = 17.38, Supplementary Fig. 12). Based on the number of household members (for which we have available data for the three missing plots), we then predicted expected pollination labour input for the three missing plots, using the estimated relationship of the two from the mentioned regression.

For all models assessing the environmental and management-related covariates of species richness, we used likelihood ratio tests using maximum-likelihood estimation¹⁰⁹ to assess the statistical significance of individual variables. Specifically, we compared the full model versus a reduced model without the individual variable (single term deletion). The final model included all explanatory variables with residual deviance χ² < 0.05 (“chi-squared value”) in the model comparison with the full model.

We assessed differences in species richness or normalized (endemic) richness between forest, forest fragment, forest-derived vanilla agroforests, fallow-derived vanilla agroforest, and fallows by fitting a glmmTMB model¹⁰⁷ or a linear mixed-effects model¹⁰⁸, respectively. We used a glmmTMB model for all analyses with count/discrete data as a response. For mean normalized richness, we used a linear mixed-effects model due to the continuous data. Here, we treated land-use type as an explanatory variable and site as a random effect. We tested if assumptions for normality and homogeneity of variances were met. If the residuals were homoscedastic, we used the glht function from the multcomp package¹¹⁰ applying Tukey’s all-pair comparisons with Bonferroni correction to assess differences in species richness between land-use types. If the residuals were heteroscedastic (withingroup deviation from uniformity significant with DHARMa quantile test¹¹¹), we used the nonparametric Kruskal–Wallis test, followed by pairwise Wilcoxon test. For all our models we used simulation plots implemented in the DHARMa package to validate our model fit¹¹¹. If our glmmTMB models were under- or over-dispersed, we changed the model family to compois (Conway–Maxwell Poisson distribution) and negative binominal (nbinom2), respectively. We assumed normal distribution for the models with mean normalized (endemic) richness as a response. To assess marginal and conditional R², we used the delta method using the rsquaredGLMM function of the package MuMIn¹¹².

We plotted all model results using RBase (Figs. 1, 2)¹⁰³. The line inside the boxplot represents the median of each land-use type. We plotted scatterplots by using estimated regression lines from our fitted models and removed the vanilla agroforest type (forest- or fallow-derived vanilla agroforests) from the model if solely vanilla yield had a significant effect on species richness (Fig. 1 panels 5, 7). We extracted the model fits with the allEffects function from the effect package^113,114. In all other cases (if no predictors were significant (e.g., Fig. 1 panel 3) or solely land-use history was significant as an additive or interactive effect (e.g., Fig. 1 panel 2) we retained the original model for plotting. In our graphs (Figs. 1 and 2), dashed horizontal lines are intercept-only linear models (lines are based on the mean of the distribution) and solid lines show statistically significant generalized linear regressions (P < 0.05). We show two colored lines as dashed lines if land-use history was significant as an additive term but no relationship of species richness with vanilla yield existed, or show them as solid lines if the effect of vanilla yield was moderated significantly by land-use history (Fig. 3).

To investigate differences in species composition among land-use types, we computed a permutational multivariate analysis of variance (PERMANOVA) using the adonis function of the vegan package¹¹⁵. Also, we used the pairwise.adonis function of the pairwiseAdonis package with false discovery rate correction to test for differences between land-use types¹¹⁶. Prior to PERMANOVA, we tested if homogenous dispersion existed among land-use types by using the betadisp function and permutest function of the vegan package (PERMDISP test)¹¹⁵. Our results show that heterogeneous dispersion significantly affected the differences in species composition for trees, birds, amphibians, butterflies, and ants (Table S14), thus differences are not only explained by the location of centroids but also by differences in dispersion between land-use types. We computed species composition by using non-metric multidimensional scaling (NMDS) with Jaccard dissimilarity distance. To compute the NMDS for amphibians, we excluded all plots where no amphibians occurred.

We generated sample-size-based rarefaction and extrapolation curves of species diversity of each taxon across land-use types to assess whether the curves reached an asymptote, indicating sampling completeness. To do so, we used ggiNEXT function of the iNext package¹¹⁷ with the following settings: type = 1, datatype = incidence_raw, q = 0 (species richness), endpoint = 20. Furthermore, we computed the sampling coverage in the percentage of each taxon across land-use types by using the function iNext with the same settings. We calculated gamma species diversity for each land-use type by summing up all species found across 10 plots and therefore excluded 10 out of 20 fallow-derived vanilla agroforests randomly (Table S18).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com