Spatiotemporal patterns in the ecological community of the nearshore Mid-Atlantic Bight

: Recognition of the need for a more holistic, ecosystem approach to the assessment and management of living marine resources has renewed interest in quantitative community eco - logy and fueled efforts to develop ecosystem metrics to gain insight into system status. This investigation utilized 12 years (2008 to 2019) of fisheries-independent bottom trawl survey data to quantify and synthesize the spatiotemporal patterns of species assemblages inhabiting the near-shore Mid-Atlantic Bight (MAB). Assemblages were delineated by ecomorphotype (EMT), and all species collected by the survey were allocated among 9 EMTs: demersal fishes; pelagic fishes; flatfishes; skates; rays; dogfishes; other sharks; cephalopods; and benthic arthropods. Annual time series and seasonal spatial distributions of relative aggregate biomass were quantified for each EMT using delta-generalized additive models. Dynamic factor analysis (DFA) revealed that the information content of the 9 annual time series was effectively summarized by 3 common trends, and DFA model fits to each EMT time series represented a new suite of ecosystem indicators for this system. Mean sea surface temperature during winter in the MAB was included in the selected DFA model, suggesting that winter environmental conditions influence the structure of this sys-tem at an annual scale. Principal component analysis uncovered a north-to-south gradient in the seasonal spatial distributions of these EMTs and identified a distinct area of elevated biomass for several assemblages along the south shore of Long Island, NY. Taken together, these results characterize the community structure of the nearshore MAB and yield requisite information to support ongoing ecosystem-scale assessment and management activities for this region.


INTRODUCTION
Societies depend on marine ecosystems as a source of both market goods and non-market services, and the magnitude and diversity of the demands on these systems have increased as the human population and associated economies continue to expand (Holland et al. 2010, Tam et al. 2017. Concomitant with this growth, rising anthropogenic pressures have altered the structure and function of living marine resource communities inhabiting these systems, often jeopardizing sustained delivery of desired economic out-puts (Link et al. 2020). For example, land use practices have increased turbidity, contaminants, and the extent of hypoxia in many adjacent marine ecosystems, disrupting systemic habitat utilization patterns and predator−prey dynamics (Craig 2012, Buchheister et al. 2013, Ortega et al. 2020. Climate change has altered the spatiotemporal distributions of marine species at broader scales as taxa shift poleward and to deeper waters in response to warming, which has likely impacted population vital rates, phenologies, and interactions (Nye et al. 2009, Doney et al. 2012, Langan et al. 2021. Meanwhile, harvest activ-ities have had perhaps the most direct and demonstrable effects on ecological communities through overfishing and sequential exploitation of target stocks (Essington et al. 2006), bycatch of non-target species (Kennelly 1995), fisheries-induced evolution (Walsh et al. 2006), and habitat alteration due to fishing practices (Collie et al. 1997). In response, governments have implemented fisheries management frameworks to reduce harvest and promote sustainable utilization of these resources.
Efforts to curb overfishing and rebuild overfished stocks in the USA have achieved measurable successes since the passage of the Fishery Conservation and Management Act in 1976, and have relied primarily on single-species approaches to assess the status of living marine resources and develop regulations for the associated fisheries (Methot et al. 2014, National Marine Fisheries Service 2022. A notable shortcoming of this reductionist approach, however, is that while stock productivity is molded by an integrated series of drivers that broadly represent biophysical forcing, trophodynamics, and exploitation, single-species assessment and management have typically only considered the influence of exploitation (Link 2010, Gaichas et al. 2012. As harvest is constrained, natural environmental processes, trophic interactions, and non-fishing human impacts that include the aforementioned examples play a greater role in shaping population dynamics (Tyrrell et al. 2011, Fogarty 2014. Further, the breadth of current assessment and management activities primarily encompasses species targeted by fisheries, such that the structure, function, and dynamics of the broader ecological community remain unresolved even in many well-managed systems. The recognized need for a more holistic, integrated approach to the stewardship of fisheries resources has fueled a renewed interest in quantitative community ecology designed to support development of Integrated Ecosystem Assessments (IEA) and operationalization of Ecosystem Approaches to Fishery Management (EAFM) for Large Marine Ecosystems (LMEs) under US jurisdiction (Levin et al. 2014, Link et al. 2020, Muffley et al. 2021. The Mid-Atlantic Bight (MAB) is the southernmost Ecological Production Unit (EPU) in the Northeast USA LME, and includes the continental shelf and estuarine ecosystems bounded by Cape Cod, MA, and Cape Hatteras, NC (Lucey & Fogarty 2013). Given the broad and increasing intra-annual temperature ranges that occur in this region (Forsyth et al. 2015) and rather complex hydrography defined primarily by the Gulf Stream, remnants of the Labrador Current, and freshwater inputs to the estuaries, this EPU is a dynamic and productive system that supports a diverse ecological community and numerous fisheries. There is a wealth of data on living marine re sources inhabiting the MAB, given the array of fisheries-independent monitoring surveys that have been sampling the shelf and estuaries of this system for decades (Howell & Auster 2012, Politis et al. 2014, Langan et al. 2021. Early efforts to quantify community dynamics in the offshore waters (i.e. > 27 m isobath) of the shelf revealed distinct species assemblages that exhibited notable seasonal and spatial dynamics (Gabriel 1992), and these assemblages have been shifting primarily northeastward since the 1960s (Lucey & Nye 2010, Kleisner et al. 2016. Many faunal groups inhabiting the shelf have also experienced appreciable changes in relative abundance through time (Methratta & Link 2006, Northeast Fisheries Science Center [NEFSC] 2022), while assemblages in Narragansett Bay (Collie et al. 2008, Langan et al. 2021), Long Island Sound (Howell & Auster 2012), and Chesapeake Bay (Jung & Houde 2003, Buchheister et al. 2013) have exhibited distinct spatiotemporal patterns shaped by natural and anthropogenic drivers. Taken together, these data on shelf and estuarine communities have formed the basis for the development of ecological metrics, particularly time-series indicators, that yield insight into the status of system structure and productivity, which in turn facilitate development of IEAs for the region (Muffley et al. 2021).
While the ecological communities of the estuarine and outer shelf regions of the MAB EPU have been the focus of intense research efforts, the structure and dynamics of species assemblages in the coastal ocean (i.e. < 27 m isobath) have yet to be explored in a holistic, synthetic manner, given the lack of directed and consistent sampling throughout this zone prior to 2008. The inner shelf represents a valuable foraging, spawning, and nursery area for an array of species, and serves as a transitional ecosystem between offshore and estuarine environments (Malek et al. 2014). The physical oceanography of this region is dynamic, influenced by estuarine outflows, seasonal upwelling events, and rapid temperature changes characteristic of shallow water environments, and this system experiences additional direct and indirect anthropogenic disturbances comparable with those observed in estuaries, due to the proximity to land (Farr et al. 2021). Given the unique nature of this coastal ocean, documented spatial variability of community-based, ecosystem indicators within the MAB EPU (Heim et al. 2021), and the recent estab-lishment of a fisheries-independent trawl survey in this nearshore zone, an evaluation of the community dynamics in this ecosystem could yield valuable information needed to support ecosystem assessment efforts for the MAB.
As such, this investigation used 12 years of trawl survey data collected from the nearshore MAB to (1) quantify relative annual, seasonal, and spatial trends in aggregate biomass for several multispecies groups; (2) answer relevant ecological questions on the degree of coherence among the annual time series, trends in productivity of these assemblages, and the influence of biophysical forcing and exploitation in shaping these trends; and (3) assess patterns in ecosystem usage through synthesis of the seasonal spatial distributions of these groups. This approach provides insight into the dynamics of the ecological community in the coastal ocean of this EPU, yields several ecological indicators that reflect temporal patterns of system status, and generates baseline spatiotemporal information against which to evaluate future system states.

Field sampling
Data for this investigation were collected by the Northeast Area Monitoring and Assessment Program (NEAMAP) Nearshore Bottom Trawl Survey from 2008 to 2019. NEAMAP is a fisheries-independent monitoring program designed to sample the late juvenile and adult stages of the living marine resources inhabiting the MAB coastal ocean ( Fig. 1) through annual spring (April to May) and fall (September to November) cruises (Atlantic States Marine Fisheries Commission 2009). The sampling frame encompasses 12 097 km 2 and is bounded by the 6.1 and 18.3 m depth contours between Montauk, NY and Cape Hatteras; although sampling in the vicinity of Delaware Bay extends to 27.4 m to maintain spatial continuity. Survey activities in Block Island Sound and Rhode Island Sound occur from 18.3 to 36.6 m depth. Sites are selected from a grid of 772 ha sampling cells using a stratified-random design, where stratification is delineated by latitudinal or longitudinal regions and depth (6.1−12.2, 12.2−18.3, 18.3−27.4, and 27.4−36.6 m), and effort is allocated in proportion to the surface area of each stratum. NEAMAP sampled at 150 sites during each cruise, except during spring 2017 when only 63 sites were sampled.
A 3-bridle, 4-seam bottom trawl with a 48 m (circumference) fishing circle and 2.54 cm knotless lined codend is towed during daylight hours at 1.5 m s −1 for 20 min at each sampling location. NEAMAP uses a net mensuration system to collect trawl geometry data and ensure consistent gear performance on each tow. Vessel position is also recorded throughout trawl hauls and used to calculate tow distance, which is coupled with wingspread data to estimate the area swept by the trawl at a site (range: 0.012 to 0.033 km 2 ). While a suite of additional variables is recorded synoptically with each haul, the most relevant for this investigation were sample time of day and water depth (to the nearest 0.05 m). Catch data are collected on each species captured at a site, and include aggregate weight (to the nearest 2 g) and count, as well as length data (to the nearest 0.5 cm)

Species assemblages
Species assemblages can be defined using an array of approaches, but often have been structured based on taxonomy, habitat associations, life history, morpho logy, or trophic mode. Here, species were grouped based on coarse similarities in taxonomy, mor pho logy, and habitat usage, and thus generally represent ecomorphotypes (EMT; sensu Compagno 1990). In designating these EMTs, consideration was given to maintaining alignment with the primary modes of fishing in this ecosystem and consistency with assemblage definitions used in previous studies of the Northeast USA LME. All taxa sampled by NEAMAP were allocated among 9 EMT assemblages that included 3 bony fish, 4 elasmobranch, and 2 invertebrate groups (Table S1 in the Supplement at www.int-res.com/ articles/suppl/m704p015_supp.pdf), and each species was assigned to an assemblage based on information in the literature (Murdy et al. 1997, Collette & Klein-MacPhee 2002. This qualitative approach was used to delineate EMTs in lieu of quantitative methods, such as cluster analysis or multivariate ordination, based on concerns that using NEAMAP catch data to both partition these assemblages and subsequently quantify the spatiotemporal patterns of these groups would represent a circular analysis. All bony fishes collected by the survey were partitioned among the demersal fish, pelagic fish, and flatfish categories. These EMTs separated fishes that primarily inhabit the water column from those typically associated with the seafloor, divided the latter taxa based on general body form, and broadly aligned with the fishery categories in the MAB (Gaichas et al. 2014). Further, the assemblage designations for these bony fishes mirrored groupings used in prior studies that focused on quantifying the community structure of and developing ecosystem indicators for the offshore waters of the Northeast USA LME , Methratta & Link 2006.
Elasmobranch species were grouped into the skate, ray, dogfish, and other shark EMTs. Note that while skates and rays share morphological similarities in that component species of both groups are dorso ventrally compressed, several skate taxa are harvested by directed commercial fisheries while rays are usually encountered as bycatch in this eco-system (NEFSC 2020). The decision to evaluate skates and rays as separate EMTs was based on these fishery differences, and because temporal patterns of the skate assemblage had been previously assessed for Georges Bank, an adjacent EPU north of the MAB (Fogarty & Murawski 1998). Similar reasoning supported the decision to separate the dogfish and other shark EMTs. Dogfish included smaller-bodied sharks and, as with skates, the historical dynamics of this assemblage have been evaluated for Georges Bank (Fogarty & Murawski 1998). The species composing this EMT are also subjected to directed fisheries (NEFSC 2006, Southeast Data, Assessment, andReview [SEDAR] 2015). In contrast, the other shark category encompassed medium-to-larger bodied ani mals that are either targeted by smaller-scale directed harvest, encountered as bycatch, or for which possession is prohibited (Mandelman et al. 2008). Finally, the invertebrate taxa collected by NEAMAP were allocated among 2 EMTs that separated softbodied cephalopods from benthic arthropods ).

Spatiotemporal modeling of ecomorphotype relative aggregate biomass
Annual trends and seasonal spatial distributions of relative aggregate biomass were quantified for each EMT assemblage by employing model-based approaches commonly used to standardize fisheriesindependent survey catch-per-unit-effort (CPUE) data (Maunder & Punt 2004). Site-specific CPUE for an assemblage was given by: (1) where B i,k is the biomass (kg) of the i th species of the EMT collected at sampling site k, S represents the number of taxa included in that assemblage, and a k is the area swept by the trawl.
Generalized additive models (GAMs) were used to relate CPUE data for each group to a suite of covariates recorded synoptically at the NEAMAP sampling sites. The GAM framework represents a flexible univariate approach that allows a response variable to be modeled as a function of both parametric (similar to generalized linear models) and non-parametric components, where the latter incorporate smoothing functions for continuous covariates and thus permit a nonlinear relationship between the response and predictor variable that is informed by the available data (Wood 2017). Covariates considered in these GAMs included: year (categorical), season (categorical), depth (m, continuous), solar zenith (angle, continuous), and a season−coastal position interaction (km, continuous) that accounted for the seasonally varying spatial distributions of each EMT. The coastal position variable was created given the need for a spatial metric, and because latitude and longitude are highly colinear within the NEAMAP survey area (Fig. 1). The first step in creating this variable was to generate a near-continuous contour that mirrored the MAB coastal shoreline. Given that the inshore boundary of the NEAMAP sampling frame is comprised of sampling cells that follow this shoreline, the landward 3 positions of each of these inshore cells (i.e. south-, central-, and northwest for cells south of New York Harbor, and northwest, -central, and -east otherwise) were selected from the full survey sampling grid. Linear regression models were then fitted to consecutive pairs of these points, 1000 location predictions were generated between each pair, and predictions were concatenated to yield a near-continuous contour of approximately 1.4 m resolution at the inshore margin of the sampling frame that paralleled the coastal shoreline ( Fig. 1). This approach was akin to applying a LOESS smoothing spline to these inshore positions using a very small span width such that a curve (i.e. straight line) is fitted between each pair of points, resulting in a contour that tracks the shoreline (Zuur et al. 2007). Next, all sites sampled by NEAMAP were assigned to the nearest location on this contour. Finally, the distances along the contour from Cape Hatteras (the 0 km reference) to these assigned locations were calculated and used to represent the coastal positions (distances) of the sampling sites (sensu Nye et al. 2009). All distance calculations were performed using the 'geosphere' package in R (v3.6.3, R Core Team 2020).
Some trawl hauls resulted in no catch for a given assemblage, and thus most of the EMT datasets contained several CPUE observations equal to zero. Initial attempts to fit GAM models for each group using the Tweedie distribution, which has positive mass at zero, were unsuccessful based on an evaluation of diagnostic plots, which was likely due to an insufficient number of zero observations in the EMT datasets (Shono 2008). A delta-GAM framework was therefore used to model CPUE of these assemblages (Lo et al. 1992). Delta-GAMs are hurdle models suitable for continuous response data that include zero values and contain 2 components: a binomial model fitted to presence/absence data to estimate the probability of encountering the assemblage and a conditional component fitted to CPUE data for the subset of sites where the EMT was encountered. Delta-GAM models were expressed as: ( 2) where represents the probability of encountering the EMT, E(y) is the expected value of the response vector CPUE for the subset of sites where the assemblage was captured, X is the fixed-effects model matrix containing observations of the categorical covariates, is the vector of fixed-effects coefficients, s j is the smoothing function for continuous covariate j, and g is the monotonic link function. Representatives of the demersal fishes and skate EMT assemblages were collected at approximately 99% of NEAMAP sites and thus were modeled using the conditional component of Eq.
(2) only, where a small quantity (0.2 kg for demersal fishes, 10 kg for skate) was added to all catches given that the domains of the available response distributions were defined on the positive real line.
The delta-GAM models were fitted for each EMT using generalized additive models for location, scale, and shape (GAMLSS) regressions (Rigby & Stasinopoulos 2005). GAMLSS extends traditional GAM models by incorporating distributions of the response variable beyond the exponential family. For each EMT, the most supported distribution of the CPUE data in the conditional component of the model (or overall model for demersal fishes and skates) was identified by fitting a model with all candidate covariates (i.e. full model) using each available continuous distribution defined on the positive real line and evaluating competing forms using Akaike's information criterion (AIC; Akaike 1973, Burnham & Anderson 2002. Once the most supported conditional distribution was identified, 4 forms of both the binomial and conditional models were fitted (Table 1), and the most supported parameterization within each component was identified by AIC and evaluation of model diagnostic plots. Note that analyses of residuals yielded no evidence of appreciable spatial (correlogram) or temporal (partial autocorrelation function) autocorrelation in these EMT datasets (Zuur et al. 2009).
While the covariates included in the delta-GAM models for each assemblage were allowed to vary be tween the binomial and conditional components, given that the processes that influence probability of encounter and conditional relative abundance may differ (Lo et al. 1992, Rubec et al. 2016, the year covariate and the season-coastal position interaction were included in all candidate models, since an objective of this investigation was to evaluate annual trends in and the seasonal spatial distributions of relative biomass for each EMT. Predictions of annual and seasonal−spatial probability of encounter and conditional CPUE were generated using marginal means (Searle et al. 1980), and these component values were multiplied to yield both an annual index of relative abundance and to estimate the seasonal spatial distribution for each EMT. For the latter metric, CPUE was predicted at 1000 points spaced evenly over the domain of the coastal position contour (~1 km resolution) for each season. Coefficients of variation (CVs) for these annual and seasonal spatial predictions were quantified from 1000 nonparametric bootstrapped samples (Efron & Tibshirani 1993). Resampling was conducted with replacement and stratified by year to characterize uncertainty in annual CPUE indices and by season and survey region for seasonal spatial relative abundance.

Coherence among and drivers of annual trends
Dynamic factor analysis (DFA) is a multivariate, dimension reduction technique designed to extract common trends (i.e. latent variables) from a suite of relatively short, nonstationary time series, evaluate coherence among the time series by illuminating relationships with these trends, and quantify the influ-ence of measured covariates on the collection of time series (Zuur et al. 2003). Because the DFA frame work is structured to incorporate time series information, survey catch data and hypothesized covariates must be summarized on some time interval (often annual) prior to analysis , Peterson et al. 2017). Thus, DFA was used in this investigation to synthesize the commonalities among the time series of annual CPUE predictions generated for the 9 EMTs from the delta-GAMs, and to evaluate the influence of biophysical forcing and exploitation covariates on these time series. The DFA model was given by: ( 3) where y t is the vector (9 × 1) of the standardized (zscored) CPUE predictions for the 9 EMTs in year t, the vector (r × 1) of r common trends (r < 9) modeled as stochastic random walks, the matrix (9 × r) of EMT-specific loadings, x t the vector (q × 1) of q covariates, D the matrix (9 × q) of covariate effects, and R and Q the variance-covariance matrices associated with the observation error vector t (9 × 1) and process error vector t (r × 1), respectively, both of which were assumed to follow a multivariate normal (MVN) distribution with mean zero. Q was set to the identity matrix so that the model parameterizations were identifiable (Zuur et al. 2003).
Twelve annualized covariates representing phy sical forcing, biological processes, and fishery re movals were evaluated within the DFA models. Specifically, the hypothesized physical drivers of the MAB ecological community included the principal components-based winter (December to March) North Atlantic Oscillation (NAO) index of atmospheric pressure differences between Iceland and the Azores (Hurrell et al. 2003), the Atlantic Multidecadal Oscillation which quantifies decadal-scale variability in sea surface temperature, and the Gulf Stream Index which measures the position of the north wall of this boundary current and has been shown to impact physical processes on the shelf (Joyce et al. 2019 (Reynolds et al. 2007). The primary production anomaly ratio measures annual fluctuations in primary productivity (NEFSC 2022), while the small:large copepod ratio characterizes the size struc- ture and quality of the zooplankton community (Perretti et al. 2017), and these 2 metrics were used to explore bottom-up biological controls on the MAB assemblages. Finally, total commercial landings served as a measure of fishing pressure and was used to evaluate top-down controls (NEFSC 2022). Given that each of these processes could impact faunal recruitment and availability to the survey, and that NEAMAP primarily samples the late juvenile and adult stages of species in the MAB, all physical and biological covariates were lagged 0, 1, and 2 yr, and exploitation was lagged 0 and 1 yr. These predictor variables were z-scored prior to inclusion in the DFA models.
In the DFA framework, the R matrix specifies the variance and covariance structure among time series (Zuur et al. 2003). The mean of the annual CV estimates generated from the bootstrapped delta-GAMs was calculated for each EMT and incorporated as the diagonal elements of the R matrix to propagate time series uncertainty (Peterson et al. 2021). Models including 1, 2, and 3 common trends were fitted and included either no covariate, a single covariate, or a lagged covariate. This approach yielded 114 candidate DFA models. Model selection occurred by evaluating AIC corrected for small sample size (AIC c ) to identify empirically supported forms, where those with ΔAIC c < 10 were considered further ). The FitRatio, which is the ratio of the sumof-squared residuals and sum-of-squared observations ( ) for the DFA model fits to the EMT time series, was evaluated for each aggregate group as well as across these EMTs (i.e. Mean Fit) for the remaining models and coupled with a visual evaluation of residuals to identify the 'final' model. Statistical significance of factor loadings and covariate coefficients was assessed for this model using 95% confidence intervals (CIs), and the fits of this DFA model to each EMT time series represented a suite of ecosystem indicators for the nearshore MAB (e.g. Nye et al. 2010).

Synthesis of seasonal spatial distributions
Principal component analysis (PCA) was applied to the predicted seasonal spatial distributions of CPUE generated for the 9 EMTs from the delta-GAM models to synthesize spatiotemporal patterns of these assemblages. PCA is an ordination technique routinely used to achieve dimension reduction and signal isolation for complex multivariate datasets. It relies on singular value decomposition to quantify principal components (PCs) as linear combinations of the original variables, such that a subset of these PCs capture an appreciable proportion of the variability in the dataset (Zuur et al. 2007). While the application of PCA to observed (i.e. 'raw') species abundance information is usually not recommended given the typical dominance of zero observations and violation of multivariate normality in these datasets (Clarke et al. 2014), model-derived, ecological community data are much less prone to these issues. As such, the use of PCA for exploring patterns in ecological metrics generated using various analytical approaches is appropriate and has proven valuable when assessing system status in the Northeast USA LME .
In this investigation, the input dataset for the PCA was a 1000 × 16 site-by-assemblage matrix, where the 16 columns represented a given EMT-season combination (spring data on rays and other sharks excluded, see Section 3.1. for details), and rows were predictions of CPUE at each of the 1000 locations along the coastal contour. The relationships among factor loadings (representing EMT-season combinations) and scores (representing locations along the contour) were visualized using a correlation biplot (Zuur et al. 2007). Uncertainty in the factor loading vectors for each EMT-season combination, which represented gradients in relative aggregate biomass for each assemblage, were characterized using a nonparametric approach that incorporated bootstrapped predictions from the GAMs (Mehlman et al. 1995). All statistical analyses were performed using the R software program (v3.6.3, R Core Team 2020). Packages 'gamlss' (Rigby & Stasinopoulos 2005), 'MARSS' (Holmes et al. 2018), and 'stats' were ac ces sed to fit the delta-GAMs, DFAs, and PCA, respectively.

Ecomorphotype catch summaries
NEAMAP completed 24 survey cruises and collected data at 3532 sites from 2008 to 2019. A total of 178 species was encountered, and the aggregate catch was comprised of 13 436 602 individuals weighing 851.7 t. The demersal fishes category included 82 species and was encountered at 99.5% of the sites sampled. This EMT encompassed the largest number of taxa and had the highest frequency of occurrence among all groups, which was not unexpected given the bottom-tending, mobile nature of the sampling gear. Although previous studies of marine communi-ˆ t 2 / y t ties based on bottom trawl survey data have excluded information on pelagic fishes due to catchability concerns (e.g. Buchheister et al. 2013), the high-rise nature of the NEAMAP trawl and sampling of relatively shallow depths combined to yield appreciable, consistent catches of this assemblage. A total of 37 pelagic fish species was collected across 96.4% of NEAMAP trawls. Although only 4 species of skate were sampled, this group was quite ubiquitous as the frequency of encounter was 98.8% and second only to demersal fishes. Flatfishes were collected at 86.1% of sites and included 11 species. The 5 taxa that composed the cephalo pod EMT were captured at a frequency (88.8%) similar to that of the flatfishes. Benthic arthro pods were a relatively diverse EMT comprised of 17 taxa and encountered in 73.9% of hauls, while dogfishes included only 2 species that were caught at 68.3% of sites. Encounter frequencies of rays (9 species) and other sharks (11 species) were relatively low, as these EMTs were collected in 28.2% and 14.8% of trawl hauls, respectively. Both of these EMTs were rarely encountered during NEAMAP spring surveys, and thus the early season catch data on these assemblages were excluded from subsequent modeling efforts. Rays were captured in 46.8% of tows that occurred in the fall, while other sharks were sampled at 23.8% of these sites, and the fall species compositions for each of these EMTs mirrored those observed at the annual scale.

Spatiotemporal modeling of ecomorphotype relative aggregate biomass
GAM models were successfully fitted for each of the 9 EMTs, as no convergence or estimation problems occurred for any of the candidate parameterizations. The generalized beta type-2 (GB2) was the most supported distribution for modeling non-zero CPUE data associated with demersal fishes, pelagic fishes, benthic arthropods, and skates, while the Box-Cox t (BCTo) distribution was chosen for the flatfish, cephalopod, dogfish, and ray EMTs (Table S2). The non-zero CPUE data for other sharks was assumed to follow the generalized inverse Gaussian (GIG) distribution. The most supported parameterizations for each assemblage were considered reliable for estimating the annual trends and seasonal spatial distributions of these EMTs based on a thorough evaluation of the associated model diagnostics, which in cluded worm plots (detrended quantile-quantile plots), histograms of residuals, and relationships be -tween residuals and fitted values as well as between residuals and the covariates.
Model M 1 , the full model (Table 1), was the most supported parameterization for both the binomial and conditional models for the flatfish, benthic arthropod, and dogfish assemblages, as well as for the binomial component of pelagic fishes and the conditional model for cephalopods (Table S3). This formulation also yielded the lowest AIC values for both the demersal fish and skate EMTs. In each of these cases, depth and solar zenith explained appreciable variation in the encounter or non-zero CPUE of these EMTs. The selected parameterization for the binomial component of the cephalopod model and conditional component of the pelagic fishes model included the depth covariate, but not solar zenith (i.e. M 2 ). Given that the datasets for the ray and other shark EMTs were restricted to the fall cruises, model M 6 received the most empirical support for both components of the ray model and for the binomial component for other sharks, while M 8 was chosen for the conditional model of the latter group.
CVs (not shown) associated with predicted annual CPUE estimates were less than 0.33, indicative of good precision, for all of the EMTs except pelagic fishes, rays, and other sharks. For these latter 3 groups, maximum an nual CVs were 0.40, 0.41, and 0.51, respectively. The CVs of the spatial predictions during spring ranged from 0.09 to 0.63 and were between 0.08 and 0.91 for most EMTs in fall. The CVs associated with the other sharks CPUE during fall were as high as 9.63 in the northern extent of the survey range, likely due to the sporadic catches of relatively few large animals that occur in this region during fall.

Coherence among and drivers of annual trends
The DFA model receiving the most empirical support and thus used for inference included 3 common trends and SSTWI (i.e. MAB mean sea surface temperature during winter; Fig. 2) as a co variate. Although a competing model with 3 common trends and NAO had a ΔAIC c value of 0.4 and a slightly lower Mean Fit (Table S4), FitRatios of DFA model fits for 7 out of 9 EMTs were larger than those from the selected model, and thus this model was not considered further. The first common trend indicated a notable shift in CPUE from slightly above average prior to 2014 to a period of lower values through 2017 (Fig. 3a). Relative ag gregate biomass then increased throughout the remainder of the time series. The factor loading of demersal fishes on this common trend was statistically significant and positive, while those of the cephalopod, ray, and other shark assemblages were significant and negative, revealing contrasting responses of relative aggregate biomass over the 2008−2019 period between the demersal fishes and these latter 3 groups (Fig. 3b). The SSTWI covariate had a significant, positive effect on the CPUE of the demersal fishes, and the influence on the ray group was also positive but not significant (ray 95% lower confidence limit [LCL]: −0.06; Table S5). The DFA model fit to demersal fishes showed a gradual de cline early in the time series, a peak in 2012 that was followed by a sharp decrease from 2013 to 2014, and then a gradual increase through 2019 with a secondary peak during 2016 (Fig. 3c). Fits to both the cephalopod and other shark EMTs exhibited increasing trends in relative abundance between 2008 and 2019 and, in contrast with demersal fishes, the CPUE of both assemblages The relative ag gregate biomass of these 2 groups was also high in 2017. The abundance of the ray group varied through out the time series with little discernable trend. These model fits represent a new suite of ecosystem indicators for these EMTs in the nearshore MAB ecosystem. Fits to the demersal fish, cephalopod, and other shark EMTs were good, as FitRatios ranged from 0.12 to 0.21, while that of the rays was relatively poor at 0.66 and likely responsible for the equivocal trend. The second common trend reflected belowaverage CPUE values during 2008 and 2009, a somewhat parabolic trend between 2009 and 2015 where relative aggregate biomass increased sharply and then declined gradually, and a subsequent stable period through 2019 (Fig. 4a). The 3 EMTs with significant factor loadings on this common trend again displayed conflicting patterns in their time series of CPUE. The factor loading of pelagic fishes was positive, and those of the benthic arthropods and dogfishes were negative (Fig. 4b). The DFA model fit for pelagic fishes exhibited a cuneate trend from 2009 to 2015, and the highest CPUE for this assemblage occurred during 2012 (Fig. 4c). Similar to the demersal fishes, pelagics were abundant during 2016 and this group was significantly and positively influenced by SSTWI. Unlike the demersal fishes, however, the abundance of this EMT mostly declined after 2012. Model fits indicated that the highest CPUE values for dogfishes occurred during the early and late periods of the time series, with below-average values during most of the intervening years. The relative aggregate biomass of benthic arthropods peaked in 2009, and in contrast with pelagic fishes, exhibited consistently low CPUE between 2010 and 2013. A secondary peak occurred in 2014 and was followed by a steady decline. The FitRatios of the ecosystem indicators for these 3 EMTs ranged from 0.07 to 0.41, indicative of good model fits.
The third common trend showed a slight increase at the beginning of the time series, followed by a consistent decline thereafter (Fig. 5a), and both the flatfish and skate EMTs were significantly and positively coherent with this trend (Fig. 5b). Although the DFA model fits for these groups displayed notable declines through time, the exact patterns exhibited by these assemblages varied (Fig. 5c) Standardized CPUE Fig. 4. As for Fig. 3, but for the second common trend covariate on this group, which was positive but not significant (skate 95% LCL: −0.01; Table S5). DFA model fits for both of these EMTs were good, as the FitRatios associated with the flatfishes and skates were 0.07 and 0.11, respectively.

Synthesis of seasonal spatial distributions
The PCA applied to the seasonal spatial predictions of CPUE for the 9 EMTs revealed appreciable differences in the locations of the nearshore MAB associated with the largest relative aggregate biomasses of these groups as well as seasonal contrasts within most assemblages. The first principal component (PC1) accounted for 37.6% of the variation in the seasonal spatial distributions of these EMTs, and effectively divided the northern and southern areas of this ecosystem at approximately Delaware Bay (Fig. 6). PC2, which delineated a gradient that separated the boundaries of this EPU from its central region, captured 28.7% of the variance in the seasonal distributions of the groups. During fall, the CPUE of demersal fishes, rays, and other sharks increased toward the southern zone of the nearshore MAB and generally peaked in the vicinity of Oregon Inlet, NC. Relative aggregate biomass of demersal fishes was highest near the northern and southern boundaries of this ecosystem in spring as indicated by the central position of, and relatively large uncertainty associated with, the factor loading of this group on PC1, along with the strong negative loading on PC2.
Spatial patterns of relative abundance for pelagic fishes during fall mirrored those of the demersal category in spring, and while spring catches of the pelagic assemblage were also elevated in the northern and southern extents of this EPU, relative abundance of this EMT tended to be greater in RIS (i.e. between Block Island, RI, and Martha's Vineyard, MA). The CPUE of benthic arthropods peaked during both seasons in the vicinity of Delaware Bay, and the dogfish group was also most abundant in this area during spring. Cephalopod relative aggregate biomass was greatest between central Long Island and about Delaware Bay in spring, while the abundance of dogfishes was highest in this area during fall, indicating that these EMTs exhibit alternating seasonal use of this area. The south shore of Long Island, and particularly the region between Moriches Inlet, NY and Montauk, represented a hotspot of elevated biomass for skates, flatfishes, and cephalo -

DISCUSSION
This investigation quantified and synthesized temporal and spatial patterns of the ecological community inhabiting the nearshore MAB, a system that serves as a vital spawning, nursery, foraging, and re fuge habitat for an array of economically and ecologically important living marine resources (Wuenschel et al. 2013, Malek et al. 2014). The annual time series and seasonal spatial distributions of relative aggregate biomass quantified for the 9 EMTs using delta-GAM modeling provided insights into spatiotemporal trends in the productivity and partitioning of bio mass among these groups (Buchheister et al. 2013). These patterns both illuminated the community structure of the EMTs and complimented prior ecological investigations that focused on the offshore and estuarine environments of this EPU. Identifying the commonalities among and drivers of the annual trends using DFA offered a key framework for summarizing the oftenoverwhelming number of time series available for consideration when evaluating ecosystem status. This approach effectively reconciled seemingly disparate temporal trends among multiple species assemblages while yielding a new suite of ecosystem indicators for this system. The synthesis of EMT seasonal spatial distributions using PCA represented a novel, albeit somewhat coarse, system-level approach to quantifying community spatial gradients that could be used by decision makers charged with mitigating risks to important marine habitats and balancing trade-offs among multiple use sectors in marine spatial planning activities.
Evaluating the spatiotemporal patterns of species assemblages in an ecosystem first requires the implementation of a classification scheme to group taxa based on similarities across traits of interest. The ap proach used to define these assemblages represents a vitally important decision, as it both influences the perceived structure of the ecosystem and shapes inference on system processes (Frost et al. 1992). To our knowledge, none of the myriad approaches commonly used to delineate assemblages has been identified as superior to the others  and thus, much like the EAFM process itself (Link 2010), selecting an appropriate grouping framework involves the consideration of trade-offs. However, because these assemblages often serve as the basis for developing ecosystem metrics, including time series indicators that underpin ecosystem assessments, it is critical that the approach used to define these groups maintains coherence with as much of the established guidance for indicator development as possible (e.g. Dale & Beyeler 2001, Rice & Rochet 2005. The common themes of this guidance are that the measures be: (T1) based in ecological the-   ory; (T2) relatively inexpensive to produce; (T3) generated consistently over time; (T4) responsive to abiotic, biotic, and anthropogenic drivers; and (T5) easily understood by the scientific community, fishery managers, and stakeholders, among others. All species collected by NEAMAP from 2008 to 2019 were assigned to 1 of 9 EMTs, which effectively mirrored assemblage definitions used previously to evaluate community status in the broader Northeast USA LME, although the groupings in these prior studies were not described as EMTs (Fogarty & Murawski 1998, Methratta & Link 2006. This EMT framework follows many of the themes associated with the development of informative and reliable ecosystem metrics. First, this classification system allowed for the incorporation of all taxa sampled by NEAMAP which, although the contributions of some species to the aggregate biomass of their assigned EMT was relatively small, was critical to support a holistic evaluation of the spatiotemporal patterns of the entire ecological community as represented by the survey trawl (T1 listed above; Meth ratta & Link 2006). Assigning taxa to the 9 EMTs was based on information available in the literature and thus was low-cost (T2), and unlike with other modes of classification, species EMT assignments likely are robust to systemic changes (e.g. demersal fishes remain demersal) which should promote longer-term consistency in the species composition of the resulting metrics (T3). These EMTs also generally align with the fisheries operating in the MAB ecosystem (Gaichas et al. 2014), and therefore are likely to be responsive to fisheries management measures (T4) and intuitive to fishery managers, stakeholders, and the public (T5). Note that further subdividing these 9 EMTs to gain greater alignment of assemblages with MAB fisheries (e.g. divide demersal fishes into warm-and cold-adapted species), perhaps by applying DFA at the specieslevel to identify these groupings, represents a potentially fruitful area of future research. Finally, and perhaps most importantly, this relatively simple EMT classification approach can be implemented for any ecosystem with at least a modest data collection program. As the value of EAFM gains greater appreciation and efforts to quantify ecosystem dynamics become more widespread, developing standard metrics (and first, assemblage delineations) based on relatively simple criteria and minimal data requirements will help to facilitate cross-system comparisons meant to illuminate common themes in ecosystem structure and function (Link & Marshak 2019).
As noted above, the choice of an assemblage framework to serve as the basis for quantifying ecological communities and developing ecosystem indicators involves trade-offs. Given similarities in habitat usage within each assemblage, spatiotemporal abundance patterns of EMTs are expected to be responsive to environmental cues and harvest pressures. Note, however, that these groupings likely would mask any trends in the ecological community driven primarily by changing trophic interactions, since most of the EMTs contain a myriad of both predators and prey as well as potential competitors (Garrison & Link 2000, Link & Auster 2013. The classification framework currently underpinning aggregate biomass indicator development for the offshore waters of the MAB shelf relies on the trophic guild concept as the basis for assemblage definition (Garrison & Link 2000, NEFSC 2022, which offers the ad vantage of being more likely to detect trends driven by changes in these species interactions. Drawbacks of the trophic guild approach include increased costs to generate the diet data required to delineate these groups and, at times, misalignment of these assemblages with established fisheries and associated management. The question as to which of these as semblage frameworks is preferable for quantifying ecological communities in support of EAFM is not new . Because these 2 ap proaches focus on different aspects of an ecosystem, we contend that both should be considered in the development and maintenance of EAFM, with re search directed toward synthesis of the metrics generated from each. The annual trends and seasonal spatial distributions of relative aggregate biomass were quantified for the 9 EMTs using delta-GAMs. Given the prevalence of nonlinear relationships between community measures and hypothesized synoptic localized drivers, recognition of the value of GAMs in standardizing and quantifying ecosystem metrics is increasing (Buchheister et al. 2013, Large et al. 2015. The availability of a broad array of candidate distributions in GAMLSS represented a key advancement over frame works that limit choices to the exponential family of distributions (Rigby & Stasinopoulos 2005), as the ability to model CPUE data using the GB2, BCTo, and GIG distributions yielded marked improvements in model diagnostics, and presumably greater accuracy and precision in predicted measures of central tendency. This framework also permits the modeling of scale (variance) and shape (skewness and kurtosis) as a function of covariates which, although not necessary in this study based on model diagnostics, could prove useful in other applications and thus represents a valuable area of future research in the development of ecosystem indicators.
Evaluating trends in and patterns of relative aggregate biomass for species assemblages has received criticism following the revelation that these metrics are prone to bias and can be hypersensitive when catchability of the component species is not constant (Kleiber & Maunder 2008). The catch data underpinning this study were derived from a fisheries-independent monitoring program that has maintained consistency in both the fishing system and sampling protocols since its inception, such that the assumption of constant species-specific catchability over space and time was reasonable, and thus any bias or hypersensitivity in these metrics was likely minimal. Indeed, preliminary analyses showed that the spatiotemporal patterns in the aggregate biomass indicators were consistent with those of many component species within these EMTs, and the value of such metrics for quantifying ecosystem structure and function has been well documented , Methratta & Link 2006, Buchheister et al. 2013, Latour & Gartland 2020, meaning that the aforementioned concerns were outweighed by the need to gain insight into the spatiotemporal patterns of the community assemblage in this system. Likewise, the restricted sampling frame and short time series of available data represented potential constraints on the applicability of this investigation. While generating information on ecosystem status at the LME (or at least EPU) level had been assumed to be an appropriate scale to support IEAs, delineating system metrics such as assemblage indicators at multiple spatial scales congruent with the hierarchal structure of the ecosystem (e.g. estuary vs. coastal ocean vs. offshore shelf) may be more relevant to these IEA efforts (Heim et al. 2021). Following guidance to avoid inference based on statistically significant annual trends in time series of less than 30 yr (Hardison et al. 2019), this study instead employed an approach designed for use with short, nonstationary time series typically associated with relatively new fisheries-independent monitoring programs.
The time series of relative aggregate biomass for the EMTs inhabiting the nearshore MAB were characterized by 3 common trends, which effectively represented latent, unmeasured variables that influenced the relative abundances of these assemblages (Zuur et al. 2003). The first common trend was characterized by an abrupt shift from positive to negative values during the 2013 to 2014 timeframe, which likely reflected a severe, system-level disturbance that impacted nearly half of the assemblages. Because SSTWI was included in the model that received the most empirical support, the effect of this co variate was removed from all of the common trends, and as such the source of this perturbation remains unresolved. An ideal DFA model would include multiple identifiable covariates that explain nearly all of the variation in the time series of interest, yielding a single common trend with no prominent features (Zuur et al. 2003). Attempts to include more than one covariate in a DFA model resulted in convergence issues that likely were due to the relatively short time series of available data, and thus represented a limitation of this analysis. The second common trend was approximately parabolic during the first half of the time series, while the third reflected a steady decline. Each of these latent variables was associated with multiple EMTs, meaning that the annual indices of these 9 assemblages primarily reflected 3 relatively simple trends.
To be clear, the DFA model fits to the EMT time series, and not these common trends, constituted the new suite of ecosystem indicators generated for the nearshore MAB from this investigation. Nevertheless, DFA common trends could prove valuable to assessments of ecosystem status in support of EAFM, as they group ecosystem indicators by identifying coherence (or conflict) among these metrics and provide synthetic information on the underlying patterns in the groups of indicators. Efforts to assess system status and implement EAFM are becoming overwhelmed by the myriad of indicators that must be considered (Bundy et al. 2019). The ability to synthesize the information content of multiple indicators based on a smaller number of latent measures (common trends) therefore represents a critical step towards promoting efficient ecosystem assessment and management activities, and one that likely will become more important as the number of indicators available for consideration continues to increase. Further, although this investigation was unable to identify the 3 latent covariates driving the annual trends in EMT relative aggregate biomass, information provided by the common trends could be used in a qualitative manner during ecosystem assessments to identify events such as perturbations (e.g. first common trend) or other concerning patterns (e.g. third common trend) that should be prioritized for further exploration. Process-oriented studies designed to unmask the covariates underlying these common trends represent a very important future direction for investigations focused on understanding community dynamics in the MAB.
The factor loadings of the EMTs associated with these 3 common trends provided insight on the responses of assemblage time series to system-level forcing (Zuur et al. 2003. Notably, the 7 groups associated with the first and second common trends showed conflicting patterns in relative abundance during the study period. Trends in relative aggregate biomass for the cephalopod, ray, and other shark EMTs generally opposed that of the demersal fish category, while the pelagic fish trend contrasted with those of the benthic arthropod and dogfish groups. A central tenet of ecosystem-based management is that trade-offs must be evaluated to sustain an optimized delivery of desirable market goods and nonmarket services from these systems (Link 2010). Quantifying these conflicting responses among EMTs in the nearshore MAB illuminated some of the potential trade-offs, and as such highlighted the value of DFA as a synthetic tool for suites of ecosystem indicators, even when these indicator time series are relatively short. The annual patterns of the flatfish and skate assemblages were not in conflict, but instead both mirrored the third common trend, and the general decline in the relative abundance of these 2 EMTs may have reflected actual trends in the abundances of the component species (Perretti & Thorson 2019, NEFSC 2020, changing availability to the NEAMAP survey as these taxa respond to climate change (Nye et al. 2009, Morley et al. 2018, or a combination of these hypotheses. The statistically significant influence of SSTWI on the demersal and pelagic fish EMTs, and near significance for skates and rays, was somewhat unexpected given that NEAMAP samples these assemblages during spring and fall. While it is plausible that warmer winter temperatures trigger earlier migrations into the nearshore MAB that in turn yield increased catches during spring, trends in aggregate biomass were generated annually for most EMTs and were limited to the fall for rays, indicating that winter temperatures may drive observed community patterns at an annual scale in this coastal zone. As noted previously, the DFA model fits to the time series of EMT relative aggregate biomass represent a new suite of ecosystem indicators for the nearshore MAB meant to support the continued development and implementation of EAFM for this EPU. An array of indicators has been quantified for this ecosystem and are updated annually in Mid-Atlantic State of the Ecosystem (SOE) reports (NEFSC 2022). Following the IEA framework for EAFM, the indicators in this report are used to inform risk assessments, where the highest risk elements are explored further through the development of conceptual models, and the performance of measures designed to mitigate these risks is tested via management strategy evalu-ation (Gaichas et al. 2016, Muffley et al. 2021). While the time series of many of the current ecosystem indicators span decades, several measures are reported on time scales that are equal to or shorter than those of the EMT indicators developed in this study, and have been included in risk assessments. As such, it is anticipated that these EMT ecosystem indicators could be incorporated into the SOE reports, where they would then be used during the annual risk assessment process specifically to evaluate ecosystem status relative to the objective of maintaining system biomass (NEFSC 2022). A 'high risk' designation could result in the development of conceptual models and testing of management strategies to minimize the identified risk. Continued monitoring of these indicators would provide feedback on the efficacy of implemented management measures, and thus aid in advancing EAFM in the MAB. Given that risk assessments are currently based on somewhat qualitative evaluations of indicator trends, developing target (i.e. desired state) and threshold (i.e. state to be avoided) reference points for these metrics, including those generated from this investigation, is an important area of future research (Large et al. 2015).
Much of the research associated with the development of ecosystem metrics for species assemblages has focused on evaluating inter-annual trends in these measures, whereas quantifying the seasonal patterns and spatial variability of these groups has received less attention (Buchheister et al. 2013). A clear north-to-south gradient was evident in the seasonal community structure of the nearshore MAB, with the nuance that assemblage compositions at the ecosystem boundaries were generally more similar relative to the central regions. Interestingly, the relative aggregate biomass of most EMTs appeared to cluster near the entrances to estuaries (i.e. Oregon Inlet for Pamlico Sound, NC; the mouth of Delaware Bay, and Moriches Inlet associated with the Moriches Bay estuary). A host of species inhabiting the MAB utilize estuaries during the warmer months (Murdy et al. 1997, Able & Fahay 2010, such that these results may reflect the composite distributions of EMT component species as they undertook seasonal migrations to (spring) or from (fall) these systems. Alternatively, estuarine outflows may create areas of elevated habitat suitability in the coastal ocean, and the implementation of multispecies ecological niche modeling to explore the distributions uncovered here from a mechanistic perspective represents a promising area of future research (Roberts et al. 2022). It is worth noting, however, that none of the factor loadings of the EMT seasonal-spatial distributions were associated with Chesapeake Bay, one of the largest estuaries in the USA. These findings are consistent with documented, appreciable declines in the relative usage of this system by several seasonally resident taxa in recent years (Schonfeld et al. 2022).
The synthesis of the seasonal spatial distributions of the EMTs was based on fisheries-independent bottom trawl data, which necessitated a somewhat coarse treatment of these distributions given that each trawl haul likely sampled over a range of unique bottom types and other environmental conditions (Sullivan et al. 2006). Nevertheless, this information on the general distributions of these assemblages could be used to inform the development of more-refined future investigations of multispecies or assemblage essential fish habitat either by directing the focus of these studies toward areas associated with large relative abundances of several EMTs or those where alternating seasonal use of a region oc curs. Further, given the expected continuing impacts of climate change and increased uses of the MAB ecosystem for non-fishing human activities (Meth ratta 2020), these seasonal spatial distributions could be valuable to future climate vulnerability assess ments and in the context of marine spatial planning activities, both of which contribute to advancing EAFM (Malek et al. 2014, Farr et al. 2021. For example, the region along the south shore of Long Island, NY, and specifically along the eastern half of the island, appeared to represent a distinct subunit of the nearshore MAB ecosystem that supported elevated biomass of 3 EMTs across both spring and fall. This aggregation of multiple assemblages, coupled with the declining trends in the ecosystem indicator time series of 2 of these 3 EMTs, suggests that this zone could be especially sensitive to the impacts of climate change, and should receive particular consideration by decision makers tasked with managing current and planned use sectors in this area. Fishery managers in the MAB have expressed the desire to move toward EAFM at a rate that is commensurate with the available science (Muffley et al. 2021). Overall, this investigation provided a classification scheme for generating ecological metrics of species assemblages that required minimal auxiliary information, built on previous efforts, and was designed to promote consistency and comparisons among ecosystems. Synthesizing the spatiotemporal patterns in the relative aggregate biomass of the resulting assemblages yielded, to our knowledge, the first integrated evaluation of the structure of the nearshore MAB ecological community, and thus much-needed baseline information against which to evaluate future system states. Further, when coupled with comparable ecosystem information derived from the offshore waters and estuaries of this EPU, these results generate a holistic characterization of this ecosystem that likely will prove instrumental to advancing EAFM in the MAB.