The biological sampling was conducted monthly from March 2017 to February 2018 in Finike Bay in the north-western Levant Basin (Fig. 2). Scientific fishing operations were conducted using longlines six to seven times per month to collect fish for the investigation. Additionally, extra pufferfish caught by commercial gill net and recreational handline fishers in the region were also collected. Each sampling day, the samples were transferred to the laboratory using coolers. First the total length and total weight were initially recorded, then individuals were dissected, the gonad weights and liver weights were measured, sex was determined, and the complete digestive tracts were removed and fixed with 70% alcohol for further stomach content analyses. For the biological analysis, a total of 1310 pufferfish were collected, including four species: 751 L. sceleratus, 325 L. suezensis, 18 L. guentheri and 216 Torquigener hypselogeneion (syn. T. flavimaculosus). Relative abundances of the four pufferfish species in total catches were 57.3% for L. sceleratus, 24.8% for L. suezensis, 16.4% for T. hypselogeneion and 1.7% for L. guentheri. Length/weight, age, growth, mortality, condition factor, reproduction and diet were investigated for the biological studies.

Total fish lengths (TL) were measured to the nearest 0.1 cm and total weights (TW) were recorded to the nearest 0.01 gram using a digital scale. The relationship between TL and TW was calculated using the equation TW=aTL ^b . Parameters (a) and (b) were estimated by applying the ordinary least squares method with natural logarithms of TW and TL. Next, to determine fish shape with growth, 95% confidence intervals of b were calculated to test whether it was significantly different from 3, as fish become thinner as they grow if b<3 and plumper if b>3 (Pauly 1984). Then, the b values of the males and females were compared using ANCOVA to determine whether there was a difference between the allometry coefficient for both sexes.

A sclerochronologic method was employed to determine the age of the specimens (Lewis 1949). First, the 2^nd to 7^th vertebrae were removed from the fish by ventral dissection. Samples were taken by scalpel, slicing a portion of the vertebrae beginning above the pectoral arch and separating large flesh portions. Residual flesh was removed by boiling the samples in hot water for three to five minutes and cleaning them with a brush. Finally, the samples were washed with distilled water and air-dried. Age was determined by counting the annual age rings on the vertebrae, which are mostly clearly seen under a stereomicroscope. The vertebrae that were unclear were sliced vertically into two pieces and examined under a stereomicroscope. Figure 3 provides an example of the annual growth rings from a vertebra sample from this study.

Length at age data was fitted to the von Bertalanffy growth function (VBGF); L_t = L_∞ (1-e^-K(t-t0)), where L_∞ is the asymptotic fish length if the fish were to grow indefinitely, K is the growth coefficient used in the VBGF expressing the rate (1/year) at which the asymptotic length is approached (Pauly 1984) and t₀ is the hypothetical age of the fish when the length was equal to zero (Sparre and Venema 1998). Growth parameters and 95% confidence intervals were estimated separately for males, females and all individuals by using the non-linear least squares method using TropFishR library (Mildenberger et al. 2017) in R statistical language (R Core Team 2021). The difference between the growth parameters of males and females was investigated using confidence intervals. Weight-specific growth parameters, Wt = W_∞ (1-e^-K(t-to))^b, were then calculated using the methods detailed in Sparre and Venema (1998). Here W_∞ was calculated as L_∞ ^b.

Total mortality rate (Z) was calculated following Beverton and Holt (1993), which estimates Z from mean fish lengths as Z = K∗(L_∞ -`L / `L - L’). L was the average fish length, and L’ the smallest fish length (Beverton and Holt 1993). Natural mortality was estimated using “Composite M: method weighting” from The Barefoot Ecologist’s Toolbox, which allows weighting of the contribution of different methods. This is an online tool compiling primary life-history-based predictors and weighting them according to input data sets. Since more than one estimator may use the same life-history data, the results are divided by the tool into the number of estimators involved. For example, there were four estimators based on maximum age, so each was weighted by 0.25 as the composite value. The tool then compiles all estimates to make a composite prior density using the empirical cumulative distribution of the point estimates. Equations and details of the composite M method can be found on the toolbox website (http://barefootecologist.com.au/shiny_m.html). Afterwards, the equation Z = M + F was used to estimate the fishing mortality rate (F). Then, the exploitation ratio was calculated from E = F/Z.

After length and weight measurements were taken, the livers and gonads were weighed to the nearest 0.01 gram. The somatic condition factor was determined using the equation SCF= (TW-GW)/TL³, where GW is gonad weight. The hepatosomatic index (HSI) is used to show the stored liver energy level of individuals, which is later transferred for reproduction purposes. It is calculated using the equation HSI = (LW/TW) 100, where LW is liver weight. These indices were calculated separately for each sex.

To determine sex, the gonads of each specimen were examined macroscopically to identify mature specimens and stereomicroscopically to identify immature specimens. Gonads containing granular structures and plump-looking tube-shaped gonads were identified as females, and those containing a milky-white substance were identified as males.

The monthly average gonadosomatic index (GSI) values were used to understand the peak reproductive period, which coincides with high GSI values. The GSI values were calculated using the equation GSI: GW/(TW-GW)100 (Gibson and Ezzi 1978). GSI provides insights into changes in ovary size and reproductive activity through gonadal development (Lowerre-Barbieri et al. 2011). Determination of reproductive timing is recommended by assessing a combination of macroscopic-microscopic evaluation and histological analysis, but high GSI is also a useful indicator and a less-costly method than others (Lowerre-Barbieri et al. 2011).

Length at first maturity (L₅₀) is defined here as the size group at which 50% of the specimens are sexually mature. To calculate L₅₀, a binomial regression model was fitted with total lengths and sex interaction. Since sex-TL interaction was not significant, sex was dropped from the model and L₅₀ was calculated for combined sexes. 95% confidence limits of L₅₀ were calculated employing 10000 bootstrap resampling.

Entire digestive tracts were removed from the samples. Then, total weights of a particular food category were measured using a digital scale (0.01 g) after absorbent paper was used to soak up extra moisture for one minute. Next, the prey items in the stomach content were examined macroscopically and identified to the lowest possible taxonomic level. To analyse seasonal variations of feeding intensity, the degree of stomach fullness was assigned using an ordinal scale between 0 and 4, where 0 = empty, 1 = 1/4 full, 2 = 1/2 full, 3 = 3/4 full and 4 = full.

Numerous numerical indices have been developed to explain the importance of different food groups. Here the following indices were used: A) frequency of occurrence (%FO) from Berg (1979), B) contribution of prey items in numbers (%N) by Hyslop (1980) and C) contribution of prey items by weight (%W) from Cortés (1997). Please note that %FO is the most robust method for prey content analysis due to fundamental problems resulting from gravimetric approaches (Becker et al. 2013).

Finally, after the above measures were derived, the index of relative importance (IRI) was calculated for each prey category using the equation IRI_i = (%N_i + %W_i) ∗ %FO_i from Hacunda (1981) to determine the main food groups. Then, percentage IRI values of each prey category were calculated with pIRI_i=IRI_i/ΣIRI.

The trophic level of L. sceleratus was then calculated based on the diet composition expressed as percentage pIRI values of each prey category using the following equation;

Here trop_i represents the trophic level of each prey item (i=1,2,..,n), which were derived from Fishbase or Sealifebase, if available (Froese and Pauly 2022). The average values of family, order or group were inputted if no information existed about the taxon’s trophic level. A bootstrap resampling procedure (9999) was performed to calculated 95% percentile confidence intervals for trophic levels using the R library boot (Canty and Ripley 2021).

A longline experiment was used for the improved catch efficiency study. Two longline sets composed of 100 hooks each were used to examine CPUE differences by season and site. The main body of the longline was a 100-gauge line with 9 m distance between branch lines, each branch line was 90 cm in length, and a straight #9/0 fishing hook was attached to each branch line.

A second experiment was performed using the standard longline type and two types of modified longlines to try to improve L. sceleratus catch efficiencies to optimize targeted removals. The first longline gear type was the standard type as described above, using a 1 mm monofilament fishing line, which is commonly used in longline fishing, as the lead line. The catches of this standard longline were compared with two types of modified longlines to compare catch efficiencies:

Then, the three gear types were trialled in fishing experiments and the results were compared for catch efficiencies. The catch (in kg) and CPUE (kg/100 hooks) were recorded along with the number of missing hooks (missing hooks/100 hooks). Number of missing hooks by different longline types were analysed using a one-way analysis of variance (ANOVA). As variances were not homogenous in analysing CPUE values, a Welch correction was applied. Pairwise comparisons were performed with Bonferroni correction. R package “one-way tests” were used in these analyses (Dag et al. 2018).

Experimental longline fishing operations were carried out specifically at the 11 main fishing sites in Finike. A total of 89 standard longline fishing operations were performed over the year-long study to test CPUE. Additionally, three longline gear type modifications were tested using ten trials each to determine whether the modifications could improve catch efficiencies for pufferfish. All longline operations were carried out between 7:00 and 14:00 during the day, with soak times ranging between two and three hours. In the study, the effect of three different snood types on the CPUE was examined. For this purpose, three different longlines, whose snood characteristics are given below, were prepared. A total of 100 snoods were used in each of these longlines. All technical features except the snood type were kept the same in the longlines (main line thickness, 1 mm; distance between snoods, 9 m; snood length, 90 cm; fishing hook No, 9/0). With these three longlines, ten operations were performed in March. CPUE (kg/hook) was calculated for each operation using the weights of individuals obtained in these ten operations. Then, it was statistically tested whether the CPUE obtained for different snood types differed. Also, the number of lost (broken) snoods in each operation was recorded. Thus, it was tested whether there was a difference in the number of lost snoods per operation according to snood types.

For each longline fishing operation, the CPUE index was used as the number of caught pufferfish per 100 hooks from each longline vessel. The CPUEs from the different fishing sites were compared using the Kruskal-Wallis test. Monthly changes in CPUE values were also compared using the Kruskal-Wallis test. The seasons were organized by assigning the last month of each season to the subsequent seasonal group (1^st group, November-December-January; 2^nd group, February-March-April; 3^rd group, May-June-July; 4^th Group, August-September-October). Of the 89 operations carried out during the 12 months, 12 were carried out in the first group, 27 in the second, 25 in the third and 25 in the fourth. The Tamhane and Dunnett T3 statistical post-hoc tests were used to examine variations by season.

Of the 751 L. sceleratus samples, 414 were male, 328 female and 8 unknown (juveniles). The M/F ratio was 1.25:1, which was significantly different from 1:1 (χ²= 9.97, p <0.01). Total length of specimens ranged from 10.2 cm to 73.4 cm. The median lengths (Wilcoxon W= 65420, p=0.39) and weights (Wilcoxon W=65091, p=0.33) of males and females were not significantly different. The mean length of individuals was 27.30±11.40, 32.34±14.90 and 49.65±8.88 cm (±standard deviation) for set nets, handlines and longlines, respectively. Based on a Kruskal-Wallis test and pairwise comparison tests with Bonferroni correction, longliners caught significantly larger specimens than set nets and handlines (p<0.05). This is because the #9/0 size hook used in this study is generally used to target medium to large gamefish, as smaller fish cannot bite larger hooks, so a standard size selectivity was applied when the hook sizes were chosen.

From Figure 4, which shows the length-frequency distribution of the samples caught by different fishing gear methods, two distinct peaks are shown, one 20 cm and the second one around 48 cm length. Clearly, the 20 cm peak is from samples caught by handlines and set nets, and the 48 cm peak is from samples caught by longlines.

Based on the ANCOVA results, the length-weight relationship parameters of each sex were not found to be significantly different. The coefficients of allometry were not significantly different from 3 in males and females, suggesting that this species retains the same body shape as it grows, but a very weak trend was detected when both sexes were combined, suggesting that some get thinner (Table 1). Somatic condition factor gradually began to decrease from April to July, coinciding with the spawning season. Except for t₀, males and females had similar von Bertalanffy growth parameters based on 95% confidence intervals (Table 1, Fig. 5). Natural mortality estimations ranged between 0.165 (Jensen method) and 0.703 (Roff method), with a mean score of 0.34 using the weighted values of both methods.

A total of nine year classes were found in this study from the vertebrae readings, ranging from 0 to 8 years. The frequency of year classes were 0.6% from the 0-age class, 17% in the first year, 15% in the second year, 10% in the third year, 26% in the fourth year (the dominant age class), 19% in the fifth year, 7% in the sixth year, 5% in the 7^th year and only 0.4% in the eighth year. The frequency of older individuals present in the samples decreased after four years in age. Longlines caught much larger fish on average (see Fig. 4), while handlines and set nets caught mostly immature ones. This is likely due to longliners targeting larger fish, hence using larger hooks able to catch larger specimens.

GSI values started to increase in March while somatic (SCF) and hepatic condition (HSI) and feeding index (FI) started to decrease at the same time (Fig. 6). There were strong negative correlations between GSI and HSI (Pearson’s r=-0.64, p<0.01) and SCF (Pearson’s r=-0.81, p<0.01). In addition, a week negative correlation with FI (Pearson’s r=-0.36, p=0.08) was found. These negative correlations are strong arguments that L. sceleratus is a capital breeder. GSI peaked in July and then rapidly declined until September, suggesting that the main reproductive period in Finike Bay is between June and August, peaking in July (Fig. 6). Total length-maturity relationship parameters were not significantly different between sexes, and length at first maturity (L₅₀) was calculated as 41.86 cm (39.73-44.08 cm; 95% confidence intervals) for combined sexes, corresponding to an age of 3.5 years (Table 1).

Of 751 individuals, 287 (38%) underwent stomach content analysis. A total of 39 different species belonging to 29 families, including 15 bony fishes, 7 crustaceans, 6 molluscs and 1 echinoderm were found as prey items (Table 2). The most important prey group was crustaceans (54%), followed by bony fish (44%). The proportion of molluscs and echinoderms was negligible compared with crustaceans and fish. There were considerable seasonal differences in prey content in terms of %IRI; crustaceans had a higher value (83%) than bony fish (15%) in spring, but this pattern reversed in summer months when more bony fish (73%) than crustaceans (19%) were found. The %IRI difference between fish and crustaceans was reduced in autumn and almost stable in winter (Table 3). The bony fish with the highest IRI were represented by the Clupeidae, Mullidae, Nemipteridae and Tetraodontidae families. The two most common decapod families were Penaeidae and Portunidae.

The stomach fullness rate was highest in March, gradually decreased to its lowest rate in July, and then increased until September. This indicates that the fish entered an intense feeding period before the spawning period; the feeding intensity decreased with the development of the gonads and intensified again after the spawning period (Fig. 6).

Table 3 shows that there is a clearly distinguishable change in either seasonal prey preferences or availability. Crustaceans were the dominant prey group in spring (March-May), but bony fish were dominant in summer, with fish and crustaceans having similar contributions in winter. Molluscs played a much more minimal role in L. sceleratus diets than fish and crustaceans.

Regarding the target fish groups, interestingly, L. sceleratus was found to prey on both pelagic (Clupeidae) and benthic families (Mullidae; Table 2). Among the fish species which could be classified to species level (n=112), predation on other pufferfish occurred in 29 cases, or 26% of the identified consumed fish prey, while cannibalism of other L. sceleratus was only found in two cases.

L. sceleratus was found to prey on 39 different families/species, clearly showing that it is an opportunistic carnivorous predator with a very wide ecological niche. The trophic level of L. sceleratus was 4.41, demonstrating that it has assumed the role of a top predator carnivore in its invaded Mediterranean Sea region.

Although there was no significant difference between the number of broken lines between the three types of longlines tested here, the highest CPUE values were attained using the second modification of the steel wire branch lines with the steel swivel hook system (p<0.05) (Table 4).

To examine spatial differences of CPUE from the 89 longline operations trialled at 11 fishing sites throughout one year, only the first fishing site had higher CPUE values (Aktaşlar with 3.6 pufferfish per 100 hooks) than the other sites, while the other ten sites had similar CPUE values (mean 2.12 pufferfish per 100 hooks; p>0.05). When the CPUE values were compared by month, catch rates were highest from November to January and May to July, and lowest from February to April and August to October (Fig. 7).