Signals of resilience and change in tidepool fish communities on the Pacific coast of Vancouver Island, Canada

Ocean warming and marine heatwaves have emerged as strong predictors of biodiversity change in coastal systems. Here, we test for systematic change in tidepool fish communities between recent and historical data sets, which include the North Pacific marine heatwave of 2014–2016.


| INTRODUC TI ON
Shifts in the abundance and location of species are restructuring life on the Earth. Biodiversity change is being driven by a combination of immigration of warm species, loss of cool species and changing abundances as temperatures approach or move away from species' thermal optima that are tightly coupled to warming (Antão et al., 2020). A fundamental question is how species will respond to warming that already persists along steep environmental gradients.
Indeed, habitats that expose animals to relatively high temperatures and other extremes can select for species with resilient physiologies (Somero, 2010) that may allow individuals to resist or recover following disturbances (Hodgson et al., 2015;Holling, 1973).
Intertidal fishes residing in tidepools display remarkable heat tolerance, but also capacity to cope with changes in oxygen and pH (Horn et al., 1999). This group of fishes lives at the interface between terrestrial and marine environments, where they are exposed to both changes in ocean surface temperature, as well as heating during aerial exposure. The taxonomic diversity of these fishes, and the morphological, physiological and behavioural adaptations to survive a dynamic environment are of great interest. Indeed, tidepools may be sensitive and convenient indicators of environmental changes in our nearshore environments (Arakaki et al., 2014;Poloczanska et al., 2011).
The communities of fishes found within tidepools on rocky waveexposed shores along the Pacific coast of North America are diverse, including both resident and transient species (Allen et al., 2006;Gibson & Yoshiyama, 1999). The habitats of tidepool fishes in the north-eastern Pacific are undergoing warming, increases in pH and other changes resulting from human activity (Harley et al., 2006;Okey et al., 2014Okey et al., , 2015. Associated with these changes are reports of climate-related impacts on coastal organisms and communities (Harley & Paine, 2009;Menge et al., 2009Menge et al., , 2019Riche et al., 2014;Thom et al., 2014), and predictions of anticipated impacts (Harley et al., 2006;Okey et al., 2015).
Here, we take advantage of a resampling effort aiming to collect empirical evidence to identify change in species, which are little influenced by local and regional human activities. We sampled a series of identifiable tidepools encompassing most of the vertical extent of the intertidal zone for fishes on multiple occasions over 55 years (historical: from 1966 to 1993 and recent: from 2012 to 2016). Recent sampling spanned an extended marine heatwave documented across the North Pacific from 2014 to 2016 that leads to the formation of a warm temperature "Blob" and strong El Nino year in 2016 (Morgan et al., 2019). The event lead to widespread impacts on seabirds, fishes, coastal marine invertebrates and plankton in systems from California to Alaska (e.g. Laurel & Roger, 2019;Morgan et al., 2019;Peña et al., 2019;Piatt et al., 2020;Shanks et al., 2020).
Yet, tidepool fish communities can be resilient (Hodgson et al., 2015;Holling, 1973) by resisting or recovering following climate events (Almada & Faria, 2004). Thus, the long-term monitoring of intertidal fish communities has special significance at a time of changing ocean conditions in coastal waters.
Here, we ask whether the distribution and abundance patterns of tidepool fishes at a study site are distinctive in the years 2012-2016 in comparison with historical data. Specifically, we test for changes in the species density, abundance and community composition between the historical and recent sampling periods. We predict a decrease in the number of species and abundance of fish communities if recent warming has challenged intertidal fish communities. Given the magnitude and extent of the marine heatwave in the North Pacific, which occurred during our sampling programme, we expected signals of geographic range shifts-for warm-affinity species, abundance increases and increases in shore height, and the opposite for cool-affinity species.

| Study site
The study was carried out at a remote and rarely visited site approximately 12 km south of Bamfield on the west coast of Vancouver Island, British Columbia ( Figure 1). Here, a gravel beach extends for about 200 m between two cliffs, which prevent shore access to the beach except at low tide and under calm sea conditions. The site is not subject to, or close to, human activities and in 1970 was included in what became Canada's Pacific Rim National Park. The beach was once accessible by a path down the steep cliff separating it from the forested land above, but because of land slips, the path no longer exists and access to the shore is difficult. The rocky, wave-exposed shore has no safe landing site for recreational water crafts, and no fishing occurs near shore.
Shoreward the beach merges with dense coastal vegetation.
Seaward, a sandstone shelf with numerous tidepools extends up to 50 m offshore. Vertically, the shelf extends from below mean lower low water (MLLW) to several metres above mean high water (MHW). The intertidal zone is directly exposed to the Pacific Ocean, and has a fauna and flora typical of a wave-exposed rocky shore (Carefoot, 1977;Green, 1971;Ricketts & Calvin, 1966), including stands of the sea palm, Postelsia palmaeformis.  Tidepool area was determined from an aerial photograph taken by a drone positioned above a pool. Prior to taking a photograph, a one-metre reference rod was placed along the edge of the pool.

| Tidepool characteristics
Surface areas were later calculated using mapping tools in ArcView GIS software. The same software was used to produce maps of the study area from aerial photographs. Pool volumes were either determined by bailing water from a pool, and measuring the amount removed or estimated from their areas and average depths. The latter was estimated by measuring water depth at a minimum of four locations in the pool.

| Fish sampling
Fish sampling was conducted on 13 occasions between 1966 and 2016, whereby sampling can be grouped as eight historical events  and five recent sampling events (annually between 2012 and 2016). While most tidepools were sampled at each event, in 2016 only two pools were sampled (Table 1). Sampling captured all fish in each pool and consisted of dispersing rotenone into a pool when it was isolated at low tide followed by a thorough search of the pool during which all fish were recovered with small dip nets (Gibson, 1999;Green, 1971). This approach is unique in collecting rare and cryptic species (the full species list for the study period is reported in Table S1). Tidepools were sampled on the same low tide on each sampling date, led by the first author (JMG). On several occasions, the lowest pool (7) could not be sampled because of waves inundating the pool at low tide. Sampling in 1966 and 1971 was with a two-person team, while all other sampling events were conducted by a team of six or more field assistants. Fishes from different pools were preserved separately in 10% buffered formalin and later transferred to 45% isopropyl alcohol for identification.
We considered how our methods might have influenced our results. First, to understand the impact of repeated sampling, we also resampled six of the tidepools after ~2 weeks (1983) and ~4 weeks (1985). Short-term recruitment results are reported in the supporting material and indicate rapid recovery ( Figure S1). Second, these intertidal fishes all have a pelagic larval stage, but both the timing of spawning and timing of settlement of larvae into pools vary between species (Green per. obs.). The timing of our sampling was not always during the same time of the year (May, June and July); thus, we disregarded recruits from the year of sampling in our analyses.
Our results are therefore not influenced by years when sampling was done later in the summer, such as in 1993 when surveys were in July.

| Data analysis
Here, we calculated responses based on the fish assemblages sampled in each tidepool: species density (total number of species observed tidepool -1 ) and abundance (total number of individuals of all species observed tidepool -1 ). We also quantified compositional similarity between the first sampling event (considered the baseline) in each tidepool with all subsequent sampling events using the Bray-Curtis similarity (using the vegdist function in the package "vegan" with method=bray: Oksanen et al., 2017).
To test for change in each of the community data responses between the historical and more recent sampling periods (termed "Sampling Period") and with vertical shore height of each pool, we used a generalized linear mixed-effects model (GLMM, link= "log," family=Poisson) or linear mixed-effects model (LME; Pinheiro et al., 2016) modelling approach implemented with the package "nlme" (Pinheiro et al., 2016). The Poisson provided a good fit to the total species and abundance count data (which were integrated across each tidepool community and thus did not contain zero values). Model coefficients and 95% confidence intervals for each term are reported in Table S2. The reference was set to 0.42 m shore height (= lowest tidepool) and the historical period.
To test for shifts in the abundance of the most abundant species (the upper 75th percentile of the abundance distribution, where a threshold of species with more than a maximum count of 3 was identified) through time, we also used a GLMM (link = "log," using glmmTMB that allowed a negative binomial distribution for the species abundances where absences occurred (Brooks et al., 2017)).
Note the results are robust to the distribution assumption and modelling approach, and were similar to results returned from a general linear mixed-effects model fitted using MCMC: MCMCglmm

| RE SULTS
Overall, declines in species density, abundance and similarity occurred in the tidepools from the lowest vertical shore heights ( Figure 2 and Table S2: "Sampling Period" coefficient is negative and provides the estimate for the reference, which represents the shore height of the lowest pool (0.42 m) and the difference between the recent and historical sampling period). The higher pools, by contrast, did not display the same magnitude of decrease in species density, abundance and community similarity in the recent versus the historical sampling period (Figure 2 and Table S2: "Sampling Period: Shore Height" coefficient is positive and predicts the change in slope for every 1-m increase in shore height starting from the low shore height pool, which is the reference ("Sampling Period" as described in the previous sentence)). Figure S2 provides the raw data plots for each tidepool through time.
Of the most common species (n = 18 species, which reached maximum counts of 4 or greater), 12 displayed non-significant abundance changes across the study period ( Figure 3, Table S3). Eight rare species that were observed in historical sampling efforts were not detected in recent surveys (Figure 4). The only species that was gained at the study site was Gibbonsia montereyensis, a non-resident tidepool species. This species has a Pacific distribution from British Columbia to Mexico at depths to 20 metres, and was detected in 2015 and 2016 during the marine heatwave.
We further tested for shifts in the vertical distribution of the most abundant species surveyed (Table S4). We found that Phytichthys chirus (a tidepool specialist with a distribution in the eastern Pacific, from the Bering Sea coast of Alaska to southern California, USA) increased its vertical shore height in the more recent surveys (2012-2016), compared with the historical time period . Shore height distributions for the species observed in this study are displayed in Figure S3. of species, as well as replacements by more southern invasive and geographic range-shifting species in response to ocean warming (Harley et al., 2006;Okey et al., 2015). Yet, of the 27 fish species identified during this period, none were indicative of a more southern species extending its range into British Columbia, nor did those resident tidepool species with the most southern ranges exhibit systematic increases in abundance.

| D ISCUSS I ON
Our analyses also indicated that resident tidepool fish communities from higher shore locations were persistent and resilient through time, and during a marine heatwave, that persisted from 2014 to 2016. These results are consistent with earlier studies of north-eastern Pacific tidepool fishes that covered shorter timescales (Grossman, 1986;Yoshiyama et al., 1986), and of studies of rocky intertidal resident fish assemblages generally (Alameda

F I G U R E 2
In each tidepool, the number of fish species was counted (a. species density), total abundance of all fish present (b. total abundance), and Bray-Curtis similarity (c. community similarity for each year sampling compared with the first year each tidepool was surveyed). Shore height is the vertical distance pools were located from mean low low water (MLLW) (see Methods), and "Sampling Period" represents the historical (8 independent events from 1966 to 1993) compared with 5 recent sampling events in all years from 2012 to 2016. Patterns were analysed using generalized linear mixed-effects model (a & b) and linear mixed-effects model (c), as described in Methods with "Tidepool" included as a random effect. Coefficients from the model outputs for the fixed effects are plotted with significance that was supported when the 95% confidence limits did not cross zero (dotted grey line =0) and are also supported by p-values represented in the model results summary table (Table S2) F I G U R E 3 Total abundance counts per tidepool for the most common species recorded over the study duration (dots). Boxplots to the right of each time period represent the data spanning the 5th to the 95th percentiles. Blue shading indicates changes in abundance between the sampling periods that can be interpreted as significant (see Table S3: Artedius fenestralis, Liparis sp., Apodichthys fucorum and Anoplarchus purpurescens), and the arrow direction indicates a decline or increase (modelled using the package "glmmTMB" and corresponding function glmmTMB with family =nbinom1). Arrows were included for two borderline species (where the confidence window just crossed zero: Artedius lateralis and Oligocottus snyderi). The coefficient estimates and 95% confidence limits are reported in Table S3 & Indeed, A. lateralis is a predator on fishes (Boyle & Horn, 2006) and has frequently been collected with juvenile O. snyderi in their mouths F I G U R E 4 Proportion of tidepools in which relatively rare species were observed for each year in which sampling occurred over the study period. Rare is defined as the lower 25th percentile of the rank abundance counts across all species sampled (Gaston, 1994). The threshold maximum count according to this definition was three individuals (per tidepool). Eight of the nine species were detected in the historical sampling period , but not in recent years (shaded grey). Gibbonsia montereyensis was only observed at the study site in 2015 and 2016 during the marine heatwave [Colour figure can be viewed at wileyonlinelibrary.com] (Green per. obs.), but whether the recent trend in their abundances is driven by environmental factors or this predatory relationship is unknown.
Recruitment processes in tidepool cottids may also explain abundance variation from year to year and have been the subject of a number of studies along the north-east coast of North America (e.g. Pfister, 1996, Ritter, 2009, Ritter & Preisler, 2006, Shanks & Pfister, 2009, Webster et al., 2007. For instance, these studies have documented significant annual and regional variation in recruitment for the four most abundant cottid species in the current study, and both pre-and post-recruitment processes have been identified as determinants of the number of adults present in a pool. Thus, while the observed differences in the abundances of species from one sampling period to another were expected, what was surprising was the similarity in the numbers of some species between the historical and recent sampling periods. This was particularly true for some species, such as C. globiceps (Figure 3). In several years, populations of some sculpins were well above average and suggest that sampling followed a year of above-average larval recruitment for that species. The fact that spawning times vary between species and that the larvae of different species utilize different parts of the pelagic zone likely accounts for the fact that good recruitment conditions for one species may not correspond to that of another (Asch, 2015;Richardson et al., 1980).
The cottid genera found in the upper pools belong to a clade of related species which while likely evolving from a subtidal ancestor are thought to have speciated in the intertidal zone (Ramon & Knope, 2008) where they evolved a suite of adaptations to intertidal life (Buser et al., 2017;Knope, 2013), including wide temperature tolerances, capability for aerial respiration and homing behaviour.
Given these physiological adaptations to a harsh and high-energy environment, it is not surprising that they have not yet exhibited changes in either vertical or latitudinal distribution patterns near the centres of their geographic ranges.
Aside from the Family Cottidae, the other major component of the tidepool fish fauna was eel-shaped fishes of the families Stichaeidae (pricklebacks) and Pholidae (gunnels). Intertidal representatives of these families are also adapted to cope with wide variations in temperature and oxygen concentration (Allen et al., 2006;Horn et al., 1999). Their body shape, physiological adaptations for aerial respiration and broad temperature tolerances enable them to occupy both the tidepool and under-rock habitats of the littoral zone (Horn & Riegle, 1981). Although species from these families were always present in samples, just two species, Phytichthyes chirus and Xiphister atropurpureus, were consistently abundant in some pools. While the abundance of the former species did not change over the course of the study, it was the only species to show an increase in its vertical distribution. It is regularly present and the most abundant stichaeid in mid-to high tidepools at wave-exposed sites in British Columbia (Green, per. obs.). Both Apodichthyes flavidus (significant) and Xiphister atropurpureus (borderline significance) increased in abundance during the most recent collections, suggesting that favourable recruitment conditions are similar for both species.
Most of the stichaeids and pholids collected during the current study occur as far south as Baja California and are perhaps the most likely species to show increases in their populations with coastal warming.
These taxa are not as dependent upon the tidepool habitat specifically, and most are present in larger numbers in less wave-exposed intertidal sites where eel grass beds or under rock habitats are available (Green per. obs, Allen et al., 2006;Yoshiyama et al., 1986).
Numerous species of fishes with more southern distributions are now occurring in Canadian coastal waters, and there is concern for climate change effects on north-eastern Pacific commercial fish populations (e.g. Auth et al., 2018;Cheung & Frolicher, 2020). Even so, the higher intertidal tidepool fish fauna at our study location are not yet experiencing significant changes in distribution patterns or invasive species. The robust physiologies of the current tidepool inhabitants may position them to be among the last faunal component to see the changes that are currently being seen in other marine systems. Even so, our final year of sampling was in 2016, and lagged effects from the marine heatwave may have unfolded and compound with ongoing ocean warming in the region which has occurred since.
Continuing monitoring of these habitats will be important as significant changes may indicate of an ecological tipping point as these robust species are pushed beyond their physiological limits.

ACK N OWLED G EM ENTS
We thank the many individuals who assisted with field sampling.

CO N FLI C T O F I NTE R E S T
The authors do not have any potential sources of conflict of interest to disclose.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw species count data supporting the results are provided in the Supporting Information (Table S5)