Long term changes (1990-2016) in carabid beetle assemblages (Coleoptera: Carabidae) in protected forests on Dinaric Karst on Mountain Risnjak, Croatia

Carabids, as well-known bioindicators, have been used to study the long term changes that have occurred in their communities in the Dinaric Alps. This study involved eight sites in the protected forests of the Risnjak National Park in the years 2015 and 2016 of which three were previously studied in 1990 and 1991. A total of 9, 521 individual ground beetles belonging to 17 genera and 33 species were collected. Species diversity and community composition, including percentages of species grouped according to their habitat preferences, body size, wing morphology, preferred moisture and temperature were used to compare the sites sampled in 1990 and 1991 and resampled in 2015 and 2016. Even though this study was carried out in protected forests within the National Park with minimal anthropogenic pressure and the fact that available climatic data didn't show any significant change in climate over the last 25 years, there was a reduction in the abundance of specialist species and increase in the spread of generalist species. Furthermore, the lower abundance of a mountain specialist and endemic species, Pterostichus variolatus, and the lack of mountain specialists Molops alpestris, Pterostichus unctulatus and Trechus croaticus in the catches indicate the importance of further monitoring of these mountain forest ecosystems and for a well-timed and appropriate conservation approach.

Long-term data could help reveal the extent of changes in species diversity, ranges, phenology of life cycles and interactions following environmental changes (such as cli-Long term changes  in carabid beetle assemblages (Coleoptera: Carabidae) in protected forests on Dinaric Karst on Mountain Risnjak, Croatia ŽELJKA JAMBROŠIĆ VLADIĆ 1 and LUCIJA ŠERIĆ JELASKA 2

INTRODUCTION
The rapid loss of biodiversity is a threat to the stability and the existence of the ecosystems we know today. While some claim that we are in the middle of the largest extinction in the Earth's history -Sixth Mass Extinction (Ceballos et al., 2015;Payne et al., 2016), others hope that the situation is overstated (Leather, 2017). Nevertheless, the loss of any form of life in some way affects the functioning of ecosystems (Bradley et al., 2012).
Of the eight sites sampled (L 1 -L 8 ), three were surveyed in 2015 and 2016 (L 3 , L 5 , L 8 ) and fi ve additional sites were chosen in 2016. Sites were at altitudes ranging from 704 m a.s.l. to 1277 m a.s.l. (Fig. 1, Table 1), with annual rainfall from 1919.3 mm yr -1 to 3708 mm yr -1 (values for period between 1990 and 2016), mean maximum annual temperatures from 20.07°C to 22.2°C, and mean minimum annual temperatures from -7.29°C to -2.34°C (values for period between 2006 and 2016).
The forest here is dominated by beech and includes mixed beech-silver fi r forests at the lower altitudes and subalpine beech forest at the higher altitudes. Dwarf mountain pine (Pinus mugo) bushes occur at the treeline, above which there are subalpine calcareous grasslands.
The experimental sites (L 3 , L 5 , L 8 ) were selected within forest habitats on the upper half of the mountain right up to the top of the tree line. These three sites were at the same location as those used in the research on carabid communities carried out in 1990-1991 (Table 1). Since we did not have information on the exact location of pitfall traps at each location (e.g. geographical coordinates), traps were positioned within forest based on the descriptions of the previous fi eld study. The sampling methods were the same as in 1990-1991.

Sampling protocol
Ground beetles were sampled between May and October, during 2015 and 2016, using pitfall traps at all eight sites (three sites in 2015 and 2016 plus an additional fi ve sites in 2016). Plastic vessels with a diameter of 9.5 cm and a depth of 14.0 cm (0.5 L) were used as pitfall traps. In 2015 six pitfall traps per site were set (site codes: L 5 , L 6 , L 7 ), and in 2016 nine pitfall traps were set at each of the eight sites (site codes: L 1 -L 8 ). Thus, a total of 90 traps were used during two years. At each site, traps were arranged in two rows with three traps in each row in 2015, and in a grid (3 × 3) in 2016, with about a 10-meter distance between each trap. The traps were fi lled up to one-third of their volume with a mixture of ethanol (96%), acetic acid (9%) and water in equal proportions (as was used in the study in 1990-1991) and reset every two to three weeks continuously throughout the sampling period (2 May 2015-10 Oct. 201512 May 2016-9 Oct. 2016. Pieces of tree bark were placed above the traps to protect them from rain. All specimens were identifi ed to species level (Trautner & Geigenmüller, 1987;Hůrka, 1996;Freude et al., 2004; Mueller collection in NHM Trieste) and stored either in a 70% alcohol solution or in a dry state.

Habitat variables
For each site, we recorded altitude (m a.s.l.), slope (degrees), aspect (i.e. the direction in which a slope faces; calculated as mate change, habitat destruction, etc.). Therefore, in this paper, we aimed to compare carabid beetle assemblages sampled in the years 1990 and 1991 with those sampled 25 years later (in 2015 and 2016), at the same sites. In addition to species richness, we analysed mean individual biomass (MIB) and changes in species traits (wing morphology, preferences for particular levels of moisture affi nity and temperature) assuming that the percentage of some traits within the community could have changed over the time indicating certain environmental modifi cations.
In this the fi rst study on carabids over a long period in the Dinaric Alps in Croatia, in addition to changes in carabid assemblages, we checked for a possible displacement of some species along an altitudinal gradient that could point to an increased level of environmental stress due to climate change.

Study area
This study was carried out in the region of the Risnjak mountain massif (1528 m) located in the Gorski Kotar area in western Croatia, in the northern part of the Dinaric Alps (Fig. 1). Its central part along with the Croatian Snježnik mountain range (1506 m) and the spring of the Kupa River (290 m a.s.l.) are in the Risnjak National Park. The forest communities within the Park are a  absolute number of degrees from south), and land cover (forest types: mixed Dinaric silver fi r and beech forest, and subalpine beech forest), ( Table 1). The fi rst three variables were calculated in a GIS environment, whereas the forest types were identifi ed in the fi eld based on the presence of particular species of plants.
The climate in the study area and how it has changed over the last few decades were analysed based on the data available from the Croatian meteorological and hydrological service (DHMZ). Precipitation data for the whole period  were available for Crni Lug (Risnjak National Park, at cca 800 m a.s.l.) while temperatures in the National Park were not measured before the end of 2003. Therefore, data for the nearest town Delnice (at 700 m a.s.l.) were used (Figs S1, S2). A non-parametric Mann-Kendall test was used to determine trends in the maximum and minimum monthly temperatures and total monthly precipitation ( Table 2).

Data analyses
In this study Margalef's diversity (DMg) and Menhinick's indices (DMn) were calculated as indicators of species richness, Simpsons' index as a dominance index, Shannon-Wiener's (H') as a diversity index and Pielou's (E) as an evenness index (Table  4).
To compare species composition between the sites in all the years sampled (1990,1991,2015,2016), we calculated Jaccard dissimilarity coeffi cient and grouped the data using Euclidean  distance as distance measure and complete linkage as linkage rule for constructing dendrograms (Krebs, 1989), (Fig. 2). Shapiro-Wilk and Levene's tests were used to check for normal distribution of data and homogeneity of variance. In addition, nonparametric Kruskal Wallis one-way ANOVA and Mann-Whitney U tests were used. Since there were no signifi cant differences in species composition and abundance among years (except for site L 3 in 1991 and 2016), we pooled data for 1990 and 1991 in one group and data for 2015 and 2016 in another group for further trait analyses. Similarities in species composition according to the percentages of tested species traits (wing morphology, moisture and temperature affi nity) recorded at the same sites in 1990-1991 and resampled in 2015-2016 were evaluated using repeated measures or one-way ANOVA. Principal Components Analysis (PCA) was used to explore variation in the distribution of the species traits among all sites. Pearson's correlation coeffi cient was used to measure the strength of the linear relationship between the percentages of the species traits with altitude. Multivariate linear regression analyses were used to determine the altitudinal distribution of the observed species traits.
Mean individual biomass (MIB) was calculated using a formula that describes the relationship between body length and live biomass of a carabid beetle: ln y = -8.928 + 2.555 × ln x (Szyszko, 1983). Data on the body length of the species caught were obtained from the Carabids.org -online datebase (Homburg et al., 2013). Differences in carabid biomass were standardise per trap per day for both of the periods 1990-1991 and 2015-2016, before further analyses.
Rarefaction curves were estimated for three of the sites in 1990-1991 and all of the sites in 2015-2016 (Fig. 3A, B).
The IBM SPSS v21; EstimateS v9.1. and Past 3 program packages were used for statistical calculations.

Changes in ground beetle assemblages 1990-2016
A total of 40 species were caught in 1990-1991 and 2015-2016 (sites L 3 , L 5  Analysing the composition and abundance of species in different years and sites, a statistically signifi cant difference was recorded between 1991 and 2016 for site L 3 (Shapiro-Wilk W: 0.2595, Levene's test for homogeneity of variance p < 0.01, Kruskal-Wallis p = 0.2799, Mann-Withney pairwise, p = 0.009), while in other years and sites they did not differ signifi cantly. The activity density (the total number of trapped specimens) of the most abundant species in 1990-1991 was for Nebria dahlii, which made up 48.99% of all the carabids caught. Based on their abundance the following species: Pterostichus variolatus (14.12%), Pterostichus unctulatus (11.97%), Pterostichus burmeisteri (6.58%) and Abax ovalis (6.38%), made up 88.04% of all the carabids caught.
In 2015-2016 the top fi ve of the most numerous species (sites L 3 , L 5 and L 8 ) were: Nebria dahlii (41.68%), Aptinus bombarda (34.09%), Pterostichus burmeisteri (8.2%), Comparing each of the three sites that were sampled in 1990-1991 and again in 2015-2016 all of the diversity indices (except DMg for site L 5 ) indicate that the diversity was higher in 1990-1991 than in 2015-2016 (Table 4), and based on species composition using Jaccard similarity index the sites sampled in 1990-1991 (L 3 , L 5 , L 8 ) formed one cluster and those sampled in 2015 and 2016 (L 1 -L 8 ) another cluster (Fig. 2).

Species trait analyses of the assemblages recorded between 1990-2016 along an altitudinal gradient
The analyses of variance (repeated measures ANOVA) did not reveal statistically signifi cant differences in carabid assemblages based on wing morphology in the two periods compared (1990-1991 and 2015-2016). Because of the zero variance between some sites, the difference in moisture preference in the two periods sampled could not be tested using factorial ANOVA. Instead, one-way ANOVA was used and did not reveal a statistically signifi cant differences [F (1,13) = 0.29, p = 0.5946], but some trends were recorded. The most numerous species, in both periods, were brachypterous and hygrophilous or moisture indifferent. Differences were recorded in temperature preference; in 1990-1991 the fi rst four of the most numerous species were low-temperature species (N. dahlii, P. variolatus, P. unctulatus, P. burmeisteri) while in 2015-2016 two of the four the most numerous species (N. dahlii, A. bombarda, P. burmeisteri, C. violaceus) were thermophilous, including the second most abundant. Site L 3 in 1990Site L 3 in -1991Site L 3 in and 2016 Based on the ecological traits of the different species, three macropterous species of which two were xerophilic (Anisodactylus intermedius, Amara eurynota, Leistus spinibarbis) were caught in 2015-2016, whereas in 1990-1991 no macropterous or xerophilic species were caught (Table  3).
A big increase in the abundance of Aptinus bombarda was recorded, from 4.71% in 1990-1991 to 43.72% in 2016, which is an increase of 42 times more specimens caught per day. Site L 5 in 1990Site L 5 in -1991Site L 5 in and 2015Site L 5 in -2016 In 2015-2016, among the top fi ve most numerous carabids caught there were two thermophilous species (Carabus violaceus, Molops striolatus), whereas in 1990-1991 only low-temperature species were dominant (Pterostichus unctulatus, Nebria dahlii, Pterostichus variolatus, Pterostichus burmeisteri), (Table 3). Site L 8 in 1990Site L 8 in -1991Site L 8 in and 2016 In 1990-1991 the top fi ve most numerous species were low-temperature species whereas in 2016 the most numerous species was thermophilous (Aptinus bombarda), (Table  3). In addition, the endemic species Pterostichus variolatus decreased in abundance from 18.06% to 1. 7% between 1990-1991 and 2016. A high positive correlation between the percentage of low-temperature species and increase in altitude was recorded (r = 0.76). Pearson's correlation coeffi cients also showed a positive correlation between large (r = 0.54), brachypterous (r = 0.33) and moisture-indifferent species (r = 0.47) with increase in altitude, indicating that large, brachypterous, low-temperature and moisture-indifferent species are more likely to occur at high altitudes. The fi rst two axes of the PCA graph (Fig. 4) represent 66.95% of the total variance (44.44% and 22.5%, respectively). The fi rst axis separates sites L 1 , L 2 and L 3 (sampled in 2015 and 2016) at low altitudes on the right side of the ordination space based on their higher percentage of thermophilic, xerophilic and macropterous species, and sites sampled in 1990 and 1991 (L 3 , L 5 , L 8 ) with smaller and polypterous species. On the left side of the ordination biplot, there are the sites L 4 , L 5 , L 6 , L 7 and L 8 (sampled in 2015 and 2016) at high altitudes, with a high percentage of low-temperature, moisture-indifferent and large brachypterous species. Along the second ordination axis, there are sites L 5 and L 8 (sampled both in 1990-1991 and 2015-2016) that are separated on the basis of their higher percentage of low-temperature species from those sampled in 1990-1991 with more polypterous and smaller species.
Distribution of the traits analysed (Fig. 5) along the altitudinal gradient revealed some general trends in terms of brachypterous, large, cold-adapted species being more abundant at high altitudes.
Analysing mea n individual biomass recorded at the sites that were compared revealed a statistically signifi cant difference for site L 3 (MIB 1990(MIB -1991  Based on the percentage of different species at given altitudes we noticed a possible displacement of some species along the altitudinal gradient (Fig. 6). Notiophilus biguttatus, a mesothermophilous and hygrophilous species whose altitudinal range in 1990-1991 was from 742 m to 1199 m, in 2015-2016 was not recorded above 1068 m; Licinus hoffmannseggi, a low-temperature and hygrophilous species, in 2015-2016 seemed to expand its range at both low and high altitudes (from 742 m to 1277 m) recorded in 1990-1991, but at only one site L 6 (1060 m); Stomis rostratus, a mesothermophilous and moisture indifferent species, seemed to reduce its range in 2015-2016 and was caught only at site L 8 (1199 m) while 25 years ago it was recorded at altitudes from 742 m to 1199 m (sites L 3 , L 5 , L 8 ). The expected increase in abundance of cold-adapted endemic species with altitude was recorded for the species C. creutzeri and P. variolatus. Both these species were recorded at altitudes from 1068 m to 1277 m and were the most numerous at the highest site (L 7 , 1277 m). This was also recorded for another endemic mesothermophilous species, C. croaticus. It was recorded at altitudes from 742 m to 1277 m and was by far the most numerous (20% of the total catch) at the highest site (L 7 ).

DISCUSSION
Analysing the species composition of carabids assemblages and distribution of species traits between communities sampled at the same sites in 1990-1991 and again 25    Table 1, and abbreviations of carabid names correspond to those in Table S1. years later, we revealed an increase in the number of thermophilous species and decrease in the number of species preferring low temperatures at some of the sites, among them some mountain specialists (e.g. P. unctulatus) and species typical of the Alpine-Dinaric region (e.g. P. variolatus). This is an expected consequence of global warming and supported by the fact that there was a slight increase in the minimum temperature recorded on Mt. Risnjak and in the Gorski Kotar region over the last 25 years (Crni Lug: m = 0.61.; Delnice: m = 0.26). The most obvious changes were recorded at site L 8, which is located in a subalpine beech forest, where in 1990-1991, more than 60% of the total catch consisted of species prefering low-temperatures and high humidity, such as Nebria dahlii, while the most numerous species at the same site 25 years later was Aptinus bombarda (70.1% of the total catch), which is a thermophilous and mesohigrophilic species. In addition, a big decrease in the low-temperature prefering species Pterostichus variolatus was also recorded at this site (an 11 fold decrease in the number of specimens caught per day). The lack of the mountain specialist, Pterostichus unctulatus, in the catches in 2015-2016 and the fact that this species made almost 12% of the total catch in 1990-1991 also indicate a decline in the abundance of low-temperature preferring species in this protected area. A similar decreasing trend in the abundance of low-temperature preferring species is recorded in long term studies in the Alps (Pizzoloto et al., 2014) and High Tatra Mts, where in 2004 after a windstorm, at some sites, less tolerant species, including P. unctulatus, disappeared (Šustek, 2013).
Certain changes can be associated with ice storms, which are characterized by freezing rain, which in February 2014 destroyed more than 680,000 m 3 of timber in the National Park (Editorial Board Šumarski list, 2014). This natural disturbance and the removal of some of the damaged trees could have stressed some parts of the ecosystem. Thinning of the forest and changes in primary habitats might account for the increase in the number of thermophilous species given that there were no signifi cant changes in the climate in this area. But the existence of only one climatological station in the National Park is insuffi cient for precise monitoring of microclimatic changes in this area.
One of the most important morphological traits of carabids used for bioindication purpose, wing morphology, clearly defi nes the ecology and evolution occuring in a certain area (Darlington, 1943;Wagner & Liebherr, 1992). Defi ning the percentage of species of carabids that are wingless has proved to be a powerful method for determining the stability of an area (Dhuyvetter et al., 2007;Wagner & Liebherr, 1992;Szyszko et al., 2000;Gobbi et al., 2007;Šerić Jelaska & Durbešić, 2009;Pizzolotto et al., 2016). Therefore, we used it to compare the sites sampled in 1990-1991 and again 25 years later. The high abundance and percentage of brachypterous species (62.5% of the species, 99.5% of the individuals) caught in 2015-2016 and the absence of differences in wing morphology recorded at these three sites then and in 1990-1991 (brachypterous species: 1990-1991: 79.3% of the species, 96.8% of the individuals; 2015-2016: 75.9% of the species, 99.2% of the individuals) may indicate that the forest ecosystems in the National Park in general are stable.
Risnjak National Park as a part of the Dinaric Mountains chain has many endemic species some of which were recorded in this study (N. dahlii, C. catenulatus, C. croaticus, C. creutzeri, M. ovipennis, P. unctulatus, P. variolatus, S. rostratus). Most of them are mountain habitat specialists and low-temperature preferring species. Although some species may benefi t from increases in availability of habitats by moving to higher altitudes, the Dinarides are socalled "mountain pyramids" on which there is a decrease in the availability of habitats with increase in altitude (Elsen & Tingley, 2015). Ultimately this can cause the local extinction of affected species (Colwell et al., 2008;Sekercioglu et al., 2008Sekercioglu et al., , 2012La Sorte & Jetz, 2010).
The displacement of low-temperature species to higher and lower altitudes and that of mesothermophilous species may be due not only to extrinsic characteristics that affect the displacements' paths of certain species, but also species' intrinsic characteristics on which the response to environmental disturbances depends .
It is expected that most endemic high mountain coldadapted specialist will tend to move upwards with increase in the average annual temperature (Wilson et al., 2005;Pizzoloto et al., 2014;Scalercio et al., 2014). Although we do not have enough data to verify a possible change in temperature over the last 25 years, we noticed that two endemic low temperature species (C. creutzeri, P. variolatus) and one endemic mesothermophilous species (C. croaticus) are now the most numerous at the site located at the highest altitude.
According to Nolte et al. (2019), mountain habitat species are at greater risk of extinction than forest species. Thus, vegetation cover may have a mitigating role and maybe extend the time species have to adapt to new conditions. Consequently, any destruction of forest ecosystems, either by man (thinning) or indirectly by unpredictable weather conditions due to climate change, will impair habitat quality and diminishes its protective role in species conservation. Species at greatest risk of extinction include specialist predatory species, those with latitudinal and altitudinal restrictions and brachypterous species (Feehan et al., 2009). Due to stenovalence and inability to cope with sudden landscape changes (Gaston & Fuller, 2009) specialists are considered to be the most vulnerable at times of evident climate change (Kotze & O'Hara, 2003;Terzopoulou et al., 2015). Penev (1996) points out that the composition of carabid assemblages is one of the most revealing ways of detecting changes in ecosystems. It contributes to the discovery of the effects of climate change and indicates the importance of regular and long-term monitoring (Kerr et al., 2007;Vaibhao et al., 2013). Current methods of monitoring insects assemblages require a lot of effort, so it is not surprising that there is only a small number of such long-term studies (Shortall et al., 2009;Fuentes-Montemayor et al., 2011;Dirzo et al., 2014). This is unfortunate as it makes it dif-fi cult to study the changes occurring in insect assemblages with changes in the environment. This is the case in the Risnjak National Park, which represents a natural link between the Alps and the Balkans and the most signifi cant example of the separation of Croatia (np-risnjak.hr). The fi rst extensive research in the fi eld in this National Park was carried out in 1963-1964(Durbešić, 1967 at six sites, then almost 27 years later at four different sites (Vujčić-Karlo, 1999) and then this study at eight sites.
Prediction of species-specifi c responses to new circumstances isn't possible without intensive and regular monitoring of carabid assemblages, especially bearing in mind the increase in anthropogenic infl uences and unpredictable weather conditions as a consequence of climate change. The protection of an area without further monitoring doesn't mean much in terms of protecting high mountain endemic species as they have already suffered to some extent the consequences of quick climatic changes as mountain ecosystems are among the fi rst to suffer the consequence of global warming. Thus, despite the study of protected forests within the National Park, there was a decrease in the number of specialist species and the spread of generalist species. The necessity for continuous monitoring, in the current period of obvious climate change, is becoming unquestionable (Saunders, 2017). The high number of endemic and threatened species in mountain ecosystems such as the Dinarides requires monitoring as it is the only way to recognise and mitigate changes that can lead to the reduction of biodiversity and permanent loss of some endemic species.