Precipitous Declines in Northern Gulf of Mexico Invasive Lionfish Populations Following the Emergence of an Ulcerative Skin Disease

Invasive Indo-Pacific lionfish Pterois volitans/miles have become well-established in many western Atlantic marine habitats and regions. However, high densities and low genetic diversity could make their populations susceptible to disease. We examined changes in northern Gulf of Mexico (nGOM) lionfish populations following the emergence of an ulcerative skin disease in August 2017, when estimated disease prevalence was as high as 40%. Ulcerated female lionfish had 9% lower relative condition compared to non-ulcerated females. Changes in lionfish size composition indicated a potential recruitment failure in early summer 2018, when the proportion of new recruits declined by >80%. Remotely operated vehicle surveys during 2016–2018 indicated lionfish population density declined in 2018 by 75% on natural reefs. The strongest declines (77–79%) in lionfish density were on high-density (>25 lionfish per 100 m2) artificial reefs, which declined to similar levels as low-density (<15 lionfish per 100 m2) artificial reefs that had prior lionfish removals. Fisheries-dependent sampling indicated lionfish commercial spearfishing landings, commercial catch per unit effort (CPUE), and lionfish tournament CPUE also declined approximately 50% in 2018. Collectively, these results provide correlative evidence for density-dependent epizootic population control, have implications for managing lionfish and impacted native species, and improve our understanding of biological invasions.

Relative condition was lower in ulcerated lionfish. Sex-specific lionfish relative condition factor (K n ) was estimated for lionfish sampled in 2017 as the ratio of an individual's measured weight to its sex-specific predicted weight-at-length (  W). These  W allometric parameters were computed from 282 males and 316 females sampled from nGOM NRs prior to disease emergence during 2013-2016. ANOVA results indicated a significant difference between male and female scaling b values (F (1, 594) = 9.14; P = 0.003). Thus, female  W values were calculated with a = 4.67e-10 and b = 3.62 and male  W with a = 1.66e-10 and b = 3.42. K n was then evaluated for 338 lionfish sampled from nGOM NRs in October-November 2017, of which 39 lionfish (20 males and 19 females) presented active skin ulcers. Sampled lionfish size in total length (TL) ranged from 118-367 mm. Ulcerated lionfish TL ranged from 212-352 mm and showed no apparent trend in ulcer occurrence by size. A two-way analysis of variance (ANOVA) and Tukey honest significant differences (HSD) test was used to test the effect of ulcer presence, sex, and their interaction on differences in K n.
ANOVA results indicated K n was significantly lower for ulcerated lionfish compared to non-ulcerated lionfish (F (1, 334) = 7.81; P = 0.006). Mean K n was 0.93 (95% CI = 0.88-0.97) in ulcerated males and 0.90 (CI = 0. 84-0.95) in ulcerated females (Fig. 4, hatched bars). Multiple comparison tests indicated this significant decline was driven by the 8.8% lower K n in ulcerated females (Tukey HSD, P = 0.025). No other pairwise comparisons were significant in the Tukey HSD test. Mean K n for apparently healthy (non-ulcerated) lionfish collected in 2017 was 0.96 (CI = 0.94-0.98) in males and 0.99 (CI = 0.96-1.00) in females (Fig. 4, empty bars). ANOVA results indicated K n was not significantly different based on sex (F (1, 334) = 2.23; P = 0.140) nor the interaction between ulcer prevalence and sex (F (1, 334) = 1.46; P = 0.230). It is notable that K n in apparently healthy lionfish was 4% lower in males and 1% lower in females compared to what would be predicted by the sex-specific weight-at-length models. Potential explanations for this difference include seasonal changes in gonadosomatic index 17 or density-dependent competition for higher density populations in 2013-2016 19 . Also, lionfish collections prior to the disease were from NR habitat like the 2017 collections, but from a broader area (Fig. 1).

Size composition changed and the proportion of new recruits declined in 2018.
Northern GOM lionfish size composition was assessed for 31,073 lionfish harvested via spearfishing from 201 sampling events (collection days) during 2014-2018. Length frequency density plots of lionfish collected during 2014-2017 show bimodal distributions with local minima c. 175-200 mm, similar to the estimated demarcation between newly recruited (age-0) and older lionfish (Fig. 5A,C, solid plots). In contrast, density plots for lionfish collected in early summer 2018 show a relatively unimodal distribution (Fig. 5A,C, hatched plots). Changes in size composition and recruitment were examined by differences in the mean proportion of age-0 lionfish (P) among monthly population samples. Differences in P were evaluated with a nested binomial generalized linear model (GLM) with sampling events as the experimental unit to test differences 1) between months within the same year and 2) between years within the same month (Table 1).
GLM results indicated P during May 2014-2017 was 0.17-0.22 (Fig. 5B) and 0.23-0.28 during May-Oct 2017 ( Fig. 5D empty bars). In 2014, P decreased throughout the summer to be significantly different (P < 0.017) in October compared to May that same year, presumably as the 2014 age-0 cohort grew larger and surpassed the cutoff size. In contrast to 2014-2017, P during May-Aug 2018 declined to 0.05-0.09 (Fig. 5D). GLM results indicate P was 0.06 in May 2018 and 84% lower compared to May 2014 (Table 1, P < 0.001). Model results also indi- www.nature.com/scientificreports www.nature.com/scientificreports/ cated mean P was 0.05-0.09 during June-August 2018 and 70-86% lower than corresponding months in 2014 ( Fig. 5D hatched vs. empty bars). These differences were not statistically significant in the GLM due to the low sample size using sampling events (n = 5, 2, and 3 respectively) per month (Table 1). However, the size of the effect along with number of lionfish sampled (n = 662, 417, and 923) suggest the difference may be biologically important. Mean P increased throughout 2018, indicating a potential return of new recruits later that year. By October, P was 0.17 and similar to what it was in October 2014 when P was 0.12 (P = 0.629).
Lionfish population density declines observed by ROV sampling. Lionfish population densities were estimated during 2016-2018 via ROV surveys (n = 338 sampling events) at 18 Okaloosa ARs, 27 Escambia ARs, and 15 nGOM NRs (Fig. 1). ARs included low-density (<15 lionfish per 100 m 2 ) reefs that had experimental  lionfish removals conducted in 2014 (Escambia) 22 and 2016 (Okaloosa) 51 . High-density (>25 lionfish per 100 m 2 ) ARs were control reefs during those experiments with no lionfish removed. Changes in lionfish densities over time (by sampling date) were evaluated with generalized linear mixed models (GLMMs) with repeated sites incorporated as a random effect. GLMMs were computed to test for differences in population densities within and between density treatments (i.e., low-versus high-density) over time for Okaloosa and Escambia ARs and over time for nGOM NRs (Table 2).
During 2016 and early 2017, mean lionfish population density on high-density ARs was 21-27 lionfish per 100 m 2 on Okaloosa ARs and 19-31 lionfish per 100 m 2 on Escambia ARs (Fig. 6A,B). Lionfish density on high-density ARs were significantly higher (c. 2-5X) than mean densities on low-density ARs in the same area ( Table 2). Results from the GLMM analyses show lionfish density declines began in late 2017 and continued through 2018. Mean lionfish density on high-density Okaloosa ARs declined 79% between October 2016 (27 lionfish per 100 m 2 ) and October 2018 (6 lionfish per 100 m 2 ). By October 2017, mean lionfish density was no longer significantly different between initially high-and low-density Okaloosa ARs (Table 2). In fact, by October 2018 lionfish density on Okaloosa high-density ARs (6 lionfish per 100 m 2 ) was less than the lionfish density on corresponding low-density ARs (10 lionfish per 100 m 2 ), although the difference was not significant in the GLMM (P = 0.324). Similarly, lionfish densities at initially high-density Escambia ARs declined 77% from November 2017 (31 lionfish per 100 m 2 ) to December 2018 (7 lionfish per 100 m 2 ), by which time they were not significantly different from low-density ARs (5 lionfish per 100 m 2 , P = 0.208). Concurrent declines were also observed on low-density ARs, but were smaller in magnitude and proportion. Mean lionfish density on low-density Okaloosa ARs, which had relatively recent removals in October 2016, increased in density during early 2017 from 4 lionfish per 100 m 2 in February up to 15 lionfish per 100 m 2 in October 2017, then declined 23% by October 2018 (11 lionfish per 100 m 2 ). Mean density on low-density Escambia ARs (lionfish removed in 2014) declined between June 2016 (9 lionfish per 100 m 2 ) and December 2018 (6 fish per 100 m 2 ). Lionfish densities on NRs were approximately two-orders of magnitude lower than ARs, similar to what has been reported in other studies in the region 19,51,52 . Mean lionfish density on NRs (Fig. 6C) increased significantly between 2016-2017 (P = 0.048, Table 2). Similar in timing and proportion to the declines on high-density ARs, densities then declined 75% between October 2017 (0.35 lionfish per 100 m 2 ) and October 2018 (0.07 lionfish per 100 m 2 , P = 0.002).
Regional declines observed in lionfish fisheries catch per unit effort. Regional changes in lionfish populations were examined with fisheries-dependent CPUE information from nGOM commercial spearfishing landings and lionfish tournaments. Commercial lionfish spearfishing CPUE was assessed from Florida Fish and Wildlife Conservation Commission (FWC) trip ticket data. Catches were measured as kg lionfish landed with effort measured per trip. During 2014-2018, a total of 44,731 kg of lionfish were landed in the GOM from 1,529 trips. Over 90% of lionfish were landed in the Florida nGOM ports of Pensacola, Destin, Panama City, and Apalachicola (Fig. 7A). Lionfish commercial landings increased from near-0 in 2014 to almost 20,000 kg in 2017 (Fig. 7A). In 2018, landings then declined by 52% and commercial trips declined by 38% compared to the previous year (Table 3). www.nature.com/scientificreports www.nature.com/scientificreports/ Changes in commercial spearfishing landings CPUE were assessed with lognormal GLMs (log linked) to test differences by year (Table 3) (Table 3, P = 0.582).
Lionfish tournament CPUE was evaluated with 21,789 lionfish speared from 1,159 nGOM ARs and NRs during 2017-2019 (Table 3). Catches were measured as the number of lionfish and effort measured as number of reef sites spearfished. Repeated measures from returning teams were incorporated in the tournament mixed model as a random effect and included six teams from 2017 participating in 2018, four teams in 2018 participating in 2019, and two teams that participated in all three years. Lionfish tournament catches and effort increased each year from 2017-2019. In contrast, mean tournament CPUE peaked in 2017 then declined through 2019 (Table 3). GLMM results indicate CPUE declined by 44% in 2018 compared to 2017 (P < 0.001) with all six repeating teams having lower CPUE (Fig. 7C). Mean tournament CPUE declined an additional 18% in 2019 (P < 0.001) with six of seven repeating teams having lower CPUE.

Discussion
Rapid and substantial changes in invasive lionfish populations were observed in 2018 over a broad area of the nGOM continental shelf. Population densities declined >77% on initially high-density ARs in our study areas, which previously had the highest reported lionfish densities within the invaded range 19,51 . Similar levels of decline were observed on NRs, which had two-orders of magnitude lower lionfish densities but make up >99% of the region's reef habitat 53 . Concurrently, an approximate 80% drop in the proportion of age-0 lionfish occurred in early summer 2018, which would be the expected period of high recruitment for subtropical lionfish 54 . These population and recruitment reductions likely drove regional declines observed in total lionfish commercial spearfishing landings (52% decline), commercial spearfishing CPUE (48% decline), and tournament CPUE (62% decline). Collectively, these changes indicate total mortality in the region increased by orders of magnitude compared to prior years 27,55,56 and that nGOM lionfish populations may have experienced a population crash 4,5 .
Population declines and fluctuations are common in invasive species 1 , but their causes are often not understood 4 . In nGOM lionfish, the rapid emergence of the ulcerative skin disease about 1 y prior to the population declines reported here appears to implicate the disease as a causative mechanism. Given that disease etiology and its mortality rates remain unknown 50 , the disease-driver hypothesis for the declines is based on correlative evidence. It is apparent that population declines were density-dependent, suggesting disease transmission was also density-dependent. The first reports of diseased lionfish were from nGOM high-density ARs 49 and these had by far the greatest magnitude of population declines. By the end of 2018, lionfish density on initially high-density ARs was not significantly different than lionfish density on initially low-density ARs that had undergone prior lionfish removals. Moreover, disease prevalence was as high as 40% in late summer 2017 and the greatest declines in relative condition were observed in ulcerated female lionfish compared to non-ulcerated females. This high disease prevalence coincided with expected seasonal highs in nGOM lionfish reproductive output when waters are warmest. Northern GOM lionfish gonadosomatic index (i.e., proportion of gonad weight to fish weight) peaks www.nature.com/scientificreports www.nature.com/scientificreports/ in August at c. 6% 17 , suggesting reallocation of energetic resources from egg production to combat infection may have decreased fitness and egg production 57 . The rebound in the proportion of age-0 fish observed in fall 2018 to levels similar to fall 2014 may then be explained by the subsided disease prevalence observed earlier that year.
Other factors also may have caused or contributed to the lionfish population declines, including control by native predators, fisheries removals, and/or severe weather disturbances. Although quantifying indirect effects of competition by native predators on lionfish is challenging 20 , there is evidence that they compete with lionfish for space 58 and prey 27 . Invasive lionfish have also been occasionally reported in stomachs of native Caribbean predators 59,60 and grouper biomass has been shown to be inversely correlate to lionfish biomass 61 . However, subsequent studies with larger and more comprehensive sampling indicated no correlation between lionfish and native predators 32,34 . Furthermore, lionfish have not been documented in stomach samples of nGOM predators 27 . The current general consensus is that lionfish venomous spines deter predators 20,21 and inter-and intraspecific (i.e., cannibalism) predation rates are too low to control lionfish populations 20,31,33 . Thus, fishery removals are widely encouraged as the primary means for controlling lionfish densities 20,21 . High removal effort can control lionfish densities on local reefs 22,62,63 , but population models have indicated regional nGOM fishing intensity has been well below the levels estimated to cause recruitment overfishing 27 www.nature.com/scientificreports www.nature.com/scientificreports/ landings increased c. 8X between 2015 to 2017, it was estimated that an approximate 100X increase in fishery removals from 2015 levels would be necessary to achieve a fishing mortality rate that would substantially reduce regional lionfish biomass 27 . Moreover, the strong declines on high-density reefs indicate population effects from removals may be mitigated by density-dependent population compensation. Analogous compensatory effects were also observed in The Bahamas after the passage of a major hurricane resulted in substantial increases in lionfish density on reefs treated with the highest removal effort prior to the storm 65 .
Strong tropical weather systems may have affected lionfish larval transport, caused mortality, or redistributed lionfish (e.g., from high-density to low density reefs). The 2017-2018 tropical weather seasons included the passage of three hurricanes in the nGOM in 2017 (Harvey, Irma, and Nate) and, in 2018, Hurricane Michael was the strongest (category 5) storm on record to make landfall in northwest Florida. Previously, hurricanes have appeared to help lionfish spread and establish lionfish in the Western Atlantic. Larval transport to southeast Florida and The Bahamas was accelerated by hurricanes in the 2000s 66 and lionfish density on reefs monitored in The Bahamas increased 4-7X following the passage of a category 3 hurricane Irene 65 . However, biophysical models show lionfish larval transport in the nGOM differs from the Caribbean 67 . Relatively weak currents in the northeastern GOM, and its proximity to the GOM Loop Current, indicate the region is likely both a lionfish source and sink 68 . These oceanographic dynamics, which likely facilitated a rapid nGOM invasion and record high densities observed during 2010-2015 67,68 , may have been disrupted by hurricanes during 2017 and/or 2018. Notably, the passage of category 3 hurricane Irma through the nGOM in September 2017, about 1 month after the first record of the ulcerative skin disease, resulted in a cold-water upwelling that decreased bottom temperatures by >11°C (monitoring location 30.14°N, 86.60°W, K. Dahl, Univ. of Florida, unpubl. data). While the minimum bottom temperature (c. 16°C) remained within the tolerance range of lionfish 13 , we do not know the potential interaction effects from thermal stress, the disease, and other potential stressors.
If the disease caused the observed lionfish population declines, it would indicate a potential end of their release from parasitic or pathogenic controls in the western Atlantic. Results from multiple studies conducted earlier in the invasion indicated lionfish had low susceptibility to native-range harmful micro-and macro-parasites [35][36][37][38][39] , possibly because lionfish-associated bacterial communities are considerably different from that of native species 69 . Furthermore, these communities exhibit disease resistance to known fish pathogens 70 . Such release from native-range infectious organisms is common for established non-native species 71 and often contributes to invasion success for introduced species 72 . However, the colonization time theory predicts invaders will acquire new parasites or pathogens within their invaded range as time since introduction increases 71,73 . For example, European round goby Neogobius menalostomus in the Great Lakes had considerably lower parasite loads 1 y post establishment as compared to either native N. menalostomus populations or indigenous Great Lakes gobies, but had comparable parasite loads 15 y later 74 . At this time, the role of enemy release (i.e., from predators, parasites, and pathogens) followed by enemy accumulation as a determinant of invasive species success and population dynamics is under active research and debate 5 .
More time will be necessary to understand how native communities will be affected by changes in the lionfish populations or the lionfish epizootic. Recent experimental lionfish removal studies have indicated declines in lionfish density result in a mixture of positive 62,75 , negative 65 , and undetectable 22,63 effects on native prey and competitor fish biomass. Possible explanations for the lack of expected, positive results include the following: recruitment and succession dynamics operate on longer time scales than could be observed in the study period; ecological damages persist after lionfish are removed 4,76 ; the density of undetected lionfish in the system is greater than the critical density to achieve community-level benefits 40,51 ; or, ecological benefits are confounded by natural 65   www.nature.com/scientificreports www.nature.com/scientificreports/ outbreaks, severe storms, or other unforeseen causes will likely have unpredictable higher-order effects in marine consumer-resource interactions. Food web modeling of the nGOM ecosystem indicates complex trophic dynamics exist between native species and lionfish 27 due to lionfish occupying a mid-level trophic position and having a generalist and adaptable diet 11,22 . Future research on lionfish population dynamics, carrying capacity, and the effectiveness of mitigation efforts will be required to understand the interactions of lionfish removals and natural controls in achieving target lionfish densities, with the ultimate goal of benefiting native species. There is also a potential of pathogen spillback (i.e., when non-indigenous species act as new hosts for native infectious agents to cause new infections in native fauna) of the lionfish ulcerative disease to native species 77 . While the FWC Fish Kill hotline has received reports of ulcerated native fishes with scale loss and visible muscle (pers. comm., Florida Fish and Wildlife Research Institute), these records predate the lionfish skin disease and a low prevalence of ulcers on native fishes is expected from diseases already present 78 . An increase in skin ulcer prevalence was also observed in some nGOM species following the 2010 Deepwater Horizon oil spill 79 . Given the undetermined etiology for the lionfish ulcerative skin disease, the potential for spillback to native fishes is currently unknown 49,50 .
Long-term scientific monitoring will continue to be critical in our understanding of lionfish and invasive species population dynamics. Invasive species boom-bust cycles can be solitary with permanent effects (e.g., St. Matthew Island reindeer Rangifer tarandus 76 ) or recurrent 5 (e.g., Hudson river zebra mussels Dreissena polymorpha 80 ). The history of Australia's invasive rabbits Oryctolagus cuniculus provides a well-documented example of invasive species populations driven by epizootic control. Intentional introductions of the myxoma virus in the 1950s 81 and rabbit hemorrhagic disease virus in 1995 82 each decreased rabbit densities >90%. Each population decline was then followed by host-pathogen coevolution 83 and rabbit population recoveries 84 . Today, naturally recurring outbreaks of both diseases persist with considerable variation in geography, season, and mortality 85 . For nGOM lionfish, it is unclear if populations will increase to previous levels and thus complete a boom-bust cycle. We postulate that a solitary population crash is unlikely given their regional spread [18][19][20] and ability for non-local replenishment 67,86 . By fall 2018, lionfish size composition and the proportion of new recruits returned to similar levels observed in 2014, prior to the population declines. This recruitment signal could indicate the initial signs of a population rebound. A lag period would be expected 4 as new recruits grow and mature, which could explain the continued decline in 2019 tournament CPUE. A lionfish population recovery would be facilitated by the same physiological attributes that contributed to their invasion success [11][12][13][14] and opportunistic life history characteristics that facilitate rapid population growth, including early sexual maturity 16,17 and high reproductive output 16,17 . That said, invasive species outcomes are notoriously difficult to forecast 4 . A lionfish population recovery could be prevented by high fishing pressure 62,75,87 or if reef systems are now less hospitable to lionfish after their initial invasion depleted prey species 23,24,42 , restructured reef communities [22][23][24] , or altered ecosystem processes [25][26][27] .

Methods
Documenting ulcer prevalence on sampled nGOM lionfish. Following the first documented report of the ulcerative skin disease on August 5, 2017, gross examination for skin ulcers was conducted for all sampled lionfish to estimate the proportion of visually identified diseased individuals present within a given sample. Ulcers were macroscopically identified by examining both sides of a given lionfish for characteristic scale loss, necrotic tissue, or muscle visible through holes in the skin (Fig. 2). No partial ulcer development was observed; however, ulcer staging was not attempted. The proportion of ulcerated lionfish was estimated opportunistically from nGOM commercial spearfishing landings during August, October, and December 2017. No other data (e.g., length, sex) were collected from these samples. From January 2018 to June 2019, the proportion of ulcerated lionfish was estimated during monthly size composition sampling, described below.
Relative condition of ulcerated lionfish compared to apparently healthy and previously collected specimens. Lionfish (n = 338) were sampled by researchers via spearfishing in October and November 2017 from nGOM NRs approximately 20 km south of Destin, FL, USA (Fig. 1) www.nature.com/scientificreports www.nature.com/scientificreports/ followed humane sampling protocol with euthanasia via pithing the brain case, as reviewed and approved by the University of Florida's Institutional Animal Care and Use Committee (UF IACUC Protocol #201810225). Lionfish were macroscopically examined for evidence of skin ulcers, measured to nearest mm total length (TL), weighed to the nearest 0.1 g, and sexed by visual identification of internal gonads.
Predicted weight-at-length (  W) was computed for each lionfish with a non-linear length-weight model (Eq. 1) Allometric parameters a and b were estimated by fitting a log-linear length-weight equation (Eq. 2) with 598 lionfish sampled prior to the appearance of the lionfish disease. These lionfish were sampled during 2013-2016 from similar nGOM NRs (Fig. 1) and with the same methods described above 19 i Differences between the  log W i models for males and females were tested with a two-factor ANOVA (Eq. 3) 88 .
i ANOVA results found a significant interaction between log TL and sex (F (1, 594) = 9.14; P = 0.003), indicating a biologically significant difference between the male and female allometric scaling b values. Therefore, sex-specific  W i models were computed. Female  W i was calculated with a = 4.67e-10 and b = 3.62; male  W i was calculated with a = 1.66e-10 and b = 3.42.
Sex-specific relative condition factor (K n ) was then estimated for each lionfish sampled in 2017 as the ratio of an individual's measured weight (W i ) to its sex-specific predicted weight-at-length,  W i (Eq. 4) 89 , as is recommended for within a sample when comparing relative weights for fishes 90 Quantile-quantile (QQ) plots were examined to determine if errors were best fit with a normal, lognormal, Poisson, or negative binomial distribution. QQ plots indicated K n had a normal error distribution, thus differences in mean K n were tested with a two-factor ANOVA. Factors included ulcer presence (ulcerated or non-ulcerated), sex (male or female), and the interaction between ulcer presence and sex (Eq. 5). Group means were then tested with a Tukey HSD test. Analyses were conducted in R (version 3.5.1) with the base stats package. Data manipulation and visualization here and below were performed with R tidyverse packages DPLYR and GGPLOT2. See supplemental material for R code, raw data, and diagnostic plots. Lionfish (n = 31,073) were sampled via spearfishing during lionfish tournaments and during non-tournament periods (n = 201 sampling events). Harvesters were instructed to spear all lionfish observed on a given reef regardless of size, which was incentivized via prizes for total number of lionfish harvested. Lionfish were measured to the nearest mm TL. Harvesters were informed and consented to information about their catch and effort being used for research purposes. The number of sampling events per month ranged from 2 to 20, the number of lionfish collected per sampling event ranged from 31 to 2,402, and the number of lionfish collected per month ranged from 502 to 8,500. Length-frequency data were visually examined by month via density plots. Changes in size composition were quantitively evaluated with a logistic regression to evaluate changes in the proportion (P) of age-0 lionfish (N age-0 ) among the population sampled (N t ) per sampling event t, (i.e., collection day). Using sampling event as the experimental unit prevented overly sensitive significance findings due to large sample size and is recommended for use in length-frequency fisheries data 91  Sensitivity analysis indicated that changes in cutoff values between 120 and 178 mm for length at age-0 lionfish did not alter significance of 2018 P . Sex was not evaluated as lionfish <180 mm are generally not sexually identifiable via visual examination of gonad 19 .
Differences in P were tested with a binomial GLM (logit link) (Eq. 6). Factors included month (categorical with 6 levels, May-October) nested within year (categorical with 2 levels, 2014 and 2018) and weighted by number of lionfish measured per sampling event with a log link [i.e., log (N t )]. Nesting allowed hypothesis testing 1) between months within the same year and 2) between years within the same month. Models were fit using Laplace approximation and maximum likelihood 92 (Fig. 1). No lionfish removals were conducted at NR sites. Lionfish density surveys were conducted with a VideoRay Pro4 mini ROV equipped with a 570-line color camera with wide-angle (116°) lens and a GoPro Hero5 high-definition digital camera mounted to the front of the ROV. The ROV was controlled with an integrated control box via the ROV's tether and real-time ROV movement was observed on a high-resolution monitor with a live feed from the ROV's video camera. Lighting was provided by twin 20-watt high efficiency halogen lights mounted on the ROV. Digital video of ROV surveys was recorded then read in the laboratory to produce fish counts. ARs were surveyed with a point-count method that utilizes a 15-m wide cylinder with the reef at the center (survey area = 177 m 2 ) with fish counts made on opposite sides of the reef, above the reef, then inside of the reef 93 . NRs were surveyed with four orthogonal 25-m long transects flown from a common central point 94 . Transect width was calculated based on ROV height off the bottom (c. 1 m) and the 45° angle of the camera relative and used to compute survey area.
GLMMs were computed to test changes in lionfish density over time on Okaloosa ARs, Escambia ARs, and nGOM NRs. Lionfish counts (positive integers) were offset in the model by survey area (m 2 ) × 100 to calculate lionfish density per 100 m 2 . Fixed effects included sample date (treated categorically) nested within treatment group (categorical with two levels for AR GLMMs: low-or high-density) (Eq. 7).
OkaloosaARs EscambiaARs nGOMNRs 2 QQ plots indicated error for Okaloosa and Escambia GLMMs was best fit with negative binomial error distributions, which is typical for over-dispersed count data 92 . Error for the NR GLMM was best fit with a Poisson error distribution. To incorporate dependency of observations for repeated measurements of the same reef, individual reef was included as a random effect (random intercept) and assumed to be normally distributed with a mean of 0 and variance σ 2 . Laplace approximation was used to estimate likelihood and information criteria based on GLMM fitting and inference protocols 92 . Models were built in R (version 3.5.1) and with the LME4 and MASS packages. See supplemental material for R code, raw data, and QQ plots.
Lionfish fisheries catches and cpUe from commercial spearfishing landings and tournaments. Commercial lionfish spearfishing landings data were provided by trip ticket records managed by the FWC for Florida waters 95 . Trips were subset for lionfish food harvest (as opposed to lionfish collected for the aquarium trade) by spearfishing in GOM waters and excluded the Florida Keys. Commercial spearfishing CPUE was computed as kg lionfish landed (44,731 total) per commercial harvest trip (1,529 total). QQ plots indicated commercial CPUE error had a lognormal distribution. Differences in commercial CPUE by year during 2015-2018 (treated categorically) were tested with a log-linked GLM (Eq. 8). The first year of reported commercial harvest in 2014 was excluded due to the low sample size (n = 25) and high variance in catches. Models were fit using Laplace approximation and maximum likelihood 92 in R (version 3.5.1) with the LME4 and MASS packages. See supplemental material for code, raw data, and QQ plots.

Commercial
Lionfish tournament catches and effort were recorded from annual 2-day lionfish spearfishing tournaments conducted May 2017-2019 from teams that participated in nGOM waters. Tournament divers were informed and consented to information about their catch and spearfishing effort being used for research purposes. Catches (number of lionfish) and effort (number of reef sites) were recorded from 11 teams (44 divers) in 2017, 12 teams (48 divers) in 2018, and 16 teams (64 divers) in 2019. Tournament CPUE was computed as the number of lionfish www.nature.com/scientificreports www.nature.com/scientificreports/ caught (n = 21,789) per reef site (n = 1,159 total). Removal of all lionfish at a reef site was incentivized with tournament prizes for the most, largest, and smallest lionfish. Consistent catchability among teams was assumed based on high removal efficiency on low complexity nGOM reefs 51 . QQ plots indicated tournament CPUE had lognormal error, similar to commercial spearfishing CPUE. Changes in tournament CPUE during 2017-2019 (treated categorically) were tested with a log-linked GLMM (Eq. 9). Team was included as a random effect (random intercept) given some teams participated in tournaments in multiple years: six teams from 2017 participating in 2018, four teams in 2018 participated in 2019, and two teams that participated in all three years. Models were estimated with Laplace approximation and maximum likelihood 92 and built in R (version 3.5.1) with the LME4 and MASS packages. See supplemental material for code, raw data, and QQ plots.