Eelgrass Meadow Edge Habitat Heterogeneity Enhances Fish Diversity on the Pacific Coast of Canada

Eelgrass (Zostera marina) meadows are important fish habitats in temperate coastal areas. Understanding the relationships between seascape patterns—the spatial and temporal variability of biological and physiochemical drivers—and fish diversity in eelgrass meadows is crucial to conserving and managing these important habitats. The main objective of this study was to determine the environmental variables that influence the diversity of fish in eelgrass meadows in British Columbia, and whether a rich mosaic of edge habitats is positively associated with species richness and diversity, owing to the increased niche dimensionality and foraging opportunities provided by heterogeneous adjacent habitats. Using a spatiotemporal multispecies model based on long-term eelgrass fish diversity monitoring program data (2004–2020), we found that seascape variables, particularly those derived from unmanned aerial vehicles (meadow area, edge habitat heterogeneity), explained the most variation in species occurrence and abundance. We also found a positive effect of edge habitat heterogeneity on species richness in small and medium-sized meadows, with higher species richness and diversity in small and medium-sized meadows with high edge habitat heterogeneity. The relationship between edge habitat heterogeneity and species richness and diversity in large meadows was less clear. We also found that species richness has declined through time while diversity has been variable through time, remaining relatively stable in one region and generally decreasing in the other region. This analysis provides key insights into how seascape variables influence the distribution of species and the diversity of fish assemblages in nearshore eelgrass habitats in British Columbia.


Introduction
Eelgrass (Zostera marina) meadows form highly productive and valuable temperate coastal habitats, providing essential three-dimensional structure, foraging and growing opportunities for juvenile fishes, diverse faunal communities (Robinson et al. 2011;Siegle et al. 2014;Kennedy et al. 2018;Chalifour et al. 2019), and invertebrate prey resources (Huang et al. 2015;Stark et al. 2020). Eelgrass also provides important ecosystem services such as erosion control (Barbier et al. 2011) and carbon sequestration (Postlethwaite et al. 2018;Röhr et al. 2018;Prentice et al. 2019). Despite their ecological importance, eelgrass habitat continues to be threatened by a variety of anthropogenic and climate stressors that lead to habitat loss and fragmentation (Waycott et al. 2009;Iacarella et al. 2018;Howard et al. 2019). Understanding how characteristics of eelgrass habitats contribute to the diversity of fish in an estuary or coastal seascape and the ecological consequences of variability in the spatial structure of the seascape is critical for guiding and developing conservation and ecosystem-based management strategies.
Seascape ecology-the application of landscape ecology concepts to the marine environment-focuses on understanding the causes and consequences of spatial patterns in marine environments (Boström et al. 2011;Pittman et al. 2011;Pittman 2017) and provides a valuable lens for exploring the fish assemblages in eelgrass meadows. Seascape patterns are the spatial and temporal variability of biological and physical drivers that influence the distribution of resources (e.g., prey, shelter) and species distributions (Pittman 2017). Here, we focus on the aspects of the seascape that relate to the structure and orientation of the habitats in space. Understanding the linkages between seascape patterns and the ecological processes that influence the spatial distributions of species and biodiversity, as well as the functional implications of these patterns, can inform effective ecosystem-based management strategies and marine stewardship activities. Olson et al. (2019) used a seascape ecology approach to investigate seagrass meadow connectivity and found that seagrass nursery function for rockfish was strongly influenced by adjacent habitats and their associated subsidies. Olson et al. (2019) recommended that managing these important nearshore habitats should involve the consideration of the linkages among different habitats across the seascape.
The habitat heterogeneity hypothesis (MacArthur and MacArthur 1961) predicts that increased heterogeneity supports a greater number of species by providing a higher density of niches compared with more homogeneous habitats. Many studies have identified an important link between eelgrass edge and fish diversity (Macreadie et al. 2010;Yarnall et al. 2022), including an associated increase in fish species richness and abundance along the edges of eelgrass meadows compared to interior spaces (Smith et al. 2010). The distribution of resources along meadow edges plays an important role in how organisms distribute themselves within the meadow; if there is higher meadow edge habitat heterogeneity, species occupying the habitat will have access to a broader variety of resources, leading to higher species richness and abundance (Macreadie et al. 2010). Although previous studies allude to the positive influence of edge habitats on eelgrass fish diversity, most do not test the specific mechanism through which edge heterogeneity affects fish diversity (but see Smith et al. (2011), Macreadie et al. (2010), Carroll et al. (2019) (and references therein) for examples of studies that investigate mechanisms driving edge effects). Nevertheless, edge habitat heterogeneity is hypothesized to have the greatest impacts on small eelgrass meadows where more of the meadow is in close proximity to a habitat edge (Smith et al. 2010). Better assessing eelgrass habitat-edge heterogeneity and meadow size interactions in relation to fish diversity will ultimately better inform a wide range of habitat management decisions as well as conservation, restoration, and monitoring activities. Furthermore, understanding species-level trends, such as what species are driving the temporal variation in assemblage composition is important, as dominant species can have a strong influence on abundance and diversity estimates. These species-level trends provide important information in a conservation and management context and also provide information for future hypothesis testing, such as those related to predator-prey dynamics. Robinson et al. (2011) examined fish diversity in eelgrass meadows (Z. marina) in coastal British Columbia at an intra-regional spatial scale (1-10 km separation) and found that the location of meadows within an estuary was an important determinant of fish species richness variability. Furthermore, Robinson and Yakimishyn (2013) identified the importance of temperature and salinity gradients, as well as the distance between meadows on fish diversity, but a lack of detailed information about habitats surrounding the eelgrass meadows precluded a more fulsome assessment of the importance of edge habitat heterogeneity.
While more recent studies have focused on using underwater video or satellite imagery to map habitats surrounding eelgrass (e.g., Reshitnyk et al. 2014), recent developments in unmanned aerial vehicle (UAV) technology and geospatial mapping have presented novel opportunities to further spatially document how seascape variables such as meadow size and edge habitat heterogeneity might influence fish diversity in estuarine habitats like eelgrass (Nahirnick et al. 2019). Other tools capable of capturing high-resolution imagery of nearshore habitats and their associated seascape parameters, such as high-resolution satellite imagery and underwater videography, are also valuable for quantifying seascape patterns. Owing to their high resolution and ability to capture imagery at the scale relevant to our study, UAVs are valuable tools for identifying and quantifying the composition of habitats along the edges of eelgrass meadows in our study area. However, other tools, such as high-resolution satellite imagery, might be more appropriate for capturing spatial data across larger spatial extents (e.g., coastwide).
The main objective of this study was to assess how seascape metrics, particularly those derived from UAVs, influence the species richness, abundance, and diversity of fish in eelgrass meadows in Clayoquot and Barkley Sounds, British Columbia, Canada. An additional objective was to determine what species are driving the temporal and spatial variation in assemblage composition to better understand the dynamics of the eelgrass fish assemblage. To address these objectives, we have used a hierarchical joint species distribution modeling approach (Ovaskainen et al. 2017;Ovaskainen and Abrego 2020) to investigate fish species and assemblage patterns in response to the complex and dynamic spatial patterning that exists in the nearshore marine environment in Clayoquot and Barkley Sounds. A benefit of this approach is that we can assess and monitor changes in fish species diversity through time and along environmental gradients, and also investigate how individual species respond to the environment and what environmental variables are the most important drivers of species distributions. We specifically aim to address the questions: (1) What environmental 1 3 variables (seascape, physical) influence fish species richness and abundance in Clayoquot and Barkley Sound eelgrass meadows? (2) Does fish species richness and/or abundance increase as the degree of eelgrass meadow edge habitat heterogeneity increases, and how does this response vary with meadow size? In addition to questions and objectives related to the spatial drivers of eelgrass fish diversity, we are also interested in addressing the temporal question: How has species richness and abundance changed through time in Clayoquot and Barkley Sounds, and to what extent is this change driven by changes in temperature and salinity compared to unmeasured temporal variables?
We hypothesize that if the environmental variables related to the seascape (e.g., meadow size, edge type, edge habitat heterogeneity) are the primary drivers of fish diversity across space, then these variables will explain more variation in the fish assemblage composition than physical variables (i.e., temperature and salinity). Secondly, we hypothesize that if a rich mosaic of edge habitats provides more shelter, refuge, and niche dimensionality for fish species inhabiting eelgrass meadows, then species diversity (richness and abundance) should be positively associated with increases in edge habitat heterogeneity (Tews et al. 2004;Frid et al. 2018), but this relationship may be less pronounced as the amount of edge to interior habitat decreases. Thirdly, we hypothesize that if changes in temperature and salinity are the primary drivers of temporal changes in the fish community, then these variables will explain the majority of the observed temporal variation in richness and abundance.

Fish Assemblage Data
Data were collected through a long-term eelgrass fish diversity monitoring program conducted by Parks Canada Agency in Pacific Rim National Park Reserve (PRNPR), British Columbia, Canada (Yakimishyn et al. 2004;Fig. 1). A component of the monitoring program included fish surveys that were conducted at core eelgrass meadows on low morning tides in June and July in Clayoquot Sound (CS) and Barkley Sound (BS), respectively. Sampling at core eelgrass monitoring meadows (n = 20) occurred annually from 2004 to 2020. Surveys occurred within the 2-h window before and after the lowest early morning tide. At each meadow, the survey was comprised of triplicate beach seines conducted using a 9.2-m long beach seine with 4 mm stretch mesh and a 3.1-m drop in the center. Sets were positioned approximately 10 m apart along the deep edge of the meadow (Yakimishyn et al. 2004). After each seine, fish were removed and held in seawaterfilled rubber totes. After the third seine was completed, all fish collected from the three seines were identified to species, counted, and returned to the ocean. The abundance of each fish species in each meadow and year was estimated using catch-per-unit-effort (CPUE) and was calculated by dividing the total number of each species by the number of beach seine sets completed at each meadow: typically 3 sets and rarely 4 sets (n = 9 of the 314 surveys had four sets, with additional sets accounted for in the CPUE calculation). In situ sea surface temperature (C) and sea surface salinity (ppt) measurements were collected during each seine using a handheld salinity, conductivity, and temperature meter (YSI Model 30), and measurements were averaged to produce meadow-level measurements of salinity and temperature. Sea surface temperature and sea surface salinity measurements were taken between 50 cm and 1 m depths. Refer to Robinson et al. (2011) and Robinson and Yakimishyn (2013) for more details about the eelgrass fish sampling program. The protocol was reviewed by the Parks Canada animal care committee as part of the Parks Canada Research and Collection Permit process.
To ensure that fish species included in the analyses were well sampled, we excluded species that were present in less than 5% of surveys. The final data set used in the analysis was composed of 43 species (Table S1) caught in 314 surveys across 20 meadows between 2004 and 2020.

Multispecies Hierarchical Model
The spatiotemporal dynamics of the eelgrass fish assemblage in BS and CS were analyzed using the hierarchical modeling of species communities (HMSC) framework and package (Tikhonov et al. 2021) which is a model-based approach for analyzing assemblage data that produces inferences at both the species and assemblage levels (Ovaskainen et al. 2017). The framework uses Bayesian inference to simultaneously fit multivariate hierarchical generalized linear mixed models for each of the 43 species in our study. We utilized the approach to integrate data on species abundances, environmental covariates, and the spatio-temporal context (Ovaskainen et al. 2017). Our analysis utilized a hurdle model approach (Shelton et al. 2014;Thompson et al. 2022), consisting of two sub-models: a presence-absence model and an abundance model that was conditional on presence. For the conditional abundance model, the log CPUE of the species sampled at each meadow in each year was used as the response variable.
For each species, the presence-absence model described the presence of species j in each meadow i as: y ij ∼ Bernoulli Φ L ij , using a probit link function with a cumulative distribution Φ L ij . The log abundance (CPUE) of species j in meadow i was the response variable in the conditional abundance model and modeled as: y ij ∼ Normal L ij , 2 . The linear predictor ( L ij ) was derived from the fixed (F) and random (R) effects in the model: L ij = L F ij + L R ij . Further details and descriptions of model fit and assessment are provided in the Supplementary Information. For the presence-absence model, explanatory power was evaluated by computing Tjur's R 2 (the coefficient of discrimination; Tjur 2009) and AUC (area under the receiver operator characteristic curve; Pearce and Ferrier 2000) for each species. For the conditional abundance model, explanatory power was assessed using R 2 . The overall explanatory power of the presence-absence and conditional abundance models was calculated as the mean Tjur R 2 , mean AUC (presence-absence model), and mean R 2 (conditional abundance model) across all species. Predictive power was assessed using the same metrics and fivefold random cross-validation.

Environmental Fixed Effects
The fixed effects in both the presence-absence and conditional abundance models are listed in Table 1, and values are provided in Table S2. Seascape variables included the eelgrass meadow area, fetch (the distance over which wind-driven waves can build and a proxy for wind or wave exposure and delivery of small fish to the meadow), edge habitat information, and distance to the nearest estuary. Physical variables considered in this study that are known to have an influence on fish physiology and ultimately distribution included in situ sea surface temperature (°C) and sea surface salinity (parts per thousand) measurements collected at the time of sampling with a handheld YSI meter. It should be noted that while temperature and salinity are not considered "seascape" variables in this study, they (and other physical variables) are often described as seascape variables that represent spatial variability in the marine environment as a spatial gradient (Pittman 2017). The surface areal extent of each eelgrass meadow was quantified using UAV and towed video surveys conducted by Parks Canada and the Hakai Institute (Reshitnyk and Yakimishyn 2021). UAV surveys occurred in 2017 and 2018 in BS and CS, respectively. Eelgrass meadow surface area (m 2 ) was quantified from the UAV imagery by creating spatial polygons representing the meadows and calculating the metrics in a Geographic Information System (GIS). By including a time-invariant size for each meadow as a fixed effect, we are able to estimate how differences in size between meadows influence their composition, independent of any changes in size through time. Edge habitat information was quantified by creating a 3 m buffer around the eelgrass meadow using ArcMap (Esri Inc. 2019) and quantifying the proportion of the outer edge of the buffer that intersected each of the following 6 classes: vegetated rock, non-vegetated rock, vegetated sand, non-vegetated sand, vegetated cobble, and non-vegetated cobble. The 3 m buffer width was chosen based on the extent of the drone imagery that was captured for each meadow (i.e., how far the imagery extended beyond the edge of each eelgrass meadow). Attempts at delineating edge habitat within a 1 m buffer were also made; however, discerning edge habitats in the imagery within 1 m of the eelgrass meadow was difficult. Adjacent habitats derived from a 5 m buffer scenario were also created and found to be very similar to the adjacent habitats derived from the 3 m buffer scenario. As such, the 3 m buffer scenario was selected because it was within the extent of the drone imagery for all eelgrass meadows and was deemed to be a reasonable distance that would influence the sampled fish assemblage.
For each meadow, an edge habitat diversity metric was calculated using the proportion of the six edge habitat types quantified for each meadow. For each meadow, Hill-Shannon diversity, which is the exponent of Shannon diversity (Jost 2006), was calculated using the edge habitat information. The metric was calculated using the Vegan package (Oksanen et al. 2020) in R v.4.1.1 (R Core Team 2021). Several of the edge habitats were strongly correlated (Spearman's correlation coefficient (0.60-0.80)) with each other and were subsequently binned into two edge habitat classes. The "sand/vegetated cobble" class was comprised of vegetated sand, non-vegetated sand, and vegetated cobble, and the "rock/non-vegetated cobble" class was comprised of vegetated rock, non-vegetated rock, and non-vegetated cobble. These two classes were highly correlated with each other, and only the rock/non-vegetated cobble was retained in the analysis.
The distance of an eelgrass meadow to the nearest estuary provided a broader spatial perspective of the influence of temperature and salinity on fish diversity (Jelbart et al. 2007b), compared to physical measurements taken at the time of sampling. Estuaries are partially enclosed bodies of brackish water where one or more rivers or streams meet the ocean and are known to be important for a variety of marine species (Carr-Harris et al. 2015;Levings 2016). Estuaries function as an important transition area between a watershed's river network and the open ocean, providing an opportunity for anadromous fish to adjust to changes in salinity during migration (Clarke and Jamieson 2006;Carr-Harris et al. 2015;Levings 2016). Distance to the nearest estuary was calculated in ArcMap (Esri Inc. 2019) using the Euclidean Distance Tool and the location of Pacific Estuary Conservation Program (PECP) estuaries (PECP 2019) using land as a barrier feature. The PECP estuary classification system ranks estuaries based on size, species rarity, waterbird density, herring spawn index, and fish escapement (PECP 2019). The PECP estuaries include large, biologically important estuaries across the BC Coast and do not include the outflows of all small creeks and streams in an area. Of the ten estuaries located in Barkley Sound in the PECP dataset, two are in close proximity (< 10 km) from the nearest eelgrass monitoring sites. These include the Toquart River Estuary (9,884,902 m 2 ) and the Maggie River Estuary (422,557 m 2 ), which are both located on the northwest side of Barkley Sound. In Clayoquot Sound, the Cypre River Estuary, located on the north end of the study area, is less than 5 km from the nearest eelgrass monitoring location and is 763,143 m 2 . In both these regions, small streams and creeks that do not drain into estuaries mapped by the PECP are not included in the analysis. Efforts are currently underway to quantify the freshwater input from these smaller creeks and streams; however, this information was not available at the time this analysis was undertaken.
Fetch, a useful proxy for wind-wave exposure in the absence of a wind-wave model for the entire study area, is known to be an indicator of small-bodied fish distributions (Jordaan 2010;Robinson et al. 2011). Fetch is an important factor that influences the distribution of eelgrass and their associated fish assemblages (Grech and Coles 2010;Robinson et al. 2011;Oreska et al. 2021). Additionally, wave exposure affects the size, shape, and spatial configuration of eelgrass patches, with eelgrass meadows exposed to higher wave energy typically being smaller and more elongated (fringing) in shape (Frederiksen et al. 2004;Uhrin and Turner 2018). Fetch is also an important driver of the resuspension of sediment in shallow areas (Lawson et al. 2007). Higher wind and wave action can lead to greater sediment suspension and reduced light availability for primary production (Lawson et al. 2007). Fetch was calculated by Fields et al. (2020) for gridded point data at a 50 m resolution by generating fetch lines at bearing intervals of 5° (resulting in 72 fetch lines per point). The sum of the 72 fetch values was calculated for each point (Fields et al. 2020) to produce fetch estimates along the shoreline.
In situ temperature and salinity measurements collected at the time of sampling were included as environmental fixed effects. Linear temporal changes not associated with temporal changes in temperature and salinity were captured by the linear effect of the year that the monitoring survey occurred (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). Region (BS or CS) was included as a fixed effect. Environmental fixed effects were classified as either seascape variables or physical variables (Table 1). All fixed effects were included additively with the exception of meadow area and edge habitat heterogeneity, which were allowed to interact to permit the effect of edge habitat heterogeneity to vary with meadow area. This interaction was pertinent because the influence of meadow edge habitat heterogeneity on the fish assemblage may be less pronounced in large meadows where edge habitats may not be in close proximity to where fish sampling occurred.
During initial data exploration, seagrass biomass and shoot density variables were also included as environmental fixed effects. However, it was found that these variables did not improve the explanatory and predictive power of the model and were thus excluded from future iterations. Furthermore, these variables are captured in previous analyses that have investigated the influence of eelgrass structural complexity (as measured by eelgrass biomass and epiphyte load) on the fish assemblage in the study area (Robinson and Yakimishyn 2013;Iacarella et al. 2018).

Species Trait Data-Taxonomy
The HMSC framework permits the inclusion of species trait information to model how species' responses to environmental covariates depend on species traits and allows predictions to be made at the level of the trait. For this study, trait information was included to help with model fitting as opposed to testing a specific hypothesis related to 1 3 family-level responses to fixed effects or how species traits influence species assemblage processes as this was outside the scope of this study. By including trait information, model fit can be improved by allowing covariation within species sharing similar traits (Ovaskainen and Abrego 2020). For this analysis, we opted to use taxonomic information (family) for species traits, as phylogeny is often highly correlated with functionally important traits (Murillo et al. 2020).

Random Effects
We used random effects and a fixed effect for year to model variation in the eelgrass fish assemblage that is not associated with the environmental covariates in the model. A temporal random effect, specified by the year of collection, as well as the fixed effect for year was used to capture temporal variation in the eelgrass assemblage that is not associated with the temporally varying environmental covariates (temperature and salinity). Year was used as both a fixed and random effect because random effects are able to capture variation at smaller time scales (< 5 years), and the fixed effect is better at capturing long-term change (Thompson et al. 2022). An advantage of this approach is that the model can account for unmeasured temporal change in the seascape, such as change in meadow size or other unmeasured temporal change in seascape attributes. Although it is not possible to elucidate what variables are captured by year fixed effect and the temporal random effect, the model is still accounting for the change and so is able to more accurately estimate temporal changes in species occurrences and abundances. This approach is particularly valuable when seascape parameters derived from UAVs are utilized in models because the mapping technology has only recently become available, and temporally varying fixed effects may not be available.

Model Estimates
Species richness, diversity, and assemblage abundance were calculated for all 1000 posterior draws of the model, and median values were reported. For each species, we estimated the linear slope for occurrence, conditional abundance, and abundance across gradients of habitat heterogeneity and through time (i.e., change per year). This was repeated for each posterior draw to obtain the full posterior distribution of linear trends for each species and metric (occurrence, conditional abundance, abundance).

Species
The probability that an individual species was caught in a meadow was estimated by the presence-absence model. The abundance of species expected in a meadow was estimated by the conditional abundance model, if that species was present in the meadow (i.e., conditional on presence). A nonconditional estimate of fish abundance in eelgrass meadows was calculated by multiplying these two estimates together.

Assemblage
Species richness was calculated as the sum of the estimated occurrence probabilities from the presence-absence model (Ovaskainen and Abrego 2020;Thompson et al. 2022). Diversity was calculated as Hill-Shannon diversity (Jost 2006) using the estimated abundance (CPUE) for all species. Assemblage abundance is the sum of the estimated species abundance values and can be used to estimate the total number of individuals that would be present in a meadow. We estimated assemblage abundance by calculating the sum of the non-conditional estimated species abundance values (Ovaskainen and Abrego 2020; Thompson et al. 2022).

Variance Partitioning
To determine the degree to which the environmental variables drive the abundance and diversity of eelgrass fish species in CS and BS, we determined whether each species had a positive, negative, or neutral response to each environmental variable in the model. Significant responses were those where the estimated coefficients differed from zero with a greater than 95% probability. To determine which of the explanatory variables are most important in contributing to model predictions, we partitioned the variance explained by the fixed effects and random effects. The variance partitioning was calculated using the computeVariancePartitioning function in the Hmsc package (Tikhonov et al. 2021), which estimates the variance at the level of the linear predictor by summing the variance and covariances of variable contributions within groups of variables and ignoring the covariance among the groups of explanatory variables (Ovaskainen and Abrego 2020). We kept predictors ungrouped, except for the fixed and random effects of Year, which were grouped into a "Time" variable. Variance partitioning was done for each species, with values scaled by the total variation explained by the model (i.e., explanatory evaluation metrics Tjur R 2 and R 2 ). Assemblage level variation was quantified by calculating the mean variation explained by each variable across all species.

Predicting Across Environmental Gradients
Variation in species richness, diversity, and abundance across gradients of meadow area and edge habitat heterogeneitytwo variables identified as explaining the most variation in species composition in the variance partitioning plots-was quantified by using the model to make conditional estimates across gradients of edge habitat heterogeneity and meadow 1 3 area. This was done while holding all other variables at their most likely values given the edge habitat heterogeneity value, based on a linear relationship (Ovaskainen and Abrego 2020). Estimates were produced for BS and CS separately, and we used the year the UAV surveys were done in each region to make predictions (2017 in BS and 2018 in CS).

Predicting Across Time
Richness, abundance, and diversity estimates across all years in the monitoring program were produced using the model. All other variables were held at their most likely values given the year, based on a linear relationship (Ovaskainen and Abrego 2020), and predictions were produced for BS and CS separately. Predictions were also produced with and without Shiner Perch (a ubiquitous species caught in high numbers in 99% of surveys) to explore the influence that a single, dominant species had on diversity and abundance predictions. Estimates of change through time for each individual species were also produced, with the goal of understanding which species were potentially driving assemblage-level changes through time. Additionally, the credible intervals on these variables are provided by the analyses, which allowed for the assessment of whether the predicted temporal changes are significant (i.e., whether or not the credible interval overlapped zero).

Environmental Variables Driving the Variation in Eelgrass Fish Assemblage Composition in Clayoquot and Barkley Sounds
Seascape variables (eelgrass meadow area, distance to estuary, fetch, edge rock/non-vegetated cobble habitat, and edge habitat heterogeneity) together explained an average of 20.7% of variation in species occurrence and 22.6% of the variation in species conditional abundance ( Fig. 2; Table 2). Physical environmental variables (temperature and salinity) explained an average of 2.3% of variation in species occurrence and 2.6% of the variation in species conditional abundance ( Table 2). The random effect associated with the meadow explained an average of 4.1% and 4.9% of variation in species occurrence and conditional abundance, respectively (Table 2). For some species (e.g., Penpoint Gunnel, Copper Rockfish), the region (BS or CS) explained a substantial amount of the variation in species occurrence and conditional abundance. Species' responses to the fixed effects (Fig. S1) and model performance for all modeled species (Table S3) are illustrated in the Supplementary Information. Many species showed increases in occurrence probability in response to increases in fetch (wave exposure; Fig. S1) and decreases in occurrence probability through time.
Many species also showed differences in occurrence probability in response to region, suggesting that variability in the fish assemblage occurred even across relatively small geographic distances. In terms of conditional abundance responses to fixed effects, the patterns were less clear, with many species showing decreases in conditional abundance through time and variable responses to the remaining fixed effects. The variation explained for each species by each of the fixed and random effects in both models are provided in Table S4.

What Explains Temporal Variation in Assemblage Composition?
Variables that change temporally (temperature and salinity) together explained an average of 2.3% of the variation in species occurrence and 2.6% of the variation in species abundance. Time (fixed and random effects combined) explained 1.0% of the variation in occurrences and 2.2% of the variation in abundance. These results suggest that temperature and salinity are the primary drivers of temporal changes in the eelgrass fish assemblage, as more variation is attributed to temperature and salinity as opposed to any unmeasured temporally dynamic variables attributed with the fixed and random effects of year. Although the difference in variation explained by temperature and salinity vs. year is small, this variation is meaningful as it is responsible for the majority of the changes in richness and diversity through time. Figure S2 illustrates how the eelgrass fish community responds to gradients in salinity. In both regions, abundance is estimated to decline as salinity increases, while species diversity generally increases. In contrast, species richness is estimated to decline towards intermediate levels of salinity and then sharply increase at around 24-26 ppt. A similar trend in species richness and abundance is estimated along temperature gradients ( Fig. S3) with richness and abundance declining as temperatures increase to around 16° Celsius, and then sharply increasing. Many species respond positively to the temporally dynamic variables, including seaperch species (Shiner Perch, Kelp Perch, Pile Seaperch, Striped Seaperch) occurrence probability, which are positively associated with increases in temperature (Fig. S1). Crescent and Penpoint Gunnel occurrence probability are positively associated with increases in salinity (Fig. S1).

Impact of Meadow Size and Edge Habitat Heterogeneity on Diversity and Abundance
Species richness is estimated to increase as the degree of edge habitat heterogeneity increases in small and mediumsized meadows (Fig. 3a). Similarly, species diversity is estimated to increase slightly as the degree of edge habitat heterogeneity increases in small and medium-sized meadows in BS, except in the the smallest meadows, which are estimated to have the highest diversity at intermediate/high levels of edge habitat heterogeneity (Fig. 3c). The relationship between edge habitat heterogeneity and species diversity in the larger meadows in CS is less clear (Fig. 3f), with species diversity declining as the degree of edge habitat heterogeneity increases, except for the largest meadows which show a peak in species diversity at intermediate levels of edge habitat heterogeneity. The model predicts that fish abundance in both regions is relatively stable across different levels of edge habitat heterogeneity, except in large meadows, which show a steep decline in abundance as the degree of edge habitat heterogeneity increases. Additionally, a comparison of predicted and observed species richness is provided in Fig. S8 to illustrate that predicted richness values are conservative at the extremes, which occurs because the predicted species richness incorporates uncertainty about the sites (via the random effects).
The relationship between species abundance and edge habitat heterogeneity for the twelve most abundant species (excluding Shiner Perch) is illustrated in Fig. 4, which allows us to identify the species that are potentially driving changes in abundance and diversity at varying levels of edge habitat heterogeneity. The decline in Threespine Stickleback and Staghorn Sculpin abundance (species caught in high Variance explained is Tjur's R 2 for the presence-absence model and R 2 for the conditional abundance model. Species are arranged on the x-axis from highest to low-est variation explained in the presence-absence model. Variables that capture similar aspects of the environment colored similarly (i.e., seascape variables are in shades of pink/peach, physical variables are in shades of blue) 1 3 abundances) at high levels of edge habitat heterogeneity is likely contributing to the corresponding increases in diversity in BS and potentially the declines in abundance in the larger meadows. In CS, Threespine Stickleback also declines at high levels of edge habitat heterogeneity and is associated with declines in diversity because the species is not as dominant in CS as in BS, which leads to less even communities at higher levels of edge habitat heterogeneity. This influence may be what is driving the increase in diversity and declines in abundance in the larger meadows in CS. However, many of the common and moderately abundant species are also increasing as the degree of edge habitat increases (e.g., plainfin midshipman, black rockfish, and copper rockfish in BS and black rockfish, copper rockfish in CS) suggesting that they are also contributing to the increasing richness at high levels of edge habitat heterogeneity at some meadow sizes but perhaps also responsible for the decreases in diversity at the highest levels of edge habitat heterogeneity due the diverging abundances of many species.

Changes in Abundance and Diversity of the Assemblage Through Time
Model results indicate that species richness is declining through time (Fig. 5a). Hill-Shannon diversity is relatively stable through time, though on a declining trajectory. When Shiner Perch, a ubiquitous species (frequency of occurrence = 99%) caught in high numbers, is included in the diversity calculations, the model estimates a steep decline in Hill-Shannon diversity throughout the time series (Fig. 5c) and an increase in abundance beginning approximately in 2010 (Fig. 5b). This suggests that this decline in Hill-Shannon diversity is influenced by a substantial increase in Shiner Perch abundance through time. Figure 5 illustrates predictions for BS and CS, with predictions for BS showing similar trends to CS, but with a more pronounced decrease in Hill-Shannon diversity through time when Shiner Perch are excluded from the calculations.

What Species Are Driving the Temporal and Spatial Variation in Assemblage Composition?
A wide range of temporal trends were estimated for the 43 individual species in the model (Fig. S4). Generally, decreases in occurrence and abundance were more common than increases, with only seven species (Pile Seaperch, Blackeye Goby, Copper Rockfish, Plainfin Midshipman, Cabezon, Speckled Sandab, and Rockeweed Gunnel) estimated to increase in occurrence probability through time (Fig. S5). Twelve species are estimated to have decreased in occurrence probability through time.

3
The estimated decreases in species richness through time (Fig. 5a) are caused by the substantial number of species (28%) that are estimated to have decreased in occurrence probability through time, while only 16% were estimated to have increased (Fig. S5). The trend showing increasing assemblage abundance when Shiner Perch are included in the calculations (Fig. 5b) influences the estimates of change in diversity through time. When Shiner Perch are excluded, diversity is relatively stable through time (Fig. 5c) but on a downward trajectory in the later years of the survey. This suggests that this trajectory is likely not driven by a single dominant species, but rather caused by declines in the Shiner perch were removed from calculations because they are ubiquitous (frequency of occurrence = 99%) and caught in high numbers, which impacts diversity calculations abundance and occurrence probability of multiple species in the assemblage, leading to a less diverse assemblage with less even composition and fewer species.

Discussion
The primary objective of this study was to determine the environmental variables that influence the diversity of eelgrass fish in BS and CS. Furthermore, we were interested in the value of seascape metrics derived from UAV imagery, particularly those that describe the habitats at the eelgrass meadow edge, in understanding the spatial dynamics of eelgrass fish assemblages in our study area. We anticipated that seascape variables describing the eelgrass meadow and the nature of the surrounding seascape (i.e., meadow area, edge habitat, edge habitat heterogeneity, fetch, distance to estuaries) were the primary drivers of fish assemblage composition and variability across space. As hypothesized, our results indicate that seascape variables explain more variation in the fish assemblage composition than physical variables (i.e., temperature and salinity), time, and the unmeasured variability between meadows (random effect) and the two regions (fixed effect). Studies have shown that these and other seascape variables (such as exposure, mean wind velocity, type, and proximity of neighboring habitat) are important in explaining differences in benthic community (Turner et al. 1999) and fish assemblage (Gilby et al. 2018) composition. We also expected species diversity and abundance to be positively associated with increases in edge habitat heterogeneity because a rich mosaic of edge habitats would provide more shelter, refuge, and niche dimensionality for fish species inhabiting eelgrass meadows. While our results do show a clear pattern where the positive effect of edge habitat heterogeneity on species richness is strongest for small and medium-sized meadows (particularly in BS where meadows are generally smaller) the relationship between edge habitat and species richness and diversity in large meadows is less clear. This could be due to a mismatch between the spatial scale of the sampling event and the larger eelgrass meadows. For example, the edge habitat on the opposite side of a large meadow from where the fish sampling occurred may not influence the fish assemblage as much as edge habitat in closer proximity to the sampling event, such as in a smaller meadow. In smaller meadows, the majority of the meadow is in close proximity to the habitat edge and the structural habitat and resources associated with it, which likely leads to a stronger influence of edge habitat heterogeneity in smaller meadows. Jelbart et al. (2006) found consistently higher species richness in small versus large eelgrass meadows and also found that larger eelgrass meadows had lower species richness in edge regions compared to edge regions of smaller meadows. They postulated that smaller meadows may not be large enough to support high abundances of predatory species, which could lead to greater numbers of juvenile fish evading predation in smaller eelgrass meadows. While model results indicate no statistically significant responses to meadow size in the abundance of the following predatory species: Buffalo Sculpin, Lingcod, and Great Sculpin, species richness is generally estimated to be higher in small meadows (except in meadows with low edge habitat heterogeneity). Further investigation is required to determine how predator-prey interactions influence species richness and abundance across meadow sizes. Furthermore, movement patterns and the abundance of estuarine fish are known to be affected by the quality and proximity of edge habitats (Irlandi and Crawford 1997); however, this relationship is most often demonstrated in adjacent or nearby biogenic habitats (e.g., saltmarsh (Irlandi and Crawford 1997;Baillie et al. 2015), mangrove forests (Jelbart et al. 2007a;Skilleter et al. 2017), kelp (Olson et al. 2019), and reefs (Irlandi and Crawford 1997;Swadling et al. 2019)). Less attention has been paid to the role of non-biogenic edge habitats such as rock and cobble fields in terms of their impact on fish assemblages in adjacent eelgrass meadows. These habitats comprise a substantial portion of the temperate nearshore seascape, and it is likely that the composition of the edge habitat (e.g., vegetated rock vs. unvegetated sand) would attract and support different species with different life history traits, feeding modes, and preferences for eelgrass and other substrates. For example, the Plainfin Midshipman (Porichthys notatus) lays its eggs in cavities underneath cobbles and rock (Arora 1948), and the newly hatched young-of-the-year (YOY) are frequently caught in adjacent eelgrass meadows. Our model suggests that the occurrence of Plainfin Midshipman is positively associated with increases in edge habitat heterogeneity. Furthermore, Cowles et al. (2009) found that invertebrate prey abundance and productivity were variable across a range of coastal habitats and that productivity was highest in hard, structurally complex habitats (e.g., subtidal and intertidal rocky reefs). Smith et al. (2008) found that food availability is also variable along eelgrass patch edges and that fish utilized adjacent sand habitats extensively. They propose that adjacent sand habitat provides complementary habitat where food availability and foraging success are greater for eelgrass fish species, while still being in close proximity to the shelter that the eelgrass provides. Future research should focus on the mechanisms that connect the various habitats across the seascape (e.g., connectivity via predator-prey interactions, larval flow, active movement, nutrient flow).
Contrary to expectations, the relationship between meadow area, edge habitat heterogeneity, and fish abundance was less clear, with few clear trends estimated from the models. Fish abundance in small and mediumsized meadows increased slightly in response to edge habitat heterogeneity, whereas large meadows showed a decline in abundance as edge habitat heterogeneity increased. This decline appears to be driven by decreasing catches in Threespine Stickleback at high levels of edge habitat heterogeneity. This again illustrates the strong influence that a single species can have on the abundance (and subsequent diversity) calculations and the importance of unpacking the diversity trends to better understand what species are driving the variability across space and through time.
Many species showed a positive response to increases in fetch, and for some species, fetch explained a substantial amount of the variance in presence and/or abundance (e.g., Pacific snake prickleback, Pacific herring, kelp greenling; Table S4). Fetch and exposure have been shown to be important predictors of fish species composition at both local and regional scales (Jordaan 2010). Jordaan (2010) found that wave energy functions to create barriers to some species while simultaneously retaining other species within a specific location. They also found that local variability in wave energy results in different fish assemblages at the site (bay) level. Wave energy and exposure are known to influence the delivery and recruitment of small fish to meadows, as demonstrated by Markel et al. (2017), who showed that recruitment of Black Rockfish and a complex of Copper, Quillback, and Brown Rockfish was higher at sites with higher fetch.
The clearest signal of temporal change in the eelgrass fish assemblage is that species richness and diversity are declining through time. However, a major contributor to the decline in diversity is the increase in Shiner Perch abundance through time. When Shiner Perch are removed from the diversity, richness, and abundance calculations, diversity is estimated to have remained stable (CS) or be declining more gradually (BS). In both regions, species richness is estimated to be declining through time and can be further explored by investigating the temporal trends in occurrence probability (Fig. S5). More species are estimated to be decreasing than increasing, which is likely causing the decline in species richness through time. That being said, in both regions, the majority of species are estimated to be stable through time (shown in gray in Fig. S5), which suggests that despite the declining species richness through time, the eelgrass meadows in BS and CS continue to function as important nursery habitat for a substantial number of species. These results emphasize the importance of not only quantifying temporal changes in diversity, but also understanding how individual species are driving the changes. In this case, if the influence of a single dominant species (Shiner Perch) was not investigated, we might assume that species diversity is on a serious downward trajectory, when really this is due to the overwhelming abundances of Shiner Perch masking the underlying temporal trends in the fish assemblage.
We hypothesized that changes in temperature and salinity would explain more variation in fish diversity, abundance, and richness than time. Model results support this hypothesis; however, a relatively small amount of variance was attributed to temporally varying variables (temperature and salinity), which together explain an average of 2.3% and 2.6% of the variation in species occurrence and conditional abundance, respectively. This is a small amount of variance relative to the spatially-varying seascape variables, which suggests that the fish assemblages vary more across space than through time and that the spatial variability across the seascape is the important driver that shapes the fish assemblages in the study area. Fig. S2 illustrates how the eelgrass fish community responds to gradients in salinity. In both regions, abundance is estimated to decline as salinity increases, while species diversity generally increases, suggesting that declines in some species at higher salinities lead to a more even assemblage. In contrast, species richness is estimated to decline towards intermediate levels of salinity and then sharply increase at around 24-26 ppt. Future work could explore hypotheses related to these trends to elucidate the drivers of this variability. While salinity is known to influence fish community composition (e.g., Iacarella et al. (2018); Robinson et al. (2011)), the responses are variable, with some studies reporting increases in diversity at higher salinities (Nicolas et al. 2010) while others reporting decreases (Sosa-López et al. 2007). A similar trend in species richness and abundance is estimated along temperature gradients ( Fig. S3) with richness and abundance declining as temperatures increase to around 15° Celsius, and then sharply increasing. In CS, species diversity is estimated to decline as temperature increases, while in BS, species diversity is estimated to increase as temperatures increase until approximately 19 °C when it starts to decline. This decline in diversity is potentially due to increasing abundances of some species, such as shiner perch, at higher temperatures, leading to a less even fish assemblage. Highly mobile and abundant seaperches, such as shiner and pile perch, occupy a range of habitats in addition to eelgrass and are known to occupy warmer shallow-water habitats in the summer and move to deeper water in the winter (Lane and Secretariat 2002). As such, it is possible that these species are occupying a range of shallow habitats in addition to eelgrass in the summer months when temperatures are favorable. A relatively small amount of variance was attributed to the effect of Region (BS or CS), which explained 2.2% and 2.3% of variation in species occurrence and abundance, respectively. Some species are more common in one region compared to another. For example, Copper Rockfish and Black Rockfish are typically more common and abundant in CS than BS, while Threespine Stickleback is more common and abundant in BS. This variability could potentially be driven by unmeasured regional differences that are not captured by any of the fixed effects included in the model, such as differences in exposure (Barkley Sound Meadows are generally located adjacent to smaller, more exposed islands relative to Clayoquot Sound) and tidal current. However, these regional differences would be captured in our analysis by the random effect for meadow and the fixed effect for region. In general, in situ temperature and salinity measurements are typically higher in BS relative to CS; however, there is substantial overlap in temperature and salinity between the two regions. This variability could lead to differences in the timing of spawning events (e.g., leading to faster egg and larval development; Beyer et al. 2021) and productivity (Genner et al. 2010) that could be driving regional differences in species composition. Gradients in salinity and temperature are known to influence the physiological tolerances and in turn the distributions of different fish species (Martino and Able 2003;Harrison and Whitfield 2006;Robinson et al. 2011). Future research is necessary to understand the regional differences that may be driving the variation in fish assemblages in CS and BS.
Eelgrass is often studied as single habitat surrounded by a matrix of non-habitat (patch matrix model of seascape structure) (Joseph et al. 2006;MacArthur and Wilson 2016). However, many multi-habitat organisms, including mobile nearshore fish species, do not exhibit a binary response to habitat patch type but rather utilize the dynamic and heterogeneous habitat surrounding the eelgrass which in turn influences ecological function (e.g. species movement, foraging activities, predator-prey interaction) of the focal eelgrass nursery habitat. Analyzing eelgrass communities in the context of their surrounding habitats is becoming increasingly common (e.g., Olson et al. (2019); Gilby et al. (2018); Whippo et al. (2018); Swadling et al. (2019)) with results producing valuable insights into the ecology of eelgrass meadows and nearshore areas in general. This suggests that eelgrass meadows, including those in CS and BS, do not function in isolation, but rather are functionally connected to the broader mosaic of habitats across the seascape. Monitoring nearshore fish communities across a wider range of habitats and habitat gradients would be valuable for understanding the ecology of nearshore fish communities in terms of their habitat use and abundance outside of eelgrass meadows-which represent a small fraction of nearshore habitat.
Meadow area and edge habitat heterogeneity (two metrics derived from UAV imagery) were the most important predictors of assemblage composition, together explaining the majority of the variation in the model. These results highlight the value of UAV imagery for understanding the spatial dynamics of eelgrass fish communities. Because UAVs are capable of collecting high-resolution data across broad spatial areas, future work could focus on understanding the 1 3 broad-scale implications of the spatial arrangement of eelgrass and habitats beyond the immediate edge in terms of their influence on species utilizing eelgrass habitat. Many studies have investigated the influence of edge habitat type on eelgrass fish assemblages (Smith et al. 2008), but these studies often incorporate a single edge habitat type (except see Olson et al. (2019) and Gilby et al. (2018) for examples of including multiple adjacent habitats), and do not typically include maps or spatial representations of all edge habitats surrounding the seagrass meadow. UAVs could be used to address this gap to quantify the size of adjacent edge habitat patches and the spatial arrangement and shape of patches across the seascape to explore how this may impact species movement and resource use across the broader seascape, beyond the habitat edges. For this study, we only investigated the influence of edge habitats within 3 m of the eelgrass meadow, yet the response in species richness and diversity to edge habitat heterogeneity is still strong, which suggests that even looking at the habitat directly adjacent to the meadow can help us better understand how the composition and configuration of the broader seascape influence eelgrass fish assemblages.
Future work should also focus on the role of nearby kelp habitat in terms of its impact on fish diversity and abundance as kelp forests are highly productive ecosystems that provide important three-dimensional structural habitat for many species (Graham 2004;Teagle et al. 2017). For example, Olson et al. (2019) found YOY rockfish inhabiting eelgrass in close proximity to kelp forests consumed higher quality prey and had better body condition relative to YOY rockfish in eelgrass sites adjacent to unvegetated sand habitat. In the tropics, there is substantial evidence that coral reefs adjacent to seagrass and mangrove habitats support higher abundances and diversity of reef fish (Olds et al. 2012;Gilby et al. 2018;Berkström et al. 2020). While the evidence that connectivity between biogenic habitats (e.g. kelp and eelgrass) is well established, additional work is required to understand the functional relationships between eelgrass and other edge habitats, such as vegetated rock, sand flats, and cobble fields-habitats that can also provide important three-dimensional habitat, foraging opportunities and nesting sites for species occupying eelgrass meadows.
Because UAV data was collected only once throughout the time series, we account for unmeasured changes in meadow area and the quantity and composition of the edge habitats by including a fixed effect for year and a temporal random effect. Eelgrass meadows are known to be dynamic and fluctuate in size and shape (Ward et al. 1997;Frederiksen et al. 2004;Duarte et al. 2007), and in many areas, meadows are decreasing in coverage through time due to anthropogenic activities such as invasive species introductions (Howard et al. 2019) and coastal development (Nahirnick et al. 2020). In our study area, however, the meadows have remained of similar extent for the past 20 years of fish sampling (Robinson and Yakimishyn personal observations), and there are minimal anthropogenic impacts or upland disturbances because the majority of meadows are located within Pacific Rim National Park Reserve (PRNPR) or in areas with minimal anthropogenic activities. Of course, this does not protect the eelgrass from widespread impacts due to climate change that could be causing changes or declines in meadow area and the nature of the surrounding seascape. Future studies on fish diversity in eelgrass meadows as well as eelgrass monitoring in general would benefit from long-term UAV monitoring data collected annually at the core monitoring meadows in PRNPR. Other tools capable of capturing high-resolution nearshore habitat imagery (e.g., satellite imagery, seafloor video surveys) could also be used to investigate seascape patterns across broader spatial extents to support long-term monitoring studies. These tools could be particularly valuable to overcome some of the limitations of UAVs, such as suboptimal weather conditions and water clarity conditions that can decrease the optical depth of the water, making eelgrass delineation challenging (Nahirnick et al. 2019).
This analysis provides key insight into how seascape variables influence the distribution of species and the diversity of fish assemblages in nearshore eelgrass habitats in British Columbia. Our estimates of how individual species and the fish assemblage as a whole responds to seascape and physical variables can help optimize eelgrass restoration strategies to conserve biodiversity. Furthermore, by understanding the relationships between seascape patterns and fish diversity and productivity, we can better prioritize and manage these areas for conservation or restoration action. These results also highlight the value of small and medium-sized meadows for supporting diverse communities, which is valuable information for inclusion in conservation and monitoring activities that occur at a coast-wide scale. In particular, the results of this work are valuable for informing the monitoring and management of eelgrass within new and existing marine protected areas (MPAs), and for informing marine spatial planning exercises that are currently underway in BC, which are guided by the identification of ecologically and culturally important habitats. This information could also aid habitat managers in assessing project proposals that threaten, damage, or remove even small eelgrass meadows and help ensure that any mitigation projects that are potentially required also promote diversity by considering the habitat edges and the surrounding seascape. Effective ecosystembased management approaches are necessary for protecting eelgrass habitats and their associated fish communities. By considering the nature of the surrounding seascape, we can better understand, conserve, and manage these crucial habitats, support their nursery function, and minimize impacts to eelgrass meadows of all sizes. Data Availability Data sets generated and used in this study may be available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.