DOI: 10.5061/dryad.dz08kprw0
title: Data and R code for: An experimental approach to assessing the impact of ecosystem engineers on biodiversity and ecosystem functions

description: Plants acting as ecosystem engineers create habitats and facilitate
 biodiversity maintenance within plant communities. Furthermore,
 biodiversity research has demonstrated that plant diversity enhances the
 productivity and functioning of ecosystems. However, these two fields of
 research developed in parallel and independent from one another, with the
 consequence that little is known about the role of ecosystem engineers in
 the relationship between biodiversity and ecosystem functioning across
 trophic levels. Here, we present an experimental framework to study this
 relationship. We combine facilitation by plants acting as ecosystem
 engineers with plant–insect interaction analysis and variance partitioning
 of biodiversity effects. We present a case-study experiment in which
 facilitation by a cushion-plant species and a dwarf-shrub species as
 ecosystem engineers increases positive effects of plant functional
 diversity (ecosystem engineers and associated plants) on ecosystem
 functioning (flower visitation rate). The experiment, conducted in the
 field during a single alpine flowering season, included
 the following treatments: 1) removal of plant species associated with
 ecosystem engineers, 2) exclusion (covering) of ecosystem engineer
 flowers, and 3) control, i.e. natural patches of ecosystem engineers and
 associated plant species. We found both positive and negative
 associational effects between plants depending on ecosystem engineer
 identity, indicating both pollination facilitation and interference. In
 both cases patches supported by ecosystem engineers increased phylogenetic
 and functional diversity of flower visitors. Furthermore, complementarity
 effects between engineers and associated plants were positive for flower
 visitation rates. Our study reveals that plant facilitation can enhance
 the strength of biodiversity–ecosystem functioning relationships, with
 complementarity between plants for attracting more and diverse flower
 visitors being the likely driver. A potential mechanism is that synergy
 and complementarity between engineers and associated plants increase
 attractiveness for shared visitors and widen pollination niches. In
 synthesis, facilitation among plants can scale up to a full network,
 supporting ecosystem functioning both directly via microhabitat
 amelioration and indirectly via diversity effects.
description: <i>Study site and ecosystem engineers</i> The study was carried out in a dry, Mediterranean alpine ecosystem in Sierra Nevada, Spain (Loma del Mulhacén, 3200 m a.s.l.). This is a well-studied model ecosystem characterized by vegetation patches that are structured by two sparsely-distributed plants that act as ecosystem engineers: <i>Arenaria tetraquetra</i> ssp. <i>amabilis</i> (Bory) H.Lindb. (Caryophyllaceae) and <i>Hormathophylla spinosa</i> L. (Brassicaceae). The engineer <i>Arenaria tetraquetra</i> is a cushion plant with compact, small and dense branches (Blanca et al. 2011, Schöb et al. 2012) and white actinomorphic flowers visited by more than twenty-five pollinator species (Gómez et al. 1996; Losapio et al. 2019). Associated plant species grow and flower on top of its canopy (Fig. 1). <i>Hormathophylla spinosa</i> is a dwarf shrub with spiny branches and small laves (Blanca et al. 2011) and has copious actinomorphic purple flowers visited by more than twenty pollinator species (Losapio et al. 2019). Associated plant species grow and flower underneath its canopy (Fig. 1). <i>Experimental design</i> After surveying the plant community (Appendix S1: Method S1), we selected a pool of plant species for carrying out the removal experiment, which consisted of the eight most-frequently associated plant species: <i>Anthyllis vulneraria</i> ssp. <i>pseudoarundana</i> H. Lindb. (Fabaceae), <i>Chaenorrhinum glareosum</i> (Boiss. & Reut.) Willk. (Plantaginaceae), <i>Crepis oporinoides</i> Boiss. (Asteraceae), <i>Jasione amethystina</i> Lag. & Rodr (Campanulaceae), <i>Leontodon boryi</i> Boiss. (Asteraceae), <i>Leucanthemopsis pectinata</i> (L.) López & Jarvis (Asteraceae), <i>Lotus corniculatus</i> subsp. <i>glacialis</i> (Boiss.) Valdés (Fabaceae) and <i>Sideritis glacialis</i> Boiss (Lamiaceae). The removal experiment consisted of the three following treatments: 1) removal (clipping aboveground parts) of all associated plants within and around randomly selected individuals of ecosystem engineers, 2) flower exclusion (covering with stones) of ecosystem engineers nearby associated plants, and 3) no removal (i.e., patches with naturally co-occurring ecosystem engineers and associated plants). In this way we ensured to remove the effects of ecosystem engineer flowers, ignoring other abiotic effects, without compromising the integrity of examined organisms and ecosystem (a protected National Park). Each treatment consisted of 20 x 20 cm standard plots, a size that reflects the fine spatial scale of plant facilitation in this alpine environment (Schöb et al. 2012, Schöb et al. 2013) as well as the intimacy of plant–pollinator interactions (Thomson & Chittka 2001, Losapio et al. 2019). Given the fine spatial resolution of our experiment, we were interested in assessing the insect response at the community scale (i.e., flower visitors) rather than at the landscape scale (i.e., species pool). We followed a randomized-block design, where each block included the three treatments for each ecosystem engineer species. The distance between plots within blocks was 50–100 cm, which was considered sufficient for undisturbed insect choice (Thomson & Chittka 2001). We kept associated plant species composition and diversity as similar as possible across blocks (Losapio et al. 2019). In total, 28 blocks, 14 blocks for <i>A. tetraquetra </i>and 14 blocks for <i>H. spinosa</i>,<i> </i>were established within a relatively homogeneous area of about 1 ha (<i>n </i>= 84). <i>Plant–visitor interactions</i> We surveyed flower visitors over the course of an entire alpine flowering season, from beginning of flowering to petal falling, during July 2015 (Losapio et al. 2019). Species interactions were surveyed by directly observing and sampling flower visitors to flowers (plant species) in each plot. Blocks were randomly sampled between 10 am and 5.30 pm for 20 min (5 min break between rounds). We followed a matched ‘pairs’ design and sampled three plots within each block at the same time in order to exclude variation within blocks due to changing weather conditions and flower visitor activity. Then, blocks were sampled following a stratified random design. Each day we sampled 14 blocks, completing sampling of all 28 blocks over two consecutive days each time. We ensured that each block was sampled at each time of the day, in a random order, for 6–9 times (blocks were dismissed once petals started falling). In total, we performed 204 sampling rounds over 15 sampling days (from 6 July 2015 – 27 July 2015), covering the whole flowering phase of examined alpine plants (Appendix S1: Tab. S1). Data were subsequently pooled for each plant species at the plot level (<i>n</i> = 84). Flower visitors were identified at the species level whenever possible, otherwise to genus or family (Appendix S1: Tab. S1). Afterwards, we classified flower visitors on the basis of life history by expert knowledge into the following trophic guilds (i.e., functional groups): pollinators, herbivores, scavengers, predators and parasitoids (Appendix S1: Tab. S2). These functional groups reflect the feeding behavior of the adult stage, the one that is observed interacting directly with flowers, rather than the entire life cycle. For instance, butterflies were considered as pollinators (adult stage) rather than as herbivores (larval stage). Since we exclusively observed visits to flowers, individuals and species classified as pollinators may not all do pollination successfully and non-pollinators may do some pollination. <i>Data analysis</i> We first tested the effects of the two plant species acting as ecosystem engineers on (1) the density and (2) diversity of associated plant species. To model the density of associated plants across microhabitats, we used a generalized linear model (GLM) with Poisson distribution with the abundance of individual plant species as dependent variable and the type of microhabitat as independent variable (categorical with open, <i>A. tetraquetra </i>and<i> H. spinosa</i>). To model the diversity of associated plants across microhabitats, we first fitted species richness as a function of the number of individuals using a Poisson GLM. Then, we modelled the residuals of this model as a function of microhabitat using a linear model (LM), therefore accounting for and removing direct effects of density on diversity. To answer the second question, we tested the effects of the presence of ecosystem engineers on the number of flower visitors to associated plant species. We used a Poisson mixed-effects GLM with visitor abundance as a function of ecosystem engineer presence (categorical variable with two levels), ecosystem engineer identity (categorical variable with two levels) and their interaction; the abundance of flowers of associated plants was fitted first as further covariate and blocks were random effects. We then tested the effects of plant removal treatments on (1) visitor phylogenetic diversity, and (2) visitor functional diversity. To do so, we first constructed a phylogenetic tree with the relationships among all 75 sampled visitors taxa (Appendix S1: Fig. S1) using NCBI taxonomy (https://www.ncbi.nlm.nih.gov/Taxonomy/CommonTree/wwwcmt.cgi). We then measured phylogenetic structure using abundance-weighted mean phylogenetic distance (MPD) among all visitor in a plot (Kembel et al. 2010). We considered the standardized effect size of MPD by comparing the observed values with the mean and standard deviation of MPD from 99 random communities generated by randomizing the identity of visitor taxa occurring in the plots (Kembel et al. 2010). We used linear mixed models (LMMs) to test the response of phylogenetic diversity to removal treatments (categorical variable with three levels), ecosystem engineer species (categorical variable with two levels), and their statistical interaction (fixed effects); plant species richness (i.e., number of subordinate plant species) and flower density (i.e., number of flowers per plot) were included as fixed effects and fitted first; blocks were random effects. Furthermore, we tested whether visitor functional diversity differed among removal treatments and ecosystem engineer species. We calculated the diversity of visitor functional groups using the Shannon index of guild abundance. Visitor diversity was standardized by flower abundance per plot. We used the same LMM structure to test the response of functional diversity to removal treatment as we did for the previous analysis of phylogenetic diversity. We proceeded with measuring the diversity effects of plant–plant facilitation on ecosystem functioning. Flower visitation rate was used as a measure of ecosystem functioning given the biological functions of flower visitors for pollination and hence for fruit production (Winfree et al. 2018; Woodcock et al. 2019). Flower visitation rate was calculated as the number of individual visitors on flowers divided by number of flowers at the plot level over the entire season. We used the additive partitioning of biodiversity effects (Loreau and Hector 2001), an approach that decomposes diversity effects on the basis of proportional deviation from expected values into two components: complementarity effects and selection effects. The complementarity effect measures the deviation in the average flower visitation rate of ecosystem engineers and associated plant species when blooming alone (i.e., treatments 1 and 2) or in association with each other (i.e., treatment 3). Positive complementarity effects occur with partitioning of flower visitors and indirect facilitation, while negative complementarity effects indicate indirect competition or chemical interference (Loreau and Hector 2001; Barry et al. 2019). The selection effect measures the covariation in the average flower visitation rate of ecosystem engineers and associated plant species between blooming alone (i.e., treatments 1 and 2) and in association (i.e., treatment 3). Positive selection effects occur if species with higher-than-average flower visitation rates are also dominating community-level flower visitation rates in mixtures, and vice versa for negative selection effects. We considered ecosystem engineers and associated plant species as two different functional groups (i.e., two functional diversity levels) in light of their distinct roles in the ecosystem (Ellison et al. 2005; Schöb et al. 2012; Losapio et al. 2019). We therefore measured the complementarity effect and the selection effect. Y<sub>EE</sub> and Y<sub>AP</sub> are flower visitation rates of ecosystem engineers and associated plants when blooming in association (i.e., treatment 3), and M<sub>EE</sub> and M<sub>AP</sub> are flower visitation rates of ecosystem engineers and associated plants when blooming separately (i.e., treatments 1 and 2). When flower visitation rate was equal to zero (n = 4), in order to avoid meaningless results (i.e., diversity effects equal to infinite) we considered those infinite values as the maximum diversity effect measured times two (i.e., CE = 0.160 and SE = 0.656; Appendix S1: Table S3). Excluding those data did not change the results. We tested each effect separately with LMs using the measure of diversity effect as response and the engineer species as predictor variables. Finally, we tested the effects of ecosystem engineers on the proportion of visitor species gained on associated plants in the presence of ecosystem engineers. This measure of species gain was calculated separately for each ecosystem engineer assemblages according to Diamond (1969) as species gain = appearance<sub>engineer present</sub> / total species Then, we used an LM with species gain (continuous variable) as a function of ecosystem engineer presence (categorical variable with two levels), ecosystem engineer identity (categorical variable with two levels) and the interaction term between them. The significance of predictors was tested both in terms of model fit and in terms of explained variance using type-II ANOVA (Fox & Weisberg 2011). Data analysis was conducted in R 3.6.0 (R Core Team 2019) using the packages ‘vegan’ (Oksanen et al. 2018) for functional diversity, ‘picante’ (Kembel et al. 2010) for phylogenetic diversity, ‘lme4’ (Bates et al. 2015) for mixed-effects models, ‘car’ (Fox and Weisberg 2011) for ANOVA tests and ‘emmeans’ (Lenth 2019) for least-square means contrasts and Cohen’s standardized effect size.
description: R script for the study "An experimental approach to
 assessing the impact of ecosystem engineers on biodiversity and ecosystem
 functions" published in Ecology by Gianalberto Losapio and
 colleagues losapio_ecology2020.R = Code to reproduce
 the analyses and plots losapio_ecology20.RData = R
 working space #######
 insect.record.csv = Sampling dataset of plant-insect
 interactions Fs = microhabitat species identity / Tr =
 treatment  ####### sndata.csv =
 Dataframe with summary data of plant-insect interactiomns
 ####### diversity_3200.csv = Sampling dataset
 of plant-plant facilitation PlotID = ordinal plot id
 number / Microhabitat = cushion plant species / cushion_length = length of
 cushion plant / cushion_width = width of cushion plant / species =
 associated plant species / indID = ordinal plant id number / cover =
 estimated coverage / demography = life history stage / phenology =
 phenological stage / nr.buds = number of buds / nr.flowers = nuber of
 flowers / nr.fruits = number of fruits
 ####### divdf.csv = Dataframe with summary
 data of plant-plant facilitation plotid = ordinal plot
 number / microh = microhabitat / nsp = plant species richness / nind =
 adult abundance / njuv = juvenile abundance / nsedl = seedling
 abundance ####### phyltree.phy =
 Phylogenetic tree of sampled insects