Nestedness structure of bird assemblages in a fragmented forest in Central Argentina : the role of selective extinction and colonization processes

Nestedness structure of bird assemblages in a fragmented forest in Central Argentina: the role of selective extinction and colonization processes. Nestedness analysis constitutes an important tool to understand the processes that shape wildlife communities. It also allows a quick first evaluation of species extinction proneness in fragmented landscapes. Here, we tested whether avian assemblages in the fragmented Espinal forest exhibited nested subset patterns. Furthermore, we examined whether selective extinction or selective colonization are driving nested subset patterns. We studied avian assemblages in 13 forest fragments in central Argentina during breeding and non–breeding seasons. We completed partial Spearman rank correlations to explore the relationship between nestedness rank order and habitat patch variables and species life history traits related to species extinction proneness and colonization rate. Bird species showed strong nestedness patterns, both for the total incidence matrix and for forest fragments and species separately. Nestedness patterns were similar during the breeding and non–breeding seasons. The nested rank order of forest fragments correlated with area and distance to nearest fragment, both of which are patch characteristics known to increase the probabilities of species extinction. The nested rank order of species was correlated with the minimum area of species requirement, trophic guild, and range size, traits that are linked to extinction risk. Selective extinction processes rather than selective colonization appear to be driving nestedness patterns of bird assemblages in fragmented Espinal forest. The most effective way to preserve forest bird species in the Espinal forest seems to be by protecting the larger fragments of this relictual forest.


Introduction
Habitat loss and fragmentation are among the most important threats to biodiversity worldwide (Wilcove et al., 1998;Sala et al., 2000;Haddad et al., 2015). Broadscale habitat fragmentation gives rise to archipelagos of natural habitat fragments or islands immersed in a matrix of anthropogenic open habitat (Matthews et al., 2015). Since species sensitivity to habitat fragmentation in a particular region is variable, species loss in those remaining habitat islands does not necessarily occur at random but may occur in a nested pattern (Patterson and Atmar, 1986;Patterson, 1993, 1995). In nested assemblages, poorer communities constitute proper subsets of increasingly richer communities (Patterson and Atmar, 1986). Therefore, less widespread species occur on sites with relatively large species assemblages while poorer assemblages are mostly composed of ubiquitous species (Cutler, 1991;Soga and Koike, 2012). Consequently, in archipelagos with 'perfect' nestedness structure, it is possible to predict the order of disappearance of the less ubiquitous species from the poorer sites in response to environmental gradients (Atmar and Patterson, 1993) as the species that are present only in the richer fragments are more likely to become extinct as environmental disturbances increase (Nupp and Swihart, 2000).
Nestedness analysis is an important tool to understand the processes that shape communities and to reveal the ecological and evolutionary limits of the species. Furthermore, it has valuable implications for conservation (Wright et al., 1998;Martinez-Morales, 2005). Nestedness analysis is attractive because it allows a quick first evaluation of species extinction proneness in species assemblages of fragmented landscapes (Ganzhorn and Eisenbei, 2001). Although this approach alone is insufficient to evaluate strategies to preserve biodiversity in fragmented biotas (Cutler, 1994) it could be highly useful as predicting species loss can be used to make informed land-use decisions and to effectively protect species that will disappear first in a determined fragmentation scenario (Fleishman et al., 2007).
Four main processes have been proposed to explain nestedness patterns: (1) selective extinction of species with large spatial requirements in relation to fragment area (Wang et al., 2012;Matthews et al., 2015); (2) selective colonization of species with low dispersal ability in relation to fragments isolation (Kadmon, 1995); (3) random, passive sampling from a common species pool, which can result in a nested pattern if sites are more likely to be occupied by species that are regionally more abundant (Cook and Quinn, 1995;Wright et al., 1998); and (4) selective occupation of hierarchically nested habitats (Honnay et al., 1999). However, studies of the mechanisms explaining nestedness structure on archipelagos resulting from habitat fragmentation have found that, in most cases, nestedness structure is driven by selective extinction and, to a lesser extent, to selective colonization process (Watling and Donnelly, 2006;Matthews et al., 2015).
The selective extinction hypothesis is related to the concept of faunal 'relaxation' (Brown, 1978;Wilcox et al., 1986). It states that fragment area is the main driver of communities' structure as species loss is predictable and follows gradients of species sensitivity to habitat size (Wright et al., 1998;Watling and Donnelly, 2006;Matthews et al., 2015). Under this mechanism animal species with large area requirements, high trophic guild (i.e. carnivorous and insectivorous), small range size and large body mass will be the first to become extinct when the area of the fragments is reduced (Matthews et al., 2015;Keinath et al., 2017;Li et al., 2019). On the other hand, the selective colonization hypothesis states that the habitat isolation would be the main mechanism behind nestedness structure of an assemblage (Watling and Donnelly, 2006;Meyer and Kalko, 2008). Under this mechanism, species with low dispersal ability -such as those with a low dispersal ratio or small body mass-will colonize only the less isolated fragments and will fail to colonize those that are more isolated (Loo et al., 2002;McAbendroth et al., 2005;Frick et al., 2009).
Although explanation of nestedness structure under selective extinction and selective colonization implies the combination of site variables with species traits (Ulrich et al., 2009), few studies have tried to analyze their roles in generating nestedness simultaneously (Wang et al., 2012;Li et al., 2019).
The Espinal xerophytic forest in central Argentina provides a suitable scenario to address the effects of habitat fragmentation in species assemblage structure. Here, open forests historically used for cattle grazing have been converted to row crop production (Baldi and Paruelo, 2008). The expansion of cultivated land has been related to a combination of climate change (increasing precipitation), increasing global demand for agricultural products, national economic policies, and new technologies (genetically modified seeds, agrochemicals, machinery) (Grau and Aide, 2008;Zak et al., 2008).
At present, the Espinal xerophytic forest is an extremely degraded lowland forest with less than 5% of the original forest area Lewis et al., 2009;Morello et al., 2012;Noy-Meir et al., 2012). Because of this severe fragmentation and habitat loss, avian diversity has been negatively affected (Dardanelli and Nores, 2001;Dardanelli, 2006;Dardanelli et al., 2006). At least eight species appear to have become extinct in this forest in the province of Córdoba, Central Argentina, and another nine species are sensitive to fragmentation . However, fragmentation effects on the species composition and nestedness structure of avian assemblages have not been assessed. Studying drivers behind nestedness structure of avian assemblages in fragmented Espinal forest would provide insights that could help avian conservation. The design of effective management plans in poorly studied and highly fragmented habitats, such as the Espinal forests of Córdoba, Argentina, could take advantage of nestedness analyses, especially in a place where there is no time or resources to undertake long-term studies and when decisions for conservation action are urgent (Ganzhorn and Eisenbei, 2001;Fleishman et al., 2007). Here, we tested whether avian assemblages in the fragmented Espinal forest exhibit nestedness patterns for winter and summer assemblages. Furthermore, we examined the mechanisms underlying the nestedness structure, particularly focusing on whether selective extinction, selective colonization, or passive sampling are driving nestedness patterns for winter and summer assemblages.

Study area
Our study was conducted in the Espinal forest fragments in the eastern lowlands of Córdoba Province, Argentina ( fig. 1). Forest fragments are located in private lands since there are no protected areas in the region. The mean annual precipitation of about 700-800 mm falls mostly in late-spring and summer, from October to March; the rest of the year is the dry season. The mean annual temperature is 16 ºC, with a maximum peak of 44 ºC and minimum temperature of -9 ºC (Morello et al., 2012). This region is regarded a semiarid environment because of the high potential of evapotranspiration that generates a water deficit for 11 months of the year (Morello et al., 2012).
All fragments in our study had three well-developed vegetation strata (tree, shrub, and herbaceous), were completely isolated (no corridors or rivers connecting any fragment, fig 1), and were embedded in a matrix of croplands, mostly soybean during the austral summer and wheat or fallow fields during the winter. Thus, we considered the contrast between forest fragments and the matrix to be high (Lindenmayer and Fischer, 2006) for forest birds.

Bird sampling
We surveyed bird species in 13 forest fragments (ranged from 0.25 to 217.4 ha). Nocturnal species (Strigidae and Caprimulgidae) and species that only flew over the fragments were not considered. Surveys were conducted during the austral winter (June-August 2001) and austral summer (December2001-March 2002) seasons. One observer (SD) carried out all the surveys by intensive searches recording all bird species seen or heard while walking slowly through the whole fragment from pre-dawn to 11:00 and from 14:00 until sunset. We surveyed each fragment until no new species were added in 4-8 additional sampling days . To adjust for differences in species detectability we compared species richness among forest patches and between seasons using rarefaction curves. Rarefaction analysis calculates species richness after standardizing differences in abundance among samples by estimating the expected number of species of each sample if all samples are reduced to a standard size (Magurran, 2004). Rarefaction curves were performed using iNext (Chao et al., 2016).
We distinguished two types of birds occurring in the fragments: forest species (species that inhabit only xerophytic forests in the study area), and habitat generalists that use both forest and open areas (table 1s in supplementary material). Because the focus of this study is on patch level effects, we centred our investigation on species for which xerophytic forest is a primary habitat. Therefore, prior to analysis, we removed all species for which xerophytic forest is not considered primary habitat (Cook et al., 2002;Watson, 2003; table 1s in supplementary material). We also removed migratory species (Nores, 1996;Barnett and Pearman, 2001), considering only year-round residents as they necessarily colonize the fragments at the beginning of the breeding season and leave (disappear) at the end of the breeding season (Watson, 2003). However, we acknowledge the response of generalist and migratory bird species to fragmentation could have some relevance and would need to be considered when designing conservation measures at regional scales for the Espinal forests in Argentina. It is important to mention that most fragments and all species analyzed in this study persist in the study area (Verga et al., 2019;eBird, 2020). Thus, we consider that the results of our study could be applied to the current scenario, as both the fragmented forest and the bird species have remained constant.
The order of the families and the generic and specific names of bird species follow the South American Classification Committee (Remsen et al., 2020).

Species traits
To analyse the influence of extinction and the colonization process in structuring species occurrences, we selected species life history traits commonly associated with avian species extinction proneness or dispersal ability (table 2s in supplementary material). Geographic range size, trophic guild and natural abundance are life history traits related to extinction proneness (Davidar et al., 2002;Henle et al., 2004;Feeley et al., 2007;Wang et al., 2010). On the other hand, body mass and dispersal ratio area are life history traits that are usually linked to species ability to colonize new sites (Schoener and Schoener, 1984;Cook and Quinn, 1995;Henle et al., 2004;Jenkins et al., 2007).
We obtained distributional range size from Birdlife International species factsheets (BirdLife International, 2020). Trophic guild data were constructed by extracting data of local species diet (Zotta, 1940;Del Hoyo et al., 1992;Alessio et al., 2005;Salvador et al., 2017) and creating four categories: 1, herbivores; 2, omnivores; 3, insectivores, and 4, carnivores. Mini-  Dardanelli et al. (2006). Body mass data were obtained from Dunning (2008). The dispersal ratio was calculated by dividing each species mean wing length (mm) by the cube root of its mean mass in grams (Fischer and Lindenmayer, 2005;Li et al., 2019). The relationship of this ratio with dispersal ability is positive so that species with higher ratios will disperse longer distances and species with lower ratios will disperse shorter distances and will consequently be poor dispersers (Fischer and Lindenmayer, 2005;Li et al., 2019). Species traits were not correlated among them (Pearson r < 0.4).

Site variables
We selected different landscape variables to characterize spatial configuration of the forest fragments: area (Area; ha), perimeter (m), two isolation variables: distance to the nearest fragment (DNF; meters) and  B proximity index in a 2 km-buffer area (PI), and shape index (SI) (table 3s in supplementary material). Area, Perimeter and Isolation metrics were calculated using Quantum GIS (QGIS) software. We estimated shape index as SI = Pm/Pc, where Pm is the measured perimeter of the fragment and Pc is the perimeter of a circular fragment of the same area. This SI index was used in similar studies and has been found to be less correlated to the area than other shape indices (Hinsley et al., 1995;Santos et al., 2002;Watson et al., 2004). We found that the Perimeter and Shape index were highly correlated with Area (Pearson r ≥ 0.7). For this reason, we excluded these variables, as they were dependent on area.

Data analysis
Matrices of presence-absence were assembled for both seasons. We used the metric based on overlap and decreasing fill 'NODF' to evaluate nestedness (Almeida-Neto et al., 2008). Through the online interface NeD (http://ecosoft.alwaysdata.net/) developed by , nestedness can be calculated for the whole incidence matrix and independently for species (NODF between rows) and sites (NODF between columns). We ran five null models to estimate the significance level of nestedness: equiprobable row and column totals (EE), equiprobable row totals-Fixed column totals (EF), fixed row totals-equiprobable column totals (FE), and fixed-column and fixed-row totals (FF) algorithms. FF algorithm has shown to be highly restrictive and EE poorly restrictive (Ulrich and Gotelli, 2012;Matthews et al., 2015;Si et al., 2015). All these null models have strengths and weaknesses (Ulrich and Gotelli, 2012;Matthews et al., 2015;Si et al., 2015). However, PP and FF null models were described as less biased than the others (Ulrich and Gotelli, 2012). Furthermore, PP has been found to be the preferred model when research systems contain relatively small islands, when the scale of analysis is small, and because it is considered more ecologically meaningful (Ulrich and Gotelli, 2012;Matthews et al., 2015;Si et al., 2015). This model provides an unbiased proportional resampling of matrix incidences proportional to row and column totals (Almeida-Neto and Ulrich, 2011). Expected nestedness metrics and related parameters were generated for winter and summer assemblages by running 1,000 Monte Carlo simulations. The passive sampling hypothesis can be tested using the Coleman's (Coleman, 1981;Coleman et al., 1982) random placement model (Calme and Desrochers, 1999;González-Oreja et al., 2012;Wang et al., 2012;Li et al., 2019). The random placement model was used to verify whether passive sampling from species abundance distributions was driving the nestedness structure of bird communities. Coleman et al. (1982) state that the number of species ŝ (α) to be found residing in a given site depends on this site relative area, α (which equals the ratio of the area of a particular fragment to the summed area of all frag-ments), and the overall abundances n 1 , n 2 ,…, n s of the S species represented in C, which is a collection of N individuals from S species (Coleman, 1981): If the hypothesis of random placement holds roughly two-thirds of the points should fall within the band bounded by ± one standard deviation of the expected curve, or if less than two-thirds of the points fall within the bands, it should be rejected (Coleman et al., 1982).
To check for spatial autocorrelation in the data (i.e. figures of variables sampled at nearby locations tend to have more similar values than would be expected by chance) we fitted a semivariogram randomisation analysis based on 99 Monte Carlo permutations (Isaaks and Srivastava, 1989). Spatial autocorrelation in the response variable (species richness) violates the assumption of independently and identically distributed errors and hence inflates type I errors (Dormann, 2007).
The order in which sites and species are organized by NODF can be compared with several independent variables to evaluate their possible roles in producing nestedness (Patterson and Atmar 2000). To test the effects of forest fragment traits on nestedness, we performed Spearman rank correlations between the forest fragments rank orders in the maximally packed matrix and ranked traits of the forest fragments (table  3s in supplementary material). Similarly, to assess the role of species traits in driving nestedness patterns, we calculated Spearman rank correlations between the species rank orders in the maximally packed matrix and ranked species traits (table 2s in supplementary material). Statistical significance was established at P < 0.05. Partial Spearman rank correlations and semivariograms were performed using R (R Development Core Team, 2016).
The semivariogram did not show a significant association between the spatial distribution of the forest fragments and species richness in winter or summer ( fig. 1s in supplementary material).
The bird assemblages were significantly nested in both seasons (table 1) for all null models except the very restrictive FF (fixed-fixed) model. Both NODF values and matrix structure were similar between seasons for resident birds (table 1). Our results show a high temporal constancy in the nested pattern for resident bird assemblages in Espinal forest fragments in Central Argentina.
Spearman's rank correlations showed that the remnant order that maximized nestedness in both winter and summer was correlated with remnants ordered according to the area and distance to the nearest fragment (table 2). The proximity index was not significantly correlated with remnant order.
Species order in matrices packed for maximum nestedness showed a significant relationship with minimum area requirement, trophic guild, and range size both in winter and summer, and with species distributional range size in summer (table 2). Species traits commonly related to colonization ability like body mass and dispersal ratio were not significantly correlated with the species order.
The nestedness of forest birds' assemblages was not caused by passive sampling in summer or winter assemblages. Only one data point in summer and four out of 13 in winter data points fell within ± 1 SD of the expected Coleman's species/relative area curves ( fig. 3, 4), which means that it did not follow expectations from the random placement hypothesis.

Discussion
Bird assemblages in fragmented Espinal forest in Central Argentina showed a non-random structure, with species aggregation consistent with the nested subset model across seasons, for the whole matrix, and for columns (forest fragments) and rows (bird species) separately. This nested structure did not follow the random placement hypothesis (Coleman, 1981;Coleman et al., 1982).
The nestedness structure in our studied system showed a structure consistent with the selective extinction hypothesis as nestedness was related to fragment area and species traits associated with extinction proneness such as trophic guild, minimum area requirement and distributional range size. The correlation of fragment area and species traits with nested rank indicated that bird assemblages on  Table 1. Comparative analyses of nestedness for resident forest birds in forest fragments between seasons, in Córdoba, Argentina. Nestedness metrics and related parameters are provided for two seasons: winter and summer. P-values were generated by 1,000 Monte Carlo simulations: EE, equiprobable-eqiprobable null model; PP, proportional-proportional null model; FF, fixed-fixed null model; SD, standard deviation; * significant nestedness (P < 0.05); matrix, nestedness estimator for the whole presence-absence matrix; species, row nestedness estimator among species (based on species incidence); fragments, column nestedness estimator among fragments (based on species composition). Tabla 1. Análisis comparativos del anidamiento de las aves residentes en fragmentos de bosque entre estaciones, en Córdoba (Argentina). Se proporcionan los valores de anidamiento y los parámetros relacionados para las dos estaciones: invierno y verano. Los valores P se generaron a partir de 1.000 simulaciones de Monte Carlo: EE, modelo nulo equiprobable-equiprobable; PP, modelo nulo proporcional-proporcional; FF, modelo nulo fijo-fijo; SD, desviación estándar; * anidamiento significativo (P < 0,05); matrix, estimador de anidamiento para toda la matriz de presencia-ausencia; species, estimador de anidamiento de fila entre especies (basado en la incidencia de especies); fragments, estimador de anidamiento de columna entre fragmentos (basado en la composición de especies).

Fig. 3. Comparación de los datos observados y esperados según el modelo de ubicación aleatoria para aves forestales residentes en fragmentos de bosque en el verano austral, en Córdoba (Argentina). Las líneas continuas representan los valores esperados y las líneas discontinuas representan las desviaciones estándar asociadas (± 1 DE; líneas discontinuas). Los círculos representan la riqueza de especies observada.
These results are in agreement with similar studies in fragmented habitats where selective extinction arises as the most common driver of nestedness structure (Wright et al., 1998;Matthews et al., 2015;García-Quintas and Parada, 2017;De la Hera, 2019). The importance of extinction driven processes in shaping community assembly has been found in many fragmented landscapes (Martinez-Morales, 2005). It has been suggested that the trigger for this kind of patterns is a faunal relaxation process, which is characteristic of highly fragmented or relictual forest ecosystems (Brooks et al., 1999;Ferraz et al., 2007).
The prevalence of colonization driven patterns is less frequent in fragmented terrestrial habitats (Wright et al., 1998;Watling and Donnelly, 2006), and it appears to be an important driver for other isolated habitats such as mountaintops, land-bridge islands and oceanic islands (Cook and Quinn, 1995;Wright et al., 1998;Watling and Donnelly, 2006;Meyer and Kalko, 2008;García-Quintas and Parada, 2017). Our results showed that selective colonization seems to have some influence on community assembly as distance to nearest fragment was correlated with fragments nested order. However, species traits commonly associated with dispersal ability such as body mass and dispersal ratio were not related to species order. In this regard, it is possible that dispersal ratio and body mass were not good indicators of dispersal ability for birds in this study. The other possibility is that because colonization has marginal importance in driving nestedness structure in our system, it does not express any significant relationship with our dispersal ability proxies. Consequently, we could venture to say that the selective colonization hypothesis only partially explains birds' nestedness structure in Espinal forest of Central Argentina. One possible explanation for the low influence of selective colonization is that species with low dispersal ability have already become extinct in the study area. We documented this in a previous study in the same area where, for example, most large birds have disappeared from fragmented forests . Moreover, as has been demonstrated in other studies (Watling and Donnelly, 2006;Matthews et al., 2015), it is very difficult to find biologically meaningful isolation effects on assemblage structure. In the case of bird communities in South America, it is even more challenging considering knowledge of colonization ability or dispersal rate of species is scarce (Faaborq et al., 2010;Jahn et al., 2017). It is therefore challenging to assess the role of selective colonization hypothesis in the assemblage structure of Espinal forest birds. Nevethless, it seems to have secondary importance as a driver of nestedness aggregation. This highly fragmented forest has almost disappeared from this region and the few remaining fragments have undergone faunal relaxation for many years, giving rise to extinction driven biotas. Consequently, we did not find any relationship between proximity indexes or any species trait related to dispersal ability with nestedness order.
Nestedness structure did not vary between seasons. Consequently, seasonality does not seem to influence nestedness in our system. These results are consistent with the results of Seoane et al. (2013), García-Quintas and Parada (2014), Zhou et al. (2014) and De la Hera (2019) who found no seasonality effects on nestedness structure of birds in isolated woodlots in Spain and birds of urban parks in Hong Kong and Spain. Our results, however, contradict the results of Murgui (2010) who found small but significant differences in nestedness structure and species-area relationships between seasons in urban parks in Spain. This author ruled out an increase in mortality outside the breeding seasons as they have mild winters. He considered that the use of alternative habitats outside of parks during autumn and winter is the most likely explanation for the observed patterns. The winters in the Espinal forest fragments in Central Argentina are mild and bird species are probably less prone to use alternative habitats than birds in urban parks. Another difference is that specialist birds analyzed in our system are generally more sensitive to disturbances and less adaptable than generalist species in urban parks in Spain.
Protecting the larger and less isolated forest fragments would be the most effective way to preserve resident birds in this relictual Espinal forests. The preservation of large and less isolated fragments would help to protect resident birds with large area requirements, small distribution range size, and high trophic guilds (i.e. carnivorous and insectivorous species). For example, by protecting the two largest forest fragments it is possible to maintain most species (97.7 % in summer; 95.3 % in winter) of forest birds in the dataset. As mentioned by other authors, nestedness analysis can be used in combination with other approaches to provide valuable recommendations for decision-making when long-term data are not available. Based on the results of the present study, future landscape management of Espinal forest should ensure the protection of large fragments as they preserve the largest populations of resident forest species throughout the year.

Acknowledgements
We thank G. and M. Esmóris, R. Parra, A. Varselotti and F. Mansilla for providing access to their properties. We also thank D. A. Serra and M. Nores for assistance with the bird surveys. Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) provided partial funding.