Multi-Scale Sampling to Evaluate Assemblage Dynamics in an Oceanic Marine Reserve

To resolve the capacity of Marine Protected Areas (MPA) to enhance fish productivity it is first necessary to understand how environmental conditions affect the distribution and abundance of fishes independent of potential reserve effects. Baseline fish production was examined from 2002–2004 through ichthyoplankton sampling in a large (10,878 km2) Southern Californian oceanic marine reserve, the Cowcod Conservation Area (CCA) that was established in 2001, and the Southern California Bight as a whole (238,000 km2 CalCOFI sampling domain). The CCA assemblage changed through time as the importance of oceanic-pelagic species decreased between 2002 (La Niña) and 2003 (El Niño) and then increased in 2004 (El Niño), while oceanic species and rockfishes displayed the opposite pattern. By contrast, the CalCOFI assemblage was relatively stable through time. Depth, temperature, and zooplankton explained more of the variability in assemblage structure at the CalCOFI scale than they did at the CCA scale. CalCOFI sampling revealed that oceanic species impinged upon the CCA between 2002 and 2003 in association with warmer offshore waters, thus explaining the increased influence of these species in the CCA during the El Nino years. Multi-scale, spatially explicit sampling and analysis was necessary to interpret assemblage dynamics in the CCA and likely will be needed to evaluate other focal oceanic marine reserves throughout the world.


Introduction
Marine protected areas (MPAs) potentially influence the dynamics of multiple species because protection from anthropogenic activity is extended, at least to some degree, to all species residing within a geographically bound region. A first step towards understanding the role of MPAs in an ecosystem context is to characterize baseline species assemblages that utilize a reserve and to evaluate how assemblages of species change through time.
An overarching objective of most MPA-based management plans is to augment regional fisheries productivity [1]. Given the diversity of habitats utilized by fishes within and around many MPAs, collecting quantitative, replicable data on multiple species is a major challenge for evaluating the dynamics of species assemblages. This constraint can be overcome to a large extent by collecting fishery independent data such as abundance of early life stages of fishes (i.e., ichthyoplankton). Because ichthyoplankton reflect fish spawning stock biomass [2,3] this method of sampling has the potential to assess directly whether MPAs impact fishery production from local to regional scales. Determining whether MPAs impact larval output can be complicated, however, because fluctuating environmental conditions are known to induce variability in fish spawning activity independent of reserve effects. To disentangle the impact of underlying environmental dynamics from the influence of MPAs on reproductive output it is necessary to first elucidate how environmental variability affects ichthyoplankton.
Oceanic MPAs are by definition nested within an open, advective system. Although the extent of MPA boundaries are static the spatial extent of water masses can change dramatically among years [4,5] introducing distinct assemblages of pelagic fishes [6] and zooplankton [7]. It is therefore likely that large-scale oceanographic processes will influence local assemblages within an oceanic MPA in a way that is not reflected solely at the local scale. Thus, a potentially fruitful approach for monitoring assemblage dynamics in an oceanic MPA is to couple detailed, focused sampling within the MPA with broader sampling to monitor the larger region encompassing the MPA.
The Cowcod Conservation Area (CCA; Fig. 1A) is the largest MPA in Southern California and the only one that includes extensive oceanic and nearshore habitats [8]. The CCA is embedded within the California Cooperative Oceanic Fisheries Investigations (CalCOFI) sampling domain where ichthyoplankton have been sampled regularly since 1951 at sixty-six stations including five within the CCA (Fig. 1B). Although analysis of CalCOFI data provided insight on ichthyoplankton dynamics at the scale of the SCB over decadal time periods (e.g., [9,10]) understanding of spatial and temporal patterns at finer scales within the CalCOFI domain is limited (but see [11,12]). To characterize baseline ichthyoplankton assemblages within the CCA, fine-scale surveys were conducted in the winters of the three years (2002)(2003)(2004) following establishment.
Here, we analyze the CCA data and accompanying CalCOFI data to determine how assemblage dynamics at the larger scale compared with the smaller CCA scale and discern how large scale dynamics affected the fish assemblage at the smaller scale. This study occurred during a transition from La Niñ a (2002) to El Niño (2003,2004) conditions [13] which enabled us to assess how fluctuating ocean conditions affected within-and between-year assemblage structure at both spatial scales.

Ethics Statement
Collections were made pursuant to a Memorandum of Understanding (MOU) by and between the California Depart-ment of Fish and Game, National Oceanic and Atmospheric Administration (NOAA)/National Marine Fisheries Service (NMFS)/Southwest Region, and NOAA/NMFS/Southwest Fisheries Science Center dated 11 February 1986 in a manner that ameliorated suffering of specimens.

Study area and sample design
The CCA is a restricted-fishing zone that was established in 2001 in response to population declines of the cowcod rockfish (Sebastes levis) [8]. It consists of an approximately 10,878 km 2 western (Figs. 1A, B) and 260 km 2 eastern area (Fig. 1B). Our analysis focused on the spatially-continuous western CCA (henceforth CCA; Fig. 1A). Within the CCA it is illegal to fish in waters deeper than 36 m. Although the CCA is entirely on the continental shelf, bottom depths are quite heterogeneous, ranging from sea level at the shores of San Nicholas and Santa Barbara Islands to near 2000 m in multiple basins.
Ichthyoplankton and oceanographic data were collected at two spatial scales: the focal CCA (Fig. 1A), and the larger, approximately 238,000 km 2 , California Cooperative Oceanic Fisheries Investigations (CalCOFI; see [14] for a description of the CalCOFI program) area that encompasses the CCA (Fig. 1B). At both scales sample sites were arranged in a grid with longitudinally separated lines running roughly perpendicular to shore. Sixty-six locations were sampled within the CCA each year with adjacent sample sites separated by approximately 9.5 and 18 km in longitudinal and latitudinal directions, respectively. Measurements were also taken from sixty-six locations at the CalCOFI scale. These sites constitute the ''core'' CalCOFI stations that have been sampled regularly since 1951 [15]. Adjacent CalCOFI lines were 74-km apart. Within lines, stations were unequally spaced as adjacent sample locations on the continental shelf were closer to one another (10 to 37 km) than seaward stations (74 km) (Fig. 1B).
The standard CalCOFI bongo net (71-cm diameter openings, 0.505 mm mesh nets, detachable 0.333 mm mesh cod ends; [11]) protocol was used to sample ichthyoplankton at all site [17]. A flowmeter attached to the net measured volume of filtered water. Samples were preserved in 5% formalin buffered with sodium borate, and all fish larvae were identified to the lowest possible taxonomic level in the laboratory. Most taxa were identified to species (Table S1) but, with the exception of Sebastes jordani, S. paucispinis, and S. levis, rockfish larvae are morphologically indistinguishable and were identified to genus (Sebastes spp.). The number of fish larvae under 10 m 2 surface area (to the maximum depth of a haul) was calculated for each taxon following the standard CalCOFI methodology [18].
Four environmental covariates were sampled at each site. Sea surface temperature (SST) and bottom depth (depth) were taken directly from ship instruments at the initiation of each haul. Small macrozooplankton displacement volume (henceforth ''zooplankton'') was calculated for each sample (ml of plankton displacement per 1000 m 3 of filtered water) [19]. We also estimated chlorophyll a values based on satellite imagery for the 30-day period encompassing each survey. This variable, however, was consistently highly correlated with at least one of the other covariates, and we therefore excluded chlorophyll a from further analysis.

Environmental Variability
We expected that the ichthyoplankton assemblage would be affected by oceanographic conditions that changed between the 2002 La Niñ a and 2003-2004 El Niñ o [13]. We determined if SST and zooplankton within our sample frames reflected oceanographic conditions expected during La Niñ a (cool SST, high zooplankton) and El Niñ o (high SST, low zooplankton) years in the Southern California region [15] and whether the sampling perspective affected the patterns.

Change in assemblage structure
To help visualize assemblage dynamics among years we first performed separate principle components analyses (PCA) on the site by species matrices from the CCA and CalCOFI scales and plotted site scores from the first two PC axes (PC1 and PC2). We then used redundancy analysis (RDA) to determine if there were differences in assemblage structure between years at the CCA and CalCOFI sample domains. Sample year was treated as a discrete variable and overall significance was assessed at a = 0.05 based on 1000 permutations of the data. If an overall difference was detected, comparisons between pairs of continuous years were made with significance levels adjusted to account for multiple comparisons (a = 0.017). Adjusted R 2 [20] values were also reported to quantify the unbiased coefficient of determination for each test. Prior to this and all other analyses on assemblages (see below) taxa were removed that were not found in at least 5% of the samples and stations were taken out that contained less than three individual larvae. One station was removed from the 2003 CalCOFI data set and nine from the 2004 CalCOFI analysis due to the low (0 or 1) number of captured larvae. Fish abundances were Hellinger-transformed prior to ordination because this transformation has been shown to produce unbiased results in Euclidian-based multivariate analysis such as PCA and RDA when zero species counts are prevalent [21].

Effect of covariates on interannual dynamics
We discerned the individual and combined effects of SST, zooplankton, depth, and year on CCA and CalCOFI assemblage dynamics. First, to evaluate the effects of the covariates on an entire assemblage, we conducted variance partitioning of scalespecific RDAs that modeled the covariates SST, zooplankton, depth and time against assemblage structure [20]. We included time as a distinct, categorical covariate to ascertain how much variance was explained by year of sampling that was not attributable to the sampled environmental parameters. Second, to better elucidate the dynamics of particular groups of species, we extracted the first two PC eigenvectors from the unconstrained, scale-specific PCAs that were based on the site by species matrices. We then used linear models to calculate how much of the variation in PC1 and PC2 was explained by the covariates. Finally, we used variance partitioning to ascertain the unique and shared contribution of each covariate to the explained variation of PC 1 and PC2 at each scale.

Effect of covariates on within-year distributions
We analyzed the horizontal distribution of assemblages within each year at both scales. We again evaluated how well covariates explained year-specific assemblage structure first using RDA and second linear models of PC1 and PC2 from PCAs of the site by species matrix from each assemblage at both scales. We then used variance partitioning to determine the amount of variation that was explained by each covariate. In addition to the environmental covariates (depth, SST and zooplankton) we included in these analyses spatial covariates. Spatial covariates provide insight towards processes affecting species distribution because they identify nonrandom distribution patterns that are not fully explained by the measured environmental covariates [21]. Such spatial effects indicate that unmeasured endogenous (i.e., behavioral or ecological) or exogenous (i.e., environmental parameters) processes affect an assemblage, and/or there is a mismatch between the scale of sampling and the response of species to environmental variables [22][23][24]. Thus, spatial covariates provide  Table 2  valuable clues regarding processes that affect assemblages that are not identified by the sampled environmental variables [24]. Spatial covariates were generated specifically for the CCA and CalCOFI domains using the principle coordinates of neighbor matrices (PCNM) method. PCNM produces eigenvectors based on a connectivity matrix (minimum spanning tree) of the geographical coordinates of sample sites (see Figure 1 in [25] for an illustration that describes the derivation of PCNM eigenvectors). The resultant spatial covariates (PCNM variables) depict spatial relationships among sample sites at a continuum of spatial scales within the bounds of the sample area [21]. PCNM variables with large eigenvalues describe broad-scale spatial relationships and vice versa.
We utilized a conservative approach for including PCNM variables in the analyses to avoid overfitting [26]. We reduced the number of PCNM variables by using a conservative forward selection technique based on adjusted R 2 values and a cutoff of p = 0.05 [27] and limited to five the number of PCNM variables in any one analysis [28].

Environmental dynamics
Interannual patterns of SST and zooplankton conformed to previous findings of low SST

Change in assemblage structure
The larval assemblage was much more dynamic at the smaller CCA than the larger CalCOFI scale over the three year period (Fig. 3). Within the CCA the overall RDA model of assemblage structure versus sample year was highly significant (R 2 = 0.  Table 1; Figs. S1, S2, S3, S4, 3A, 4; see [31] for the definition of species-specific habitat associations). From 2003 to 2004 the influence of the oceanic species declined while that of the coastal-oceanic species increased (Table 1; Figs. S1, S2, S3, S4, 3A, 4). Although the overall model at the CalCOFI scale was significant (p,0.01) the coefficient of determination (R 2 = 0.02)

Effect of covariates on interannual dynamics
The response of the assemblage to environmental fluctuation was much more predictable at the larger than the smaller scale. Variance partitioning of RDA models and linear models of individual PC axes from the PCAs of the original site by species matrices ( Table 2) showed that environmental parameters explained much more of the variation in assemblage structure at the CalCOFI than the CCA scale (Table 3). Further, in the CCA the amount of variation explained purely by time (i.e., time with other variables factored out) was greater than that characterized by the environmental covariates for both the RDA and linear model of PC1, which characterized a dichotomy between stations with benthic versus coastal-oceanic taxa ( Table 2, 3). At the CCA scale only SST was relatively effective in describing the variability of a linear model of PC2 (coastal-oceanic and benthic versus oceanic taxa; Tables 2, 3). At the CalCOFI scale, by contrast, pure time was unimportant for all analyses, and the environmental covariates were particularly effective in explaining the interannual variability in PC1 (72%), which separated oceanic versus coastaloceanic and benthic taxa (Tables 2, 3).

Effect of covariates on within-year distributions
Variance partitioning of models relating environmental and spatial covariates to assemblage structure at each scale within years revealed two main findings. First, the amount of variation explained by the environmental covariates (i.e., depth+SST+zooplankton) measured in the CCA was quite different among years, ranging from 7 to 22% for the RDA models and 14 to 45% for the linear models of PC1 (Table 4). At the CalCOFI scale, by contrast, the amount of variation explained by the environment was consistently higher as RDA models explained between 17 and 30% and linear models of PC1 between 66 and 79% of the variance in assemblage structure (Table 4). Second, the amount of variation explained purely by the spatial covariates (i.e., space with other covariates factored out) in the RDA models was consistently similar to the environmental covariates at the CCA scale. At the CalCOFI scale, however, spatial covariates explained much less of the variation than the environment for both RDA models and linear models of PC1 and PC2 (Table 4).

Discussion
The ichthyoplankton assemblage was much more dynamic within the CCA than the CalCOFI area over the three-year survey period that was characterized by a transition from La Niñ a to El Niñ o oceanographic conditions. During the La Niñ a in 2002 the CCA was dominated by species with coastal-oceanic habitat affinities (as described by [31]) such as E. mordax, M. productus, and L. stilbius as well as benthic rockfishes. During the 2003 El Niño, however, there was a decline in the influence of the coastal-oceanic species and an increase in rockfishes and several oceanic species (N. ritteri, D. atlanticus, S. californiensis). In 2004, E. mordax, L. stilbius and, especially, M. productus, increased (although not to 2002 levels). These results suggest that the coastal-oceanic species moved out of the CCA between 2002 and 2003 as the oceanic Table 2. Taxa loadings on the first two principle component axes from a PCA of the site by species matrix of samples from the CCA and CalCOFI sampling domains using data from all years.  species moved in, and then returned to the area in 2004. They also highlight how variable the ichthyoplankton assemblage can be in this region at an annual time scale. Previous ichthyoplankton research in this area characterized similar patterns but averaged over longer time frames. Hsieh et al. [9] examined the distribution and abundance of oceanic (e.g., N. ritteri, P. crockeri, L. ochotensis, D. atlanticus) and coastal-neritic (e.g., E. mordax, M. productus, S. sagax) species in association with a regime shift [4,5] from cool  to warm  conditions. They found that the abundance of oceanic species expanded and encroached shoreward during the warm period whereas the coastal-neritic species retreated shoreward and northward at this time. Similarly, ichthyoplankton abundances of species with subtropical affinities increased and moved shoreward throughout Table 3. Adjusted R 2 values from partitioning of the amount of variation explained by depth (D), sea surface temperature (SST), zooplankton volume (Z) and year of sampling (time) for interannual analysis of assemblage structure at the CCA and CalCOFI scales.  California and Baja California, Mexico, following a transition from cool, La Niñ a conditions in 1954-56 to warm, El Niño conditions in 1958-1959 [32]. Our results build on these analyses of longer time series by demonstrating how abruptly the ichthyoplankton assemblages can change over a relatively short time frame. Comparison of the relationship between covariates and assemblage structure at the CCA and CalCOFI scales revealed two main points. First, although the CCA assemblage changed through time, the sampled environmental variables explained relatively little of the variation. At the CalCOFI scale, however, environmental covariates better explained the dynamics of the assemblage. Second, purely temporal (year) and spatial (PCNM) covariates consistently explained more of the variation in assemblage structure than the measured environmental covariates at the CCA scale, but the opposite was true at the CalCOFI scale. These results may be caused by multiple nonexclusive factors [33]. One possibility is that the variation captured by year and unexplained spatial covariates within the CCA was induced by unmeasured environmental factors such as salinity [34], mixed layer depth [35] or dissolved oxygen concentration [10]. Another possibility is that the unexplained time and space variables were proxies for biological interactions. Biological processes such as schooling behavior, competition, and migration/dispersal can induce spatial clustering independent of environmental effects either by concentrating species within a small portion of suitable habitat or by causing species to occupy suboptimal habitat [33]. Because most of the larvae captured in our study were in the very early stages of development and were poor swimmers it is unlikely that behavioral processes had a large effect on horizontal distribution. However, if behavior influenced the distribution of spawning adults then biological interactions could be reflected in larval distributions.
A third possibility is that there was a mismatch between the scale of sampling within the CCA and the response of species to environmental variability. This explanation is supported by recent theoretical work indicating that if species respond to environmental fluctuation at a scale larger or smaller than a sample frame, then there is a tendency for the amount of variation that is explained by spatial covariates to increase relative to environmental parameters [23]. In our case, it is possible that at least the pelagic species responded to environmental variability at a scale that more closely matched the CalCOFI than the CCA sampling domain. For example, the range of SST across the CalCOFI region was 11.5-17.2uC during the study period but across the CCA it was only between 12.7 and 16.2uC. The relatively narrow range of temperatures encountered in the CCA likely had more subtle effects on fish distributions than did the large range of temperatures across the CalCOFI region. In particular, examination of the assemblage from the CalCOFI scale suggests that the assemblage was affected by shifting water mass boundaries as PC1 scores correlated highly with SST. CalCOFI PC1 separated stations characterized by taxa with oceanic affinities (e.g., N. ritteri, P. crockeri, and S. californiensis) from those with coastal-pelagic and benthic habitat association (e.g., Sebastses spp., E. mordax, S. leucopsarus and M. productus), and the distribution of these taxa shifted among years. In 2002 stations with high PC1 scores (i.e., oceanic species) clustered in the southwest portion of the sample domain. In addition, stations west of the shelf in the northern part of the CalCOFI domain contained coastal species, suggesting offshelf transport in 2002, a pattern also identified for S. sagax in spring of that year [13] Results from the CalCOFI scale demonstrate that ichthyoplankton assemblage structure is affected by the distribution of large-scale water masses that can have spatially dynamic boundaries [32]. Similar large-scale processes were documented in other parts of the world. For example, ichthyoplankton assemblages in the western Mediterranean were partitioned by two surface water masses, and the distribution of the water masses and assemblages varied seasonally [36]. Similarly, the spatial structure of the ichthyoplankton assemblage throughout the equatorial Pacific varied in response to extended periods of warmth over a 13-year period [37]. In addition, ichthyoplankton assemblages off the coast of Oregon varied significantly in response to warming and cooling of surface waters between 1999 and 2006 [38,39]. Together, these findings demonstrate that although ichthyoplankton assemblages are often spatially segregated, the locations where particular assemblages are found can vary dramatically through time in response to environmental forcing [40][41][42].
Although our results suggest that broad-scale processes affected the CCA assemblage in a way that was not fully apparent by examining local conditions, there is evidence that local processes strongly impacted the CCA assemblage in 2002 as environmental covariates explained more than twice as much variation in this year than in other years. This result might be driven by interannual variability in mesoscale oceanographic structure. The Southern California Bight is bathymetrically heterogeneous and contains several deep basins with steep vertical drops that can induce the formation of local fronts (Fig. 1). Nishimoto and Washburn [43], for example, documented a mesoscale eddy just north of our study area that affected significantly the distribution of late-stage fish larvae in 1998 but found no evidence of the eddy or spatially-structured fish distributions the following year. Working within the CCA region, McClatchie et al. [12] documented a strong north-south front that hugged the ridge south from San Nicholas Island (Fig. 1A) in 2010. This front separated warm, saline water in the east and cool, fresher water in the west, and there were different zooplankton, larval fish, and egg (fish and squid) assemblages on either side of and within the front. In our study there was a clear east-west gradient in both species distributions and SST within the CCA in 2002 but less distinct transitions in other years (Fig. 4, 5). These results suggest that local oceanographic structure contributed more to assemblage structure within the CCA in 2002 than other years and imply that both local and broad-scale oceanographic dynamics can combine to affect the distribution of ichthyoplankton within the CCA.
These results stress the importance of being cognizant of sampling scale when describing species-environment relationships [44]. Although ichthyoplankton studies from around the world documented differences in assemblage structure from very small (e.g., ,1-m; [45]) to very large spatial scales (e.g, 4000-km; [6]), explicit evaluations of assemblage-environment relationships across sampling scales within a study are rare. Among the existing multi-scale investigations, Catalan et al. [46] quantified ichthyoplankton assemblage structure from a mesoscale sampling grid (stations separated by 18 km) that was embedded within a macroscale grid (stations separated by 40 km) off the coast of Spain and Portugal. Similar to our results they found consistent differences in assemblage composition between onshelf and offshelf stations at the macroscale. In addition, they discovered through mesoscale sampling that the spatial boundary separating assemblages was porous as taxa with primarily onshelf affinities were advected to offshore areas. Another multi-scale study from the Gulf of Mexico demonstrated that different oceanographic variables best explained ichthyoplankton distribution patterns at fine (1-km), meso (10-km) and coarse (100-km) scale sampling perspectives [47]. These results stress how examining ichthyoplankton from multiple perspectives is important given the inherently dynamic nature of fish distribution.
Ultimately, the aim of future/ongoing research is to compare icthyoplankton assemblage data reported here with future samples to ascertain whether the MPA has an effect on the local production of fishes. Given the findings of this initial effort several recommendations can be made to augment the efficacy of future sampling in this and other ichthyoplankton assemblage studies. First, we found that the CCA assemblage was highly variable at an annual time scale. However, our results are based on only three years of sampling and additional annual sampling is necessary to evaluate how these fluctuations compare to long-term trends. Second, we showed that assemblages were likely affected by unmeasured environmental covariates, endogenous behavior by the species, and/or a mismatch between the scale of sampling and the response of species to the environment. To better resolve the relative importance of these factors future sampling needs to quantify additional oceanographic covariates such as salinity, current velocity and oxygen in concert with species sampling. Third, we obtained insight into processes that affected the CCA assemblage by also sampling at the larger CalCOFI scale. This result emphasizes the importance of recognizing that oceanic MPAs are nested in a broader system where forces external to the reserve can affect assemblage structure within the reserve. Hence, our results stress the need for multi-scale monitoring to elucidate the cause of species fluctuation in MPAs. Our findings set the stage for documentation of the effect of MPAs on fisheries production, which is a critical, yet largely unresolved question in the study of MPAs.