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.

Study site and ecosystem engineers

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: Arenaria tetraquetra ssp. amabilis (Bory) H.Lindb. (Caryophyllaceae) and Hormathophylla spinosa L. (Brassicaceae). The engineer Arenaria tetraquetra 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). Hormathophylla spinosa 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).

Experimental design

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: Anthyllis vulneraria ssp. pseudoarundana H. Lindb. (Fabaceae), Chaenorrhinum glareosum (Boiss. & Reut.) Willk. (Plantaginaceae), Crepis oporinoides Boiss. (Asteraceae), Jasione amethystina Lag. & Rodr (Campanulaceae), Leontodon boryi Boiss. (Asteraceae), Leucanthemopsis pectinata (L.) López & Jarvis (Asteraceae), Lotus corniculatus subsp. glacialis (Boiss.) Valdés (Fabaceae) and Sideritis glacialis 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 A. tetraquetra and 14 blocks for H. spinosa, were established within a relatively homogeneous area of about 1 ha (n = 84).

Plant–visitor interactions

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 (n = 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.

Data analysis

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, A. tetraquetra and H. spinosa). 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_EE and Y_AP are flower visitation rates of ecosystem engineers and associated plants when blooming in association (i.e., treatment 3), and M_EE and M_AP 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_{engineer present} / 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.

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