Spatio-temporal patterns of the oceanic conditions and nearshore marine community in the Mid-Atlantic Bight (New Jersey, USA)

Study area and survey design

To evaluate marine community dynamics, the survey area (Fig. 1) was divided into 15 sampling strata (12–26) by depth and geographical location (latitude and longitude). To be consistent with established federal marine resource field-sampling programs, the designers of the New Jersey Ocean Stock Assessment (OSA) survey incorporated the same latitudinal boundaries defined by the National Marine Fisheries Service (NMFS), Northeast Fishery Science Center (NEFSC), Northeast Atlantic Groundfish Survey Program; exceptions were those strata at the northern and southern ends of the New Jersey coastline where NMFS extended its survey into New York and Delaware waters (ASMFC, 1994). The boundaries were also truncated in the northern and southern strata to include only the waters adjacent to the New Jersey coastline and the ocean waters off Delaware Bay. The longitudinal boundaries consisted of the 9.1 (30 ft.), 18.3 (60 ft.), and 27.4 m (90 ft.) isobaths. The bottom contours were somewhat irregular, so the stratum boundaries were smoothed using GIS techniques (ASMFC, 1994).

To reduce potential sampling bias, each sampling area was divided into smaller blocks. Mid-shore blocks (9.1–18.3 m) and offshore (18.3–27.4 m) blocks were 2.0 min longitude by 2.5 min latitude, whereas nearshore (5.5–9.1 m) blocks were 1.0 min longitude by 1.0 min latitude. Nearshore block dimensions were smaller because the strata were narrower and encompassed a smaller area than the mid and offshore strata; thus, the smaller block size permitted a greater number of potential sampling sites than would be possible with larger dimensions. It should be noted the blocks truncated by stratum boundaries encompassed a smaller area (>50%) than the whole blocks (Byrne, 1994; Byrne, 2008).

Experimental field sampling approach

Field sampling was conducted bimonthly (every two months: February, April, June, August, October, and December) from 1988 to 1989. From 1990 to date, the December and February surveys were replaced by a single winter survey in January, followed by surveys in April, June, August, and October (ASMFC, 1994) The annual sampling survey effort during 1988 through 1990 varied slightly because of the limited budget (e.g., high charter vessel costs), but annual sampling generally consisted of two hauls per stratum. Overall, the sampling effort averaged 39 hauls (i.e., two samples from each strata plus one additional haul in each of the nine larger strata) per survey. The average number of stations sampled each year was around 182.

Station selection

Constrained randomization was used to select unique sampling stations for each survey trip. Sampling stations (survey site location) were randomly selected by the NJDEP program leader during 1988 through 1991, but this method was replaced in 1992 by a computer generated random number selection program. Because stratum shapes were elongate and the sampling effort was limited, a station selection procedure was used to reduce any spatial distribution sampling bias. The station selection procedure consisted of limiting the first station to only the top half of the block numbers and the second station to the bottom half; however, if a third station was selected then these limitations were not imposed in the procedure process. For instance, haul one would be selected from blocks 1 to 25, haul two from blocks 26 to 50, and haul three from blocks 1 to 50 for a stratum with 50 blocks. For each station, three additional alternate sites were also pre-selected using the same procedures described above to account for any fixed fishing gear (e.g., traps or nets), bottom obstructions, or other impediments that prevented sampling at the initial station Byrne, 1994; Byrne, 2008.

Field sampling gear

Field sampling was conducted with a tapered (forward to rear netting) three-in-one bottom otter trawl. The forward (i.e., wings and belly) and rear netting was constructed with #3/20 twisted polyethylene twine. The forward and rear netting was constructed with 12 cm and eight cm stretch mesh, respectively. The cod-end was constructed with 7.6 cm stretch mesh four mm polyethylene twine and lined with #147 style (6.4 mm stretch mesh) white knitted nylon netting. The round and drop-mesh corners were hung to yield a 50 percent hanging ratio. The trawl doors were 2.44 m × 1.27 m (3.1 m²) constructed with marine grade wood and steel shoes that weighed approximately 453.5 kg. Based on hydroacoustic sensors, the average wing spread was around 13 m. The headrope was 25 m and the footrope was 30.5 m. The top and bottom bridles were 36.6 m and constructed with a 1.27 and 1.91 cm wire rope, respectively. The groundwire was 18.3 m. The bottom bridle and groundwire were covered with 6.03 cm rubber cookies. The rigging (sweep) consisted of 7.62 cm rubber cookies on 14.3 mm wire rope, 9.5 mm long drop chains (six link long), and 32 center-hole floats (20.32 cm). The extension chain was a 1.27 cm Trawlex chain Byrne, 1994; Byrne, 2008.

Field sampling procedures

All tows were conducted between sunrise and sunset. The trawl tow duration was standardized at 20 min and the vessel speed was between 4.7 and 5.6 kilometers per hour (2.5 and 3.0 knots). The swept area (a) was estimated with the following equation: a = D∗hr∗X2, D = V∗t; where V is the velocity of the trawl over the ground when trawling, hr is the length of the head-rope, and t is the time spent trawling. X2 is that fraction of the head-rope length, hr, which is equal to the width of the path swept by the trawl, the “wing spread”, hr∗X2. Based on vessel speed, one 20 min tow generally covered a distance of 1.85 km. Given the trawl dimensions and distance towed, the total swept area was around 24,050 m². The cardinal direction of the tow was determined by the oceanic conditions (wind, waves, and current) at the time of deployment; tows were generally made in the direction of the waves, wind, and current. A 91.5 m wire was used to maintain a tow depth ratio of approximately 3:1 Byrne, 1994; Byrne, 2008. Survey replicates (tow) were considered independent of one another given the random station selection process, and the distance between sampling sites and time between tows; each tow was considered a random sample of the population.

Data collection protocol

At each sampling station the surface and bottom environmental conditions (water temperature (degrees Celsius (°C), salinity (parts per thousand [psu], and dissolved oxygen (milligrams per liter (mg/l)) were measured with a CTD and recorded before deploying the trawl. After the 20 min tow was completed, the trawl was retrieved and the catch (fish and macroinvertebrates) was rough sorted into plastic buckets. Afterwards, the entire catch was identified to species, enumerated, and the length (fork and/or total, as appropriate) was measured to the nearest cm for fish (20 individuals randomly selected); the disk width (cm) was measured for skates and rays. All species were identified to the lowest taxa. The total weight was taken using either a hanging or floor scale. The individual weights of every species were determined by weighing individual baskets (total weight) of every species collected and dividing by the total count of individuals in the basket. Various other measurements were recorded depending on the macroinvertebrate species. For example, the carapace width (mm) was measured for crabs, the carapace length (mm) for lobsters and mantle length (mm) for squids. Because some catches were too large to sort in the field, a representative thoroughly mixed sub-sample was randomly selected and weighed. After the sub-sample was sorted, species composition was extrapolated to determine the total catch (Byrne, 1994; Byrne, 2008).

Data treatment/processing

To minimize any potential spatial non-independence, data were pooled among stations within each individual sampling area. Before initiating statistical hypothesis tests, environmental and biological data were transformed (e.g., logarithmic, square root, fourth root, or arcsine) to meet normality assumptions, and down-weight the statistical effects (i.e., reduce skewness) of abundant taxa, while allowing less common taxa to contribute to sample discernment (Thorne, Williams & Cao, 1999; Korsman, 2013). Normal probability plots were examined, and Kolmogorov–Smirnov and Bartlett tests were used to assess normality and homoscedacity (Zar, 1999). Outlier observations were investigated to determine whether the outlier occurred by chance; all outliers were retained for these analyses.

To evaluate the nearshore marine community and oceanic conditions, 28 years (1988–2015) of fishery-independent monitoring data (environmental and biological) were compiled, sorted (time [year and month] and space [area and zone]), and summarized. After pooling the data by stations sampled within each area, the marine community (catch characteristics [total number, estimated abundance, and estimated biomass]) was evaluated using two approaches: a single dataset (pooling all the data) and segregating the data in various time-series datasets. The data were segregated into six 5-year time-series periods to help discern patterns and test for potential differences among time and space. This approach was driven from the perspective that most available time-series data for nearshore/offshore fisheries are only two to five years in duration. Data were pooled by stations and segregated by individual sampling areas (12–26) and geographical north/south zones defined as the following: 1 (sampling areas 12–17), 2 (sampling areas 18–23), and 3 (sampling areas 24–26). Data was also segregated by west/east zones defined as the following: 1 (sampling areas 12, 15, 18, 21, and 24), 2 (sampling areas 13, 16, 19, 22, and 25), and 3 (sampling areas 14, 17, 20, 23, and 26). It should be noted the amount of area for each designated geographical zone was a different size in terms of km².

Following Howell & Auster (2012), marine species were classified a priori as coldwater-adapted species (primarily distributed in cold temperate regions), warmwater-adapted species (primarily distributed in warm temperate regions), or subtropic-adapted species (primarily distributed in subtropical and tropical regions). Classification followed Froese & Pauly (2018) and published life-history literature (e.g.,  Murdy, Baker & Musick, 1997; Collette & Klein-MacPhee, 2002; Able & Fahay, 2010) describing a species’ distribution relative to the MAB, water temperature tolerance (minimum and maximum), preferred water temperature range, and preferred spawning water temperature. In general, the mean preferred water temperature was used to select the best water temperature preference group for each species. Species preferring water temperature <15 °C were generally classified as coldwater-adapted, while those preferring water temperatures 15–29 °C were classified as warmwater-adapted. Species preferring temperatures >30 °C were classified as subtropic-adapted.

The statistical significance level was defined as P < 0.05. In the presence of significance at the 95 percent confidence level, a Tukey post-hoc multiple pairwise comparison tests were used to differentiate the specific differences among the population means. Data were evaluated using various software, including Microsoft Access^®, Microsoft Excel^®, and Statgraphics Centurion XVI^®.

Statistical analyses

Physicochemical conditions

Salinity

The mean annual surface salinity off the coast of New Jersey within 15 strata (areas: 12–26) during 1988 through 2015 ranged from 29.94 in 1996 to 32.11 psu in 1998 with a mean of 30.75 psu (±2.06 psu). Overall, the mean salinity decreased about 0.13 psu per year or 1.3 psu per decade. The surface salinity varied significantly among years, and there was a weak negative (non-significant) association between the surface salinity and time (Table 1). Averaging over years, the mean monthly surface salinity ranged from 30.25 in April to 31.48 psu in November with a mean of 30.73 psu (±2.07 psu). The surface salinity varied significantly among months (F [11, 5084] = 31.46, P < 0.05).

The mean surface salinity ranged from 28.92 in sampling area 14 to 31.45 psu in sampling area 26 with a mean of 30.74 psu (±2.07 psu). The surface salinity varied significantly among sampling areas (Table S1), and there was a strong positive (significant) association between the surface salinity and sampling area (Table S2). The surface salinity was generally lower in the northern sampling areas than the southern sampling areas (Table S1 ). The surface salinity ranged from 28.13 psu in sampling area 12 to 31.49 psu in sampling area 26. The lowest surface salinity levels were found in the areas closest to shore, and the highest surface salinity levels were detected in areas furthest from shore.

A paired t-test showed the mean annual bottom salinity level was higher than the mean annual surface salinity (t (5095) = −34.25; P < 0.05). The lowest (30.52 psu) bottom salinity was in 2005 and the highest (32.87 psu) was in 1988 with a mean of 31.66 psu (±1.4 psu). The bottom salinity varied significantly with time (F [27, 5069] = 50.12, P < 0.05), and there was a weak negative (significant) association between the surface salinity and time (Table 1). The bottom salinity ranged from 31.09 in December to 31.94 psu in November. Averaging over years, the bottom salinity also varied significantly among months (F [11, 5085] = 5.13, P < 0.05).

The mean bottom salinity ranged from 30.33 in sampling area 12 to 32.32 in sampling area 23 with a mean of 31.66 psu (±1.4 psu). The bottom salinity varied significantly among sampling areas (Table S2), and there was a weak negative association between the bottom salinity and the sampling area (Table S2). Segregating the sampling areas into zones (See Methods), the mean bottom salinity varied significantly among northern and southern zones (F [2, 5094] = 14.31, P < 0.05); a post-hoc test showed mean bottom salinity varied significantly between zones 1 and 2 (−0.18), and 2 and 3 (0.15). The mean bottom salinity level increased from northern to middle sampling areas, and then decreased in the southern sampling areas. The mean bottom salinity varied significantly among eastern and western zones (F [2, 4829] = 373.5, P < 0.05); a post-hoc test showed mean bottom salinity varied significantly between zones 1 and 2 (−0.54), 1 and 3 (−1.07), and 2 and 3 (−0.52). Mean bottom salinity increased from eastern to western zones.

Dissolved oxygen

The mean surface DO off the coast of New Jersey within 15 strata (areas: 12–26) during 1988 through 2015 ranged from 7.9 in 1995 to 8.98 in 1998 mg/L with a mean of 8.54 mg/L (±1.52 mg/L), and there was a weak positive association between the surface DO and time (Table 1). Overall, the mean DO increased about 0.009 mg/L per year or 0.09 mg/L per decade. Averaging over years, the mean monthly surface DO ranged from 7.17 in September to 10.27 mg/L in December. The surface DO varied significantly among months (F [11, 5080] = 608.75, P < 0.05).

The mean surface DO ranged from 8.20 in sampling area 24 to 9.06 mg/L in sampling area 14. The surface DO varied significantly among sampling areas (Table S2), and there was a strong negative association between the surface DO and the sampling area (Table S2). In general, the surface DO decreased from northern to southern sampling areas, but no pattern was evident for the sampling areas closest and furthest from shore (Table S1).

A paired t-test showed the annual bottom DO levels were significantly lower (t (5091) = 132.17; P <0.05) than the annual surface DO levels suggesting a strong water column stratification. The mean bottom DO off the coast of New Jersey within 15 strata (sampling areas: 12–26) during 1988 through 2015 ranged from 6.64 in 1988 to 8.35 mg/L in 1993 mg/L with a mean of 7.58 mg/L (±1.84 mg/L), and there was a weak positive association between the bottom DO and time (Table 1). The mean monthly bottom DO ranged from 5.55 in August to 9.87 mg/L in January. Averaging over years, the bottom DO varied significantly among months (F [11, 5081] = 1121.85, P < 0.05).

The mean bottom DO ranged from 7.09 in sampling area 16 to 8.13 mg/L in sampling area 24. The surface DO varied significantly among sampling areas (Table S2), and there was a strong positive association between the bottom DO and sampling area (Table S2). Segregating the sampling areas into zones (See Methods), the mean bottom DO varied significantly among northern and southern zones (F [2, 5090] = 54.78, P <0.05); a post-hoc test showed mean bottom DO varied significantly between zones 1 and 2 (−0.42), 1 and 3 (−0.69), and 2 and 3 (−0.27). The bottom DO generally increased from northern to southern sampling areas. The mean bottom DO also varied significantly among western and eastern zones (F [2, 4825] = 2.77, P < 0.05); a post-hoc test showed mean bottom DO varied significantly between zones 1 and 3 (0.15). The mean DO generally decreased from nearshore to offshore sampling areas.

Marine fauna community

Annual and spatial dynamics

Over 20.4 million fish and invertebrates (1,338.3 mt) representing 214 (water temperature preference classified) species (not including unidentified species) were collected within 15 strata (areas: 12–26) off the coast of New Jersey from 1988 to 2015. Three marine fauna water temperature preference groups (coldwater-adapted, warmwater-adapted, and subtropic-adapted) were identified in the study area (Fig. S1). The total number of individuals collected ranged from 4.0 million (coldwater-adapted) to 6.8 million (subtropic-adapted). In each temperature preference category, three species represented the majority of the catch. The main coldwater-adapted species collected were longfin squid (Loligo pealei) (n = 2,225,975), Atlantic herring (Clupea harengus) (n = 544,032), and little skate (Leucoraja erinacea) (n = 316,356), while Atlantic butterfish (Peprilus triacanthus) (n = 2,873,138), scup (Stenotomus chrysops) (n = 1,318,569), and northern searobin (Prionotus carolinus) (n = 503,230) represented the warmwater-adapted group. Bay anchovy (Anchoa mitchilli) (n = 9,227,960), striped anchovy (Anchoa hepsetus) (n = 245,214), and Atlantic moonfish (Vomer setapinnis) (n = 38,691) denoted the subtropic-adapted group.

The most abundant water temperature preference group was the subtropic-adapted and the coldwater-adapted group was the least abundant group in the study area. The number of individuals collected per group varied over time with the coldwater-adapted group slightly decreasing since 2004 (Fig. 2). In terms of percent composition, coldwater-adapted group declined over time, but in some years (1989, 1990, 1991, 1993, 1995, 1997, 2006, and 2012) it was the second most dominant group (Fig. S2). In 1998, the coldwater-adapted group was the most abundant group. Overall, the coldwater to warmwater-adapted ratio declined over time, and there was a weak negative association between the coldwater: warmwater-adapted ratio and time (F [1, 26] = 10.71, P = 0.003; r² = 29.2%) (Fig. 3)

Pooling all samples, the lowest (0.1784 individuals per m²) mean estimated abundance (# individuals/m²) was the warmwater-adapted group and highest (0.2340 individuals per m²) was the subtropic-adapted group (F [2, 100201] = 618.47, P < 0.05). The lowest mean estimated biomass (kg/m²) was the subtropic-adapted group and the highest was the coldwater-adapted group (F [2, 100201] = 687.35, P < 0.05).

The mean estimated abundance of the coldwater-adapted group was consistently lower than the estimated abundance of the warmwater-adapted group over time (Fig. 4). The mean annual estimated abundance varied significantly by water temperature preference category (F [2, 100203] = 587.86, P < 0.05), time (F [27, 100203] = 11.93, P < 0.05), and the interaction between the water temperature group and time (F [54, 100203] = 9.33, P < 0.05). Overall, pooling the six 5-year time-series, the mean annual estimated abundance of the coldwater-adapted (F [5, 39039] = 8.82, P < 0.05) and subtropic-adapted (F [5, 7382] = 8.16, P < 0.05) groups decreased, and the warmwater-adapted group slightly increased (F [5, 39039] = 14.10, P < 0.05) over time.

The mean estimated abundance of the water temperature groups significantly varied by time (F [2, 100203] = 600.08, P < 0.05), space (F [14, 100203] = 3.93, P < 0.05), and the interaction between time and space (F [28, 100203] = 3.82, P < 0.05). The highest mean estimated abundance for the coldwater-adapted group was in the northern sampling areas (sampling areas 13, 14, and 12), while the highest mean estimated abundance for the warmwater-adapted (sampling areas 20, 18, and 19) and subtropic-adapted (sampling areas 22, 20, and 19) groups was in mid and southern sampling areas.

Similarly, the mean annual estimated biomass varied significantly by water temperature preference category (F [2, 100203] = 705.8, P < 0.05), time (F [27, 100203] = 6.49, P < 0.05), space (F [14, 100203] = 2.92, P < 0.05), and the interaction between time and space (F [54, 100203] = 3.79, P < 0.05). The mean annual estimated biomass for all three water temperature preference groups (coldwater-adapted [F [5, 39039] = 38.84, P < 0.05]; warmwater-adapted [F [5, 39039] = 44.03, P <0.05]; and subtropic-adapted [F [5, 7382] = 10.35, P < 0.05] increased with time, but the mean annual estimated biomass of the subtropic-adapted group increased the most and coldwater-adapted group increased the least. The interaction pattern between the mean annual estimated biomass and space was less clear than the mean annual estimated abundance and time; there was no pattern in the mean annual estimated biomass of the coldwater-adapted [sampling areas 18, 13, and 23], warmwater-adapted [sampling areas 20, 12, and 19], and subtropic-adapted [sampling areas 18, 12, and 22] groups.

Pooling the data (1988–2015), GLMs showed that year and sampling area were significant predictors of the total number and estimated abundance of the coldwater and warmwater-adapted groups. The total number and estimated abundance of the subtropic-adapted group were significantly predicted by the year, month, and sampling area, whereas the estimated abundance was only significantly predicted by the year and sampling area (Tables 2–4).

Table 2:

General Linear Model and associated ANOVA Type III Sums of Squares.

Coldwater-adapted species Temperature Preference. Five-yr Time-series (1988–2015). CW, Coldwater-adapted species.

Dependent variable	Source	Sum of squares	Df	Mean square	F-ratio	P-value	Fitted model	R2
Total number (Coldwater-adapted)	Time Series (CW)	11.3453	1	11.3453	7.99	0.0047	√√ Total Number (CW) = 6.3662–2.13748E–7*Time Series (CW) + 0.0192003*Month (CW) - 0.0104509*Area (CW)	0.41%
	Month (CW)	136.547	1	136.547	96.15	0.0000
	Area (CW)	74.3637	1	74.3637	52.36	0.0000
	Residual	55435.6	39036	1.42011
	Total (corrected)	55663.8	39039
Estimated abundance (Coldwater-adapted)	Time series (CW)	0.0292921	1	0.0292921	2.38	0.1229	√√ Estimated Abundance (CW) = 0.413077–1.0861E–8*Time Series (CW) + 0.00168064*Month (CW)–0.000971917*Area (CW)	0.36%
	Month (CW)	1.0462	1	1.0462	85.01	0.0000
	Area (CW)	0.643144	1	0.643144	52.26	0.0000
	Residual	480.408	39036	0.0123068
	Total (corrected)	482.161	39039
Estimated biomass (coldwater-adapted)	Time Series (CW)	0.0336629	1	0.0336629	6.63	0.0100	√√ Estimated Biomass (CW) = −0.0790099 + 1.16431E–8*Time series (CW) - 0.00145997*Month (CW) − 0.00139156*Area (CW)	1.06%
	Month (CW)	0.789501	1	0.789501	155.47	0.0000
	Area (CW)	1.31841	1	1.31841	259.62	0.0000
	Residual	198.236	39,036	0.00507829
	Total (corrected)	200.357	39,039

DOI: 10.7717/peerj.7927/table-2

Table 3:

General Linear Model and associated ANOVA Type III Sums of Squares.

Warmwater-adapted species Temperature Preference. Five-yr Time-series (1988–2015). WW, Warmwater-adapted species.

Dependent variable	Source	Sum of squares	Df	Mean square	F-ratio	P-value	Fitted model	R2
Total number (Warmwater-adapted)	Time Series (WW)	60.3819	1	60.3819	40.24	0.0000	√√ Total Number (WW) = −8.22361 + 4.92926E-7*Time Series (WW) + 0.0550105*Month (WW) −0.00515665*Area (WW)	1.66%
	Month (WW)	922.387	1	922.387	614.75	0.0000
	Area (WW)	18.0308	1	18.0308	12.02	0.0005
	Residual	58,570.3	39,036	1.50042
	Total (corrected)	59,559.4	39,039
Estimated abundance (warmwater-adapted)	Time Series (WW)	0.780301	1	0.780301	60.04	0.0000	√√ Estimated Abundance (WW) = −0.967161 + 5.6035E-8*Time Series (WW) + 0.00501172*Month (WW)–0.000490182*Area (WW)	1.64%
	Month (WW)	7.65589	1	7.65589	589.07	0.0000
	Area (WW)	0.162928	1	0.162928	12.54	0.0004
	Residual	507.331	39036	0.0129965
	Total (corrected)	515.812	39,039
Estimated biomass (warmwater-adapted)	Time Series (WW)	0.800626	1	0.800626	164.43	0.0000	√√ Estimated Biomass (WW) = −1.0566 + 5.67601E–8*Time Series (WW) + 0.00400921*Month (WW) + 0.000216194*Area (WW)	2.91%
	Month (WW)	4.89938	1	4.89938	1006.22	0.0000
	Area (WW)	0.0316933	1	0.0316933	6.51	0.0107
	Residual	190.069	39,036	0.00486908
	Total (corrected)	195.747	39,039

DOI: 10.7717/peerj.7927/table-3

Segregated spatial dynamics

West/east spatial dynamics

The mean annual estimated abundance of subtropic-adapted species was highest and the coldwater-adapted species were lowest in every west/east zone; west/east zones were designated as the following: west (sampling areas 12, 15, 18, 21, and 24); mid (sampling areas 13, 16, 19, 22, and 25), and east (sampling areas 14, 17, 20, 23, and 26). The mean annual estimated abundance of subtropic-adapted species decreased from west (Zone 1) to east (Zone 2) Zones, and coldwater and warmwater-adapted species increased from west to east zones.

The mean annual estimated abundance of coldwater-adapted species varied significantly by time (F [27, 40368] = 2.07, P = 0.0009) and space (F [2, 40368] = 4.81, P = 0.008). The mean annual estimated abundance in Zone 1 was the highest in 1993 and 2004. In Zone 3, the mean annual estimated abundance peaked in 1990 and 2011.

The mean annual estimated abundance of warmwater-adapted species varied significantly by time (F [27, 40368] = 2.09, P = 0.002) and interaction between time and space (F [42, 40368] = 1.64, P = 0.005). The mean annual estimated abundance peaked in Zone 3 during 1992. In Zone 2, the mean annual estimated abundance peaked in 1994, 1996, 1999, 2002, 2004, and 2006.

The mean annual estimated abundance of subtropic-adapted species varied significantly over time (F [27, 7471] = 1.90, P = 0.003). The mean annual estimated abundance in Zone 1 peaked in 1990, 2000, 2002, and 2014. The mean annual estimated abundance in Zone 2 peaked 1992, 1995, 2008, and 2013. In Zone 3, the mean annual estimated abundance peaked in 1993, 2009, and 2011.

Overall, the highest mean annual estimated biomass was warmwater-adapted species, but there was no evidence to suggest the mean annual estimated biomass among temperature preference categories significantly changed from western to eastern areas.

The mean annual estimated biomass of coldwater-adapted species did not vary significantly by time (F [27, 40368] = 1.21, P = 0.21) or space (F [2, 40368] = 2.90, P = 0.0552). The mean annual estimated biomass in Zone 3 peaked in 2000, and in Zone 2 it peaked during 2007; the mean annual estimated biomass in Zone 1 was relatively low throughout the time period.

Similarly, the mean annual estimated biomass of warmwater-adapted species did not vary significantly by time (F [27, 40368] = 1.09, P = 0.35) or space (F [2, 40368] = 0.35, P = 0.89). In Zone 1, the mean annual estimated biomass peaked in 2012, and peaked in Zone 2 in 2005. In Zone 3, the mean annual estimated biomass peaked in 1992, 2001, and 2007.

The mean annual estimated biomass of subtropic-adapted species did not vary significantly by time (F [27, 7471] = 0.31, P = 0.99) or space (F [2, 7471] = 0.03, P = 0.97). In Zone 1, the estimated biomass peaked in 1991, and increased from 2012 to 2015. In Zone 2, the mean annual estimated biomass peaked in 2004 and 2012. The mean annual estimated biomass in Zone 3 was relatively low during most of the years, but it peaked in 1997 and 2014.

Physicochemical conditions

Identifying the annual and seasonal variability in the oceanic conditions and the response of community and component populations within the ecosystem is critical for predicting long-term community dynamics, trends, and evaluating disturbance. Given the broad oceanographic hydrodynamics off New Jersey (Kohut, Glenn & Chant, 2004), the water temperature, salinity, and DO levels significantly varied over time and space in the study area. The mean oceanic conditions were highly variable over the 28-year period, but various alternating or cyclic patterns were evident, along with increasing or decreasing trends. In general, water temperature increased and the salinity decreased over time. The rising water temperature and falling salinity trends echoed previous studies for the region (Howell & Auster, 2012; Geiger et al., 2013; Thomas et al., 2017). The trend rates for mean (0.06 °C per year) and maximum (0.02 °C per year) surface water temperature were similar to published rates (Thomas et al., 2017).

The nearshore waters off New Jersey are a dynamic hydrological system influenced by summer stratification and winter mixing, which are associated with the prevailing wind and buoyancy factors (Glenn et al., 2004; Kohut, Glenn & Chant, 2004). Northeast wind often causes downwelling, while southwest wind causes upwelling (Kohut, Glenn & Chant, 2004). The physicochemical conditions varied not only by year and season, but among specific sampling areas within the study area. Besides the southeast corner of the study area (sampling areas 21 and 24), water temperature (surface and bottom), salinity (surface), and DO increased from northern to southern sampling areas, and salinity increased from western to eastern sampling areas. Bottom salinity decreased from middle to southern areas, and the surface DO decreased from northern to southern areas. Bottom DO increased from northern to southern sampling areas, and decreased from eastern to western sampling areas.

Upwelling and downwelling off New Jersey are frequent events given the wide continental shelf and gently sloping topography; the continental shelf extends about 200 km off the coast of New Jersey (Song, Haidvogel & Glenn, 2001; Kohut, Glenn & Chant, 2004). It is possible the relatively warmer water temperature and higher salinity in sampling areas 21 and 24 was influenced by the nearby underlying topography (i.e., topographic bump and oblique sand ridges) or seasonal upwelling events. The nearshore waters (∼40 m) are more stratified than the offshore waters in June because of the influence of lower salinity water from the Hudson River (Schofield et al., 2008), which could be causing the decreasing salinity trend off the New Jersey coast, noted here. The warmwater period in the MAB is longer now than in the past and this environmental change has influenced the stratified period (June–September) thereby altering the timing of spring and fall phytoplankton blooms. Over the past decade, the warming period is beginning earlier each year (Thomas et al., 2017).

The surface and bottom DO varied significantly over time and the positive association between DO and time was explained adequately by regression, but a low correlation coefficient value (0.08) indicated a weak association. It is difficult to explain why the DO slightly increased (0.009 mg/L per year) with time given research has shown several recurrent hypoxia events (1994, 1996, and 2001) occurring along the southern New Jersey coast (Glenn et al., 2004). Hypoxia is relatively common along the New Jersey coast, but these events are somewhat short-term and related to coastal upwelling, which sometimes occurs in summer when the wind is from the southwest (Glenn et al., 2004). The hypoxia centers (∼150 km²) are spatially isolated in duration (∼1 week), frequency (∼5 times in 9 years), and space (Barnegat Inlet, Mullica River Estuary, and Townsend/Hereford Inlets). The hypoxia locations are downstream of a series of topographic highs associated with ancient river deltas in the southern waters off New Jersey. The most significant upwelling events occurred after the most severe cooling seasons in 1994, 1996, and 2001 (Glenn et al., 2004). Researchers hypothesized that severe cooling seasons often causes colder and larger Cold Pools, which produce more significant summer upwelling events by summertime wind-driven forces. Upwelling also depends upon wind, precipitation, and storm frequency. Given that water temperature is increasing in the study area, it is possible that summer upwelling events could be less severe than in the past, which is reducing the magnitude of hypoxia events in specific areas within the study area. This hypothesis might explain the slight increase in DO over time. Then again, it might be related to the ongoing water quality improvements (less pollution) in nearby New York and New Jersey waterbodies (HydroQual, 2010).

Marine fauna community

The nearshore coastal waters off New Jersey provide year-round and seasonal habitat for three temperature preference groups (subtropic-adapted, warmwater-adapted, and coldwater-adapted), which complemented previous findings in the MAB (Wood, Collie & Hare, 2009; Howell & Auster, 2012). However, different than Howell & Auster (2012), subtropic-adapted species were the most abundant and coldwater-adapted species were the least abundant in the study area. Bay anchovy was the most abundant subtropic-adapted species and butterfish was the most abundant warmwater-adapted species. These findings did not agree with Howell & Auster (2012) who reported more coldwater-adapted species than warmwater-adapted and subtropic-adapted species in Long Island Sound. New Jersey waters are not only located south of Long Island Sound, but anecdotal information from local fishermen suggest warmwater eddies and whorls spinning off the Gulf Stream Current can provide habitat to subtropic-adapted marine fauna, which could explain why there were more subtropic-adapted species in New Jersey than New York.

The estimated abundance of warmwater-adapted species is increasing, coldwater and subtropic-adapted species is decreasing, and the coldwater to warmwater-adapted ratio is decreasing over time. This ongoing pattern for marine fauna seems to be becoming more common in the MAB (Nye et al., 2009; Wood, Collie & Hare, 2009; Howell & Auster, 2012), and throughout the world (e.g.,  Polovina et al., 2011). Most researchers attribute this shift in species distribution and composition to climate variability (Cheung et al., 2011; Jang et al., 2011; Polovina et al., 2011; Brander, 2013). The abundance of warmwater-adapted species is also increasing with time in Narragansett Bay (Wood, Collie & Hare, 2009). In nearby Long Island Sound, Howell & Auster (2012) reported a shift from a coldwater to warmwater-adapted dominated species community, and an increase in subtropic-adapted species over a 25-year duration (1984–2008). It is possible the present study did not show this same trend for subtropic-adapted species because the environmental conditions off New Jersey are more variable than in the Long Island Sound (a semi-enclosed estuary) given the frequency of upwelling and downwelling events (Glenn et al., 2004; Kohut, Glenn & Chant, 2004). Though the mean water temperature is rising, the frequency of upwelling/downwelling is also increasing in the MAB (Kohut, Glenn & Chant, 2004), which could explain the decline in estimated abundance of subtropic-adapted species over time. It is also possible subtropic-adapted species cannot tolerate abrupt changes in water temperatures caused by these events. Thus, these oceanographic events could be negatively impacting some subtropic-adapted species since the influx of cold water/warm water (±1 −4 °C) can often occur in a short (∼3 or 4 weeks) period (Kohut, Glenn & Chant, 2004).

The highest estimated abundance of coldwater-adapted species was in northern sampling areas (12–17), and the highest estimated abundance for warmwater and subtropic-adapted species was in mid and southern sampling areas (18–26). The total number of coldwater and warmwater-adapted species decreased from north to south, but the estimated abundance did not vary significantly among the study area. Explaining the distribution of species is not straightforward in terms of water temperature (Wood, Collie & Hare, 2009), as it appears to vary by time, space and season given the oceanographic dynamics in the study area. Despite this variable pattern, the estimated abundance of coldwater-adapted species was generally higher in the southern sampling areas (24–26) during earlier years (1990s) and higher in the northern sampling areas (12–14) during later years (2000s). The estimated abundance of warmwater-adapted species was more variable, but it was higher in the southern sampling areas during a few earlier years (1990s) and higher in the northern sampling areas during a few later years (2000s). The estimated abundance of subtropic-adapted species was highest in the middle sampling areas (18–23), followed by the northern (12–17) and southern sampling (24–26) areas. In general, the highest estimated abundance was higher in the southern sampling areas (24–26) during a few earlier years, and higher in the northern sampling areas (12–17) during later years. In some ways, the findings suggest the distribution of species, based on their temperature preference, is currently in a transitional phase; the period when some individuals can still tolerate the mean water temperature. For instance, the water temperature preference classification for this study was based on the mean preferred water temperature, which means 50 percent of the individuals can tolerate either a lower or higher water temperature. Assuming the water temperature continues to rise with time, then a full transition from a coldwater dominated community to a warmwater community will occur over time.

The estimated abundance of coldwater and warmwater-adapted species increased from nearshore to offshore and subtropic-adapted species decreased from nearshore (12, 15, 18, 21, 24) to offshore (14, 17, 20, 23, 26) sampling areas. The increase in the estimated abundance of warmwater-adapted species from nearshore to offshore sampling areas suggests the environmental conditions are ideal for a geographical range expansion within the study area given the rising water temperature. In general, the estimated abundance of coldwater-adapted species was higher in the nearshore sampling areas in a few earlier years (1990s) and higher in the offshore sampling areas in a few later years (2000s). These findings suggest the distribution of coldwater and warmwater-adapted species is shifting north, which agrees with other findings in the MAB (Nye et al., 2009; Howell & Auster, 2012).

The estimated abundance for all three temperature preference categories has changed over time. The estimated abundance of subtropic-adapted species decreased from nearshore to offshore sampling areas, and coldwater and warmwater-adapted species increased from nearshore to offshore sampling areas. Despite these general patterns, inter-annual patterns were challenging to decipher given their high annual population and habitat selection variability. For instance, coldwater-adapted species had high abundance in the nearshore sampling areas during 1993 and 2004, and high abundance in the offshore zone during 1990 and 2011. The findings partially support the hypotheses that the estimated abundance and biomass of warmwater and subtropic-adapted assemblages increased from offshore to nearshore sampling areas, and the estimated abundance and biomass of coldwater-adapted assemblage increased from nearshore to offshore sampling areas with time. Overall, it appears that biomass trends are idiosyncratic to specific areas within the overall study area, but these findings could be more related to how the sampling was spatially segregated rather than to biological reasons.

The estimated biomass for all three water temperature preference categories is increasing with time in the study area; however, the estimated biomass of coldwater-adapted species is increasing the most even though they are decreasing in abundance throughout the study area. Larger individuals seem to be replacing smaller individuals, or less abundant larger individuals are more common in the study area, such as bullnose sting ray (Dasyatis sayi). It is difficult to explain why the estimated biomass of subtropic-adapted species is increasing the most since their estimated abundance is decreasing with time. Smaller or younger (juvenile life-stage) coldwater-adapted individuals could be declining, moving away, or maybe the larger less abundant subtropic-adapted species are moving into the study area to feed. Actually, it is possible that juvenile/sub-adults are moving offshore to deeper colder waters, while adults are moving nearshore in pursue of prey. Another potential explanation is that rising water temperature is causing coldwater-adapted species to grow faster (Duffy et al., 2016). It is often difficult to explain or generalize the response of fish populations to climate change given the number of influential factors and individual species broad responses (Rijnsdorp et al., 2009).

The estimated biomass of coldwater and warmwater-adapted species was influenced by month and sampling area, whereas subtropic-adapted species were generally more influenced by the sampling area. This observation seems reasonable since the estimated abundance of most species in the coldwater and warmwater-adapted groups is correlated with season (i.e., water temperature). However, it was somewhat surprising that the subtropic-adapted group was not influenced by season, but it is possible that certain sampling areas have ideal and less fluctuating water temperatures given the oceanic dynamics (upwelling/downwelling and eddies) in the study area (Kohut, Glenn & Chant, 2004). The total number and estimated abundance of coldwater-adapted and warmwater-adapted species in recent years (2013–2015) were influenced more by month and sampling area than the previous years (2008–2012) suggesting the seasonal water temperature could be rising and falling faster depending on the sampling area. Actually, the warming period is occurring earlier and lasting longer (Thomas et al., 2017), which could explain the decreasing and increasing abundance in coldwater-adapted and warmwater-adapted groups, respectively.