2.5 Statistical analysis
Changes in pollinator species richness were evaluated in response to green area fragmentation (i.e. , the variable Edge Density) and flower richness (i.e. , the number of flowering species per site). To do this, a Generalised linear mixed model (GLMM) regression (glmmTMB R package; Magnusson et al., 2017) with Poisson distribution was used, with island included as a random effect. The flower richness was included as a predictor along with the edge density, since it could represent an important local driver of pollinator richness (Blüthgen & Klein 2011). The same variables (edge density, and flower richness) were used along with the network size as predictors of change of the Connectance network index. Network size, calculated as the product between the number of insects and plants included in the networks for each site, was also included as a predictor in the model to account for its effect on Connectance variation (as in Biella et al., 2020). In this case, a GLMM with beta distribution and island included as a random effect was used. Changes in individual pollinator degree were evaluated in response to green area fragmentation and flower richness. The effects of these covariates were evaluated in interaction with the pollinator species identity to highlight differences among the considered pollinator species. A GLMM with Poisson distribution was used, with sites nested in the island as a random effect. Variation in the pollination efficiency was evaluated using the pollinator richness and the Connectance as covariates. Moreover, the plant degree (mean number of pollinator species interacting with each plant species considered in pollination efficiency analysis) was calculated from DNA metabarcoding data to estimate the mean plant generalism and included as model covariate. The role of these covariates was evaluated in interaction with the plant species identity, to highlight differences among the investigated plants. A GLMM with negative binomial distribution was used to account for overdispersion. Also in this case, the site nested in island was included as a random effect.
All the analyses were performed with R (version 3.6.1; R CoreTeam 2019). Predictor significance was evaluated through a log likelihood ratio test (P < 0.05). The Vif function of the car package (i.e. , Variance Inflation Factor with an exclusion threshold of 3 (Zuur, Ieno, & Smith, 2007)) was used to exclude collinearity among variables. In all cases, the final models were obtained by removing the variables that did not improve the model fit through backward stepwise regression based on second-order Akaike Information Criterion (AIC) (Zuur et al., 2009) calculated with the package MuMIn.