Study areas
We studied the bird communities of inland wetlands located in the pre-Alpine belt of Lombardy Region (Northern Italy, Fig. 1), one of the most densely populated area of Europe, with a mean population density of 421 people/km2 (ISTAT 2018). Inland wetlands of Lombardy are severely threatened by anthropic pressure and measures to foster their long-term conservation should be a priority for regional and national governments, also considering the EU requirements. Our study was carried out in 21 inland wetlands within the ‘Natura 2000’ network (Fig. 1, ‘study sites’ thereafter). Most of the sites are rather small wetlands (mean extension in ha ± SD 371.58 ± 441.56; min 20.23, max 1610.64) of post-glacial origin (sites 4–14 and 19–20 in Fig. 1), while others originate as shallow water areas along major rivers (i.e. river Adda, sites 16–18) or lakes (i.e. Lake Maggiore, sites 1–3; Como Lake site 15). The study sites are scattered over an area of ~ 7800 km2, ranging between 46°12′–45°38′N and 8°36′–10°01′E (Fig. 1), at an average elevation of 257.95 (± 82.89 SD) m a.s.l. (range 171–537). The climate is temperate continental, with a mean annual temperature of 11.1 °C and mean annual precipitation of 1236 mm (values referred to Montorfano Lake, nearly occupying the mid-point of the study area, source: Merkel 2017). The habitat of the study areas is mainly characterised by the occurrence of patches of Phragmites australis, in some cases forming wide continuous reedbeds up to 290 ha (site 5, Fig. 1), but generally reduced to thin strips of a few meters of width bordering open water extents. Together with the reedbeds, the study areas host patches of wetland habitats protected by the European Union’s Habitat Directive as the “alluvial forests with Alnus glutinosa and Fraxinus excelsior (Alno-Padion, Alnion incanae, Salicion albae)” (EU code: 91E0) and the “Natural eutrophic lakes with Magnopotamion or Hydrocharition-type vegetation” (EU code: E3150).
Target species
We performed our study by surveying different sets of species at two different spatial scales, i.e. the ‘study site’ scale (with Natura 2000 sites being the sampling units, n = 21), and the ‘point count’ scale (with sampling points surveyed during point counts as sampling units, n = 75). At the site scale, we checked the occurrence as breeding taxa of 10 species of conservation priority in the EU (i.e. included in the Annex I of the Birds Directive) and known to breed in some of the target sites: Western Marsh Harrier (Circus aeruginosus), Purple Heron (Ardea purpurea), Ferruginous Duck (Aythya nyroca), Kingfisher (Alcedo atthis), Black-crowned Night Heron (Nycticorax nycticorax), Little Bittern (Ixobrychus minutus), Black Kite (Milvus migrans), European Honey Buzzard (Pernis apivorus), Black Woodpecker (Dryocopus martius) and Lesser Spotted Woodpecker (Dendrocopos minor; the last not included in Annex I of the Birds Directive but of conservation concern at the local scale). We also checked for the occurrence of further three species listed in Annex I of the Birds Directive, but these eventually resulted completely absent (i.e. Moustached Warbler Acrocephalus melanopogon) or too scarce to be analysed in the study (Little Crake Porzana parva and Spotted Crake Porzana porzana occurred with only 1–2 pairs in the entire study area).
At the point count scale, we checked the occurrence of four bird species, not listed in the Annex I of the Birds Directive: Eurasian Reed Warbler (Acrocephalus scirpaceus, referred to as Reed Warbler hereafter in the text), Marsh Warbler (Acrocephalus palustris), Great Reed Warbler (Acrocephalus arundinaceus), passerines belonging to Acrocephalidae; Water Rail (Rallus aquaticus), belonging to Rallidae. The three passerines are insectivorous, long-distance migrants with diurnal habits and all breed in the reedbed, where they occur along a habitat gradient, with Great Reed Warbler occupying permanently flooded and monospecific reedbeds, and with Marsh Warbler at the other end of the gradient, inhabiting drier reedbeds mixed with bushes (Cramp and Simmons 1998; Kennerley and Pearson 2010). The Water Rail is a mainly nocturnal bird, which feeds on invertebrates found in the mud (Cramp and Simmons 1998) and year-round occurring in the area, although the migratory behaviour of the local populations is still poorly known. Overall, the habitat requirements of these four species taken together represent fairly well the broader spectrum of ecological niches occurring in reedbeds. Therefore, their presence can be taken as an indicator of the overall ecological status of a reedbed (see also Brambilla and Jenkins 2009; Mortelliti et al. 2012).
Environmental variables
Environmental variables were collected at the two different spatial scales of Natura 2000 sites and point counts respectively.
At the Natura 2000 sites scale, we estimated land-cover within each study site, thanks to a land-use/land-cover (LULC) map available from Lombardy Region and updated to 2015 (DUSAF 2015; resolution 1:10,000). Even if this LULC map, based on photointerpretation, is very accurate and updated, we field-validated the boundaries of the habitat types and manually corrected them when needed in QGis (v 2.4.0, Quantum GIS Development Core Team 2018). Fifty-two LULC categories were included in our study areas: we then merged them into 9 habitat categories (see Additional file 1: Table S1 for the complete list of original habitat types): (1) ‘reedbeds/mires’ (mainly reedbeds of P. australis, but including mires and Magnocaricetum assemblages); (2) ‘inland open waters’ (natural and artificial ones); (3) ‘non-vegetated areas’ (including buildings); (4) ‘arable lands’; (5) ‘meadows’; (6) ‘riparian vegetation’ (low shrubby vegetation, at least occasionally flooded); (7) ‘transitional shrubland’ (high shrubs); (8) ‘deciduous woods’ (mainly Salix alba and A. glutinosa forest patches); (9) ‘wood crops’ (poplar plantations, orchards).
At the point count scale, we mapped the habitat composition within a 100 m-buffer from each point (measured in QGis). Successively, we estimated the proportion of each habitat type within the buffer and eventually converted this proportion into an extent of habitat cover (ha). At this step, we considered 8 habitat categories: Open shallow water; P. australis reedbed; trees (> 5 mt height); shrubs (< 5 mt); meadows (< 1 mt); Magnocaricetum complex (typical wetland assemblages of tall grasses, up to 1.5 m, dominated by Carex spp., mostly C. elata, C. pendula and C. acuta within the study area); bare soil; other LULC types. Furthermore, since each point was visited for the environmental survey at least twice in the season, we added a categorical variable indicating if the reedbed fraction of the buffer was permanently inundated during the season (2), occasionally inundated (1), or permanently dry (0).
Survey methods
Surveys aimed to detect the species occurrence were also performed at the two different spatial scales of sites and point count, respectively (the sets of target species were mutually exclusive for the two spatial scales). To check the occurrence of the 10 species of conservation concerns, we followed the species-specific guidelines described in the ‘Region Lombardy Plan for monitoring of avian species of Annex I’ (Fondazione Lombardia per l'Ambiente 2015), adopting the adequate survey method for each target species and adapting the methodology to the two survey scales.
Natura 2000 site survey scale
At the site scale, the presence/absence of Lesser-spotted and Black Woodpecker was assessed by means of linear transects carried out along existing pathways: within each wetland, an observer walked at a slow pace along the transect in early morning in March–April, when woodpeckers’ territorial activity peaks (Cramp and Simmons 1998). The presence/absence of Marsh Harrier, Purple Heron, Night Heron, Black Kite and Honey Buzzard was assessed by means of morning observation sessions (duration: 2–4 h according to site characteristics), starting 1 h after sunrise, in May–June. Whenever possible, vantage points allowing a comprehensive view of the wetland were chosen for these observational sessions. Wood patches of the studied areas are often of limited extent, which makes easier—compared to other environmental contexts—to assess the occurrence of breeding tree-nesting raptors. Furthermore, to assess the occurrence as breeding species of Black Kite and Honey Buzzard, we also collected anecdotal observations from local birders that frequently visited the areas during our survey period. The occurrence of breeding Ferruginous Ducks was assessed in June, when chicks are still unable to fly, by means of observations from shorelines, coupled with kayak surveys of reedbed borders. Each site was monitored with the help of 20 × 60 telescopes and 10 × 42 binoculars, at least twice during the period. Opportunistic observations of Kingfishers, realized during the previously described dedicated sessions, were registered from March to June for all the surveyed sites and the presence of the species as a breeder was therefore established based on these data. The occurrence of breeding Little Bitterns was assessed in May–June, by means of night-transects (from dusk to the 2–4 following hours) along the pathways closest to reedbeds, listening to males’ territorial call.
Point count survey scale
At the point count scale, the occurrence of the three warblers species and of Water Rail was established by the means of point counts realized within 2 h after sunrise for warblers and at dusk (from 1 h prior to sunset to 2 h after) for Water Rail. Point counts have been reported as the most time-efficient way to estimate presence/absence and population density in elusive species (Sutherland et al. 2004). Census based on listening is the common method used to estimate presence/absence and abundance for reedbed-dwelling warbler species (i.e. Mortelliti et al. 2012; Sozio et al. 2012; Ceresa et al. 2016) and proved to be effective for rallids too (Polak 2005; Brambilla and Jenkins 2009; Jedlikowski et al. 2014). In the point counts aimed to detect water rail presence, we broadcast a registered call (1-min length) to stimulate a possible territorial response (e.g. Brambilla and Jenkins 2009), followed by 3 min of listening, again a 1-min stimulation and 3 further minutes of listening, for a total length of 8 min. Increasing the duration of a listening point from the standard 5 min up to 8 should significantly enhance the detectability of elusive species, especially in censuses that are not maintained over consecutive years (Leu et al. 2017). The trace used for playback stimulation was built using the sounds records by Schultze and Dingler (2007). Point counts targeted at warblers were censused with no acoustic stimulation and lasted for 5 min. Point counts for Water Rail were surveyed in the period 24 March–20 April, when the territorial activity of the species is at its peak in northern Italy (Brambilla and Rubolini 2004). Point counts focused on warblers were instead surveyed in the period 20 April–10 June, with two repetitions separated by at least 10 days.
Statistical analyses
We investigated the relationships between habitat variables (either collected at study site scale or at point count scale) and the probability of species occurrence by means of multivariate adaptive regression splines (MARS models). MARS is a non-parametric, machine-learning technique (Friedman 1991; Hastie et al. 2009) that well suits to investigate non-linear relationships between habitat features and occurrence (Elith and Leathwick 2007). Due to its effectiveness and flexibility, the use of MARS models is growing, and it has been recently applied in several studies based on survey data, analogous to those collected in the current study (e.g. Heinanen and von Numers 2009; Brambilla and Gobbi 2014; Brambilla and Pedrini 2016; Assandri et al. 2017). MARS models allow for multi-response tests, thus giving the opportunity to explore which factors have the main importance in determining species occurrence over an entire set of species, which are analysed at once. This MARS feature is particularly interesting for our set of data since some of our target species occurred at low frequency, and multi-response models have been shown to outperform single-species models for scarce taxa (Elith and Leathwick 2007; Brambilla and Gobbi 2014; Assandri et al. 2017). The core result of the final multi-response model is, therefore, the selection of a set of habitat variables that are likely important for all the species; at the same time, multi-response models individually describe the effect of each habitat variable on the occurrence probability of each species.
We fitted MARS models with the earth package, version 3.2-1 (http://cran.rproject. org/web/packages/earth/index.html), in R 3.1.2 (R Development Core Team 2018), using a binomial distribution of the response variable (Milborrow 2011a). We used the following settings for model selection: threshold = 0.01, penalty = 3. In all the cases, we tested models both allowing for the possibility of interaction between two variables (two-way interaction) and excluding the possibility of variable interaction (additive model; degree of interactions = 1, no interaction allowed among variables). We evaluated the variable importance by means of the evimp command (Milborrow 2011a; Brambilla et al. 2013; Jedlikowski et al. 2014). The evimp estimates variable importance in MARS models reporting (1) the number of models in which, after the pruning pass, is included a given variable; (2) the change in the residual sum of squares (RSS), scaled to 100, among each subset and the previous one thus indicating the change in the proportion of not-explained variance among the two models; (3) the generalized cross validation (GCV) of each model, which considers the change in the GCV once the reference variable is entered or removed from the model. GCV is calculated using the penalty argument and scaled to 100 (Milborrow 2011a). The plotmo package version 1.3-1 (http://cran.r-project.org/web/packages/plotmo/index.html) was used to plot the fitted functions (Milborrow 2011b). Results obtained with machine-learning methods are robust to collinearity in the set of the predictors (Garg and Tai 2013), but we applied a conservative approach and checked by means of VIF (variance inflation factor) that the set of the habitat variables presented a small degree of autocorrelation (VIF < 3; Salmerón Gómez et al. 2016). At this step, we found collinearity in both set of predictors and we therefore excluded the extent of meadows and those of transitional shrubland areas from further analyses at the site-level scale (max VIF value after exclusion: 2.01), and the extent of tree cover from the predictors used in the analyses at the point count scale (max VIF value after exclusion: 2.08). Based on these predictor sets, we run two multi-response MARS models, one to explore the habitat variables predicting the probability of occurrence of the target species at the site level, and a second one exploring the analogous relation at the point count level. For the point count scale, we also run single-species MARS models to compare their performance with those of multi-response ones. We did not run single-species models for data collected at the site level, since the sample size was considered to be too reduced to have reliable results, as some preliminary tests confirmed (details not shown).
Graphically, the output of MARS models is presented in form of plots representing with a spline the non-linear function linking the extension of a given habitat type and the occurrence probability of a given taxon. In most cases, these functions have a threshold point: above or below this value the spline flattens, thus indicating that further increases (or decreases) in habitat extension has no longer effects on occurrence probability. The exact values of this threshold point, if present, are indicated in the output generated with the summary function of the MARS models, and we took it as a threshold also for practical management recommendations.