Could Ilıca Stream (Antalya-Türkiye) become a new stronghold for Garra turcica (Teleostei: Cyprinidae)? An evaluation of ecological and population parameters

Cyprinids belonging to the genus Garra are widespread in sub-Saharan Africa, the Middle East, South Asia, and Southeast Asia, and a total of 39 species have been described (Bayçelebi et al., 2021). There are nine species of the genus Garra in Türkiye (Kaya et al., 2026). Species of the genus Garra inhabit fast-flowing waters, often with fine gravel bottoms, and feed almost entirely on algae. They inhabit a wide range of water temperatures (15–40°C). They spawn in groups and breed in late spring and early summer (Krupp & Schneider, 1989).

Garra turcica Karaman, 1971, one of these nine species, is a fish species that naturally lives in the Kızıl (Mersin), Seyhan, and Ceyhan rivers and in small coastal streams south of Ceyhan in Türkiye (Bayçelebi et al., 2018). This species is known worldwide for its use in ichthyotherapy. Both species (Garra rufa and G. turcica) are known as “doctor fish” and form the basis of ichthyotherapy, which is used to alleviate the symptoms of skin conditions such as psoriasis and eczema (Gorzel & Kiełtyka-Dadasiewicz, 2019). “Doctor fish” can pose a potential threat in their habitats due to their wide water temperature tolerance (15–40°C) and high climatic adaptability (Krupp & Schneider, 1989; U.S. Fish and Wildlife Service 2019). A recent morphological and molecular study reported the species for the first time from the western Antalya basin (Kapız Creek), Türkiye, outside its natural range, suggesting that the individuals may have been released into the local aquatic environment after being used as “doctor fish” in tourist spas (Güçlü et al., 2025).

There are only two studies on the biology of G. turcica in its natural distribution area. One of these is a study solely on the length–weight relationship data (Seyhan Reservoir/Seyhan River), and the other is a study on some growth characteristics (Aksu Stream/Ceyhan River) (Çiçek et al., 2021; Ergüden, 2015). Thus, information on the growth characteristics of the species in its natural area is very limited, while there is no information on its reproductive characteristics. This study reports, for the first time, the growth and reproductive characteristics of G. turcica outside its natural range in a population established in the Ilıca Stream (Manavgat) in the eastern Antalya basin. Based on the baseline growth and reproductive data obtained, this study may also inform future research evaluating the ecology and potential spread of the species.

The presence of G. turcica (Fig. 1) in Ilıca Stream (Fig. 2) was documented through monthly sampling conducted from February to October 2025, during which a total of 505 individuals were collected using electrofishing gear. In each sampling event, a 600 m² stream area was systematically surveyed for 45 min to ensure standardized sampling effort. The study determined the length and weight distributions, length–weight relationships, condition factors, Von Bertalanffy growth equation (VBGF) parameters, fecundity, egg diameter, egg number, and age at first reproduction of the Ilıca Stream population. Additionally, seine and dragnet sampling was conducted in the dune-like areas of the Ilıca Stream to determine the species’ population density (area/time covered). The total length (TL) of each specimen was measured to the nearest 0.05 mm using a digital caliper, and body weight was determined to the nearest 0.01 g with a digital scale under laboratory conditions. Age determination was based on the examination of scales taken from the area between the left lateral line and the base of the dorsal fin.

Figure 1

(A) Sampling point-Ilıca Stream, (B) Garra turcica, (C) Mouth morphology of G. turcica.

Figure 2

Distribution areas and sampling point of G. turcica.

Female and male distributions were determined according to age. The overall ratio of males to females was evaluated with χ² test (α = 0.05) (Zar, 2010). TL and weight frequency distributions for all specimens were calculated. To construct the length–frequency distribution, the TL data of all fish specimens were grouped into class intervals. The class width (h) was calculated by dividing the range of observed lengths by the number of class intervals (k) using the following formula: h=Xmax−Xmink,h = {{{X_{max}} - {X_{min}}} \over k}, where X_max and X_min represent the maximum and minimum observed standard lengths, respectively, and k is the number of class intervals. The number of intervals was determined according to Sturges’ formula based on the sample size (n); k = 1 + 3.22* log10(n).

This method allows for a balanced representation of the dataset by minimizing the oversimplification or excessive fragmentation of the frequency distribution (Sturges, 1926; Zar, 2010).

The relation of weight to TL was established using the exponential regression equation W = a × TL^b, where W was the body weight in g, TL in cm, “a” is intercept, and “b” is the regression coefficient. R² (the coefficient of determination) was also estimated. In regression analysis, the R² value ranges from 0 to 1 and represents the portion of the total variance observed in the dependent variable that is explained by the independent variables. An R² value closer to 1 indicates that the model explains a large portion of the variability; a value closer to 0 indicates weaker explanatory power (Ricker, 1975). The growth of the G. turcica population was estimated with the following VBGFs: Lt=L∞(1−e{ −k(t−t0) }),Lt = {L_\infty }\left( {1 - {e^{\left\{ { - k\left( {t - {t_0}} \right)} \right\}}}} \right), where L_t is the total length in centimeters at age “t,” L_∞ the average asymptotic length in centimeters, k the body growth coefficient, “t_o” the hypothetical age, and “a” and “b” are the constants (Cren, 1951). The statistical significance level of the coefficient of determination (R²) and 95% confidence intervals (95% CI) of b were also estimated (Zar, 2010). The deviation of slope values from b = 3 (isometric growth) was assessed for all species, and Pauly’s t-test was performed (Pauly, 1984).

Pauly’s t-test was calculated as: t=SdlogTLSdlogW|b−3|1−r2n−2,t = {{S{d_{\log \pi L}}} \over {S{d_{\log W}}}}{{|b - 3|} \over {\sqrt {1 - {r^2}} }}\sqrt {n - 2} , where Sd_logTL is the standard deviation of the logTL values, Sd_logW is the standard deviation of the logW values, and n is the number of fish species used in the computation. The value of b is different from 3 if the t value is greater than the tabled t values for n – 2 degrees of freedom (Pauly, 1984). Measured TL and calculated TL in VBGF were evaluated with t-test (0.05). The average growth performance (Φ, phi prime) was calculated using the following formula (Gayanilo et al., 1988): Φ=logk+2*logL∞\Phi = log\;k + {2^*}log\;{L_\infty }

Fulton’s coefficient of condition factor was calculated by Sparre & Venema (1998) as Cf=(W/TL3)*100{C_f} = \left( {W/T{L^3}} \right)*100

The gonadosomatic index (GSI) was calculated as: GSI=(GW/W)*100,{\rm{GSI}} = {({\rm{GW}}/{\rm{W}})^*}100, where GW was the gonad weight and W is the total body weight of the fish (Gibson & Ezzi, 1980). Ovulation period was determined by GSI assessment. Average (absolute) fecundity per female was estimated gravimetrically from the number of mature oocytes in ripe females (spawning stage). The diameters of eggs were measured by means of a microscopic micrometer (Nikolsky, 1980). The method used to estimate the maturity level of mature females was based on fitting of the sigmoid, logistic curve. Estimation of the length at first sexual maturity is as follows (De Martini et al., 2000): First, it was plotted L against ln(1−PxPx)\ln \left( {{{1 - {P_x}} \over {{P_x}}}} \right) using simple linear regression to estimate values for “a” and “b,” where Px is the observed proportion of maturity at length x. Second, the mean lengths at 50% maturity was calculated by Lm=−abLm = {\raise0.7ex\hbox{{ - a}} \!\mathord{\left/ {\vphantom {{ - a} b}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{b}}.

Thereafter, the estimated proportions of maturity at length L were calculated using the equation and the first sexual maturity was drawn as PL=1001+eb(L−Lm),{P_L} = {{100} \over {1 + {e^{b\left( {L - {L_m}} \right)}}}}, where P_L is the estimated proportion of maturity at length L (De Martini et al., 2000).

The presence of age groups 0–5 years in Ilıca Stream indicates successful reproduction and the continuous addition of new individuals. While this age structure points to the existence of a locally established population, the possible introduction of individuals from outside areas cannot be entirely ruled out (Table 1). It is quite striking that 84.7% of the population consists of young individuals (ages 0 and 1) (Table 1). A “young and dynamic” age structure can be considered as a classic indicator of a population’s high growth and expansion potential (Gotelli, 2001). This suggests that, particularly in the early stages of colonization of a new habitat, resources are being allocated for rapid population growth, and the population has not yet reached equilibrium. This young and productive structure may continue to exert propagule pressure on both the Ilıca Stream and other water bodies in the basin in the near future (Sakai et al., 2001). The marked imbalance in the sex ratio (female:male = 0.32:1 or 1:3.13) provides important clues about the reproductive dynamics of the population (Table 1). The expected male:female ratio in fish populations is 1:1 (Nikolsky, 1980). The electrofishing method used may have increased the likelihood of male capture due to different behaviors of the different sexes (e.g., more aggressive male nest-site defense) or habitat preferences. The population may have started with a small founding group (e.g., a few “doctor fish” individuals) that included very few females. This genetic bottleneck may have caused the initial sex ratio to be reflected to the present day (Allendorf & Luikart, 2007).

The length–frequency distribution shows a unimodal structure. This indicates that the population is concentrated within a specific length range. Most of the data are clustered around the median length values (Fig. 3). Deviations may be due to growth stages or environmental factors. The weight distribution is generally skewed to the right (positive skew). This indicates that the population is composed of small/heavy individuals (Fig. 3). This skewness may be due to the abundance of young individuals. The significant negative allometric growth observed in the length–weight relationship (b = 1.756) is a significant deviation from the values found in the species’ natural range (b = 2.917 in the Aksu Stream-Ceyhan River; b = 3.202 in the Seyhan River) (Çiçek et al., 2021; Ergüden, 2015) (Fig. 4, Table 3). The observed negative allometry indicates that individuals gain less weight relative to their length, resulting in a more slender form. These measurements reflect the growth patterns of the studied population under the current environmental conditions. The value of R² = 0.768 indicates that length explains 76.8% of the variation in weight (Fig. 4). Compared with the >0.95 values reported in two natural populations (Çiçek et al., 2021; Ergüden, 2015) (Table 3), the current value appears lower, suggesting that growth and condition may be more variable within this population. However, as only two reference populations are available, this interpretation should be considered preliminary. The Fulton condition factor of the Ilıca Stream population (1.68) was higher than that of the Aksu Stream population in its native range (1.23) (Çiçek et al., 2021) (Table 3). This difference reflects the observed variation in the average body condition between the populations. The negative allometry observed in the Ilıca Stream population indicates that individuals gain less weight relative to length, resulting in a more slender form. Together, these measurements provide insights into the growth and condition patterns of the population under the current environmental conditions. However, no further inferences regarding energy allocation or reproductive investment can be drawn from these results, and future studies are needed to explore the relationship between body condition and reproductive status. The observed mean length values do not show a statistically significant difference from the model predictions (p > 0.05), demonstrating that the Von Bertalanffy model is valid for this population (Table 2). This high agreement means that the obtained L_∞, k, and t₀ parameters are reliable in describing the population (Zar, 2010). The theoretical maximum length of the Ilıca Stream population (L_∞ = 16.45 cm) is significantly lower than the population in its natural distribution area (Aksu Stream, L_∞ = 21.77 cm) (Çiçek et al., 2021), indicating that the individuals here have adopted a life strategy to reach a smaller final length (Fig. 5, Table 3). The growth performance index calculated for the Ilıca Stream population (Φ = 1.90) is lower than the natural population of the Aksu Stream (Φ = 2.17) (Çiçek et al., 2021) (Tables 2 and 3), confirming that the overall growth performance of the population is lower. Evaluation of the von Bertalanffy growth parameters indicates that the Ilıca Stream population exhibits relatively low L_∞ and early maturation compared with available data from the natural populations. These traits may suggest a tendency toward faster growth and earlier reproduction; however, further studies are needed to confirm any potential shift in the life-history strategy. This plasticity is one of the most important characteristics determining the success of non-native species in new environments (Richards et al., 2006; Winemiller & Rose, 1992). These findings may indicate a degree of ecological flexibility in G. turcica, particularly in terms of growth and reproductive allocation. Such flexibility could potentially facilitate its establishment and persistence in the Antalya basin. However, further ecological and long-term monitoring studies are needed to evaluate its invasive capacity and ecological impact. The GSI values begin to rise in both females (F) and males (M) by the end of February, peak in April–May, and then decline back to baseline values towards the end of summer (Fig. 6). This suggests that G. turcica has a “seasonally synchronous” reproductive strategy. The maximum GSI values in females are significantly higher than in males (Fig. 6). Male GSI values begin to rise earlier and remain relatively high throughout the breeding season. The seasonal dynamics of GSI differed between sexes, with females showing a sharper and more temporally restricted peak, while males exhibited an earlier increase and more prolonged elevation (Fig. 6). Sex-specific differences in reproductive timing were observed during the breeding season. Absolute fecundity in the studied population averaged 1533.44 ± 133.02 eggs. As no other fecundity data are available for this species, no conclusions can be drawn regarding the reproductive strategy or relative fecundity compared with other populations based on these data, and further studies are needed to investigate these aspects. The average diameter of a mature egg was determined to be 1.24 ± 0.11 mm, while the average diameter of an immature egg was 0.65 ± 0.12 mm. Although the oocyte size–frequency distribution has not been analyzed in detail, the distinct difference between immature and mature oocyte diameters, along with the significant increase and subsequent decrease in GSI values during the breeding season, supports the existence of synchronous oocyte development. This suggests that the species exhibits mass (single-time) ovulation rather than partial (batch) ovulation.

Figure 3

TL (A) and weight (B) in the frequency distributions of G. turcica. TL, total length.

Figure 4

Length–weight relationships of the G. turcica population.

Figure 5

L_∞ value according to VBGF of G. turcica. VBGF, Von Bertalanffy growth equation.

Figure 6

GSI values by months of G.turcica population in Ilıca Stream. GSI, gonadosomatic index.

The characteristic sigmoid (S-shaped) maturation curves observed (Figures 7 and 8) indicate that the length at which individuals first reach sexual maturity (Lm₅₀) can be reliably estimated in this fish population. The points corresponding to the 50% maturity rate, Lm₅₀ (69.47 mm) and tm₅₀ (1.03 years), represent the mean first reproductive threshold for this population. A particularly low level of these values is one of the most critical factors in shortening a population’s generation time. If a fish can produce offspring just 1 year after birth, its population size has the potential to increase geometrically rapidly. This significantly increases the population’s intrinsic growth rate (r) and reinforces its invasive potential (García-Berthou, 2007). The ability of an individual to contribute to the population size within its first year of life is key to its rapid dominance in new habitats, and this trait has been reported in many successful invasive fish species (e.g., G. holbrooki) (García-Berthou, 2007). This strategic adaptation could significantly increase the threat it poses to the local biodiversity. Furthermore, this trait means the species could spread rapidly to other water bodies within the Antalya basin. Therefore, developing preventive measures and early warning systems to address the “doctor fish” trade is critical.

Figure 7

Length at first maturity of males and females of G. turcica in Ilıca Stream.

Figure 8

Age at first maturity of males and females of G. turcica in Ilıca Stream.

The study area hosts four endemic (Alburnus baliki (VU), Capoeta antalyensis (VU), Pseudophoxinus alii (VU), Salariopsis burcuae (LC)) (Fig. 9), one native (Anguilla anguilla-CR) (Kaya et al., 2026), and four non-native (Alburnus escherichii, Pseudorasbora parva, Gambusia holbrooki and Carassius gibelio) species. The habitat characteristics where sampling was conducted in Ilıca Stream include partially sandy areas with marginal zones covered with reeds and vegetation, and a bottom structure dominated by gravel. In samplings conducted over a 600 m² area, it was determined that the species composition was largely dominated by G. turcica. Of the ~200 individuals captured using electrofishing, 90%–95% (~180–190 individuals) of G. turcica were captured. Following this were the non-native species A. escherichi and the endemic species C. antalyensis. In seine net sampling, which showed lower efficiency compared with electrofishing, the majority of individuals caught were again G. turcica. This difference in efficiency is thought to be due to the relatively slow drag speed of the seine net, which allows fish to escape. Other species in the community were found in lower abundance. The non-native species P. parva (3 individuals), G. holbrooki (1 individual), and C. gibelio (3 individuals), as well as the endemic species C. antalyensis, A. baliki, P. alii, S. burcuae, and the native species A. anguilla, were all recorded as rare (<5 individuals). In summary, G. turcica exhibited higher relative abundance in electrofishing samples than native, endemic, and other non-native species in the surveyed area. However, this numerical predominance alone does not constitute evidence of ecological displacement. The rapid establishment of a dense population indicates a strong capacity for successful establishment within the system and warrants further investigation into its potential ecological interactions. In the future, such competition could potentially put pressure on endemic species by limiting access to food sources, affecting growth and reproductive performance, and altering habitat use. If sustained over time, these interactions could contribute to reductions in local population abundance, but the long-term demographic consequences have not yet been demonstrated. G. turcica exhibits life history characteristics such as broad temperature tolerance, early sexual maturation, and the presence of multiple age classes; these characteristics have been associated with ecological resilience and successful establishment in other systems. Such features may facilitate population persistence under favorable environmental conditions and promote potential spread. While the present study does not directly assess the invasive risk outside of Ilıca Stream, these attributes may necessitate a more detailed investigation of the species’ ecological performance in the neighboring Mediterranean basins. It is clear that commercial activities and indiscriminate releases under the guise of “doctor fish” are the primary drivers of this invasion (Güçlü et al., 2025). Therefore, strict regulation of trade and transport of this species in fisheries regulations, awareness campaigns for the tourism sector, and the development of early detection and rapid response protocols in potential distribution areas are urgently needed. Such measures are considered the most effective and cost-effective strategies for invasive species management (Simberloff et al., 2013).

Figure 9

Endemic Ichthyofauna of the Ilıca Stream; (A) S. burcuae, (B) A. baliki, (C) P. alii, (D) C. antalyensis.