Comparative Abundance, Species Composition, and Demographics of Continental Shelf Fish Assemblages throughout the Gulf of Mexico

We analyzed the results of the first comprehensive, systematic, fishery-independent survey of Gulf of Mexico (GoM) continental shelves using data collected from demersal longline sampling off the United States, Mexico, and Cuba. In total, 166 species were sampled from 343 longline sets during 2011–2017, which deployed 153,146 baited hooks, catching 14,938 fish. Abundance, species richness, and Shannon–Wiener diversity indices by station were highest in mid-shelf depths (~100 m), declining by about half in deeper waters. Six spatial assemblages were identified by testing the results of cluster analysis using similarity profile analysis and then plotting the geographic location of identified station clusters. A high degree of depth-related and horizontal zonation was evident for demersal fish species. Multispecies CPUE (number per 1,000 hookhours) was highest off the north-central (NC) and northwestern (NW) GoM and lower on the West Florida Shelf (WFS), Cuba (CUB), Yucatan Peninsula (YP), and southwestern (SW) GoM. Snappers and groupers were most abundant in the WFS and CUB, while elasmobranchs were the dominant taxa in the NC and NWGoM. Pelagic species were relatively rare everywhere (owing to the use of demersal longline gear), but were most dense off CUB. Species richness was highest in the NC and WFS subareas and lowest in the NW and CUB. Slopes of multispecies size spectra, which integrated mortality, recruitment, growth, and species interactions among size-groups, were shallowest in the NW and NC GoM and steepest off the WFS and YP. These results provide a basis for evaluating the relative resiliency potential of species assemblages across the continental shelves of the GoM, and thus for identifying subareas that are most vulnerable to acute and chronic perturbations from cumulative effects of fishing, climate change, pollution (including oil spills), habitat loss, and invasive species. Subject editor: Patrick J. Sullivan, Cornell University, Ithaca, New York *Corresponding author: smurawski@usf.edu Received December 26, 2017; accepted April 9, 2018 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Marine and Coastal Fisheries: Dynamics, Management, and Ecosystem Science 10:325–346, 2018

Productivity, species diversity, and demographics of fishes occupying large marine ecosystems (LMEs), such as the Gulf of Mexico (GoM; Figure 1), are affected by a wide range of both chronic ("press") and acute ("pulse") perturbations (May 1974). These perturbations include fishing (commercial, recreational, and subsistence), climate change and variability, habitat loss, pollution (both chronic [e.g., nutrient inputs] and acute [e.g., oil spills]), invasive species, disruptive natural phenomena such as hurricanes, and other factors. Fishery and environmental management regimes, which may vary widely among nations sharing contiguous LME waters (Tunnell 2017), will also influence the abundance and demographics of species both within territorial jurisdictions and across national boundaries, especially for species with high dispersal potential. Effects of intensive perturbations can be manifested in reduced abundance of species and whole assemblages, truncated size distributions, and higher variability in recruitment, all of which undermine ecosystem resilience. Responding to the diverse challenges presented by chronic and acute pressures on ecosystems requires information on fish assemblages within LMEs and determination of their relative resilience potential to existing threats (Fogarty and Murawski 1998;Dulvy et al. 2008;Stefansdottir et al. 2010;Johannesen et al. 2012). Management actions, such as controlling excessive fishing mortality and protecting vulnerable habitats and rare species, can mitigate some threats to ecosystem resilience, thereby reducing the risk of catastrophic regime shifts (deYoung et al. 2008; Barange et al. 2010).
Defining and evaluating fish assemblages require systematic data collected using consistent methods that sample multiple species simultaneously, employing an appropriate statistical design that is amenable to robust comparisons over space and time. Assemblage structure and species abundances have been evaluated for GoM fishes, but previous studies have been restricted to spatial subsets or have been used to evaluate regionally specific environmental management issues. For example, assemblages of fishes have been determined for nearshore waters of the northern GoM (Lewis et al. 2007;Jordan et al. 2010), on the northwestern (Monk et al. 2015) and West Florida continental shelves (Darcy and Gutherz 1984), and the northwest continental slope (Wei et al. 2012) by using bottom trawling. Ichthyoplankton surveys have been used to define and analyze assemblages of larval fishes independently in the northern (Muhling et al. 2012) and southern GoM (Flores-Coto et al. 2014). Using fisherydependent data, species assemblages were evaluated for U.S. waters (Scott-Denton et al. 2011;Farmer et al. 2016;Pulver et al. 2016). Assemblages of fishes have also been evaluated on artificial versus natural reef habitats (Streich et al. 2017) to assess the effects of hypoxic waters (Craig and Bosman 2013), to determine impacts from the 2010 Deepwater Horizon (DWH) oil spill in nearshore waters (Schaefer et al. 2016), and for impact assessment from the 1979-1980 Ixtoc-I oil spill off Mexico (Soto et al. 2014;Amezcua-Linares et al. 2015). In a limited number of cases, multinational surveys in the GoM have been conducted, but these either sampled spatial subsets of the shelf (e.g., Nelson and Carpenter 1968) or were directed at particular species groups (e.g., coastal sharks; NOAA Fisheries 1998). Thus, while a number of fishery-independent surveys and results from commercial fisheries data have been used for assessment of fish assemblage characterization prior to our study, there have been no comprehensive, systematic, and multinational (e.g., USA, Mexico, and Cuba) surveys of a broad spectrum of GoM fish resources, thus obviating GoM-wide fish assemblage evaluations.
In this study, we analyze the results of the first comprehensive, systematic, fishery-independent survey of GoM continental shelves ( Figure 1). This study was motivated by the DWH oil spill (Lubchenco et al. 2012) and its likely pulse impacts on the resiliency of continental shelf fish populations in the GoM (DWH Natural Resource Damage Assessment Trustees 2016; Murawski et al. 2016). Given the ever-expanding quest for oil and gas resources in all three national exclusive economic zones of the GoM, it is prudent to assess the vulnerability of fish assemblages at the LME scale. An important consideration (and hypothesis) is the extent to which fish communities are structured in multiple dimensions (i.e., depth, horizontal, time). If there is significant spatial modularity, then identifying the relative sensitivity of spatially structured fish assemblages to stressors becomes an important factor in designing appropriate risk-averse conservation strategies. Using data we collected from demersal longline sampling off the United States, Mexico, and Cuba during 2011-2017, spatial assemblages of fishes were identified, the abundance and diversity of fishes in each assemblage were compared, and differences in community-level demography metrics (multispecies size and diversity spectra) were evaluated. These results thus provide the basis for evaluating the relative resiliency potential of species assemblages across the continental shelves of the GoM and for identifying particularly vulnerable subareas. Implications for resource management strategies that may strengthen resilience to perturbations are discussed.

STUDY AREA
The GoM is a large (158 × 10 6 ha), deep (average depth = 1,485 m), semi-enclosed sea (Figure 1; Darnell 2015) that supports highly diverse fish communities (1,500+ species; McEachran 2009) and valuable fisheries (Murawski et al. 2016). Continental shelf waters of the 326 GoM exhibit a diversity of benthic habitat types and depth profiles (Figure 1). The West Florida Shelf (WFS) and Yucatan Peninsula (YP) shelf are broad, gently sloping, sandy areas with carbonate underlayment interspersed with carbonate rock outcroppings and, in more southern areas, coral reefs (Uchupi and Emery 1968;Tunnell et al. 2007;Darnell 2015). At the western edge of the YP, the bathymetry is extremely steep (Figure 1). The north-central (NC), northwestern (NW), and southwestern (SW) GoM shelves mainly comprise fine-grained mud and clay sediments of terrigenous origin, interspersed with hardbottom banks (NW GoM; Dennis and Bright 1988) and shallow coral reefs (SW GoM; Tunnell et al. 2007). The shelf is extremely narrow in the SW GoM, whereas it broadens in the NW and NC areas. The area off northwestern Cuba (CUB) is steeply sloping and consists of both shallow and mesophotic coral reefs interspersed with sand and mud (Claro et al. 2001). Substantial "built" fish habitat, consisting mainly of oil and gas infrastructure (wells and pipelines; Murawski and Hogarth 2013) and other artificial habitats, exist in all subareas of the GoM except off northwest CUB.
In the wake of the DWH oil spill, we undertook GoMwide sampling of demersal fish communities off the USA, Mexico, and CUB primarily to monitor aspects of fish health (Murawski et al. 2014), to identify potential FIGURE 1. Geographic locations and relative abundance (number per 1,000 hook-hours) of continental shelf fishes sampled in the Gulf of Mexico (GoM), 2011-2017 (WFS = West Florida Shelf; NC = north-central GoM; NW = northwest GoM; SW = southwest GoM; YP = Yucatan Peninsula; CUB = Cuba). Fishes were sampled with demersal longline gear along defined shallow-to-deep transects. Locations of the Deepwater Horizon (DWH) and Ixtoc-I oil spills are indicated by black triangles. Solid black lines are national exclusive economic zone (EEZ) boundaries; dashed black lines are the putative boundaries of the GoM (Felder et al. 2009). GULF OF MEXICO FISH ASSEMBLAGES changes in populations and vital rates (Herdter et al. 2017), to obtain tissues for contaminant analyses (Murawski et al. 2014;Snyder et al. 2015;Granneman et al. 2017), and to complete associated baseline studies.

METHODS
Field sampling procedures.-Fishes were collected using demersal longline sampling gear and a transect survey design extending throughout the GoM continental shelves (Figure 1). Nominal station placement was in continental shelf waters from 40 to 300 m deep. We used longline sampling (baited hooks) as the standardized sampling methodology primarily because this gear is selective for large juvenile and adult fishes of exploitable sizes (Scott-Denton et al. 2011) that occupy relatively high trophic levels, and because the gear can be deployed in complex bottom habitats where trawls and other bottom-tending mobile gears are not efficient and may either destroy or be destroyed by the habitat. Because of their economic importance and central roles in food webs, "reef" fishes (e.g., snappers, groupers, and associated species) were our primary sampling targets, although our methods caught many nontarget demersal teleosts as well as elasmobranchs and some pelagic species (Table 1).
Longline sets were generally deployed at six stations along predefined transects that extended from relatively shallow to deep continental shelf areas ( Figure 1). The nominal depths sampled along each transect were 37, 73, 110, 146, 183, and 274 m. Along several transects, the bathymetric slope was so steep that six unique stations could not be effectively sampled, and depth control of those stations that were successfully sampled was difficult since the beginning and end of sets were 8 km apart (especially off northwestern CUB and the western YP; Figure 1). Thus, for some stations, the average depth exceeded the nominal set depth, accounting for a few outliers in the depth distributions of species ( Figure 2).
We used three chartered commercial fishing vessels (2011-2012) and the RV Weatherbird II (2012-2017) to deploy longlines. At each predetermined sampling location, we searched for suitable habitat for our target species groups. This search process involved using the ship's echosounder to locate "hard-bottom" habitat typically associated with concentrations of reef fishes. The vessel ranged up to 9 km from the center line of the sampling transect in search of suitable habitat. At each station, 8 km of 3.2-mm galvanized steel (2011)(2012) or 544-kg-test monofilament (2013-2017) main line were deployed, with a mean of 446 baited hooks per set. We used 136-kg-test leaders (2.4 m long) clipped to the main line and attached to size-13/0 circle hooks. Bait was cut fish (Atlantic Mackerel Scomber scombrus) and squid (primarily Humboldt squid Dosidicus gigas wings), switching baits haphazardly from hook to hook during deployments. The number of hooks in each set was counted as they were deployed. At the beginning and end of each set, we deployed a "high-flyer" buoy and attached Star:Oddi CDST Centi-TD temperature/time/depth (TTD) recorders to the main line with sufficient scope to reach the bottom and to record bottom time, bottom temperature (°C), and fished depth (m). The recording interval of these instruments was 5 min. At set-out and haul-back, we recorded latitude and longitude, time, depth, and the unique code numbers of the TTD instruments deployed at either end of the main line. Once the longline was deployed, the vessel returned to the start high-flyer buoy, and haul-back began. Fishing was accomplished only during daylight hours.
At retrieval, we determined species and recorded the SL, FL, and TL (cm), as appropriate for each specimen caught. Each specimen was weighed to the nearest gram on a Marel motion-compensated scale; large fish (>6 kg) were weighed with a hand scale (nearest 0.1 kg). For large sharks (e.g., ≥2 m), species were identified at the rail (not put aboard the vessel) and, to the extent feasible, the TL was estimated to the nearest 0.3 m prior to releasing the fish alive. Some captured specimens were incomplete (e.g., partially consumed by predators) or otherwise could not be accurately measured or weighed or identified to species. A total of 134 specimens (0.9% of the catch) was only identified to genus or were otherwise unknown (one specimen). In the cases where identifications were to genus, we included them in size-based analyses for species groups (e.g., snappers/groupers, etc.; see below).
Data standardization.-Abundance data obtained from each longline set were standardized to account for variations in the number of hooks deployed and the total soak time of each set of the gear. The standardization procedure adjusted the nominal catches to CPUE, defined as the number of fish caught per 1,000 hook-hours fished: CPUE i,j = individuals caught i,j × [(1,000/number of hooks j )/average hours of soak time j ] for species i and set j. Average soak time was calculated as [(B e − B s ) + (E e − E s )]/2, where B s , B e , E s , and E e are the times that the beginning (B) and end (E) "high flyers" were set (s) and retrieved (e). Average soak time of the gear was 2.08 h. The average station standardization coefficient (accounting both for the numbers of hooks fished and the set duration) was 1.24; thus, the adjusted catches were similar, on average, to the nominal catches obtained at each station. For the field work in CUB, the number of hooks deployed per set was halved from normal procedures to minimize gear loss and damage owing to the extremely rough bottom encountered there, which also shortened average soak times. On average, the calibration coefficients for sets in U.S. and Mexican waters were 1.05 times the average number of individuals caught, whereas the average calibration for CUB sets was 3.18 times the number of individuals caught.    The distribution of selected species in relation to average bottom temperature and average water depth (recorded at the beginning and end of sets) of capture was evaluated ( Figure 2). We also determined station-by-station abundance (CPUE), species richness (R; number of species caught), and the Shannon-Wiener diversity index H′ s ¼ À∑ R i¼1 ðp i;s log e p i;s Þ, where H′ s is the diversity index for station s; and p i,s is the proportion of species i at station s (Figure 3).
Each specimen was classified into one of four nominal species groupings (i.e., snappers/groupers, elasmobranchs, other demersal species, and pelagic species; Table 1) based on taxonomy and life history (Page et al. 2013). The use of species groups versus individual species data in some analyses was intended to (1) increase sample size and reduce station-to-station variability within and among the identified subareas and (2) reduce the impacts of relatively rare species for comparing subareas (Table 1). These species groups likely share a high degree of trophic overlap among areas.
Spatial classification.-Spatial patterns of species distributions were assessed using Bray-Curtis similarity for all pairwise station comparisons; the similarities were clustered using the group-average method, and the statistical significance of the resulting clusters was tested using similarity profile analysis (SIMPROF) in PRIMER software (Clarke et al. 2008;Clarke and Gorley 2015). The null hypothesis in this analysis was that the data sets are unstructured and the similarity profiles are random (Clarke et al. 2008). Individual longline sets (stations) were classified into site groups with similar species compositions based on cluster analysis using Bray-Curtis similarity (Figures 4, 5;Clarke et al. 2008;Clarke and Gorley 2015). Standardized catches (CPUE) for each species at each station were used as the basis for simultaneous classification of station groups and species associations (Figure 4; Clarke et al. 2008;Clarke and Gorley 2015). Square-root transformation of the data was used to minimize the effect of aberrant large catches. Similarity profiles were used to assess the statistical significance of clusters of stations (Figures 4, 5).
Abundance, multispecies size spectra, and diversity spectra.-Average abundances (CPUEs) by species group and subarea were calculated to evaluate differences in productivity and species dominance by subarea and to index the relative regional effects of differences in productivity, fishing, and other cumulative stressors ( Figure 6A). Statistical comparisons of average longline CPUE (log e [{multispecies catch per 1,000 hook-hours} + 0.1]) by GoM subarea were conducted with ANOVA and associated multiple range tests of subarea means (Bonferroni t-test).

334
For size spectra, we plotted log e standardized CPUE versus length (cm), since this depiction produced consistent and readily comparable indices (Figures 6B, 7; Tables 2, 3). Multispecies size, R, and H′ at length (Figure 6B-D; Tables 2, 3) were computed for each 3-cm length interval for each of the six identified subareas (Figure 7) and for the combined northern GoM (NGoM), southern GoM (SGoM), and Cuban "super-areas" (Figure 6). These super-areas reflect areas where, presumably, there are consistent policies of resource management affecting the population abundance and demography of managed resources. To combine data across taxa, we developed a length rule that used the most appropriate length measurements for different species and groups. All sharks were measured to TL, eels and groupers were measured to TL (the latter due to the shapes of their tails), and most other teleosts were measured to FL. Initial testing with 1-cm intervals indicated the same general patterns of size spectrum slopes, but 3-cm length groups showed lower variability and consequently better regression fits to the data. Similarly, R and H′ showed lower variability using 3-cm groupings. Linear regression models of size (length L) and diversity spectra were fitted to the descending limbs of abundance and H′ at length data, censored for incomplete gear selection; for regression analyses, we used 48 cm ≤ L ≤ 125 cm for size spectra (Figures 6B, 7; Tables 2, 3) and 48 cm ≤ L ≤ 150 cm for diversity spectra ( Figure 6D) to minimize the effects of small numbers of animals caught at larger sizes and to avoid the bias of incomplete gear selectivity of fish less than approximately 50 cm. We used ANCOVA to test the slopes, and, where appropriate, the adjusted means (intercepts) of regressions with common slopes for size spectra (Table 3; Figure 7).
Because of the disparity in total numbers of stations sampled per subarea, we computed species accumulation curves for each subarea (Figure 8; Gotelli and Cowell 2001;Cowell et al. 2004) to evaluate whether variations in sampling intensity confounded the comparisons of R and H′ among areas in the GoM. Two methods were used to compute species accumulation. First, deterministic species accumulation curves were computed without randomizing the selection of field samples. Alternatively, we computed stochastic accumulation curves (Figure 8, bottom panel) by using the "specaccum" function in the VEGAN package implemented in R (Oksanen et al. 2017). This procedure used 1,000 replicate sampling experiments employing the station-based field data in a jackknife procedure of station selection without replacement. The stochastic approach calculates the average distributions of the number of species encountered in 1, 2, 3,…, N sequential stations ( Figure 8, bottom panel).

RESULTS
We deployed a total of 153,146 baited hooks at 343 station locations and caught 14,938 specimens (Table 1) for a mean (unstandardized) catch of 44 fish/station and an average success rate of 10% (percentage of hooks with a retained fish). Depths and bottom temperatures associated with the captures of some of the most prevalent species showed consistent structuring along these environmental gradients (Figure 2). Given the negative relationship between depth (m) and bottom temperature (bt;°C: bt = 25.31 × exp[−0.0023 × depth], F = 846.7, P < 0.001, r 2 = 0.71), it is difficult to deconvolve which environmental variable is more important in determining this structure. However, in general, variability in temperature of occurrence was less than the variability in depth of occurrence, especially for deeper-dwelling species (Figure 2).

336
Given the wide range of depth and temperature conditions sampled (first row in Figure 2), the obvious species structuring along these gradients is not an artifact of the choice of station locations (Figure 1). Relatively shallow, warmwater species (e.g., Red Grouper, Atlantic Sharpnose Shark, and Red Snapper) were most abundant where the median depths of occurrence were less than the median depths sampled and where median bottom temperatures of occurrence were greater than the median sampled (first row in Figure 2). Conversely, the deep, coldwater species (e.g., Gulf Hake, Shortspine Dogfish, and Tilefish) were found at median depths of occurrence that were greater than the median depth for all stations and at median temperatures of occurrence that were less than the median for all stations. Species abundance, R, and H′ showed considerable variation, with peaks (especially in R and H′) around 100 m and lower levels at greater depths (Figure 3), especially beyond the photic zone (i.e., ≥200 m). In general, trends in CPUE, R, and H′ by depth were similar for the NGoM, SGoM, and CUB regions ( Figure 3); H′ at depth for the CUB region was somewhat greater than values for the other two regions.
Results of SIMPROF identified statistically significant groups of stations with similar species compositions (columns in Figure 4), and a simultaneous cluster analysis (without SIMPROF) identified species associations (rows in Figure 4). Although some significant station clusters were defined by unique occurrences of rare species (small numbers of stations per cluster), 16 of the clusters accounted for 91% of all stations (Figure 4). With some  Figure 4; not all station groupings are plotted (where small numbers of stations form unique groups).

GULF OF MEXICO FISH ASSEMBLAGES
exceptions, mapping of the station clusters revealed a series of semi-discrete subareas along GoM continental shelves ( Figure 5). Within these subareas, there were depth-related clines in species dominance, resulting in different shallow-water and deepwater species assemblages (Figures 2, 5). There were thus strong associations between identified species clusters ( Figure 4) and their respective depths and temperatures of occurrence (Figure 2). Some species ranged across multiple station clusters (e.g., Red Snapper, King Snake Eel, and Tilefish; Figure 4), while others were confined mostly to one cluster (e.g., Red Grouper, Little Gulper Shark, and Red Hind). Combining the station clusters into defined subregions was an important step and was based on the relative consistency of stations within the subregions versus variation across subarea boundaries ( Figure 5). In some cases, identified station clusters ranged across subarea boundaries (e.g., station group "i" in Figure 5, which occurred both in the NW and SW GoM). However, in splitting the SW and NW regions, other station groups were not common to both (e.g., station group "aa"). Some previous studies clustering fish species and sampling stations into regions have either prespecified station membership to sampling   Figure 5). Data are standardized catch rates (log e CPUE, where CPUE = number per 1,000 hook-hours) for 3-cm fish length intervals for all species captured. Linear regression lines are plotted for the descending limb (48-125 cm) for each subarea. Slopes and statistical significance of the resulting regression lines are provided in Table 2. 338 regions and/or restricted the analyses to a subset of species occurring in an arbitrary percentage of the stations (e.g., Monk et al. 2015). These pre-analysis decisions potentially increase the uniqueness of station clusters but may underemphasize the diversity of species distributional differences observed in the raw data. We chose subarea boundaries that minimized the number of shared station clusters; however, particularly in the vicinity of subarea boundaries, there may be more conflicts than between the centroids of the defined subareas ( Figure 5). Based on the results of these clustering procedures, we identified six subareas of the GoM continental shelf ( Figure 5) for which we contrast relative fish abundance, species dominance, and resilience metrics.
There was a statistically significant subarea effect in log-transformed multispecies abundance data (log e [CPUE + 0.1]; ANOVA: F = 9.114, P < 0.001), with five significant comparisons of subarea means: NC versus SW (t = 5.449, P < 0.001), NC versus CUB (t = 4.550, P < 0.001), NC versus WFS (t = 3.325, P = 0.023), NW versus SW (t = 3.444, P = 0.01), and NW versus CUB (t = 3.131, P = 0.028). Elasmobranch species (sharks and rays) were most abundant in NC and NW; snappers and groupers were densest on the WFS and off CUB; and the abundance of "other demersal" fishes was similar across the GoM, with the exception of CUB (Figure 6A). Pelagic species (e.g., tunas, Great Barracuda, etc.) were not abundant in the catches of the demersal     Catch-at-length data were characterized by a steep ascending (left) limb to a maximum value of about 40-50 cm (Figures 6B, 7). The ascending limb of catch at length reflects incomplete size selectivity due to the relatively large fish hooks used and also reflects the tendency for low-trophic-level fishes to be relatively small. Maximum densities occurred at similar sizes for all subareas (Figure 7), but the slope of the descending limb for the SGoM was significantly steeper (Tables 2, 3) than those for the NGoM and CUB. The NW and NC subareas had shallower slopes than the other four subareas (Figure 7; Tables 2, 3). Despite having relatively high densities of snappers and groupers ( Figure 6A), the WFS had a steeper size spectrum slope than the NC and NW subareas of the GoM, primarily due to the relatively low numbers of small-bodied sharks (e.g.,~1 m in length; Figure 6A; Tables 2, 3) in the WFS subarea. Adjusted means of size spectrum regressions were not different for any of the pairs of subregions having common slopes (e.g., NC versus NW: P = 0.506; WFS versus YP: P = 0.284; WFS versus CUB: P = 0.173; YP versus CUB: P = 0.832; Figure 7). These differences suggest lower overall productivity in some subareas, which was consistent with patterns in overall abundance ( Figure 6A).
Species richness R was calculated for each subarea ( Figure 8) and for 3-cm length intervals within the three super-areas ( Figure 6C). Numbers of species at length should be greatest at small fish lengths because although all teleosts transition from a few millimeters to larger asymptotic lengths, few species actually reach large asymptotic lengths. The ascending limb of the species richness curve ( Figure 6C) was due to the size and species selectivity of the sampling gear for high-trophic-level predators with relatively large asymptotic sizes. Paradoxically, the number of species at similar lengths (betweeñ 50 and 150 cm) was much higher for the NGoM despite the more tropical environment in the SGoM and CUB ( Figure 6C). Even though the number of samples obtained from the NGoM was substantially greater than the number from southern areas, regional differences in R are unlikely to be attributable to sampling intensity, as all species accumulation curves reached their asymptotes well before the maximum number of stations per subarea was reached (Figure 8, top panel). For example, in the NC region, the total number of species encountered was observed in only 34 of the 144 stations sampled. Deterministic and stochastic species accumulation curves showed the same endpoints (total number of species encountered in each subarea) but in some cases differed in slope to projected or realized asymptotes (Figure 8), typical of such analyses (Gotelli and Cowell 2001).
Besides R (Figure 6C), other metrics of diversity were more reflective of numerical evenness across species, including H′ ( Figure 6D), which also showed an ascending limb to a maximum, with a rapid reduction to low richness at the greatest sizes. The slopes of the diversity spectra were also significantly steeper for comparisons of SGoM versus NGoM, SGoM versus CUB, and CUB versus NGoM ( Figure 6D); the slopes of the regression models were ranked as follows: SGoM > CUB > NGoM, with all results significant at P < 0.05. The H′ values for SGoM and NGoM were similar at moderate size intervals (~50-80 cm; Figure 6D), but with steeper slopes, the H′ values diverged at greater sizes. Since H′ incorporates both species evenness and absolute richness, this may reflect similar relative levels of trophic redundancy for SGoM and NGoM at these moderate sizes.

DISCUSSION
The 166 unique species identified in this study (Table 1; based on the nomenclature of Page et al. 2013) represent 11% of the 1,541 known fish species from the GoM (McEachran 2009). The remainder are estuarine or nearshore, occupy the deep GoM, or are otherwise not available to the demersal longline gear we used for sampling (e.g., highly migratory pelagic species were only sampled incidentally as the gear was being set or retrieved). Our comparative analyses of species richness used only the species we sampled as opposed to all known species in various subareas of the GoM (McEachran 2009). Thus, these analyses should be considered in relative, not absolute, terms (e.g., for R and H′ calculations). Nevertheless, because we used standardized sampling methods, our results do reflect relative abundance and thus sampling probability of occurrences consistently across sampling areas.
Are the relative species richness calculations (Figure 8) reflective of the true diversity differences among subareas? The deterministic and probabilistic approaches to species accumulation (Gotelli and Cowell 2001) yielded the same righthand endpoints (numbers of species encountered or estimated by subarea), but there are some differences in slopes, and the stochastically derived curves approached but did not reach asymptotes (Figure 8). The most obvious discrepancy in interpretation of relative R among subareas between the two methods was for the YP. The slope of the accumulation curve for the YP subarea was steeper than those for the NC and WFS subareas, which exhibited higher cumulative R than the YP (Figure 8). The steep slope for the YP region reflects the high R occurring at a few stations-the maximum R per station in the entire program was 18, occurring in the YP (Figure 3, middle panel). A total of 30 of the 69 species encountered in the YP subarea occurred at just two stations (Figure 8, bottom panel). The stochastic averages for the YP for the first and second samples randomly taken were only 8.0 and 14.9 species, respectively-less than half the numbers actually sampled in the field. Thus, even though the actual field sampling data showed an asymptote after 18 of the 24 sampled stations (Figure 8, top panel), the stochastic curve appeared steeper since average numbers of species accumulated more slowly as remaining species were discovered. Results from stochastic models also should be interpreted cautiously because of sampling design considerations. In our case, simple random selection of stations within subareas for jackknife analyses assumes that the individual samples have an equal probability of species richness. However, as observed in Figure 3 (middle panel), species richness depends partially on sample depth, with lower average richness at deeper depths. The potential for biases is unknown, but random deals of stations may be inconsistent with the sampling design that systematically sampled along transects. Despite differences in the numbers of stations between subareas, interpretations of relative R among subareas appear to be robust, recognizing that there are additional relatively rare species to be encountered in all regions.
Our sampling gear was similar to that used in the demersal (bottom) longline fishery in the NGoM (Scott-Denton et al. 2011). In the Scott-Denton et al. (2011) study, at-sea observers recorded the catch, disposition, fishing characteristics, and effort from 1,503 longline deployments sampled from 2006 to 2009 but confined to the NGoM (with concentration on the WFS). Scott-Denton et al. (2011) sampled about 6 times more longline sets, deployed 13 times more hooks, and caught 6 times more fish than our study for comparable areas (i.e., the NGoM). However, even with this large disparity in sampling effort, the number of unique species we observed in the NGoM (129; Table 1) was 86% of that recorded by Scott-Denton et al. (2011;150 unique species). Most of the disparity in R between the two studies was due to incidental catches of rare species at the extreme righthand side of the species accumulation curves (Figure 8). The observed success rate of the bottom longline fishery (percentage of hooks with fish) was only 5%-less than half that observed by us in the NGoM (11%), probably because our sampling effort was not directed to specific fishery targets (e.g., Red Grouper).
Comparisons between our data and those of Scott- Denton et al. (2011) are important because they document that species assemblage determinations and associated species compositions (Figure 4) are relatively insensitive to sampling effort, even at levels 13 times greater than in our study. Furthermore, Scott-Denton et al. (2011) presented the size composition of Red Grouper, which were the dominant target of bottom longline fishing in the NGoM. The mode (46 cm [18 in]) and length distribution of the 40,992 Red Grouper they sampled were nearly identical to those derived from the 849 Red Grouper we caught (Table 1), thus confirming that the size data we obtained from fishery-independent surveys are representative of the size selectivity of commercial fishing practices (at least in this case). Therefore, we are confident that the patterns in species composition and size distributions and the inferences drawn from them based on our data are robust to sampling variability. Since the data reported by Scott-Denton et al. (2011) were collected prior to the DWH oil spill and ours were collected afterwards, this comparison between surveys also shows that the general compositions of fish assemblages in the NGoM were similar before and up to 7 years after the spill.
The relatively high numbers of stations sampled in the NC region (144 of 343 stations, or 42%; Figures 1, 5) reflect a time series of six surveys conducted in the NC subarea from 2011 to 2017. Although not reviewed in GULF OF MEXICO FISH ASSEMBLAGES detail here, the issue of spatial comparisons (e.g., across subareas) using a multi-year sampling campaign raises the issue of interannual variability in community composition and relative abundance. The overall CPUE of total species catches in the NC was stable over the six surveys, although some species showed modest increases and, conversely, some species decreased in abundance (e.g., Red Snapper and Southern Hake). Importantly, aggregating species into the taxonomic groups resulted in relatively consistent catches without significant trend. The other five subregions doubtlessly have trends in species abundances as well; however, as a snapshot of the relative compositions among regions, there is no reason to believe that our interpretations from multi-year sampling are not robust. This interpretation is supported by the consistency of our catches in the NGoM with those of Scott-Denton et al. (2011), who sampled longline catches in the NGoM prior to the DWH oil spill.
Our data demonstrate classic patterns of alpha (stationto-station; Figures 3, 4), beta (along environmental gradients; Figures 2, 3), and gamma (GoM-wide) species diversity (Figures 4, 5;Whittaker 1972;Anderson et al. 2011). We observed distinct patterns of species zonation by depth (inversely correlated with bottom temperature) throughout the GoM (Figure 2). This is similar to observed depthrelated patterns of demersal fish zonation in deep waters from the shelf break to the abyss in the NW and NC GoM (Wei et al. 2012). Although the patterns of species zonation we observed along the mid-to outer shelf were similarly depth structured (unlike Wei et al. 2012), we also documented significant horizontal patterns in zonation over the LME scale ( Figure 5). Additionally, we observed significant differences in average fish abundance (CPUE) among some subareas and by species group ( Figure 6A). The lack of horizontal contrast in species dominance in deeper waters (Wei et al. 2012) may be because the deep GoM is characterized by more uniform temperature conditions, thus contributing to more similar species distributions, at least in the NC and NW GoM. Along the inner continental shelf of the NW and NC GoM (18-55 m deep), Monk et al. (2015) documented horizontal structuring of these shallow shelf demersal fish communities, with a strong seasonal component and boundaries roughly equivalent to those we defined in overlapping areas along the mid-and outer shelf ( Figure 5). Darcy and Gutherz (1984) also identified both significant depth and spatial (north versus south) differences in fish abundance from 9 to 193 m on the WFS. Thus, fishes on the continental shelves of the GoM are structured both vertically (by depth) and horizontally, whereas deeper-dwelling demersal fish assemblages appear more horizontally homogeneous.
In addition to depth-related patterns in richness and diversity, there were also substantial differences among subareas (Figures 4, 5, 8). The NW and CUB subareas had significantly fewer species than the other areas (Figures 6, 8). Relatively low species richness for reef fishes on the NW shelf has previously been reported (Dennis and Bright 1988;Streich et al. 2017). Several theories potentially explain the relatively low diversity there, including the presence of a widespread nepheloid layer (Dennis and Bright 1988;Streich et al. 2017) that may limit benthic productivity, low habitat diversity (Dennis and Bright 1988), and the generally lower levels of primary productivity in the NW GoM (Benway and Coble 2014). The low number of species encountered off CUB is a paradox given that the area we sampled contains abundant shallow and mesophotic coral reefs (Claro et al. 2001), which are generally thought to be highly diverse. However, while coral reefs have highly diverse fish assemblages, only a portion of the fishes encountered along the northwest CUB coast were vulnerable to our longline sampling gear (e.g., animals ≥ 50 cm and high-trophic-level predators; Figure 6B). Furthermore, subsistence fisheries off northwest CUB are intensive (Claro et al. 2001;Tunnell 2017;S. A. Murawski and M. Armenteros, personal observations). Baisre (2018) concluded that 79% of Cuban marine fishery stocks were overfished or collapsed. This may limit the availability of large, economically valuable fishes to the sampling gear. The relatively steep slope of the size spectrum off CUB ( Figure 6B) is consistent with this argument, as is the negative offset of the data and regression slope for H′, which are below the data for SGoM and NGoM ( Figure 6D).
There are strong theoretical and empirical bases for using the multispecies size spectrum of fish communities to index both within-ecosystem temporal responses to drivers such as fishing effort (Murawski and Idoine 1992;Rice and Gislason 1996;Shin and Shannon 2010) and to compare status across ecosystems (Bianchi et al. 2000; Figure 6B). The slope of the multispecies size spectrum integrates the abundance, growth, and mortality of species and potential compensatory responses, such as competitive release and niche replacement of depleted populations (Murawski and Idoine 1992). If the productivity or mortality of ecosystem components differs among areas, then these should be reflected in differing slopes and perhaps adjusted means of the size spectra.
What do differences in size and diversity spectra tell us about the relative resilience potential of fish assemblages occupying GoM subareas? Overfishing and collapse of keystone species have long been recognized as having destabilizing effects on marine ecosystems (Fogarty and Murawski 1998;Steneck et al. 2002;Hughes et al. 2007;Newton et al. 2007;Mumby and Steneck 2008). Ecosystems in which biomass is distributed broadly across varying animal sizes are considered to be more resilient than ecosystems in which biomass is concentrated at the lowest trophic levels (Sprules and 342 Munawar 1986;Pope et al. 1988;Jennings et al. 2001;Mumby and Steneck 2008). Differences in the steepness of size spectra reflect both differences in fishery management outcomes and species dominance differences among subareas ( Figures 6D, 7). The absolute values of spectrum slopes are determined by species abundance and are related both to low-trophic-level productivity and the total mortality on the species complex (Rice and Gislason 1996). To a point, conservation of the size spectrum slope may indicate potential compensatory feedbacks, such as increased individual growth rates and lowered natural mortality as a function of reduced densities of selected species, and thus may serve as a metric of inherent ecosystem resilience.
Slopes of the size and diversity spectra in the NGoM were significantly shallower than those of other two superareas ( Figure 6B, D; Tables 2, 3). A number of factors likely contribute to this result. Since 1999, the proportion of exploited species in the NGoM undergoing overfishing has declined from 40% to less than 5% (Karnauskas et al. 2017), with many species increasing in biomass and extending size distributions as a function of reduced fishing mortality. The high abundance of small-bodied elasmobranchs (~1 m TL), particularly in the NW and NC subareas ( Figure 6A), is also likely a major contributor to the shallower size spectrum slope in the NGoM. Sharks in the NGoM have been more intensively regulated in recent years through restrictive catch limits, closed areas, and, in some cases, prohibitions on landings (NOAA Fisheries 2017). Overall, it is likely that the effectiveness of fishery management regimes extant in the three super-areas scales the relative slopes of the size spectra and, to an extent, the diversity spectra.
The NC subarea had the shallowest size spectrum slope ( Figure 7) and the highest R (Figure 8) of all subareas. Ironically, this is also the subarea in which the DWH spill occurred (Figure 1), where much of the U.S. oil and gas offshore infrastructure exists (Murawski and Hogarth 2013;BOEM 2017), and where the highest frequency of GoM hurricane disturbance has occurred over the past 150+ years (NOAA 2017). In contrast, the NW subarea, although having a shallow slope for its size spectrum (Figure 7), exhibited a low level of R (Figure 8). Thus, even though these subareas are adjacent, the NW subarea appears less resilient overall to the effects of pulse-type perturbations, such as a catastrophic oil spills, primarily due to less functional redundancy in the event of the collapse of one or more important species. Like the NC subarea, a large proportion of U.S. oil and gas wells are located in the NW subarea (Murawski and Hogarth 2013;BOEM 2017), and therefore this subarea may be particularly vulnerable to the effects of pulse-type perturbations. The two subareas in the SGoM had relatively steep size spectra ( Figure 7) and moderate levels of R (Figure 8).
Hurricane frequency is lower in the SW subarea than in the YP (NOAA 2017), but the SW also supports most of the offshore Mexican oil and gas industry. In the SW subarea, a large exclusion zone (17,500 km 2 ) was set around the Campeche oil platform area, restricting all fish activities and non-oil industry-related navigation (DOF 2003) from 2003 until the zone was partially modified in 2017. However, the reserve effect does not seem to be reflected in the SW size spectrum (Figure 7), which was the steepest observed among all subareas. This may be due to the high total fishing effort (artisanal plus industrial) on mobile resources there that are in fully exploited or overexploited conditions (D ıaz de Le on et al. 2004). The CUB subarea had relatively low R and a steep size spectrum (Figure 7; Table 2), indicating that it too may be vulnerable to acute resource perturbations. Because our sampling occurred after the DWH and Ixtoc-I oil spills, it is difficult to ascertain whether the NC and SW subareas, respectively, were significantly more resilient prior to those incidents.
Spatial comparisons of fish biodiversity across the breadth of the GoM are possible because we used consistent gears and sampling designs in all subareas. Obviously, all fish sampling gears (e.g., longlines, midwater and bottom trawls, gill nets, hook and line, camera-based surveys, ichthyoplankton nets, etc.) have biases. This is clear with respect to the size selectivity ( Figure 6B, C) and species selectivity ( Figure 6A) of demersal longlines. How would the interpretation of abundance, size structure, and diversity differ if we used any of these other technologies in a comprehensive regionwide survey? There is some information for such gear comparisons from systematic sampling programs in the NGoM using trawls (Monk et al. 2015) and alternative commercial fishing gears (Scott-Denton et al. 2011), among other published studies. Trawl catches are dominated by small-bodied, lowtrophic-level fishes and invertebrates, which are generally underrepresented by demersal longlines (Monk et al. 2015; Table 1). Conversely, trawling in the NGoM captures relatively few of the large-bodied, fast-swimming fishes that are the primary targets of commercial fisheries (e.g., snappers, groupers, and other species). This is partially due to the inappropriateness of trawl gear for the high-relief habitats where some of these species congregate. These and other sampling techniques are thus complementary, and the choice of the most appropriate technology depends upon what questions are being asked. A systematic, multi-gear sampling scheme across the GoM-for example, using bottom and midwater trawls, longlines, camera systems, and ichthyoplankton gears -would provide for a more complete, complementary, and potentially powerful set of data with which to understand the biodiversity, connectivity, and ecosystem dynamics of fishes in the entire GoM.
Our study compared continental shelf fish communities in the GoM, emphasizing economically and ecologically important demersal species (Murawski et al. 2016). Due GULF OF MEXICO FISH ASSEMBLAGES to the importance of these target species in all GoM subareas, there is concern that the cumulative impacts of an array of stressors may trigger communitywide regime shifts in abundance and species dominance. The three countries bordering the GoM all have differing fishery management regimes, resulting in varying degrees of success in achieving sustainable fisheries (Tunnell 2017). Harmonizing management goals and strategies among countries for interconnected species and communities now occurs only for those stocks that are managed under regional fishery management organizations (e.g., the International Commission for the Conservation of Atlantic Tunas). Many species occurring in GoM countries connected by animal migrations or larval dispersal are not subject to formal management coordination. Additionally, all GoM countries are actively pursuing offshore and, increasingly, ultra-deep (≥1,500-m) oil drilling, which may be problematic if oil spills on the scale of the Ixtoc-I (1979Ixtoc-I ( -1980 and DWH (2010) spills occur again, particularly in areas of the GoM with low fish community resilience potential (i.e., deepwater communities; Koslow et al. 2000) and where connectivity among subareas of the GoM is also low. Understanding the relative sensitivity of subareas to perturbations and strengthening the resilience of fish assemblages through conservative fishery and environmental management policies can help mitigate risks to ecosystem stability and coastal economies at both the subarea scale and the interconnected LME scale.