We adapted a modeling framework predicting age-structured population dynamics for sex-changing species (Kindsvater et al., 2017) to address the demographic consequences of fishing for PGG and bacalao, given their life histories and the specifics of their fisheries in the tropical eastern Pacific. The model structure is very similar to the general age-structured model used in fisheries stock assessments (Mangel, 2006; Haddon, 2011), but, unlike standard demographic models, follows the fates of males in addition to females. We used the published values of life history parameters for each species from our literature review, or inferred parameters from ecological principles or related species (Table 1). Individual length and age are related by data on species' growth, allowing us to model their vital rates - maturation, sex-change, fecundity, and fishing mortality—as a function of both size and age. For each species, we used discrete versions of the von Bertalanffy growth curve, which depend on its maximum length L_∞ and growth coefficient k:

This function provides the length-at-age, given a length recruitment (a = 1), until the end of life (A_max).

For both species, we assumed all individuals were born and mature as females. It is possible that some individuals mature as males, but the evidence for this is lacking and has only been shown definitively for a few grouper species (Provost et al., 2017). We assumed that maturation is a function of length (which in turn is a function of age), where the probability an individual is mature at each size, p_mat, is a logistic curve (ogive) (Figures 1C,D). The curve depends on the steepness parameter q and the parameter L_mat, which determines the length at which half of females of that length are mature (the inflection point):

Mature females transition to the terminal male class according to their length, relative to the average length of individuals in the population. We again represent the probability of sex change p_SC with a logistic function of length, which depends on a steepness parameter ρ, and the mean length of individuals in the population, L¯. Because the mean length incorporates immature individuals, we also included an offset parameter, ΔL, which adjusts the sensitivity of the sex change function to changes in mean size:

These assumptions mean that the probability an individual changes sex varies according to its length relative to the length of other mature individuals in the population. The probability an individual transitions at smaller sizes increases if large individuals are rare, and can vary as a function of time t. While little is known about the degree of plasticity in sex change rules in general, for many species there is evidence that individuals adapt to fishing by decreasing their size at sex change (Alonzo and Mangel, 2005). We also assumed that the probability to transition from female to male never exceeds the threshold of 0.8 (Figures 1G,H), so that some large individuals remain female.

We assumed that egg production is related to female mass. As is commonly assumed in fisheries science, individual mass M is related to length with a simple cubic function M(a) = vL(a)³, where v scales the mass to the appropriate units (in this case, kilograms). We can then use another power function to relate mass to egg production E at each age, with a shape coefficient c and exponent b that can be adjusted according to measurements for each species. The number of viable eggs is further dependent on the specific fertilization probability f (where f is a number between 0 and 1):

Parameter values for Equation (4) can be found in Table 1. The parameter f is of interest because it allows us to address different rates of fertilization, representing potential scenarios of male depletion. Fertilization success is difficult to study in natural populations of sex-changing species, due to our limited knowledge of their spawning behaviors. The key unknown is how changes in the sex ratio of adults affect the chance that there are enough males, or large enough males, to fertilize all eggs produced in a spawning aggregation. The numbers and sizes of adult females and males change due to direct removal by fishing, as well as indirect effects on sex change (males could transition earlier as large individuals become more rare) and female egg production (fewer eggs are produced as large females are fished out or transition to males).

The pool, or year class, of larvae P in year t can be calculated as the sum of all individuals alive at time t, multiplied by the probability they are mature, their size-specific egg production, fertilization rates, and the probability they spawn as females or transition to males (which depends on the size structure of the spawning population).

Survival from the larval stage to “recruitment” to the age-1 year class N(1) in year t+1 depends on population density. As is customary in fish population models, we assumed a non-linear relationship between larval abundance and recruitment that was characterized by the Beverton-Holt recruitment function (Beverton and Holt, 1957; Kindsvater et al., 2017). This saturating function is characterized by the parameter α, which determines the slope near the origin, and β, which, along with α, determines the asymptote (Kindsvater et al., 2016):

Because there is very little empirical information on values of α and β, we varied the values of these parameters in sensitivity analyses, but we found they did not have a strong effect on reference points.

For simplicity, after age one, we modeled an average natural mortality rate μ that applied to all ages. While it is certain that natural mortality rates vary with size, for these long-lived species natural mortality is low overall, so we approximate it with a constant rate (Kindsvater et al., 2017). As very little is known about fishing effort (e.g., number of fishers or boats), and because fishing effort certainly varies spatially for these species, we considered a range of fishing scenarios, S. The overall intensity of fishing varied with the coefficient F_max, but was distributed across sizes (and ages) according to a selectivity function F[L(a)], which were chosen to represent different types of size-selective fishing scenarios. We focused on two scenarios, based on published information about fisheries for PGG and bacalao. First, we considered a preference for the largest individuals, known as “asymptotic” selectivity and second, a fishery with “plate-sized” selectivity, where fishers targeted only small fish, between 59 and 90 cm for PGG (Baos et al., 2015), and between 35 and 55 cm for bacalao [derived from Usseglio et al. (2016) and personal observations]. The asymptotic selectivity function was:

Note the parameter j is the steepness of curve and D determines the inflection point of the asymptotic selectivity function for each species. The parameters were chosen as follows: j = 0.05 and D = 90 for PGG and j = 0.1, D = 60 for bacalao. The selectivity functions (Figure 2) were adjusted by eye to match the distribution of sizes observed in fish markets (Usseglio et al., 2016; Castellanos-Galindo et al., 2018). By comparing the demographic structure of populations in each fishing scenario, and calculating the respective biological reference points, we are able to characterize the relative effects of fishing intensity and selectivity on population biomass and spawning potential.

We simulated the population dynamics through time, given the above functions governing births, recruitment, sex change, and mortality. The population N was initialized at an arbitrary starting size. The dynamics of each age class, in annual time steps, followed

We simulated populations until they reached a steady state, from which we were able to calculate the stable age distribution. We then “fished” our simulated population until it re-equilibrated to a new steady state. The unfished steady state (the carrying capacity) of the population emerges from the balance of births, deaths, and our assumptions about density-dependent recruitment. As very little information is known about recruitment function in tropical fisheries, an alternative is to focus on the relative effects of fishing on population biomass and larval production (spawning potential) under different fishing and fertilization scenarios (Kindsvater et al., 2016). Therefore, we calculated the spawning potential ratio (SPR) for each fishing scenario, S and rate of fertilization success, f.

The spawning potential of a cohort of females is equivalent to the expected lifetime larval production of a female in a population at equilibrium (O'Farrell and Botsford, 2005). The larval production is given in Equation (5). To calculate the spawning potential ratio, the larval production of a population that has reached a fished steady state is compared to the larval production of a population that is unfished.

We also calculated the biomass reference point, where in each fertilization scenario f, the steady-state fished biomass B_{S, f} is compared with the unfished biomass B_{0, f}. The steady-state biomass in each case is the product of the number of individuals in each age class a multiplied by their mass.

This ratio reflects the degree to which the population has been depleted by fishing, assuming constant recruitment and mortality rates, which represent average conditions over the long term. It is a useful metric that integrates changes in abundance and size structure of each sex caused by direct fishing mortality and the indirect effects of fishing on sex change.

Our literature review revealed that documented growth rates [the parameter k in Equation (1)] differed within and between our focal species. The reported growth rate for goliath grouper varied across its range. We assumed that growth rates among Atlantic and Pacific species would be similar at the same latitude. Therefore, we used the estimate of k = 0.2 from Artero et al. (2015) for AGG from French Guiana (approx. latitude 4° S) as a proxy for the growth rate of PGG from Pacific Colombia (approx. latitude 4.5° S). For bacalao, growth was estimated as k = 0.11. Goliath groupers had larger reported asymptotic lengths [L_∞ = 192 cm, Artero et al. (2015)] than bacalao L_∞ = 110 cm (Usseglio et al., 2015, 2016).

Maximum age for PGG was again estimated from studies of AGG to be approximately 37 years (Bullock et al., 1992). For bacalao, the maximum age was estimated to be 30 years (Usseglio et al., 2015). Figures 1A,B shows the resulting von Bertalanffy growth function (Equation 1) for PGG and bacalao, respectively. As these species both live longer than 30 years, natural mortality is likely to be relatively low. We assumed an average mortality rate of M = 0.1 for both species, after recruitment, which is similar to estimates used in other studies of long-lived teleosts (Kindsvater et al., 2016).

We found no explicit studies of PGG reproductive biology and spawning behavior. Hence, we used the estimated L_mat = 110 cm from Bullock et al. (1992) in the Gulf of Mexico for AGG. For bacalao, Usseglio et al. (2015) found that 50% of the female population reaches sexual maturity at approximately 65.4 cm in total length (Figures 1C,D). In terms of sex-change, the probability of females transitioning to males increases as a function of length, as seen in Equation (3). In the absence of explicit information on sex change for both PGG and bacalao [e.g., only one male found in Usseglio et al. (2016)], we manually adjusted the parameters Lc and ρ until the function approximated observations on male size structure in fish markets (Figures 1G,H).

In summary, we compiled parameters derived from empirical reports and inferred from related species to produce comprehensive descriptions of growth, maturation, egg production, and sex change as a function of age for each species (Figure 1).

We next evaluated the effects of two types of fishing selectivity (i.e., asymptotic and plate-sized selectivity) on the age-structure and sex ratio of PGG and bacalao (Figure 3), given our parameters for the life history of each species. In the unfished scenario, for both species females vastly outnumbered males, a common pattern in female-first fish populations (Allsop and West, 2004). When assessing the effects of fishing (F_max = 0.3) on the population age structure, asymptotic selectivity had a greater effect, depleting the largest and oldest individuals on the population for both species (Figures 3B,F). Targeting larger individuals of PGG significantly reduced the number of fish older than 15 years in the population (Figure 3B). As in PGG, a fishery with asymptotic selectivity also depleted the oldest individuals of the bacalao population, with fewer individuals reaching ages above 20 years old (Figure 3F). As expected, fishing that targeted large and old individuals resulted in populations characterized by younger individuals and a highly female-biased sex-ratio. All else being equal, the plate-sized selectivity scenario had less effect on the steady-state population structure for both species (Figures 3C,F). Targeting individuals only within size-specific slots (PGG, 59–90 cm TL; bacalao, 35–55 cm TL) allowed some old (large) males to persist. In both selectivity scenarios, due to plasticity in the size threshold for sex change, females transitioned to males earlier in life.

We assessed the effects of varying levels of fishing mortality (F_max) and fertilization success (f) on the spawning potential ratio (SPR) and relative biomass (B_F/B₀), for PGG and bacalao populations in each selectivity and fertilization scenario (Figure 4). We also varied overall fishing effort. Across our scenarios, targeting the largest individuals (asymptotic selectivity) had the greatest impact on the biological reference points (both SPR and B_F/B₀ were lowest) for populations of both species. Reduced fertilization success affected larval production: when fertilization rates dropped from 100 to 50% the SPR values declined by half regardless of the fishing selectivity. By contrast, the relative biomass (B_F/B₀) of both species was not affected (Figure 4), because of density-dependent compensation in larval recruitment. This result would likely change with stochastic variation in the recruitment function and warrants further investigation.

Despite the lack of species-specific life history information for PGG and limited data for bacalao, we inferred realistic biological and fishery parameters based on ecological principles and field observations of both species. We found the relative vulnerability of both of these species populations depended primarily on fishery selectivity. The PGG and bacalao populations were the most vulnerable to fishing scenarios targeting large individuals. Asymptotic selectivity led to age structure truncation, heavily skewed sex ratios, and severe declines in both reference points (SPR and relative biomass).

When we considered management scenarios that could improve the sustainability of fisheries for sex-changing species, plate-size selectivity was less detrimental (given the reference points of both species) when fishing effort was low or moderate (Figures 4B,D). This suggests that fisheries management that concentrates harvest on individuals within a limited size range could allow females to survive long enough to change sex and maintain sex ratios ensuring all eggs are fertilized. However, we still do not have sufficient information for either species to gauge what sex ratio is sufficient to maintain fertilization rates. Our model predicted that a fishery targeting smaller individuals will still lead to rapid overfishing (SPR < 0.4) if fertilization success rates are reduced by 50%, for either species (Supplementary Tables 1, 2). Furthermore, in contrast to the asymptotic selectivity scenario, the slower growing bacalao had a lower SPR and relative biomass than PGG in fishing scenarios with plate-sized selectivity. This is because slower growing fish (such as bacalao) spent more time at sizes vulnerable to fishing mortality, while fast growing fish (such as PGG) outgrow this range more quickly. This is an underappreciated source of variation in the response-of species to size-based management strategies. For example, Kindsvater et al. (2017) predicted that plate-sized selectivity that targeted sub-adults decreased the spawning potential of large sex-changing species more than other types of harvest selectivity. The differences between that analysis and the results of this paper likely arise from the particular biology of each species, including the somatic growth rate, the fishery selectivity relative to maturation size, and the details of the sex change function. The analyses presented here and in Kindsvater et al. (2017) are in agreement that management actions that preserve population size-structure are critically important for sex-changing species.

Our model predicted rapid overfishing of PGG (SPR < 0.4) when fishing selects the largest individuals of the population, even at relatively low levels of fishing mortality. This prediction is consistent with observed dynamics of goliath grouper fisheries. Similar population trends have been documented for the Atlantic goliath grouper (AGG, Epinephelus itajara) and Red Hind grouper (Epinephelus guttatus), where sport fisheries targeting the largest individuals (i.e., trophy fish) of the population contributed to sperm limitation and ultimately led to reproductive failure and population collapse (Bullock et al., 1992; Beets and Friedlander, 1999; Provost and Jensen, 2015). Across most of the range of PGG, there is little to no fisheries-dependent or independent data on population trends over time (Castellanos-Galindo et al., 2018; Erisman et al., 2018). However, in Pacific Panama and Ecuador, where large individuals are routinely harvested, it is likely that PGG have experienced severe disruption in their ability to recruit, due to low numbers of spawning males in the older cohorts.

The potential for a plate-sized fishery to be sustainable for PGG is consistent with the fact that in Pacific Colombia, where “plate-size” fish fetch a premium price and thus are selectivity targeted by fishers, larger individuals of PGG are still seen (Castellanos-Galindo et al., 2018). These larger fish potentially contribute to the pool of viable gametes, sustaining recruitment. Additional factors that are unique to this area might contribute to the health of the PGG population, including low human population density, relatively low fishing pressure, and the existence of large mangrove areas providing suitable nursery habitat for juvenile PGG (Castellanos-Galindo et al., 2018).

Our model predicted when bacalao were fished with low intensity and with asymptotic selectivity, the estimated SPR value was close to the fully exploited threshold (SPR = 0.44) at the lowest level of fishing mortality (assuming all eggs are fertilized; Supplementary Table 2). In reality, the bacalao fishery in the Galapagos Islands occurs year-round and targets a wide range of sizes, including the largest individuals (Nicolaides et al., 2002). We do not know the level of fishing effort. Usseglio et al. (2016) found a highly skewed female population, and signs of overfishing (SPR < 0.5). These findings are consistent with our predictions, and suggests that if depletion of large individuals persists, reproductive failure and collapse can occur even at low levels of fishing mortality (Alonzo and Mangel, 2004).

Currently there are no catch size limits for any finish species in the Galapagos Islands, including bacalao. In the Galapagos National Park jurisdiction, fisheries management strategies include a zonation scheme (e.g., no-take zones), gear restrictions, and a licensing system (Castrejn et al., 2014). Research has suggested the potential benefits of marine reserves for protogynous species populations (e.g, Alonzo and Mangel, 2004; Heppell et al., 2006; Molloy et al., 2008; Easter et al., 2020) if they successfully protect enough of the male population to maintain healthy sex ratios (Allsop and West, 2004). However, the time it takes to observe population recovery will depend on the demography of the fished population and whether the population is demographically closed or open [i.e., the larval pool comes from inside or outside the protected area (Easter et al., 2020)].

Evidence has shown that protogynous hermaphrodite fish species are more susceptible to overfishing than gonochoristic species (Provost and Jensen, 2015; Kindsvater et al., 2017). In a sex-changing population, fishing mortality can disproportionately affect one sex, whereas in a gonochoristic population it is expected that fishing will impact both sexes equally (Provost and Jensen, 2015). When fishing mortality affects one sex more than the other, it can induce instability in the sex ratio, growth, and fertilization success of the targeted population (Alonzo et al., 2008). Therefore, sustaining exploited protogynous populations requires specific fisheries management strategies that address the details of sex-changing behaviors (Alonzo et al., 2008; Provost and Jensen, 2015).

We presented evidence suggesting that fisheries management targeting small (plate-sized) fish would be more sustainable for PGG and bacalao than current fishing practices. However, it is essential to understand the exact sex-changing rules governing these species to inform species-specific management strategies. In particular, further investigations on the size or age at which sex change in the population occurs, natural mortality, and recruitment dynamics are needed to improve model accuracy. Our population model was sensitive to our assumptions regarding natural mortality and the likelihood of transition (size at which females change sex to males). We assumed a constant natural mortality across all ages for each species, as is common in fisheries population models, but it is likely that natural mortality varies across different age groups and between sexes. Moreover, maturation and transition probabilities are size-specific and may vary among populations and locations (Alonzo and Mangel, 2004).

Finally, complex social structures of hermaphroditic fish can affect a species' susceptibility to fishing in ways that we could not address in our model (Alonzo and Mangel, 2004; Provost and Jensen, 2015). For example, during spawning aggregation events, a single large dominant male can fertilize the eggs of multiple females. Therefore, the removal of large males during a spawning period can disrupt complex social structures, causing fertilization failure (Alonzo and Mangel, 2004). However, there is no information about the mating behaviors of PGG and bacalao, and whether their reproductive patterns would differ from other protogynous grouper species (e.g., Red Hind; Sadovy et al., 1994). We recommend further investigation into the complex social structures, mating behaviors, and age-specific changes in survival and fecundity for bacalao and PGG to provide particular fisheries management strategies that would allow the recovery of the populations.

While PGG and bacalao represent commercially and culturally valuable fish species in Pacific Colombia and the Galapagos islands, there is little empirical information on the sex-changing rules governing these populations. Our modeling framework provided theoretical insights into the relative vulnerability to overfishing of these two data-poor fisheries, using life history characteristics gathered from the literature. This modeling framework has the flexibility to be applied to other fish species (i.e., protogynous and gonochoristic fish) of the Tropical Eastern Pacific, and to a variety of fishing scenarios. This approach presents an opportunity to prioritize species and fisheries for which sustainable management could be improved in developing, data-limited fisheries. For the fisheries considered here, addressing knowledge gaps with further investigation into the sex-changing rules of these species will improve our understanding of the potential consequences of fishing on population dynamics, and inform the best management strategies for the sustainability of these valuable fisheries.