Study system and site
Antirrhinum majus L. (wild snapdragon; Plantaginaceae) is a short-lived perennial herb (Andalo et al. 2010). Snapdragons live in rocky cliffs and disturbed areas in the French and Spanish Pyrenees Mountains (Khimoun et al. 2013). Individuals often aggregate in small persistent patches (Jaworski et al. 2016), likely due to passive seed dispersal (Andalo et al. 2010; Khimoun et al. 2013).
Only large bumblebees (Hymenoptera) have enough strength to force their way through the flower lips to access snapdragon nectar (Guzmán et al. 2017) and pollinate snapdragons (Tastard et al. 2011, 2014). Bumblebees provide resources to their offspring by foraging on flowers in the area surrounding their nests (Goulson 2010). In our focal population, (Andalo et al. 2019) found that two of the three bumblebee species that pollinate snapdragons - Bombus terrestris L. and Bombus muscorum L. - shift their behaviour to nectar robbing across the flowering season. Specifically, they found that nearly 90% of visits produced between 05/25/2016 and 07/12/2016 were nectar robbing. Nectar robbing is the foraging behaviour in which a nectar reward is obtained without providing any benefit to the plant. Its effects on plant fitness can be neutral (Morris 1996) or detrimental, reducing the seed set (Goulson 2010; Irwin & Maloof 2002). Primary nectar robbing is defined as the first time a bumblebee robs nectar on a flower. It is easily traceable by observing the holes that bumblebees chew in the basal half of the snapdragon corolla (Andalo et al. 2019; Goulson 2010).
In our study system, we called the behaviour of feeding on flower ovules “frugivory”. This type of frugivory can have highly detrimental effects on plant fitness (Whitney & Stanton 2004). We observed larvae of Brachypterolus vestitus Kiesenwetter (snapdragon beetle; Coleoptera) inside 19% of sampled flowers, feeding on the ovules and possibly nectar. Often multiple individuals were found in each flower (up to 8). Adult B. vestitus lay their eggs inside snapdragon flower buds, where both adults and larvae feed on inner floral tissues of A. majus (Butler-Stoney 1988; Jelíken 2007; Wagner 1994). To identify the species, we sequenced the standard mtDNA COI barcode region from tissue of two larvae collected inside A. majus flowers, and six adult coleopterans collected in or on A. majus at our study site using standard methods (standard mtDNA COI sequences analysis performed by Brent Emerson’s laboratory at CSIC and IPNA in La Laguna, Tenerife, Spain). Then we used BLAST analysis against GenBank to identify the insects (NCBI 2022). The two frugivorous larvae and one adult specimen had > 99% match to B. vestitus (Supplementary Table 1 in Online Resource). Voucher specimens are deposited at Museo Nacional de Ciencias Naturales, Spain (larvae ID: MNCN_Ent 344534, MNCN_Ent 344535, adult ID: MNCN_Ent 356420).
Our study site is the Vall de Ribes, Catalonia, Spain. The study area is 3 km wide, between Planoles (42.3166N, 2.1042E) and Fornells de la Muntanya (42.3240N, 2.0472E) (Fig. 1). The site has been used in past research because it contains a flower colour cline, where yellow- and magenta-flowered A. majus meet and interbreed (Andalo et al. 2010; Whibley et al. 2006). Here, we ignore the effects of flower colour on plant-insect interactions because preliminary analysis indicates that nectar robbing and frugivory are not affected by flower colour (C. Baskett, unpublished data).
Data collection
We quantified the prevalence of nectar robbing and frugivory, two trophic interactions (i.e., interactions based on resource exploitation) involving snapdragons. We conducted a field survey from 06/16/2020 to 07/07/2020, during the peak flowering season. Censuses were taken between 9 am and 6 pm and conducted in a balanced manner to avoid temporal biases between locations. For 944 individuals (a haphazard subset of the population, due to limited time) we collected the uppermost flower and stored it in the fridge (4ºC) for a maximum of 48h before quantifying the trophic interactions. Primary nectar robbing sensu Goulson (2010) was indicated by a hole in the basal half of the flower corolla. Frugivory was indicated by the presence of B. vestitus larvae inside flowers.
In order to perform spatial analysis, we collected precise location data on the focal plants and all additional snapdragons we were able to access in the study area, for a total of 2333 individuals. Location data (accurate to within 3.7m at our site, Surendranadh et al. 2022) was recorded using GeoXT handheld GPS units (Trimble, Sunnydale, CA, USA). Using QGIS software (v. 3.10.14) (QGIS Development Team 2021), we manually moved the few points that were accidentally situated inside the road - due to Trimble measuring errors - to the nearest site outside it. We incorporated the point dataset inside an Observational Window (sensu Baddeley et al. 2016), which comprises the sampled area next to roads and paths where snapdragons flourish (14.86 ha), informed by 12 years of sampling snapdragons at this study site (N. Barton, unpublished data). We reprojected the spatial objects from the geographic coordinates system (EPSG: 4326) to the projected coordinates system (EPSG: 25831) in R (v. 3.6.3) (R Core Team 2021) with the package sf.
Aggregation patterns
Point patterns analysis (PPA) quantifies spatial aggregation, which can provide insight into how resources are exploited (Lancaster & Downes 2004; Rodríguez-Pérez et al. 2012) and if coexistence mechanisms exist (Brown et al. 2011). We used first and second-order PPAs to quantify the spatial patterns of the resource (snapdragons) and both trophic interactions (frugivory and nectar robbing). The test statistics of PPA were conducted with the R package Spatstat (Baddeley et al. 2016).
Before performing PPAs to test our hypotheses, we tested for homogeneity of the data, a crucial first step for deciding the appropriate PPAs. The intensity of points can be constant through space - called homogeneity - or vary by an intensity function, called inhomogeneity (Baddeley et al. 2016; Law et al. 2009). We tested for homogeneity using the quadrat counting and the Kolmogorov-Smirnov techniques. The quadrat counting technique is a chi-squared test defined from Pearson goodness-of-fit sensu Baddeley et al. (2016, pp. 163–168). We decided the optimal grid to perform the quadrat counting test was 18 x 12 cells. Each cell had a horizontal length of 185 m and a vertical length of 130 m. This resolution reflects a trade-off between snapdragon cluster size and bias. The Kolmogorov-Smirnov technique is a measure of discrepancy between cumulative distribution functions (Perry et al. 2006) sensu Baddeley et al. (2016, pp. 184–186). Both techniques showed significant intensity differences for the snapdragons and trophic interactions (Supplementary Table 2 in Online Resource), so our subsequent analyses assumed inhomogeneous data.
We performed first-order PPA, which asks whether the snapdragon aggregation is caused by spatially structured features and ecological gradients (Gimond 2021; Law et al. 2009). Additionally, we tested if the aggregation of the trophic interactions follows from the degree of aggregation of the snapdragons, which act in this case as the spatially structured features. We used the Hopkins-Skellam analysis, as recommended for inhomogeneous PPA (Baddeley et al. 2016, p. 259; Rubak 2019). Its formula is defined in Table 2a. Hopkins-Skellam index (HSI) = 1 is consistent with completely spatially randomness (CSR) pattern, HSI < 1 indicates clustering, and HSI > 1 indicates overdispersion. The Hopkins-Skellam analysis was performed using the hopskel.test function from Spatstat in R. To estimate the P-values, we proceed with a bootstrapping test. We ran 999 Monte-Carlo simulations of the Poisson null model, which is a CSR process inside the same Observational Window. The significance threshold was set at 0.05. Edge correction was not used because snapdragons were distributed next to roads and could not grow inside dense forests. These actual habitat boundaries are not an artefact of sampling, so edge correction is not needed (Law et al. 2009).
Out of 944 sampled plants, we found evidence of frugivory in 175 and of nectar robbing in 685. To estimate the degree of aggregation of the interactions independently from that of the resource itself, we proceeded to reshuffle the interaction presence - also called random labelling - across the 944 sampled plants 999 times and obtain significance with bootstrapping.
Next, we performed second-order PPA, which asks whether the exploitation patterns are explained by the resource aggregation and behavioural differences. Second-order analyses quantify the influence between pairs of individuals (such as exclusion by competition or aggregation by mutualistic interactions) (Gimond 2021; Law et al. 2009). The method used is the Pairwise Distances with Ripley’s K function (Baddeley et al. 2016, pp. 242–246; Dale 1999). The modification of the K function that applies to inhomogeneous data is the Kinhom (Baddeley et al. 2016, p. 243; Law et al. 2009; Perry et al. 2006). The inhomogeneous L function is analogous to the inhomogeneous K function. However, we used the L function because it stabilizes the variance and helps with visual assessment of the graph. It was computed with the Linhom function from Spatstat (Baddeley et al. 2016, p. 245) and it is defined in Table 2b. Then the L function is compared with the theoretical curve of the Poisson null model, defined in Table 2c (Baddeley et al. 2016, p. 206). The intensity was estimated with the density function from Spatstat (Baddeley et al. 2016, pp. 242–246). Some judgement is needed in choosing the appropriate bandwidth to smooth the intensity (Law et al. 2009; Rubak 2019). We chose a bandwidth value of 100 m because it created a density map large enough to not interfere with the natural aggregation processes of snapdragons and small enough to distinguish the local intensity differences. Then we proceed to bootstrap, generating 999 simulations of the Poisson null model of the L function (Baddeley et al. 2016; Swift et al. 2017). We then created simulation confidence envelopes (SCE) with an alpha level of 0.05. If the empirical L function lies outside the SCE, we rejected the CSR null hypothesis (Baddeley et al. 2016, pp. 390–394). When the L function lies above the SCE, it indicates clustering, and if it is below, it suggests an overdispersed pattern (Baddeley et al. 2016, pp. 390–394). We did the simulations for all plants and each plant-insect interaction.
In order to ask whether frugivory and nectar robbing occur closer together than expected, we used another form of second-order analysis to compare the spatial relationship between two point patterns (bivariate) and infer biological meaning such as coexistence or dependence between two species (Dixon 2014). To compare bivariate data, we used the inhomogeneous Lcross function. We transformed the Kcross function to Lcross function as described in Baddeley et al. (2016, pp. 594–596) to check how many type j events are nearby the i type (Table 2d). We used the Lcross function to examine whether the non-trophic interaction is more or less aggregated than expected. The Lcross function can also compare a point pattern with a subsample of itself with a reshuffling procedure (Baddeley et al. 2016; Fletcher & Fortin 2018). We used this feature to demonstrate that the interactions were spatially randomly sampled to avoid bias in the spatial analyses (Supplementary Fig. 1 in Online Resource), and to find the spatial patterns of trophic interactions isolating the resource effect.
Density of trophic interactions v. habitat patch size
To test whether resource exploitation depends on resource patch size, we used quantile regression analysis of interaction frequency in discrete resource patches. We defined patches with the Kernel Density Estimation (KDE) map tool in QGIS. KDE is a commonly used tool to determine hotspots - or patches in our case (Nelson & Boots 2008). No edge correction was utilised because (1) the KDE applies to the intrinsic dataset, not to the individuals that could lie outside the Observational Window and (2) we clipped the area of the KDE inside the Observational Window. We chose the Triweight method (Table 2e) because it creates sharper patches (Guidoum 2015). The KDE bandwidth estimation differs from that chosen in the second-order analyses. In this case, the objective was to create a sharp KDE to identify the patches of resources. The likelihood cross-validation bandwidth selection method, with the bw.ppl function from Spatstat, provided the optimal bandwidth for our dataset: 3.2m. Based on our visual perception, we considered that one plant per square meter is the minimum density to create a continuous snapdragon patch for the bumblebees and adult coleopterans. Therefore, we created the KDE contour map with the QGIS tool contour at a density of 1 ind./m2 and a pixel size of 0.1 m.
We compared the resource patch area with the number of trophic interactions detected inside. The abundance was transformed to decimal logarithm after adding 1 (Connor et al. 1997; Schooley & Wiens 2005). We then used quantile regression analysis because it supports heteroscedastic data - like ours - and is equivariant to monotonic transformations (Cade et al. 1999). We proceed as sensu Schooley and Wiens (2005), indicating the 5% regression quantile as the lower constraint, the 50% as the median and the 95% as the upper one. The slope and confidence intervals (CI) of the 95% quantile were compared with the slope of the Null Linear Regression (NLR). We calculated the NLR slope by dividing the number of snapdragons surveyed for plant-insect interactions by the total number of snapdragons across the dataset. If the NLR is between the 5% and 95% quantile regression lines, then our data is equivalent to the ideal free distribution hypothesis. If the NLR is significantly above the 95% quantile, it corresponds with the undermatching hypothesis. When the NLR is significantly below the 5% quantile, it matches the resource concentration hypothesis (Schooley & Wiens 2005). Simultaneously, if the 95% quantile had a significantly greater slope (CI 95%) than the 50% quantile, then the resource patch size is the active constraint of the trophic interaction among other factors, called a limiting-factor relationship (Cade et al. 1999; Schooley & Wiens 2005).
Resource accessibility in a heterogeneous landscape
To assess the accessibility to resources in a heterogeneous landscape, we computed the minimum distance between the location of each trophic interaction and six landscape features. Initially, we verified the spatial autocorrelation of the point-pattern data across four scales, employing the joint count test (Stevens & Jenkins 2000) conducted with the R package spdep. The scales used were 5, 20, 40, and 100 nearest neighbor points. This preliminary analysis confirmed that the trophic interactions were spatially autocorrelated.
To mitigate the influence of spatial autocorrelation and ensure minimal data loss while examining the association between trophic interactions and landscape features, we implemented a thinning process (AKA pruning) on the point spatial data, separately for frugivory and nectar robbing. We followed the procedure outlined by Assis (2020) to determine the minimum thinning required to eliminate spatial autocorrelation. First, we thinned the complete dataset to retain one individual every X meters in the surrounding area, where X ranged from 1 m to 30 m in 0.5 m increments (yielding 59 subsamples). Next, we evaluated spatial autocorrelation in the subsamples at the same four spatial scales (5, 20, 40, and 100 nearest neighbor points) using the joint count test. For subsequent analysis of landscape features, unambiguously avoiding autocorrelation across multiple scales, we selected the minimum thinning interval without significant spatial autocorrelation at any of these four scales. To account for potential variation introduced by randomness in the thinning process, we conducted bootstrapping, repeating the thinning and autocorrelation tests one hundred times. Finally, we calculated the mean of the minimum thinning distances obtained through bootstrapping. The resulting minimum thinning distance necessary to eliminate significant autocorrelation was 18.5 m for nectar robbing and 3 m for frugivory, reducing the datasets to 155 and 466 plants, respectively.
We extracted the landscape features from the soil use layer of the Institut Cartogràfic i Geològic de Catalunya (Cobertes del sòl. ICGC 2018) and categorized them as bare soil, fields, farming, forests, gardens and urban areas (Supplementary Table 4 in Online Resource ). We computed the minimum distance between the location of the trophic interactions and the landscape features using the QGIS geoalgorithm v.distance. We then generated a logistic GLM to compare the presence or absence of each trophic interaction (previously thinned to eliminate spatial autocorrelation) with the minimum distance to each landscape feature, thus it is independent of plant location.
Table 1
Summary of the hypotheses and spatial analyses used
spatial scale (m) | interactions | prediction: exploitation patterns are explained by... | analysis |
0-120 | trophic | behaviour (tendency of insects to aggregate) | point patterns |
non-trophic | the spatial pattern of the other interaction |
1–220 | trophic | the resource patch size | quantile regression |
0–1500 | trophic | the landscape structure | distance GLM |
Table 2
Summary of the formulas used
| metric | description | formula |
a. | HSI | Hopkins-Skellam index* | \(={\sum }_{i = 1}^{m} { d}_{i}^{ 2} / {\sum }_{j = 1}^{m} { e}_{j}^{ 2}\) |
b. | \({\text{L̂} }_{inhom}\left(r\right)\) | Transformation of the K function to the inhomogeneous L function * | \(=\sqrt{ \mathbb{E}\mathbb{ }\left[\mathbb{ }{\sum }_{{x}_{j}\in \varvec{X}}^{ } \frac{1}{\lambda \left({x}_{j}\right)} 1\right\{0<\left|\right|u - x{ }_{j} \left|\right| \le r\left\} \right| u \in \varvec{X} ] /\pi }\) |
c. | \(\text{L}{ }_{inhom, pois}\left(r\right)\) | Poisson null model of the L function * | \(\equiv r\) |
d. | \({\text{L̂}}_{cross inhom}\left(r\right)\) | Inhomogeneous Lcross function * | \(=\sqrt{ {{\lambda { }_{j}}^{-1}}_{ }\mathbb{E} \left[t ( u, r, {X}^{\left(j\right)} ) 丨u ϵ {X}^{\left(i\right)} )\right] / \pi }\) |
e. | \(K(x;6)\) | Triweight method for the Kernel Density Map ✝ | \(= \frac{35}{32} {(1-{x}^{2}{)}^{3}{1}_{\left(\right|x|\le 1)} }^{ }\) |
HSI: i and j are the index of a vector of the nearest-neighbour points, di is the nearest-neighbour distances for m randomly sampled data points, and the empty-space distances ej for an equal number m of randomly sampled spatial locations. \({\text{L̂}}_{inhom}\left(r\right)\): each point xj from the point process X is weighted depending on the function of the intensity λ that is subject to the location u on a specific radius r. \(L{ }_{inhom, pois}\left(r\right)\): \(\equiv\) represents identical to. \({\text{L̂}}_{cross inhom}\left(r\right)\): \({X}^{\left(j\right)}\)is the sub-process of points of type j, with intensity λj, and t(u, r, X) is the count of r-neighbours. \(K(x;6)\): x in a given point and 6 is the maximum derivative of Kernel for the Triweight method.
* Baddeley et al. (2016)
✝ Guidoum (2015)