Chorological and Ecological Differentiation of the Commonest Leech Species from the Suborder Erpobdelliformes (Arhynchobdellida, Hirudinea) on the Balkan Peninsula

This study is the result of extensive investigations of leeches on the Balkan Peninsula. Our aim was to detect actual and potential (modeled) distributions of common Erpobdellidae species, and to identify their ecological differentiation with respect to the altitudinal and waterbody type gradient. Although widespread, these species rarely live together. Intense competition is avoided by preferences for different types of habitats. This was confirmed by Pearson correlation analyses that yielded negative results. Differentiation of these species was clarified by the results of logistic Gaussian regression analyses. While Erpobdella octoculata and Dina lineata have a similar distribution along the altitudinal gradient, they prefer different waterbody types. Erpobdella vilnensis prefers higher altitudes than the other two species. Its preferred habitats are smaller rivers and streams located at altitudes from 400 to 1000 m a.s.l. Although present in all waterbody types, large lowland rivers and standing waterbodies are the preferred habitats of E. octoculata. Fast-flowing springs and streams are mostly inhabited by D. lineata. While the distribution of the species overlaps to a large degree, the ecological preferences of species differ significantly and thus they can be used as confident typological descriptors and indicators of ecological status.


Introduction
Despite a long history of leech investigations on the Balkan Peninsula , important data on their distribution are scarce.
According to Sket and Trontelj [25], two families of Hirudinea (Erpobdellidae and Glossiphoniidae), dominate in the Palearctic biogeographical region. The family Erpobdellidae is represented by 25 taxa in the Balkan and neighboring areas. Among them, a large number of species are endemic and have a narrow distribution [4,11,13,[15][16][17][18]26]. Erpobdelliformes are dwellers of freshwater habitats and are macrophagous predators, feasting on molluscs, arthropods and other annelids, although cases of cannibalism have also been reported [27,28]. With their role in the regulation of a number of prey

Materials and Methods
Extensive fieldwork was conducted from 2010 to 2018. The leeches were collected using a benthological hand net (mesh size 500 µm); additionally, individuals were collected using tweezers, from hard substratum and vegetation. Each animal was relaxed in 10% ethanol, and then transferred to 70% ethanol for further analysis. Identification of the leeches was done according to Nesemann and Neubert [30], using two stereomicroscopes: a Nikon SMZ800N (Nikon Corp., Tokyo, Japan) (magnification 10-80×) and a Zeiss Stemi 2000-C (Carl Zeiss Microscopy GmbH, Göttingen, Germany) (magnification 6.5-50×). Small, juvenile or damaged individuals were not taken for analysis, since their identification at the species level was not possible. Material was deposited in the collection of the Institute for Biological Research "Siniša Stanković", University of Belgrade.

Study Area
The study area, extending from the Sava River (CRO) in the northwest to the Dojran Lake (NMCD) in the southeast, covers a large portion of the Balkan Peninsula, as shown in Figure 1. The waters of this area drain to three watersheds: the Zeta, Morača, Bojana (MNE) and Black Drim (NMCD, AL) rivers that drain into the Adriatic Sea (the Skadar Lake and numerous karstic springs in its basin also belong to this watershed). The waters of the Vardar River basin, along with the Dojran Lake, (NMCD) and Maritza River (BLG), belong to the Aegean watershed. Most of the studied area drains drain into the Black Sea through the River Danube. The Sava and Velika Morava rivers are two of the biggest tributaries of the Danube, draining major parts of Croatia, Serbia, and Bosnia and Herzegovina. Larger rivers, such as the Bosna, Una, Sana (BIH), Tara (MNE) and Drina (SRB), drain the waters from the western parts of the investigated area into the Sava River. In the investigated area, a variety of waterbodies were studied, from springs and small rivers to large waterbodies such as the Lakes Ohrid and Skadar. Aside from the geographical parameters (recorded with GPS), the hydromorphological properties of the sampling sites were assessed using six categories of waterbodies, according to the modified national typology of surface waters of Serbia [35]. Reservoirs, lakes, ponds and other standing waterbodies were marked as type 1 (T1); large lowland rivers (the Danube, lower stretches of the Sava, Velika Morava and Vardar) were marked as type 2 (T2); lower stretches of their bigger tributaries (the Bosna, Sana, Una, Drina, Zapadna Morava, Južna Morava, Kolubara and Tara) were assigned to type 3 (T3); type 4 (T4) incorporated various watercourses of medium size (wadeable rivers) at elevations below 500 m a.s.l.; those at higher elevations were categorized as type 5 (T5); and small waterbodies, like springs and upper stretches of streams, belonged to type 6 (T6). Each site was categorized to the appropriate waterbody type according to its characteristics. The basic properties of these waterbody types are represented in Table 1.

Statistical Analyses
Spatial data analysis was performed with DIVA-GIS 7.5 software [36]. The climate parameters (mean temperature and mean precipitation of sampling sites) were extracted from the WorldClim database [37] at a resolution of 30 arc-seconds (~1 km 2 ). To model the potential distribution of the analyzed species, the Maximum Entropy (MaxEnt) method [38][39][40] was used. This method belongs to a family of species distribution models that relates field observations to environmental predictor variables [41]. The MaxEnt method operates using binary presence/absence data. Distributions of analyzed species were modeled using MaxEnt software (open source software available at https: //biodiversityinformatics.amnh.org/open_source/maxent/), version 3.3.e [38].
Both altitude and hydromorphological types may be considered as complex gradients (sensu Whittaker) [42]. A complex gradient represents an assemblage of environmental factors that change together. A set of climate factors (e.g., temperature, moisture and precipitation) change predictably along the altitudinal gradient. Waterbody types that range from 1 to 6 as shown in Table 1 represent a complex hydromorphological gradient. A set of environmental factors (flow velocity, dissolved oxygen concentrations, load of nutrients, and bottom properties) change predictably from waterbody type 1 to waterbody type 6 [43,44]. The logistic Gaussian regression [45,46] was used to detect the response and ecological preferences of the analyzed species, along the altitudinal gradient and gradient of waterbody types.
The logistic Gaussian regression was performed using FLORA software (software freely available upon request at branko@ibiss.bg.ac.rs), updated version [47].

Results
During nine years of field investigations, 2202 individuals of E. octoculata, E. vilnensis and D. lineata were collected from 229 localities (Table S1). The water from 57 of the sampling sites drained into the Adriatic Sea; 14 sites belonged to the Aegean watershed; and the majority of sampling sites (158) drained into the Black Sea through the Danube River, as shown in Table 2. E. octoculata was recorded on 99 sites; 23 sites were inhabited by E. vilnensis, and D. lineata was the most frequently detected species with 142 records, as shown in Table 3.  T1  3  3  18  24  T2  0  1  7  8  T3  2  0  15  17  T4  2  5  33  40  T5  4  2  42  48  T6  46  3  43  92  Σ  57  14  158  229 T-Type of waterbody; Ws-Watershed.

The Actual Distribution of the Analyzed Species
The distribution of species in each investigated area is presented on the map in Figure 1. All three species had a similar distribution along latitudinal and longitudinal gradients, as shown in Table 4; however, due to their widely overlapping distributions, the chorological separation of these species was obscured. The analyzed species differed with respect to their distribution along the altitudinal gradient. E. octoculata and D. lineata occured at similar altitudes, whereas E. vilnensis prefered higher altitudes ( Figure 2). Most of the sampling sites belonged to the Black Sea and Adriatic watersheds. D. lineata dominated in sites that belonged to the Adriatic Sea watershed (Figure 3). A high percentage of the sites inhabited by E. octoculata belonged to the Aegean Sea. The distribution of the analyzed species with respect to waterbody type is presented in Figure 4.

The Potential Distribution of the Analyzed Species
In addition to the actual distribution, we analyzed the potential (i.e., the most probable) distribution of the investigated species. Results of the MaxEnt analysis indicated that E. octoculata preferred large waters in the lowlands, as shown in Figure 5. D. lineata had a narrower range, as shown in Figure 4. Compared to E. octoculata, D. lineata occupied higher altitudes. In lowland areas, D. lineata avoided large waterbodies. Instead of large rivers and lakes (the rivers of the Danubian and the Thracian plains, the Danube, Maritza (BLG) and Vardar (NMCD) rivers), D. lineata preferred fast streams (in mountainous regions) and springs in lowlands. Compared to E. octoculata and D. lineata, E. vilnensis preferred habitats at higher altitudes (>400 m a.s.l.). However, this species avoided altimontane regions ( Figure 4). Such potential distributions correspond to actual distributions of the analyzed species.

The Ecological Differentiation of Analyzed Species
Despite their broadly overlapping distribution, the analyzed species were negatively correlated, as shown in Table 5. E. vilnensis was almost uncorrelated with E. octoculata and D. lineata (r = −0.145 and r = −0.247, respectively). However, a negative value of correlation coefficient between E. octoculata and D. lineata (r = −0.625) indicated a strong competitive interaction among these taxa. Gaussian logistic regression was used to assess the ecological similarity among the analyzed species. The distribution overlap of the species along the altitudinal gradient is presented in Figure 6. The optimal altitude for E. octoculata and D. lineata was about 450 m a.s.l. E. vilnensis had a preference for higher altitudes when compared to the aforementioned species. The ecological tolerance (expressed as the standard deviation of a logistic Gaussian curve) of E. octoculata was greater than that of D. lineata.
The greatest ecological tolerance with respect to altitude was observed for E. vilnensis. The response curves of E. octoculata and D. lineata along the altitudinal gradient were similar. Due to overlapping distributions, competition between these two species was intense, and such a situation was confirmed by a negative Pearson correlation. Gaussian logistic regression of the analyzed species with respect to the waterbody gradient is presented in Figure 7. The ecological optimum of E. octoculata with respect to waterbody type (3) indicates that this species preferred the large tributaries of large lowland rivers, with a slow to medium flow, and a riverbed with fine to medium-sized sediments. E. vilnensis preferred smaller and faster waters (the optimum was 4.5), as shown in Figure 7. Finally, D. lineata dominated in fast-flowing rivers and both lowland and mountainous streams and springs (the optimum was 5), as shown in Figure 7.
E. octoculata and D. lineata occured in all types of waterbody, as shown in Figure 7. Therefore, their ecological tolerance with respect to waterbody type was wide. A narrow ecological tolerance was recorded for E. vilnensis. This species was absent in large rivers. Moreover, E. vilnensis was recorded in only one lake (Lake Prespa).

Discussion
Despite their overlapping distributions, the analyzed species differ with respect to their geographic and ecological preferences. E. octoculata is widely distributed (from the upper stretch of the Sava River to the Maritza River along the longitudinal gradient, and from the Sava and the Danube Rivers to the Prespa and Dojran Lakes latitudinally). E. vilnensis and D. lineata have narrower distributions. All three species span a wide altitudinal range. Higher altitudes are preferred by E. vilnensis, while E. octoculata and D. lineata prefer lower altitudes compared to aforementioned species.
Potential distributions of the analyzed species obtained by MaxEnt analyses were in agreement with actual distributions recorded in the field. E. octoculata prefers large waterbodies in lowlands. This finding is in accordance with several articles that describe its ecology and preference to α-mesosaprobic, β-mesosaprobic and polysaprobic waters with a slower flow [30,33,34,48,49].
Although D. lineata occurs in lowland areas, this species avoids large rivers and standing waters. At lower altitudes, it is usually found in springs (e.g., numerous springs in the Skadar Lake basin). The preference of D. lineata for fast-flowing waters presented in this paper corresponds to our previous findings [24]. These habitats are not typical for populations of D. lineata in northern and central Europe [50][51][52]. Sket et al. [4] described three subspecies of D. lineata in the Balkans: D. lineata lacustris inhabits glacial lakes in western North Macedonia, D. lineata montana is usually found in the alpine areas of the Prokletije and Komovi mountains and D. lineata dinarica is widespread in the Dinaric Alps. Grosser et al. [15,18]  Compared to E. octoculata and D. lineata, E. vilnensis prefers habitats at higher altitudes. However, this species avoids both lowlands and altimontane regions. It was recorded in all waterbody types except large lowland rivers. Considering the lake habitats, E. vilnensis was found only in Lake Prespa, located at 840 m a.s.l. Records of E. vilnensis at altitudes lower than 400 m a.s.l. are for rivers that flow through large urban settlements (SRB-Južna Morava, Nišava and BLG-Maritza). These rivers are under anthropogenic pressures, and are often exposed to water quality degradation [53,54]. Tolerance to mesosaprobic conditions is known in the case of E. vilnensis [48,55,56]. It can be concluded that in these cases, the combination of the faster-flowing and higher saprobic state has favored E. vilnensis in competition with other species. Many articles support the conclusion that E. vilnensis prefers smaller water bodies located at medium-high and higher altitudes [55,57,58].
Despite their broadly overlapping distribution, the analyzed species are negatively correlated. This indicates strong competitive interactions between the species. Due to the overlapping distribution along the altitudinal gradient, competition between D. lineata and E. octoculata is intense. However, these species are clearly separated along the gradient of waterbody types. E. octoculata prefers lowland rivers and their large tributaries with slow to medium flow rates and a bottom with fine to medium-sized substrate. D. lineata dominates in fast rivers and both lowland and mountainous streams and springs with a rocky bottom.
Papers dealing with macrozoobenthic fauna as a whole rarely show the presence of E. vilnensis in the Balkans. Only when authors deal with leeches as the main focus of their investigations [4,6,17,23] does this species appear. This can be attributed to the misidentification of E. vilnensis as some similar species with dark paramedian stripes on the dorsum, such as D. lineata. A similar problem in the misidentification of this species has been reported by several authors throughout Europe [50,59,60].
The Assessment System for the Ecological Quality of Streams and Rivers throughout Europe using Benthic Macroinvertebrates (AQEM protocol) [61] prevails in assessments of water quality and implementation of the Water Framework Directive [62] throughout Europe. Many countries in the Balkans implement this protocol in the routine monitoring of water quality, but unlike many EU countries, they do not have modified protocols for local/national use [63]. Autecological preferences of D. lineata in the AQEM protocol indicate that it prefers the epi-and metapotamal sectors of rivers. These sectors are described as the slow-flowing upper and middle parts of lowland streams and rivers, with a fine sediment on the riverbed [64]. In relation to the saprobic index of Zelinka and Marvan, it prefers α-, βmesosaprobic and polysaprobic waters [65]. However, this contradicts the recorded habitats and preferences of D. lineata given in the present study. The utilization of these indicator values for D. lineata in the Balkans, in the assessment of water quality, could give false results. It is recommended that authorities dealing with the implementation of the WFD and the modification of protocols for assessing water quality make the necessary changes with regard to this species and its use as a bioindicator.
In light of the above, it is clear that these common leech species, with significantly overlapping ranges in the Balkans, rarely coexist in the same site. Each species has preferences for certain types of waterbodies, with specific conditions in which it outcompetes the others. On the rare occasions when they live in sympatry, this can be attributed to the presence of favorable microhabitats that provide optimal conditions for each species.