Fish Population Trends in the Pontchartrain Basin

Author

Group 6: Mariana, Aleja & Isabella

Introduction

Monitoring fish populations is critical for detecting declines, understanding the drivers of change, and informing sustainable management. Globally, fish populations have experienced long-term declines driven by overfishing, habitat loss, and degradation (Sumalia et al. 2020). In the United States, many fish populations have shown long-term declines, and 18% of federally managed stocks are currently classified as overfished (NOAA 2023). The main source of fish mortality in US coastal waters is recreational fishing (Coleman et al. 2004). Therefore, implementing long-term, standardized, and effective monitoring programs that can track population trajectories and evaluate the resilience of heavily exploited species becomes especially important for coastal estuaries.  

As a major estuarine system in the northern Gulf of Mexico, Lake Pontchartrain in Louisiana is the type of coastal environment where long-term monitoring is of importance. This lake basin has multiple freshwater sources; this mix of fresh and salt water creates nutrient-rich conditions that sustain abundant prey resources and species diversity (Louisiana Department of Wildlife & Fisheries 2022). Strong water circulation and connectivity with Lake Borgne, the Mississippi Sound, and the surrounding marshes promote the movement of baitfish and juvenile estuarine species through the system (Vasconcelos et al. 2011). Due to its high productivity, this area is considered one of the most important recreational fishing destinations in the northern Gulf of Mexico, providing millions of dollars in income through recreational fishing and tourism (Keddy et al. 2007, Louisiana Department of Wildlife & Fisheries 2022). The lake’s conditions support year-round fisheries that target both seasonal and resident species such as spotted seatrout (Cynoscion nebulosus) and red drum (Sciaenops ocellatus), with black drum (Pogonias cromis) and sheepshead (Archosargus probatocephalus, Louisiana Department of Wildlife & Fisheries 2022).  However, fluctuations in abundance have been observed for many species due to variable recruitment, fishing pressure, agricultural runoff, and periodic environmental disturbances, including hurricanes, salinity shifts, and harmful algal blooms (O’Connell et al. 2004). Therefore, understanding these population trends is particularly important in this area, where both recreational demand and environmental variability can influence populations from year to year.

Given the high fishing pressure in Lake Pontchartrain and the ecological and economic importance of its resident species, evaluating long-term population trends is essential for effective management. Therefore, this study aims to (i) assess the population trends for several commonly found species in Lake Pontchartrain and (ii) examine changes in community structure over time to better understand the status, resilience, and dynamics of the lake’s fish assemblages.

Methods

Fish survey

From 1992 to 2019, monthly gillnet surveys were conducted across the lake across 11 sites. These sites are shown in a map below. In most years, sampling occurred in all twelve months, with one to eleven sites sampled per month. However, some years had lower sampling efforts: only five months were sampled in 2019; 2017, 2015, 2014, 2013, 2005, 2001, 1997, and 1993 were sampled for only eight to eleven months.

fish = read_csv("Data/group6.csv")

fish <- fish %>%
  mutate(date = as.Date(date))

fish <- fish %>%
  mutate(site = paste(lat, lon, sep = "_")) %>%
  mutate(site = as.factor(site),
    site = paste0(as.numeric(site)))

# Total abundance and % of total per species 
species_summary <- fish %>%
  group_by(species) %>%
  summarise(total_n = sum(n, na.rm = TRUE),
    .groups = "drop") %>%
  mutate(percentage = 100 * total_n / sum(total_n))

# All unique coordinates (lat, lon) and add site column to fish df
coords_df <- fish %>%
  distinct(lat, lon)

# Sites df
sites <- fish %>%
  distinct(lat, lon, site)

# Map with coordinates
leaflet(sites) |>
  addTiles() |>  
  addCircleMarkers(lng = ~lon,
    lat = ~lat,
    label = ~site,
    radius = 5,
    stroke = FALSE,
    color = "black",
    fillOpacity = 0.8)
# Sites / surveys  per month
sites_per_month <- fish %>%
  mutate(year = year(date),
    month = month(date)) %>%
  group_by(year, month) %>%
  summarise(n_sites = n_distinct(site),
    .groups = "drop")

sites_per_month
# A tibble: 312 × 3
    year month n_sites
   <dbl> <dbl>   <int>
 1  1992     1       1
 2  1992     2       5
 3  1992     3       3
 4  1992     4       9
 5  1992     5       4
 6  1992     6       8
 7  1992     7       7
 8  1992     8       9
 9  1992     9       3
10  1992    10       5
# ℹ 302 more rows
# Months per year
months_per_year <- sites_per_month %>%
  group_by(year) %>%
  summarise(n_months = n_distinct(month),
    .groups = "drop")

months_per_year
# A tibble: 28 × 2
    year n_months
   <dbl>    <int>
 1  1992       12
 2  1993       11
 3  1994       12
 4  1995       12
 5  1996       12
 6  1997       10
 7  1998       12
 8  1999       12
 9  2000       12
10  2001       11
# ℹ 18 more rows

Population Models

To construct the full sampling matrix, we added zeros for all sampling events in which a given species was not detected, ensuring that absence was explicitly represented in the data. For population model development, we first calculated the mean monthly abundance for each species and then summed these monthly values to obtain annual abundance estimates. Our initial approach, in which we averaged abundances abundances across all samples within each year, produced unrealistically low population estimates (often fewer than five individuals per year). Moreover, because sampling effort varied among months and years, annual averaging would disproportionately weight months with more sampling occasions and could bias abundance estimates in species exhibiting strong seasonal movement patterns. For these reasons, we adopted the monthly average-sum approach described above, which gave us more reasonable annual abundance estimates. Nevertheless, we acknowledge that these values may be much lower than true abundances and should be interpreted as relative rather than absolute population sizes. For instance, a technical report by the Louisiana Fish and Wildlife Department estimated Red Drum abundances in the millions from 1982 - 2021 (West et al. 2022). Our estimates of red drum abundances peaked at about 80. However, a study also using gill nets (O’Connell et al. 2014) had estimated abundances of some of the same fish species (e.g. Ariopsis felis) at the same magnitudes of abundances as our data.

Lastly, although the sampling locations represent distinct sites within the Pontchartrain basin, we did not analyze fish abundance at the site level for several reasons. First, sites were not sampled at equal intervals across the time series, with some sites being sampled more frequently than others in certain years. This inconsistency prevents sites from functioning as true replicates. Second, based on their locations, the sites probably differ substantially in habitat type, salinity regime, and proximity to passes or marshes, meaning that spatial variation would confound rather than clarify temporal abundance trends. Third, using sites to calculate annual abundances would be biased, as years with more sampled sites would appear to have more abundance than years with fewer sampled sites. For these reasons, we focused on month–year sampling units, which provided a more consistent basis for estimating annual population trends.

To estimate changes in population abundances over time, we created a density independent and a density dependent model for three species: Spotted seatrout (Cyanoscion nebulosus), hardhead catfish (Ariopsis felis) and Gafftopsail catfish (Bagre marinus). These three species were selected based on a combination of ecological and fisheries importance, their suitability as good candidates for gillnet surveys (i.e., their abundance estimates obtained from gillnets are less likely to be too biased compared to small, benthic or structure associated species), and the fact that they exhibited fewer outliers in the abundance dataset compared to the other 10 most abundant species.

Spotted seatrout represented 10% of all captured individuals throughout the whole study period across all species. This species is a batch spawner with high fecundity, has a long spawning season and early maturity (Brown-Peterson & Warren 2001). These characteristics make them highly responsive to environmental conditions, and due to their moderate productivity, density-dependent models are more likely to describe their population dynamics. This species is also one of the most heavily targeted species by recreational fisheries in Louisiana, making its population trends highly relevant for management (West et al. 2021).

The two catfishes are not considered high value for fisheries (Bohannon et al. 2015), however, they represented about 22% of all captured individuals, and comprise a large portion of shrimp-trawl bycatch in the Gulf (Eustis 2011). Both catfish species are known for having low fecundity rates, small brood size, late maturity ages and high parental investment, meaning they are what is considered an equilibrium strategist or K-selected species (Flinn et al. 2020, Pensinger et al. 2021). These characteristics make it more likely to show density dependent population regulation but also create a population less resilient to stressors.

Two population models were carried out to assess which best fit each species’ dynamics: density dependent and density independent. Selection of the best fitting model was based on the relationship between lambda and abundances throughout time, the projected model values compared to the real values from the dataset, and life-history traits for each species. We also performed a simulation (100 replicates) of (density dependent) population growth beyond the year 2000 for each of the selected species. We built a model of the distribution of lambdas from the beginning of the dataset (1992) until 1999. In each run of the simulation, a random lambda from this distribution was selected for each time step and the population was modeled until 2018. This was done 100 times, and we then determined the mean and the standard error to calculate the 95% confidence interval for the model. Evaluation of the performance of the models was done visually, based on its overlap with real abundance values. These results are discussed later in this report.

Community Analysis

To assess changes in the fish community throughout the study period, we first evaluated whether sampling effort was sufficient to capture most of the species diversity by generating a species accumulation curve. We then examined temporal changes in species diversity by identifying the species that comprised 75% of total abundance, and used these to estimate annual trends of alpha and beta diversity, as well as the inverse Simpson diversity index. Beta diversity was calculated two ways: using additive beta diversity to study annual changes in species composition, and using Jaccard dissimilarity with NMDS to visually examine how community composition differed among years. Finally, to evaluate the annual stability of the fish community across sites, we used a community trajectory analysis.

Results & Discussion

Population models

1.1 Use a discrete density-dependent and -independent population model to assess the population growth dynamics of any three species within the 10 most abundant species. What can you infer from both models (Model i | Species i)? For each species, what is the most proper model to assess their growth? Why?

First, we identified the top 10 most abundant species:

fish_top_species = fish %>%
  group_by(species) %>%
  summarise(total = sum(n, na.rm = TRUE), 
            .groups = 'drop') %>%  
  slice_max(total, n = 10) %>%
  ungroup()

fish_top_species
# A tibble: 10 × 2
   species             total
   <chr>               <dbl>
 1 Red Drum             2885
 2 Black Drum           2839
 3 Sea Catfish          2588
 4 Striped Mullet       2469
 5 Gulf Menhaden        1830
 6 Spotted Seatrout     1809
 7 Gafftopsail Catfish  1236
 8 Blue Crab             544
 9 Gizzard Shad          356
10 Sheepshead            143

Here we explored the abundance data for these 10 species over the years. This helped us select the species based on having fewer extreme outliers. We averaged per month and then summed those together to get an estimate of yearly abundance.

Note: A significant outlier was removed from the Spotted seatrout data (filter Nt < 400) to avoid overestimations.

abundance = fish %>%
  unite("lat_long", lat, lon, sep = "_", remove = TRUE)|>
  group_by(lat_long) |>
  mutate(site = as.factor(cur_group_id()))|> #create new column for site
  ungroup()|>
  mutate(year=year(as.Date(date)), 
         month = month(as.Date(date)))|>  # add column for year and month
  select(-sci_name)|>
  pivot_wider(names_from = species, 
              values_from = n, values_fill = 0)|> #fill in missing values
    pivot_longer("Redear Sunfish":"White Bass", names_to = "species",
               values_to = "Nt")

abundance_monthavg <- abundance |>
  filter(species %in% fish_top_species$species)|>
  group_by(species, year, month)|>
  summarise(Nt = mean(Nt), .groups = 'drop')|>
  group_by(species, year) |>
  summarise(Nt = sum(Nt), .groups = 'drop') |>
  filter(Nt <400) |>
  ungroup()
Figure 1. Observed relative annual abundance (Nt) for the top 10 most abundant species in the dataset.
ggplot(abundance_monthavg, aes(x = year, y = Nt, group =species))+
  geom_line()+
  facet_wrap(~species, scales = "free_y")

Our 3 species species selection: Spotted Seatrout, Hardhead Catfish (listed as Sea Catfish in the dataset), and Gafftopsail Catfish.

Lambda calculation at a yearly time-step, grouping all sites together.

abundance_monthavg <- abundance_monthavg|>
  filter(species %in% c('Spotted Seatrout', 'Sea Catfish', 'Gafftopsail Catfish'))

observed_density = abundance_monthavg |>
  arrange(species, year) |>
  group_by(species) |>
  mutate(lambda=(lead(Nt) / Nt)) |>
  ungroup()
Figure 2. Annual population size (Nt) plotted against lambda to visualize the effect of population size on growth.
#plotting Nt against lambda to see the per-capita growth rates
ggplot(aes(x=Nt, y=lambda, color=species, group = species), 
              data=observed_density)+
  geom_point()+
  geom_line()+
  geom_smooth(method = "lm")+
  xlab("Relative annual population size (Nt)")+
  ylab("Lambda")+
  theme_bw()+
  facet_wrap (~species)

The relationship between the relative annual population size (Nt) and the population growth rate (lambda) was negative for all three species. In general, lambda declined as population size increased (Figure 2). This negative relationship is consistent with density-dependent population regulation. This pattern suggests that density dependent models are more appropiate for describing the dynamics of these species. The sea catfish (hardheaded catfish) exhibited the strongest relationship between lambda and Nt, suggesting its population exhibited the most density-dependence. For Spotted seatrout, there was more variability in lambda, which suggests that while density dependence may still be present, other processes such as recruitment fluctuations, fishing pressure or environmental variability may play an important role in this species’ abundance.

Figure 3. Distribution of annual population growth rates (Lambda) among the three species.
#boxplot of lambda values
ggplot(observed_density, aes(x=species, y=lambda, color = species))+
  geom_boxplot()+
  stat_summary(fun = "mean", geom = "point", shape = 8,
               size = 2, color = "black")+
  scale_y_continuous(limits = c(0, 8), breaks = c(0, 2, 4, 6, 8))+
  xlab("Species")+
  theme_bw()+
  theme(legend.position = "none")

Mean lambda values were similar across species, generally being above 1, which indicates that on average, none of the species exhibited strong long-term decline, and exhibited growth on average (Figure 3). The lambda values are skewed due to a large number of extreme increases in the population sizes. However, the variability of lambda differed among species. For both catfish species, variability was relatively low suggesting a stable population (Figure 3). Spotted seatrout showed the highest variability, including three high-lambda outliers. However, spotted seatrout exhibits strong seasonal movements in and out of Lake Pontchartrain or between areas in response to salinity and temperature (Baer 2019). Threfore, their annual abundance estimates derived from gillnet data could potentially reflect seasonal availability rather than true population size, which may increase interannual variability in lambda.

Calculation of density-dependent model parameters: r, r(d), and K.

# Calculating r and rd
densityd = observed_density %>%
  arrange(species, year) %>%
  group_by(species) %>%
  mutate(r=log(lambda),
         rd=lambda-1) %>%
  na.omit()%>%
  ungroup()

densityd$r[is.infinite(densityd$r)] <- NA

#Calculating K-estimate from linear regression of r and Nt
dd_results <- densityd %>%
group_by(species) %>%
  summarise(fit = list(lm(r ~ Nt, data = cur_data())),
            .groups = "drop") %>%
  mutate(a = map_dbl(fit, ~ coef(.x)[1]),
         b = map_dbl(fit, ~ coef(.x)[2]),
         K = -a / b)

dd_long <- dd_results %>%
  select(species, K)
Figure 4. Annual population size (Nt) for Gafftopsail catfish, hardhead catfish and spotted seatrout from 1992 to 2019, with estimated carrying capacity (K) shown as a horizontal line for each species.
#Plot of real data over time with K estimates for each species
ggplot(abundance_monthavg, aes(x = year, y = Nt, color=species)) +
  geom_line() +
  geom_hline(data = dd_long, aes(yintercept= as.numeric(K), color = species))+
  facet_wrap(~species, scales = "free_y")+
  scale_color_discrete(name = "Species")+
  xlab("Year")+ ylab("Observed relative annual population size (Nt)")+
  theme(legend.position = "none")

Annual population size fluctuated for all species, with several pronounced peaks through the time series. All species showed years with sharp increases in estimated abundance, which were then followed by rapid declines. Those patterns may be influenced by gillnet encouters with larger schools or recruitment pulses, however, all species showed a strong decline in their abundance after or around 2010. This decline could be attributed changes in environmental conditions. During 2008 and 2011, the Bonnet Carre spillway was opened to mitigate flooding in New Orleans. This spillway diverted river water into Lake Pontchartrain, carrying large nutrient loads (Mize et al. 2009, Roy et al. 2013). High nutrient loads resulted in multiple algal blooms which are known to affect fish survival and recruitment, potentially contributing to the observed declines.

When viewed in relation to carrying capacities (horizontal lines, Figure 4), the populations of Gafftopsail catfish remain mostly below K, suggesting their population may be more regulated by other external drivers rather than by density dependence. In contrast, Hardhead (Sea) catfish and Spotted seatrout had several years approaching or exceeding the estimated K, which suggests these populations are not living at equilibrium and is subject to the effects of episodic environmental events, recruitment pulses, temporary migrations or a bias in sampling method rather than a population regulated solely by density dependence.

Let’s make a projection of the density-independent and density-dependent models of the three species to further evaluate which model best fits the observed population trends. For both models, the same initial conditions were used: the geometric mean of lambda or r(d), and an initial population of 2. The initial population was chosen as it matched the initial estimated abundance values at the beginning of the dataset. Each model was run for 150 years.

Density-independent model calculation:

#find geometric mean
lambdas <- observed_density |>
  group_by(species) |>
  filter(lambda>0)|>
  summarise(geom.mean = exp(mean(log(lambda), na.rm = TRUE)),
            sd = sd(lambda, na.rm = TRUE))

#using these values, lets project with both types of models. we will use the same starting values for both. 

#lets use 10 as our starting value as our populations are relatively small
years = 150
N0 = 2
species = 3
popt <- matrix(NA, nrow = years, ncol = species)

for (h in 1:nrow(lambdas)) {
    N <- numeric(years)
    N[1] <- N0
  
     for (t in 2:years) {
        N[t] <- N[t-1] * lambdas$geom.mean[h]
  }
  popt[, h] <- N
  }

diproj <- as.data.frame(popt) %>%
  setNames(c("Gafftopsail Catfish", "Sea Catfish", "Spotted Seatrout")) %>%
    mutate(year = 1992:(years+1991)) %>%
  pivot_longer(cols = "Gafftopsail Catfish": "Spotted Seatrout", 
               names_to = "species", values_to = "Nt")
Figure 5. Density-independent population projection for 150 years, starting in year 1992.
ggplot(diproj, aes(x= year, y =Nt, color = species)) +
  geom_line(aes(x = year, y = Nt, color = species))+
  theme_bw()+
  scale_color_discrete(name = "Species")+
  labs(title = "Density Independent Model Projection")+
  theme(plot.title = element_text(hjust = 0.5))+
  xlab("Year")

Density dependent growth calculation:

#now lets do the density dependent projection, same lambda and initial values
dlogisticD = function(K, rd, N0, years) {
  logis = tibble(Nt = NA, time = 1:years)
  logis$Nt[logis$time == 1] = N0
  
  for(i in 1:(length(logis$time)-1)){
    logis$Nt[i+1] = logis$Nt[i]  + rd*logis$Nt[i]*(1 - logis$Nt[i] / K)
  }
  return(logis)
}

#bind K and calculate rd values from geometric mean of lambda
rd <- lambdas$geom.mean -1
dd_long$rd <- rd

#create vectors
Ks <- dd_long$K
rds <- dd_long$rd

#apply function
ddproj <- dd_long |>
  mutate(pop = map2(Ks, rds, \(Ks, rds) dlogisticD(K = as.numeric(Ks), 
                                                   rd = as.numeric(rds), 
                                                   N0 = 2, years = 150)))|>
  unnest(pop)|>
  mutate(year = time + 1991)
Figure 6. Density-dependent population projection for 150 years, starting in year 1992.
ggplot(ddproj, aes(x = year, y = Nt, color = species))+
  geom_line()+
   xlab("Year")+
   theme_bw()+
  scale_color_discrete(name = "Species")+
  labs(title = "Density Dependent Model Projection")+
  theme(plot.title = element_text(hjust = 0.5))

We then evaluated how the models behaved compared with the estimated population data, which took into account only the first 37 years of the two models.

Figure 7. Density dependent and independent model projections from year 1992 to year 2019 represented in lines. Gillnet-derived population size data shown as points.
ddproj2 <- ddproj |>
  select(species, Nt, year)|>
  filter( year < 2020)

observed_density2<- observed_density |>
  select(species, Nt, year)

diproj2 <- diproj |>
  filter( year < 2020)

ggplot(data = observed_density2, aes(x = year, y = Nt, group = species))+
  geom_point()+
  geom_line()+
  geom_line(data = ddproj2, aes(x = year, y = Nt, color = 'red'))+
  geom_line(data = diproj2, aes(x = year, y = Nt, color = 'blue'))+
  ylab("Observed relative annual abundance Nt")+
  xlab("Year")+
  facet_wrap(~species, scales = 'free')+
  scale_color_discrete(name = "Model",
                       labels=c('Density Dependent', 'Density Independent'))

For all three species, density-dependent and density-independent models produced similar projected population trajectories for the 37 years that coincided with the dataset (Figure 7). Both models showed slight positive trends over the study period, and were generally low and stable (flat) relative to the observed annual abundances. For Hardhead catfish (Sea catfish), observed abundances fluctuated from year to year, but were generally high. However. both models again predicted only small increases in population size overtime. Overall, both models do not capture the magnitude or fluctuations presernt in the observed data, and appeared remarkably similar to eachother. For the purpose of this study, based on Figure 2 and life history traits of these species, we selected the density-dependent model for further analysis.

1.2 Based on the selected models, how does lambda change over time?

Lambda projection calculation for each model:

ddproj <- ddproj |>  
  arrange(species, year) |>
  group_by(species) |>
  mutate(lambda=(lead(Nt) / Nt))|>
  na.omit() |>
  ungroup()

diproj <- diproj |>
  arrange(species, year) |>
  group_by(species) |>
  mutate(lambda=(lead(Nt) / Nt))|>
  na.omit() |>
  ungroup()
Figure 8. Annual variation in population growth rate (Lambda) for Gafftopsail catfish, Hardhead (sea) catfish and Spotted seatrout from 1992 to 2019. Black lines show the observed lambda values derived from year to year changes in estimated abundance: N(t+1) / N(t). Horizontal colored lines represent the projected lambda for both density-independent and density-dependent models for each species.
ggplot(observed_density, aes(x = year, y = lambda, group = species))+
  geom_line()+
  geom_line(data = ddproj, aes(x = year, y = lambda, color = "blue"),
            linewidth = .5)+
  geom_line(data = diproj, aes(x = year, y = lambda, color = "red"),
            linewidth =.5)+
  theme_bw()+
  scale_x_continuous(limits = c(1990, 2020), 
                     breaks = c(1990, 1995, 2000, 2005, 2010, 2015, 2020))+
  xlab("Year")+
  ylab("Lambda")+
  facet_wrap(~species)+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))+
  scale_color_discrete(name = "Model",
                       labels=c('Density Independent', 'Density Dependent'))

Population growth rates varied largely from year to year for all species, with several sharp spikes and crashes. These large changes indicates high interannual varability in population dynamics that could be driven by environmental fluctuations, variable recruitment or sampling noise due to gillnets capturing large schools. Despite the high variability in observed lambda, both model-estimated lambda values (horizontal lines) were constant remained close to 1. We can infer that our modeled populations (both density-dependent and density independent) were largely inaccurate due to a non-varying lambda.

1.3 Starting with the average abundance in 2000, perform a simulation (N = 100 replicates) of population growth using the growth variability based on the abundance values from start to 1999 (i.e., use 100 samples of lambda from 1980 to 1999s frequency distribution to project population beyond 2000). How many times are the real abundance values intercepted by the uncertainty bounds of the simulation (i.e., values within the 95% confidence intervals)? Briefly discussed the performance of the models and potential reasons for the results.

Model simulation calculations:

#using density dependent model
grvar <- observed_density |>
  group_by(species) |>
  filter(year < 2000)|>
  mutate(rd = lambda - 1)|>
  summarise(mean = mean(rd, na.rm = TRUE),
            sd = sd(rd, na.rm = TRUE))|>
  column_to_rownames(var = "species") 
  

#100 replicates
years = 20
reps = 100
N0 <- observed_density |>
      filter(year == 2000) |>
      pull(Nt)

set.seed(224)
#in my for loop i made any negative values equal zero because negative pop values are not realistic
popt <- matrix(NA, nrow = years, ncol = reps)
list.pop <- vector("list",3)
for (h in 1:nrow(grvar)) {
    
  for (r in 1:reps) {
    N <- numeric(years)
    N[1] <- N0[h]
    rd <- rtruncnorm(28, a = -1, mean = grvar[h,1], sd = grvar[h,2])
    for (t in 2:nrow(popt)) {
     N[t] = N[t-1]  + rd[t-1]*N[t-1]*(1 - (N[t-1] / Ks[h]))
     N[t] <- ifelse(N[t] <0 | is.nan(N[t]) | is.infinite(N[t]), 
                    N[t] == 0, N[t] <- N[t])
    }
  popt[, r] <- N
  }
  meanpop <- apply(popt, 1, mean, na.rm = TRUE)
  sds <- apply(popt, 1, sd, na.rm = TRUE)
  list.pop[[h]] <- data.frame(meanpop, sds)
} 

grvar <- grvar |>
  rownames_to_column("species")

names(list.pop) <- grvar$species

#plot
dd.plot <- as.data.frame(do.call(rbind, list.pop))|>
  rownames_to_column(var = "spy")|>
  separate_wider_delim(spy, ".", names = c("species", "time"))
  

dd.plot$time <- as.numeric(dd.plot$time)
dd.plot$year <- dd.plot$time + 1999

#calculate confidence interval
dd.plot <- dd.plot |>
  mutate(error = (qnorm(0.975)*sds)/(sqrt(years)))
Figure 9. Density dependent population model simulations for Gafftopsail catfish, Hardhead (Sea) catfish and Spotted seatrout from 2000 to 2019. Solid black lines represent the predicted annual population trajectories derived from the fitted density-dependent model. Shaded regions show 95% confidence intervals of model predictions.
ggplot(dd.plot, aes(x = year, y = meanpop))+
  geom_line(linewidth = .5)+
  geom_point(size =.8) +
  geom_ribbon(aes(ymin=meanpop -error, ymax=meanpop+error),alpha=0.3)+
  xlab("Year")+
  ylab("Population Estimate")+
  labs(title = "Population Models: 2000 - 2019")+
   theme(plot.title = element_text(hjust = 0.5))+
  facet_wrap(~species, scales ='free_y')

Figure 10. Comparison of density-dependent model simulations (Black line and uncertainty shaded area) with observed gillnet-derived abundance data (Colored points) for Gafftopsail catfish, Hardhead (Sea) catfish and Spotted seatrout from 2000 to 2019.
observed_density.model <- observed_density |>
  filter(year > 1999)


ggplot(dd.plot, aes(x = year, y = meanpop)) +
  geom_point()+
  geom_ribbon(aes(ymin=meanpop -error, ymax=meanpop+error),alpha=0.3)+
  geom_point(data = observed_density.model, aes(x=year, y = Nt, color = species))+
  geom_line(data = observed_density.model, aes(x=year, y = Nt, color = species))+
  facet_wrap(~species) +
  ylab("Population Size")+
  xlab("Year")+
  theme_bw()+
  labs(title = "Population Models v. Sample Data: 2000 - 2019")+
   theme(plot.title = element_text(hjust = 0.5))+
  theme(legend.position = "none")

Density dependent model simulations (Figure 9) showed large uncertainty, and observed population trajectories from 2000 to 2019 fell outside the range predicted by the models most of the time. Overall, even with varying lambdas at each time step, the modeled population dynamics did not capture all of the variablility shown in the estimated abundance data. Gafftopsail catfish showed a gradual increase through time in the simulations, while the observed data showed high fluctuations. Only three years fell within the model’s predicted range (Figure 10). For Hardhead (Sea) catfish, the simulations predicted a relatively stable population, with the observed data showing high fluctuation. Observed data from only two years fell within the model’s boundaries. Spotted seatrout simulations indicated long-term decline, with observed values falling withing the model’s boundaries five times. The spotted seatrout model may have shown greater agreement with the data (based on intersections) only because of the larger confidence intervals calculated by the model.

Even though the models did not capture the full magnitude of year to year variation present in the gillnet-derived data (Figure 9 and Figure 10), the density-dependent model seemed to capture the general direction and trend of the population sizes. This pattern is expected, as simple density-dependent models do not incorporate environmental stochasticity, sampling noise or species specific movements, therefore, they are not intended to replicate every annual observation but capture the general population dynamics. The general trends captured by the model are supported by other published papers. Spotted seatrout populations are known to be overfished and declining (West et al. 2021). Catfish populations in the Gulf are stable or growing due to diminished bycatch as a result of reduced shrimping effort (Walters et al. 2006, Eustis 2011).

Community Analysis

2.1 Does the sampling assigned to your group reach the asymptote of the species accumulation curve? Discuss results.

# Create site x date matrix
fish <- fish %>%
  mutate(sample_id = paste(site, date, sep = "_"))

comm_mat <- fish %>%
  group_by(sample_id, species) %>%
  summarise(N = sum(n), .groups = "drop") %>%
  pivot_wider(names_from = species,
    values_from = N,
    values_fill = 0 ) %>%
  column_to_rownames("sample_id")

# Run species acc curve
sac <- specaccum(comm_mat, method = "random", permutations = 100)
sac_df <- tibble(samples   = sac$sites,
  richness  = sac$richness,
  sd = sac$sd) %>%
  mutate(lower = richness - sd,
    upper = richness + sd)
Figure 11. Species accumulation curve for all gillnet sampling events conducted from 1992 to 2019 in the Lake Pontchartrain basin. The solid line represents the cumulative number of fish species detected as number of samples increased, the shaded area represents 95% confidence interval.
ggplot(sac_df, aes(x = samples, y = richness)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), 
              fill = "lightblue", alpha = 0.5) +
  geom_line(size = 0.5, color = "black") +
  labs(x = "Number of Samples",
    y = "Cumulative Species Richness",
    title = "Species Accumulation Curve") +
  theme_classic(base_size = 14)

The species accumulation curve shows a rapid rise in species richness during the first portion of the surveys, followed by a level off as additional samples contributed fewer species (Figure 11). This pattern indicates that the gillnet survey adequately captured the majority of the species detectable with this gear in Lake Pontchartrain. Although the final plateau doesn’t fully flatten, the slow rate of species added by the end indicates that the fish community is well represented by the data. In addition, the accumulation curve may not reach the asymptote due to species turnover throughout the time series.

2.3 - Using a community trajectory analysis, assess the annual stability of the fish community across the sites from the assigned basin and gear type. Describe and discuss the main results.

# Prep data
fish <- fish %>%
  mutate(date = as.Date(date),
    year = year(date),
    site = as.factor(site))

fish$site <-factor(fish$site, levels=c("1", "2", "3", "4", "5", "6", "7", "8",
                             "9", "10", "11"))
# Build community data: one row per site-year, one column per species
comm_df <- fish %>%
  group_by(site, year, species) %>%
  summarise(n = sum(n), .groups = "drop") %>%
  arrange(site, year, species) %>%
  pivot_wider(names_from  = species,
              values_from = n,
              values_fill = 0)

# Extract site and survey (year) vectors
sites   <- comm_df$site 
surveys <- comm_df$year

# Community matrix (only species columns)
comm_matrix <- comm_df %>%
  select(-site, -year) %>%
  as.data.frame()

# Bray-Curtis & trajectories 
# Dissimilarity matrix
b_dist <- vegdist(comm_matrix, method = "bray")

# Define trajectories: each site moves through years
d <- defineTrajectories(b_dist,
                        sites   = sites,
                        surveys = surveys)
Figure 16. Annual community change expressed as total trajectory length (path length) for each sampling site in Lake Pontchartrain. Bars represent the cumulative amount of year to year change at each site, calculated as the sum of dissimilarity distances through time.
## Trajectory lengths (stability) 
tl <- trajectoryLengths(d) %>%
  rownames_to_column(var = "site") %>%
  as_tibble()

tl$site <-factor(tl$site, levels=c("1", "2", "3", "4", "5", "6", "7", "8",
                             "9", "10", "11"))


# Total path length per site (overall temporal change; shorter = more stable)
ggplot(tl, aes(site, Path, fill = site)) +
  geom_col() +
  labs(x = "Site",
       y = "Total trajectory length (Path)",
       title = "Annual community change")+ #lower = more stable
  scale_fill_viridis_d() +
  theme_bw() +
  theme(legend.position = "none")

Figure 17. Distribution of annual trajectory length (path lengths) for each sampling site in Lake Pontchartrain basin. Boxplots show the mean and interquartile range of path lengths, with individual points representing yearly values.
# Segment lengths (month-to-month / year-to-year change) 
# Automatically detect segment columns S1, S2, ..., Sk
seg_cols <- grep("^S[0-9]+", names(tl), value = TRUE)

tll <- tl %>%
  pivot_longer(cols = all_of(seg_cols),
               names_to = "seg",
               values_to = "len") %>%
  mutate(seg_num = as.numeric(sub("S", "", seg))) %>%
  arrange(site, seg_num)

# Summary per site (mean and SD of segment lengths)
tll %>%
  group_by(site) %>%
  summarize(mean_seg = mean(len),
    sd_seg   = sd(len),
    .groups  = "drop")
# A tibble: 11 × 3
   site  mean_seg sd_seg
   <fct>    <dbl>  <dbl>
 1 1        0.717  0.240
 2 2        0.770  0.192
 3 3        0.705  0.209
 4 4       NA     NA    
 5 5        0.729  0.186
 6 6       NA     NA    
 7 7       NA     NA    
 8 8        0.574  0.237
 9 9        0.764  0.201
10 10      NA     NA    
11 11       0.623  0.181
# Boxplot of segment lengths (how variable change is within site)
ggplot(tll, aes(site, len, fill = site)) +
  geom_point(aes(color = site),
             size = 1,
             position = position_jitterdodge()) +
  geom_boxplot(outliers = FALSE, alpha = 0.6) +
  labs(x = "Site",
    y = "Segment length") +
  scale_fill_viridis_d() +
  theme_bw() +
  theme(legend.position = "none")

Figure 18. Speed of community change through time across sampling sites. Each colored line represents a different site.
## Speed of change vs segment index 
ggplot(tll, aes(seg_num, len, group = site, color = site)) +
  geom_point(size = 2) +
  geom_line(linewidth = 1) +
  labs(x = "Time (years)",
    y = "Segment length (speed of change)") +
  scale_color_viridis_d() +
  facet_wrap(~site)+
  theme_bw() +
  theme(legend.position = "none")

Figure 19. Distance from initial community state over time for all sites. Temporal trajectories show how far each site’s community composition deviated from the first year sampled across the monitoring period.
# Distance from baseline (first year) 
tl_init <- trajectoryLengths(d, relativeToInitial = TRUE) %>%
  rownames_to_column(var = "site") %>%
  as_tibble()

# Columns like Lt1_t2, Lt1_t3, ...
base_cols <- grep("^Lt1_t", names(tl_init), value = TRUE)

tl_init_long <- tl_init %>%
  pivot_longer(cols = all_of(base_cols),
    names_to  = "step",
    values_to = "len") %>%
  mutate(step_num = as.numeric(sub("Lt1_t", "", step))) %>%
  arrange(site, step_num)

ggplot(tl_init_long, aes(factor(step_num), len, group = site, color = site)) +
  geom_point(size = 2.5) +
  geom_line(linewidth = 1) +
  labs(x = "Time (years)",
    y = "Distance from baseline") +
  scale_color_viridis_d() +
  facet_wrap(~site)+
  theme_bw() +
  theme(legend.position = "none")

Figure 20. Accumulated community change over time.
# Cumulative distance (accumulated change)
tll <- tll %>%
  group_by(site) %>%
  mutate(dist_cuml = cumsum(len)) %>%
  ungroup()

ggplot(tll, aes(seg_num, dist_cuml, group = site, color = site)) +
  geom_point(size = 2.5) +
  geom_line(linewidth = 1) +
  labs(x = "Time (years)",
    y = "Cumulative distance") +
  scale_color_viridis_d() +
  theme_bw() +
  theme(legend.position = "none")

Community-trajectory analysis revealed substantial variation in fish community stability across sites, suggesting an episodic rather than a gradual change. Community trajectories for each site (Figure 15) showed that all locations experienced frequent directional shifts in species composition. Some sites (site 2, 3 and 7) showed long distance shifts in particular years, indicating major community shifts during that time. Other sites, such as 9 had a wide and elongated trajectory, indicating frequent, recurrent shifts, making it appear the most vulnerable site to disturbances. In contrast, points in sites 4, 8, 10, and 11 are much closer together and have shorter paths, reflecting a smaller annual turnover and a relatively more stable community.

Despite these site-level differences, overall annual community change, quantified as total trajectory length, was broadly similar among sites (Figure 16). Site 2 and 9 showed the highest levels of overall temporal change, but were not too far apart from the other sites. The distribution of segment lengths (Figure 17) further supports this conclusion: each site showed a wide range of year to year change, some sites having smaller shifts than others, consistent with environmental disturbances, recruitment pulses or rapid changes in dominant species. Even the sites with the smaller boxes still contained extreme large or small values, indicating that all sites experience at least some years with rapid atypical shifts. However, in general, all sites have overlapping values, suggesting no diffferences in their compositional changes.

Temporal dynamics were also evident when examining the speed of community change. Segment lengths plotted over time (Figure 18) showed that high-magnitude shifts occurred intermittently throughout the study period, rather than clustering in a particular decade. Distance from the initial community state (Figure 19) showed no consistent pattern, instead, communities oscillated around moderate to high distances from the baseline composition.

Finally, the accumulated community change (Figure 20) showed that total compositional displacement increased steadily through time across all sites. Although slopes varied, all sites had an upward trajectory, suggesting that all changes accumulated into a long-term differentiation from earlier community states.

Together, these results demonstrate that the Lake Pontchartrain fish community is characterized by high variability and large shifts, and has been experiencing a long-term cumulative change, rather than reaching stability. This is likely due to the effects of multiple disturbances such as environmental variability, human pressures or species turnover.

Conclusion

Across all analyses, Lake Pontchartrain’s fish assemblage showed evidence of substantial temporal variability and long-term change. Population models revealed strong interannual fluctuations in abundance for all three focal species, with pronounced peaks and declines and a marked drop in population size around 2010, likely linked to large-scale environmental events such as nutrient-driven algal blooms associated with the Bonnet Carre spillway openings. Although simple density-dependent models did not capture the full magnitude of year-to-year variability, they reproduced the general population trends, consistent with expectations for species whose dynamics are influenced by environmental variability and sampling noise. Diversity metrics similarly indicated long-term shifts: alpha diversity declined after 2010, beta diversity showed a progressive homogenization of communities, and inverse Simpson diversity revealed decreased species richness and reduced evenness. Community-trajectory analyses further demonstrated that fish assemblages experienced recurrent directional shifts at all sites, with episodic spikes in the speed of change and steadily increasing cumulative displacement from initial community states. Together, these results show that the Lake Pontchartrain fish community has undergone sustained, long-term alteration driven by a combination of environmental disturbances and species-specific dynamics, and has not returned to a stable or previous state over the nearly three-decade monitoring period.

References

Baer, A. M. 2019.?Population dynamics and demographics of spotted seatrout, Cynoscion nebulosus, through spatial analysis: Towards an integrative management approach. [Master’s thesis, Louisiana State University and Agricultural & Mechanical College]

Bohannon, C., Esslinger, J., & Wagner, T. 2015. Trends in Texas commercial fishery landings, 1994–2012 (Management Data Series No. 290) [Data set]. Texas Parks and Wildlife Department, Coastal Fisheries Division.

Brown-Peterson, N. J., & Warren, J. W. 2001) The reproductive biology of spotted seatrout, Cynoscion nebulosus, along the Mississippi Gulf Coast. Gulf of Mexico Science, 19(1), 7–19. https://doi.org/10.18785/goms.1901.07

Cates L, Hale S, Olsen Z. 2023 Gafftopsail catfish in Texas estuaries: population trends and ecosystem implications. N Am J Fish Manag. https://doi.org/10.1002/nafm.10967

Coleman, F. C., Figueira, W. F., Ueland, J. S., & Crowder, L. B. 2004. The impact of United States recreational fisheries on marine fish populations. Science, 305(5692), 1958-1960. https://doi.org/10.1126/science.1100397


Eustis, S. P. 2011. Bycatch of the Lake Pontchartrain Basin inshore shrimp fishery an its effect on two sea catfish species: The Gafftopsail Catfish (Bagre marinus) and the Hardhead Catfish (Ariopsis felis) [Master’s thesis, University of New Orleans].

Flinn, S., Midway, S., & Ostrowski, A. 2020. Age and Growth of Hardhead Catfish and Gafftopsail Catfish in Coastal Louisiana, USA. Marine and Coastal Fisheries: Dynamics, Management, and Ecosystem Science, 11(5), 362-371. https://doi.org/10.1002/mcf2.10089

Keddy, P. A., D. Campbell, T. McFalls, G. P. Shaffer, R. Moreau, C. Dranguet, and R. Heleniak. 2007. The wetlands of Lakes Pontchartrain and Maurepas: Past, present and future. Environmental Reviews 15:43–77.?

Louisiana Department of Wildlife & Fisheries. 2022. Lower Pontchartrain Sub-Basin Waterbody Management Plan. Office of Fisheries, Inland Fisheries Section, Baton Rouge, Louisiana.

McCorquodale, J. A., Roblin, R. J., Georgiou, I. Y., & Haralampides, K. A. 2009. Salinity, nutrient, and sediment dynamics in the Pontchartrain Estuary. Journal of Coastal Research, (10054), 71-87.

Mize, S. V., & Demcheck, D. K. 2009. Water quality and phytoplankton communities in Lake Pontchartrain during and after the Bonnet Carr? Spillway opening, April to October 2008, in Louisiana, USA. Geo-Marine Letters, 29, 431–440. https://doi.org/10.1007/s00367-009-0157-3

National Oceanic and Atmospheric Administration, Fisheries. 2023. Status of Stocks 2023: The year in review. U.S. Department of Commerce.https://www.fisheries.noaa.gov/national/sustainable-fisheries/status-stocks-2023

O’Connell, M. T., R. C. Cashner, and C. S. Schieble. 2004. Fish assemblage stability over fifty years in the Lake Pontchartrain estuary; comparisons among habitats using canonical correspondence analysis. Estuaries 27:807–817.?

Martin T. O’Connell, Ann MU O’Connell, and Christopher S. Schieble. (post print) “Response of Lake Pontchartrain Fish Assemblages to Hurricanes Katrina and Rita” Estuaries and Coasts 37.2 (2014): 461-475.

Roy, Eric Daniel. 2013. A multi-scale investigation of nutrient dynamics in the Lake Pontchartrain Estuary and Basin. [Doctoral Thesis, Louisiana State University and Agricultural & Mechanical College]

Skrobialowski, S. C., Green, W. R., & Galloway, J. M. 2007. Water quality of Lake Pontchartrain and outlets to the Gulf of Mexico following Hurricanes Katrina and Rita: Chapter 7E in Science and the storms-the USGS response to the hurricanes of 2005 (No. 1306-7E, pp. 221-224). US Geological Survey.

Sumaila, U. R., and T. C. Tai. 2020. End overfishing and increase the resilience of the ocean to climate change. Frontiers in Marine Science 7.?

Vasconcelos, R. P., P. Reis-Santos, M. J. Costa, and H. N. Cabral. 2011. Connectivity between estuaries and marine environment: Integrating metrics to assess estuarine nursery function. Ecological Indicators 11:1123–1133.

Walters, C., Martell, S.J.D., Mahmoudi, B., 2006. An Ecosim model for exploring ecosystem management options for the Gulf of Mexico: implications of including multistanza life history models for policy predictions. In: Presentation for Mote Symposium #6, Florida State University, November 2006, p. 23.

West J., X. Zhang, T. Allgood, and J. Adriance. 2021. Update assessment of Spotted Seatrout (Cynoscion nebulosus) in Louisiana Waters 2021 Report. Louisiana Department of Wildlife and Fisheries.

West J., E. Lang, X. Zhang, J. Adriance. 2022. Assessment of Red Drum Sciaenops ocellatus in Lousiana 2022 Report. Louisiana Department of Wildlife and Fisheries.