Monitoring the abundance of saproxylic red-listed species in a managed beech forest by landsat temporal metrics

Frncesco Prisi, Eli Vngi,,＊, Sverio Frncini,c, Gherro Chirici,c, Dvie Trvglini,Mrco Mrchetti, Roerto Tognetti

a GeoLAB-Laboratorio di Geomatica Forestale,Dipartimento di Scienze e Tecnologie Agrarie,Alimentari,Ambientali e Forestali,Universit`a degli Studi di Firenze,Via San Bonaventura 13, 50145, Firenze, Italy

b Dipartimento di Bioscienze e Territorio, Universit`a degli Studi del Molise, Contrada Fonte Lappone, 86090 Pesche, Isernia, Italy

c Fondazione per il Futuro delle Citt`a, Firenze, Italy

d Dipartimento Agricoltura, Ambiente e Alimenti, Universit`a degli Studi del Molise, Via De Sanctis, I-86100, Campobasso, Italy

Keywords:Alpha-diversity Beetle assemblages Biodiversity hotspots Normalized difference vegetation index Temporal metrics

ABSTRACT

1. Introduction

Saproxylic beetles, i.e., organisms that depend upon wounded or decaying woody material, during at least a stage of their life cycle,represent the most threatened animal assemblages in European forests(Nieto and Alexander,2010;Stokland et al.,2012),represent the highest percentage of world forest diversity(about 34%is related to deadwood),and act as a keystone in forest dynamics.Indeed,since they are involved in the degradation of wood, they contribute to soil fertility and the incorporation of nutrients in ecosystems (Buse et al., 2009; Mic′o et al.,2011;Parisi et al.,2018).Finally,influencing deadwood abundance and quality,saproxylic beetles have a significant impact on biodiversity and other saproxylic organisms that depend on deadwood.On the other hand,saproxylic beetles have been widely used to study changes in the quantity and quality of deadwood(Alexander,2010;Stokland et al.,2012)and to evaluate the forest management effects on compensating the loss of biodiversity and ecosystem services(Lindenmayer et al.,2000).

Forest structure, tree species composition, deadwood volume, abundance, and diversity of microhabitats are essential factors influencing biodiversity in forest ecosystems(Paillet et al.,2010;Müller et al.,2012;Regnery et al., 2013; Parisi et al., 2020a). Indeed, many saproxylic species depend on tree microhabitats for food,shelter,and breeding(Vuidot et al.,2011;Larrieu et al.,2018;Parisi et al.,2020b).Although ecological studies on saproxylic beetle communities in southern Europe have increased in recent years(Gossner et al.,2016;Lachat et al.,2016;Lelli et al.,2019;Parisi et al.,2016,2019,2020a,2020b;Vogel et al.,2020a,2020b), there are relatively few studies on saproxylic organisms in managed mountain forest ecosystems (Siitonen, 2001; Parisi et al.,2020a). Most research is concentrated in unmanaged forests and protected areas(Parisi et al.,2021;Sabatini et al.,2016).

While saproxylic beetles are relevant indicators for studying different aspects of biodiversity, most are small (many species of beetles do not exceed 2 mm), cryptic, and difficult to sample (Bouget et al., 2008).Together with the high number of species and the multiple trophic roles they occupy, these characteristics make their study remarkably challenging.Indeed,monitoring these groups of organisms has high costs.For this reason,in general,studies on saproxylic beetles focus on small areas,which are inappropriate to derive information on sustainable management and conservation practices to the spatial scale of forest properties.Furthermore, a multi-year sampling would be appropriate to follow population trends over time, with all the repercussions in terms of operational cost and performance.

In this context, remote sensing (RS) represents a powerful tool,capable of providing open access(Woodcock et al.,2008)and up-to-date(Francini et al.,2020,2022b;Vaglio et al.,2021;Francini et al.,2022a)data at a spatial resolution that is in line with the scale at which human interacts with forest ecosystems (Wulder et al., 2018; Francini et al.,2021). In particular, Landsat has been providing consistent data since 1985 with a frequency of 15 d and at 30-m spatial resolution. Recent changes to the open and free data policy have further revolutionized the use of Landsat data, allowing the creation of robust standard products and science applications (Wulder et al., 2020). In addition,cloud-computing platforms have drastically increased the ability to process these remotely sensed images, allowing the development of cloud-based processing tools,thus simplifying large-scale global projects(Gomes et al., 2020). A key example of a cloud computing platform is Google Earth Engine (GEE), which brings massive computational capabilities to bear on a variety of high-impact societal issues, including deforestation, drought, disaster, disease, food security, water management,climate monitoring,and environmental protection(Gorelick et al.,2017).

RS-based estimates of the quantity and spatial distribution of deadwood and biodiversity indicators(e.g.,microhabitats)at different spatial scales reflect the availability of food resources for saproxylic beetles, as suggested by previous studies evaluating foraging behavior,abundance,and space use models for this group of animals related to deadwood(Alaniz et al., 2021; Bombi et al., 2019; Müller and Brandl, 2009). RS provides objective and time-efficient measures of forest structure,health,and management, helpful in concentrating biodiversity analyses in specific hotspots for monitoring purposes (Alaniz et al., 2021). RS techniques to assess forest biodiversity from landscape to stand scale have increased rapidly (Udali et al., 2021), with a significant focus on vegetation components and the macro-fauna.Fassnacht et al.(2016)and Yu et al. (2015) provide comprehensive reviews of RS applications on tree species classifications and forest structural variable estimations.Several studies have used optical and LiDAR data to model breeding birds' diversity and habitat distribution (Foody, 2003; Culbert et al.,2012;Vierling et al.,2013),while most studies on saproxylic beetles have focused on the distribution of endangered species related to environmental and climatic variables,such as growing stock volume,net primary production, land use, and precipitation (Hermosilla et al., 2015a; Chen et al.,2020;Della Rocca and Milanesi,2020).Just a few studies tested the relation between RS structural metrics and saproxylic abundance,even if good results were reported.For example,Airborne Laser Scanning(ALS),Terrestrial Laser Scanning (TLS), and satellite RADAR data were successfully used to test the relationship between canopy structure and arthropod diversity (Müller and Brandl, 2009; Bae et al., 2019; Knuff et al.,2020).

Although Landsat Time Series (TS) analysis has a high potential for studying and monitoring the distribution of the saproxylic community,the use of remotely sensed data for assessing the species diversity of saproxylic beetles is a new technique.

Here, we aimed to verify the hypothesis that Landsat TS-derived variables can be used as predictors to estimate the abundance of saproxylic red-listed species and, consequently, help the selection of biodiversity hotspots over large areas. By supporting the sample selection,Landsat-derived variables may decrease the effort needed for data acquisition during field sampling,while increasing the number of insects available for the analysis.

2. Materials and methods

2.1. Overview

We selected four monitoring sectors in a Mediterranean-managed beech forest located in the Apennines (Molise, Italy). First, we characterized each sector concerning on-site species richness and alphadiversity and assessed differences between sectors to compare biodiversity patterns within the beech forest (Hsieh et al., 2016). Second, for a 37-year study period (1984–2020), we constructed Landsat Best Available Pixel composites (White et al., 2014) from which an NDVI time series was obtained and used to calculate eight Temporal Metrics(TMs)that summarise the NDVI trend over time. Then, we investigated the relationship - in terms of Pearson's product-moment correlation - between Landsat TMs and the abundance of saproxylic key species at the family and trophic category level.

2.2. Study area

The study was conducted in a beech forest of the Apennines (Italy),the Roccamadolfi forest,which is included within the Site of Community Importance (SCI, http://natura2000.eea.europa.eu) “La Gallinola -Monte Miletto - Monti del Matese” (Cod. IT 7222287) and the National Park of Matese(Fig.1).The experimental area covers about 400 ha,and its altitude ranges from 1,180 to 1,737 m a.s.l.The site is characterized by volcanic soils on calcareous substrates (aluandic and silandic andosols).According to the K¨oppen-Geiger classification, the climate is warmsummer Mediterranean (Peel et al., 2007). The mean annual temperature is 12.2°C,and the mean total annual precipitation is 1,802 mm,with the maximum in autumn (283 mm) and the minimum in summer (52 mm). Forest covers 70% of the total area and is represented by nine different forest types(EEA,2006),the most common being beech forest of about 8,000 ha, mostly coppices resulting from past forest management(Vizzarri et al., 2015;Parisi et al.,2020c).

In this area, the priority habitat 9210 (Apennine beech forests with Taxus and Ilex) is mainly extended, especially at a higher elevation.Associated trees in the beech forests include Taxus baccata (L.), Acer pseudoplatanus (L.), Ilex aquifolium (L.), Sorbus aria (L.) Crantz and S. aucuparia (L.), while the herbaceous layer consists mainly of Daphne laureola (L.), Lathyrus venetus (Mill.) Wohlf., Melica uniflora Retz, Geranium versicolor(L.),Potentilla micrantha Ramond ex DC.

2.3. Experimental design

The study was conducted in the same areas as Proietti et al. (2020).Four monitoring sectors were selected: high altitude and north aspect(HN - Range: 1,395-1,509 m a.s.l), low altitude and north aspect (LN -Range:1,180-1,273 m a.s.l),high altitude and south aspect(HS-Range:1,556–1,737 m a.s.l), and low altitude and south aspect (LS - Range: 1,190–1,220 m a.s.l) (Proietti et al., 2020). Within each sector, 15 permanent circular sample plots (radius = 13 m and area = 530 m2) were established.

The Roccamandolfi forest,regularly managed for timber production,is characterized mainly by mountainous beech forest type and, secondarily, by sub-mountainous beech forest type at lower altitudes (Vizzarri et al.,2015;Parisi et al.,2020c).The forest is predominantly even-aged of gamic and agamic origin.Tree density ranges between 701 and 2,174 trees·ha-1(Proietti et al.,2020).

UTM-WGS84 coordinates(Zone 32T)and elevation were recorded in each plot using a Juno SB Global Positioning System (GPS) (Trimble,Sunnyvale,California).

2.4. Deadwood survey

The sampling protocol followed Lombardi et al.(2015):dead downed trees, snags, coarse woody debris (CWD), and stumps were sampled,measuring their length/height, minimum (≥5 cm), and maximum diameter, recording the species, when possible. The volume of snags,CWD, and stumps was calculated through the cone trunk formula(Lombardi et al.,2012):

where V = volume (m3), h = height or length (m), D = maximum diameter(m),and d=minimum diameter(m).

2.5. Saproxylic beetle trapping

In the same 60 plots where deadwood was surveyed, the capture of saproxylic adult beetles was carried out, using window flight traps for flying beetles. Traps were checked approximately every 30 d for five surveys in 2018(from May to October).All the monitoring systems were then removed during winter.

Systematics and nomenclature followed Bouchard et al. (2011) and Audisio et al.(2015).All the taxa collected during the field activities are alphabetically listed in Additional file: Appendix 1. Species strictly considered as saproxylic(sensu Carpaneto et al.,2015)are also reported in Additional file: Appendix 1, together with their risk category at the Italian level(see Audisio et al.,2015;Carpaneto et al.,2015).The beetle assemblages were grouped according to the prevalent trophic categories,defined by Audisio et al.(2015)and Carpaneto et al.(2015):i) Xylophagous (XY) (organisms feeding exclusively or mainly on wood), ii) Saproxylophagous (SX) (organisms feeding exclusively or largely on fungus-infected wood), iii) Mycophagous (MY) (organisms feeding exclusively or mainly on fungi), iv) Mycetobiontic (MB) (organisms feeding on carpophores of large Polyporales and other fungi living on old trees and stumps), v) Sap-feeder (SF) (sap-feeders on trees attacked by Xylophagous), and vii) Predator (PR) (organisms that primarily obtain food,killing and consuming other organisms).

2.6. Landsat time-series and TMs calculation

Candidate images were selected from the Landsat archive available on GEE (https://developers.google.com/earth-engine/datasets/catalog/landsat). The dataset consists of Landsat surface reflectance imagery atmospherically corrected using LEDAPS (Wolfe et al., 2004) and acquired with a solar zenith angle smaller than 76°. Images contain three visible bands(blue,green,and red),a near-infrared band(nir),and two short-wave infrared bands(swir1,swir2)pre-processed to orthorectified surface reflectance, and one thermal infrared (TIR) band processed to orthorectified brightness temperature. Images include masks informing clouds, shadow, water, and snow produced using CFMASK (Foga et al.,2017).

Selected images are covered by clouds less than 50%and are acquired between Jun-15 and Aug-15 during the 1984–2020 period.Images with cloud cover greater than 50%were removed because they are more prone to geographical location errors due to the challenges of performing geometrical corrections when ground control points are obscured(White et al.,2017).

For each year in the study period(1984–2020),we calculated cloudfree composites of the study area. When more than one cloud-free observation was available for a specific year, the “best” pixel was selected using the Best Available Pixel procedure (BAP), for which a comprehensive description was given by Griffiths et al. (2013), White et al.(2014)and Hermosilla et al. (2015a,b).

Fig. 1. Study sites and distribution of sectors for each study site. LN = Low altitude and North aspect (41.48° N; 14.34° E); HN = High altitude and North aspect(41.56° N; 14.34° E); HS = High altitude and South aspect (41.45° N; 14.35° E); LS = Low altitude and South aspect (41.50° N; 14.29° E).

When any observation was available for a specific year (data gaps),we obtained a valid synthetic observation by linear interpolating the first two valid observations in previous and subsequent years. As a result of this step,we obtained a BAP-collection of 37 cloud-free composites,one for each year between 1984 and 2020.

We calculated the NDVI as the normalized difference of NIR and RED Landsat spectral bands for each image.We extracted the NDVI time series over the plot in our reference dataset by calculating for each image the mean of the pixels included in each plot.When the plot was included in a single Landsat pixel,that pixel was selected.To summarise the NDVI 37-year time-series and to obtain more generalizable variables we calculated a set of eight Temporal Metrics TMs similar to those proposed by Potapov et al. (2020):(i)the slope and(ii) the Pearson correlation coefficient of the regression line obtained by linear interpolating the NDVI values over time;(iii)the RMSE calculated between the regression line and the NDVI time series; (iv) the median, (v) the mean, (vi) the minimum, (vii) the maximum, and (viii) the standard deviation value of the NDVI time series.

2.7. Diversity indices

This study assessed species diversity for each assemblage using a unified family of diversity measures called Hill numbers or effective numbers of species, incorporating relative abundance and richness. Hill numbers are defined, for q/=1 as

where S is the number of species in the assemblage,and piis the relative abundance of ith species,for i=1,2,…,S.The parameter q determines the sensitivity of the measure to the relative frequencies.When q=0,D is simply the species richness, where only species presences are counted without regard to their relative abundances. For q = 1, D is undefined,but its limit as q tends to 1 is the exponential of the familiar Shannon index,referred to here as Shannon diversity(Chao et al.,2014):

for q = 1, D weighs species in proportion to their frequency and can be interpreted as the effective number of common species in the assemblage.

The measure for q=2,referred to as Simpson diversity,discounts all but the dominant species placing more weight on the frequencies of abundant species and discounts rare species and can be interpreted as the effective number of dominant species in the assemblage.

This study reports two of the most widely used members of the Hill numbers family(q=1 and 2).

We also compared species diversity across assemblages, using rarefaction and extrapolation curves,with functions provided by Hsieh et al.(2016)for the two used members of the Hill number.

To characterize the saproxylic, we compared the four monitoring sectors in terms of species diversity, the amount of deadwood (stumps and snags),and TMs.We searched for significant differences between all different pairs of monitoring sectors using the Wilcoxon-Mann-Whitney test.

2.8. Relationship between TMs and saproxylic communities

Using our reference data, we grouped the species according to their family and trophic category to verify possible relationships between the TMs and the saproxylic community. Pearson's product-moment correlation(r) matrix between each TMs and the number of trapped individuals for each family and the trophic category was calculated.Results obtained using data from families and trophic categories with less than five individuals trapped over the whole study area were considered not statistically reliable and excluded from the analysis.

3. Results

3.1. Deadwood abundance

Deadwood volume showed an irregular distribution among the sample plots, with a mean value of 248.1 m3·ha-1and a standard deviation of 638.9 m3·ha-1. Stumps were the most abundant deadwood component(235.5 m3·ha-1on average among the sample plots),followed by CWD(6.8 m3·ha-1on average)and snags(5.7 m3·ha-1).

3.2. Saproxylic communities

We collected 2,334 specimens belonging to 64 species referring to 27 families of saproxylic Coleoptera (see Additional file: Appendix 1). The most abundant species were Ernoporicus fagi(Curculionidae family)with 1,232 specimens and Hemicoelus costatus (Ptinidae family) with 214 specimens, representing 62% of the total sampled beetles. The Curculionidae family represented 63%of the whole beetles captured,followed by Ptinidae(9.8%)and Melyridae (7.7%).

Regarding the IUCN risk categories, the sampled saproxylic beetles were classified as follows: Vulnerable (VU; 2 species), Near Threatened(NT;8 species),Data Deficient(DD;1 species),and Least Concern(LC;53 species).All the collected taxa are alphabetically listed in Additional file:Appendix 1. The proportion of the collected species,grouped according to the prevalent trophic categories,is reported.Xylophagous represented 28%of the total sampled families,followed by Saproxylophagous(25%),Mycophagous(21.8%),Predators(11%),and Mycetobiontic(6.2%).Sapfeeder,only one specimen was collected.

3.3. Diversity indices

The values of diversity indices in the four monitoring sectors are shown in Fig. 2, for the whole saproxylic community and grouped by trophic categories.

The diversity indices analyzed were generally higher for the whole community than for the single trophic categories,except for the Simpson index for the mycophagous and predator.

Regarding the whole community, the diversity indices were always significantly different for the monitoring site pairs HS-LN and LN-LS(Table 1). Considering the Shannon index, no significant differences were found among the monitoring sites for mycophagous and predator,while the Simpson index showed a more significant difference among sites for all the trophic categories(Fig.2).

Fig. 3 showed the rarefaction and extrapolated curves for the two Hill's numbers used in the study.For the Shannon diversity(q =1), the number of equally-common species was higher in HN sites than in all other sites,reaching an observed diversity of 8.68,against the minimum of 5.10 in the HS sites.From the extrapolated curve it can be stated that the sampling effort was adequate to describe the species diversity in all sites. The Simpson diversity (q = 2) showed the same pattern as the Shannon diversity in which the HN sites reached an observed diversity of 5.80, against the minimum of 2.44 in the HS sites. The analysis of the Wilcoxon-Mann-Whitney test revealed that monitoring sites differed primarily in NDVI TMs and the amount of deadwood(stumps and snags).Specifically, only the RMSE and the mean NDVI metrics were significantly different among the six pairwise combinations of monitored sites,followed by the median,slope,and R2,with five combinations out of six(see Additional file: Appendix 2). Regarding the deadwood, the volume of snags had the strongest effect in determining the difference among sites,followed by the stumps and coarse woody debris(Table 2).

Fig. 2. Boxplot of diversity indices for each monitoring sector for the whole saproxylic community and the main trophic categories. LN = Low altitude and North aspect;HN=High altitude and North aspect;HS=High altitude and South aspect;LS=Low altitude and South aspect.XY=Xylophagous;SX=Saproxylophagous;MY = Mycophagous; PR = Predator.

Table 1 Results of the Wilcoxon-Mann-Whitney test and related effect size of the diversity indices.

The other variables had less influence in determining differences among the four monitoring sites.

3.4. Relationship between NDVI TMs and saproxylic communities

Fig.4 shows the distribution of the four most correlated TMs with the species abundance for each monitoring sector.

Based on the relationships between the NDVI TMs and the saproxylic community, the strongest correlation was reached by the Monotomidae family (represented by one species with 12 individuals), followed by Cerambycidae and Curculionidae(eight species with 128 specimens and eight species with 1,470 specimens,respectively),with a mean absolute correlation (r) of 0.46, 0.31, and 0.25, respectively. The strongest correlation was reached between the Monotomidae family and the RMSE metrics, with an R-value of 0.66. TMs with the strongest correlation among all families were the slope(maximum correlation of 0.5,between the Curculionidae family),followed by the minimum of the NDVI,with a maximum correlation of 0.52(between the Monotomidae family).RMSE was the only metric in which all the correlations were positive for all families, while minimum and mean were always negative. Fig. 5 shows the correlation between the abundance of each family of saproxylic insects and the four NDVI TMs that reached the average strongest correlations among families in terms of R.

Fig.3. Rarefaction curves with 95%confidence intervals(shaded areas)for the two Hill's numbers(q=1 Shannon diversity and q=2 Simpson diversity),considering the whole saproxylic species.

Table 2 Results of the Wilcoxon-Mann-Whitney test and related effect size of the deadwood variables.

Fig.6 shows correlation coefficients between the abundance of taxa at the trophic category level and the four NDVI TMs that reached the average strongest correlations among trophic categories. TMs with the strongest correlation among all trophic categories was the slope, with a maximum correlation of 0.5 between the Xylophagous group(28%of the insect species), followed by the minimum, with a maximum correlation of-0.4(between the Mycetophagous group,representing the 25%of the total number of species). The Xylophagous group reached the strongest correlation against the slope,with a value of 0.5.Moreover,Xylophagous and Mycetophagous groups reached the strongest mean absolute correlation among the metrics(0.36 and 0.33,respectively).The only metrics with all positive correlations among all trophic categories was the RMSE,while consistently negative correlations were obtained by minimum and mean metrics. These results are consistent with those obtained for the saproxylic beetle families.

4. Discussion

4.1. Diversity of saproxylic communities

In the investigated area we found as many as ten threatened species included in the IUCN red list(two Vulnerable and eight Near Threatened species) present in the different proportions of all forest sectors (Additional file: Appendix 1). The analysis of alpha diversity of the whole saproxylics community and the trophic categories shows a certain uniformity across the four monitored sectors. The LN sector tends to have greater diversity, considering communities (Fig. 2). This aspect can be attributable to the trophic category of xylophagous, which reaches its maximum expression in LN (Fig. 2), probably due to the conspicuous presence of stumps and insect holes, relatively abundant in this sector compared to the others. Also, in HS, xylophagous species are wellrepresented thanks to the presence of stumps(Parisi et al.,2021).

Fig.4. Boxplot of TMs values distribution for each monitoring sector.LN=Low altitude and North aspect; HN = High altitude and North aspect; HS = High altitude and South aspect; LS = Low altitude and South aspect.

Furthermore,the differences between the four sectors are confirmed by analyzing the number of individuals per trophic category. The rarefaction curve proves that a sufficient number of individuals were sampled to describe the diversity in the various forest parts monitored (Fig. 3).Exposure is a determining factor for beetle diversity in the higher elevation monitoring sectors(HS and HN).While in the lower sectors(LS and LN),the exposure does not seem to affect the diversity of the species,probably because of minor differences in temperature. It is a matter of fact that, in this Mediterranean mountain environment, stronger winds and lower temperatures at higher elevations may result in a thresholdtype response for beetle species in HS and HN, which is less evident in LS and LN.In HS,although a greater diversity than HN can be expected,due to higher temperature, higher irradiance may favor the establishment of ecologically more specialized species, resulting in a lower diversity.

Significant variations in TMs distribution among the monitoring sectors are visible for the slope and minimum, which show the most significant differences(Fig.4).Mean and minimum are metrics of direct photosynthetic activity and tend to have an opposite trend to RMSE and slope in the monitoring sectors. Slope reaches the highest values in the sector with the highest photosynthetic activity (HS). Indeed, HS is characterized by relatively low coverage(701 vs.2,174 trees·ha-1for H and L, respectively) compared to other sectors, and has higher photosynthetic variations over the years than the others.

4.2. Relationship between NDVI TMs and saproxylic communities

The best correlations were obtained with four out of eight TMs(mean,minimum,slope,RMSE).Of these,the mean and the minimum show only negative correlations and the rest only positive, both at the level of families and trophic categories. The first metrics (mean, minimum) are directly related to the degree of forest cover and thus influence the solar radiation that filters into the forest,as well as the soil temperature and its seasonal variation (Horak, 2014; Henneron et al., 2017) (Fig. 5). Radiation exposure determines different ecological niches that differentiate the food substrates and, therefore, the trophic categories of insects.Because some families of beetles need direct solar radiation and, therefore,low tree cover(low NDVI)to carry out their biological functions,an inverse correlation is expected to exist between mean and min TMs and species abundance. Indeed, solar radiation plays a fundamental role in deadwood decomposition processes.On the other hand,RMSE and slope are related to the variability in photosynthetic activity, which may directly or indirectly affect the diversity and abundance of saproxylic species favoring the biological activity of adult beetles and allowing pollination(Thorn et al.,2020).

Regarding the beetle families found in the study area,Monotomidae reached the greatest correlation among three out of the four considered TMs(mean = -0.64,min = -0.52,RMSE= +0.66).The abundance of these beetles depends on light radiation, and, possibly, lower NDVI values indicate that a low canopy cover may favor their biological cycle.For example, the Monotomidae family includes various predatory and rhizophagous species, most living on decaying organic material, under bark or in wood. Several species are found in the galleries of other saproxylic species (Parisi et al., 2018). The beetle family that obtained the lowest correlations in all metrics is the Elateridae, probably due to the great diversity of trophic categories occupied by the species of this group of insects. In the study area, the sampled elaterids are represented by predators of other saproxylics,hence,not directly affected by stationary and climatic conditions. These results are confirmed by the analysis of trophic categories,in which predators show a low correlation in all TMs.

Fig.5. Correlation coefficients(R)between the abundance of each family of saproxylic insects and the four TMs returning the strongest correlations(mean,minimum,RMSE, slope).

Fig. 6. Correlation coefficients (R) between the abundance of species in each trophic category of saproxylic insects and the four TMs returning larger correlations(mean, minimum, RMSE, slope). XY = Xylophagous; SX = Saproxylophagous; MY = Mycophagous; MB = Mycetobiontic; PR = Predator.

The abundance of insects in some families (i.e., Curculionidae,Lucanidae, and Cerambycidae family) is directly related to the slope,representing the trend of photosynthetic activity over the past 37 years.These results are also confirmed by correlations between the abundance of species in each trophic category (Xylophagous (XY) and Saproxylophagous (SX)). The Ernoporicus fagi mainly represents the Curculionidae family in all monitoring sectors,whose presence is strictly related to mature beech forests. Slope can be considered a proxy of forest maturity and management since steep terrain makes it difficult for forest operations.

Among the various ecological parameters that influence the diversity and abundance of species, the available solar radiation represents an element of discrimination in beech forests(Salmone et al.,2008).Dense cover promotes, in particular, the diversity of saproxylophagous and xylophagous,while secondarily supporting low abundances(Fahrig and Merriam, 1994). Differences in solar radiation determine microclimatic conditions that allow a more diversified community of beetles (Parisi et al.,2020b;Thorn et al.,2020).For example,Horak(2014)argued that high species richness is determined by an open structure, recalling the conditions of solitary trees. Furthermore, Henneron et al. (2017) found that high canopy coverage, typical of even-aged beech forests, reduces predator diversity.

Good solar radiation is, for example, essential for several species of Cerambycidae and Elateridae.In addition to deadwood,the development of various saproxylic larvae also requires other sources of nourishment,such as pollen and herbaceous plants,whose presence is greater in forests with canopy gaps or forest margins.These stand structures and substrates need to be available and well distributed within the forest(Horak,2014).Moreover,the Monotomidae family shows a significant correlation with the RMSE metrics.Indeed,the RMSE metric can be considered a measure of the global variation occurring in the time series.

Regarding trophic categories, the four NDVI TMs follow the same pattern of beetle families (Fig. 6). The Xylophagous shows the best correlations in all sectors, while the Mycetobiontic has the worst correlations, probably due to the low presence of wood-inhabiting fungi attributable to the management practices in this forest stand with relatively young trees (Proietti et al., 2020). Several authors have observed that wood-inhabiting fungi are more abundant in unmanaged or old-growth forests than in young managed forests (Siitonen, 2001; De Zan et al.,2014).

4.3. Opportunities and limitations of RS in monitoring saproxylic communities

These observations may allow the design of conservation strategies tailored for threatened saproxylic beetles and the development of management practices to address the ecological needs of beetle communities and overall biodiversity decline. In particular, as already suggested in other studies (Koivula, 2011; Lange et al., 2014), management options aimed at preserving saproxylic beetles should promote larger trees,favor habitat trees, and reduce forest density, in short diversifying the forest structure. These characteristics are typical of mature forests, the achievement of which should be a major objective of sustainable forest management in the studied site, which is under conservation management(i.e.,Natura,2000).In practice,a mature/aging forest structure can be achieved by selecting several dominant trees,eliminating competitors,and reducing overall forest density (Timonen et al., 2010). This is partially in contrast with the forest management approach currently adopted in the present Roccamandolfi forest, as well as in many Mediterranean mountain beech forests. Structural variables are not the only habitat drivers for animal species,which are also determined by multiple interspecific and intraspecific biotic interactions (Parisi et al., 2018).Incorporating such interactions into modeling and prioritizing approaches is an important challenge(e.g.,Araujo and Guisan,2006;Wisz et al., 2013). Therefore, our integrated approach may allow identifying hotspots in the forest with the greatest need for management interventions aimed at improving the conservation of saproxylic and non-saproxylic red-list species.

The ability to identify conservation priorities, combined with the resolution of Landsat data,is linked to the use of appropriate predictors,as well as the availability of detailed ground data (Bombi et al., 2019;Campanaro and Parisi,2021;Parisi et al.,2021).In the present case,we paid particular attention to the design of a sampling protocol that requires a relatively high field monitoring effort, obtaining an exhaustive analysis of the saproxylic beetle communities (see Additional file1: Appendix 1). Integrating a high-resolution modeling approach of various biodiversity indicators (e.g., deadwood and microhabitat) allows identifying general and local management priorities. On the one hand,large-scale analyses make it possible to overcome local errors (Groves,2003).On the other hand,local analyses may,at the same time,allow the planning of conservation initiatives for practical implementation(Bombi et al., 2019; Huber et al., 2010). Although there are many knowledge gaps,especially on the response of beetles to forest succession in natural conditions and to anthropogenic disturbances in managed stands, saproxylic beetles are considered local proxies for the sustainability of forest management (Pearce and Venier, 2006). Our analyses show that NDVI TMs are related to the abundance of red-listed saproxylic beetles in these Mediterranean mountain forests.

While the fieldwork in this study involved three people for six months(excluding the long phase of preparation and taxonomic determination of the species), the processing of the satellite images was done by one person in a few days.This study highlights the potential of applying RS tools in conservation planning and decision-making in the forest landscape. These tools can be used at different spatial scales, are highly repeatable,and are helpful for monitoring purposes,as data can be easily compared over time.The recent advancement in the availability of highresolution satellite images and global Landsat images now freely available on the web (NASA Geocover dataset; Tucker et al., 2004) allows estimating productivity using vegetation indices and contemporaneously examining the relationships between these estimates and biodiversity models (Turner et al., 2003). Despite the great opportunities offered by RS and,in particular,the long Landsat time-series data,some limitations have to be accounted for.First,the spatial resolution of Landsat could not be adequate to capture micro-spatial variations in the distribution of saproxylic species, which have poor dispersion capacity, making the Landsat sensor suitable only for large-scale monitoring. Second, due to the well-known saturation effect of multispectral data,the Landsat sensor could not be sensitive to multilayer canopy forests, dense forests or complex topographic features (Chirici et al., 2020; Vangi et al., 2021;D'Amico et al., 2022), affecting NDVI values. Authors were aware that satellite data, by providing canopy-level information, cannot fully explain the biodiversity richness of the ecosystem. However, the integration of RS approaches and on-site monitoring activities in forest management guidelines may help design and optimize conservation strategies for this taxon, tuning forestry practices to meet habitat preferences and spatial requirements of threatened beetle species.

Although RS cannot replace fieldwork or identify individual insect species, rarity, and composition, we suggest that, given the valuable products of such images, investment in processing and analysing these data will be highly affordable.In this regard,the availability of the GEE cloud platform allows an unprecedented view of forest areas worldwide.

5. Conclusions

In this study, we demonstrate that Landsat time series are effective variables for predicting the presence and abundance of saproxylic insects.Although referred to a relatively small dataset,we suggest that saproxylic trophic categories in the present Mediterranean mountain forest ecosystem can be predicted with R up to 0.5, while saproxylic families with R up to 0.7.We believe that RS data cannot replace on-site sampling and analysis, as they are needed for a reliable assessment of population dynamics and activity patterns, identifying community structure and function, rare and threatened species,and understanding the ecological factors on which targeted interventions need to be planned. Indeed,suitable RS imagery can hardly be available,currently,with the desired resolution in mountain environments, which are characterized by persistent cloud cover and extensive snow mantle.On the other hand,the herein proposed BAP-derived TMs can be considered a valuable tool to identify insects’ hotspots – or at least areas where very few insects are expected to be.In this sense,advances in remote sensing technology will probably open new avenues to account for temporal and spatial changes in species distribution. These observations have important implications for studying species distribution and may help the selection of monitoring areas for concentrating field analysis, therefore, increasing the amount of information acquired and decreasing the level of effort required.Accordingly,wall-to-wall maps of the proposed metrics need to be further investigated in future research as a valid indication of areas where more saproxylic insects are expected.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

All data generated or analyzed during this study are included in this published article.

Competing interests

The authors declare that they have no competing interests.

Funding

This work received no external funding.

Authors’contributions

Conceived and designed the experiment:Francesco Parisi,Elia Vangi,Saverio Francini; data and samples in the field: Francesco Parisi. Processed samples in the lab:Elia Vangi and Saverio Francini.Analyzed the data and wrote the manuscript: Francesco Parisi, Elia Vangi, Saverio Francini.All authors read and approved the final manuscript.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors thank the specialists of the various taxonomic groups:Paolo Audisio (Nitidulidae), Alessandro Bruno Biscaccianti (Biphylidae,Cerylonidae, Erotylidae, Lucanidae, Mycetophagidae, Salpingidae, Tenebrionidae, Trogossitidae), Enzo Colonnelli (Curculionidae), Massimo Faccoli(Curculionidae Scolytinae), Fabrizio Fanti(Lycidae),Gianfranco Liberti (Melyridae), Emanuele Piattella (Scarabaeidae), Giuseppe Platia(Elateridae), Pierpaolo Rapuzzi (Cerambycidae (pars)), Adriano Zanetti(Staphylinidae (pars)). We are grateful to Bruno Lasserre (Universit`a degli Studi del Molise, Italy)for technical support.

Abbreviations

TS Time Series

RS Remote Sensing

GEE Google Earth Engine

ALS Airborne Laser Scanning

TLS Terrestrial Laser Scanning

HN High altitude and North aspect

HS High altitude and South aspect

LN Low altitude and North aspect

LS Low altitude and South aspect

XY Xylophagous

SX Saproxylophagous

MY Mycophagous

MB Mycetobiontic

SF Sap-feeder

PR Predator

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://do i.org/10.1016/j.fecs.2022.100050.