Integrating Landscape Disturbance and Indicator Species in Conservation Studies

Successful conservation plans are conditioned by our ability to detect anthropogenic change in space and time and various statistical analyses have been developed to handle this critical issue. The main objective of this paper is to illustrate a new approach for spatial analysis in conservation biology. Here, we propose a two-step protocol. First, we introduce a new disturbance metric which provides a continuous measure of disturbance for any focal communities on the basis of the surrounding landscape matrix. Second, we use this new gradient to estimate species and community disturbance thresholds by implementing a recently developed method called Threshold Indicator Taxa ANalysis (TITAN). TITAN detects changes in species distributions along environmental gradients using indicators species analysis and assesses synchrony among species change points as evidence for community thresholds. We demonstrate our method with soil arthropod assemblages along a disturbance gradient in Terceira Island (Azores, Portugal). We show that our new disturbance metric realistically reflects disturbance patterns, especially in buffer zones (ecotones) between land use categories. By estimating species disturbance thresholds with TITAN along the disturbance gradient in Terceira, we show that species significantly associated with low disturbance differ from those associated with high disturbance in their biogeographical origin (endemics, non-endemic natives and exotics) and taxonomy (order). Finally, we suggest that mapping the disturbance community thresholds may reveal areas of primary interest for conservation, since these may host indigenous species sensitive to high disturbance levels. This new framework may be useful when: (1) both local and regional processes are to be reflected on single disturbance measures; (2) these are better quantified in a continuous gradient; (3) mapping disturbance of large regions using fine scales is necessary; (4) indicator species for disturbance are searched for and; (5) community thresholds are useful to understand the global dynamics of habitats.


Introduction
Habitat destruction converts previously homogeneous landscapes into a mosaic of new land use types that can be unfavorable to species persistence [1][2][3]. This mosaic usually includes areas of intensive agriculture, exotic species plantations, non-managed alien vegetation and remnant natural areas. However, border areas between different habitats create edge effects, often promoting the extinction of many disturbance sensitive species in native habitats [4]. Such transition areas are called ''ecotones'' and often represent diverse habitats in the ecological boundary between native vegetation and managed habitats [5][6][7]. In addition, ecotones may promote the dispersal of species through landscapes (e.g. invasive species), modify species interactions and trophic web structure and promote ecosystem change [8].
Often, only a few habitat patches and/or sites are kept relatively safe from human disturbance and it is becoming increasingly important to identify which areas still keep mostly intact communities. Detecting ecological community thresholds, i.e., points along continuous gradients where there are major changes in community composition [9][10][11], is a way to identify where and how to protect native species and their biotic interactions [12]. This is particularly important when clear anthropogenic gradients are identified, because such gradients may be explicitly modified in a conservation planning framework, for example, by selecting reserve areas where certain uses are forbidden and others encouraged. It is important to note that two different approaches can be followed when studying edge effects [13]: (1) a discrete approach where 'edge' is compared with 'interior'; (2) assessing how communities change along a continuum.
A number of methods have been proposed to identify ecological thresholds [13][14][15][16]. Most are aggregated methods, in the sense that species-specific responses are often lost, especially for the least abundant species, which often are of special conservation concern [17]. Recently, Baker and King [18,19] introduced a new method called Threshold Indicator Taxa ANalysis (TITAN) to identify ecological thresholds at both individual species and community level along environmental gradients. This approach revealed to be useful for different taxa and environmental gradients at different regions. However, to our knowledge, it was never tested along an anthropogenic disturbance gradient that explicitly considers a mosaic landscape configuration. In this paper, we present a new approach to identify ecological community thresholds along disturbance gradients, combining TITAN with a new disturbance metric that explicitly considers landscape configuration. To test our approach, we studied Azorean arthropods, which represent a particularly interesting model system for both theoretical and practical reasons.
Oceanic islands are home to large numbers of endemic species, either through in-situ speciation (neoendemics) or after the extinction of populations outside the islands (palaeoendemics). Such species, together with other native species (which occur anywhere else and reached the islands through their own means) create very unique communities [20]. These communities are often more prone to be disrupted by the introduction of exotic species than their continental counterparts, as they may not be able to cope with the introduction of certain, previously inexistent, ecological traits into the ecosystem [21][22][23]. Additionally, being isolated, many populations are unable to recover from past disturbance events, when rescue effects are impossible. As such, oceanic islands have been stage to most recorded extinctions worldwide, mainly driven by habitat destruction and species introductions [20,24,25].
The Azorean archipelago, which was mostly covered by Laurisilva forest prior to human settlement, has undergone drastic land use changes since the first inhabitants arrived almost 600 years ago [26][27][28]. Such changes are thought to have caused the extinction of numerous endemic species, particularly in the most disturbed islands, where few and minute native forest patches remain [29]. However, while some human-modified habitats, such as exotic forest plantations, may harbor some endemic and native species [26,30], others (like intensively managed pastures) are known to be mainly dominated by exotic species [31,32]. In particular, there is a clear dissimilarity gradient in community composition according to the disturbance level of a particular land use type [26,33]. If and where any threshold in community composition along a disturbance gradient can be found is nevertheless still unanswered.
Arthropods are the most diverse group in the Azores as mostly elsewhere [34]. Although they provide multiple ecosystem services, they are often neglected in conservation policies and programs, even if known to be prone to the same threats as any other organisms [35]. This is caused by different reasons, the most obvious being the lack of information on which species are living (the Linnean shortfall), where they are persisting (the Wallacean shortfall), how their abundances change in space and time (the Prestonian shortfall) and how sensitive they are to ecological change (the Hutchinsonian shortfall) [35]. Because of intensive sampling during the latter decades [36], the Azorean Islands are exceptions to the rule and much is now known about their arthropod fauna, including how different species abundances change in response to anthropogenic disturbance.
As seen, the arthropods of the Azorean Islands are a particularly well suited case study to illustrate the general aims of our approach, namely, (1) to quantify a continuous disturbance gradient according to the different land uses, their perceived disturbance level and the surrounding matrix of different uses; (2) to find the species that respond to the gradient according to their abundance and spatial distribution; (3) to determine the shift point for each species responding significantly to the gradient (i.e. uneven distribution along the gradient); (4) to assign species to either negative or positive response groups (i.e. the ones whose abundances are negatively or positively correlated with disturbance); (5) to combine species information to compute community threshold(s) and; (6) to investigate whether some particular features of the selected species may be related to the direction of their response.

Study area
The Azorean archipelago is located in the North Atlantic, roughly between 37u to 40uN latitude and 25u to 31uW longitude. It comprises nine main islands and some small islets aligned along a roughly WNW-ESE trend. They are fully oceanic, that is, they are totally volcanic islands of recent origin (8.12-0.25 Myr for the oldest areas of the main islands). In this study we focused on Terceira Island (402 km 2 ), as this island had the most information about all arthropod taxa and land use types [32,36]. Terceira is formed by four main volcanic polygenetic complexes (Cinco Picos, Guilherme Moniz, Pico Alto and Serra de Santa Bárbara). The highest point (Serra de S. Bárbara, 1023 m) is also the most recent (0.025 Myr B.P.) of the three major island complexes. The climate is temperate oceanic, i.e. strongly influenced by the ocean and by its topography, which produces high levels of relative atmospheric humidity that can reach 95% of annual average in the native high altitude semi-tropical evergreen laurel forest, while it restricts temperature fluctuations throughout the year.

Fieldwork
We used a dataset of arthropod species distribution and abundance in Terceira Island previously published [26]. A total of 72 sites were sampled between 2000 and 2007, between June and September, once per site: 36 in natural forests and 36 in nonnative habitats, the latter including 9 in exotic forests, 11 in seminatural pastures and 16 in intensively managed pastures. We tried to spread the sampling sites all over the island independently of the surrounding land use matrix, although intensively managed pastures tend to be located in peripheral low-altitude areas, while natural forests tend to be present only in central high-altitude areas. At each site, epigean soil fauna was captured along 150 m long transects. Thirty pitfall traps, consisting of plastic cups with a top diameter of 42 mm and 78 mm deep, were dug into the soil so that the rim of the cup was level with the surface. Half of the traps were filled with approximately 60 ml of diluted ethylene glycol (anti-freeze liquid) and the other half with the same volume of a general attractive solution (Turquin). Traps were spaced 5 m from one another, alternating Turquin and ethylene glycol traps, and were left open for two-weeks at each site [26]. All Araneae, Opiliones, Pseudoscorpiones, Diplopoda, Chilopoda and Hexapoda (excluding Collembola, Diplura, Diptera and Hymenoptera) were initially sorted into morphospecies by students using vouchered specimens under supervision of a trained taxonomist (PAVB). All unknown morphospecies were subsequently sent to several taxonomists for species identification. All species were classified as endemic (E), native non-endemic (N) or introduced/ exotic (I). Endemic species refer to species found only in the Azores. Native non-endemic species arrived by long-distance dispersal to the Azores, cannot be associated with human activities (intentional or accidental human introduction) and are also known from other regions. Exotic species are those believed to have arrived in the archipelago as a consequence of human activities and often have a cosmopolitan distribution. All species were also classified as predators, herbivores, fungivores or saprophagous, according to the criteria of Moran and Southwood [37]. Fungivores are a specialized guild of beetles living in forest litter or associated with dead wood or under bark, many of the species found in Azores being endemic (e.g. Zopheridae), and although very few species were present, we kept them separate from herbivores.

Disturbance index
We designed a new index of ''landscape disturbance'' (D) that reflects an anthropogenic disturbance gradient by explicitly considering landscape configuration. A land use map of 1006100 m resolution depicting the location of all land use types was built based on aerial photography [38], with native forests further delimited and confirmed by fieldwork (C. Gaspar, unpublished data). Based on previous work [26], we knew the proportion of endemic, native and exotic species typical to each land use type present in the island. With such data, it was possible to infer the disturbance level of each land use relative to an undisturbed native forest. This was used to rank the different land uses and each was given a value of ''local disturbance'' (L) as follows : Natural forests = 0, Natural(ized) vegetation or rocky outcrops = 1, Exotic forests = 2, Semi-natural pastures = 3; Intensively managed pastures = 4; Orchards/agriculture areas = 5; Urban/industrial areas = 6. Different scales for L values were tested with similar results and we chose to present and discuss the simplest case. To the ocean we attributed the value of ''no data''. The landscape disturbance index of each 1006100 m cell in the island was then calculated as: where: D i,j is the final index value of the cell in row i and column j; L is the local disturbance value of each cell (as defined above); r is number of rows in the map; c is number of columns in the map; d is the distance between the centroids of each two cells; max is the maximum theoretical value of disturbance each cell may take (in this case max = 6, corresponding to urban/industrial areas). Thus, the influence of each cell surrounding the focal cell is inversely proportional to the square of the distance between the two cells.
That is, a cell next to the focal cell (d = 1) has 4 times more influence that a cell two rows apart (d = 2). Although all cells in the island were taken into account, the ones far away from the focal cell had an almost negligible individual influence. The L value of the focal cell was multiplied by 2 to guarantee that the land use of the focal cell (d = 0) is double weighted compared with the immediately surrounding cells (d = 1). Also, the division by the maximum value of each cell was necessary to guarantee that the presence of a land use without L values at the border of the study area (in our case the ocean) would not make all border (coastal) cells having low D (as if the ocean was equivalent to natural forest (L = 0)). As a beneficial side-effect, this division by the theoretical maximum also guarantees that D is not open-scaled, since D is limited between 0 (no disturbance at all, only possible if all cells in the island had natural forest) and 100 (maximum possible disturbance, only possible if all cells in the island were occupied by urban/industrial areas).

Identifying species response to disturbance
To identify species response to disturbance, we use Threshold Indicator Taxa ANalysis (TITAN). TITAN is a recently developed method [18] allowing the identification of threshold(s) or change point(s) along environmental gradients for each taxon. Basically, TITAN detects changes in species distributions along an environmental gradient over space or time, and assesses synchrony among species change points as evidence for community thresholds. TITAN uses IndVal (Indicator Value) scores [39] to identify change points along a continuous environmental gradient. Indval is a method based on the comparison of relative abundances and relative frequencies of species in different groups of sites chosen a priori. A good indicator species is one which occurs at all sites in a given group and never in any other groups [39]. Here, midpoints between observed values of disturbance are candidate change points (x i ) used to iteratively split observations into two groups, and thus produce two IndVal scores at each split. The relative magnitude of IndVal scores for groups on each side of a candidate change point reflects whether a species shows greater association with the left (negative response with respect to x = associated with low disturbance) or the right (positive response = associated with high disturbance) side of each split. The greater the difference in association created by a specific split, the greater the IndVal score for one of the two groups (i.e. low or high disturbance). The greatest IndVal score at each split and the direction of the split (i.e. toward low disturbance or toward high disturbance) on which it occurs are kept for comparison with those at other candidate change points. The probability of obtaining an equal or larger IndVal score from random data is estimated by comparing the magnitude of each observed IndVal score with those generated by randomly assigning group membership via permutations [39]. TITAN uses bootstraps to compute the confidence interval of the change point location along the gradient for each taxon. Additionally, this bootstrap procedure is also used to estimate two new measures for each taxon: the ''purity'' (proportion of the bootstrap replicates with the same response direction, i.e., negative or positive), and ''reliability'' (proportion of the replicates with p-values for the indicator value score at the change point below a specified probability). Therefore, a species may be considered significantly associated to either low or high disturbance if IndVal ,0.05, purity .0.95 and reliability .0.95.

Identifying community disturbance threshold
By using the disturbance thresholds identified for each species (see section above), subsequently TITAN is also able to identify two community thresholds both associated with low and high disturbance. To estimate these two thresholds, the IndVal scores are first standardized to z scores (i.e. by subtracting the mean of the randomization from the observed IndVal and dividing by the SD of the randomization) to describe the magnitude of the response relative to each species' abundance distribution. Standardization is required to allow rare species with small indicator value scores to have high z scores if they show a sharp response at a particular point along the disturbance gradient. Then, species are grouped according to the direction of the response into z2 species, with a negative response (low disturbance affinity), and z+ species, with a positive response (high disturbance affinity), along the disturbance gradient. Finally, the z2 and z+ values at each point along the gradient are summed. Then, negative and positive community thresholds correspond to the disturbance value where the sum (z2) or sum (z+) scores show a peak, respectively. These two community thresholds can be subsequently mapped to visually identify which sectors of the study areas (in our case, the study island) they delimited.

Investigating relationship between species response and species attributes
The TITAN analysis distinguishes two species groups: (1) species significantly associated with low disturbance (i.e. negative z2 response; significant IndVal scores ,0.05, purity .0.95 and reliability .0.95), and (2) species significantly associated to high disturbance (i.e. positive z+ response; significant IndVal scores ,0.05, purity .0.95 and reliability .0.95). We therefore used TITAN results to test for potential association between species attributes and species response to disturbance. We considered four main characteristics: biogeographic categories (endemics, natives non-endemics and exotics), taxonomic order, feeding guild (predator, herbivore, fungivore and saprophagous) and body size. Since our data were greatly unbalanced (i.e. sometimes one species per characteristics' level), we tested each factor independently by using a series of x 2 tests and non-parametric tests. First, we tested whether biogeographic categories, taxonomic orders, feeding guilds or body sizes differed between the two groups of species (i.e. between species with negative response and species with positive response). x 2 tests were performed to assess differences in species frequency of the first three characteristics, while differences in body size were assessed by performing Kruskal-Wallis tests. Second, we assessed whether taxonomic order, feeding guild or body size differed among biogeographic categories within each group of species by using x 2 tests and Kruskal-Wallis tests. To account for potential expected frequencies below 5, all p-values for the x 2 tests were computed by permutations tests using Monte Carlo simulations. For Kruskal-Wallis, if significant differences were detected, pairwise Mann Whitney-U post hoc tests were implemented.

Statistical analysis software
All the analyses were performed using R [40]. Particularly, we used the R code provided by Baker & King [18] to implement TITAN in R 2.13.0. Abundance data were transformed (log 10 (x+1)) to reduce the influence of highly abundant species on IndVal scores [18]. Following the recommendations of Baker & King [18], taxa occurring in less than 3 sites across the gradient were deleted to remove outliers representing a potential bias. In our study, we also use a minimum group size of five observations to compute IndVal for TITAN analyses. We reran TITAN across 500 bootstrap replicates to compute purity and reliability of individual threshold indicator taxa and uncertainty surrounding thresholds based on the distribution of maximum TITAN sum(z2) (individual taxa) and TITAN sum (z+) values.

Results
The landscape disturbance index was calculated for all the 1006100 m cells in Terceira Island and mapped (Fig. 1). Although it closely follows the distribution of different land uses, it allowed obtaining a continuous landscape of disturbance values for the entire island. For the 72 sites sampled across the different land uses, the index ranged from D = 14.38 (the lowest disturbed site) to D = 75.58 (the highest disturbed site).
The TITAN framework identified 72 species out of the total 140 analysed (51%) with significant indicator value for disturbance (Supporting Information S1). Twenty-eight species composed the group significantly associated with low disturbed sites (i.e. corresponding to 38% of the 72 species identified). The disturbance threshold of these species ranged from 17.17 to 60.17 but 24 of the 28 species had a threshold below 40 (Fig. 2). Forty-four species composed the group associated with high disturbance (i.e. corresponding to 62% of the 72 species identified) with a threshold ranging from 23.94 to 70.80.
Additionally, the TITAN analysis allowed estimating two community-level-thresholds associated respectively with low and high disturbance by identifying peaks in distribution of the sum (z2) and sum (z+) along the gradient ( Table 1). The distribution of both sum (z2) and sum (z+) did not show abrupt peaks meaning that the existence of clear disturbance thresholds is not evident. Therefore, these thresholds should be interpreted with caution. This uncertainty is also reflected in the large confidence limits associated with them ( Table 1). The threshold associated with low disturbance roughly followed the distribution of native forests (Fig. 3a), although other habitats connecting some of the forest fragments were also included and many forest margins were excluded. The threshold associated with high disturbance roughly followed the distribution of intensively managed pastures, orchard/agriculture areas and urban/industrial areas (Fig. 3b), although marginal areas were either included or excluded depending on the surrounding landscape.
Investigation of the link between species response groups and species characteristics revealed significant differences in biogeographic categories and taxonomic order composition (Table 2; Fig. 4). As expected, most of the endemic species are significantly associated with low disturbance areas and most of the exotic species are significantly associated with high disturbance areas. However, within the group of species associated with low disturbance, the three biogeographic categories are almost equally represented with 11 endemic, 9 native and interestingly, 8 exotic species. Exotics clearly dominate the group associated with high  disturbance areas with 32 exotic, 10 native and only two endemic species. Species associated with low disturbance are taxonomically diverse, but Araneae (n = 8), Coleoptera (n = 5) and Hemiptera (n = 5) are the most species rich. While the Coleoptera are almost all endemics, there are a few native and even exotic spiders and Hemiptera associated with low disturbed sites. Among the species associated with high disturbance, Coleoptera (n = 26) and, in a lesser extent, Araneae (n = 10), dominate. All spiders and most beetles are exotics. Feeding guild and body size did not show any significant association with species response groups. Within the group associated with low disturbance, no significant differences were found between endemics, natives and exotics in their taxonomic order composition, feeding guild or body size. Within the group associated with high disturbance, only body size was slightly different between the three biogeographic categories. However, the Mann Whitney-U post hoc tests did not show significant pairwise comparisons (All p.0.05).

Discussion
The general goal of this study was to illustrate a new approach for spatial analysis in conservation biology. This was based on a two-step protocol: (1) the creation of a new metric that provides a continuous measure of disturbance for any focal community on the basis of the surrounding landscape matrix (disturbance index D) and; (2) the use of a recently developed analysis of taxa responses to specific gradients to estimate species and communitylevel disturbance thresholds (TITAN analysis). Spatial detection of anthropogenic disturbance effects on species communities is a cornerstone concern in conservation biology and we believe that the approach proposed here may be useful for future researches.

Disturbance index
Ecologists often consider disturbance levels as being well represented by habitat categories (e.g. native vs. anthropogenic habitats [41,42]). Although this approach is useful in many contexts and often the only possible option due to the lack of data, disturbance is more likely to be a continuous phenomenon than a discrete one [43]. A number of metrics have been previously developed to measure disturbance on a continuous scale. These may be based on multiple factors that are thought to disrupt community composition and may be intrinsic to the focal sites, such as urbanization, logging, trampling, livestock grazing and land degradation or dependent on the surrounding landscape matrix, such as fragmentation or connectivity between patches [44][45][46][47]. All these may be studied either in isolation or through aggregated indices. Taking the surrounding landscape into account is, however, critical, given the meta-population and source-sink dynamics that necessarily occur among different habitats [32,48], which may determine that habitat fragments with similar local disturbance levels may present different species compositions depending on the surrounding matrix [45,49]. In a previous work on Terceira Island, Borges et al. [32,48] found that a source-sink dynamics between different habitats was, in fact, critical to understand community composition, particularly in marginal areas of forest fragments. Our newly developed disturbance measure reflects the effects of both local and regional processes in a single, continuous metric. Consequently, it can be viewed as a simple univariate gradient, easy to integrate in many statistical frameworks and with a wide application in conservation biology studies.
Transforming a simple ordinal scale of disturbance intensity into a continuous measure depicts disturbance patterns in a more realistic way, especially in buffer zones between different land uses. In our application of the new index D, L-values of local disturbance were arbitrarily established by the authors according to the knowledge of the studied system and the aim of the study. In general, L values can be determined a priori to reflect environmental disturbance, independently of the characteristics of the group of concern as done here. For this, the researcher identifies a number of basic land use types or habitats that reflect different levels of anthropogenic impact and assigns a score to each category [50]. For example, completely artificial habitats, such as urban or industrial areas, are assumed to be the most disturbed, thus receiving the highest score, whereas virtually pristine habitats, such as primary forests, receive the lowest score. What constitutes Figure 3. Maps of Terceira Island with (a) areas below the community threshold for the negative response (sum z2) and (b) areas above the community thresholds for the positive response (sum z+). (a) Areas below the community thresholds for the negative responses (sum z2) of species to disturbance gradient are native forests (dark green) and other habitats (yellow). Areas of native forest that are above this threshold are also represented (light green). (b) Area above the community thresholds for the positive (sum z+) responses of species to disturbance gradient are the intensively managed pastures, orchards/agriculture areas or urban/industrial areas (dark red) and other habitats (yellow). Intensively managed pastures, orchards/agriculture areas or urban/industrial areas that are below this threshold are also represented (light red). doi:10.1371/journal.pone.0063294.g003 Table 1. Community-thresholds and its associated percentiles. the extremes of the scale will depend on the particular context. For example, in certain ecosystems, such as most areas in Europe, where even well-preserved biotopes suffered some kind of disturbance starting in prehistoric times, the currently 'most natural' settings are always the result of a certain disturbance, thus the lowest score will be given to habitats that are in fact not pristine. Regarding the scale to adopt, the most obvious solution is a linear one, and this can be particularly appropriate if transitions among habitats are assumed to be sufficiently smooth, as in our case study. However, other scales, such as a geometric series, might be appropriate if one would want the values to reflect strongest differences between natural or semi-natural habitats vs. highly impacted or completely artificial habitats. When L values are intended to express 'absolute' disturbance, the index can be used to identify which species respond negatively, and which respond positively, to increasing disturbance, and this is the approach illustrated in the case study discussed here. We can also imagine a reverse situation, in which the researcher uses a particular indicator group to evaluate the disturbance level of the areas or habitat types. In this case, the response of the indicator group to disturbance can be used to assign the scores to the habitat/land use types. For example, species known to increase their population densities with increasing disturbance (as commonly found for invasive plants or opportunistic animals, such as many birds associated with urban areas) can be used to assign scores to habitats according to their frequencies or abundances in each habitat.
Besides integrating both local and regional processes and being a continuous measure, the index we propose has the additional advantage that it is particularly adequate for mapping disturbance of large regions, as it does not require measuring a number of variables at all single sites/cells/habitat patches in a region of interest. This allows, among other applications, the use of this metric to model predicted species distributions [51] accounting for the disturbance value of each site/cell. In fact, we used our disturbance index here presented to model 47 arthropod species in Terceira Island and it had a high explanatory power for the distribution of most species, being one of the main contributing variables to most final models [52].

Indicator species analysis
As explanatory variables, such as diverse disturbance sources, may be studied in isolation or in combined metrics, response variables, usually the presence/absence or the abundance of taxa, may also be studied either in isolation or as a combined community metric. Using the newly developed TITAN protocol [18] allows both options under the same framework. This analysis allowed us to simultaneously look for individual indicator species and analyse their characteristics and discover community thresholds common to most of the species.
As previously suggested [26], local composition of arthropod communities in Terceira is mostly determined by the presence of nearby natural forests or intensively managed pastures. To the latter land use we may add agricultural and urban/industrial areas. Native forests are the main source for endemic species, whereas intensively managed pastures and even more highly disturbed land uses are the main source for exotic species [26,32]. Therefore, their borders are found by the TITAN analysis as community thresholds for species responding negatively and positively to disturbance respectively. Native non-endemic species as a group do not show clear patterns, as although they are naturally occurring in the islands, they tend to be less specialized in habitat type. In fact, the Azorean Laurissilva is a very particular forest, completely different from, for example, the Madeira Laurissilva, which has both a different plant species composition and a radically different structure [53]. Forest trees in Madeira present a high stature and the ground is covered by leaf litter, whereas the trees in the Azores are contrastingly low in stature and the ground is mostly covered by mosses and ferns. Native species, which occur in other regions besides the Azores, must therefore be able to cope with different habitats, and this is evident from our results. The communities present in semi-natural pastures and exotic forests are a mixture of the other land uses, although they seem to play an important role as corridors between natural forests for both endemic and native species or even as a refuge for some endemic species [26,30].
Analyzing the two orders with most indicator species, spiders and beetles, the differences are particularly evident. Spiders  biogeographic categories (a-b), taxonomic order (c-d), feeding guild (e-f) and body size (g-h) for the species group associated with low and high disturbance (n = 28 and n = 44 respectively). The three biogeographic categories (endemics, natives and exotics) are also distinguished in the panels c to h for correspondence with the contingency table and body size distributions used in within groups tests (See table 1). Note that for the species group associated with low disturbance, 43 species out of the 44 were considered since feeding guild was undetermined for one species (see Supporting Information S1). Boxplots indicate the distribution of the body size per biogeographic category. For graphical convenience, body size was log-transformed. doi:10.1371/journal.pone.0063294.g004 identified as indicators of low disturbance sites are mostly endemic or native although two exotic species, Agyneta decora and Ero furcata mostly occur in native forests. It is hard to know the impact of these species on the native and especially the endemic assemblages because closely related and therefore competing species may have been driven to extinction already, either in Terceira Island only or in the archipelago [29]. The first species has a close endemic relative in the Azores, Agyneta rugosa. Interestingly, although both species occupy different islands from West to East of the archipelago, they only co-occur on São Jorge, an island that, while presenting few and small native forest patches, has some of the most undisturbed forests in the archipelago, namely Topo, which is a small fragment with two single island endemic spider species [44,52,54]. Ero furcata is a specialist spider-hunter, deceiving web-builders in their own snares as most members of the family Mimetidae [55]. As no other araneophagic spiders live in the archipelago except for three exotic Ero, these species may have just occupied an empty niche. Their impact on their spider prey is however unknown. All spiders species identified as indicators of disturbed sites are exotic species, reflecting the lack of particular adaptation of native and endemic species to such habitats. Beetles indicators of low disturbance sites are mostly endemic, and include species known previously to be dependent of high altitude pristine native forests [48,56]. Many endemic species only occur in native core areas, distant to the nearest ecotone, confirming previous findings [48]. On the contrary, generalist species experiencing landscape source-sink dynamics can cope with local disturbance in a dynamic fragmented landscape [57]. One endemic beetle, Heteroderes azoricus, is indicator of highly disturbed habitats, a fact also demonstrated in previous studies [30].
It should be noted, however, that the TITAN approach uses a step-function model of community threshold, one which tries to identify steep changes in the frequency or abundance of taxa along environmental gradients. In many cases the change may be gradual, such as in broken-stick or dose-response models [58,59]. In such cases great care should be taken in the interpretation of results. In our study, the community thresholds were correctly identified, as the disturbance index is based on the spatial distribution of discrete land use types. Moreover, results corroborate previous studies [26] and match our empirical knowledge of the ecology of the studied communities and individual species.

Conclusions
Concluding, the landscape disturbance index coupled with the TITAN analysis allows perceiving patterns that would probably go unnoticed otherwise. This framework may be useful in many different situations, namely: (1) when both local and regional processes are to be reflected on single disturbance measures; (2) when these are better quantified in a continuous gradient; (3) when mapping disturbance of large regions using fine scales is necessary; (4) when indicator species for disturbance are searched for and; (5) when community thresholds are useful to understand the global dynamics of habitats.

Supporting Information
Supporting Information S1 List of the 72 species having significant response with either low or high disturbance and their corresponding threshold (with associated 90% confidence limits) and the z-score (See material and methods). Taxonomic information as well as biogeographic categories, feeding guild and body size are also given. (PDF)