Time-Varying Epipelagic Community Seascapes: Assessing and Predicting Species Composition in the Northeastern Pacific Ocean

The vast spatial extent of the ocean presents a major challenge for monitoring changes in marine biodiversity and connecting those changes to management practices. Remote-sensing offers promise for overcoming this problem in a cost-effective, tractable way, but requires interdisciplinary expertise to identify robust approaches. In this study, we use generalized additive mixed models to evaluate the relationship between an epipelagic fish community in the Northeastern Pacific Ocean and oceanographic predictor variables, quantified in situ as well as via remote-sensing. We demonstrate the utility of using MODIS Rrs555 fields at monthly and interannual timescales to better understand how freshwater input into the Northern California Current region affects higher trophic level biology. These relationships also allow us to identify a gradient in community composition characteristic of warmer, offshore areas and cooler, nearshore areas over the period 2003–2012, and predict community characteristics outside of sampled species data from 2013 to 2015. These spatial maps therefore represent a new, temporally and spatially explicit index of community differences, potentially useful for filling gaps in regional ecosystem status reports and is germane to the broader ecosystem-based fisheries management context.


INTRODUCTION
Major initiatives worldwide have recognized the importance of measuring diversity and community structure as indicators of ecosystem condition (Skidmore and Pettorelli, 2015). Satellite-based remote sensing (RS) is a tool that provides habitat information across large extents at high spatial and temporal resolutions, having the potential to describe the distribution of multiple facets of biodiversity, including species distributions, alpha diversity, predator-prey overlap, and community composition (Skidmore and Pettorelli, 2015;Pittman, 2017;Wallis et al., 2017;Walters and Scholes, 2017;Muller-Karger et al., 2018). As a greater variety of RS-based oceanographic information becomes increasingly available and accessible, scientists have coupled these data with biological datasets collected in situ to develop knowledge in key areas of seascape ecology research, which are in turn useful for fisheries management and dynamic ocean management (Zwolinski et al., 2011;Hazen et al., 2013;Scales et al., 2017). Three key topic areas of seascape ecology research [i.e., research linking oceanography with landscape ecology (Pittman et al., 2011)] include: (1) a better understanding of how terrestrial landscapes affect adjacent seascape patterns and processes, (2) determining which structural attributes of seascapes (defined in Pittman et al., 2011, as: 'wholly or partially submerged marine landscapes') drive biotic assemblages and distribution of biodiversity, and (3) quantifying the impacts of climate change on seascape patterns (Pittman et al., 2011). In this study, we address the first two areas and demonstrate potential applicability of our approach to the third.
Pelagic fish and invertebrate communities are among the most ecologically, culturally, and economically important components of marine ecosystems worldwide (Pikitch et al., 2014). The species we focus on here (salmon, sardines, anchovies, squid and mackerel) are at once the object of dedicated and emerging fisheries and critical links connecting lower and higher trophic levels in the coastal ocean (Pikitch et al., 2014;Szoboszlai et al., 2015). As such, understanding the dynamics of pelagic communities is fundamental to ecosystem-based fisheries management efforts. As in other regions, the pelagic community composition of the California Current Large Marine Ecosystem (CCLME) has varied substantially in response to changes in environmental conditions (Brodeur et al., 2006;Ralston et al., 2015;Peterson et al., 2017;Morgan et al., 2019). One of the major environmental drivers influencing the coastal ocean off Washington, Oregon, and northern California is the input of freshwater from the Columbia River (e.g., Hickey et al., 2005;Henderikx Freitas et al., 2018) and the small coastal rivers along the coast (Mazzini et al., 2014;Saldías et al., 2020). These freshwater outflows are a significant source of nutrients, sediments, organic matter and other constituents for the coastal ocean (Sigleo and Frick, 2007;Goñi et al., 2013). The Columbia River plume has been identified as a crucial environment providing nutrients and enhancing the coastal productivity to sustain the ecosystem during periods of delayed upwelling (Hickey and Banas, 2008;Kudela et al., 2010). Thus, the plume can modulate plankton distribution and the aggregation of zooplankton around plume fronts (e.g., Morgan et al., 2005;Hickey et al., 2010). However, the use of RS data products related to coastal freshwater input in the study of epipelagic community dynamics is very limited.
Moving toward ecosystem-based fisheries management of the CCLME requires consideration of the effects of environmental variability and climate on the biota (Field et al., 2006), including the development and evaluation of ecological indicators that directly and indirectly measure environmental impacts on marine communities (Levin et al., 2009;Samhouri et al., 2014;Thompson et al., 2019). Under current marine management policies, predictions of species and ecosystem distributions are necessary for ongoing integrated ecosystem assessments and indicator development . Leveraging RS data, predictive statistical models can allow for near real-time and hind-casted distributions of individual species as well as community metrics onto multi-layered RS fields. Innovative methods are being used to couple high-resolution RS data with field survey data via statistical modeling techniques to generate fine-resolution predictive maps (Manderson et al., 2011;Hobday and Hartog, 2014;Thorson et al., 2020). However, research that integrates satellite RS data with species diversity, has largely focused on alpha-diversity metrics (e.g., Pittman et al., 2007;Hazen et al., 2013) and less on other community metrics, such as ordination scores representative of assemblage similarities/dissimilarities. Ultimately, forecasts of individual species distributions and community composition metrics can provide an important foundational layer to be used in marine spatial planning, such as dynamic ocean management Maxwell et al., 2015 and references therein).
In this study, we use a pelagic seascape approach and model the epipelagic fish and macro-invertebrate community structure off the coasts of Oregon and Washington using RS and in situ oceanographic biophysical data. First, we assess the relationships among species across different months and community groups. To understand the relationship between community structure and the environment, we use generalized additive mixed effects models (GAMMs) and couple ordination scores with both local in situ collected environmental data and remotely sensed oceanographic data fields. Using the results of these community seascape models, we construct predictive spatial maps of community composition onto the RS layers as hind-casts and predictions at yearly and seasonal scales. From these community structure prediction maps, we extract time series of predicted community gradients. We propose that by coupling readily available RS data with community metrics, we can predict community species composition in the pelagic marine realm, gaining information that can potentially be of use for dynamic regional ocean management.

Survey Specifics
From 1998 to 2012 NOAA's Northwest Fisheries Science Center's Fish Ecology Division regularly conducted pelagic surface trawl surveys within the Northern California Current (NCC) epipelagic ecosystem with the primary goal of quantifying the marine survival and habitat associations of juvenile salmonids [Juvenile Salmon and Ocean Ecosystem Survey (JSOES); . These surveys were conducted regularly each year in spring-early summer, and late summer, off the coasts of Oregon and Washington between 44.25 and 48.23 • N and −125.61 and −123.97 • W ( Figure 1A). Here, we analyze survey data collected between 2003 to 2012 in order to synchronize with MODIS data availability. The sampling methodology of these surveys was altered after 2012 with respect to the geographical and temporal extent of sampling, and modifications to the gear, impacting catchability and abundance estimates of certain species (Wainwright et al., 2019). However, while comparable ecosystem data (with similar gear configuration) from 2013 to 2015 were collected, they were not made available to us for this analysis/validation purposes. Approximately, 10 transects consisting of between 6 and 8 stations extend from the nearshore to off the continental shelf (approximately 50 km from shore). Thirty-minute surface tows were conducted at each station with a Nordic 264 rope trawl (Nor'Eastern Trawl Systems, Inc., Bainbridge Island, WA, United States) at a speed over ground of approximately 6 km/h. The length of the trawl is 108 m long with a mouth opening of 30 m × 20 m, and the head rope at approximately 1 m from the surface. The mesh size varies along the length of the rope trawl, ranging between 162.6 cm at the mouth to 8.9 cm at the cod end, which had a 0.8 cm liner inside.

Species Data
For this study, we restricted the hauls to only those conducted exclusively during daylight hours of the months of May, June, and September. The cruises occurred during the last half of the cruise month, however, on few occasions, some trawls were conducted into the first few days of the next month (into July for example for the June cruise, similarly for May and September). Due to the mesh size used in the trawl, we restricted our community analysis to those species quantitatively retained by the trawl, including both invertebrate and fish species. All species captured in the trawl were identified, counted and measured, with the exception of large catches, where a random sub-sample of fish from each species was measured, counted and weighed and then total abundance re-calculated based on total weight. Catches were standardized by calculating catch-per-unit-effort (CPUE, number per km 2 ) using the distance fished (geographic distance between beginning and end of trawl) ( Figure 1A; Daly et al., 2012). Of the total quantitatively sampled species captured between 2003 and 2012 in May, June, and September, we restricted the epipelagic community analysis presented here to species that were found in 3% or more of the sampled stations.

In situ and Remotely Sensed Oceanographic Data
In situ (IS) ocean temperature, salinity, and pressure were recorded at each meter of the water column down to 5 m from the bottom, or to 100 m depth, either prior to or immediately after conducting each surface trawl at a station using a Sea-Bird SBE 19 SeaCAT conductivity/temperature/depth profiler (CTD). At each trawl station (Figure 1A), water samples were collected with a Niskin bottle at a depth of 3 m, and filtered at sea. Pigments were later extracted and analyzed using acetone (90% HPLC grade). Sample fluorescence was measured with a Turner Designs fluorometer (Arar and Collins, 1997). Following Lentz (1992), we calculated mixed layer depth as the depth where the temperature difference relative to the sea surface exceeded 0.02 • C.
Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua ocean temperature and multi-spectral color fields for the region off Oregon and Washington were acquired from http: //ocean.color.gsfc.nasa.gov/. Specifically, sea surface temperature (SST), chlorophyll-a (CHL), and Remote sensing reflectance at 555 nm (Rrs555) fields were used. SST and satellite derived metrics of primary productivity have been commonly used to define marine habitat for pelagic species as they describe both the thermal conditions as well as potential prey availability (Suryan et al., 2012). Rrs555 fields have been previously used to map the surface expansion of the Columbia River plume throughout the annual cycle (Thomas and Weatherbee, 2006;Mazzini et al., 2015;Saldías et al., 2016) and is likely a good descriptor of community differences due to biological responses to plumeinfluenced waters (Saldías et al., 2016). All MODIS fields were smoothed using a two-dimensional 3 × 3 grid cell median filter to reduce noise associated with cloud edges and to enhance frontal features in satellite images (Wall et al., 2008;Saldías et al., 2016).
Sea level anomaly (SLA) data were obtained by combining gridded, daily AVISO altimeter fields (Pujol et al., 2016) with low-pass filtered coastal tide gauge data (Risien and Strub, 2016). This 0.25 • × 0.25 • , blended dataset improves SLA estimates along the coast by removing altimeter data approximately 55-70 km from the coastline and replacing it with a linear interpolation between the tide gauge and remaining offshore altimeter data (Risien and Strub, 2016). In order to have a consistent grid resolution for all RS data fields, we interpolated the blended SLA data to the 4 km MODIS grid using a Barnes objective analysis scheme (Barnes, 1964;Emery and Thomson, 2004). As we were interested in quantifying the distinct water properties associated with distributions of species in the local community, we used both SST and SLA in this study (as done in other species distribution studies, see Becker et al., 2016). SLA and SST are correlated at large scales (Emery et al., 2011) and in the data used here. While MODIS SST is limited to the skin of the sea surface (top 10 microns), SLA is integrative, providing information about the density structure of the water column. As SLA is related to pycnocline depth it may impact habitat quality and availability for distinct species (e.g., a higher SLA is indicative of a deeper pycnocline and therefore potentially more diluted prey availability for pelagic fish and invertebrates, whereas a lower SLA is indicative of a shallower pycnocline which may in turn concentrate prey for pelagic species).
We collocated sampling dates (using actual trawl dates, rather than assigned cruise month) and locations (midpoint of the start and end locations) of surface trawls (Figure 1) with coincident 8-day MODIS and SLA data to the nearest 4 km grid cell for all fields. As we were interested in understanding the fluctuations of the community associated with the regional features (upwelling front, Columbia River plume, eddies) at monthly and interannual time scales, we used 8-day MODIS composites, which McKibben et al. (2012) show to be adequate for comparing IS with RS data off central Oregon. While the sampled species and in situ time series begin in 1998, the time series analyzed here begins in 2003 to allow for the use of consistent RS fields and temporally congruent IS data. Prior to 2003, AVHRR SST and SeaWiFS ocean color data were collected by two distinct satellites with different orbits and timings. The MODIS-Aqua instrument collects SST and ocean color information simultaneously, thus avoiding the issue of mismatched data in time and space. Additionally, SeaWiFS data are only available through mid-December 2010, which does not cover the in situ sampling time series (up to 2012) or the prediction period (up to 2015). To better understand the correlations of IS and collocated RS fields across all stations, all months, and all years, we ran Pearson correlations and report the correlation values (Supplementary Figure 1). In order to create maps of predicted community gradients and to maximize the number of cloudfree grid cells,

Community Analysis of the Epipelagic
Non-metric multidimensional scaling (NMDS) ordination (Clarke and Warwick, 2001;McCune and Grace, 2002) using Bray-Curtis dissimilarities was used to quantify the variability in community composition by year and season of the epipelagic juvenile and adult nekton and associated macro-invertebrate community. We specifically chose to use NMDS ordination to quantify the community composition as we do not assume that species are linearly related to each other [a central assumption of principal component analysis (PCAs)/empirical orthogonal functions (EOFs)] (McCune and Grace, 2002).
To explore the spatial changes of the epipelagic nekton and associated macro-invertebrate community by year and season, we performed a NMDS ordination of hauls in species space with a random starting configuration, calculating a dissimilarity matrix using the Bray-Curtis distance measure. The community dataset was composed of 23 relative species abundances (CPUE) by 1215 individual hauls as sample units. The 23 species consisted of salmon, forage fish, and some top predator species, as well as several gelatinous macro-zooplankton and squid representing the macroinvertebrate community. The relative species abundances (CPUE) were log 10 (x + 1) transformed to reduce the variation between different species abundances. In order to visualize the relationships of in situ and RS data in relation to axes of the NMDS ordination -we used vegan's envfit (Wood and Scheipl, 2014) function to first obtain a vector overlay on the ordination representing a linear regression of the environmental variables with each NMDS axes. Next, to visualize potential non-linear relationships between each environmental variable and the ordination, we used the ordisurf function that fits a generalized additive model with each variable individually, and produces response surfaces of oceanographic conditions overlaid on the ordination using 2-D thin-plate regression splines. The degree of smoothing is automatically selected by cross validation (Wood and Scheipl, 2014).
We performed a hierarchical cluster analysis to define two community clusters (refer to dendrogram of hauls in Supplementary Figure 6), henceforth referred to as cold and warm communities as they were distributed along the temperature gradient associated with NMDS axis 1. Multiple Response Permutation Procedure (MRPP) was used to determine if the two clusters were significantly different from one another (refer to McCune and Grace, 2002 for a methodological description of MRPP) (A-statistic = 0.056, p < 0.01). MRPP was also conducted by year and month to determine whether these two clusters were significantly different from one another across all years and months (A-statistic range: 0.046 to 0.057, p < 0.01). To assess the exclusivity and fidelity of the species to particular clusters, as well as the statistical significance of the relationship between species abundance and groups of sites in each month, we performed an indicator species analysis (ISA) (McCune and Grace, 2002) on the two groups identified a priori by the hierarchical cluster analysis for each season of data. All non-parametric multivariate statistical analyses were done in R v. 3.4.0 (R Core Team, 2019) using the vegan and indicspecies libraries (Oksanen et al., 2019; R Core Team, 2019; De Caceres and Jansen, 2020).

Relating Community Composition to in situ or Remotely Sensed Oceanographic Data
We used the dimensionless NMDS axis 1 scores from the ordination of hauls in species space as response variables in GAMMs (Zuur et al., 2009) to explore the non-linear effects of oceanographic conditions on the community structure. The GAMM analyses focus on the cold to warm community gradients and how the Columbia River plume affected species communities. We constructed two model set ups, using NMDS axis 1 as the response variables and either in situ and remotely sensed variables as co-variates (i.e., NMDS1∼IS covariates and NMDS1∼RS covariates) ( Table 3). We selected the best fit models for each based on minimization of the AIC. All final models and model selection details for NMDS1 are reported in Table 3.
Only the community seascape model (NMDS1∼RS) results are discussed here, the other model (NMDS1∼IS) and their resulting functional relationships with environmental co-variates are presented in the Supplemental Information.
Using the linear and non-linear functional relationships of RS environmental data with NMDS axis 1 derived from the GAMM, we predicted the spatial distribution of the response variable (NMDS axis 1) for the different years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) and month (May, June, and September) combinations, using monthly composited satellite images. We also used the relationships of NMDS axis 1 to RS covariates, to predict outside the range of the data (2013-2015 in May, June, and September) used here. Spatial predictions were restricted to within the minimum convex hull surrounding the sampled locations.
Generalized additive mixed effects models are a non-linear regression technique where the covariates (environmental variables) are modeled with non-parametric smoothing functions (Hastie and Tibshirani, 1990;Wood, 2015) in the gamm4 library of R (Wood and Scheipl, 2014). GAMMs do not require an a priori assumption of the type of relationship between the covariates and response variables in a model. The environmental co-variates were used as the fixed effects in the models and individual stations as the random effect to account for within station correlation. Variable selection was achieved by minimizing the Akaike information criterion (AIC) as well as the genuine cross validation (gCV) score for each set of models. To calculate gCV, we fit each model to a randomly selected training dataset (with 90% of the total observations), generated predictions for a validation dataset (with the remaining 10% of the observations) and then calculated the prediction error. This procedure was repeated a total of 500 times and a model gCV criterion score was computed by taking the mean squared prediction error (Ciannelli et al., 2007(Ciannelli et al., , 2012. The candidate model with the lowest combined AIC (weighted more heavily) and gCV criterion scores was determined to be the best model for determining local IS and RS data relationships with NMDS1.

Oceanographic Correlates With Individual NCC Epipelagic Species
For each species (at stations with positive catch only), we provide the median and range values of remotely sensed oceanographic data (SLA, SST, CHL, Rrs555) in Table 1. The median and range of in situ sampled oceanographic variables (temperature, salinity, density at 10 m depth, and chlorophyll-a (CHL) concentration sampled at 3 m depth), for all 23 species included in this analysis are included in Supplementary Table 1.

Interannual and Monthly Community Variability
The NMDS 2-D ordination of the survey hauls (1215 total hauls and 23 species) had a final stress of 0.42 (Figures 2A-D). There were significant differences in the communities of fish and macroinvertebrates sampled off Oregon and Washington amongst years (data from months and stations for a given year) ( Table 2). However, the degree of differentiation between years was low given the large variability associated with spatial variables and ocean conditions (e.g., Oregon vs. Washington and cold vs. warm conditions) ( Figure 2B). The overall yearly MRPP A-statistic was 0.013 (p < 0.05) for year-to-year comparisons when using sample units from the full sampling grid in each month.

Indicator Species of Monthly Community Clusters
Indicator species of monthly communities are presented in Table 2. Pacific saury was the only consistent warm community indicator across all three cruise months. The other warm indicator species differed by month ( Table 2). The months of May and June included egg-yolk jelly and sablefish as indicators of the warm community, while in June, Jack mackerel (Trachurus symmetricus) was also an indicator of this warm  Component 'A' is the probability that the surveyed sample unit belongs to the target site group given the fact that the species has been found (or specificity) and 'B' is the probability of finding the species in sites belonging to a particular cluster (or fidelity) and 'stat' is the indicator value (composed of both components).
community cluster. The warm community in September included Pacific saury and Jack mackerel as well as ocean sunfish. The indicator species of the cold community cluster were more variable than those of the warm cluster.  Table 2). Market squid was a significant cold indicator species in September but not in other months. Steelhead was also found to be an indicator species for the cold community only for the month of May. Species differed by month in their specificity, probability and fidelity to a particular community cluster ( Table 2).

In situ Environment Relationships With Community Composition
With respect to the in situ sampled environment, the ordination had a notable cold-warm gradient along NMDS axis 1 and a salinity gradient along NMDS axis 2 (henceforth referred to as NMDS1 and NMDS2). Simple linear regressions indicated that NMDS1 was positively correlated with in situ temperature at 10 m depth and the seafloor depth of the sampling location, which was also partially negatively correlated with NMDS2 (Figure 2A). In situ salinity was positively correlated along NMDS2 while mixed layer depth was positively correlated along both NMDS1 and NMDS2. In situ chlorophyll (at 3 m depth) was negatively correlated with NMDS1 and NMDS2 (Figures 2A, 3A-H). The surfaces of each in situ physical and biological variable overlaid on the NMDS ordination presents the non-linear component of each variable in relation to the community composition (with both NMDS axes) of each haul the ordination (Figures 3A-H). Temperature at 10 m depth had higher temperatures at positive NMDS1 scores and lower temperatures at negative NMDS1, while having both high and low values of temperature at both extremes of NMDS2 (Figures 3A-H). Salinity was more linearly related to NMDS2, indicating a gradient of fresh to warm waters with a positive correlation along NMDS2 (Figures 2A, 3A-H). Chlorophyll had a negative relationship with both NMDS axes (Figure 2A). While oxygen was only sampled from 2006 to 2012, we overlay the non-linear relationship with the community (Figure 3D), but do not use it in the analyses. Communities sampled at shallower stations were concentrated in the upper left quadrant of the ordination (low values of NMDS axis 1 and high NMDS2) whereas communities sampled the deeper water column stations were on the right side of the NMDS ordination (all NMDS2, and NMDS1 values > 0) ( Figure 3E). Mixed layer depth largely showed a positive linear relationship with both NMDS axes ( Figure 3F). Relationships with latitude, longitude are also presented (Figures 3G,H).

Remotely Sensed Environment Relationships With Community Composition
The relationships of the RS data fields showed slightly different patterns than the IS data as expected, given the different data sources and depth at which the variables were measured relative to satellite data measurements (see Figures 3, 4). MODIS-Aqua SST had a different pattern compared to the in situ data at 10 m depth. Relatively warm SST (>14 • C) ranged between −0.5 to 1.5 values along NMDS1 and along almost the entire NMDS2 ( Figure 4A). SST of less than 13 • C was restricted to the far-right side of the ordination (NMDS1 values less than −0.5) ( Figure 4A). Rrs555 showed a similar pattern to in situ salinity in relation to the ordination, with high values in the lower right quadrant of the NMDS (Figure 4B). MODIS Chl-a demonstrated an almost identical pattern to that observed with in situ Chl-a collected at 3 m depth ( Figure 4C).

Community Seascape Model Results
The final RS NMDS axis 1 model explained 45.8% of the variability in NMDS axis 1 and included as covariates: (1) an interaction of latitude and longitude, (2) an interaction of SLA, SST, and (3) Rrs555 ( Table 3). The average spatial pattern of NMDS axis 1, with low NMDS axis 1 values distributed along most of Washington and higher NMDS axis 1 values distributed off Oregon (Figure 5A). The GAMM interaction biplot of SST and SLA showed a positive relationship of NMDS axis 1 with both SST and SLA ( Figure 5B). There was a non-linear relationship of NMDS1 to the Rrs555 data field, with high NMDS1 at low values of Rrs555 (<0.004 sr −1 ). At Rrs555 values >0.004 sr −1 there was a quasi-linear relationship with negative NMDS axis 1 scores ( Figure 5C).

Remote Sensing-Based Prediction of Ordination Scores
The best fit NMDS1∼RS model (community seascape model) had significant correlations with the interaction term of latitude and longitude, interaction of MODIS SST with AVISO+TG SLA, and with Rrs555 ( Figure 5). Positive values of NMDS1 were generally found offshore, and at higher SST and SLA values and vicea-versa for negative NMDS1 values. While a formal threshold analysis was not done here, the Rrs555 value where there is no change in NMDS1 occurs at approximately 0.004 sr −1 . Based on the community seascape model (NMDS1∼RS; Figure 5), spatial predictions of NMDS axis 1 values (indicating the gradient from cold to warm communities) indicated differences across years and months (Figure 6). Based on boxplots of extracted predicted NMDS axis 1 values for the Oregon and Washington regions as well as the whole sampled region, Oregon had a higher presence of warm community during all months-years (Figure 7), while Washington had a greater presence of the cold community across all months-years ( Figure 5). During some months and years (May 2005, May 2008, June 2008, June 2012, June 2006, September 2006, September 2012, June 2014, the difference between Oregon and Washington was pronounced (Figure 7).
In May, the period between 2003 and 2005 had a higher presence of the warm community across the entire sampled extent, with a greater presence of the warm community in Oregon than for Washington or the full region (Figure 7). For this month, from 2006 to 2012, there was a lower presence of cold community across the whole region (Figure 7). For May 2013-2015, a significant increase in the presence of the warm community was predicted for all regions, with 2015 having the highest presence of the warm-offshore community within our hind-casted and fore-casted time series. In June, between 2003 and 2005 there was a greater incursion of the warm community than other years (Figure 7). Years with little to no presence of the warm community included 2007 and 2008 (Figures 5, 6). For both Oregon and Washington, September of 2004 and 2009 had the highest predicted warm community presence.
Between 2013 and 2015, for May and September there was a pronounced predicted increase in the presence of the warmoffshore community, especially off Oregon (high positive values of NMDS1). June and September of 2014 both had notably higher predicted presence of the warm community relative to the full area and to the Washington area. For all regions, June 2015 was similar to the long-term mean predicted NMDS1 scores. However, in both May and September, the predicted NMDS1 scores were higher than any previous year and had values greater than the most positive of the observed NMDS1 scores (up to 2012).

DISCUSSION
The field of pelagic seascape ecology so far has dealt mainly with the prediction of species distributions onto RS data fields (e.g., Alvarez-Berastegui et al., 2014, as well as the long-standing tradition of categorizing ocean waters based on RS data (Jerlov, 1968;Kavanaugh et al., 2014). More recently, predicting individual species distributions onto clustered RS data fields has also been attempted. For example, Breece et al. (2016) classified MODIS ocean color and SST fields, and then predicted Atlantic sturgeon habitat (using acoustic receiver tracked movements of individuals) onto categorized images based on a cluster analysis of multiple RS data fields. Here, we chose a different seascape approach: prediction of community gradients onto RS data by linking ordination axes describing community composition onto statistically amalgamated RS-data fields given the relationships defined by non-linear regressions (GAMMs). The relationships between marine species distributions and the environment are often non-linear (McCune, 2006;Lintz et al., 2011), as are the relationships of gradations in community composition to the environment. Our use of a continuous community gradient response variable (ordination axes) avoids drawbacks associated with classified community clusters, which includes unrealistic crisp borders and potential misclassification of species to distinct clusters (Lintz et al., 2011). It also allows for the description of non-linearities between individual RS data fields and community gradients (e.g., Figures 4A-D). Our approach of using a community metric has been instrumental in uncovering useful ecological information that other approaches could not have revealed. Notably, we found that at the community level, spatial and temporal dynamics are well represented by and correlated with remotely sensed and in situ oceanographic variables, indicating a high degree of spatio-temporal variability of each identified assemblage (e.g., warm-cold) in the NCC. This is the first study in this region that assesses Rrs555 in relation to the nekton and macro-invertebrate community. MODIS Rrs555 data are known to be related to pelagic biology in this region (Thomas and Weatherbee, 2006;Saldías et al., 2016), but this is the first study in this region to demonstrate the utility of using MODIS Rrs555 fields for biological studies. Furthermore, the species associated with the highest Rrs555 median values are generally understood to be related to the plume or relatively fresh waters, while those with the lowest Rrs555 median values given their distributions during this sampling period are indicators for the 'warm community' . In addition, we demonstrate the utility of using the interaction of both SST and SLA, as together they capture the local physical expression of ocean indices such as PDO and the North Pacific Gyre Oscillation in species distribution models, rather than individually or together as additive co-variates (as is often done in species distribution modeling of organisms).
SST in this region is at times strongly influenced by the warmer Columbia River plume waters (Saldías et al., 2016) and relatively cool upwelled waters during the spring-summer season. From our correlations of collocated MODIS SST and SLA at sampling locations, SST is positively correlated with SLA, indicating that these variables capture complementary yet distinct aspects of the marine environment (SST is only representative of the surface skin conditions of the ocean and SLA provides a more depth integrative understanding of the surface ocean capturing both density and regional circulation features).
The GAMM models predicted community differences based on RS oceanographic variables, capturing predicted NMDS1 scores outside the range observed during the sampled 10 years during all months. This is likely indicative of a distinct offshore community configuration occurring in the Oregon and Washington shelf and slope waters during the 2014-2015 (variable by season) period, due to anomalously warm water temperatures during 2014-2015 heatwave years (Gentemann et al., 2017;Peterson et al., 2017;Brodeur et al., 2019;Morgan et al., 2019). There are limitations with respect to the inferences we can make about the spatio-temporal community predictions past the extent of our empirical sampled data. That said, our results do suggest an anomalous offshore community in 2014-2015 which agrees with results from other studies in the region and adjacent that analyzed empirical data from that time period (e.g., Brodeur et al., 2019;Morgan et al., 2019). For example, Auth et al. (2018) recently identified novel ichthyoplankton assemblages in this region after 2014, due to the Northeastern Pacific (NEP) marine heatwave of 2014-2016 influencing the early and more coastal spawning of northern anchovy and Pacific sardine, and the more northerly and coastal spawning of Pacific hake (Merluccius productus). In the central California Current System, Santora et al. (2017) identified anomalously high diversity during this 2014-2016 NEP marine heatwave. Our prediction of community composition onto RS-data fields allowed for the spatial distribution of the compositional gradient, and thus has broader ecological implications. Almost all years and months had a significantly lower presence of the warm community along the Washington shelf than off Oregon (Figures 5, 6). This region is similar to that identified by Barceló et al. (2018) as a cold water refugia and provides evidence that the coastal region off southern Washington has a consistent onshore affinity community that is consistent across multiple months during the spring-summer upwelling season.
The two community clusters identified in this study, represent what we describe as the 'cold community' or the 'warm community.' These communities are represented by pelagic species that have either warm (offshore and/or southern) or cold-water (northern and/or nearshore) affinities (Brodeur et al., 2003). The 'warm community' might be further subdivided into a Columbia River plume vs. the warm, but non-plume community (Supplementary Figure 6). However, this separation was not distinguished in this study, and possibly not clearly defined for all months due to the variable spatial extent of the plume as well as the use of the plume by a wide range of both warm and cold affinity, nearshore and offshore species. Northern anchovy and Pacific sardine clustered together with sablefish in the ordination (Figure 2A) in a different region of the ordination than either Pacific saury or ocean sunfish, and are species known to be associated with the Columbia River Plume during the spring-summer upwelling season Litz et al., 2008). The indicator species in each month for the monthly warm clusters ( Table 1) consistently had Pacific saury as an indicator. This species is often captured at offshore stations; however, its schooling behavior leads to a highly variable catch predictability in these stations (Brodeur and Pearcy, 1986;. Sablefish (age-0) was also an indicator of the 'warm community' in May and June, but given median values of CTD salinity and Rrs555 it may also be a Columbia River plume associated species (Brodeur and Pearcy, 1986).
The ISA indicated that September also included Jack mackerel and ocean sunfish, likely due to an increased presence of warm offshore water along the shelf. With respect to the shared indicator species among the 3 months, May and June shared two of the species, and June and September shared two different species, while all 3 months shared Pacific saury. This is indicative of an expected, compositional continuum (vs. completely different warm communities between May to September) throughout the spring-summer season of makeup of the offshore community composition.
The fluctuation of the warm vs. cold community gradient has implications from a food-web and fisheries perspective. The closer to shore the preferred habitat for species of warm water affinity is, the higher the degree of competition and/or potential habitat compression for colder affinity species (such as juvenile and sub-adult salmonids). This concept is analogous to vertical habitat compression of the mesopelagic community due to the vertical expansion of oxygen minimum zones (Stewart et al., 2014;Checkley et al., 2017), but instead, from a horizontal perspective. A concentrated presence of the cold community into the nearshore implies potentially higher competition among similar trophic level species, leading to reduced success, as well as higher predation by coastal predators (colonial seabirds and pinnipeds), potentially further reducing already low populations of forage fish and salmon species. On the other hand, the near-shore intrusion and increased presence of the offshore community may introduce higher trophic level species (e.g., Santora et al., 2020), such as albacore, making these species more easily available to the fishery (less travel time to fishing grounds for fishers).
Further, community seascape derived indices can add to integrated ecosystem assessments, by filling in gaps in survey coverage (e.g., resulting from global pandemics) and providing information about the state of the ecosystem that can then be used to relate to other trophic levels and human activities. Similar to the cold-water lipid-rich copepod community versus warmwater, lipid-poor copepods in the northeastern Pacific Ocean (Peterson et al., 2017), the warm-cold community gradient metrics described here (Figures 5, 6) have both spatial and temporal qualities. As such, they measure both the 'flavor' of the regional community, as well as the spatial presence of different community groups. This metric of warm-cold epipelagic fish fluctuations could be added to the suite of regional biological condition indicators of the upper ocean that have been related to higher trophic levels, such as salmon returns (Burke et al., 2013), or as a predictor for the yearly distribution of albacore habitat for fishermen. This metric also could contribute to the suite of diversity indicators used in the California Current Integrated Ecosystem Assessments as a new, community composition-based indicator of upper ocean ecosystem state. Furthermore, these data and community predictions could be integrated into seasonal forecasting products, such as J-SCOPE 1 (Siedlecki et al., 2016;Malick et al., 2020) and extended to use Regional Ocean Model System (ROMS) data to provide useful information on coastal pelagic species' essential fish habitat.
In this study, we combined analytical techniques for understanding community structure with in situ and RS oceanographic observations to develop a better understanding of the community seascape. This analysis determined the habitat associations and community gradients in this upwelling and river-plume influenced pelagic marine ecosystem. These results demonstrated a statistically significant relationship between in situ and remotely sensed environmental variables and community composition that can have major benefits in monitoring and managing multi-species communities regionally, ranging from zooplankton to mammal and seabird communities. Daily satellite passes, and calculated satellite data composites (weekly, monthly) allow for continually updating community 1 www.nanoos.org/products/j-scope/home.php seascape maps that can be easily visualized by managers and stakeholders. The information gained by this approach allows for the identification of distinct communities and associated habitat characteristics for the whole community as well as distinct species and is directly applicable to ecosystem-based management as this metric could be linked to salmon returns data (Burke et al., 2013) and landings data of other species. This index, however, does not replace the important need for maintaining long time series of biological fishery independent surveys in the ocean so as to be able to "seatruth" (Barth et al., 2019) these remotely sensed derived indices.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These data are available through a data request to NOAA-NWFSC's Estuarine and Ocean Ecology Program [contact: David Huff (david.huff@noaa.gov) or Brian Burke (brian.burke@noaa.gov)].

ETHICS STATEMENT
All animal work was conducted according to relevant national guidelines. Fish were collected under the Endangered Species Act (ESA) Section 10 permit #1410-7A, which is the federal procedure for research directed by NOAA that includes ESAlisted species. Neither NOAA nor OSU collections of fishes required a separate review by an Institutional Animal Care and Use Committee.

AUTHOR CONTRIBUTIONS
CB with input from all authors, developed concept. RB and ED contributed biological and oceanographic data. LC provided base generalized additive model and gCV code. GS and CR contributed satellite data. CB conducted analysis and developed all figures and tables. All authors discussed the results and edited the manuscript.