Globally, freshwater gastropods inhabit a wide range of aquatic environments, from lotic systems such as mountain streams to lentic habitats like wetlands (Johnson et al. 2013). Although they have adapted to diverse physical environmental changes, most freshwater gastropods exhibit limited mobility, and many species are endemic with restricted geographic ranges (Clements et al. 2006). Koreoleptoxis globus (Martens 1886), commonly known as the Korean freshwater snail, is an endemic species of South Korea. It inhabits streams in the northern regions of Gangwon State, Chungcheongbuk-do, and Yeoncheon-gun in Gyeonggi-do (National Institute of Biological Resources 2022). K. globus is primarily found in relatively deep, fast-flowing sections of upper reaches of streams (Min and Lee 2005). However, this information is based on qualitative observations, and quantitative studies are lacking. K. globus is predominantly herbivorous, primarily consuming periphytic algae, and unlike terrestrial and marine gastropods, it does not exhibit predatory behavior (Burch 1989). As a result, it has a significant influence on algal primary production in river ecosystems, playing a vital role in aquatic food webs and nutrient cycling (Covich et al. 1999). Major threats to the survival of K. globus include overharvesting for food, habitat alterations due to dam and weir construction, and water pollution. Consequently, the National Institute of Biological Resources (NIBR) of the Ministry of Environment (MOE) of South Korea has classified it as “Vulnerable (VU)” in the “Red Data Book of Republic of Korea”. The MOE designated it as an endangered candidate species for the first time in 2022.

The aim of this study was to quantitatively determine a physically suitable habitat range for K. globus using probability density functions (PDF). Unlike discrete analysis based on occurrence data, continuous probability distribution analysis facilitates a continuous interpretation of habitat characteristics across all values of environmental factors and allows for the presentation of habitat suitability as a probabilistic function (Kong and Kang 2023). Previous applications of Weibull probability distribution models for analyzing benthic macroinvertebrates have included interpretations of the relationship between the survey area and the number of species (Kong and Kim 2015), calculation of physical habitat suitability indices (Kim and Kong 2018; Kong and Kim 2017), and studies on the physical habitat characteristics of Koreoleptoxis nodifila (Kim et al. 2022) and Cambaroides similis (Kim et al. 2023). Kong and Song (2023) used various distribution models to assess the suitability of the physical microhabitats for Ephemera. In this study, we applied 12 different distribution models to quantitatively analyze the habitat suitability for K. globus for the first time.

As the frequency of concentrated rainfall events increases due to climate change, the risk of rapid fluctuations in stream discharge increases, which may compromise the stability of physical habitats (Kong and Kim 2016). Additionally, threats from artificial habitat disturbances, such as stream construction activities, are increasing. Hydrological changes, such as alterations in water flow, can disrupt the entire food chain of freshwater ecosystems (Joo et al. 2008). Temporary increases in discharge can lead to disturbances that reduce benthic algal populations, which, in turn, may decrease the population of benthic macroinvertebrates that depend on benthic algae for food (Gafner and Robinson 2007). Due to their limited mobility, gastropods have a restricted ability to escape from threats, making it difficult for aquatic ecosystems to recover once they are disrupted (Heino et al. 2002; Poff 1997). The aim of this study was to provide essential information for habitat conservation by delineating the physical habitat characteristics of K. globus.

Field survey

A total of 34 sites were sampled to determine the density of K. globus (individuals m^–2) from April 2024 to May 2024. These survey sites were selected based on the locations where K. globus was previously identified in the “National Ecosystem Survey (2020–2022)” and the “Inland Wetland Survey (2022)” conducted by the National Institute of Ecology, as well as the “Stream/River Ecosystem Survey and Health Assessment (2008–2022)” conducted by the National Institute of Environmental Research (NIER). Surveys were conducted at 340 sampling units across these 34 sites (Fig. 1), where the presence of K. globus was confirmed, targeting various physical environments, including riffles, runs, and pools. Following the guidelines for surveying and assessing benthic macroinvertebrates outlined in NIER Announcement No. 2016-372, K. globus was quantitatively sampled using a Surber sampler (30 cm × 30 cm). After counting the number of individuals in the field, the population density (individuals m^–2) were calculated for quantitative analysis. Physical environmental factors considered in this study included current velocity (U, cm s^–1), water depth (D, cm), and streambed substrate particle size (Φm). Water depth was measured usingh a 100 cm steel ruler. The current velocity was calculated using the difference in water height observed when a ruler was placed perpendicular to the direction of stream flow, following the formula of Craig (1987). Substrate composition was visually assessed by measuring the area ratio of particle sizes on the streambed surface according to the system outlined by Cummins (1962) and converted into phi (Φ) values as per Kong and Song (2023).

Figure 1. Map of South Korea, the southern half of the Korean Peninsula bordering the East Sea and Yellow Sea, Northeast Asia; and survey sites of Koreoleptoxis globus (presence: 11 red circles, absence: 23 yellow circles).

Combined relative abundance analysis

Environmental factors were categorized into intervals (water depth: 15 cm, current velocity: 10 cm s^-1, streambed substrate (Φμ): 1), and the combined relative abundance (C_i) was calculated for each interval. The relative mean number of individuals (A¯i) for each interval was defined as specificity, while the relative frequency of occurrence (R_i) was defined as fidelity. The product of these two values is referred to as combined relative abundance (Dufrêne and Legendre 1997). C_i considers R_i and A¯i as independent factors while integrating them comprehensively. Consequently, when both R_i and A¯i increase within a specific interval, the corresponding C_i for that interval increases significantly. Therefore, C_i tends to more accuratetly reflect the occurrence patterns across habitat intervals based on the distribution characteristics of the organisms.

The relative frequency of occurrence (R_i) in interval (i) of K. globus was the ratio of the number of sampling units where a taxon occured in each class interval of the environmental factor (n_i) to the total number of sampling units in that class interval (N_i) (Eq. 1).

The relative mean number of individuals (A¯i) in interval (i) was calculated by dividing the sum of the number of individuals of K. globus (A_ij) observed at each sampling unit (j) within the class interval of an environmental factor by the total number of sampling units (N_i) in that class interval (Eq. 2).

The combined relative abundance (C_i) was defined as the product of the relative frequency of occurrence and the mean number of individuals (Eq. 3).

The combined relative abundance must first be converted into a discrete probability mass function (PMF) and then into a continuous PDF to continuously interpret the occurrence of taxa in relation to changes in physical environmental factors. Thus, the PMF value was determined by multiplying combined relative abundance by the range value of class interval (Δxi) and summing these products, and then dividing by the combined relative abundance for the corresponding class (Eq. 4).

According to Equation 4, the cumulative mass function (CMF) up to class μ was calculated using Equation 5. In Equation 5, the value of the CMF for all classes (m = k) equals 1.

Selection of probability distribution model

The probability distribution models applied in this study were the same as those used by Kong and Kang (2023), comprising eight types divided into threshold values. Theses two-parameter models, i.e., exponential distribution (type 1), truncated normal distribution (type 2), logistic distribution (type 3), and Gumbel distribution (type 4); three-parameter models, i.e., log-normal distribution (type 5), Weibull distribution (type 6), and gamma distribution (type 7); and a four-parameter model, i.e., beta distribution (type 8). Additionally, four inverted models were included to interpret the negatively skewed distributions. These were inverted exponential distribution (type 1-1), inverted log-normal distribution (type 5-1), inverted Weibull distribution (type 6-1), and inverted gamma distribution (type 7-1). Thus, a total of 12 probability distribution models were used to analyze the physical habitat suitability for K. globus.

The domain of random variables in each probability distribution model may differ from the actual range in which organisms are distributed. For example, while the domain of random variables in a normal or logistic distribution extends from -∞ to ∞, the mean and variance derived from these probability distribution models should be calculated based on the actual range of environmental factors where the organisms occur. Organism distribution patterns can be classified into the following cases: those with a dummy distribution when the range of environmental factors is from 0 to ∞ or 0 to a (Fig. 2A), those without a dummy distribution (Fig. 2B), those in the positive domain ranging from a to ∞ or a to b (Fig. 2C), and those that include negative domains, such as a to ∞, -∞ to b, or a to b (Fig. 2D). Since organisms can appear even at a current velocity of zero, their distribution may follow the pattern shown in Figure 2A, or the threshold-based pattern in Figure 2B. In contrast, since organisms cannot occur at a water depth of zero, their distribution may follow the pattern in Figure 2C, or the threshold-based pattern in Figure 2B. Furthermore, since the mean grain size of the substrate, when converted to phi values (Φm), can take negative values, the distribution may follow the pattern shown in Figure 2D. The commonly used PDF, cumulative distribution function (CDF), and statistical functions can be applied to the domain of random variables for types that exclude dummy distributions. However, for dummy distributions, the PDF and CDF were adjusted to ensure that the domain of the random variables was 0 or greater, and statistical functions were derived accordingly (Kong and Song 2023). Detailed information on the calculation methods for the CDF, PDF, mode, mean, and variance of each probability distribution model can be found in Kong and Song (2023).

Figure 2. Distribution patterns of aquatic organisms based on environmental factors. (A) With a dummy variable, (B) Without a dummy variable and non-shifted, (C) Without a dummy variable and shifted, (D) Transformed, without a dummy variable and shifted.

Parameter estimation of the probability distribution

The parameters of the probability distribution were derived using the solver function in Microsoft Excel, by minimizing the normalized root mean squared error (NRMSE) between the CMF and CDF (Eq. 6).

Physical habitat suitability assessment

The habitat suitability index (HSI) assessment method of the ‘United States Instream Flow and Aquatic Systems Group (Bovee 1986)’ was adopted. Suitability values of 1.0, 0.5, 0.1, and 0.05 were assigned to ranges of the probability distribution of 50%, 75%, 90%, and 95%, respectively. The optimal habitat range was determined on the basis of a threshold of 50%.

Water depth

The NRMSEs between the CDF and CMF of each probability distribution model for the water depth of K. globus are shown in Table 1. Water depth exhibited a weak negative skew, with an inverted log-normal distribution model providing the best fit (NRMSE, 2.41%). The inverted gamma distribution (NRMSE, 2.48%), beta distribution (NRMSE, 2.49%), and Weibull distribution (NRMSE, 2.94%) models also demonstrated good fits. Figure 3A illustrates the PMF of water depth for K. globus, and the PDF of the model showed the best fit. Central values were as follows: mean value of 35.8 cm and mode value of 44.4 cm, with a standard deviation of 15.8 cm in Table 2. The depth ranges for K. globus based on the HSI are shown in Table 3. Considering the range corresponding to an HSI value of 1 as the optimal habitat range, the lower and upper bounds of the 50% range in Figure 3B were considered the boundaries of the optimal habitat range. Accordingly, the optimal habitat range for the water depth of K. globus was 31.2 to 53.7 cm.

Table 1 . Normalized root mean squared error (%) between cumulative mass function derived from combined relative abundance based on water depth, current velocity, and mean diameter gradient of substrate grain and cumulative distribution function of each model.

Model	Water depth	Current velocity	Streambed substrate (phi value)
Exponential	26.10	17.52	4.12
N-exponential	14.94	9.11	20.25
Normal	4.52	1.12	3.06
Lognorma	4.94	1.02	0.76
N-lognormal	2.41	1.52	3.47
Logistic	5.53	1.60	2.52
Weibull	2.94	0.86	0.55
N-Weibull	3.12	0.98	0.89
Gamma	4.76	1.06	0.64
N-gamma	2.48	1.98	3.34
Beta	2.49	0.71	0.96
Gumbel	7.99	1.41	0.81

Table 2 . Best-fit probability distribution models and corresponding statistical metrics for water depth, current velocity, and mean diameter gradient of substrate grain.

	Water depth (cm)	Current velocity (cm s-1)	Streambed substrate (phi value)
Model	N-lognormal	Beta	Weibull
Mean	35.8	63.0	–5.7
Mode	44.4	75.6	–6.7
Standard deviation	15.8	21.4	1.5

Table 3 . Habitat suitability range of Koreoleptoxis globus based on water depth, current velocity, and mean diameter gradient of substrate grains.

	Range of habitat suitability
	50%	75%	90%	95%
Water depth (cm)	31.2–53.7	20.2–58.0	9.5–60.9	4.3–61.9
Current velocity (cm s-1)	57.8–88.8	44.3–94.7	31.1–98.0	23.7–99.1
Streambed substrate (phi value)	–7.4 to –5.6	–7.7 to –4.7	–7.8 to –3.5	–7.8 to –2.5

Figure 3. Water depth. (A) Probability mass function (PMF) and probability density function (PDF) using inverted lognormal distribution, and (B) adjusted habitat suitability index (HSI) for Koreoleptoxis globus.

Current velocity

The NRMSE between the CDF and CMF for the current velocity is listed in Table 1. The current velocity exhibited a negatively skewed distribution, with the beta distribution model providing the best fit (NRMSE = 0.71%). The Weibull distribution (NRMSE = 0.86%) and inverted Weibull distribution (NRMSE = 0.98%) models also demonstrated an excellent fit. Figure 4A illustrates the PMF of the current velocity for K. globus and the PDF of the model with the best fit. Central values were as follows: mean value of 63.0 cm s^-1 and mode value of 75.6 cm s^-1 with standard deviation of 21.4 cm s^-1 in Table 2. The composite frequency of K. globus was very low in sections without water flow and decreased drastically at current velocities above 90 cm^-1. Figure 4B shows the HSI for the current velocity of K. globus. The optimal habitat range for current velocity was between 57.8 and 88.8 cm^-1 in Table 3.

Figure 4. Current velocity. (A) Probability mass function (PMF) and probability density function (PDF) using beta distribution, and (B) adjusted habitat suitability index (HSI) for Koreoleptoxis globus.

Streambed substrate

The NRMSE between CDF and CMF for the mean particle size of the substrate is listed in Table 1. The streambed substrate exhibited a positively skewed distribution, with the Weibull distribution model (NRMSE of 0.55%) providing the best fit. The gamma (NRMSE of 0.64%) and log-normal (NRMSE of 0.76%) models also demonstrated an excellent fit. Figure 5A shows the PMF of the mean particle size of the substrate and the PDF of the best-fitting probability distribution model. Central values were as follows: mean value of −5.7 and mode value of −6.7 with a standard deviation of 1.5 in Table 2. Figure 5B shows the HSI for the streambed substrate of K. globus. The optimal habitat range for the mean particle size of substrate was between −7.4 and −5.6 in Table 3. Based on substrate data measured at sampling points for K. globus in this study, the average substrate type was strongly characterized as ‘pebbly cobble’ (Fig. 6).

Figure 5. Median diameter. (A) probability mass function (PMF) and probability density function (PDF) using weibull distribution, and (B) adjusted habitat suitability index (HSI) for Koreoleptoxis globus.

Figure 6. (A) Textural classification and (B) composition of stream substrates of Koreoleptoxis globus.

As a stream flows from headwaters to large rivers, geomorphic and physical variables such as stream flow and channel morphology influence the composition of biological communities (Vannote et al. 1980). The key physical factors that affect benthic macroinvertebrates’ microhabitats include current velocity, substrate type, and water depth (Kim 2014; Kong and Kim 2017; Orth and Maughan 1983; Pan et al. 2015). Coarse-particle substrates such as rocks and gravel are predominant in areas with high current velocities (Church 2002). Incontrast, fine sediments such as sand are more prevalent in areas with low current velocities (Colby 1964). Clingers that attach to the surface of the streambed substrate are dominant in coarse substrates, whereas the habitat orientation group (HOG) shifts towards burrowers that dig into and inhabit the sediments in fine substrates (Kim et al. 2017). Kwon et al. (2013) classified the HOG of K. globus into clingers and sprawlers. Burrower populations, which prefer slow current velocities, are most strongly correlated with the substrate type, whereas clingers are most influenced by current velocity (Kim et al. 2017). However, because these findings were based on the average data for each HOG, we quantitatively analyzed and evaluated the habitat suitability of K. globus, a clinger species, in relation to current velocity, water depth, and substrate type in this study.

Using 12 probability distribution models to evaluate the physical habitat suitability of K. globus, we found that the beta distribution, which is flexible and can be expressed in various forms (Kong and Song 2023), along with weibull distribution, provided the highest overall suitability, with NRMSE of less than 3%. The beta distribution models were particularly well suited for cases exhibiting negative skewness, whereas weibull distribution models provided a better fit for cases with positive skewness. The inverted Weibull distribution, gamma distribution, inverted gamma distribution, log-normal distribution, inverted log-normal distribution, and normal distribution demonstrated acceptable suitability, with NRMSE values of less than 5%. However, in this study, the model with the smallest NRMSE value for each environmental factor was selected as the representative model for the following results, i.e., inverted log-normal distribution for water depth, Weibull distribution for current velocity, and beta distribution for the streambed substrate. The distribution patterns of organisms in relation to environmental factors can be predicted using probability distribution models, with commonly applied models including the exponential, normal, log-normal, and logistic distributions (Fisher et al. 1943; Godsoe et al. 2016; Matthews and Whittaker 2015; Preston 1948). However, the appropriateness and flexibility of these distribution models may vary based on the range of environmental variables considered (Ahmadi‐Nedushan et al. 2006). As the diversity of environmental variables increases, the number of applicable models may also expand (Golestani and Gras 2013; Lenart and Missov 2016). Consequently, selecting the most suitable probability distribution model is essential in habitat suitability analysis. Failure to account for the specific distribution characteristics of the target species when calculating the HSI can result in skewed predictions, and thus requires careful consideration.

In a study by Kim (2014) on benthic macroinvertebrates in Gapyeong Stream, the primary occurrence range for clingers was reported as a current velocity of 6 cm s^-1 to 15 cm s^-1, a water depth of 20 cm to 50 cm, and a cobble substrate. Similarly, in a study focused on the Han River by Kim (2017), median values of environmental factors for the clinger genus Semisulcospira were a current velocity of 29 cm s^-1, water depth of 44 cm, and substrate particle size of −3.84 (gravel). These findings suggested that the physical habitat of clingers generally consists of coarse substrates in shallow-flowing water. However, in a study by Kim et al. (2017) on microhabitat classification of benthic macroinvertebrates, clingers were reported to occur even at current velocities of up to 120 cm s^-1. This indicated that the preferred range of current velocities varies among clinger species.

In this study, mean values of environmental factors for K. globus were a current velocity of 63.0 cm s^-1, a water depth of 35.8 cm, and a substrate particle size −5.7 (pebble). The microhabitat of K. globus was similar to that of clinger species, which prefer shallow, wadeable depths with well-developed coarse-particle substrates. However, K. globus prefered higher current velocities, even among the clinger species. Given that mode current velocity for K. globus was 75.6 cm s^-1, but considering that the observed maximum current velocity range is 90 cm s^-1 to 100 cm s^-1, it can be inferred that the species has expanded its ecological niche to adapt to habitats with fast current velocities. Conversely, K. globus was occasionally found at current velocities below 15 cm s^-1, all occurring in riffle-adjacent runs or pool sections, where the current velocity temporarily decreased. Although current velocity was low in these areas, the physical habitat conditions were not significantly different from those in the riffle sections. Kim et al. (2017) reported that the number of individuals of clinger and burrower species may be reversed at water depths of 40–50 cm, and our results also confirmed a drastic decrease above 60 cm. Kong and Kim (2016) classified streambed substrate characteristics into ten categories using Shepard's (1954) method, applying a triangular diagram based on the composition ratios of sand, silt, and clay. They proposed five types of substrate preferences based on the average substrate particle size: lithophilous, psephophilous, moderate, psammophilous, and pelophilous (Table 4). By substrate type, psephophilous substrates predominated, with 72% of the total, followed by moderate substrates at 16%, lithophilous substrates at 10%, and psammophilous substrates at 2%. Kim et al. (2017) classified microhabitats into nine types based on the distribution characteristics of benthic macroinvertebrate HOG in terms of water depth, current velocity, and substrate type in streams. Considering the distribution patterns observed in this study, the mode value of environmental factors for K. globus was applied to define the habitat type as a riffle (Zone IV). This provided quantitative evidence that K. globus prefers riffle sections and coarse- particle substrates.

Table 4 . Number and percentage of sampling units, mean diameter gradient of substrate grains and classification of lithophilic series according to substrate types.

No.	Substrate type	Average of mean particle diameter (Φμ)	Number and percentage (%) of sampling units of Koreoleptoxis globus		Lithophilic series
1	Bouldery cobble	–7.1	6	72	Lithophilous
2	Pebbly cobble	–6.1	66
3	Gravelly cobble	–5.4	6	10	Psephophilous
4	Cobbly pebble	–5.2	4
5	Pebble	–4.7	-	16	Moderate
6	Copegra	–4.5	15
7	Gravelly pebble	–4.1	1
8	Cobbly gravel	–3.5	1	2	Psammophilous
9	Pebbly gravel	–3.0	1
10	Gravelly clay	–0.9	-	-	Pelophilous

Revised from the article of Kong and Kim (J Korean Soc Water Environ. 2016;32(1):1-14).

When selecting habitats, food source availability is one of the most critical factors for organisms. Variations in food sources across different habitats can be inffered from fatty acid composition extracted from tissues of target organisms (Shimma and Taguchi 1964), Fatty acids such as 20:5ω3, 18:2ω6, and 18:3ω3 in benthic macroinvertebrates can serve as biomarkers for algae and phytoplankton (Shin et al. 2012). According to Lim et al. (2009), the ratio of these fatty acids in K. globus (12.8%), which inhabits riffle sections, is approximately 1.7 times higher compared to that in species of the genus Semisulcospira (7.4%), which inhabit run sections. Additionally, chlorophyll a and chlorophyll b contents in K. globus (1.2 mg/100 g) was about twice that found in Semisulcospira coreana and Semisulcospira tegulata (0.6 mg/100 g), indicating a relatively stronger preference for algae attached to rocks than in species of the genus Semisulcospira. Generally, riffle sections experience less deposition of fine organic particles because of the high current velocity, whereas coarse-particle substrates that are essential for the growth of diatoms make up a significant portion of the habitat. This suggests that the primary reason K. globus prefers riffles is the availability of food sources. Increasing the frequency of intense rainfall events due to climate change can significantly reduce the diatom biomass (Lee et al. 2017). Moreover, the influx of highly concentrated turbid water from civil engineering activities can block light and inhibit the growth of diatom communities (Kim et al. 2011). Future research on the conservation of K. globus habitats should include additional studies that focus on food sources and factors that threaten them.

Analysis of the suitability of the physical microhabitat for the endangered candidate species K. globus (Mollusca, Gastropoda) using probability distribution models. We identified the physical microhabitat suitability of K. globus as being a "riffle" section characterized by coarse-particle substrate, fast current velocity, and shallow water depth in streams. This is the first instance where the physical microhabitat suitability of K. globus has been statistically validated. However, we assessed the suitability of each physical environmental factor for K. globus individually and did not evaluate the interactions between these factors. Consequently, further research is necessary to understand how each environmental factor influences the survival of K. globus. Moreover, additional studies are needed to investigate the impacts of increased frequency of intense rainfall due to climate change and the influx of high-concentration turbid water resulting from civil engineering activities on the availability of food sources for the effevtive conservation of K. globus habitats.

Not applicable.

CDF: Cumulative distribution function

CMF: Cumulative mass function

PDF: Probability density function

PMF: Probability mass function

HOG: Habitat orientation group

HSI: Habitat suitability index

NRMSE: Normalized root mean squared error

JYK conducted formal analysis, visualization and writing of the original draft. JKM and YJK were responsible for data curation, as well as writing, review and editing. YSP handled funding acquisition, project administration, resource management, and contributed to writing, review and editing. DK contributed to conceptualization, supervision, formal analysis, and writing, review and editing. All authors have read and approved the final manuscript.

This study was supported by the project title ‘Survey on Endangered Candidate Species (NIE-C-2024-103)’ funded by the Ministry of Environment.

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Not applicable.

Not applicable.

The authors declare that they have no competing interests.