Age, growth and population structure of invasive lionfish (Pterois volitans/miles) in northeast Florida using a length-based, age-structured population model

Ethics statement

All lionfish used in this study were handled in strict accordance with a UNF IACUC protocol (IACUC#13-004) and tissues of opportunity waivers approved by the University of North Florida. UNF IACUC defines tissues of opportunity as samples collected: (1) during the course of another project with an approved IACUC protocol from another institution; (2) during normal veterinary care by appropriately permitted facilities; or (3) from free-ranging animals by appropriately permitted facilities. Lionfish removals are encouraged by the State of Florida and sample collection locations did not require any specific permissions. No endangered or protected species were harmed during the course of this study.

Sample collections

Lionfish were collected from locations offshore of northeast Florida by trained spearfishermen (Fig. 1). Sampling occurred primarily during several large-scale public removal events in 2013, 2014 and 2015 (April and August) and during opportunistic sampling by recreational spearfishermen in 2014 (July, September, October, November) and early 2015 (January). All lionfish were captured offshore (>15 km) on hardbottom and artificial reef habitats in approximately 20–35 m of depth. Lionfish were measured in the field for standard (SL) and total length (TL) to the nearest 1 mm. A subset of lionfish from each tournament, and all lionfish from the opportunistic samples, were transported to the University of North Florida for further processing. In the laboratory, lionfish were measured (SL and TL), weighed to the nearest 0.1 g and sexed. Some small lionfish were difficult to accurately sex and were considered immature (Morris Jr, 2009). A subset of lionfish (n = 100) had their sagittal otoliths removed, which were used to determine lionfish age directly and in model validation.

Length-based, age-structured model

Lionfish TL data from 2014 were used to construct length-frequency histograms for the observed data from 0 to 450 mm (TL) using 10 mm increments (46 length bins) for each collection month separately (Fig. 2). Length-frequency data from 2013 and 2015 were not used in model fitting because of low temporal resolution; however these data were used to assess model performance (see Model validation). Because sex was not determined for many lionfish, length-frequency data were pooled. Growth and population age structure were estimated by fitting a length-based, age-structured model to the observed length-frequency data.

The model uses length as a proxy for age and estimates the proportion of fish of a given size in each age class using a maximum likelihood approach by fitting a predicted length frequency distribution to the observed data (Pauly & David, 1981; Johnson, 2004). To generate the predicted length-frequency distribution, the mean size-at-age was first estimated using either the traditional (von Bertalanffy, 1938; Beverton & Holt, 1957) or seasonalized (Somers, 1988; García-Berthou et al., 2012) formulation of the von Bertalanffy growth function (VBGF) which expresses fish length as a function of age. The traditional formulation of the VBGF is given below (Eq. 1): (1)${\displaystyle L_t=L_{\infty }\left[1-e^{-\left(K\left(t-t_0\right)\right)}\right]}$

where L_t is the length of a fish at age t, L_∞ is the asymptotic maximum length, K is the Brody growth coefficient, and t₀ is the theoretical time at which a fish was length 0. The seasonalized VBGF (Eq. 2) extends the traditional VBGF to allow the growth rate to vary seasonally, and may better reflect the growth of fishes in temperate regions which experience pronounced seasonal temperature fluctuations. The seasonalized VBGF (Somers, 1988) is given in the series of equations below (Eqs. 2–4): (2)${\displaystyle L_t=L_{\infty }\left[1-e^{-\left(K\left(t-t_0\right)-S\left(t\right)+S\left(t_0\right)\right)}\right]}$ (3)${\displaystyle S\left(t\right)=\left(\frac{CK}{2\pi }\right)sin\pi \left(t-t_s\right)}$ (4)${\displaystyle S\left(t_0\right)=\left(\frac{CK}{2\pi }\right)sin\pi \left(t_0-t_s\right)}$

where L_t, L_∞, K, t₀ are the same as previously defined (Eq. 1), t_s sets the timing of seasonal growth oscillations relative to t₀ (t_s + 0.5 = winter point (t_w); the time of slowest growth), and C is the intensity of seasonal growth oscillation (0 ≤ C ≤ 1). A value of C = 0 indicates no seasonality in growth, while C = 1 indicates extreme seasonality and complete cessation of growth at t_w. Because the VBGF only estimates the mean size-at-age over time; variation in the size of individual fish within each age class was estimated by including variance in size-at-age (${\sigma }_a^2$) as a model parameter. The proportion of lionfish in each age class during each sampling month and year was also estimated within the model (P_a,t). The expected number of individuals of age a in size class i in month t (n_a,i,t) was then calculated: (5)${\displaystyle n_{a,i,t}=N_t\cdot P_{a,t}\cdot P\left(L_{lower}\le L_i\le L_{upper}\left|\right.N\left({\overline{L}}_{a,t},{\sigma }_a^2\right)\right)}$

where N_t is the total number of individuals captured in month t, P_a,t is the probability of a fish captured in month t being age a, L_lower and L_upper are the lower and upper bounds of a predicted 10 mm size bin (e.g., 230 and 240 mm), and $N\left({\overline{L}}_{a,t},{\sigma }_a\right)$ defines a normal probability density function with the mean length ${\overline{L}}_{a,t}$ of fish of age a in month t estimated from the VBGF (Eqs. 1 or 2–4), and model estimated standard deviation at age, σ_a. Because size distributions overlap across ages, the total number of expected fish of size i in month t was then calculated by summing the expected contributions to size class i from each age: (6)${\displaystyle N_{i,t}=\sum _{a=0}^3{n}_{i,a,t}.}$

Four different candidate models were compared based on all combinations of two possible model structures (1) seasonal versus non-seasonal growth and (2) fixed versus variable variance in size-at-age (Table 2). Model fit was assessed by varying model parameters to minimize the log-likelihood between observed and predicted (Eq. 6) monthly length-frequency data using the SOLVER optimization routine in Microsoft Excel (MS Excel 2013, Microsoft, Inc. Seattle, WA, USA). Akaike’s Information Criterion corrected for sample size (AICc) was used to select the best model from the candidate set and quantify the relative support of each model given the data (ω_i).

The key assumptions of all candidate models were as follows: (1) predicted length-at-age follows a normal distribution with mean ${\overline{L}}_a$ and standard deviationσ_a, (2) there are only four age classes present in the observed length-frequency samples (age 0, 1, 2, and 3; this assumption was verified by aging of sagittal otoliths from a subset of lionfish (n = 100) which found that only 8% of individuals were age three (despite non-random sampling that was biased to select larger individuals), and no individuals were age four or older (see Model Validation), and (3) lionfish recruitment is assumed to occur at a single point during the year and was estimated in the model by the parameter, t_b, the estimated birth date of an annual cohort. This simplifying assumption is surely violated to some degree (every fish is not born on the same day), yet model outputs are not significantly affected because t_b in essence represents the average birth date with the variance in size-at -age (${\sigma }_a^2$) reflecting a combination of variation in birth dates confounded with the variability in growth rates among individual fish. We further assumed that (4) because diver effort varied across time and the pattern of selectivity for lionfish of varying ages is unknown, the proportion of lionfish in each age class may not accurately reflect overall population structure. While this assumption does not directly impact our estimates of model parameters, it does limit the some of the conclusions that can be drawn from the data. As a result, no attempt is made to make quantitative inferences regarding relative changes in abundance of cohorts between years or over time (e.g., recruitment strength, natural mortality).

Model performance and sensitivity

Two types of analyses were conducted to examine the robustness of our model and associated parameter estimates: (1) a randomized grid search designed to evaluate the ability of the model to converge on a consistent solution from randomly generated sets of initial parameter values (±25% best fit values) and (2) the sensitivity of model outputs was assessed by fixing individual model parameters at ±10% of their best fit values, allowing the model to converge on a new constrained solution, and examining the resulting effect on model fit and parameter estimates.

Model validation

Direct aging of a 100 fish subsample using sagittal otolith analysis was performed to verify ages and provide external validation for model outputs. Otoliths were extracted by first making a transverse cut into the skull, and removing the otoliths from under the brain cavity. Otoliths were rinsed and stored dry in envelopes until processing. Aging analysis was conducted by the Florida Fish and Wildlife Research Institute (FWRI) following the procedures outlined in VanderKooy & Guindon-Tisdel (2003). Briefly, a singular otolith from each fish was embedded in casting resin and 500 µm sections were cut using a Buehler low-speed Isomet saw. Sections were then mounted on glass slides with histomount and viewed under reflected light with a dissecting microscope at 32× magnification. Marginal increment analysis was conducted and the distance from the most recent annulus to the otolith margin was scored 1–4. Two readers aged the otoliths independently. If the ages did not agree, the otolith was removed from further analysis (n = 7). Otolith ages were then reconciled with the model chronology using Eq. (7) which generates the actual age (A_a) at the time of the first annulus formation in years (time elapsed from the estimated date of birth (t_b) to deposition of the first annulus): (7)${\displaystyle A_a=t_w+1-t_b}$

where A_a is the actual fractional age in years, t_w is the winter point from the seasonal VBGF, and t_b is the model estimated birth date. This adjustment is valid assuming that annuli deposition is coincident with the winter point (t_w) when growth is slowest. To generate ages for older fish (1+) we simply added a year for each additional annulus present. To evaluate model performance, the VBGF from the best fit model was plotted against observed size-at-age from otolith analysis to examine the level of agreement in the two independent measures of growth.

As a further evaluation of model performance, we fit the best model to three sets of observed lionfish length-frequency data from northeast Florida that were not used in parameter estimation and thus external to the model (Fig. 3). These were not used in model fitting because available samples in those years (2013, 2015) lacked sufficient temporal resolution to evaluate seasonal growth models.

General

The total lengths of lionfish sampled (n = 2, 137) ranged from 41 to 448 mm in northeast Florida over the study period. Maximum length (448 mm) and minimum length (41 mm) were both recorded in August of 2014. Some fish were not sexed and some fish were immature, but of the fish that were sexed there were 466 females present and 727 males.

Model selection

There was almost complete support for model 1 as the best fit model (ω_i ≈ 1), which assumed seasonal variability in growth and age dependent variance in size-at-age (Table 1). Models that did not assume seasonal variability in growth or had fixed variance in size-at-age fit the data poorly and had essentially no support (ω_i ≈ 0). The best model fit observed length-frequency distributions well (Fig. 2), particularly in months with large sample sizes. The best model converged on biologically realistic values for life history parameters (Table 2) with the exception of L_∞ (see below). A grid search using randomly generated initial parameter values consistently converged on the same set of parameter values indicating a global minimum and robust solution were reached. Sensitivity analyses revealed that parameter estimates were generally robust with the exception of L_∞ and k which were inversely correlated (Table 2).

Growth

The best fit model allowed for seasonal growth (Tables 1 and 2; Figs. 2 and 4); traditional VBGF formulations generated poor fits to the data (Table 1). The estimated seasonalized VBGF from the model was: (8)${\displaystyle L_t=448\mathrm{mm}\left[1-e^{-0.47\left(t-0\right)+\left[\frac{\left(0.61\ast 0.47\right)}{2\pi }sin\pi \left(t-0.71\right)\right]-\left[\frac{0.61\ast 0.47}{2\pi }sin\pi \left(0.71\right)\right]}\right].}$

All growth parameters were freely estimated with the exception of L_∞ which was fixed at the maximum observed size of 448 mm TL, a common convention (Sammons & Maceina, 2009). Constraining L_∞ was required because relatively few old fish were captured and the length data contained little information about maximum size. Consequently, the model produced biologically unrealistic estimates of L_∞ in excess of 800 mm TL when this parameter was not constrained. This is a commonly reported problem since K and L_∞ are negatively correlated and adequate fits can be achieved over a broad range of parameter values (Shepherd, 1987). We also effectively fixed t₀ at 0 by removing it from the model; this nuisance parameter was not required since the time of birth (t_b) was independently estimated within the model (Gulland & Rosenberg, 1992; Taylor, Walters & Martell, 2005). The estimated value for C was 0.61 indicating relatively strong seasonality in growth with slowest growth occurring on February 13th (t_w + t_b = 0.12). The seasonal VBGF for lionfish predicted patterns in external (not used in model fitting) observed length-frequency data collected in 2013 and 2015 well, although some evidence for interannual variation in growth is apparent (Fig. 3) and were in close agreement with size-at-age determined from otolith analysis (Fig. 4).

Recruitment

Annual lionfish recruitment was estimated to occur as a single event in early summer (t_b = 0.41 = June 2nd). The distinct bimodal distribution of total lengths among early ages in all years is consistent with a brief period of annual recruitment (Figs. 2 and 3) and also with results from otolith aging where 94% and 75% of otoliths from lionfish collected in April and August, respectively, had the same marginal increment. Moreover, the observed marginal increment relative to annuli deposition was consistent with a birth date in summer.

Population age structure

The lionfish population in northeast Florida was relatively young. Aged otoliths ranged from age 0 to age 3 (Fig. 4) with age 0, 1, 2, 3 fish comprising 9%, 69%, 14% and 9% of the population, respectively (Table 3). Otolith analysis identified no fish greater than 3 years of age supporting the model assumption of only four age classes in the population (Figs. 2 and 4). Older lionfish are probably more common than suggested by our data since the primary goal of otolith analysis was model validation and the largest fish aged in this study was smaller (342 mm TL) than the maximum size captured (448 mm TL). However, this bias is likely small since very few lionfish were larger than 342 mm TL (6.1%) and only 0.8% were greater than 400 mm TL. Model estimated population age structure generally agreed with the results from otolith aging (Table 3, Fig. 4), but direct quantitative comparisons between otolith samples and model outputs are not possible because we selectively aged larger fish, for which age is more uncertain, for model validation. For all months, the highest proportions of fish were age 1 and age 2 (Table 3). The highest proportion of age 0 fish occurred in April just prior to the cohorts estimated first birthday (A_a = 0.91) on June 2nd (t_b = 0.41; see above) when these fish, although not fully recruited to the fishery, were large enough (TL = 151 mm) to be effectively targeted by divers (Barbour et al., 2011; Edwards, Frazer & Jacoby, 2014). Newly recruited age 0 fish were largely absent during July and August, likely because fish at this age are translucent and cryptic and too small (TL = 32–59 mm) to be captured by spearfishers. Model predicted proportions at age indicate that older fish were rare in all months with Age 3 fish comprising only 6% of the population on average (Table 3).