Increase of insular exotic arthropod diversity is a fundamental dimension of the current biodiversity crisis

A dramatic insect decline has been documented on the grasslands and forests of European or North American mainland. Yet, other parts of the world and other ecosystems remain much less studied with unknown patterns. Using a unique time‐series dataset, we investigate recent trends on abundance and richness of arthropods sampled in Azorean native forest over 6 years (2013–2018). We test the hypothesis that biodiversity erosion drivers are changing the diversity and relative species abundance structure (species abundance distribution, SAD) of endemics, native non‐endemics and exotic species over time. We also examine temporal trends in abundance for each individual species. In contrast with mainland studies, we observed no decline in overall arthropod diversity, but a clear increase in the diversity of exotic arthropods and some evidence of a tendency for decreasing abundance for some endemic species. We also document stronger species turnover for exotic species, but no specific changes in the SAD. We argue that many changes, particularly in unique systems such as islands, will be noticed not at the richness but mostly at compositional level. Special attention should be given to exotic species which are known to be one of the major drivers of biodiversity erosion on islands.


Introduction
Biodiversity loss is well recognised as a key challenge for this century Cardinale et al., 2012). Insect biodiversity in particular is an essential component for ecosystem functioning (Allan et al., 2015;Bennett et al., 2015), so loss in insect diversity is a risk to both ecosystem sustainability and human health and well-being (Sandifer et al., 2015). Despite their importance, insect biodiversity and population abundances tend to be overlooked by the public, policy-makers and local and national authorities . Yet, recently, insect population declines have made the news, with multiple reports describing what appear to be dramatic losses in richness, abundance and biomass (Halmann et al., 2017(Halmann et al., , 2019(Halmann et al., , 2020Habel et al., 2019;Homburg et al., 2019;Seibold et al., 2019;Wagner, 2020). There is a recent call for urgent measures to be taken to avoid massive insect declines (e.g. Cardoso et al., 2020;Harvey et al., 2020;Samways et al., 2020) and obtain more longterm standardised robust data (Cardoso & Leather, 2019;Thomas et al., 2019). The many insect population reports that are arriving from across the world do however present a strong geographical bias, focusing on the European or North American mainland. The rest of the world remains much less studied and to what point the current findings are general is unknown.
The present biodiversity crisis is hitting islands disproportionately (Kier et al., 2009), with impacts at all levels of biodiversity but also affecting human populations and economies Weiss, 2015). This biodiversity crisis on islands and elsewhere is a consequence of several biodiversity erosion drivers including land-use changes, habitat degradation, pollution, invasive species and climatic changes . In fact, historical human colonisation of most isolated oceanic archipelagos promoted whole island ecosystem extirpation, and at best, native habitats fragmentation. Despite evidence of historical losses in species diversity on islands after Human occupation (e.g. Goodfriend et al., 1994, Blackburn et al., 2004Triantis et al., 2010;Alcover et al., 2015;Kirch, 2015;Terzopoulou et al., 2015), a key gap is to investigate the recent temporal trends on species diversity and composition changes on island native faunas.
In this contribution, we use a unique time-series dataset of arthropods sampled in Azorean native forest over 6 years (2013-2018) (see Borges et al., 2017b;Matthews et al., 2019b) to investigate recent trends in arthropod diversity and composition. The Azores is one of the world's most isolated archipelagos made of nine main islands aligned on a WNW-ESE axis in the Atlantic Ocean. At the time of human colonisation, around AD 1440, the Azorean archipelago was almost entirely forest-covered. In less than 600 years, 95% of the original native forest has been destroyed  due to the development of an economy based on agriculture (intensive dairy production) and forestry (the plantation of forests of fast growing exotic species; Cryptomeria japonica). The protection of the remaining native forest fragments became recently a priority for local authorities, since those fragments play a fundamental role as refuges and source habitats for the endemic fauna and flora (Gaspar et al., 2011).
Several negative biodiversity erosion drivers are currently impacting Azorean biota . For instance, Triantis et al. (2010) investigated the possibility of an extinction debt attributable to the forest destruction in Azorean arthropod fauna and have estimated that more than half of the extant native forest-dependent arthropod species might eventually be driven to extinction. The current lack of connectivity between forest patches (Aparício et al., 2018) and potential enhanced habitat fragmentation due to climate change (Ferreira et al., 2016;Aparício et al., 2018) is reducing the habitat quality for many arthropod species. Based on this previous work and despite recent studies showing that there is a dominance of endemics (Ribeiro et al., 2005;Rego et al., 2019) and possible resistance to invasions (Florencio et al., 2016), we expect a visible decrease in endemics diversity over time most markedly in abundance. Moreover, since exotic arthropod species are not in equilibrium and consequently are more dynamic in space (Cardoso et al., 2009;Rigal et al., 2018) and time (Matthews et al., 2019b) than indigenous species, we also expect a temporal positive signal for this group of species. Given the expected loss of endemics but gain of exotics, overall richness and abundance would remain similar across time. We test these predictions using the diversity metrics (Hill numbers), change in species composition over time and species abundance distribution (SAD) for all arthropod species, endemics, native non-endemics and exotic species.

Study sites and sampling methods
The Azorean archipelago is located in the North Atlantic, roughly between 37 -40 N latitude and 25 -3 W longitude. Terceira Island (Fig. 1), a roughly circular island with an area of 402 km 2 , contains the largest area of native pristine forests in the Azores (Triantis et al., 2010), with five main fragments of native forest distributed across four main volcanic polygenetic complexes. The study was initiated in 2012 as part of the EU-NETBIOME-ISLANDBIODIV (Cicconardi et al., 2017;Borges et al., 2018a) (Fig. 1). SLAM (Sea, Land, and Air Malaise) traps (see Supporting Information Fig. S1) were setup in 10 plots within native forest for long-term monitoring. These traps are passive flight interception traps with approximately 110 × 110 × 110 cm, similar to Malaise traps, but allowing the interception of insects from all directions (arthropods hit an area of black mesh and are funnelled into a sampling bottle). SLAM trapping is a suitable technique for long-term monitoring, but these traps still only collect a more mobile set of the invertebrate fauna. The collecting bottles, containing propylene glycol as preservative, operated continually for 6 years, being collected and changed every 3 months; thus, each sample covers one season of the year. Two of the original ISLANDBIODIV plots were not considered due to the fact that some samples were missing in some of the years. In total, 192 (8 plots × 4 seasons × 6 years) samples were considered for the study. A more comprehensive outline of the study is provided by Borges et al. (2017b) and Matthews et al. (2019b).
Arthropods were grouped into three colonisation categories consulting the last updated list of Azorean arthropods (Borges et al., 2010): endemic (i.e. restricted to the Azores); native non-endemic, i.e. species that arrived naturally to the archipelago but are present both in the Azorean Islands and elsewhere; and exotic non-native species, i.e., species whose original distribution range did not include the Azores and that are believed to have been introduced in the Macaronesian region after human settlement in the 15th century. The exotic status was inferred either from historical records of detected species introductions or from their current distribution being closely associated with human activity. For unidentified species, if other species in the same genus, subfamily or family were present in the archipelago and all belonged to the same colonisation category (according to Borges et al., 2010), the unknown species were classified similarly. Otherwise, we assumed the species to be unclassified.

Data analysis
For the current study, we used data sampled over the years of 2013-2018 (inclusive) and pooled the data from the four seasons of each year to create yearly datasets for each of the eight plots considered (temporal α diversity). We further pooled the eight plot trap data to have yearly datasets at the regional scale (temporal γ diversity). Unless otherwise stated, all the following analysis were performed with temporal γ diversity (one value per year). A sensitivity analysis was performed previously by Matthews et al. (2019b) to ensure our results were robust, namely sampling completeness estimates for each year in each plot and an evaluation if a single SLAM trap was sufficient to capture the relevant community properties, evaluating the results of three SLAM traps setup in one of the plots during 1 year. In both cases, the data were considered robust.
Our first set of analyses aimed to test whether species diversity of arthropod assemblages show temporal trends over 6 years period of sampling. We calculated the Hill numbers at three different orders (q) of diversity with a q value of 0 for species richness with all species having the same weight, q = 1 for the exponential of Shannon's index with species being weighted exactly for their abundance in the community and whereas q = 2 for the inverse of Simpson's index which favours abundant species. We used a sample size-based rarefaction approach to estimate the rate of increase in Hill number with increasing number of individuals sampled and then extrapolated the observed accumulation curve using the rarefaction-extrapolation approach (Chao & Jost, 2012). Second, we estimated accumulation rates in relation to the degree of completeness of sampling effort, rather than just sample size, as recommended by Chao and Jost (2012). Degree of completeness for each year was evaluated using the coverage estimator (Ĉ n ; Chao & Jost, 2012) which estimates the proportion of the total number of individuals in an assemblage that belong to the species represented in the sample. Sample coverage is considered to be an objective measure of completeness (Chao & Jost, 2012). Hill numbers were subsequently calculated for endemic, native non-endemic and exotic species separately. Generalised linear models with quasi-Poisson error and linear regression models with log-transformed data were used to examine the relationships between diversity and years for species richness and the exponential of Shannon's index and the inverse of Simpson's index, respectively. This was carried out for all species and for endemic, native non-endemic and exotic species. All aforementioned analyses were also performed by averaging diversity across the eight plots (mean temporal α diversity) to test whether trends found at regional scale were congruent or not with patterns found at plot scales.
Our second set of analyses was to evaluate temporal trends in abundance per individual species across 6 years. We only selected species occurring in a minimum of 4 years out of 6 years and with a minimum of five individuals per year of occurrence. We used this cut off because temporal changes of very rare species could be caused by minor, chance events. Generalised linear models with quasi-Poisson error were run to test for the relation between the abundance and years.
Our third analysis aimed to examine β diversity patterns over time considering all species together but also for endemic, native non-endemic and exotic species separately. We used the Bray-Curtis index as a β diversity measure with square-root transformed abundance. Relationships between β diversity and distance in time were tested with Mantel permutation tests, based on 9999 permutations and using Spearman's correlation coefficient. We also tested for differences in β diversity values between endemic, native non-endemic and exotic species using Kruskal-Wallis tests followed by post hoc pairwise Wilcoxon tests.
Our fourth set of analyses was to evaluate temporal changes in SAD to further quantify the contribution of endemic, native nonendemic and exotic species abundance patterns over time. We used the Gambin model (Ugland et al., 2007, Matthews et al., 2014 to characterise the shape of the SAD (Ugland et al., 2007). The unimodal gambin model has a single free parameter (α), which characterises the distribution shape with low values indicating logseries-distribution and higher values indicating more lognormal-distribution. Gambin has been shown to provide good fits to a wide variety of empirical datasets (Matthews et al., 2014(Matthews et al., , 2019a. Recent developments (Matthews et al., 2019a) have derived the likelihood functions for multimodal gambin models (i.e. with one α per mode) providing a means of easily assessing multimodality in SAD datasets. Therefore, for each year, we fitted both unimodal and bimodal gambin models and used the Bayesian information criterion (BIC) values to compare both models (Burnham & Anderson, 2002). The bimodal model was considered to better fit the SAD if its BIC value was lower than the unimodal model. We also recorded the χ 2 goodness-of-fit statistic and its associated P-value for both one and bimodal models. Given that SAD model parameters are sensitive to variations in sample size (McGill, 2011) and that we were also interested in comparing parameter values across years, we used a procedure where, for each sample, we subsampled 3000 individuals (the least abundant sample having 3944 individuals), fitted the best gambin model (either uni-or bimodal) to this subsample and stored the α parameter value(s). Given that this subsampling procedure is stochastic, we repeated the process 100 times for each sample and took the mean α value. Mean α values were regressed against years using linear regressions. Finally, we calculated the mean weighted octaves (i.e. weighted with species number) for endemic, native non-endemic and exotic species per year and (i) we used weighted t-tests to examine whether mean weighted octaves differ between groups in each year and (ii) run linear regressions to test whether, for a given group, mean weighted octaves increase or decrease over time. Statistical analyses were implemented within the R programming environment (R Development Core Team, 2014) using the packages vegan (Oksanen et al., 2013), iNEXT (Chao et al., 2014), BAT (Cardoso et al., 2015) and gambin (Matthews et al., 2014).

Results
A total of 30 875 arthropod specimens were collected and 159 (morpho) species were identified, representing 17 orders. Of the 159 species, 32 were considered to be endemic to Azores, 63 to be native (non-endemics) and 57 to be exotics while seven (7) species were not assigned to any group. About 10% of the taxa were identified to morphospecies rather than species. The number of species and total abundance per year for all species and for the three species groups are presented in Table 1. Overall, we found marginal significant differences in species richness over time between the three groups (Friedman's test, P = 0.04) with exotics representing a larger fraction of diversity over time (paired Wilcoxon post hoc test P = 0.036).
Sample size-based species accumulation curves (Supporting Information Fig. S2) indicated that the arthropod fauna was well sampled for each year but did not represent a complete census of species. The majority of these unsampled species are likely to be rare (and potentially vagrants) as coverage-based rarefaction revealed a sample completeness of almost 100% for 6 years, suggesting that very few individuals of the community belonged to species not represented in our samples. Therefore, we conclude that our sampling was satisfactory for characterisation of Table 1. Degree of completeness estimated with coverage (Ĉ n ) and number of species and individuals for all species (ALL) and for the endemic (END), native (non-endemic) (NAT) and exotic (EXO) species separately. The coverage estimator estimates the proportion of the total number of individuals in an assemblage that belong to the species represented in the sample.  Fig. S2; Table 1). Overall, we only detected a significant increase of gamma diversity with time for exotic species, for both the exponential of Shannon's index (R 2 = 0.81, P = 0.01) and the inverse of Simpson's index (R 2 = 0.84, P = 0.009) ( Fig. 2; see detailed results in the Supporting Information Tables S1 and S2). Similar results were obtained with mean α diversity (exotic species: exponential of Shannon's index; R 2 = 0.86, P = 0.007 and inverse of Simpson's index; R 2 = 0.86, P = 0.007) (See detailed results in the Supporting Information Table S2).
Only 48 species were retained for the individual abundance analysis, representing 31% of the total species recorded in our study but 96% of the individuals (Fig. 3). Overall, average effect size was −0.064 ± 0.41, −0.066 ± 0.417 and − 0.066 ± 0.249 for endemic, native non-endemic and exotic species respectively. Nevertheless, 11 species showed a significant decline in abundance, while only three species showed a significant increase (Supporting Information Table S3).
We found significant and positive relationships between β diversity with distance in time for overall, endemic and exotic species (r s = 0.85; P = 0.002; r s = 0.71; P = 0.016 and r s = 0.60; P = 0.011, respectively; Fig. 4) but not for native non-endemic species (r s = 0.20; P = 0.247). On average, β diversity values were strongly and significantly different between the three species groups (Kruskal-Wallis; H = 36.25; P < 0.001) with exotic species having the highest β diversity over time and endemics the lowest (post hoc Pairwise Wilcoxon test: all pairs being significantly different, P < 0.001; Fig. 4). Figure 2. Relationships between species diversity and years. Analyses were performed with first three Hill numbers namely species richness, the exponential of Shannon and the inverse Simpson and for all species and endemic, native (non-endemic) and exotic species separately. Generalised linear model with quasi-Poisson error and linear regression models with log-transformed data were used to examine the relationships between diversity and years for species richness and the exponential of Shannon's index and the inverse of Simpson's index, respectively. Dotted red lines indicate non-significant relationships while solid red lines indicate significant ones.
For each of 6 years, SADs were best represented by a bimodal gambin model according to the χ 2 goodness-of-fit test and BIC ( Fig. 5; see Supporting Information Table S4). Non-significant relationships were detected between the mean α parameter values and time for both modes (R 2 = 0.15, P = 0.44 and R 2 = 0.04, P = 0.68). Our weighted t-tests showed that, regardless of the years, exotic species were always significantly associated to smaller octaves (rare species, Supporting Information Table S5) than endemic species that were mostly represented in higher octaves (dominant species). We did not detect any significant temporal trends in the weighted mean octaves for endemic, native nonendemic or exotic species (R 2 = 0.04, P = 0.68, R 2 = 0.01, P = 0.85 and R 2 = 0.003, P = 0.91, respectively).

Discussion
Contrasting with recent studies on the European and North American mainland, our study shows no evidence of a decline of overall species diversity for endemic and native non-endemic species, although the data are limited to 6 years of sampling. Yet, our study clearly identified three important patterns: (i) an increasing diversity of exotic species over time; (ii) the high temporal dynamism of exotic species (see also Matthews et al., 2019b) and (iii) some evidence of a tendency for decreasing abundance for some endemic species.
Therefore, all our predictions were generally confirmed, and here we clearly identify exotic species as a potential driver of Figure 3. Effect sizes of the relationships between individual species abundance with years. Effect size refers to the standardised slope of the relationship extracted from the Generalised linear models with quasi-Poisson error fitted to the data. Analyses were performed only for the species occurring in a minimum of 4 years out of 6 years and with a minimum of five individuals per year of occurrence. In total, 48 species were considered representing 31% of the species sampled across 6 years but 96% of all individuals. Horizontal lines indicate the confidence interval at 95% calculated for each effect size. Negative effect size indicates decline in abundance over 6 years while positive effect size indicates an increase in abundance. Endemic, native (non-endemic) and exotic species per octave are indicated with different colours. Species for which the effect size was significant are labelled. future biodiversity erosion on Azorean native forests. This is partly similar to the pattern found by Seibold et al. (2019) in German forests, where few initially abundant invasive (and potential pest) species increased in abundance. In previous studies we observed that exotic species are already permeating Azorean indigenous arthropod communities (Gaston et al., 2006;Rigal et al., 2013), being however only abundant in anthropogenic habitats (Cardoso et al., 2009). Rapid demographic fluctuations (Matthews et al., 2019b) and contingent spatial factors related with native forest fragmentation and land-use changes (Cardoso et al., 2009;Florencio et al., 2016) are surely involved on the increased diversity through time of exotic species in the Azores. It is worth noting that the increasing diversity of exotic species in 6 years is mostly due to juveniles. This suggests the increasing importance of a source-sink dynamics between natural forests as sinks and other habitat types as sources. The dispersion of juveniles from the surrounding matrix of semi-natural grasslands and exotic plantations of Cryptomeria japonica towards native forests seems to be more significant from year to year (see also Borges et al., 2008). Despite previous observations suggesting that Azorean canopies may act as physical barriers to the colonisation of exotic spiders (Florencio et al., 2016), the results obtained in this six-year time-series study are less optimistic. If this process continues, we may expect that in addition to the ongoing invasion of soil communities (see Cicconardi et al.;, the Azorean arthropod canopy communities can also undergo an ecological meltdown and homogenisation process. The observation that there is a tendency for a decreasing abundance of some endemic species is also of some concern. A recent assessment of the IUCN red list status of Azorean endemic forest arthropods (Borges et al., 2017a(Borges et al., 2018c; see also http://www. maiisg.com/), concluded that most species are threatened. There is no current easy explanation for the observed trend in the endemic arthropod species, but it can be associated with the current spread of invasive plants that are changing dramatically the cover and structure of forest understory, namely the cover of ferns and bryophytes in the forest floor (Borges et al., 2017a(Borges et al., 2018c. In addition, António Frias Martins (pers. comm.) observed a tendency for a decrease in the cover of bryophytes Figure 4. Results of the β diversity analyses over time for all species, endemic, native (non-endemic) and exotic species. β diversity was computed with the Bray-Curtis index. (a) Relationship between β diversity and distance in years. Relationships were tested using Mantel test with Spearman correlation. Dotted red lines indicate non-significant relationships while solid red lines indicate significant ones. (b) Distribution of the β diversity values between endemic, native and exotic species. Differences were tested with Kruskal-Wallis test. The test was significant and post-hoc pairwise Wilcoxon test showed that all groups differed between each other. on the canopy of some endemic Azorean trees at high elevations in the last 20 years with a consequent decrease in the abundance of endemic micro land snails. This observation indicates an ongoing impact of climate change on Azorean native forest, resulting in a change in habitat structure that alters the diversity and abundance of endemic invertebrates.
Continuing monitoring of island forest ecosystems is badly needed (Borges et al., 2018b), on the Azores as well as other regions to avoid unexpected declines and extinctions of native, often endemic, species.

Conclusions
In this study, no major insect diversity decline was observed in the native forests of a small Atlantic volcanic island in Azores (Portugal). In contrast, in addition to some evidence of decreasing abundances of endemic species, the increase through time of the diversity of exotic arthropods was the stronger pattern.
We can consider that our results are representative of the current situation on Azorean native forests since we included most of arthropod orders encompassing several trophic levels, even if excluding Acari, Collembola, Lepidoptera, Diptera and Hymenoptera from our diversity measurements. The impact of exotic species on islands is well known (Sax & Gaines, 2008), and as a result of this 6-year study on Azores, important management actions should be taken to tackle this problem (see also Harvey et al., 2020;Samways et al., 2020). This is therefore a different, complementary view of many studies showing insect population declines. On islands, as studied extensively in the past, invasive species are particularly relevant in a conservation context. This signal would have been missed if species were not identified and classified according to their origin, masking a worrisome pattern.