Introduction

Per-recruit analysis is one of the oldest and most widely used tools in fisheries science. It asks two fundamental questions:

  1. Yield-per-recruit (YPR): How much catch (in weight) do we get from each fish that enters the fishery?
  2. Spawning-per-recruit (SPR): How much reproductive potential is left compared to an unfished stock?

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:

  1. Load the Beverton-Holt YPR equation and compute YPR curves
  2. Build and interpret yield isopleth diagrams
  3. Understand YPR limitations and how SPR can address those
  4. Identify and interpret reference points (Fmax, F0.1, F40%)
  5. Understand the YPR–SPR trade-off
  6. Know how to apply those tools to a new species from an external dataset

1 Setup

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")))

2 The Beverton-Holt YPR Model

2.1 Example species: North Sea Cod

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")
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.

2.2 The Beverton-Holt YPR equation

The classic Beverton & Holt (1957) YPR formula:

2.3 2.2 The Beverton–Holt YPR equation

\[ 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")

3 YPR Curves and Reference Points

3.1 Compute the YPR curve

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))
)

3.2 Find Fmax and F0.1

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

3.3 Plot the YPR curve

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.

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?


4 Yield Isopleth Diagrams

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.

Yield isopleth diagram. The star marks the global YPR maximum. Increasing tc (larger mesh) shifts the optimum to lower F.

4.1 Comparing tc slices

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.

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?


5 Spawning Potential Ratio (SPR)

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)

5.1 SSB-per-recruit function

We will use a discrete age-structured model

# functions already loaded. Load again if needed
source("YPR_SPR_functions.R")

5.2 SPR curve and reference points

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.

SPR declines with increasing F. The shaded zone (SPR < 20%) indicates a depleted state.

5.3 The YPR–SPR trade-off

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.

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.


5.4 Summary table

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)")
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

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

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.)


Summary

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