Marine protected areas can increase the abundance of invasive lionfish (Pterois miles)

Marine protected areas (MPAs) can protect and restore marine biodiversity and fisheries, but there are concerns that they may also benefit invasive species. The spatial and temporal colonization of invasive lionfish (Pterois miles) in the eastern Mediterranean was compared across zones with varying fishing restrictions (no fishing, recreational and commercial fishing, and commercial fishing only), and stations where targeted removal events were conducted by volunteer SCUBA divers. Lionfish density in no fishing areas was nearly double that of areas with commercial fishing only, and over four times greater than in areas where both commercial and recreational fishing were allowed. Lionfish density increased with depth, possibly due to easier human exploitation in shallow waters (0–10 m) that are accessible to recreational spearfishers. Targeted removals by volunteer divers decreased lionfish densities by over 60%, while areas without removals had a 200%–400% increase. Along with management actions, natural and ecological processes might drive lionfish densities within MPAs, and the speed with which lionfish colonized fishery‐restricted zones, emphasized the need for a more sophisticated MPA management strategy that considers invasive species impacts and dynamics in an ecosystem‐based approach.

strategy that considers invasive species impacts and dynamics in an ecosystembased approach.

K E Y W O R D S
alien species, bioinvasions, conservation, fishery, invasive species, marine ecosystems, Mediterranean Sea, non-indigenous species, Suez Canal, targeted removals

| INTRODUCTION
Marine protected areas (MPAs) can restore depleted fisheries and degraded habitats (Blampied et al., 2022;Fraschetti et al., 2013;Guidetti & Sala, 2007;Stevens et al., 2014), and are a key management tool to conserve marine biodiversity and other marine resources (Duarte et al., 2020;O'Leary et al., 2016).From 2000 to 2020, the number of designated MPAs in the Mediterranean Sea has increased from 109 to 1209, and their coverage has increased over 30-fold (data from MedPAN & SPA/ RAC, 2021).More than 100 countries have committed to designate new MPAs and aim to achieve 30% protection of marine areas by 2030 (HAC, 2022).This goal aligns with the EU Biodiversity Strategy for 2030 (EC, 2020) and the Kunming-Montreal Global Biodiversity Framework (CBD, 2022) which calls for 30% of marine areas to be effectively conserved and protected by 2030.
Species diversity and biomass in highly protected areas, where no extractive activities are allowed (no-take MPAs), can be many-fold higher than in unprotected areas (Grorud-Colvert et al., 2021;Hall et al., 2023;Sala & Giakoumi, 2018).Increased biodiversity can, in turn, increase resilience to storms and other perturbations (Davies et al., 2022;Sheehan et al., 2013Sheehan et al., , 2021)).The establishment of highly protected areas is limited due to conflicts of conservation with fisheries and other societal goals and needs (Andradi-Brown et al., 2023).Partially protected areas, while popular for balancing socioeconomic interests by allowing certain fishing methods, often yield variable and modest conservation benefits compared to fully protected zones, with positive effects mostly observed in a few selected species targeted by fishing (Hollitzer et al., 2023;Sciberras et al., 2015;Turnbull et al., 2021).Moreover, there are many examples of ineffective or counter effective MPAs due to inappropriate management, science, and enforcement (Claudet & Fraschetti, 2010;Devillers et al., 2015;Rife et al., 2013).Even when highly protected and well-managed areas are established, a key question is whether such MPAs will be resilient to uncontrolled disturbances, such as climate change (Pettersen et al., 2022) and the establishment of invasive species (Giakoumi & Pey, 2017).
In the Mediterranean Sea, MPAs are strongly affected by climate change associated impacts such as the spread of warm-water invasive alien species (IAS) that arrive through the Suez Canal (D'Amen & Azzurro, 2020;Frid et al., 2023;Giakoumi, Pey, et al., 2019).Conversely, species that are particularly vulnerable to climate change may not benefit as expected from the protections MPAs offer, creating a 'Protection Paradox' (Bates et al., 2019).The presence of IAS can diminish the impact of fishery management measures (Corrales et al., 2018).Invasive species can synergistically interact with climate change to alter species interactions and food-web dynamics, potentially leading to the collapse of sensitive species populations that MPAs are designed to protect (Chaffin et al., 2016;Corrales et al., 2018;IPBES, 2023).There are concerns that MPAs in the Mediterranean, where IAS are a chronic and expanding issue, may not function as biodiversity conservation areas but as breeding grounds for IAS (Frid et al., 2023;Galil et al., 2017;Giakoumi, Pey, et al., 2019).Impacts of IAS are expected to worsen as other human pressures intensify (Geraldi et al., 2020).Long-term and standardized monitoring is essential to inform and synchronize management actions with the evolving understanding of IAS population dynamics and impacts across the expanding network of MPAs in the Mediterranean (Kleitou, Hall-Spencer, et al., 2021).
Our research focused on the colonization and spread of invasive Pterois miles (Bennett, 1828) (hereafter, 'lionfish') in an eastern Mediterranean MPA where fishery restrictions were recently established.Lionfish are currently expanding in two separate invasions of large marine ecosystems; the Mediterranean Sea and the Western Atlantic Ocean (Kletou et al., 2016;Ulman et al., 2022).This expansion is being favored by multiple traits, including an opportunistic diet, fast growth, high fecundity, early maturity, absence of predators, and naïve native prey (Agostino et al., 2020;Mouchlianitis et al., 2021;Savva et al., 2020;Ulman et al., 2021;Zannaki et al., 2019).Lionfish are a major management concern due to their high densities and demonstrated adverse impacts on biodiversity, fisheries, and food web processes (Albins & Hixon, 2008;Ballew et al., 2016;Benkwitt, 2015;Chagaris et al., 2017Chagaris et al., , 2020;;Green et al., 2012;Lesser & Slattery, 2011;Raymond et al., 2015).Invasive lionfish have reached high densities in the Western Atlantic, even in areas with healthy predator communities such as those found within MPAs (Hackerott et al., 2013;Valdivia et al., 2014).Limited predation on lionfish is attributed to factors such as their unusual shape, cryptic coloration, venomous spines, and the ability to extend them when threatened, along with the protective mucous covering their eggs (Hixon et al., 2016;Morris et al., 2011;Ulman et al., 2021).The efficacy of native predators as control agents is further diminished by their unfamiliarity with lionfish as prey and by the decreased populations of these predators due to overfishing, highlighting the need for targeted removal efforts (Kleitou, Crocetta, et al., 2021;Ulman et al., 2021).
Our objectives were to assess the lionfish relationship with fishery restrictions, habitat, and prey-fish density by comparing lionfish densities among three MPA exclusion zones: no-take, partially no-take, and fished.Sites with fishery restrictions were also compared with adjacent sites where targeted removals of lionfish were applied to understand the potential effect of such management interventions.

| Study site
Surveys were carried out in the Kavo Gkreko Natura 2000 area (CY3000005), which was designated as a Site of Community Importance (Habitats Directive 92/43/EEC) in 2008, and therefore forms an MPA under EU law.The MPA was mainly established to protect the seagrass Posidonia oceanica (Linnaeus) Delile, 1813 (Annex I habitat type with code 1120 of the Habitat Directive).Whole-site protection was established in 2018 (Κ.Δ.Π.115/2018) with the designation of three zones: (i) zone A where fishing is not allowed (hereafter: No Fishing), (ii) zone B where only commercial fishers are allowed (hereafter: Commercial Only), and (iii) zone C (hereafter: Fishing Allowed) where all legal types of fishing gears/practices are allowed, including both recreational and commercial activities (Figure 1a).Habitat mapping data have shown that the No Fishing zone has less hard substrate but substantially more coverage by Posidonia oceanica meadows compared to Commercial Only and Fishing Allowed zones (Figure 1a).Bottom trawling is prohibited in the entire area as it was enforced by Regulation (EC) 1626/1994 for depths shallower than 50 m.In the region, commercial fishers primarily employ gillnets, trammel nets, demersal longlines, and fishing pots, while recreational fishers often use spearfishing, traps, and shore-based rods (Moutopoulos et al., 2021).Both recreational and commercial fishers are subject to the same regulations regarding the minimum sizes of captured fish.Additionally, recreational fishers are subject to a general harvest limit of 10 kg per boat per fishing trip or day.This limit applies to all species caught, including lionfish, during any 24-h period.In contrast, commercial fishers do not have specific harvest limits for lionfish.The sale of recreational fishery catches is not allowed.More details about fishing practices in the region are available in Moutopoulos et al. (2021).
Between 2019 and 2020, nine targeted lionfish removals were conducted by volunteer SCUBA divers at sites within the MPA (zone C) as part of RELIONMED research project (Kleitou, Hall-Spencer, et al., 2022;Kleitou, Rees, et al., 2021;Savva et al., 2024).Specifically, five removal events took place from April to November 2019, and four occurred between February and July 2020, with each event involving 2 to 42 divers (median = 8 divers).Detailed information on these removal events can be found in Table S1.

| Research design
Two sets of underwater visual surveys were conducted to collect representative data on lionfish populations across time and space.First, repeated surveys (hereafter, 'temporal surveys') was conducted at nine monitoring stations (27 transects).Six of the monitoring stations were located at 20 m (±2 m) depth and three at 7 m (±2 m).All monitoring stations were separated from its nearest same-depth neighbor station by at least 0.5 km and were characterized by rocky substrata.The monitoring stations (Figure 1b) were surveyed consecutively for 3 years (2018-2020) during late July or early August.Four of the monitoring stations were in the Commercial Only or No Fishing zones.Five monitoring stations were in areas where all fishers could operate (Fishing Allowed), with two of these being regularly subject to targeted removals by SCUBA divers.The nine stations were categorized as (i) Fishing Allowed, (ii) Fishing Allowed & Targeted Removals (additional targeted removals were carried out by volunteer divers), and (iii) Protected (Commercial Only or No Fishing zones).
Second, spatial density sampling (hereafter, 'spatial surveys') was conducted in June-August 2020 by randomly surveying 45 monitoring stations (135 transects).A stratified-random design was employed to cover three depth categories (0-9, 10-19, and 20-30 m) within each enclosed zone area (Figure 1c).Nine stations were established randomly on P. oceanica seagrass or rocky substrata.Specifically, 27 transects were surveyed in No Fishing zone with 19 on rocky substrata and eight on P. oceanica meadows, 54 transects were surveyed in Commercial Only zone with 30 on rocky substrata and 24 on P. oceanica, and 54 were surveyed in Fishing Allowed zone with 33 on rocky substrata, 18 on P. oceanica, and three on sandy substrata (excluded from the statistical analyses).All stations below 10 m depth were located on rocky substrata.

| Visual census surveys
For both temporal and spatial monitoring surveys, three repetitive strip transects, 25 m Â 5 m (length Â width), were surveyed at each station to monitor the fish community as described by Katsanevakis et al. (2012).The dive observer positioned a Keson Transect Tape along a 25-m diving line, establishing the delineation for each strip's dive trajectory (Figure 1b).Within each strip, the observer recorded all fish found within a 2.5-m radius from the centerline in both directions.Individuals (other than lionfish) below 13 cm were categorized as potential prey items for lionfish (Côté et al., 2013;Green et al., 2012).Species with a trophic level greater than 3.85 (sourced from FishBase) and a maximum size exceeding 80 cm (sourced from FishBase) were categorized as predators, including certain species that Ulman et al. (2021) have identified as potential predators or competitors of lionfish.
In the temporal monitoring survey, an additional observer was following a modified technique to monitor only lionfish in three repetitive 25 m Â 20 m (length Â width) transects as described by Green (2012).Lionfish-specific monitoring was conducted by an observer swimming in a zig-zag pattern (Figure 1b) who searched under overhangs, crevices, and cracks in the substratum, using a dive light as needed.This technique was found to be more suitable for accurately assessing lionfish densities compared to the conventional (25 m Â 5 m transect) technique (Green et al., 2013;Kleitou, Hall-Spencer, et al., 2022).
For each recorded individual, an estimate of its length was made in situ to the nearest centimeter by the observer.Prior to data collection, all observers underwent calibration in measurement techniques using reference objects of known sizes to ensure consistency and reliability of their length estimates.
Total densities and sizes were quantified and compared among stations and zones and illustrated using 'ggplot2' package in R (version 4.2.0).

| Lionfish density over time between management zones
Collected data from the temporal survey had a high degree of multicollinearity (variance inflation factor, VIF > 10) among the investigated predictor variables 'Location' and 'Year'.Thus, we used permutational analysis of variance (PERMANOVA) (Anderson, 2014) which is robust to multicollinearity due to its permutation procedure that disrupts the inherent sample structure.PERMANOVA partitioning was performed using a similarity matrix based on the Euclidean distance of square root transformed lionfish density.Preliminary comparisons revealed only minor differences between zone A (No Fishing) and zone B (Commercial Only), as indicated by high p-values (p > .5).Consequently, the two zones were merged into a single category named 'Protected.'This consolidation was undertaken to address the challenges associated with low sample sizes and degrees of freedom, thereby enabling more reliable comparisons among factors.To test the effects of location (three levels: Fishing Allowed zone, Fishing Allowed zone combined with Targeted Removals, and 'Protected' located in No fishing or Commercial Only zones) year (three levels: 2018, 2019, and 2020), and depth (two levels: 7 and 20), and station (nine levels, nested in location and depth) on the lionfish density of each station (three replicates each), a total of 9999 permutations were computed using residuals under a reduced model.PER-MANOVA was also employed to assess the impact of location and year, along with their interaction, on the size of lionfish across different stations (nested in location).Furthermore, the correlations between lionfish densities, prey densities, and predator densities were examined using the Kendall's tau rank correlation coefficient.

| Drivers of spatial variation in lionfish density
For data collected from the spatial surveys, lionfish counts were over-dispersed with unequal conditional mean and conditional variance with many zero-captures of lionfish counts.Thus, we fitted a zero-inflated Poisson GLMM with a logit link for the zero-inflation part and a log link for the Poisson part (Brooks et al., 2017) to model the effect of the location (three levels: No Fishing, Commercial Only, and Fishing Allowed), depth (continuous in meters), habitat type (two levels: Posidonia meadows or rock), and prey-fish density (count per 100 m 2 ) on lionfish densities (count per 100 m 2 ).Transects were included as nested random effects within stations and model fit was assessed using maximum likelihood estimation with the 'BFGS' or 'L-BFGS-B' optimizer.Twelve models were fitted including location, habitat type, prey-density, predators, and depth.After fitting the zeroinflated Poisson GLMM, we conducted pairwise comparisons to assess the significance of differences in lionfish densities across the different fishing management measures.The correlations between lionfish densities, prey densities, and predator densities were also examined using the Kendall's tau rank correlation coefficient.The impact of the location on the size of lionfish across different stations (nested in the location) was examined using PER-MANOVA with 9999 permutations of residuals under a reduced model.Akaike's Information Criterion (AIC) (Akaike, 1973) was used to rank the GLMM models, with a decrease in AIC of ≥2 considered a significant improvement (Burnham & Anderson, 2004).The model with management zone, depth, prey density, and habitat as explanatory variables was evaluated as the most suitable with the lowest AIC score.To illustrate the spatial distribution of lionfish in a map, the study area was divided into 100 Â 100 m cells.The best-fit model was applied to predict the densities of lionfish in the Posidonia oceanica meadows and rocky habitats of the entire MPA.The lionfish density map was produced in QGIS.The GLMMs were conducted in R (version 4.2.0) using the 'glmmTMB' package.Pairwise comparisons were performed using the emmeans package.Model diagnostics were checked with 'DHARMa' and 'performance' packages.114.8 ± 149.6 in 2019, and 169.6 ± 171.0 in 2020.The increase in lionfish density was more pronounced in the two zones with fishery restrictions (No Fishing or Commercial Only).Over the 3 years, those stations had an increase of 378% in lionfish densities (Figure 2a).Specifically, stations within the Commercial Only zone experienced an increase of 370%, while those in the No Fishing zone experienced a 400% increase.A smaller increase was observed at stations where fishing was allowed with an increase of 210% in lionfish densities (Figure 2a).Stations with targeted removal events had a 64% decrease in lionfish numbers.Stations at 20 m depth had 230% more lionfish compared to those at 7 m depth (Figure S1).The interaction between location and year as well as depth and year had a significant effect on lionfish densities (Table 1).Increases in lionfish densities were not correlated with prey or predator densities (Kendall's Tau, p > .05; Figure 2a).The average estimated lionfish size at stations with frequent targeted removal events decreased from 19.8 ± 1.3 cm (SE) in 2018 to 17.7 ± 1.1 cm and 17.8 ± 0.7 cm in 2019 and 2020, indicating high exploitation.On the other hand, sizes in protected stations increased from 19.7 ± 1.7 cm (SE) in 2018 to 21.8 ± 0.7 cm in 2019 and to 22.4 ± 0.8 cm in 2020.Similarly, estimated sizes in fished areas increased from 19.0 ± 1.1 cm (SE) in 2018 to 19.2 ± 0.9 in 2019 and 19.6 ± 0.8 in 2020 (Figure 2b).Despite these trends, the PER-MANOVA analysis indicated no significant interaction effect between year and location on lionfish sizes (df = 4, Pseudo-F = 1.17, p > .05).

| Spatial surveys
The No Fishing zone had significantly higher lionfish densities (Figures 3 and 4; Table 2).The average densities of lionfish in the No Fishing zone were 225.2 ± 66.8 (SE) individuals per hectare, compared to 128.9 ± 34.4 individuals in the Commercial Only, and 51.8 ± 17.0 individuals in the Fishing Allowed (Figure 3a).The average size of the lionfish individuals in the No Fishing zone was 22.0 ± 0.8 cm (SE) compared to 20.4 ± 0.4 cm in the Commercial Only and 18.8 ± 0.6 cm in the Fishing Allowed (Figure 3b).The PERMANOVA analysis did not reveal these differences to be statistically significant (df = 2, Pseudo-F = 0.64, p > .05).All zones had similar T A B L E 1 Results of PERMANOVA (Type III, partial) regarding the effects of location, year, depth, and station on lionfish density based on a permutation of residuals under a reduced model (9999 number of permutations).prey densities (Figure 3a).There was a weak negative correlation between lionfish and prey densities (Kendall's Tau = À0.16,p < .05)and a moderate positive correlation between lionfish and predator densities (Kendall's Tau = 0.2, p < .05, Figure 3a).Model selection indicated the best fit GLM was the one with habitat, depth, protection zone, and prey density, albeit most of these variables did not exhibit significant predictive power (Table S2).

Source
Habitat type (rocky substrata or P. oceanica meadows) had a marginal significant effect on lionfish densities with preference displayed for rocky substrata (Table 2).
Most lionfish detected on P. oceanica habitat were not found in the actual meadows but mostly at the eroded edges of meadows with prominent "matte" comprised of rhizomes, roots, and the sediments that fill the interstices.

| DISCUSSION
MPAs should ideally be informed by science and tailored to local fisheries, biological, and socioeconomic contexts (Hilborn & Kaiser, 2022) to provide ecosystem-based conservation across suites of habitats (Solandt et al., 2020).
Although the sites with fishery restrictions initially recorded fewer lionfish in 2018, they witnessed a disproportionate increase in lionfish numbers between 2019 and 2020, contrasting with the trends observed in fished sites.The No Fishing zone (zone A) had nearly double the lionfish density of Commercial Only zone (zone B), and over four times higher lionfish densities than nearby areas that allowed fishing (zone C/D).These findings demonstrate the potential proliferation of invasive species within fishery-restricted zones, which can potentially undermine desired outcomes (D'Amen & Azzurro, 2020; Galil, 2017).
Our results also add to growing evidence that targeted lionfish removals can mitigate their impacts in MPAs (Davis et al., 2021;Kleitou, Hall-Spencer, et al., 2022;Kleitou, Rees, et al., 2021).Nine targeted removals with volunteer divers (median number of divers was eight) helped decrease lionfish densities at selected stations by over 60%, whereas other stations without targeted removals faced 200%-400% increase in lionfish density.Although lionfish were previously considered to exhibit high site fidelity (Akins et al., 2014;Bos et al., 2018;Jud & Layman, 2012;Tamburello & Côté, 2015), recent acoustic monitoring has revealed that they can also travel much farther than expected, with individual movements ranging from 1 to 10 km over a 10-day period (Green et al., 2021).Consequently, even successful removals may only temporarily reduce populations, as favorable habitats can quickly be recolonized by lionfish from adjacent areas (Kleitou, Rees, et al., 2021).To protect MPAs effectively, invasive species control programs must be both cost-effective and tailored to the specific site characteristics and the recovery dynamics post-removal (Kleitou, Rees, et al., 2021).
Developing ecotourism packages that incorporate targeted lionfish fishing could provide a sustainable and economically viable solution to control this invasive species while also generating positive social, economic, and environmental benefits (Clements et al., 2021;Kleitou, Rees, et al., 2021;Rahman et al., 2022).The challenge of maintaining suppression of these harmful invasive species while also limiting access to and regulating harvest of key fisheries species is a challenge considered by many territories and nations invaded by lionfish (Buddo, 2012).Jurisdictions across the Western Atlantic employ a range of strategies to effectively manage invasive species within MPAs.These strategies include issuing special licenses  S3.
T A B L E 2 Zero-inflated generalized linear mixed model outputs the changes in mean lionfish population density with depth, habitat, prey density, and protection zone as fixed effects and repeated transect line incorporated as a random effect.Values in bold indicate statistical significance at p < 0.05.

Factor
Odds Note: Depth, protection zone, and prey density were included as zero-inflation effects.The habitat was initially included but then removed due to multicollinearity with the rest of the variables.Model used the Poisson family with the log link function and odds ratio represents the odds of observing a change of lionfish density per 100 m 2 area, holding all other variables constant, on the log scale.Bold indicates significant ( p < .05)variables.
and permits for supervised removals using short pole spears, largely ineffective for other fish species (Ulman et al., 2022).Additionally, specific exceptions are granted to tourism operators, allowing them to conduct controlled activities during designated times and in specified areas to enhance sustainable tourism and contribute to invasive species management (Buddo, 2012;Ulman et al., 2022).Moreover, efforts are made to allocate extra resources to natural resource management staff, enabling them to develop and execute MPA-specific control plans (Kyne et al., 2020).Decreased lionfish densities outside the restricted areas emphasize the potential contribution of fisheries in controlling IAS (Kleitou, Crocetta, et al., 2021).Fishery investment strategies could be used to motivate fishers and facilitate sustainable pressure on invasive species (Bogdanoff et al., 2021;Kleitou, Crocetta, et al., 2021).For example, investing in community capacity (awareness, knowledge, skills, and collaboration), markets (increase demand, valorization opportunities, and development of novel products such as lionfish jewelry), and development of removal tools and selective fishing gears with limited destructive impact on the environment.Market incentives and campaigns to control lionfish populations need to challenge norms and engender a shift in consumer choice (Kleitou et al., 2019) but at the same time, foster a conservation-minded approach among fishers, which prioritizes ecological sustainability over marine environmental degradation (Kleitou, Moutopoulos, et al., 2022;Quintana et al., 2023).In addition to targeted removals, commercial spearfishing for lionfish may be a potential market-based solution to control their densities while diversifying fisher livelihoods (Burgess et al., 2016;Chapman et al., 2019;Harris et al., 2023;Kleitou, Crocetta, et al., 2021).It was notable that we observed higher lionfish densities at deeper depths.Freediving spearfishing is depth limited, and even SCUBA removals cannot control deeper lionfish populations.Lionfish-specific traps could be tested in fisheryrestricted zones and depths inaccessible for spearfishing removals (Harris et al., 2020(Harris et al., , 2023)).
In a mature ecosystem where niches are occupied, increased functional redundancy should add biotic resistance and resilience to IAS (Davies et al., 2022;Sheehan et al., 2013Sheehan et al., , 2021)).Strict fishery restrictions in Kavo Gkreko were established only 2 years before this study, and exploited populations may need decades to return to baseline levels after fishing pressure was removed (Duarte et al., 2020).Other factors can also affect ecosystem health and resilience of an MPA including enforcement level, MPA size, connectivity, and the intensity of fishing outside the MPA (Edgar et al., 2014;Ferreira et al., 2022;Halpern, 2014;Lester et al., 2009;Watson et al., 2014).Illegal fishing appears to be common in the area (Moutopoulos et al., 2021).Although these are valid considerations, the speed at which lionfish colonized the No Fishing zone of the MPA suggests that additional measures are needed to reduce the rate, and support eastern Mediterranean MPAs recovery processes especially in the early years of their designation.
MPAs are systematically established and promoted in biologically rich areas (Stevenson et al., 2020) characterized by complex seascape (Zhao et al., 2020).These biological and physical conditions favor the aggregation of lionfish (Bejarano et al., 2015;Hall et al., 2023;Hunt et al., 2019).Hence, in addition to management actions, natural and ecological processes within MPAs might also contribute to the (future) proliferation of lionfish populations in MPAs.Indeed, the No Fishing zone was characterized by more extensive Posidonia oceanica and rocky substrate coverage compared to the other two zones, which could also have played a role in the observed increase in lionfish densities.
Allowing ecosystems to naturally recover without any human support was suggested as a possible, cheap, and easy solution for invasive species management (Giakoumi, Katsanevakis, et al., 2019); however, this solution should be viewed with caution, especially in MPAs designed to produce conservational benefits for the surrounding ecosystems and fish populations in the eastern Mediterranean which is overfished, highly invaded and changing rapidly due to repeated marine heat wave effects.Climate change is expected to also favor biological invasions of warm water species in the Mediterranean Sea (D'Amen & Azzurro, 2020).Invasive species can reduce ecosystem resilience (Holling, 1973), drive regime shifts, and result in negative socioeconomic effects (Chaffin et al., 2016;Sheehan et al., 2021).
Inaction on the current invasions within the expanding network of MPAs seems ill advised while control of their population densities and adaptation through market promotion is a viable option (Kleitou, Crocetta, et al., 2021;Kleitou, Hall-Spencer, et al., 2022).Measures to increase ecosystem resilience and densities of lionfish predators may help mitigate invasive species impacts (Kleitou, Rees, et al., 2021).This could include rebuilding groupers (Epinephelus spp.) and octopus (Octopus vulgaris) populations which prey on lionfish in the Mediterranean Sea (Crocetta et al., 2021;Ulman et al., 2021).However, studies of the Western Atlantic lionfish invasion suggest limited biotic resistance to the invasion (Davis, 2019;Hackerott et al., 2013;Valdivia et al., 2014).Additional management measures that facilitate targeted removals of invasive species might be necessary to help MPAs achieve their conservation objectives and at the same create market opportunities via fisheries and dive tourism.
Finally, while our study highlights the potential impacts of MPA protection on lionfish populations, it is essential to acknowledge the influence of confounding variables that were not controlled in this research.These factors, which include but are not limited to, the human pressures, habitat health and complexity, and variability in oceanographic conditions, as well as intricate ecological interactions, could significantly influence our findings (Dunham et al., 2020;Fidler et al., 2021;Hunt et al., 2019;Schilling et al., 2024).Moreover, the inherent challenges in quantifying the effects of protection status on population structure, such as understanding spatial and temporal variability, add layers of complexity to our analysis.Recognizing these limitations and the need for a comprehensive approach, we call for future multifaceted research that involve more MPAs and longitudinal studies to track invasive species population changes over time while accounting for more confounding variables, such as habitat quality and complexity, along with other important factors like fishing pressure and the intrinsic characteristics of each invasive species.

F
I G U R E 1 (a) Kavo Gkreko MPA with three fishing zones: zone A (No Fishing), zone B (Commercial Only), and zone C (Fishing Allowed) and habitats coverage.Additionally, zone D is marked to show an adjacent area where fishing is allowed under the same regulations as zone C, but it lies outside the officially designated boundaries of the MPA.(b) Location and category of the nine sampling stations, and methodology applied as part of the temporal survey.(c) Location of the 45 sampling stations established in the three protection zones, and methodology applied as part of the spatial survey.The habitats of the area were mapped as part of the project DFMR (2022).
Average lionfish density in the nine stations increased from 77.0 ± 75.6 (SE) individuals per hectare in 2018 to F I G U R E 2 Temporal survey: (a) Boxplots showing the density per hectare of lionfish, identified predators (species with a trophic level higher than 3.85 according to FishBase and a maximum length exceeding 80 cm), and prey (individuals smaller than 13 cm) observed from 2018 to 2020 across the sampling stations of the temporal monitoring survey.The central line in each box indicates the median density, the box spans the interquartile range, and the whiskers extend to the furthest non-outlier points.The y-axis is limited to 750 lionfish/predator density and 30,000 for prey density to magnify the range of most common values.Values above this range were considered outliers and are not shown.(b) Kernel density plot with the mean (±SE) illustrating the lionfish size distribution changes from 2018 to 2020.

F
I G U R E 3 Spatial survey: (a) Boxplots showing the density per hectare of lionfish, identified predators (species with a trophic level higher than 3.85 according to FishBase and a maximum length exceeding 80 cm), and prey (individuals smaller than 13 cm) at the sampling stations of the spatial monitoring survey.The central line in each box indicates the median density, the box spans the interquartile range, and the whiskers extend to the furthest nonoutlier points.(b) Kernel density plot with the mean (±SE) illustrating lionfish size differences across the three zones (zone A, No Fishing; zone B, Commercial Only; and zone C, Fishing Allowed).F I G U R E 4 Spatial distribution of lionfish densities at the three zones (zone A, No Fishing; zone B, Commercial Only; and zone C, Fishing Allowed) as predicted with the zero-inflated Poisson generalized linear mixed model for hard substrate and Posidonia oceanica habitats.Zero lionfish correspond to the habitat 'soft substratum' which was not examined in the present study.Soft substrata are characterized by low (usually zero) lionfish densities.The outputs of the model are shown in Table