Abstract

Environmental factors can often have population effects on aquatic organisms, though studies of environmental effects are often specific to a given life stage. Stage-structured demographic models provide a means of exploring the multivariate suite of life history parameters associated with a population and can provide a tool for understanding population-wide impacts of single stage events (e.g., mortality and fecundity). Here, the unique dynamics observed for an isolated population of black drum (Pogonias cromis) and the population-wide impacts of salinity as a driver of young-of-the-year (YOY) survival were investigated. This modeling exercise revealed that the dynamics observed in the black drum population are potentially driven by increased survival in the post-settler/YOY stage and that this increased survival is at least partially a result of the higher salinities that typify the Upper Laguna Madre of Texas, a hypersaline estuary (34% increase in population growth rate across the range of salinity examined). Early maturation in this population was also shown to have the potential to provide infrequent, large pulses of fecundity to the population. Quantifying the population-wide impact of such drivers can place management decisions into the context of the environment and provide both a proof-of-concept for specific management action and realistic expectations for managers and constituents alike. Without such formal quantification, it will be difficult for habitat concepts to move from an abstract management tool to widespread application.

Environmental effects on fishes have long been studied and cited as potential drivers of population dynamics (Baird 1873; Beck et al. 2001; Peterson 2003). Researchers often attempt to identify trends or optima of growth, survival, or abundance for various environmental parameters. Scharf (2000) and Olsen (2019) sought to recognize such correlates of growth and mortality among juvenile red drum (Sciaenops ocellatus) and black drum (Pogonias cromis), respectively. Similarly, Peterson et al. (1999) examined salinity impacts on Atlantic croaker (Micropogonias undulatus) growth rates and Doerr et al. (2016) conducted an in-situ experiment to determine optimum salinity preferences for juvenile white (Litopenaeus setiferus) and brown shrimp (Farfantepenaeus aztecus). Although such studies are often specific to a given life stage (e.g., juveniles), they are potentially impactful to the full population dynamics. Exploring and formally quantifying population effects of such stage-specific impacts is necessary for the application of this information to fisheries and their respective ecosystems management.

Matrix-based stage-structured demography (Leslie 1945; Caswell 2001) provides an analytical structure for quantifying population impacts of stage-specific vital rates. In general, demographic analysis extends the focus of population modeling beyond abundance of the various life stages to the relevant vital rates of each life stage that potentially impact overall population abundance (i.e., mortality and fecundity). For example, Levin & Stunz (2005) employed this technique to examine the potential impacts of nursery habitat restoration on the red drum (Sciaenops ocellatus) population in Galveston Bay. Similarly, Wong & Dowd (2016) examined the production potential for various fish species resulting from coastal habitat restoration. Although environmental relationships are well understood for many populations, the population-level impacts of these relationships are rarely quantified. Such a formal quantification can place management decisions into the context of the environment and provide both a proof-of-concept for specific management action and realistic expectations of for managers and constituents alike.

To demonstrate the utility and importance of such population-wide analyses, I investigated the causation of unique dynamics observed in an isolated population of black drum (Anderson et al. in press). In Texas waters, black drum support a commercial fishery worth nearly 1.5 million dollars per year (ex-vessel value; Bohannon et al. 2015), with greater than half of these landings typically occurring in the Upper Laguna Madre (ULM). The ULM is a hypersaline estuary that supports an adult population of black drum with an abundance three to six times greater than that observed in other Texas bay systems (Olsen 2014). Additionally, recruitment of juveniles in this system is up to two orders of magnitude greater than that observed elsewhere (Olsen 2014). Increasing trends in both adult abundance and recruitment began in the early 1990s, which coincided with several other unique events including record drought conditions, winter freeze events, and the onset of a long-term algal bloom known as “brown-tide” that has since impacted the system (Buskey et al. 2001; Wetz et al. 2017). Calculated juvenile mortality is known to be lower in the ULM as compared to other Texas bays and density, mortality, and growth of juvenile black drum in this system has been linked to salinity (Olsen 2019). Further, black drum in the ULM are known to mature earlier than in other regions (Bumguardner et al. 1996; Olsen et al. 2018) presumably resulting in increased population fecundity. While these unique aspects of mortality and fecundity in the ULM population have been studied, the synthesis of this information into a single model is necessary to pinpoint stage-specific drivers and to assess the magnitude of stage specific and environmental control on the observed population dynamics.

Thus, an examination of fecundity and mortality, their simulated impacts on the overall population, and known environmental drivers of these vital rates will better allow us to explore and understand observed trends and further elucidate the potential links of environmental parameters and system wide perturbations. The objectives of this study were to (1) construct a stochastic stage-structured demographic model for black drum in the ULM using published and calculated vital rates, (2) assess the relative importance of each stage specific vital rate in driving the growth rate of this population, and (3) exemplify a use of this model for examining the effects of salinity on stage level vital rates.

Materials & Methods

Model construction.–The ULM is located on the south Texas coast and represents a unique hypersaline estuary in which salinities regularly exceed that of seawater (i.e., >35 ppt), and salinity gradients are often negative (i.e., highest salinities in the upper estuary – here termed the Baffin Bay complex; Tunnell & Judd 2002). These characteristics are a result of two factors: greater evaporation than freshwater inflow that typify this estuary and limited direct connection to the Gulf of Mexico (Fig. 1). A five-stage demographic model was constructed to exemplify the life history of the ULM black drum population (Fig. 2). The associated size and age for each stage are given in Figure 3. Adults were grouped into two stages to examine the potential impact of early maturation on population dynamics. The ULM population is unique in that maturation begins to occur at the adult 1 stage of the model while other black drum populations do not mature until the adult 2 stage (Olsen et al. 2018).

Figure 1.

Map of the Upper Laguna Madre showing the secondary bay system, Baffin Bay. Packery Channel, which is the only direct connection to the Gulf of Mexico, is indicated with a star.

Figure 1.

Map of the Upper Laguna Madre showing the secondary bay system, Baffin Bay. Packery Channel, which is the only direct connection to the Gulf of Mexico, is indicated with a star.

Figure 2.

Schematic of the five-stage demographic model used in this analysis. Calculations for fecundity (Fi), probability of survival to the next stage (Gi), and the probability of remaining in the current stage (Pi) are given in the Methods section.

Figure 2.

Schematic of the five-stage demographic model used in this analysis. Calculations for fecundity (Fi), probability of survival to the next stage (Gi), and the probability of remaining in the current stage (Pi) are given in the Methods section.

Figure 3.

Timeline of fish age as it relates to stages in the demographic model shown in Figure 2.

Figure 3.

Timeline of fish age as it relates to stages in the demographic model shown in Figure 2.

The core of a demographic model is a projection matrix that consists of calculated vital rate information divided into sequential life stages in the rows and columns of the matrix (Stevens 2009). All analyses are based on this projection matrix. Such analysis allows for an understanding of the sequential, multivariate interactions of the various stages in the model and for an understanding of stage specific drivers of larger scale patterns in the matrix. Following Leslie (1945) and Caswell (2001), there are three parameters for each stage in the demographic model: survival within the stage (Pi), survival and movement to the next stage (Gi), and fecundity (Fi), where i indicates the stage of interest (Fig. 2). Calculations of Pi and Gi followed that of Crouse et al. (1987):
formula
formula
where pi is the probability of survival during stage i and di is the stage duration. Probability of survival (pi) was calculated from instantaneous mortality (Z) as follows:
formula
Where e is Napier’s coefficient, Zi is instantaneous mortality at stage i, and Di is the stage duration given in the unit of time for which Zi was calculated. The demographic model was constructed such that each model iteration represents a single year. Thus, Pi is only relevant to stages with a duration >1 year (i.e., adult 1 and adult 2; Fig. 3). For the larval stage, Zi was based on published values from Cowan et al. (1992) under various levels of predation. This was an in-situ study on larval black drum mortality based on a 30-d time interval (Di) and represents the only published estimates of black drum larval mortality. These values are slightly higher, though comparable to those used in Levin & Stunz (2005) for red drum larvae in Galveston Bay, Texas and calculated for marine larvae by Comyns et al. (2003) for the northern Gulf of Mexico. For the remaining classes (post-settler/young-of-the-year (YOY), sub-adult, adult 1 and adult 2), Zi values were calculated from a fishery independent dataset collected and maintained by the Texas Parks and Wildlife Department (Martinez-Andrade 2015). This dataset encompasses the entirety of the Texas coast from 1983–present. Data used in this study are from the Upper Laguna Madre from 1983–2015 and from three commonly used sampling gears: bag seines (18.3 m long by 1.8 m deep; 19 mm and 13mm nylon mesh in the wings and bag, respectively), trawls (6.1 m wide; 31 mm nylon mesh), and gill nets (182.9 m long by 1.2 m tall; four panels of differing monofilament mesh size: 152 mm, 127 mm, 102 mm, and 76 mm). Based on the length frequency of capture (Fig. 4), data from each of these gears were assigned to the stage they best represented: bag seine data-post-settler/YOY, trawl data- sub-adult, and gillnet data-adult 1/adult 2. For each stage, Zi was calculated with a catch curve analysis (Ricker 1975). The catch curve analysis exploits sequential decreases in catch rates for the various black drum stages either seasonally (when stage duration was 1 year or less, i.e., post-settler/YOY and sub-adult stages) or annually (when stage duration was >1 year, i.e., adult 1 and adult stages) as a means of quantifying mortality. This is based on an exponential model of decline in catch:
formula
where Ni is relative abundance at time t (based on fishery independent catch rates), N0 is initial relative abundance, e is as defined above, and Zt is instantaneous mortality across the timeframe of the analysis. For the post-settler/YOY stage calculation, monthly bag seine catch rates for black drum 50–200 mm were used. The sub-adult stage calculation used monthly trawl catch rates for black drum 200–299 mm. Black drum in the adult 1 stage (300–449 mm) are not fully selected in any of the gears (i.e., they are too big to be efficiently captured in the bag seine and trawl gear and too small for the gill net gear). For this reason, Z values were assumed to be comparable for adult 1 and adult 2. The adult 2 stage calculation used annual gill net catch rates for black drum ≥450 mm. Because this stage spans multiple year classes, catch rates were calculated by year class based on published von Bertallanfy growth curves (Bumguardner et al. 1996; Olsen et al. 2018) and decreasing catch rates across years within year classes used in the catch curve analysis. Fish in year classes 5–15 were used for calculation of adult 2 mortality.
Figure 4.

Length frequency of black drum captured in fishery independent sampling gears.

Figure 4.

Length frequency of black drum captured in fishery independent sampling gears.

Fecundity (Fi) was calculated as follows:
formula
where WBFi is weighted batch fecundity at stage i, 100 represents the length (in days) of the black drum spawning season for a given year in the northern Gulf of Mexico, and 4 is the spawning frequency for black drum during the spawning season in the northern Gulf of Mexico (Nieland & Wilson 1993). Given that past work has identified variation in population size and size at maturity within the Baffin Bay complex of the ULM (Olsen et al. 2018), WBFi was calculated as a weighted average accounting for the relative proportion of the mature black drum population in the Baffin Bay complex and the remainder of the ULM:
formula
where BFi is batch fecundity at stage i based on Bumguardner et al. (1996), PMBBi is the proportion mature at stage i in the Baffin Bay complex, CPUEBB is the proportion of fishery independent catch rates of the ULM observed in the Baffin Bay complex, PMULMi is the proportion mature at stage i in the remainder of the ULM, and CPUEULM is the proportion of fishery independent catch rates in the ULM observed in the remainder of the ULM (i.e., outside of Baffin Bay). Fecundity was assumed to be zero for larval, post-settler/YOY, and sub-adult stages (Olsen et al. 2018). For adult 1 and adult 2 stages, BFi (and thus Fi) was taken from Bumguardner et al. 1996 and PPBBi and PPULMi values were taken from Olsen et al. (2018). CPUEBB and CPUEULM were calculated from TPWD fishery independent gill net catch rates. These values are 0.64 and 0.36, respectively (i.e., the average proportion of the total amount of black drum captured in the Upper Laguna Madre as a whole is 0.64 from Baffin Bay and 0.36 from the remainder of the Upper Laguna Madre).

Two separate methods were used to create the demographic model and to meet the objectives of this analysis. First, a base model was constructed using a Monte-Carlo approach with simulations of Zi (with conversions to Pi and Gi as given in equations 13 above) and batch fecundity (with conversions to Fi as given in equation 56 above) used to construct 1000 stochastic iterations of the projection matrix assuming a normal distribution (though truncated at 0) with calculated mean and standard deviations given in Table 1. The purpose of this stochastic base model was to examine general patterns in the population given natural variation about the vital rates (Zi and Fi). The second method explicitly examined the influence of salinity on black drum YOY survival and the resulting impacts on the dynamics of the full population. Olsen (2019) fit a bell curve model to the relationship between annual mean salinity and YOY black drum mortality. Thus, to model the role of the salinity – YOY Z relationship on overall population dynamics, I used this relationship to determine YOY Z across the range of observed annual salinities (Fig. 5) in the ULM and fed these outputs into the stochastic base model. This resulted in stochastic demographic models (1000 iterations of the projection matrix for each model) with varying YOY Z relationships based on salinities across the following values: 25 ppt, 30 ppt, 35 ppt, 40 ppt, 45 ppt, 50 ppt, and 55 ppt. This range of salinities was based on annual mean salinities typically observed in the ULM (Tunnell & Judd 2002).

Table 1.

Mean (standard deviation) instantaneous mortality (Z) and annual fecundity (F) values used to draw Monte-Carlo simulations for each stage in the demographic model. Size range of stage members and stage duration are given. A schematic of the demographic model is given in Figure 2 and a timeline for the various stages in this model is given in Figure 3.

Mean (standard deviation) instantaneous mortality (Z) and annual fecundity (F) values used to draw Monte-Carlo simulations for each stage in the demographic model. Size range of stage members and stage duration are given. A schematic of the demographic model is given in Figure 2 and a timeline for the various stages in this model is given in Figure 3.
Mean (standard deviation) instantaneous mortality (Z) and annual fecundity (F) values used to draw Monte-Carlo simulations for each stage in the demographic model. Size range of stage members and stage duration are given. A schematic of the demographic model is given in Figure 2 and a timeline for the various stages in this model is given in Figure 3.
Figure 5.

The modeled salinity-young-of-the-year mortality (Z) relationship calculated by Olsen (2019) and here used to model the impact of this relationship on overall population dynamics for black drum in the Upper Laguna Madre. Demographic model outputs are given for 25, 30, 35, 40, 45, 50, and 55 ppt.

Figure 5.

The modeled salinity-young-of-the-year mortality (Z) relationship calculated by Olsen (2019) and here used to model the impact of this relationship on overall population dynamics for black drum in the Upper Laguna Madre. Demographic model outputs are given for 25, 30, 35, 40, 45, 50, and 55 ppt.

Model analysis.–Eigenanalysis was used to analyze the iterated projection matrices of the demographic models such that:
formula
where A is the projection matrix, λ is an eigenvalue, and w is a right eigenvector (termed ‘right because it is solved with w on the right side of A). Essentially, eigenanalysis finds an independent set of solutions that captures all information in A (similar to a principal components analysis and, in fact, eigenanalysis is at the core of principal components analysis). While there are an infinite number of solutions, the first solution typically captures the most important features of A. In this context, the first eigenvalue (λ1) represents the intrinsic rate of population increase where values < 1 correspond to decreased population growth rates and values > 1 correspond to increased population growth rates. Additionally, the first right eigenvector (w1) gives the stable-stage distribution and is presented as the proportion of each stage under a stable intrinsic rate of population increase (i.e., λ = 1). Similarly, the eigenvector can be solved from the left such that:
formula
where A and λ are as presented above and v is the left eigenvector (termed ‘left’ because it is solved with v on the left side of A). Here, the first left eigenvalue (v1) represents the reproductive value for each stage which is presented as the proportional contribution of each stage to future stages based on fecundity and survival to that stage. Further, elasticities were calculated on λ as follows:
formula
where Eij is the elasticity of element aij in the projection matrix, λ is as presented above, and δλ and δaij represent small changes in these elements of the matrix. Elasticities represent the proportional impact of each input parameter (Pi, Gi, and Fi) on the intrinsic rate of population increase. All parameter calculation, model construction, and model analyses were conducted in R (R Core Team 2016) using the “truncnorm” package to create truncated normal distributions for the Monte-Carlo approach (Trautmann et al. 2014) and R base packages for all other model analysis (R Core Team 2016). R code for demographic analysis largely followed that given in Stevens (2009).

Results

Demographic base model.–Instantaneous mortality coefficients (Zi) calculated for this study were largely comparable to those observed for other sciaenids (e.g., Levin & Stunz 2005) and for other populations of black drum (e.g., Jones & Wells 1998; Table 1). Additionally, the general pattern of calculated Zi observed among stages is expected for most fishes (i.e., decreasing mean and standard deviation of Zi with increasing age, Haddon 2001). While the increased mean Zi from the post-settler/YOY stage (mean ± standard deviation; 0.09 ± 0.08) to the sub-adult stage (0.10 ± 0.06) does break from this pattern, standard deviations for these stages still follow this pattern and show considerable overlap. Fecundity calculations were largely based on published values and were observed to increase with age (i.e., from adult 1 to adult 2) as would be expected for most fishes (Hixon et al. 2014).

The projection matrix resulting from mean mortality and fecundity (Table 1) is given in Table 2 with the distribution for each matrix element (based on the of the variation about these mean values and 1000 iterations of this base model) given in Figure 6. Lambda values (i.e., intrinsic rate of population increase) calculated on these 1000 iterations of the base model were right skewed and ranged from 0.60 to 45.98 with a median value of 1.43 (Table 3). Trends in fishery independent catch rates for black drum in the ULM suggest extreme booms in recruitment of juveniles (bag seine surveys, Fig. 7a) and subadults (trawl surveys, Fig. 7b) with a generally positive trend in the older stages (gill net surveys, Fig. 7c) across the timeseries. These observed trends seem logical in the contexts of the observed distribution of simulated lambda values. As is expected for general life history strategies typical of many fishes (Winemiller 2005), the stable stage distribution was overwhelmingly dominated by the larval stage and each successive stage systematically decreased in abundance and the distribution of these iterations moved from highly left skewed (larval stage) to highly right skewed (YOY/post settler, sub-adult, adult 1) to relatively even (adult 2; Table 4). Median reproductive values were largest for the YOY/post-settler stage and increased for each successive stage (Table 3).

Table 2.

The projection matrix calculated from mean values of mortality and fecundity given in Table 1. Values given along the top row of the matrix represent fecundity while values along the diagonal of the matrix represent survival as calculated in the Methods section.

The projection matrix calculated from mean values of mortality and fecundity given in Table 1. Values given along the top row of the matrix represent fecundity while values along the diagonal of the matrix represent survival as calculated in the Methods section.
The projection matrix calculated from mean values of mortality and fecundity given in Table 1. Values given along the top row of the matrix represent fecundity while values along the diagonal of the matrix represent survival as calculated in the Methods section.
Figure 6.

Plotted distributions of projection matrix parameters with mean values (Table 3) in bold. These distributions are based on 1000 iterations of a normal distribution (truncated at 0) with mean and standard deviations of mortality and fecundity values given in Table 2. Calculations for the parameters of G, P, and F are given in the Methods sections. Stages included larvae, post-settler/young-of-the-year (YOY), sub-adult, adult 1, and adult 2. Stage definitions are given in Figure 3.

Figure 6.

Plotted distributions of projection matrix parameters with mean values (Table 3) in bold. These distributions are based on 1000 iterations of a normal distribution (truncated at 0) with mean and standard deviations of mortality and fecundity values given in Table 2. Calculations for the parameters of G, P, and F are given in the Methods sections. Stages included larvae, post-settler/young-of-the-year (YOY), sub-adult, adult 1, and adult 2. Stage definitions are given in Figure 3.

Table 3.

Calculated metrics from the eigenanalysis of the demographic base model of black drum in the Upper Laguna Madre, Texas. Summary statistics are given from the 1000 iterations of the projection matrices used for metric calculations. Elasticities are further given for each parameter in the projection matrix: G, P, and F.

Calculated metrics from the eigenanalysis of the demographic base model of black drum in the Upper Laguna Madre, Texas. Summary statistics are given from the 1000 iterations of the projection matrices used for metric calculations. Elasticities are further given for each parameter in the projection matrix: G, P, and F.
Calculated metrics from the eigenanalysis of the demographic base model of black drum in the Upper Laguna Madre, Texas. Summary statistics are given from the 1000 iterations of the projection matrices used for metric calculations. Elasticities are further given for each parameter in the projection matrix: G, P, and F.
Figure 7.

Black drum catch rates from fishery independent data collected in the Upper Laguna Madre. Bag seine catch rates (catch/hectare) represent the post-larval/young-of-the-year stage, trawl catch rates (catch/hour) represent the sub-adult stage, and the gill net catch rates (catch/hour) represent the adult 1 and adult 2 catch rates.

Figure 7.

Black drum catch rates from fishery independent data collected in the Upper Laguna Madre. Bag seine catch rates (catch/hectare) represent the post-larval/young-of-the-year stage, trawl catch rates (catch/hour) represent the sub-adult stage, and the gill net catch rates (catch/hour) represent the adult 1 and adult 2 catch rates.

Table 4.

A summary of lambda values calculated on 1000 iterations of a demographic model with salinity varied young-of-the-year (YOY)/post-settler mortality. The following summary statistics are given: mean, standard deviation (SD), minimum value (min.), quartile 1 (Q1), median, quartile 3 (Q3), and maximum value (max.).

A summary of lambda values calculated on 1000 iterations of a demographic model with salinity varied young-of-the-year (YOY)/post-settler mortality. The following summary statistics are given: mean, standard deviation (SD), minimum value (min.), quartile 1 (Q1), median, quartile 3 (Q3), and maximum value (max.).
A summary of lambda values calculated on 1000 iterations of a demographic model with salinity varied young-of-the-year (YOY)/post-settler mortality. The following summary statistics are given: mean, standard deviation (SD), minimum value (min.), quartile 1 (Q1), median, quartile 3 (Q3), and maximum value (max.).

The largest median elasticities calculated on iterations of the projection matrices were for larval, post-settler/YOY, and sub-adult G values. The P parameters calculated on the adult 1 and adult 2 stages were the widest ranging of calculated elasticities (with largest maximum elasticity values). Of the F parameters, the adult 2 stages had the largest median elasticity (almost four times larger than that of adult 1). However, the elasticity of the adult 1 F stage had a larger maximum value and a more skewed distribution than the elasticity of the adult 2 stage (Table 3; Fig. 8).

Figure 8.

Plotted elasticities for each parameter calculated from iterations of the demographic base model. Calculations for the parameters G, P, and F are given in the Methods section. Stages included larvae, post-settler/young-of-the-year (YOY), sub-adult, adult 1, and adult 2. Stage definitions are given in Figure 3.

Figure 8.

Plotted elasticities for each parameter calculated from iterations of the demographic base model. Calculations for the parameters G, P, and F are given in the Methods section. Stages included larvae, post-settler/young-of-the-year (YOY), sub-adult, adult 1, and adult 2. Stage definitions are given in Figure 3.

Demographic salinity model.–While the distribution of lambda values from iterations of the projection matrix spanned from <1 (i.e., decreasing population) to >1 (i.e., increasing population) for the range of salinity values examined, all these distributions were heavily skewed to the right (increasing population) with all median values >1 (Table 4). Lambda increased with increasing salinity across the range of values examined. Median values of lambda increased 34% from a low of 1.54 at 35 ppt to a high of 2.07 at 55 ppt (Fig. 9).

Figure 9.

Plotted lambda values from iterations of the demographic model resulting from changes to young of the year Z values with increasing salinity. See Figure 5 for this salinity-Z relationship. Dashed lines represent quartile 1 and quartile 3 while the solid line represents median values.

Figure 9.

Plotted lambda values from iterations of the demographic model resulting from changes to young of the year Z values with increasing salinity. See Figure 5 for this salinity-Z relationship. Dashed lines represent quartile 1 and quartile 3 while the solid line represents median values.

Discussion

Demographic base model.–The lambda values calculated from the demographic base model iterations indicate that given the input parameters (i.e., mortality and fecundity) and typical variation about these parameters, the ULM black drum population is generally increasing and only rarely decreasing in abundance at the annual timestep. In fact, the upper quartile and maximum values of lambda calculated from the iterations of the projection matrix suggest some possibility for extreme population growth rates. Such extreme swings have been observed in this population. Olsen (2014) noted both a rapid increase in adult abundance beginning in the early 1990s in addition to booms in YOY recruitment that were orders of magnitude larger than black drum recruitment observed in adjacent bay systems. These trends observed by Olsen (2014) continue into the present (Fig. 7). Further, it has been suggested that this population exists at or near carrying capacity with seasonal breaches in carrying capacity resulting in observed emaciation of individuals in this system (Olsen 2016; Ajemian et al. 2018). The asymptotic trend in the Figure 7c implies this as well. These observations suggest that model outputs regarding lambda and the population that would presumably result are at least somewhat representative of the actual population.

With regards to the reproductive values, the output of this metric for the post-settler/YOY and sub-adult stages can potentially be misleading considering that these stages exert no direct reproductive input to the population (i.e., individuals in this stage of the model are not considered reproductively mature and so are assigned a fecundity of 0). However, this reproductive metric is calculated on the full population matrix and is affected not only by stage fecundity but by the relative influence of a stage’s survival to downstream, reproductively mature stages. Thus, in the demographic base model, it is an indirect influence of these stages on the reproductive capacity of the population that results in the high proportion contribution to future stages communicated by the calculated reproductive values. Further, the distribution of reproductive values among stages suggests that while the abundance of the younger stages largely impacts the typical reproductive capacity of the population, the spikes in fecundity associated with the oldest fish in the population have the potential to impact population growth (i.e., the distribution of reproductive values in these later stages are skewed towards larger values and contain larger upper quartile and maximum values). The concept of the oldest and largest individuals in a population providing a large proportion of the population’s fecundity is common in many fish species (Nieland & Wilson 1993; Hixon et al. 2014). However, the ULM population of black drum is unique from other black drum populations in that individuals have been found to mature at age 2 or 3 rather than age 5 for most populations (Olsen et al. 2018). In this model, these ‘early maturing’ fish are represented in the adult 1 stage. Due to smaller body size, the adult 1 stage exhibits lower fecundity compared to the adult 2 stage. However, the stable stage distribution for the population suggests a greater relative abundance at these ‘early maturing’ stages compared to the adult 2 stage and we see this potential impact in the distribution of elasticities for the adult 1 F parameter.

In the context of the population scale, our metric for understanding the relative importance of each stage and their associated vital rates is the calculation of elasticities. Elasticities measure the proportional contribution of each model parameter (Gi, Pi, and Fi) to the intrinsic rate of population increase (λ). That is, they assess the population wide impact of stage-specific vital rates. It should be noted that the summary of values given in Table 3 do not sum to one because these are summaries from 1000 iterations of the input parameters and their associated distributions. As discussed above, these iterations make the model more robust to natural variation in these vital rates and thus more representative of reality. Elasticities should still be interpreted as they would be for a single projection matrix. The larger median elasticities of survival (G) at the larval, post-settler/YOY, and sub-adult stages suggest that it is the survival at these early stages that has the greatest impact on the overall population growth rate and resulting dynamics of this population. Interestingly, black drum do not enter the commercial or recreational fishery until they reach the adult 1 stage (minimum size limit of 355 mm). Similarly, Levin & Stunz (2005) found that the bulk of the red drum population growth rate could be explained by survival of the larval and juvenile stages. Levin & Stunz (2005) used this information to support restoration and protection of nursery habitats for this species.

Although median elasticities for G are largest for these first three stages, they are not overwhelmingly large (0.157) compared to that of other matrix parameters. The distributions of elasticities for adult 1 P (stationary survival within the stage) and adult 2 P along with their large maximum values (1) suggest that these parameters can have a dominating impact on the overall population growth rate under extreme scenarios, although this is probably not the norm. Regarding elasticities of the fecundity parameters (F), it is interesting to note the larger maximum values at the adult 1 stage compared to the adult 2 stage despite the larger median value at the adult 2 stage. This skew in the adult 1 F elasticities towards larger values suggests a potential for this ‘early maturation’ stage to play a large role in the dynamics of this population (though infrequently). However, the more even distribution observed in the adult 2 stage highlights the stability of reproductive input for the largest females in the population. Olsen et al. (2018) highlight the dominant role that this early maturation plays in defining the unique life history strategy of black drum in the ULM population compared to other populations throughout their range. This analysis revealed that the survival of larvae and YOY plays a major role in the unique population dynamics and absolute abundance of black drum in the ULM. Relatively uniform distribution of elasticities among the various model parameters suggest that no single characteristic of this population is overwhelmingly causative of the observed population trends. It is likely that both depressed YOY mortality and early maturation (i.e., fecundity at sub-adult and adult 1 stages) contribute to the abundance of black drum in the ULM.

Demographic salinity model.–Vital rates are known to be intimately related to environmental conditions (Beck et al. 2001). The uniquely hypersaline conditions in the ULM have been linked to YOY survival of black drum (Olsen 2019). I applied the structure of the demographic base model to changing YOY Z values across a range of salinities based on an analysis by Olsen (2019). The purpose of this application was to formally quantify the potential for measurable population-wide impacts of environmental parameters. Understanding the extent of this impact is important when managing fisheries and the ecosystems in which they exist.

Across the range of salinity–YOY Z values examined, there was a general pattern of increasing lambda at the highest salinity values though it should be noted that mortality would certainly begin to increase at higher salinities than those examined here. In their study of red drum, Levin & Stunz (2005) stated that an increase of two percent in lambda due to habitat-related changes was significant. This analysis revealed that increasing salinity resulted in a 34% increase in median lambda, which suggests that salinity mediated juvenile mortality considerably influences black drum population dynamics in the ULM. However, the ULM is notable for both hypersalinity and extreme interannual variations in salinity (Tunnell & Judd 2002). Such extreme salinity variation probably results in highly variable YOY mortality and likely dampens such long-term population growth to some extent. Recruitment of black drum in the ULM is known to be highly variable (Fig. 7a, b), although adult abundances have, nonetheless, remained stable over the past decades (Olsen 2014; Fig. 7c). This salinity relationship was previously known (Olsen 2014; Olsen 2019) but the extent of its population-wide impact had only been speculated.

In summary, it seems that the unique dynamics observed in the ULM black drum population are potentially driven by the increased survival in the YOY stage and that this increased survival is at least partially a result of the higher salinities that typify this hypersaline estuary. It should be noted, however, that the causative nature of this salinity–Z relationship is not explicitly examined in this model. Olsen (2014) suggested the possibility for an indirect salinity relationship such that high salinity periods resulted in decreased predation on black drum juveniles. Further, Street et al. (1997) suggested the potential for trophic disruptions in this system resulting from “brown-tide” algal blooms known to correlate with extreme hypersalinity. No such relationships are explicitly accounted for in this model and represent areas that are in need of further research.

Early maturation in this population may also play a role in driving observed population dynamics in some years, though with limited frequency as suggested by the highly skewed elasticities for the adult 1 F values. More consistent elasticities (i.e., less skewed iterations of the model) for the adult 2 F values suggest more consistent (though less potential for extreme) reproductive input to the population. Thus, while the adult 1 stage has the potential to provide infrequent but large pulses of fecundity into the population, the more stable adult 2 stage provides a regular, high level of fecundity that may sustain the population more consistently. It is likely some combination of these factors that drives the large and highly stable (in recent decades) black drum population in the ULM. To note, this modeling exercise did not explicitly differentiate fishing mortality versus natural mortality in calculations of Z (though Z estimates were presumed to include fishing mortality). Neither does it assume immigration and emigration of the ULM population to neighboring systems. While movement of black drum to neighboring population does likely occur to some extent, several studies suggest that black drum in the ULM are at least somewhat isolated to this region (Ajemian et al. 2018; Anderson et al. in press). However, relative to neighboring bay systems, fishing mortality (especially commercial fishing mortality) is much higher in the ULM (Bohannon et al. 2015). Given current minimum and maximum length restrictions on black drum fishing (i.e., fish must be between 355 mm and 762 mm for harvest), fishing mortality is primarily impacting adult 1 and early adult 2 stage fish. While this study did not seek to specifically evaluate the impact of the current level of fishing mortality on the population, it does suggest that the current maximum size limit is a valid strategy to preserve population fecundity. Trends in the population do indicate that current levels of fishing in the ULM are generally sustainable (Fig. 7). However, given the variability in ULM population demographics and the potential impact of salinity variation on these demographics, these trends and environmental relationships should continue to be closely monitored.

Fisheries and ecosystem managers are often tasked with both the identification of vital habitat parameters (e.g., salinity regimes, seagrass coverage, etc.) and discerning the relative influences of environmental (i.e., natural) effects. The stage-structured demographic model provides a framework for exploration of observed population trends and the potential drivers of those dynamics. In the present example, understanding the quantitative influence of salinity on the ULM black drum population will help fisheries managers better understand observed juvenile recruitment trends and resulting population dynamics. In this context, decisions regarding the health of the black drum population in the ULM will be more robust. Thus, the demographic modeling approach used here can both highlight conservation strategies that surround the protection of certain habitat parameters and modify management expectations in the context of existing habitat quality. Without such formal, population-wide quantification, it will be difficult for habitat concepts to move from an abstract management tool to widespread application.

Acknowledgments

I would like to thank past and present technicians, biologists, and support staff of the Texas Parks and Wildlife Department Upper Laguna Madre management team for the collection and maintenance of the fisheries independent dataset on which many of these vital rates calculations and supporting observations were based. I would also like to thank Mark Fisher and a number of anonymous peer reviewers for their helpful critiques and comments on early versions of this manuscript. Finally, I would like to thank my graduate committee at Texas A&M University – Corpus Christi: Gregory Stunz, James Tolan, Jennifer Pollack, and Paul Montagna.

LITERATURE CITED

LITERATURE CITED
Ajemian,
M. J.,
Mendenhall
K. M.,
Pollack
J. B.,
Wetz
M. S.,
&
Stunz
G. W.
.
2018
.
Moving forward in a reverse estuary: habitat use and movement patterns of black drum (Pogonias cromis) under distinct hydrological regimes
.
Estuar. Coasts.
41
:
1410
1421
.
Anderson,
J.,
Williford
D.
&
Olsen
Z.
.
In press.
Estuarine level genomic variation confirms demographic and life history differences among black drum Pogonias cromis populations in Texas
.
N. Am. J. Fish.
Baird,
S. F.
1873
.
Report on the condition of the sea fisheries of the south coast of New England in 1871 and 1872
.
Report of the United States Fish Commission,
vol 1
.
Government Printing Office
,
Washington, D.C.
.
Beck,
M. W.,
Heck
K. L.
Jr
,
Able
K. W.,
Childers
D. L.,
Eggleston
D. B.,
Gillanders
B. M.,
Halpern
B.,
Hays
C. G.,
Hoshino
K.,
Minello
T. J.,
Orth
R. J.,
Sheridan
P. F.
&
Weinstein
M. P.
.
2001
.
The identification, conservation, and management of estuarine and marine nurseries of fish and invertebrates
.
BioScience
51
:
633
641
.
Bohannon,
C.,
Esslinger
J.
&
Wagner
T.
.
2015
.
Trends in commercial fishery landings, 1994–2012
.
Texas Parks and Wildlife Department – Coastal Fisheries Division Management Data Series No. 290, Austin, Texas.
Bumguardner,
B. W.,
Colura
R. L.,
Young
E.,
Westbrook
D.
&
Buckley
R.
.
1996
.
Black drum life history in Texas bays with emphasis on the Upper Laguna Madre
.
Final Report for project funded through U. S. Department of Interior, Fish and Wildlife Service under DJ 15.605 (Grant F-36-R, Project 18).
Buskey,
E. J.,
Lu
H.,
Collumb
C.
&
Bersano
J. G. F.
.
2001
.
The decline and recovery of a persistent Texas brown tide algal bloom in the Laguna Madre (Texas, USA)
.
Estuaries
24
:
337
346
.
Caswell,
H.
2001
.
Matrix population models: construction, analysis, and interpretation
.
Sinauer Associates, Inc.
,
Sunderland, MA
.
Comyns,
B. H.,
Shaw
R. F.
&
Lyczkowski-Shultz
J.
.
2003
.
Small-scale spatial and temporal variability in growth and mortality of fish larvae in the subtropical northcentral Gulf of Mexico: implications for assessing recruitment success
.
Fish. Bull.
101
:
10
21
.
Cowan
J. H.,
Jr.,
Birdsong
R. S.,
Houde
E. D.,
Priest,
J. S.
Sharp
W. C.
&
Mateja
G. B.
.
1992
.
Enclosure experiments on survival and growth of black drum eggs and larvae in lower Chesapeake Bay
.
Estuaries
15
:
392
402
.
Crouse,
D. T.,
Crowder
L. B.
&
Caswell
H.
.
1987
.
A stage-based population model for loggerhead sea turtle and implications for conservation
.
Ecology
68
:
1412
1423
.
Doerr,
J. C.,
Liu
H.
&
Minello
T. J.
.
2016
.
Salinity selection by juvenile brown shrimp (Farfantepenaeus aztecus) and white shrimp (Litopenaeus setiferus) in a gradient tank
.
Estuar. Coasts.
39
:
829
838
.
Haddon,
M.
2001
.
Modelling and quantitative methods in fisheries
.
CRC Press
,
Boca Raton
.
Hixon,
M. A.,
Johnson
D. W.
&
Sogard
S. M.
.
2014
.
BOFFFFs: on the importance of conserving old-growth age structure in fishery populations
.
ICES J. Mar. Sci.
71
:
2171
2185
.
Jones,
C. M.
&
Wells.
B.
1998
.
Age, growth, and mortality of black drum, Pogonias cromis, in the Chesapeake Bay region
.
Fish. Bull.
96
:
451
461
.
Leslie,
P. H.
1945
.
On the use of matrices in certain population mathematics
.
Biometrika
33
:
183
212
.
Levin,
P. S
&
Stunz.
G. W.
2005
.
Habitat triage for exploited fishes: can we identify essential “Essential Fish Habitat?” Est
.
Coast. Shelf. Sci.
64
:
70
78
.
Martinez-Andrade,
F.
2015
.
Marine resource monitoring operations manual. Texas Parks and Wildlife Department, Austin, Texas
.
Nieland,
D. L.
&
Wilson.
C. A.
1993
.
Reproductive biology and annual variation of reproductive variables of black drum in the northern Gulf of Mexico
.
Trans. Am. Fish. Soc.
122
:
318
327
.
Olsen,
Z.
2014
.
Potential impacts of extreme salinity and surface temperature events on population dynamics of black drum, Pogonias cromis, in the Upper Laguna Madre, Texas
.
Gulf Mex. Sci.
1–2
:
60
68
.
Olsen,
Z.
2016
.
Emaciated black drum (Pogonias cromis) in the Upper Laguna Madre, Texas: tracking the recovery of the population over two years
.
Tex. J. Sci.
68
:
79
90
.
Olsen,
Z.
2019
.
Quantifying nursery habitat function: variation in habitat suitability linked to mortality and growth for juvenile black drum in a hypersaline estuary
.
Mar. Coast. Fish.
11
:
86
96
.
Olsen,
Z.,
McDonald
D.
&
Bumguardner
B.
.
2018
.
Interspecific variation in life history strategies and implications for management: a case study of black drum (Pogonias cromis) in Baffin Bay
,
Texas USA. Fish. Res.
207
:
55
62
.
Peterson,
M. S.
2003
.
A conceptual view of environment-habitat-production linkages in tidal river estuaries
.
Rev. Fish. Sci.
11
:
291
313
.
Peterson,
M. S.,
Comyns
B. H.,
Rakocinski
C. F.
&
Fulling
G. L.
.
1999
.
Does salinity affect somatic growth in early juvenile Atlantic croaker, Micropogonias undulatus, (L.)?
J. Exp. Mar. Ecol.
238
:
199
207
R Core Team.
2016
.
R: A language and environment for statistical computing
.
R Foundation for Statistical Computing
,
Vienna, Austria
.
Ricker,
W. E.
1975
.
Computation and interpretation of biological statistics of fish populations
.
Fish. Res. Board Can. Bull.
191
,
Ottawa
.
Scharf,
F. S.
2000
.
Patterns in abundance, growth, and mortality of juvenile red drum across estuaries on the Texas coast with implications for recruitment and stock enhancement
.
Trans. Am. Fish. Soc.
129
:
1207
1222
.
Stevens,
M. H. H.
2009
.
A primer of ecology with R
.
Springer Science+Business Media
,
New York, NY
.
Street
G. T.,
Montagna
P. A.
&
Parker
P. L.
.
1997
.
Incorporation of brown tide into an estuarine food web
.
Mar. Ecol. Prog. Ser.
152
:
67
78
.
Trautmann,
H.,
Steuer
D.,
Mersmann
O.
&
Bornkamp
B.
.
2014
.
Truncnorm: truncated normal distribution. R package version 1.0–7
.
Tunnell
J. W.
Jr.,
&
Judd
F. W.
(eds.).
2002
.
The Laguna Madre of Texas and Tamaulipas
.
Texas A&M University Press
,
College Station, TX
,
xxii+346
pp.
Wetz,
M. S.,
Cira
E. K.,
Sterba-Boatwright
B.,
Montagna
P. A.,
Palmer
T. A.
&
Hayes
K. C.
.
2017
.
Exceptionally high organic nitrogen concentrations in a semi-arid South Texas estuary susceptible to brown tide blooms
.
Est. Coast. Shelf. Sci.
188
:
27
37
.
Winemiller,
K. O.
2005
.
Life history strategies, population regulation, and implications for fisheries management
.
Can. J. Fish. Aquat. Sci.
62
:
872
885
.
Wong,
M. C.
&
Dowd.
M.
2016
.
A model framework to determine the production potential of fish derived from coastal habitat for use in habitat restoration
.
Estuar. Coasts.
39
:
1785
1800
.