Drivers of migrant passerine composition at stopover islands in the western Mediterranean

Study islands and bird data

Systematic ringing in spring on Mediterranean islands has been promoted by the Piccole Isole project since 1988²⁶. Standard methods of the project involve ringing between 16th April and 15th May attempting to include the peak of the spring passage of long-distance migrants. Ringing is performed from dawn to nightfall using a constant number of nets within ringing stations placed at stable sites located at representative habitats in each island (Supplementary Table S1). The use of tape-lures is not allowed. We have compiled ringing data for all the Spanish Mediterranean islands that have been applying this methodology, with the exception of Mallorca and Menorca where the ringing stations were located in wetlands and captured a large percentage of local birds (Fig. 2, Table 1). The nine study islands are spread along a south-west to north-east gradient and, with the exception of Columbrets, they are distributed in pairs of similar longitude but different latitudes (Fig. 2). Ringing stations have been operating over a variable number of years (5–27 years), with the maximum number of ringing stations operating at the same time occurring between 2003 and 2010. To include between-year variation on islands that started ringing campaigns more recently we used data from the years 2003–2018.

The ringing period within each spring also varied in most islands, owing to funding or logistic limitations; thus, to reduce the possible effects on migrant composition we only used data from the standard period of the Piccole Isole project and from years that included at least one week of ringing in the fortnight of each month within this interval. This procedure excluded the use of some years for several islands, and the final number of data years for islands ranged between 5 and 16 (Table 1).

We used only data for trans-Saharan nocturnal migrant passerines, which form the bulk of species ringed on Mediterranean islands during the standard period. The standard ringing period only covers the tail end of the short-distance migrants’ passage; thus, these species were excluded as their contribution to composition of migrants could vary mainly due to between-year variation in migration phenology. Diurnal migrants, like hirundinids and fringillids, also represent a small fraction of birds ringed and may use different cues to select stopover islands. In addition, some of these species nest in some of the islands studied and birds ringed could include breeding birds. To avoid the distorting effect of species that are captured accidentally in very small numbers, we considered only the species that were ringed in at least five separate years, or on five different islands, which limited the species considered to 35 (Supplementary Table S2). This led to the exclusion of just two species (Ficedula semitorquata with three individuals ringed in two islands and Locustella luscinioides with one individual ringed in Aire island). In addition, we only considered the number of ringed birds, since the proportion of recaptures varies among islands, likely reflecting variation in the duration of stopovers²¹, which could bias the comparison of the patterns of migrant species composition.

Island descriptors

We obtained two groups of variables describing the characteristics of the study islands (Tables 2, 3): (1) Variables related to geographical location: latitude, longitude, straight distance and minimum distance to the North African coast, minimum distance to the closest large body of land (continent or large island) in any direction and to the closest large body of land situated in a southerly angle between SW and SE. (2) Variables related to the habitat characteristics of the islands: area, maximum altitude and Normalized Difference Vegetation Index (NDVI). We estimated NDVI from Landsat 8 Images taken during the standard ringing period in the years 2015 and 2016. Pixels containing shoreline were excluded and the average NDVI was calculated for the rest of the pixels.

Continental abundance data

Abundance estimates for western Europe were obtained from the European Red List of Birds²⁷. We used the mean of the minimum and maximum number of pairs estimated for the 27 EU Member States as a measure of continental abundance (Supplementary Table S2).

Data analysis

All analyses were done using R 3.6.1²⁸. We built a matrix of island-year x species containing the number of individuals of each selected species ringed in the study period in each island and year. Average number of individuals of each species ringed at each island was calculated and was used (log-transformed) as a dependent variable in a linear model with continental abundance (log-transformed), island and their interaction as predictors. This model was simplified using AICc as criteria to identify the best model.

To analyze variation of species composition, the matrix of island-year x species was transformed using the chord transformation²⁹ with the function decostand in the vegan package³⁰.

Using the function beta.div of the adespatial package³¹ we calculated beta diversity, including temporal and between-island variability (BD_I,T), as the total variance of the aforementioned transformed matrix (BD_Total in²⁹), which varies between 0 and 1 when chord distance is used. Considering that y_ijk is the chord transformed abundance of the species j in the island i and year k and (overline{{y }_{j}}) is the mean for species j in all islands and years altogether, then:

where N is the total number of samples. The function beta.div also provides an estimation of contribution of localities (LCBD) and species (SCBD) to beta diversity (Table 3). Yearly LCBD (log transformed because of skewed distribution) of each island were averaged and compared between islands using ANOVA and a post-hoc Tukey test.

We partitioned the above sum of squares in several ways. First, we calculated a beta diversity that considered only between-island variability, excluding temporal variability (BD_I), by averaging the chord transformed abundances of each species j in each island along study years (({overline{y} }_{ij})) and applying the same procedure, but using the number of studied islands (n):

Second, we calculated a beta diversity due to inter-annual variation of migrant composition within islands (BD_T) as:

where Y is the total number of study years and n is the number of studied islands (9). We also calculated a temporal beta diversity for each island i (BD_Ti) as the sum of squares due to variation within the island divided by the number of the island study years (Yi) minus 1:

Differences in temporal variability between islands could be due to different predominance of species that are more or less variable between years. To check this, we calculated Spearman’s rank correlation between the percentage of captures of each species in the total ringed on each island and BD_Ti and LCDB indices, for species present on all islands.

We tested for the existence of differences between islands in migrant species composition using Permutational Multivariate Analysis of Variance (PERMANOVA) using the function adonis2 in the vegan package. We performed a multivariate test of homogeneity of variances using the betadisper function (vegan package) with the adjustment for small sample bias, to test if temporal variability in species composition differed between islands. We made post-hoc comparisons between islands with False Discovery Rate (FDR) correction using the function pairwise.perm.manova of the package RVAideMemoire³².

To identify gradients in migrant species composition and the island characteristics that were associated with them, we employed Redundancy Analysis using the rda function (vegan package). We used the chord transformed matrix of species x island-year as a response matrix. We used two explanatory matrices, one including variables of geographical location and the other the variables related to habitat characteristics of the islands. We evaluated the relative importance of each group of variables to explain migrant species composition by performing a variation partitioning analysis, using the varpart function (vegan package). For that analysis, we followed the steps and R scripts recommended in³³.

Variables describing island characteristics were transformed using natural logarithms and collinearity within each group was evaluated with variance inflation factor (VIF)³⁴. All the habitat variables presented VIF < 3 and were retained in subsequent analyses, but some of the geographical location variables presented VIF values much larger than 10, so we removed the variable with the highest VIF and recalculated VIF values. This procedure was repeated until all variables presented VIF < 10 and led to the exclusion of minimum distance to Africa and latitude. Maximum VIF in the remaining variables was 2.1. We then performed a separate forward selection of variables within each group using the vegan function ordistep and all variables were selected.

A second redundancy analysis was performed with all the explanatory variables together to identify the particular variables most related to gradients of composition of migrants. Because we have repeated measures for each island, we fitted linear mixed models to identify the explanatory variables significantly related to each RDA axis, using the function glmmTMB (glmmTMB R package³⁵). In these models, scores of each RDA axis were dependent variables, each island descriptor was tested as an explanatory variable and the island identity was included as a random effect. To evaluate the species whose abundance was significantly related to RDA gradients, we calculated correlation coefficients of each species abundance (chord transformed) and RDA axis scores.

Models of species response

To evaluate the response of each study species to the geographical and habitat characteristics of the islands, we calculated the average number of individuals captured daily per 100 m of mist net in each study year in each island. However, the capture index from Columbrets Island is not comparable to the other islands because most birds were trapped in just three nets placed around a clump of isolated Opuntia cactus, which strongly attracted birds landing on this island. On the other islands, nets were set in several straight lines across the available habitat, and thus sampling was more uniformly distributed. The relationship between the percentage of total captures and the number of birds captured per 100 m of net is similar in all islands except Columbrets (Supplementary Fig. S1 online) and thus it was not possible to include Columbrets in the models for the captures/100 m net.

To check if excluding Columbrets from these analyses changed the results, the species response to island variables was modeled twice: using as a dependent variable the proportion of captures, including Columbrets, and using an abundance index calculated as the average number of birds captured daily per 100 m of mist net (excluding Columbrets). In both cases, we used generalized linear mixed models (GLMM, R function glmmTMB). In the proportion of captures models, we included the natural logarithm of the total number of birds ringed of the selected species in each island and year as an offset. For the captures/100 m net models, the offset was calculated as the natural logarithm of the result of multiplying the total length of nets (in hundreds of meters) by the number of days that the nets were operating in each spring. Island identity and year, considered as a factor, were included in the models as random effects. Each of the island variables included in the RDA, after excluding colinear variables, was included in turn as a fixed effect in these GLMMs. We considered three potential family distributions: Poisson and two parameterizations of the Negative Binomial distribution (nbinom1: variance increases linearly with the mean; nbinom2: variance increases quadratically with the mean). A null model including only the random effects was run with each family distribution and the one with the lowest AICc was selected. Then, linear and quadratic models including by turn one of the island descriptors were fitted and retained if it had an AICc lower than the null model. For each species we considered as equally plausible the models with ΔAICc ≤ 2. We used the performance_aicc function from the performance package³⁶ to calculate AICc. Significance of regression coefficients in quadratic models was checked by comparing quadratic and linear models with the anova function.

Effect of wing shape on island selection

We used the Kipp index as an indicator of wing pointedness³⁷. This index is calculated as the percentage of the distance from the tip of the first secondary feather to the tip of the longest primary feather (primary projection) of the wing length (Kipp = 100*primary projection/wing length)³⁸. Primary projection was measured with a transparent rule placed over the folded wing and wing length was measured as the maximum length with a butted ruler³⁸, both with a precision of 0.5 mm. Primary projection data are only available for Tabarca Island from 2018 onwards, so these measures were lacking for some of the species less frequently ringed in the study area. The primary projection of great reed warblers (Acrocephalus arundinaceus) and wood warblers (Phylloscopus sibilatrix) were obtained from a nearby ringing station in the mainland (Hondo Natural Park). The Kipp index of individuals were averaged to obtain mean and SD for each species. We also used the wing aspect ratio (WAR) available in³⁹ for most of the species studied.

We tested if mean Kipp index and wing aspect ratio differed between groups of species whose abundance was related to different RDA axes or selected island variables using phylogenetic generalized least-squares (PGLS) models⁴⁰ constructed using the R package ape⁴¹ and nlme⁴². These models allow to account for potential non-independence among species due to shared evolutionary history by estimating a phylogenetic signal index, Pagel’s λ⁴³, that measures phylogenetic dependence of observed trait data. For that aim, following⁴⁴ 2500 trees were downloaded from Birdtree web site⁴⁵ for the set species with Kipp index available and the same quantity of trees for the species with WAR estimates available. A consensus phylogenetic tree was built for each set of species using the R package phytools⁴⁶ (Supplementary Figs. S4 and S5 online). Iduna opaca was not available in Birdtree data and was substituted by the closely related Hippolais pallida (currently named Iduna pallida), from which it was split. Bivariate PGLS models including interaction were fitted and interaction was removed if not significant. Models were fitted with estimated Pagel’s λ and with λ fixed to zero (i.e. no phylogenetic signal). Likelihood ratio test was used to compare these models and if not significant the simplest model was selected⁴⁷.

Source: Ecology - nature.com