Growth during early life history is recognized as an important determinant of recruitment success in fish populations (Houde 2008). Moreover, growth rates become increasingly important if selective processes reduce mortality in larger or faster-growing individuals relative to their smaller, slower-growing counterparts (Sogard 1997). This differential mortality may be due to decreased predation risk, increased resilience to starvation and environmental change, or some combination of these factors as individuals become larger (van der Veer et al. 2000). By these mechanisms, even minute changes in growth rates during early life history can substantially impact survivability and recruitment strength (Houde 1987).

Winter Flounder Pseudopleuronectes americanus is a historically important species that once supported a robust commercial and recreational fishery in the inshore bays of Long Island, New York, until a precipitous decline of local populations occurred in the 1990s. Winter Flounder in this region are presently considered to be at record low abundances (Nuttall et al. 2011 ; Sagarese and Frisk 2011). In spite of a commercial moratorium in federal waters and recreational catch limits in New York, Winter Flounder populations have shown little to no signs of recovery in Long Island waters. In the most recent stock assessment, spawning stock biomass estimates in the southern New England and Middle Atlantic Bight, which includes Long Island populations, ranged from 16% to 24% of target levels, and recent recruitment has been consistently lower than the long-term average (NEFSC 2011). Furthermore, microsatellite analyses of Winter Flounder from Long Island bays have revealed an alarmingly high incidence of inbreeding, which is likely the result of a very low effective population size in these bays (O'Leary et al. 2013). Previous studies have proposed that inadequate recruitment of young-of-the-year (age-0) fish is limiting recovery in this region (Socrates and Colvin 2006; Yencho 2009), although the specific mechanisms underlying recruitment limitation in Long Island bays have not been extensively studied. Moreover, daily mortality in Long Island bays has been observed to be higher than in Connecticut and Rhode Island estuaries (Yencho 2009; L.A.H., unpublished data), and identifying the factors that influence growth may provide insight into the potential causes of these high rates of mortality.

Throughout their range, adult Winter Flounder enter estuaries from the continental shelf during the winter and early spring to spawn (Klein-MacPhee 2002). Juveniles are typically retained within the estuaries where they were spawned and spend the first 2 years of life in these inshore habitats. Thus, the environmental conditions of these inshore habitats are considered to be important determinants of growth and survival in juvenile Winter Flounder (Goldberg et al. 2002). For example, previous studies in other estuarine systems have demonstrated that habitat suitability for age-0 Winter Flounder is driven primarily by dynamic and sometimes simultaneous changes in temperature, dissolved oxygen, salinity, and food availability (Meng et al. 2001; Manderson et al. 2002). The timing of settlement, wherein pelagic larvae metamorphose into demersal, laterally compressed juveniles, also influences the environmental conditions that juveniles are exposed to and thus may affect long-term growth and other vital rates (Sogard et al. 2001; Manderson et al. 2003).

Growth in Winter Flounder would be expected to decline with increasing age, and studies in both field (Sogard and Able 1992; Phelan et al. 2000) and laboratory (Keller and Klein-MacPhee 2000) settings have indicated that high temperatures (typically > 21°C) negatively influence growth in Winter Flounder. However, the effect of temperature may diminish later in the growing season as fish grow larger (Sogard et al. 2001). Long Island bays are near the southern extremity of the range of Winter Flounder, which may suggest that fish in this region are more likely to experience long periods of unfavorably high temperatures (Manderson 2008; Able et al. 2014). Many watersheds in Long Island are also highly developed and have been historically impacted by nutrient and contaminant loading (Bopp et al. 1993). Exposure to contaminants, such as metals, hydrocarbons, and synthetic organic compounds, can lead to alterations in metabolism, hormone production, condition, and liver function in fish (Polgar et al. 1985; McArdle et al. 2004). The excessive loading of nutrients from sewage or septic wastewater, particularly nitrogen, can lead to eutrophication and subsequent hypoxic conditions that degrade fish habitat, especially when aided by thermal stratification (Nixon 1995). Thus, age-0 Winter Flounder in Long Island are likely subjected to numerous local stressors and may be especially vulnerable when these stressors combine with unfavorable temperatures and dissolved oxygen levels during summer months.

The objectives of this study were threefold: (1) to identify which individual attributes significantly and consistently affected patterns in otolith increment width (hereafter referred to as daily growth) in separate models for 2010 and 2011, (2) to construct and compare multivariable models to select the best model for each year, and (3) within the selected models, to examine and compare relationships between age, temperature, and individual attributes and daily growth in each study site. Central to this framework was the testing of 14 variables related to the ontogeny (age at capture, settlement date), condition (Fulton's K, hepatosomatic index [HSI], RNA/DNA ratio), and expression of genes associated with contaminant exposure (cytochrome P4501A [CYP1A], vitellogenin [VTG]), immune response (pleurocidin, hepcidin II, complement C3, phospholipase A2), and metabolic stress (glutamate decarboxylase, glucokinase, glycerol 3-phosphate dehydrogenase [GPDH]) within individual fish. We conducted an exploratory analysis, utilizing hierarchical linear models to investigate which of these factors influenced otolith-derived estimates of daily growth while controlling for age, daily average temperature, and location. Ultimately, we sought to provide insight into how individual variation in these ontogenetic, physiological, and molecular traits control patterns of growth, and possibly recruitment processes, in these vulnerable populations of juvenile fish.

METHODS

Study sites.—The sites evaluated in this study were Jamaica, Moriches, and Shinnecock bays (Figure 1), each of which are shallow inshore bays along the south shore of Long Island, New York. Although separated by less than 200 km, these sites lay along a decreasing west-to-east gradient of human population density (LIRPC 2010) and environmental conditions, primarily associated with differences in nutrient input and loadings of organic and inorganic pollutants (Bopp et al. 1993; Hartig et al. 2002). Analysis of priority pollutants in fine-grain sediments sampled from 2000 to 2006 as part of the U.S. Environmental Protection Agency's National Coastal Assessment generally supports a west-east gradient, with average concentrations of organic (PAHs and PCBs) and inorganic (Pb, Cu, Cd, and Hg) pollutants decreasing by at least a factor of five from west to east (USEPA 2006; McElroy et al. 2015). A previous study of Atlantic Silverside Menidia menidia in Long Island detected a significant longitudinal trend in sex ratios, such that female-biased sex ratios were more prevalent in more urbanized western embayments, which reflected a west-east gradient in municipal wastewater discharge (Duffy et al. 2009). In addition, environmental data collected from our three study sites in the summer of 2011 suggest that there is a gradient in bottom-water hypoxia, with the percentage of dissolved oxygen measurements below 2.3 mg/L being 23.0, 6.0, and 0.6% at Jamaica, Moriches, and Shinnecock bays, respectively (McElroy et al. 2015). Consequently, we anticipated that growth patterns in age-0 Winter Flounder would vary in accordance with this environmental gradient, with greater reductions in growth due to stress expected in more impacted western sites.

FIGURE 1.

Map showing the locations of the three sites sampled in Long Island, New York. An arrow denoting a gradient of anthropogenic impact reflects that Jamaica and Moriches bays are more impacted sites, while Shinnecock Bay is comparatively less impacted, with no direct sewage effluent and a lower human population density in its watershed. Sample locations consisted of areas within the bays known to support juvenile Winter Flounder and were randomly sampled (Frisket al. 2013).

Fish collections and basic metrics.—Winter Flounder were randomly sampled as a part of a larger study during a biweekly beam trawl survey from May through October in 2010 and 2011 (Frisk et al. 2013). Fish were flash frozen on dry ice after capture and brought back to the laboratory and stored in freezers at -80°C until they were processed. The total length (to nearest millimeter) and weight (whole body weight; to nearest 0.01 g) were recorded for each individual, and these measures were subsequently used to calculate condition indices. Fish from Shinnecock, Moriches, and Jamaica bays were included in this analysis if they were caught during time periods that allowed for 20 or more days of increment width data and corresponding age and temperature data to be collected (see below).

Otolith processing and analysis.—The left and right sagittal otoliths of each individual were removed and affixed to clean, labeled slides using heat-activated Crystalbond 509 adhesive. The left otolith was then polished whole, in the sagittal plane, using 30-µm and 3-µm lapping paper. The left otolith was selected for age and growth estimation because it has been observed to have a stronger relationship to somatic growth than the right otolith (Sogard 1991). After polishing, each otolith was immersed in Cargille Type A immersion oil and four coincident photographs were taken at 100× magnification using a Nikon DX1200C digital camera connected to a Nikon Eclipse 80i compound microscope. These images were merged to improve the depth of field and then analyzed using Image-Pro 6.0 software (Media Cybernetics 2006). An axis was placed along the postrostrum after the method of Yencho (2009). The axis started approximately 10 increments away from the formation of the secondary growth primordium because metamorphosis is known to be rapid and typically occurs in conjunction with the formation of secondary growth centers in this species (Sogard 1991). The axis was extended to the otolith edge so that increments could be counted back to the settlement period (Figure A.1 in the appendix). Age was estimated in each individual by counting the number of daily growth rings present, and the width between each increment was measured. Daily increment width is a valid proxy for daily growth in age-0 Winter Flounder in Long Island because otolith size (total width of the postrostral axis) and fish size have been observed to be linearly and strongly correlated (n = 126, r² = 0.77, P < 0.001 in 2010; n = 119, r² = 0.81, P < 0.001 in 2011; data not shown). Settlement date was back-calculated by subtracting the number of daily rings from the date of capture in ordinal days (day number).

Environmental data.—The modeling procedure required water temperature data in order to relate growth and temperature histories; however, environmental data were only continuously collected from the end of June through August of 2011 with sondes (Hydrolab MS5; Hach Hydromet, Loveland, Colorado) in each bay. To complete the summer temperature histories for each site, the sonde data from 2011 was combined with water temperature data obtained from a variety of sources and converted to daily averages for 2010 and 2011. In Jamaica Bay, daily bottom temperatures for all of May—August were obtained by combining conductivity, temperature, and depth sensor data collected by the New York City Department of Environmental Protection (NYCDEP 2010, 2011) and YSI Model 85 probe (Yellow Springs, Ohio) data obtained during the beam trawl survey used to collect fish. Temperature in Moriches Bay were recorded every 6 min by a tide and temperature gauge (Sea Bird Electronics SBE-26) at Smith Point, and data for all of May–August during both years was extracted from LIShore ( http://www.lishore.org; Thomas Wilson, Stony Brook University, personal communication). Continuous data from Shinnecock Bay were available from the beginning of June to the end of August in both 2010 (Christopher Gobler, Stony Brook University, unpublished data) and 2011 (Bradley Peterson, Stony Book University, unpublished data). There were gaps in the Jamaica Bay data set varying from 5 to 8 d in length, with one 14-d gap in 2011, and temperature during these time periods was inferred using linear interpolation. The resulting temperature data sets for both 2010 and 2011 spanned from June to August in Shinnecock Bay and from May to August in Jamaica and Moriches bays. For each fish, daily growth measurements were subsequently aligned with the corresponding age and daily average temperature during these time periods.

Condition indices and gene expression data.—Data detailing individual condition, gene expression, and microsatellite profiles at the time of collection were developed through other components of a larger study of Long Island Winter Flounder (Frisk et al. 2013). Fulton's K was calculated by dividing the body weight by the cube of total length for each individual and multiplying by 100, while HSI was calculated as the ratio of liver weight (to nearest 0.001 g) to body weight. Differential mRNA expression was examined in liver tissue from individuals to identify cellular pathways associated with the expression of two genes related to contaminant exposure (CYP1A and VTG), four associated with immune response (pleurocidin, complement C3, hepcidin II, and phospholipase A2), and three associated with glucose and glycogen metabolism (glucokinase, glutamate decarboxylase, GPDH). The raw data from the gene expression analysis spanned several orders of magnitude, so values were log_e transformed for all subsequent analyses. Full details of RNA isolation, cDNA synthesis, and quantitative PCR are provided in McElroy et al. (2015). Descriptions of all individual-level variables can be found in Table 1.

Hierarchical linear models.—Otolith increment widths, while useful in reconstructing growth patterns (Brothers et al. 1976; Campana and Neilson 1982), typically result in unbalanced, autocorrelated repeated-measures data (Campana and Jones 1992), which are often difficult to analyze from a statistical standpoint. In juvenile fish, growth variation can be attributed to factors operating over large scales, such as climatic and hydrographic regimes, as well as responses to local environmental conditions (Houde 1989; Sogard et al. 2001). Thus, factors affecting daily growth may be considered hierarchically structured (Figure 2). The three-level hierarchical linear models (HLMs) used in this study employed iterative maximum likelihood methods to estimate parameters, as well as variance and covariance components, in nested linear regression equations (Raudenbush and Bryk 2002). The hierarchical framework allowed data to be structured across three levels: level-1 factors that affected daily growth, level-2 factors that affected individual fish, and level-3 factors that affected all fish within each bay.

In HLMs, the levels are nested (i.e., growth increments within fish within bays) so that variation at one level is taken into account when calculating the variance components at any other level. In addition, level-2 and level-3 coefficients are estimated to explain variability in group means rather than individual level-1 responses, which addresses the correlation of error within groups when working with hierarchically structured data (McMahon and Diez 2007). One previous study utilized a hierarchical nonlinear modeling framework to investigate trends in annual size-at-age data in adult fish (Schaalje et al. 2002); however, no research, to our knowledge, has applied HLMs to describe daily growth patterns in juvenile fish.

TABLE 1.

List of the 14 candidate level-2 variables, with their category and a brief description. The last two columns denote how each variable performed during the first stage of the exploratory analysis each year. The abbreviation NS means the variable was not significant in at least four of the six models tested, and PI (shown in bold italics) means that the variable was considered potentially important because it had a P-value <0.01 in at least four of the six models tested; the parentheses denote the number of times each variable was significant for that year.

Continued.

The general form of a three-level HLM is provided in this section to fully conceptualize the hierarchical structure of this analysis. The notation is modified from Raudenbush et al. (2011), and regression coefficients are indicated by α, β, and γ, independent variables by x, y, and z, and random error terms by u, v, and w. At level 1, daily incremental growth G_ijk for increment i in fish j collected in bay k is represented as follows:

where α_pjk are level-1 regression coefficients, χ_pjk are level-1 independent variables, and the random error term u_ijk is assumed to be N(0, σ²). In the current study, age and temperature associated with each otolith increment were taken to be the level-1 independent variables.

At level 2, individual attributes, such as measurements of ontogeny, condition, and gene expression, were allowed to influence the level-1 coefficients:

where β_qpk are level-2 regression coefficients, y_qjk are level-2 independent variables, and the random error terms for individual j in bay k are ν;_pjk.

The final level, level 3, was used to represent bay-level characteristics:

where γ_rqp are level-3 regression coefficients, z_rk are level-3 independent variables, and the random error terms for bay k are W_qpk. In the present study, level-3 variables were categorical variables denoting different bays.

If equation (3) is substituted into equation (2), and then the terms from the resulting equation are substituted into equation (1), the result is a linear mixed-effects model in terms of γ_rqp. This final mixed model incorporates fixed and random factors on all three levels, in which parameters are estimated simultaneously. In this study, full maximum likelihood was used for estimation so that the deviance of models with different fixed factors could be directly compared (Kreft and de Leeuw 1998).

FIGURE 2.

Conceptual diagram of the hierarchical structure of daily growth data in age-0 Winter Flounder. The data are broken up into three levels, which are nested within one another. The effects of these variables on growth move up the hierarchy (arrows) such that the characteristics of each bay influence the individual attributes of fish, which in turn affect the daily growth patterns reflected in otolith increments.

Model selection.—An exploratory analysis was conducted using HLM7 (Raudenbush et al. 2011) in order to determine which independent variables and model structures best predicted daily growth in separate data sets for 2010 and 2011. To do so, the analysis was broken up into two steps. The level-1 independent variables, age and temperature, were retained throughout the selection process. This resulted in an intercept term, an age slope, and a temperature slope affecting daily growth (as described in the previous section). An interaction term between age and temperature did not substantially improve fit during either year and was thus not used in subsequent models. In addition, all level-2 variables were scaled and centered to their grand mean throughout the selection process.

Stage 1: Initially, there were 14 candidate explanatory variables at level 2 related to ontogeny, condition, contaminant exposure, immune response, and metabolism (Table 1) that were each tested independently in order to trim the number of explanatory variables. Each model incorporated a single level-2 variable affecting one of the three level-1 coefficients and level-3 site descriptors influencing the level-2 variable in one of two possible ways (i.e., through the intercept or slope). Random error terms were assigned to the intercept terms on level 2 (β_0pk) and level 3 (γ_0pk) in every model (i.e., random intercept model). Dummy variables denoting Shinnecock and Jamaica bays were used as level-3 variables because including all three sites resulted in collinearity. Thus, six models (i.e., three level-1 coefficients × two level-2 coefficients) were tested for each level-2 variable for each year, yielding 196 models overall. Variables progressed to the next stage of the analysis if they consistently had a significant influence on growth, such that the effect on level-1 coefficients had a P-value < 0.01 in four or more out of the six model runs conducted for each variable for each year. Six variables satisfied this criterion in 2010 while four variables did so in 2011.

Stage 2: In the final stage, various model configurations were compared to select the best models for each year; these models included different combinations of the selected level-2 variables affecting the level-1 coefficients (e.g., the intercept, age slope, and temperature slope). In addition, each level-2 effect was allowed to be influenced by site on level 3. The Akaike information criterion (AIC) was used to compare the relative fit of all model structures tested until the model with the lowest AIC was determined, selecting a single model for each year without creating over-parameterized models.

Model performance.—The performance of level-2 variables throughout stage 2 was reported as the proportion of total model runs that had a P-value < 0.05. Akaike weights (AIC_w) were used to evaluate the performance of models of varying complexity (i.e., number of parameters estimated), given the data (Wagenmakers and Farrell 2004). The AIC_w is estimated by first determining the difference between the AIC of the best model (minAIC) and that of all the candidate models (i):

The AIC_w is computed for all models with the following equation:

where n is the total number of models being compared. The AIC_w can be interpreted as a conditional probability and provides an estimate of how likely a model is relative to all candidate models (Wagenmakers and Farrell 2004).

Model performance was also investigated by constructing level-1 residual box plots for the model with the smallest AIC for each year (hereafter referred to as the selected model). These residuals represented the deviations, at level 1, of daily growth from the corresponding predicted values for each individual fish in the 2010 and 2011 models. Box plots for individuals were displayed in order of decreasing median increment width, with outliers and extreme values included. Intraclass correlations, which represent the proportion of total variance that was attributed to each level, were also calculated for the selected model for each year. Finally, model performance was evaluated by comparing observed growth-at-age data to predicted growth-at-age trajectories in each bay during each year. Ideally, the predicted trajectories would provide a realistic representation of the data and would fall within the range of observed growth values at a given age.

The selected models for 2010 and 2011 were analyzed graphically to illustrate the effects of each level-1 and level-2 variable on daily growth and to evaluate how these effects varied between the three sample sites. For each variable, we plotted the intercept, slope, and range of values at each site while holding all other variables constant at their mean. Contrasts were performed to (1) examine whether each slope was significantly different from 0 and (2) test whether the slopes of each relationship were significantly different among the three study sites.

RESULTS

Daily Growth

In total, otoliths from 57 fish were analyzed in 2010, consisting of 14 fish from Jamaica Bay, 21 from Moriches Bay, and 22 from Shinnecock Bay. In 2011, otoliths from 45 individuals were examined, including 17 individuals from Jamaica Bay, 16 from Moriches Bay, and 12 from Shinnecock Bay.

However, HLMs treat the structure of the variance in a manner that allows otolith increment widths to be considered independent when estimating degrees of freedom at level 1, such that 3,851 data points were used in the 2010 model and 2,158 data points were used in 2011. Otolith increment widths during the study period ranged from 1.96 to 15.01 µm in 2010 and from 2.12 to 12.53 µm in 2011. Average monthly otolith increment widths among all sites (Figure 3) throughout the whole growing season were variable and displayed similar trends during both years. Growth was highest early in the season, from February to April, and was followed by decreases in growth during the late spring through the summer months (i.e., those evaluated in the HLM analysis), before increasing again in the fall. Fish collected in our three study sites over the entire sampling period averaged 65.4 mm in length (range: 28–135 mm) and 116 d in age (55–238 d) in 2010, while averaging 59.9 mm in length (24–99 mm) and 101 d in age (38–185 d) in 2011.

FIGURE 3.

Plot of average otolith increment widths for all three sites in 2010 (black) and 2011 (gray), with error bars representing one standard deviation. The black dashed lines denote the time period analyzed in Moriches and Jamaica bays, while the gray dashed line denotes the start of the time period analyzed in Shinnecock Bay (which also ended in August), when daily growth was declining on average.

Hierarchical Linear Models

Model selection.—Stage 1: Model selection at this stage resulted in a reduction of the potential level-2 explanatory variables. In 2010, six explanatory variables satisfied the criterion of having four or more effects on daily growth with a P-value < 0.01 among the six runs for each variable: settlement date, age at capture, Fulton's K, HSI, and expression of pleurocidin and VTG. In 2011, four variables (settlement date, age at capture, and expression of CYP1A and GPDH) met the criterion.

Stage 2: The proportion of significant effects (P < 0.05) of level-2 variables in all potential 2010 models ranged from 0.98 for pleurocidin to 0.29 for Fulton's K in 2010. In 2011, the proportion of significant level-2 variables ranged from 0.80 for GPDH to 0.31 for CYP1A (Table 2). A total of 73 models in 2010 and 43 in 2011 were compared using AIC to determine the selected model for each year. The selected model in 2010 did not include VTG expression as a level-2 variable and resulted in five variables at level 2. Level 3 (site) influenced terms on level 2 in seven places, yielding a 25-parameter model. In 2011, the best model retained all four level-2 variables that passed through stage 1 of model selection, and level-3 variables influenced level-2 terms in six places, resulting in a 22-parameter model. Akaike weights based on the number of parameters estimated in all model configurations in both years (Table 3) generally showed that the selected models in both years were much more probable than the other model structures examined.

TABLE 2.

Table summarizing the performance of each level-2 variable tested in stage 2 of model selection. The Models row shows the number of models tested, the Significant row the number of model runs that were significant at the P < 0.05 level, and the Proportion row the proportion of significant model runs. Variables listed with an asterisk did not pass the stage-1 criterion and were not analyzed in stage 2.

2010 Selected-model contrasts.—A summary of all terms in the selected model for 2010 can be found in Table A.1 in the appendix. The intercept term was highly significant and positive in the selected model for 2010 (γ₀₀₀ = 7.703; P = 0.003), and this was the only term in the model which was not influenced by site. The age slope term was significantly negative in Shinnecock (γ₁₀₁ = -0.019; P < 0.001), Moriches (γ₁₀₀ = -0.023; P < 0.001), and Jamaica (γ₁₀₂ = -0.039; P < 0.001) bays, and the relationship in Jamaica Bay was significantly different from the other bays (P < 0.020). The temperature slope term was negative in Shinnecock (γ₂₀₁ = -0.101; P < 0.001) and Moriches (γ₂₀₀ = -0.041; P < 0.001) bays but was positive in Jamaica Bay (γ₂₀₂ = 0.040; P = 0.128), and all three sites were significantly different from one another P < 0.010). Settlement date was negatively related to daily growth in Shinnecock (γ₀₁₁ = -0.005; P = 0.003), Moriches (γ₀₁₀ = -0.062; P < 0.001), and Jamaica (γ₀₁₂ = -0.002; P = 0.459) bays, and Moriches Bay was significantly different than the other sites (P < 0.001). Age at capture significantly affected daily growth in Moriches Bay (γ₁₁₀ = -1.148; P = 0.001) but not in Shinnecock (γ₁₁₁ = 0.254; P = 0.595) or Jamaica (γ₁₁₂ = 0.160; P = 0.842) bays in 2010, and Moriches Bay had a significantly different relationship than Shinnecock Bay (P = 0.020). The expression of pleurocidin had a varying relationship with daily growth in Shinnecock (γ₁₂₁ = 0.006; P = 0.098), Moriches (γ₁₂₀ = -0.003; P < 0.001), and Jamaica (γ₁₂₂ = -0.018; P < 0.001) bays, and each site was significantly different from one another (P < 0.020). The condition indicator HSI was negatively related to daily growth in Shinnecock (γ₂₁₁ = -0.010; P = 0.431), Moriches (γ₂₁₀ = -0.019; P = 0.002), and Jamaica (γ₂₁₂ = -0.087; P < 0.001) bays, and this effect was significantly different in Jamaica Bay relative to the other sites (P < 0.001). The relationships between daily growth and Fulton's K were not significant at any site (P > 0.080), and no site was significantly different from any other (P > 0.100).

TABLE 3.

Akaike weights (AIC_W) for all models from stage 2 of the exploratory analysis, organized by the number of parameters estimated (i.e., model complexity).

2011 Selected-model contrasts.—A summary of all terms in the selected model for the 2011 data can be found in Table A.1 in the appendix. In 2011, the intercept term was again highly significant and positive (γ₀₀₀₀ = 8.875; P < 0.001), and this was once again the only term in the model not influenced by site. The age slope was also significantly negative in Shinnecock (γ₁₀₁ = -0.041; P < 0.001), Moriches (γ₁₀₀ = -0.042; P < 0.001), and Jamaica (γ₁₀₂ = -0.065; P < 0.001) bays, and Jamaica Bay once again had a different slope than the other sites (P < 0.001). Similar to 2010, the temperature slope term was negative in Shinnecock (γ₂₀₁ = -0.108; P < 0.001) and Moriches (γ₂₀₀ = -0.009; P = 0.589) bays but positive in Jamaica Bay (γ₂₀₂ = 0.019; P = 0.552), but Shinnecock Bay was significantly different than the other two sites (P < 0.010). Settlement date positively affected daily growth in Shinnecock (γ₀₁₁ = 0.697; P = 0.037) and Moriches (γ₀₁₀ = 0.694; P = 0.018) bays but was negative in Jamaica Bay (γ₀₁₂ = 0.388; P = 0.487), and none of these relationships were significantly different from one another (P > 0.080). Age at capture was positively related to daily growth in Shinnecock (γ₁₁₁ = 0.013; P < 0.001), Moriches (γ₁₁₀ = 0.017; P < 0.001), and Jamaica (γ₁₁₂ = 0.026; P < 0.001) bays, but no sites differed from one another (P > 0.050). The effect of CYP1A expression on daily growth varied in its direction in Shinnecock (γ₂₁₁ = -0.012; P = 0.495), Moriches (γ₂₁₀ = 0.044; P < 0.001), and Jamaica (γ₂₁₂ = -0.025; P = 0.038) bays, and Moriches Bay was significantly different than the other two sites (P < 0.010). In 2011, GPDH differentially affected daily growth in Shinnecock (γ₁₂₁ = 0.001; P = 0.514), Moriches (γ₁₂₀ = -0.008; P = 0.025), and Jamaica (γ₁₂₂ = -0.007; P = 0.072) bays, and Shinnecock Bay had a significantly different relationship than Moriches Bay (P < 0.010).

FIGURE 4.

Level-1 residual daily growth box plots for individual fish (sample sizes for each site in parentheses) in the selected models in 2010 (left column) and 2011 (right column). The gray-shaded boxes indicate the 25th, 50th, and 75th percentiles for individual fish, the error bars represent the 10th and 90th percentiles, the white dots denote outliers, and the black dots denote extreme values. Note the different scales of the y-axes (i.e., magnitude of variation) for each bay.

TABLE 4.

Intraclass correlation (ICC) estimates denoting the proportion of total variance attributed to level 1, level 2, and level 3 in the selected models for 2010 and 2011.

Model performance.—In general, box plots of level-1 residuals did not appear to display a strong trend in relation to the median observed increment width (Figure 4), and residuals were well centered around 0. However, there were differences in the range and the number of outliers and extreme values in each model for each year. In the selected models for both years, the intraclass correlation estimates indicated that over 90% of the total variance was attributed to level 1 and level 2 (Table 4), which suggests that the majority of variability in growth occurs at the daily growth (temperature and age) and individual fish (condition, gene expression, etc.) levels. Within each model, settlement date only had a significant effect in Shinnecock and Moriches bays, while gene expression and (or) condition metrics effects were significant in Moriches and Jamaica bays, the two more impacted sites in this study.

The majority of growth-at-age trajectories performed reasonably well during both years, with predicted values generally falling within the range of observed values (Figure A.2 in the appendix). There was a tendency, however, for the model to underestimate growth at early ages (<50 d), and this was particularly evident in Shinnecock Bay in 2010. In addition, there was one individual in Jamaica Bay (ID number 63) that had a consistently low predicted growth trajectory, eventually falling to 0 at an age of ∼110 d. This individual fell well below the range of observed values, but this may be attributed to the fact that all values of its level-2 variables were below average (and were thus negative) and also because of the comparatively late age and high temperature at which it entered the model.

The graphical analyses of the selected models in 2010 and 2011 reflected the results of the contrasts, wherein relationships between increment width and variables of interest were displayed as being significantly different from 0 or not. The effect of age was significantly negative during both years, while temperature did not have the strong negative effect on daily growth that was expected (Figures 5, 6). In addition, all level-2 variables, with the exception of Fulton's K in 2010, had significant effects on increment width in at least one of the study sites.

DISCUSSION

Daily growth in age-0 Winter Flounder was analyzed using HLMs and showed significant relationships with ontogenetic factors, environmental conditions, and indices of health and condition. The selected models for 2010 and 2011 each contained a different set of several level-2 variables, thereby supporting previous studies that have demonstrated that relationships between growth and environmental conditions in this species are driven by multiple factors (Manderson et al. 2002). Interestingly, temperature did not have a strong direct effect on growth over the range of observed values (13–27°C; Figure A.3 in the appendix). Age at capture was negatively related to growth for Moriches Bay in 2010 and was positively related in all three locations in 2011. Settlement date was negatively related to daily growth in Shinnecock Bay and Moriches Bay but was positively related to daily growth in 2011 in both bays. Indicators of fish health (HSI in 2010) and immunological and metabolic stress (pleurocidin in 2010, CYP1A and GPDH in 2011) only had significant effects on daily growth in more impacted sites (Jamaica and Moriches bays). Our findings suggest that the optimum conditions for growth are spatially and temporally dynamic and that daily growth cannot be reduced to a single factor; instead, growth patterns are a complex response to a highly variable environment.

The switch from a negative to a positive relationship between settlement date and growth suggests that annual variation in environmental conditions influences the success of recruits settled during different periods of a given year. In general, early settlement corresponds with cooler temperatures, late settlement would favor higher temperatures, and both would have an impact on the bioenergetics and growth of juveniles (Bullock 1955). However, the association of daily growth with temperature was generally weak, suggesting other factors likely have an important influence on the role settlement date plays in determining growth. Previous research has explained lower growth rates in juvenile Winter Flounder with higher temperature (Meng et al. 2000) and a greater range of temperature (Meise et al. 2003); however, the former also found a positive relationship with temperature, suggesting additional factors are likely involved. Other factors that could have additive effects on growth include higher temperatures resulting in higher feeding rates (Fonds et al. 1992) or increased prey abundance (Franz and Tanacredi 1992). Alternatively, older and larger individuals that settled early in the season may have advantages in terms of resiliency to stressful environmental conditions and predation (Johnson and Hixon 2010) or in terms of outcompeting smaller conspecifics (Shepherd and Cushing 1980). These observations highlight the trade-offs that age-0 Winter Flounder encounter, and the optimal settlement date likely changes annually depending on the temperature, food resources, size reached by a given date, and ability to survive unfavorable conditions, such as hypoxia. These results also suggest that adult spawning behaviors that result in multiple cohorts of young of year in the same year-class (Yencho 2009) may have historically enabled Long Island Winter Flounder to produce viable year-classes in the face of uncertain environmental conditions.

FIGURE 5.

Plots of site-specific effects of level-1 and level-2 variables on the predicted otolith growth (y-axis) relationship in the selected 2010 model. Solid lines correspond to fish from Shinnecock Bay, hatched lines correspond to fish from Moriches Bay, and dashed lines correspond to fish from Jamaica Bay. Black lines denote slopes that are significantly different than 0 (P < 0.05), while lines in gray are not significant; PLEUR = pleurocidin.

The selected models during both years also contained five other variables: Fulton's K, HSI, GPDH, pleurocidin, and CYP1A. We expected that the condition indices (HSI, Fulton's K) would display a positive relationship with daily growth (Bolger and Connolly 1989), while the gene expression metrics (pleurocidin, CYP1A, and GPDH) would be more complex, with the potential for a mixed response depending on exposure to contaminants and potential compensatory responses, such as heightened immune response or mobilization of energy reserves. Fulton's K was included in the final 2010 model selected but did not have significant slopes for any site in 2010. The other condition index examined, HSI, had a significant negative relationship with daily growth for the more impacted Jamaica and Moriches bays in 2010. The HSI may be reflective of the energetic status and ration of individual fish, and it is known to be negatively affected by exposure to metals in age-0 Winter Flounder (Pereira et al. 1993).

Cytochrome P4501A mRNA expression, enzymatic activity, or protein levels are a common biomarker of exposure to organic pollutants, responding to many contaminants found in wastewater and urban runoff, such as many polycyclic aromatic hydrocarbons and a variety of planar chlorinated hydrocarbons, including some PCBs, dioxins, and furans (Stegeman and Hahn 1994; Schlenk et al. 2008). However, as reactive metabolites produced through CYP1A-mediated metabolism can have toxic consequences, and as induction has been linked mechanistically to a variety of developmental deficits, induction of CYP1A is sometimes considered a biomarker of toxic effect. The data presented here is among the first to link a physiological meaningful response (growth) to CYP1A induction. The fact that the relationship between CYP1A and growth is negative in Jamaica Bay Winer Flounder and positive in those from Moriches Bay, however, suggests that other factors may also be influencing CYP1A expression and its relationship to growth.

Glycerol 3-phosphate dehydrogenase is an enzyme involved in glycolysis whose expression has been associated with the mobilization of glycogen reserves in European Flounder Platichthys flesus exposed to hypoxia (Jørgensen and Mustafa 1980). The negative relationship between daily growth and GPDH observed in Moriches Bay in 2011 may therefore be the result of exposure to low dissolved oxygen in these fish, although this relationship was not significant in Jamaica Bay, another site known to experience intermittent hypoxic and hyperoxic conditions (Frisk et al. 2013). Pleurocidin is a peptide isolated from the mucous secretions of Winter Flounder that possesses antimicrobial properties (Cole et al. 1997). In 2010, pleurocidin mRNA expression was significant in the model, but the relationship was positive for Jamaica Bay and negative in Moriches Bay. The difference in the direction of this effect is consistent with adaptive (immune stimulation) and maladaptive (immune suppression) responses occurring in Jamaica and Moriches bays, respectively (Romany et al. 2015). The significant relationships between growth and pleurocidin expression in Moriches and Jamaica bays in 2010 may reflect the potential importance of immune status in controlling growth patterns in age-0 Winter Flounder, especially in more impacted sites.

FIGURE 6.

Plots of site-specific effects of level-1 and level-2 variables on the predicted otolith growth (γ-axis) relationship in the selected 2011 model. Solid lines correspond to fish from Shinnecock Bay, hatched lines correspond to fish from Moriches Bay, and dashed lines correspond to fish from Jamaica Bay. Black lines denote slopes that are significantly different than 0 (P < 0.05), while lines in gray are not significant.

Compensatory growth may explain the increased daily growth rates that were observed in fish with a reduced condition (HSI in 2010) and an increased gene expression associated with environmental stress (pleurocidin in 2010, CYP1A and GPDH in 2011). Although compensatory growth is likely an important driver of growth in wild fish, there is a paucity of field-based research on the subject (Ali et al. 2003). Laboratory experiments have confirmed that compensatory growth in age-0 Winter Flounder is possible. Bejda et al. (1992) exposed three treatments of fish to low, high, or fluctuating dissolved oxygen levels for 10–11 d, and growth rates were significantly lower in conditions of low or fluctuating oxygen. Following initial treatments, all individuals were placed in higher oxygen conditions for 5 weeks and individuals from the hypoxic treatment grew significantly faster than those from the other groups. This suggests that while prolonged exposure to hypoxic conditions can suppress growth rates, these same fish may be able to alter their consumption rates or reallocate energy reserves to increase short-term growth and compensate for adverse environmental conditions. These experiments may provide an explanation for some of the unexpected growth patterns observed in Jamaica Bay and Moriches Bay, the two most impacted sites. However, compensatory growth is not without trade-offs and it has been associated with decreased survival (Munch and Conover 2003) and various short- and long-term costs (Mangel and Stamps 2001; Mangel and Munch 2005). The populations of Winter Flounder in this study are in a depleted state, and previous research suggests they are suffering from a population bottleneck occurring during the first year of life (Yencho 2009; Frisk et al. 2013); the potential for compensatory growth may be another indicator of a highly stressed population.

Distinct responses in growth to attributes related to health (HSI) and stress (gene expression metrics) in Jamaica and Moriches bays may be due in part to population-specific responses. O'Leary et al. (2013) detected weak, but significant, fine-scale genetic differences in age-0 Winter Flounder in Long Island bays. Similarly, recent and historical studies indicate that Winter Flounder display diversity in migratory and spawning behaviors and the existence of resident individuals (Poole 1966; Sagarese and Frisk 2011; Frisk et al. 2014). Black et al. (1988) documented that PCB contamination can be passed down from mother to offspring in Winter Flounder from New England, and it may be possible for this process to elicit phenotypic and genetic changes over many generations. Over time, this may result in fish from more anthropogenic ally impacted sites becoming either adapted or acclimated to with-stand higher contaminant loadings (Wirgin and Waldman 2004). Understanding how different local populations of Winter Flounder respond and adapt to such long-term stressors is also important for rebuilding the species in Long Island.

Hierarchical linear models provided a useful framework to evaluate complex relationships between important traits, such as growth, and a myriad of potential explanatory variables. The model structure and statistical technique provided a means to analyze data measured across various scales, including longitudinal daily growth data (Chambers and Miller 1995; McMahon and Diez 2007). This type of data structure is likely to arise in many situations when investigators are interested in evaluating growth in juvenile fish, especially if repeated measurements are taken from fish collected from multiple locations. The HLMs performed reasonably well, capturing growth in age-0 Winter Flounder. In some individuals, very high growth early in the season was not always captured by the models, especially in Shinnecock Bay in 2010. Future research exploring alternative model structures and additional explanatory variables may help better capture individuals with high growth early in the season. Our study would have benefited from more consistent and longer-term measurements of temperature, dissolved oxygen, and salinity, as these and other potentially important level-1 and level-2 factors have been observed to influence growth in age-0 Winter Flounder. A combination of controlled laboratory experiments and additional field-based research is needed to explore the temporal variation in the influence of potential drivers of growth, the timing of settlement patterns, and the relationship between somatic growth rate and otolith increment width. In addition, a better understanding of compensatory growth and its triggers is important for understanding observed growth trajectories in impacted sites.

In conclusion, daily growth patterns in age-0 Winter Flounder are complex, being related to settlement patterns, condition, and gene expression associated with contaminant exposure, immune response, and energy metabolism. The dynamic response of growth with settlement date and other level-1 and level-2 factors suggests that variability in the timing of spawning within bays may also play an important role in controlling these processes. A full life cycle research approach is needed to link spawning periodicity, larval stage duration, timing of settlement, and recruitment success in order to elucidate drivers of Winter Flounder population dynamics and recovery in Long Island bays (Frisk et al. 2014).