Per-recruit analysis is one of the oldest and most widely used tools in fisheries science. It asks two fundamental questions:
These models require only life history parameters and assumptions about selectivity — no catch or abundance data. This makes them especially useful for data-limited situations and for understanding how gear regulations and fishing pressure interact.
Data file needed:
ypr_codlike_params.csv (a cod-like species used as case
study) Data file needed:
YPR_SPR_functions.R (for loading YPR and SPR funtions)
Learning objectives:
library(ggplot2)
library(dplyr)
library(tidyr)
library(viridis)
theme_set(theme_minimal(base_size = 14) +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold")))
We start with parameters to a species similar to North Sea cod — a temperate demersal species
codlike <- read.csv("ypr_codlike_params.csv")
codlike %>% knitr::kable(caption = "Life history parameters for a cod-like species from temperate waters")
| parameter | value | unit | source | notes |
|---|---|---|---|---|
| species | Cod like species | FAO | A simulated cod-like species | |
| Linf_cm | 120 | cm | FishBase | von Bertalanffy asymptotic length |
| K_per_year | 0.2 | year-1 | FishBase | von Bertalanffy growth coefficient |
| t0_years | -0.5 | years | FishBase | theoretical age at length 0 |
| M_per_year | 0.2 | year-1 | assumed | estimated from tmax using Hoenig method |
| a_WL | 0.01 | g/cm^b | FishBase | weight-length coefficient |
| b_WL | 3 | FishBase | weight-length exponent (slightly positive allometric) | |
| tr_years | 1 | years | assumed | age at recruitment to fishing grounds |
| tmax_years | 15 | years | local study | maximum observed age |
| mat50_years | 4 | years | local study | age at 50% maturity |
| mat95_years | 6 | years | local study | age at 95% maturity |
| F_current | 0.55 | year-1 | stock assessment | estimated current fishing mortality |
| tc_current | 2.5 | years | selectivity study | current age at first capture (from trawl fishery experiment) |
# Parse into a params list (with a structure that we can handle)
get_val <- function(param_name) {
codlike$value[codlike$parameter == param_name] |> as.numeric()
}
params <- list(
Linf = get_val("Linf_cm"),
K = get_val("K_per_year"),
t0 = get_val("t0_years"),
a = get_val("a_WL"),
b = get_val("b_WL"),
M = get_val("M_per_year"),
tr = get_val("tr_years"),
tmax = get_val("tmax_years"),
mat50 = get_val("mat50_years"),
mat95 = get_val("mat95_years")
)
params$Winf <- params$a * params$Linf^params$b
F_current <- get_val("F_current")
tc_current <- get_val("tc_current")
cat("Linf =", params$Linf, "cm, K =", params$K,
", M =", params$M, ", M/K =", params$M / params$K, "\n")
#> Linf = 120 cm, K = 0.2 , M = 0.2 , M/K = 1
cat("Winf =", round(params$Winf, 0), "g\n")
#> Winf = 17280 g
The M/K ratio is a key life-history invariant. For this cod-like species M/K = 1, which is at the low end of the typical range (1.0–2.0 for teleosts). Low M/K species tend to produce dome-shaped YPR curves with a clear Fmax. High M/K species (fast growth, high mortality) often produce flat curves where Fmax is very high or undefined.
The classic Beverton & Holt (1957) YPR formula:
\[ YPR = F \, e^{-M(t_c-t_r)} \, W_\infty \sum_{n=0}^{3} \frac{U_n \, e^{-nK(t_c-t_0)}}{F+M+nK} \]
Where:
| Symbol | Meaning |
|---|---|
| \(YPR\) | Yield per recruit (e.g., kg per recruit) |
| \(F\) | Fishing mortality rate (typically yr\(^{-1}\)) |
| \(M\) | Natural mortality rate (typically yr\(^{-1}\)) |
| \(t_r\) | Age at recruitment to the population / model (years) |
| \(t_c\) | Age at first capture (knife-edge selectivity) (years) |
| \(W_\infty\) | Asymptotic weight (weight at infinite age) |
| \(K\) | von Bertalanffy growth coefficient (typically yr\(^{-1}\)) |
| \(t_0\) | von Bertalanffy theoretical age at zero length (years) |
| \(n\) | Index of the expansion term (\(0,1,2,3\)) |
| \(U_n\) | Coefficients from the binomial expansion of \((1-e^{-K(t-t_0)})^3\): \(\{1,-3,3,-1\}\), which describes the weight-at-age curve when b=3 |
| (isometric growth) |
# load the Beverton-Holt YPR and SPR function we will need
source("YPR_SPR_functions.R")
# Quick test
#cat("YPR at F=0.3, tc=3:", round(beverton_holt_ypr(0.3, 3, params), 1), "g/recruit\n")
F_values <- seq(0, 1.5, by = 0.01)
tc_fixed <- 3
F_current <- 0.5
ypr_data <- data.frame(
F = F_values,
YPR = sapply(F_values, function(f) beverton_holt_ypr(f, tc_fixed, params))
)
Fmax is the fishing mortality that maximises yield-per-recruit. F0.1 is more precautionary — it is the F at which the slope of the YPR curve is 10% of the slope at the origin. At F0.1, additional effort produces very little extra yield.
# Fmax
Fmax_idx <- which.max(ypr_data$YPR)
Fmax <- ypr_data$F[Fmax_idx]
YPR_max <- ypr_data$YPR[Fmax_idx]
# F0.1
delta_F <- 0.001
slope_origin <- (beverton_holt_ypr(delta_F, tc_fixed, params) -
beverton_holt_ypr(0, tc_fixed, params)) / delta_F
target_slope <- 0.1 * slope_origin
slopes <- diff(ypr_data$YPR) / diff(ypr_data$F)
F0.1_idx <- which.min(abs(slopes - target_slope))
F0.1 <- (ypr_data$F[F0.1_idx] + ypr_data$F[F0.1_idx + 1]) / 2
YPR_F0.1 <- (ypr_data$YPR[F0.1_idx] + ypr_data$YPR[F0.1_idx + 1]) / 2
cat("Fmax =", round(Fmax, 3), " → YPR =", round(YPR_max, 1), "g/recruit\n")
#> Fmax = 0.4 → YPR = 2048.9 g/recruit
cat("F0.1 =", round(F0.1, 3), " → YPR =", round(YPR_F0.1, 1), "g/recruit\n")
#> F0.1 = 0.185 → YPR = 1861.6 g/recruit
YPR_current <- beverton_holt_ypr(F_current, tc_fixed, params)
ggplot(ypr_data, aes(F, YPR)) +
geom_line(color = "#065A82", linewidth = 1.2) +
# Current F
geom_vline(xintercept = F_current, linetype = "dotted", color = "grey40") +
annotate("text", x = F_current + 0.04, y = YPR_max * 0.15,
label = sprintf("Fcurrent = %.1f", F_current),
hjust = 0, color = "grey40", size = 3.5) +
# F0.1
geom_vline(xintercept = F0.1, linetype = "dashed", color = "#F96167") +
geom_point(aes(x = F0.1, y = YPR_F0.1), color = "#F96167", size = 4) +
annotate("text", x = F0.1 - 0.05, y = YPR_F0.1 * 0.88,
label = sprintf("F0.1 = %.2f", F0.1), hjust = 1, color = "#F96167") +
# Fmax
geom_vline(xintercept = Fmax, linetype = "dashed", color = "#00A896") +
geom_point(aes(x = Fmax, y = YPR_max), color = "#00A896", size = 4) +
annotate("text", x = Fmax + 0.05, y = YPR_max,
label = sprintf("Fmax = %.2f", Fmax), hjust = 0, color = "#00A896") +
labs(title = "Yield-Per-Recruit Curve",
subtitle = sprintf("tc = %d years, M = %.2f year⁻¹", tc_fixed, params$M),
x = expression(paste("Fishing Mortality (F, year"^{-1}*")")),
y = "Yield per Recruit (g)") +
scale_x_continuous(expand = c(0, 0), limits = c(0, 1.5)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, YPR_max * 1.1))
YPR curve for a cod-like species. F0.1 captures ~95% of the maximum yield at much lower F.
Q: Why is F0.1 always less than Fmax? What does the “10% slope” rule mean in practical terms?
Q: At Fcurrent = 0.5, is this stock above or below Fmax? Is the fishery getting close to the maximum possible yield per recruit?
Q: Notice the flat top of the curve between F0.1 and Fmax. What does this plateau tell you about the cost-effectiveness of increasing fishing effort?
The YPR curve above assumes a fixed age at first capture (tc = 3). But in reality, both F and tc are management levers. Isopleth diagrams show YPR as a surface over both dimensions.
F_range <- seq(0.01, 1.2, by = 0.02)
tc_range <- seq(1, 8, by = 0.1)
isopleth_data <- expand.grid(F = F_range, tc = tc_range)
isopleth_data$YPR <- mapply(
function(f, tc) beverton_holt_ypr(f, tc, params),
isopleth_data$F, isopleth_data$tc
)
# Find the global maximum
max_idx <- which.max(isopleth_data$YPR)
optimal_F <- isopleth_data$F[max_idx]
optimal_tc <- isopleth_data$tc[max_idx]
cat("Global maximum YPR at F =", round(optimal_F, 2),
", tc =", round(optimal_tc, 1), "years\n")
#> Global maximum YPR at F = 1.19 , tc = 5.7 years
ggplot(isopleth_data, aes(x = tc, y = F, z = YPR)) +
geom_contour_filled(bins = 15) +
geom_contour(color = "white", alpha = 0.5, bins = 10) +
geom_point(aes(x = optimal_tc, y = optimal_F),
color = "red", size = 4, shape = 8, stroke = 2,
inherit.aes = FALSE) +
annotate("text", x = optimal_tc + 0.5, y = optimal_F + 0.05,
label = sprintf("Max YPR\n(tc=%.1f, F=%.2f)", optimal_tc, optimal_F),
color = "white", fontface = "bold", size = 3.5) +
scale_fill_viridis_d(name = "YPR (g)", option = "C") +
labs(title = "Yield Isopleth Diagram",
subtitle = sprintf("M = %.2f, K = %.2f, Linf = %.0f cm",
params$M, params$K, params$Linf),
x = "Age at First Capture (tc, years)",
y = expression(paste("Fishing Mortality (F, year"^{-1}*")")))
Yield isopleth diagram. The star marks the global YPR maximum. Increasing tc (larger mesh) shifts the optimum to lower F.
To see this more clearly, we can extract YPR curves at different ages of first capture:
tc_slices <- c(2, 4, 6)
ypr_slices <- isopleth_data %>%
filter(tc %in% tc_slices) %>%
mutate(tc = factor(tc, labels = paste("tc =", tc_slices, "years")))
ggplot(ypr_slices, aes(F, YPR, color = tc)) +
geom_line(linewidth = 1.1) +
scale_color_manual(values = c("#065A82", "#00A896", "#F96167"),
name = "Age at First\nCapture") +
labs(title = "YPR Curves at Different Ages of First Capture",
x = expression(paste("Fishing Mortality (F, year"^{-1}*")")),
y = "Yield per Recruit (g)")
YPR curves at three different ages of first capture. Delaying tc increases maximum YPR and shifts Fmax lower.
Q: In the isopleth diagram, what happens to YPR if you fish very hard (high F) but catch fish very young (low tc)? Why?
Q: What does increasing tc represent in management terms? (Hint: think about mesh size or minimum landing size regulations.)
Q: Look at tc = 6. The maximum YPR is higher than at tc = 2, but it requires lower F. What is the practical implication?
SPR measures what fraction of the unfished reproductive output remains under fishing:
\[ SPR(F)=\frac{(SSB/R)_F}{(SSB/R)_{F=0}} \]
| Symbol | Meaning |
|---|---|
| \(SPR(F)\) | Spawning potential ratio at fishing mortality \(F\) |
| \(SSB\) | Spawning stock biomass |
| \(R\) | Recruitment |
| \((SSB/R)_F\) | Spawning biomass per recruit under fishing (\(F\)) |
| \((SSB/R)_{F=0}\) | Spawning biomass per recruit when unfished (\(F=0\)) |
| \(F\) | Fishing mortality rate (typically yr\(^{-1}\)) |
YPR tells you how much yield you get; SPR tells you the cost of that yield in terms of reproductive potential. Common reference points:
| Target | Interpretation |
|---|---|
| F40% | Fishing Mortality (F) where SPR = 40% (widely used as a proxy for FMSY) |
| F35% | Intermediate threshold |
| F30% | Lower limit; Values of SPR < 30% often considered risky |
| SPR < 20% | Limit / severe depletion risk (often treated as a management limit) |
We will use a discrete age-structured model
# functions already loaded. Load again if needed
source("YPR_SPR_functions.R")
ssbr_F0 <- ssb_per_recruit(0, tc_fixed, params)
spr_data <- data.frame(F = F_values) %>%
mutate(
SSBR = sapply(F, function(f) ssb_per_recruit(f, tc_fixed, params)),
SPR = SSBR / ssbr_F0,
YPR = sapply(F, function(f) beverton_holt_ypr(f, tc_fixed, params))
)
find_F_at_SPR <- function(target_spr, spr_df) {
idx <- which(diff(sign(spr_df$SPR - target_spr)) != 0)
if (length(idx) == 0) return(NA)
idx <- idx[1]
x1 <- spr_df$F[idx]; y1 <- spr_df$SPR[idx]
x2 <- spr_df$F[idx+1]; y2 <- spr_df$SPR[idx+1]
x1 + (target_spr - y1) * (x2 - x1) / (y2 - y1)
}
F40 <- find_F_at_SPR(0.40, spr_data)
F35 <- find_F_at_SPR(0.35, spr_data)
F30 <- find_F_at_SPR(0.30, spr_data)
SPR_current <- ssb_per_recruit(F_current, tc_current, params) /
ssbr_F0
cat("F40% =", round(F40, 3), "\n")
#> F40% = 0.182
cat("F35% =", round(F35, 3), "\n")
#> F35% = 0.214
cat("F30% =", round(F30, 3), "\n")
#> F30% = 0.253
cat("F0.1 =", round(F0.1, 3), " (for comparison)\n")
#> F0.1 = 0.185 (for comparison)
cat("Fmax =", round(Fmax, 3), " (for comparison)\n")
#> Fmax = 0.4 (for comparison)
cat(" SPR at Fcurrent =", round(SPR_current, 3),
sprintf("(%.0f%%)\n\n", SPR_current * 100))
#> SPR at Fcurrent = 0.108 (11%)
spr_refs <- data.frame(
F_ref = c(F40, F35, F30), SPR_ref = c(0.40, 0.35, 0.30),
label = c("F40%", "F35%", "F30%"),
color = c("#00A896", "#F5A623", "#F96167")
)
ggplot(spr_data, aes(F, SPR)) +
annotate("rect", xmin = 0, xmax = 1.5, ymin = 0, ymax = 0.20,
fill = "#FFCCCC", alpha = 0.4) +
annotate("text", x = 1.4, y = 0.10, label = "Depleted\n(SPR < 20%)",
hjust = 1, size = 3, color = "#C0392B") +
geom_line(color = "#065A82", linewidth = 1.2) +
geom_hline(yintercept = c(0.40, 0.35, 0.30),
linetype = "dashed", color = spr_refs$color, linewidth = 0.7) +
geom_vline(xintercept = c(F40, F35, F30),
linetype = "dashed", color = spr_refs$color, linewidth = 0.7) +
geom_point(data = spr_refs, aes(x = F_ref, y = SPR_ref),
color = spr_refs$color, size = 3, inherit.aes = FALSE) +
annotate("text", x = c(F40, F35, F30) + 0.05, y = c(0.43, 0.38, 0.33),
label = sprintf("%s = %.3f", spr_refs$label, spr_refs$F_ref),
hjust = 0, size = 3.3, color = spr_refs$color) +
labs(title = "Spawning Potential Ratio vs Fishing Mortality",
subtitle = sprintf("tc = %d years, mat50 = %.0f years", tc_fixed, params$mat50),
x = expression(paste("Fishing Mortality (F, year"^{-1}*")")),
y = "SPR") +
scale_x_continuous(expand = c(0, 0), limits = c(0, 1.5)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1.05),
labels = scales::percent_format(accuracy = 1))
SPR declines with increasing F. The shaded zone (SPR < 20%) indicates a depleted state.
This is one of the most important diagnostic plots in fisheries management. It shows that maximising yield (at Fmax) often comes at the cost of severely depressing reproductive output:
ref_points <- data.frame(
F = c(F0.1, F40, F35, F30, Fmax),
label = c("F0.1", "F40%", "F35%", "F30%", "Fmax"),
color = c("#E91E63", "#00A896", "#F5A623", "#F96167", "#3F51B5")
) %>%
mutate(
YPR = sapply(F, function(f) beverton_holt_ypr(f, tc_fixed, params)),
SPR = sapply(F, function(f) ssb_per_recruit(f, tc_fixed, params) / ssbr_F0)
)
ggplot(spr_data, aes(SPR, YPR)) +
annotate("rect", xmin = 0, xmax = 0.30,
ymin = 0, ymax = max(spr_data$YPR) * 1.1,
fill = "#FFE0E0", alpha = 0.5) +
annotate("text", x = 0.15, y = max(spr_data$YPR) * 1.02,
label = "SPR < 30%\n(Risk zone)", size = 3, color = "#C0392B") +
geom_line(color = "#065A82", linewidth = 1.2) +
geom_point(data = ref_points, aes(SPR, YPR),
color = ref_points$color, size = 4, inherit.aes = FALSE) +
geom_text(data = ref_points, aes(SPR, YPR, label = label),
nudge_y = max(spr_data$YPR) * 0.04,
color = ref_points$color, size = 3.5, inherit.aes = FALSE) +
geom_vline(xintercept = 0.40, linetype = "dashed", color = "#00A896",
linewidth = 0.7, alpha = 0.7) +
labs(title = "YPR vs SPR Trade-off",
subtitle = "Right = less fishing; Up = more yield per recruit",
x = "Spawning Potential Ratio (SPR)",
y = "Yield per Recruit (g)") +
scale_x_continuous(limits = c(0, 1), expand = c(0, 0),
labels = scales::percent_format(accuracy = 1)) +
scale_y_continuous(limits = c(0, max(spr_data$YPR) * 1.1), expand = c(0, 0))
The YPR–SPR trade-off. Fmax gives highest yield but pushes SPR dangerously low. F0.1 and F40% offer much better compromises.
Q: What is the SPR at Fmax? Would you consider this a safe level of reproductive output?
Q: Compare F0.1 and F40% — are they similar? Why might a manager prefer one reference point over the other?
Q: In the trade-off plot, notice that moving from F40% to Fmax gains relatively little extra yield but costs a lot of SPR. Explain this in terms you’d use when advising a manager.
tibble(
`Reference Point` = c("Fmax", "F0.1", "F40%", "F30%", "Fcurrent"),
`F (year⁻¹)` = c(Fmax, F0.1, F40, F30, F_current),
`YPR (g)` = c(YPR_max, YPR_F0.1,
beverton_holt_ypr(F40, tc_current, params),
beverton_holt_ypr(F30, tc_current, params),
YPR_current),
`SPR (%)` = c(
ssb_per_recruit(Fmax, tc_current, params) / ssbr_F0 * 100,
ssb_per_recruit(F0.1, tc_current, params) / ssbr_F0 * 100,
40, 30, SPR_current * 100
)
) %>%
mutate(across(where(is.numeric), ~round(., 2))) %>%
knitr::kable(caption = "Cod-like species (reference points summary)")
| Reference Point | F (year⁻¹) | YPR (g) | SPR (%) |
|---|---|---|---|
| Fmax | 0.40 | 2048.90 | 15.00 |
| F0.1 | 0.18 | 1861.56 | 35.96 |
| F40% | 0.18 | 1807.27 | 40.00 |
| F30% | 0.25 | 1905.92 | 30.00 |
| Fcurrent | 0.50 | 2036.47 | 10.83 |
Q: Is Fcurrent above or below F0.1? Above or below F40%? What does this tell you?
Q: Look at the isopleth diagram. Is the current operating point (red cross) near the optimal yield zone, or could a different tc improve things?
Q: At the current F, what is the SPR? Is the spawning stock at risk?
Q: Write a 3-sentence management recommendation for this stock based on all the results.
Exercise — Mesh size change for the cod-like species
The fisheries authority is considering increasing the mesh size so that tc shifts from 2.5 to 5 years. Re-compute the YPR and SPR at the current F under the new tc:
ypr_new_tc <- beverton_holt_ypr(F_current, 5, params)
spr_new_tc <- ssb_per_recruit(F_current, 5, params) / ssbr_F0
cat("Current tc:", tc_current, "→ YPR =", round(YPR_current, 1),
"g, SPR =", round(SPR_current * 100, 1), "%\n")
cat("New tc: 5.0 → YPR =", round(ypr_new_tc, 1),
"g, SPR =", round(spr_new_tc * 100, 1), "%\n")
Does yield increase? Does SPR improve? Is this a good management measure?
Class discussion — The “flat curve” problem
For some species, the the YPR curve may be quite flat (high M/K). A manager says “Since there’s barely any dome in the YPR curve, it doesn’t matter how hard we fish.” Write 3–4 sentences explaining why this conclusion is wrong. (Hint: remember that the YPR curve says nothing about recruitment, that is where SPR comes in.)
| Concept | Key message |
|---|---|
| YPR curve shape | Low M/K → clear dome with Fmax; High M/K → flat curve |
| Fmax | Maximises yield but often depletes spawning biomass |
| F0.1 | Precautionary; captures ~95% of max yield at much lower F |
| Isopleths | Both F and tc matter — mesh regulations complement effort control |
| SPR | Measures cost of fishing to reproductive output |
| F40% | Widely used proxy for FMSY; preserves 40% of unfished SSB/R |
| Trade-off | Moving from F40% to Fmax gives little extra yield at high SPR cost |