Assessment of Two Adjacent Mountainous Riparian Areas along Nestos River Tributaries of Greece

: The riparian areas of the Mediterranean are unique but face many pressures from anthropogenic and climate change impacts. They have very rich and diverse plant communities as a result of the dynamic ﬂuvio-geomorphologic conditions of the Mediterranean streams. In this study, the riparian areas of two adjacent mountainous tributaries (Arkoudorema and Vathirema) of the Nestos River were studied. To assess the condition of riparian areas holistically, diverse measurements are required. This is why ﬂuvio-geomorphologic (in the ﬁeld and with GIS), vegetation (surveys and visual protocols) and ground-dwelling insect (pitfall traps and indices) measurements were taken along an elevational gradient. The results of all three methodologies draw to similar conclusions, with Vathirema sub-watershed riparian areas being in better condition than Arkoudoreama. This was expected, since Vathirema has less anthropogenic pressures. In addition, the riparian areas in higher elevations were in better condition for the same reason. To implement integrated water resources management plans, ﬂuvio-geomorphologic and biological (e.g., vegetation and insects) datasets are required to provide a holistic view on the watershed and riparian area conditions. For the studied sub-watersheds, we recommend these measurements to continue, to record the current anthropogenic pressures and based on this information to suggest best management practices that will secure long-term sustainability.


Introduction
Mediterranean riparian areas have been inhabited for thousands of years [1]. In many cases, these areas were the center of great civilizations. This long-term interaction of riparian areas with humans in the Mediterranean is the reason why anthropogenic disturbances have been as important as natural factors in structuring their plant communities [1]. Most riparian areas are heavily degraded, especially in lowlands, due to river flow regulation, channel alterations, groundwater withdrawal, agricultural activities and urbanization [2][3][4][5]. Few natural riparian ecosystems remain in the Mediterranean, primarily located in mountainous areas with few inhabitants [6].
Riparian areas of the Mediterranean are of greater importance compared to those of mesic regions [7]. They are distinctive because of the region's climatic, topographic and hydrologic unique conditions [8] along with the frequent occurrence of natural wildfires [9]. The rivers and streams have natural high stream flow variability [10]. In the high elevation areas (precipitation > 1000 mm year −1 ), the streams have a perennial flow that peaks typically in spring after major rainfall events or snowmelt [11]. In the lowland areas, precipitation is substantially less (200-500 mm year −1 ), leading to many intermittent or ephemeral streams [11]. These stream types have a highly irregular hydrologic regime with a great amount of water and sediment transported at high velocities in short periods of time, leading to frequent flash floods [12]. The fluvio-geomorphologic conditions of Mediterranean streams influence riparian diversity and functionality, and understanding these conditions is essential to sustainably manage riparian areas [13,14].
The dynamic fluvio-geomorphologic conditions of Mediterranean streams lead to high year to year and interannual fluctuations in riparian plant species numbers and composition that can vary from woody, to shrubby, to herbaceous vegetation and mixtures [15]. This diversity is also spatial, as they have distinguishable gradients because of the more discrete changes in the groundwater level and soil moisture, and different degrees of natural disturbance impact on the riparian areas as you move away from the aquatic ecosystem. Overall, riparian vegetation assemblages experience high colonization rates and desiccation resistance mechanisms, while their life expectancies are relatively short [16,17].
The protection and conservation of riparian areas has become a priority worldwide due their severe degradation [8]. Riparian areas are ecotones (transitions zones) with biophysical gradients as you move from the terrestrial to the aquatic ecosystem [18]. These gradients lead to the many ecosystem services of riparian areas that have and continue to be utilized, mostly unsustainably, by humans.
The increase in temperatures due to climate change has also increased evapotranspiration in rivers in the summer and decreased snow accumulation in the mountains during winter [19][20][21]. This alters the flows of snow-fed mountainous streams of the Mediterranean region. The riverine ecosystem annual stream flow is expected to be reduced, while water temperature and large flood magnitudes are expected to increase and directly impact riparian areas [22,23].
Most Mediterranean riparian areas will be adjacent to intermittent and ephemeral streams due to climate change's impact [8]. These riparian areas have characteristics different to those adjacent to perennial streams [24][25][26]. Native riparian species are adapted to seasonal peak flows, and this alteration will retard their regeneration, lessen their growth and allow other species to be more competitive and move into the riparian areas (typically invasive species). Finally, longer and more intense droughts might exceed the natural resilience boundaries of native riparian flora and lead to their demise [27,28].
Assessing the condition of riparian areas requires an interdisciplinary approach that measures abiotic and biotic characteristics [28]. Such an approach will evaluate more accurately the impacts of human actions and climate change on the ecological functioning of riparian areas [29]. This is particularly true for the complex Mediterranean riparian areas that require hydro-, geo-and bio-diversity monitoring to suggest best management practices for functional riparian ecosystems [13,30]. In addition, monitoring methods need to be easy-to-use, reliable, inexpensive and efficient [28].
The aim of this study was to assess and compare the condition of riparian areas of two adjacent sub-watersheds of the Nestos River. These two watersheds have different characteristics and anthropogenic pressures but are representative of Mediterranean mountainous sub-watersheds. Many natural mountainous watersheds and, particularly, their riparian areas, including those of this study, have had increased visits due to the mainstreaming of ecotourism [28]. Specifically, hydrological, geomorphological, vegetation and insect parameters were estimated and measured. These parameters were measured with channel characteristics, vegetation surveys, visual protocols and insect collection and estimated with GIS and biodiversity indices. Spatial comparisons (based on elevation) for each watershed were also conducted. This is one of the few studies that has measured and compared the hydrologic, geomorphologic and biodiversity characteristics of mountainous riparian areas of Greece to such an extent.

Study Area
Arkoudorema (ADR) and Vathirema (VTR) were the two sub-watersheds of the Nestos Basin that were assessed ( Figure 1). Both are within the Rodopi Mountain Range National Park in Northern Greece (Figure 1). The annual precipitation ranges from 600 mm year −1 in the lowlands to 1500 mm year −1 in the mountainous areas [31]. The climate in the lowlands is characterized as Mediterranean, with relatively mild winters and dry, hot summers. It must be noted that these lowlands receive greater amounts of rainfall compared to typical Mediterranean areas. In the mountainous areas, the climate is a transition between Mediterranean and mid-continental climate. The ADR sub-watershed is primarily mountainous with 6 villages within its boundaries and Paranesti (the largest town in the region) very close. The VTR sub-watershed is one of the most sparsely-populated regions in Greece with 1 village with its boundaries and 4 villages nearby. The road network in the ADR is substantially greater than in VTR. This along with the more inhabitants in ADR than VTR is the reason for more anthropogenic activities. These activities include livestock grazing by local residents along with fishing and hunting. In the last decade, an increase in ecotourism has also been experienced in ADR because of its proximity to Paranesti and more extensive road network. It must be noted that as the elevation of both sub-watersheds increased the distance from villages and town increases, the network become sparser leading to more limited anthropogenic activities. Overall, both watersheds are very scenic and belong to the Natura 2000 Network.
Within each watershed seven sampling stations were established. Emphasis was given to the homogeneous distribution of the stations to cover all elevation ranges, vegetation zones and geomorphologic characteristics. The selection criteria were: (a) complete coverage of the river route from its springs to its outflows, (b) accessibility by the forest road network and (c) recording of all vegetation types. The distribution of sampling stations is depicted in Figure 1b,c. Table 1 indicates the type of measurements taken in each sampling location. We must mote that all type of measurements were not taken in all sampling stations. Table 1. The type of measurements taken at each sampling station of the two studied sub-watersheds. The elevation at each sampling stations is also provided, and an X indicates the type of measurements taken at each station.

Fluvio-Geomorphologic Measurements
Fluvial geomorphology studies the natural process that shape streams and adjacent landforms to understand how streams interact with the surrounding riparian and terrestrial areas [32]. This study concentrated on parameters that impact stream stability.

Watershed Scale-GIS
The following frequently used fluvio-geomorphologic parameters were estimated for each sub-watershed: area, perimeter, maximum elevation, minimum elevation, mean elevation, maximum slope, minimum slope, mean slope, relief, length, shape and roundness degree; for the stream: total length, drainage density, Shreve stream order, Strahler stream order, ruggedness, relief ratio and relative stream power [33]. All the aforementioned parameters provide insights on hydrologic and erosion processes, especially on the ability of a watershed to transport water and sediment [34]. For example, when the watershed shape is elongated, water transport capacity is slower but stream flow has a more constant response to precipitation events [35]. In contrast, circular watersheds have greater peak flows (less water concentration time to the stream channel). Another example is the maximum and minimum watershed elevation that can provide information on the source of the discharge, e.g., if it is snow-fed or not. The stream flow regime is very different in snow-fed watersheds that have higher discharges typically in spring when the snow melts [11]. Using regional maps from the Geographic Service of the Hellenic Army (scale 1:50,000) that were inserted in GIS, the contour lines, watershed boundaries and stream network were digited, and the abovementioned parameters were estimated.

Field Measurements
Four field measurements were conducted: (a) pebble count, (b) water surface width (WSW), (c) bankfull width (BW) and stream bank slope [36,37]. Stream reaches of the sampling stations were surveyed during summer baseflow conditions. Pebble counts are used to estimate the substrate particle size distribution. Stream bed material can indirectly be indicative of streams functions [36]. During the pebble count, the observer walks across the stream in a zigzag pattern within the stream channel and measures stream bed material [36]. The material was measured (mm) with a gravelometer at the intermediate axis (width). At least 100 particles were measured and classified to size classes according to Wentworth [36]. The bed particle that touched the tip of the observer's boot as he walks the stream was measured. It must be noted that riffles and pools were sampled at the same proportion.
The edge of the stream water provides the WSW. To measure BW, indicators such as changes in vegetation and the slope of the bank (slope breaks), staining lines on rocks and boulders, deposits of twigs, leaves, pine needles or garbage and the height of depositional features were utilized [33,36]. Both WSW and BW vary with stream size order and elevation [34]. Finally, the stream bank slopes measured were categorized. Specifically, as: 1 when the slope < 20 • , 2 when the slope was 20 • -45 • , 3 when the slope was 46 • -75 • and 4 when the slope > 75 • [37].

Riparian Vegetation Identification and Evaluation
This was accomplished by surveying the riparian species but also by utilizing visual protocols. Visual protocols are frequently used to assess riparian areas, as they provide accurate and rapid results [28]. Many protocols have been developed and its selection depends on the aim of the study and the experience of the assessor [28]. In this study, we selected the riparian forest quality (QBR) protocol and riparian forest evaluation (RFV) protocol that have been develop for the Mediterranean region.

Vegetation Species Survey
The vegetation of the riparian and terrestrial zone was recorded along the seven sampling stations of both sub-watersheds. Specifically, 50 m upstream and downstream from the sampling station, the dominant and subdominant tree species and shrubs within the riparian zone were recorded. Dominant trees are the largest that form the stand main canopy. Subdominant are smaller trees found close to dominant ones, with smaller crown size. The riparian zone was determined visually and its width differed from station to station. Additionally, the dominant species from the terrestrial zones were recorded. The terrestrial zone surveyed was along the same 100 m transect of the riparian zone. Specifically, a 15 m wide zone from the edge of the riparian zone of both stream banks was surveyed. The identification of the species was based on national (Greek) flora guides [38][39][40]. Based on these tree and shrub species, the habitat types were also determined [41].

Riparian Forest Quality (QBR) Protocol
The QBR evaluates riparian forest quality based on four criteria [42]: (a) Total riparian vegetation cover (TRC), which assesses the riparian areas interconnectivity with the adjacent terrestrial ecosystems. (b) Cover structure (CST), which assesses tree, shrub and understory plant cover in the riparian area. (c) Quality of the cover (CQU), which assesses the riparian habitat and the native tree species based on the geomorphologic characteristics of the site. There are three geomorphologic categories.

Riparian Forest Evaluation (RFV) Protocol
The RFV protocol [37] complies with the objectives of the WFD. It has great applicability at a broad range of riverine areas, particularly in the Mediterranean [37]. The RFV assess the spatial and temporal continuity of the riparian forest. Spatial connectivity is assessed in its three dimensions, longitudinal, transversal and vertical. The temporal connectivity is assessed by the regeneration capacity. Initially, BW and riparian forest width need to be defined, before implementing the RFV.
Longitudinal connectivity is assessed on both banks along a reach that is 10-14 times the BW. It is determined by the autochthonous tree and shrub cover. Bank stretches with rocky substrates are not considered. Transversal connectivity is assessed along 5-7 channel sections. The length of these sections needs to be equal to the width of the riparian zone. This connectivity is determined by the presence of autochthonous tree, shrub and macrophyte species. In the transversal cross-section, vertical connectivity is also assessed in regard to the tree and understory shrub cover, along with the existence of lianoid and epiphytic species. Non-native species when present indicate severe forest quality degradation. Finally, regeneration capacity is assessed along the reach of the longitudinal connectivity. Emphasis is given to the presence of young trees or samplings. Based on these four parameters each riparian stand is classified as: (a) very good, (b) good, (c) moderate, (d) poor and (e) bad.

Visual Protocols Statistical Analysis
A one-way analysis of variance (ANOVA) was conducted for the QBR and RFV protocols. The level of significance was at 5%. The total and each criteria value of the two protocols were compared between the sub-watersheds The two protocols were also compared among sampling stations with different elevation (stations below and above 450 m) and habitat types.

Pitfall Traps
Pitfall traps were used to collect ground-dwelling beetles [43,44]. Studies have shown that these beetles are sensitive to anthropogenic disturbances in riparian areas [16]. The traps were installed in the first week of August 2014. Insects were collected from the traps every 2 weeks in the two studied sub-watersheds until October 15th. The traps were left in situ for 10 weeks. Three sampling stations were chosen in each watershed in areas representative of the riparian forest structure in different elevations (see Table 1). At each station, pitfall traps were installed in three zones, the bankfull, the riparian and the terrestrial. Special effort was given to place traps away from open spaces or paths to avoid edge effects [45]. Pitfall traps consisted of a plastic cup, 10 cm in diameter and 12 cm in depth. They were filled halfway with one part of propylene glycol and three parts of water as a preservative until collection. After the collection, ground insects were conserved in 70% ethanol and were identified to species level according to taxonomic keys in the laboratory [46][47][48].

Population and Diversity Indices
Species composition for each station and zone was estimated based on species identification and the relative abundance. In addition, two different diversity indices were employed. The Simpson index [49] focuses on species dominance, with values from 0 (highest diversity)-1 (lowest diversity) and was estimated for all sampling stations. The Shannon-Wiener index [50] is based on species richness and evenness of the community. When its value is closer to 0, abundance is concentrated on one species. The Shannon-Wiener index was estimated for all stations, each elevation and each microhabitat.

Fluvio-Geomorphologic Parameters
VTR has an elongated shape while ADR is more circular (Figure 1andFigure 2). As ADR also had a larger watershed area than VTR (approximately by 80 km 2 ) ( Table 2), it is expected to have greater discharge. In addition, the ADR had greater values of drainage density, relative stream power and ruggedness and higher stream order than the VTR ( Table 2). The relief ratio [48] values of the two watersheds were almost the same.  The pebble count particle size cumulative curves are presented in Figure 3a. Four sampling stations of the ADR (3, 4, 5, 6) and four sampling sites of the VTR (2, 3, 5, 6), cumulatively, were compared. The ADR had a greater proportion of coarser (large size) particle material. This is because of the reduced stream power and the decreased slope gradient in the VTR. In addition, the variation of the cumulative curves of the different sampling stations is greater in the ADR that in the VTR (Figure 4b,c).  The WSW decreased as the elevation of the sampling stations increased ( Table 3). The ADR WSW are steadily higher than VTR at similar elevations. These findings indicate that the ADR has probably greater discharge magnitudes. BW in most cases had small values with only one station having double digits ( Table 3). The stream bank slopes were categorized as either 1, 2 or 3 (Table 3), although most ranged between 46 and 75 • . * Not possible to measure due to topography or stream flow restrictions.

Vegetation Identification and Evaluation
Along the ADR stations, elevations ranged from 140-1400 m, while along VTR stations, they ranged from 400 m to 1350 m, so differences in vegetation species and habitat types were expected. The hydrophilic Alnus glutinosa L. was the most dominant species prevailing in four stations in both watersheds (Table 4). Corylus avellana L. was dominant at two stations in VTR and subdominant at three stations in ADR. Salix species and Carpinus betulus L. were also present in many stations. The invasive species Populus × canadensis and Robinia pseudoacacia were recorded at one station at each sub-watershed. At the very high elevation zones, the riparian area was not distinguished with the same species present at the riparian and terrestrial areas. At these high elevations, Fagus sylvatica was the dominant species at one station in each sub-watershed and cold conifers such as Pinus sylvestris L., Picea abies (L.) H.Karst. and Abies borisii-regis Mattf. at two stations at VTR and one station at ADR. Salix and Rubus genera were the most common shrubs. The most common habitat type was Alno-Pandion, Alnion incanae and Salicion Albae found in eight stations of both watersheds from 400 to 900 m. The habitat type Salix alba and Populus alba galleries was found at one station at the lowest elevation of the ADR (below 150 m). The two non-riparian forest habitats were recorded at the high elevations (above 900 m). Three stations had the Asperulo-Fagetum habitat type and two had the Rhodopide and Balkan Range Scots pine forests habitat type.
Based on the QBR (Figure 4), there were four riparian quality categories at the ADR locations, with most being good (three stations). For the VTR, the riparian quality categories were only three with most being natural (three stations). The RFV ( Figure 5) had four quality categories with good and very good the most frequent (two stations). The VTR sampling stations had only good and very good conditions (each three stations). No significant difference was found between the QBR mean values of ADR and VTR. The TRC and CST values of the QBR were significantly higher in VTR than the ADR (p-value 0.034 and 0.0012, respectively). No significant differences were found in the RFV values between the ADR and VTR. The stations above 450 m had significantly higher quality values compared to stations below 450 m. The differences were significant for both protocols (QBR: p-value 0.002; RFV: p-value 0.024). Based on the forest habitat types, stations with Alno-Pandion, Alnion incanae and Salicion Albae presented lower riparian quality values compared to stations dominated by Asperulo-Fagetum beech forests and Rhodopide and Balkan Range Scots pine forests, but only for QBR (p-value 0.034).

Ground-Dwelling Insects
In total, 139 individuals were collected and identified to 20 species of ground-dwelling beetles (Coleoptera). Individuals were not equally distributed in each sub-watershed, with most collected in VTR (114 individuals). Similarly, species richness was much higher in VTR (16) than in ARD (7) (Figure 6a,b). At all three elevations, the Shannon-Wiener index was the highest at VTR (Figure 7), with the differences increasing with elevation (greatest difference between ADR 6 vs. VTR 6).
For all traps, both diversity indices showed greater diversity in VTR compared to ADR (Figure 8a,b). The VTR Shannon-Wiener index values also outperformed ADR in each micro-habitat (Figure 9). The bankfull micro-habitat exhibited the lowest Shannon-Wiener index value at ADR, and the other two micro-habitats had slightly higher values. For VTR, the riparian micro-habitat had the lowest Shannon-Wiener index value. Overall, the VTR Shannon-Wiener index microhabitat values were 2-to 3-fold higher compared to the respective values of the ADR.

Discussion
Effective integrated water resources management (IWRM) plans require the measurements or estimation of fluvio-geomorphologic parameters [51,52]. The differences in fluvio-geomorphologic parameters can be quite distinct even between adjacent watersheds. The most distinguishable difference in this study was the shape, since ADR tends to be round, while the VTR is elongated. Watershed shape can impact stream flows, especially peak flows when the other parameters are relatively similar [53]. ADR had larger values than VTR for most fluvio-geomorphologic parameters, such as watershed area and slope, drainage density, stream order, ruggedness and relative stream power. Taking in consideration these parameters and the watershed shape, ADR appears to have greater potential for erosion and flooding events [54]. The use of GIS to determine fluvio-geomorphologic parameters can provide the basic dataset for IWRM plans.
Stream bed particle size can be an indirect indicator of stream velocity, sediment transport capacity and disturbances. Stream beds that exhibit a high percentage of fine material indicate sedimentation that is associated with stream or watershed scale disturbances that cause surface erosion [55]. In Mediterranean streams, this can occur after extensive wildfires [56] that leave the watershed bare of vegetation and susceptible to surface erosion and runoff whose sediments are deposited on the stream bed as fine particles. Both watersheds contained different size particles in their stream bed including larger diameter. This was expected, since both streams are in a mountainous region with relative steep slopes and high elevational differences (ADR is 1630 m and VTR is 1490 m). The ADR had greater stream water power based on the GIS fluvio-geomorphologic parameters, something that was confirmed by the cumulative percentage particle size figures (larger particles on the stream bed). VTR streams on average had a greater percentage of finer material (≤2 mm) (~20%), while ADR had less (≤5%). The ADR cumulative particle sizes at the lower elevation stations as expected contained higher percentages of smaller diameter particles [57,58]. This was evident between ADR1 and ADR 6, but less apparent among the other four stations. This lack of differences (ADR 2-5) is likely due to localized lithological controls [59]. The differentiation of the cumulative particle size at the sampling stations along VTR is not as clear. VTR 2 had higher percentages of finer material but the cumulative lines of all stations are close to each other, especially after the medium size particles. The finer particles percentage for VTR 2 is relatively high, reaching 30%, and is always greater than 15% in all VTR stations. In contrast in all ADR stations, the percentage is always below 10%. High sedimentation can negatively affect invertebrates and fishes [60] and destabilize stream channels [61], making them susceptible to re-suspension [62] but was not a serious problem for these watersheds. In contrast, the very low percentages of finer material in all ADR stations should be a concern because it decreases water retention capacity of streams [57].
Moving downstream should normally lead to an increase in WSW and BW, as larger order streams have greater discharges with deviations due to lithologic or topographic restrictions [63]. ADR sampling stations had higher WSW values than VTR stations at similar elevations. This is a response to greater discharge magnitude events that are expected at ADR. With most sampling stations of both watersheds in semi-mountainous or mountainous areas, topographic restrictions (ravines) were present. Most banks were categorized as 3, which was something expected, since both tributaries have high stream power and steep stream slopes. The lack of bank slopes greater than 75 • is encouraging and indicates that both watersheds do not have serious stream bank erosion [64]. Overall, field measurements and GIS estimations of fluvio-geomorphologic parameters are essential, since they help understand the fluvial processes of the watershed and stream that strongly influence the health in riverine and riparian ecosystems, especially in Mediterranean environments [59].
Riparian areas are azonal [8] and can encompass a great variety of plant species, even dominant and subdominant ones. Vegetation species and habitat types changed along the elevational gradient in both sub-watersheds [18]. This is expected and due to climatic changes, that occur primarily in temperature and precipitation. In this study, the elevation difference between the highest and lowest station was greater than 1000 m in ADR and almost 1000 in VTR. ADR had one more habitat type (Salix alba and Populus alba galleries) that dominates Mediterranean riparian areas below 300 m [41]. As expected, in the elevation from 150 until 900 m due to climatological condition, Alnus glutinosa was the most abundant species, and Alno-pandion, Alnion incanae and Salicion Albae the most dominant habitat. As the temperature became colder (above 900 m), non-riparian species and habitats dominated. These two sub-watersheds showcase that the riparian areas of the Mediterranean Basin exhibit high species richness along elevation gradients that signifies the importance of finding reference riparian sites for restoration efforts [35].
The visual protocol results indicated that the riparian areas of both sub-watersheds are in good condition compared to many other studies in Greece and the Mediterranean [65][66][67]. This was expected since both sub-watersheds have limited anthropogenic impacts, especially agricultural activities. Other studies, also using visual protocols, have found that anthropogenic activities can degrade riparian areas in the Mediterranean [24,[68][69][70][71]. Still, both protocols used in this study indicated that VTR riparian areas are in better condition compared to ADR, even though difference was not statistically significant. The QBR protocol, found that ADR sampling stations showed significantly lower interconnectivity with the adjacent terrestrial ecosystems and lower coverage of trees, shrubs and understory plants in their riparian areas. The ADR sampling station was expected to have a lower score because of the greater anthropogenic activities. For both protocols, the elevation and habitat type of the sampling stations strongly impacted the assessment results. Sampling stations at an elevation higher than 450 m were in better condition. Stations with Asperulo-Fagetum beech forests and Rhodopide and Balkan Range Scots pine forests habitats were of better quality than stations with Alno-Pandion, Alnion incanae and Salicion Albae. In both cases, the key characteristics were the substantially less anthropogenic activities and fewer roads in higher elevation and non-riparian habitat sampling stations. Other studies have also found that riparian areas of lowland regions of Greece have substantially higher anthropogenic activities, particularly agricultural activities [26,69] that degrade their quality. In both ADR and VTR sub-watersheds, anthropogenic activities were more intense near or along riparian areas at stations below 450 m [69] because they are nearer to few villages and towns of the region and have a denser road network.
The results of the two protocols were further confirmed by the evaluation of the ground-dwelling beetle assemblages. More studies are combining the use of visual protocols and insect as bioindicators with good results when assessing riparian areas [70] Ground-dwelling beetles are sensitive to anthropogenic activities and changes in fluvial conditions and vegetation structure [15,[71][72][73][74]. Species richness and abundance and diversity indices of the ground beetles along VTR was considerably higher compared to ADR. ADR exhibited an uneven distribution of species, with Quedius fuliginosus accounting for almost half of the individuals trapped (44%), whereas in VTR, there were five species with a percentage higher than 10%. As anthropogenic activities increase (e.g., grazing, hunting, fishing, ecotourism), these impact the number and the diversity of beetle species. Riparian areas exhibit higher biodiversity compared to the terrestrial and bankfull areas [8]. This was reflected in the Shannon-Wiener index among these habitats in the ADR. In the VTR, terrestrial habitats exhibited the highest Shannon-Wiener index value. This might be attributed to the limited anthropogenic impacts [75] and fewer flood flows [73].
All methods seem very promising, since they were able to detect differences among the two sub-watersheds. The methods provide both a physical and biological assessment of the anthropogenic pressures and changes to riparian quality and will allow to develop a holistic framework for their sustainable management [24]. The greater rural population and more extensive road network in the ADR leads to more pressures such as grazing, hunting, fishing and other ecotouristic activities that has impacted the riparian areas, especially those that are easily accessible in the lower elevations. This is a little worrisome, since both sub-watersheds do not have extensive anthropogenic pressures. The potential increase in ecotourism should be monitored along with the other activities that take place to avoid the degradation of the riparian areas [65]. New management plans should be developed that should not be centrally implemented and be based on public participation [76]. Despite the implementation of the WFD, the current plans do not seem very effective. Alternative plans should be adopted based on ecosystem-based approaches and nature-based solutions [55,76] for the sustainable protection and conservation of the riparian areas, always developed at the basin scale [72,77].

Conclusions
Natural areas are threatened by anthropogenic disturbances and degradation, and for their sustainable management, effective assessment methods are a necessity. For mountainous natural watersheds and riparian areas, this assessment is complicated, since fluviogeomorphologic conditions interact with biological conditions and impact their natural processes. Thus, to implement IWRM plans, a holistic view is required to secure long-term sustainability. The fluvio-geomorphologic, vegetative and entomofaunistic datasets of this study provided different assessment perspectives but had the same conclusions. The leastdisturbed sub-watershed exhibited the better conditions. The fact that the methods used in this study found differences in the two sub-watersheds despite both not having intensive anthropogenic pressures is very encouraging and indicates that they should be adopted by the responsible authorities in Greece and other Mediterranean countries. Differences were also evident between the sampling stations, with higher elevation stations having the best conditions, since anthropogenic disturbances were limited. Overall, integrating different assessment methods provides a complementary view on the watershed and riparian areas conditions that should help in the development of sustainable management plans.
Ecotourism is gaining greater interest, especially in mountainous and riparian areas of Greece and the Mediterranean. Similar datasets should be collected for riparian areas throughout Greece and even the Mediterranean to be able to capture the current conditions. Such datasets will provide the necessary information to implement ecotourism plans and best management practices to mitigate anthropogenic and climate change impacts.
For the two studied sub-watersheds, although the overall condition of the riparian areas is relatively good compared to other riparian areas of Greece, some measures need to be taken. Firstly, the same or similar assessments need to be continued in both subwatersheds in order to be able to develop a long-term dataset of the conditions. Another key is to better monitor the anthropogenic pressures. Specifically, the livestock grazing, hunting, fishing and other ecotourism activities need to be checked and measured. In the past, these activities were limited, but based on this study, there are some indications that they are negatively impacting the riparian areas in the low elevations that have more residents living nearby and a more extensive road network that leads to more visitors. A new management plan should be developed to regulate such activities. Specifically, after the current pressures are recorded, the number of visitors for ecotourism might need to be regulated to avoid degradation and promote sustainability, while best management practices should look into minimizing the impacts of residents (e.g., grazing, wood harvesting) from the adjacent communities to maintain the naturalness of these mountainous sub-watersheds and their riparian areas. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.