Seasonal Dynamics of Faunal Diversity and Population Ecology in an Estuarine Seagrass Bed

Biodiversity is important for communities to be resilient to a changing world, but patterns of diversity fluctuate naturally over time. Understanding these shifts — and the species driving community dynamics — is crucial for informing future ecological research and conservation management plans. We investigated the impacts of seasonality, small-scale changes in seagrass cover, and small-scale spatial location on the epifaunal communities occupying a temperate seagrass bed in the South Island of New Zealand. By sampling epifaunal communities using a fine-mesh push net two to three times per season for 1 year, and using a combination of multivariate and hierarchical diversity analyses, we discovered that season, seagrass cover, and the location within the bay, and their interactions, explained 88.5% of the variation in community composition. Community composition and abundances, but not numbers, of species changed over seasons. The most common taxa were commercially important Caridean shrimp and juvenile flounder (Rhombosolea spp.), and both decreased in abundance in summer (shrimp: 1.40/m2 in winter to 0.80/m2 in summer; flounder: 0.15/m2 in winter to 0.01/m2 in summer). Other commercially important species were captured as juveniles, including blue cod (Parapercis colias), kahawai (Arripis trutta), and whitebait (Galaxias spp.). The only adult fish captured in the study were two pipefish species (Stigmatopora nigra and Leptonotus elevatus), which had distinctly seasonal breeding patterns, with reproductively active adults most likely to be found in the spring and fall. Our study highlights the importance of estimating biodiversity parameters based on sampling throughout the year, as some species will be overlooked. We demonstrate that the temperate estuarine seagrass-affiliated animal communities differ in response to season and fine-scale local environments, causing fluctuations in biodiversity throughout the year.


Introduction
Foundation species, such as seagrasses, mangroves, and coral reefs, provide ecosystem functions and services such as nutrient cycling, sediment trapping and stabilisation, attenuating water flow, and carbon sequestration (Costanza et al. 1997;Heck et al. 2003;Morrison et al. 2014;Nordlund et al. 2016;Ruiz-Frau et al. 2017;Orth et al. 2020). Moreover, these foundation species also facilitate biodiversity by providing habitat and supporting diverse food webs (Spalding et al. 2003;Gillanders 2006;Guannel et al. 2016). Uncovering the habitat components facilitating biodiversity is important for understanding the ecological impacts of predicted habitat loss (Waycott et al. 2009).
Estuaries with seagrass beds are known to harbour high biodiversity, in part because seagrass provides refugia to a large number of organisms (Orth and van Montfrans 1987;Orth et al. 2006;Grech et al. 2012), especially juvenile fish that rely on seagrass as a nursery habitat (Hemminga and Duarte 2000;Heck et al. 2003;Gillanders 2006;Shibuno et al. 2008;Bertelli and Unsworth 2014). Some species that use seagrass as a nursery include species important for fisheries industries, such as shrimp and prawns (Watson et al. 1993;Haywood et al. 1995;Schaffmeister et al. 2006;Taylor et al. 2017), flounder (Earl et al. 2014;Hyndes et al. 2018), and Australasian snapper . By providing shelter for prey species, seagrass has the ability to alter predator-prey relationships while simultaneously supporting species in other trophic levels such as sea turtles, Communicated by Ronald Baker dugongs, tiger sharks, and marine birds (Aragones and Marsh 1999;Beck et al. 2001;Heithaus et al. 2002;Turner and Schwarz 2006;Bertelli and Unsworth 2014). These nursery habitats might provide an influx of certain species during specific seasons, for example when juvenile fish settle out from larvae or when larger predators use nurseries as feeding grounds, thereby affecting seasonal patterns of community diversity.
Seagrass biomass can also fluctuate seasonally, impacting the species that rely on it as a habitat. Seagrass characteristics such as above-ground biomass, shoot density or height, and epiphyte biomass can positively correlate with the densities of certain species and overall community diversity (Orth et al. 1984;Ray et al. 2014). In temperate regions, seagrass traits such as density are impacted by seasonal changes in water temperature, light, day-length, salinity, nutrient concentrations, and dominant phytoplankton species (Hemminga and Duarte 2000;Jankowska et al. 2014). Temperate seagrass meadows often follow a cyclic pattern in response to these abiotic variables (Duarte 1989;Olesen and Sand-Jensen 1994), and faunal communities in seagrass beds likewise experience seasonal fluctuations (Bauer 1985;Edgar and Shaw 1995;Nakaoka et al. 2001;Ribeiro et al. 2006;Hutchinson et al. 2014;Jankowska et al. 2014;Włodarska-Kowalczuk et al. 2014). Fluctuations in seagrass and its associated communities also respond to more episodic or local stressors like algal blooms Carstensen et al. 2015), marine heat waves (Oliver et al. 2018Smale et al. 2019), and extreme weather events (Harris et al., 2020), which can be seasonally linked, but might also vary at fine spatial scales -for example, due to freshwater inputs carrying nutrients into coastal habitats; (Lapointe et al. 2015). These unpredictable environmental perturbations might counteract and complicate predictable cues in shaping estuarine biodiversity.
Furthermore, biodiversity of estuarine seagrass beds can also be impacted by seascape features such as fragmentation and loss of connectivity between biogenic habitats (Boström et al. 2011), for example created by freshwater inflows or large unvegetated areas (Boström et al. 2006), or by anthropogenic stressors like propeller scars, haul seining, and trampling (Orth et al. 2017;Uhrin and Holmquist 2003). For species that use seagrass as a nursery, proximity to habitats used as adults (e.g., rocky reefs) is an important predictor of species abundance (Rees et al. 2018). The spatial scale of the impacts of fragmentation between habitats differs based on the type of organism, and has been measured across a wide variety of scales (Boström et al. 2011). Biodiversity across interconnected seascapes, such as seagrass and mangrove habitats, can also fluctuate across seasons (Tarimo et al. 2022), as fauna transition between nursery grounds or change in response to biogenic habitat seasonality (Cuvillier et al. 2017).
Here, we quantify seasonal changes to epifaunal community diversity in a temperate seagrass bed in New Zealand. In this bioregion, only one species of seagrass occurs, Zostera muelleri (Turner and Schwarz 2006;Short et al. 2007), which has seasonal but rare sexual reproduction in New Zealand (Ramage and Schiel 1998;Dos Santos and Matheson 2017). Seagrass in New Zealand primarily occurs in intertidal mudflats of sheltered estuaries but has also been found on rocky reef platforms, coastal beaches, and in subtidal waters (Inglis 2003;Turner and Schwarz 2006;Morrison et al. 2014). Seagrass is declining throughout New Zealand due to anthropogenic impacts (Turner and Schwarz 2006;Matheson et al. 2011;Matheson and Wadhwa 2012;Morrison et al. 2014). In New Zealand, seagrass habitats have been suggested to be biodiversity hotspots, but seasonal data on epifaunal community structures supporting this notion remain sparse (Turner and Schwarz 2006;Morrison et al. 2014). More specifically, research on seasonal dynamics of epifaunal seagrass communities are lacking in terms of understanding small-scale habitat characteristics shaping biodiversity and have largely focussed on sampling in one season (e.g., van Houte-Howes et al. 2004;Morrison et al. 2014). Here, we address this research gap by testing for impacts of small-scale spatial variation, seagrass cover, and season on epifaunal community diversity in a New Zealand seagrass bed. We analyse population dynamics of three focal taxa for which we collected size or age data, providing a unique contribution to understanding the drivers of community diversity in an estuarine seagrass bed.

Methods
This study was conducted in Duvauchelle Bay (43.7504° S, 172.9334° E), located in the Akaroa Harbour on the South Island of New Zealand (Fig. 1). Duvauchelle Bay is one of the northernmost bays in Akaroa Harbour, over 1 km from the harbour mouth. The bay is approximately 800 m wide with mosaic patches of no seagrass (mud and sand flats) and sparse, medium ("patchy"), or consistent dense cover within close proximity to each other (Figs. 1 and S1). We assessed temporal variation in animal communities occupying the seagrass habitat in Duvauchelle Bay by sampling monthly from November 2019 through October 2020 (see Table S1 for detail on sampling dates and effort). We refer to the sampled animal communities as epifaunal communities, using a broad definition that encompasses invertebrates and fish found on or among the seagrass leaves (Chen et al. 2021). The full list of animals we include is listed in Table S2, as recommended by Chen et al. (2021).

Monthly Sampling Methods
Each sampling event began approximately 2 h before low tide in the intertidal zone and sampled the shallow subtidal zone as the tide receded. To sample, we used a handhauled 1.2 m × 1.2-m push net (5-mm mesh) mounted on a PVC frame, dragging the net parallel to shore in a series of non-overlapping tows through 0.5-0.75-m-deep water. This method is ideal for capturing small animals with minimal impact on the seagrass beds (Thomsen et al. 2020). Because tows were a non-uniform length, we recorded time-stamped coordinates for the start and end of each tow using a Garmin GPS (average tow length was 30.6 m). At the end of each tow, the density of the seagrass was recorded based on a visual assessment from one of four semi-quantitative categories: bare patches, in which no seagrass was observed; sparse seagrass, in which some seagrass shoots were seen but the majority of the substrate was bare; patchy seagrass, when the tow included a mix of bare areas and seagrass, but more seagrass occurred than in sparse tows; or dense seagrass, when seagrass was the primary substrate visible (see Fig. S1 for representative drone images). Tows were conducted to remain within a single seagrass cover category as much as possible. All animals captured by the net in each tow were identified to the lowest taxonomic level (species for some, family or infraorder for others), and the counts of each type were recorded. Some taxa, including all fish, were recorded in size categories to capture ontogenetic variation throughout the year.
Pipefish species (Stigmatopora nigra and Leptonotus elevatus) were also identified by sex and classified as either adults or juveniles. Adult males of both species can be recognised by the presence of a brood pouch. We identified adult females based on body size and the presence of female ornamentation. All pipefish larger than 60 mm were tagged using visible implant fluorescent elastomer tags (Northwest Marine Technologies). Using established protocols for tagging syngnathid fishes for minimally invasive population monitoring (e.g., Masonjones et al. 2010;Flanagan et al. 2017;Masonjones and Rose 2019), pipefish were anesthetised using 100 mg/L of clove oil in 500 mL of seawater and three coloured tags in unique combinations (to facilitate individual identification) were injected just below the skin. Fish recovered in buckets of fresh seawater and were gently released back into the subtidal seagrass habitats. These

Analysis of Community Diversity
All analyses were conducted in R v. 4.1.0 (R Core Team 2021), and plots were created using packages magick (Ooms 2021), leaflet (Cheng et al. 2021), and mapview (Appelhans et al. 2021). Each sampling date was lumped by southern hemisphere season (Summer = Dec-Feb; Fall = Mar-May; Winter = June-Aug; Spring = Sep-Nov). We used the GPS coordinates from the start and stop of each tow to calculate the tow distances using the R package geosphere (Hijmans 2019) and multiplied the linear distances by the width of the net (1.2 m) to calculate the tow area. A single tow was removed from the analyses because the seagrass cover was not recorded.
To analyse effects associated with fine-scale spatial variation, we used distance-based Moran's Eigenvector Mapping (dbMEM), a method which accounts for spatial autocorrelation in the analysis of community composition (Dray et al. 2006;Legendre and Legendre 2012). This approach is valuable because it accounts for interacting features of the bay impacting community composition such as seagrass cover, tidal zone, and unmeasured spatial features of the environment. This method uses eigenvalue mapping of the tow coordinates to reduce spatial data to a smaller number of axes that describe spatially autocorrelated community variation, which can then be used as predictor variables in a partial redundancy analysis of the community composition. The final model of relative densities of taxa included the effects of seagrass cover, season, the most informative axes describing small-scale spatial variation, and the overlap of each of those components. A full description of the analysis can be found in the Supplementary Methods. We used R packages adespatial (Dray et al. 2021), SoDA (Chambers 2020), vegan 2.5-7 (Oksanen et al. 2015).
A hierarchical analysis of patterns of diversity is most appropriate, given that our study design was hierarchical, with seagrass cover nested within sampling date, which itself is nested within season. The hierDiversity package (Marion et al. 2015a, b) is designed for such sampling designs, and estimates alpha, beta, and gamma diversity in a hierarchical framework for specified Hill numbers (Hill 1973;Jost 2007;Chao et al. 2014). Hill numbers are a unified set of diversity metrics, each associated with an order q, which down-weight rare species as q increases. When q = 0, all species carry equal weight (corresponding to traditional "species richness"); at q = 1, the weighted geometric mean is used, and the diversity measure is a transformed version of Shannon's entropy; for q = 2, the weighted arithmetic mean of species is used, and abundant species have more influence on the diversity metric (Hill 1973;Jost 2007;Chao et al. 2014). Calculating multiple diversity estimates using multiple Hill numbers (i.e., creating "diversity profile plots") allowed us to inspect the importance of dominant versus less abundant species in shaping diversity and turnover (Marion et al. 2015b(Marion et al. , 2021. For this analysis, we calculated pairwise turnover (scripts available on GitHub: https:// github. com/ spfla nagan/ ecolo gy-duvau chelle/ tree/ main/R/ hierD ivers ity) to avoid biasing our beta diversity estimates due to uneven sampling across hierarchical levels (Marion et al. 2017). See Supplemental Methods for further details. Rarefaction analysis, done in vegan (Oksanen et al. 2015), indicated that some of our tows did not sufficiently capture the diversity of the community. Therefore, we created 999 rarefied communities, in which we down-sampled to the number of animals caught in the shortest tow, and ran our modified hierDiversity analysis on each of those communities. Using the distributions of these diversity estimates from the rarefied communities, we assessed whether the diversity profiles from our observed community were impacted by sampling effects. This analysis used the densities (the number of individuals divided by the tow area, to account for variation among tow areas) of all observations except those with zero animals captured in the tow.
To test for associations between particular taxonomic groups and seagrass cover and season, we used the multilevel pattern analysis in the R package indicspecies (Cáceres and Legendre 2009). We ran this analysis independently for seagrass cover and for season using the function multipatt to identify the "core" community membership associated with each test factor. The results of these analyses were compared to the species loadings from the distance-based Moran's eigenvector maps analysis to identify species with strong associations with these two variables.

Population Dynamics of Focal Taxa
We further investigated the impacts of season and seagrass cover on the abundance of key taxa in Caridean shrimp, pātiki (flounder, Rhombosolea spp.), and pipefish (Stigmatopora nigra and Leptonotus elevatus). We focused on shrimp and pātiki because they were the two most abundant and frequently captured taxa in the dataset (see "Results"), and on pipefish because they were part of a larger study. Furthermore, we also measured sizes for these taxa, allowing us to test for seasonal ontogenetic changes. For each taxon, we fitted generalised linear models with the most appropriate link functions and used a model testing approach to identify the model with the lowest Akaike information criterion, corrected for small sample sizes (AICc).
Shrimp densities were normally distributed after log transformation and we therefore used a linear model (instead of a generalised linear model). The predictor variables included seagrass cover, season, and size class of the shrimp (≤ 50 mm or > 50 mm), and their interactions. Residual plots from the best-fitting model were evaluated for model fit and assumptions.
Flounder densities were not normally distributed after log-transformation nor were their residuals from linear models, and the raw counts were overdispersed (the mean counts, 1.59, were much smaller than the variance, 14.29). Data with these properties are best modelled using a zeroinflated negative binomial regression, here implemented in the pcsl package v. 1.5.5 (Zeileis et al. 2008;Jackman 2020). We modelled raw counts and included the tow area as a predictor. As with the shrimp data, we modelled counts predicted by seagrass cover, season, and size class (≤ 10 mm, 10-20 mm, 20-50 mm, > 50 mm). All four size classes are considered juveniles as we did not catch any fish > 15 cm (Earl et al. 2014). To assess the fit of the model, we visualised the residuals compared to each variable in the model to ensure residuals were centred around 0 and calculated the dispersion statistic as the Σ(residuals 2 )∕(N − p − 1) , where N is the number of individuals and p is the number of parameters (including interaction coefficients). Furthermore, we created a confusion matrix of the 0 and non-zero observations to estimate the false negative rate. Finally, we used ordinary least squares regression to investigate the relationship between the observed non-zero counts and the non-zero counts predicted by the model. Specifically, we extracted the predicted counts and used bootstrapping implemented by the boot package v. 1.3-28 (Canty and Ripley 2020) to estimate confidence intervals.
Because pipefish were rare compared to flounder and shrimp (see results), and therefore their counts overdispersed and densities highly skewed, we modelled their presence, instead of abundance, using a multinomial logistic regression. This approach enabled us to predict the probability of encountering reproductively mature pipefish, given the season, seagrass cover, and pipefish species identity (Stigmatopora nigra vs Leptonotus elevatus). Photographs and size data were used to identify individuals as adults or juveniles and assign them to the appropriate species, and count data to match the individuals with the specific tows in which they were captured. Tows without pipefish were given the status "absent" instead of "juvenile" or "adult." We then modelled the status of individual pipefish, predicted by seagrass cover, season, their interaction, and species identity using a multinomial log-linear regression implemented in nnet v. 7.3-16 (Ripley and Venables 2021). To ensure assumptions of the model were not violated, we visually inspected the residuals compared to each variable in the model. Finally, we used the effects package version 4.2-0 (Fox et al. 2020) to generate predicted probabilities and confidence intervals for each level of status across the predictor variables.

Results
A total of 403 tows covering 14,812 m 2 were sampled over 12 dates, with a total of 133 tows in summer (3 dates), 52 in fall (2 dates), 108 in winter (3 dates), and 110 in spring (4 dates). We caught 24,511 animals, represented by 29 taxa, including 10 fish, 7 crabs, and five snails and whelks (see Table S2 for total counts and mean densities). The most abundant taxa were Caridean shrimp (67.72% of the total catch; all shrimp were identified to infraorder), followed by Rhombosolea spp. (10.65% of total catch) and Diloma subrostrata (5.75% of total catch). Species accumulation curves suggest that sampling effort was insufficient for most species except Rhombosolea spp. and Caridean shrimp (Fig. S2).

Community Composition and Diversity
Partitioning the variation in the seagrass communities revealed that the spatial position of the tows, seagrass cover, season, and their overlapping effects shaped the abundances of the animals in the seagrass beds, with each factor having a significant effect in partial redundancy analyses (P ≤ 0.001; Table S3). The distance-based redundancy analysis with the interacting effects of all three types of explanatory variables (seagrass cover, season, and space) was significant (F 317,85 = 2.063, p = 0.001). The constrained axes explained 88.5% of the variation in community composition, with the first two axes accounting for ≥ 65% of variation (64.51% of total and 72.91% of constrained; Fig. 2A and B). Three of the 25 retained axes describing the spatial location of the tows had substantial loadings in the redundancy analysis, and those axes primarily corresponded to specific regions of the bay, matching bay-wide east-west and diagonal intertidal off-shore gradients ( Fig. S3; Supplementary Results). Summer samples were most differentiated from other seasons ( Fig. 2A) and were most similar to the ellipse describing samples from dense seagrass habitats ( Fig. 2A and B). The species with the largest magnitude loadings on the first axis were Caridean shrimp, Rhombosolea spp., Austrovenus stutchburyi, and Potamopyrgus antipodarum, and on the second axis Diloma subrostrata, A. stutchburyi, P. antipodarum, Rhombosolea spp., and Caridean shrimp (Fig. 2B). Flounder were associated with variation in a direction orthogonal to dense habitats (i.e., they were least associated with dense seagrass), although the four seagrass cover classes had substantial overlap in the first two multivariate axes (Fig. 2B).
The alpha diversity profile plots suggest that, across hierarchical levels, community composition was dominated by a few abundant species, with most species occurring in low numbers. Both the observed (Fig. 3, top row) and the rarefied (Fig. 3, bottom row) communities showed that alpha diversity decreased with increasing q (i.e., as rare-species are down-weighted), with values for richness (q = 0) around 5 effective taxa per tow for each seagrass cover category and 20 for each season, both decreasing to 2 to 5 effective taxa per tow for all categories at q = 2 (the inverse of Simpson's index). These patterns suggest that the communities are dominated by several highly abundant species, with most other species occurring at low abundances. However, we found no major differences in alpha diversity across seagrass cover or seasons, although estimates from the rarefied communities suggest that more species were present in the spring than in the other seasons (Fig. 3).
The pairwise turnover plots suggest that within each level of seagrass cover, there were few disparities between common versus rare species (see flat profile plots in Fig. 4),   Fig. 2 The first two constrained coordination analysis axes incorporating the effects of season, seagrass cover, and spatial covariates from a distance-based Moran's eigenvalue matrix. A and B show the same axes (explaining 65% of the constrained variation) and sample points, but where overlain ordiellipses representing seasonal patterns (A) and seagrass cover (B), respectively. Also plotted are eigenvectors, scaled by eigenvalue, for the five taxa with the highest loadings  F). These values reflect the effective number of taxa captured per tow. At q = 0, all species have the same weighting, but as q increases, rare species have less influence on the metric. The negative slopes on the profile plots indicate that common species drive patterns of spe-cies richness. A-C Observed alpha diversity values with bootstrapestimated confidence intervals, D-F Estimates of alpha diversity and associated confidence intervals estimated from 999 permutation communities rarefied down to the smallest tow distance in the dataset (note that some confidence thresholds are so narrow that they are not visible in the plots) suggesting that samples share dominant and low-abundance species within each seagrass habitat. These analyses suggest that sparse habitats had lower estimates of pairwise turnover than other seagrass types. However, at higher hierarchical levels (season and overall), communities differed in less abundant species but shared the same highly abundant species, as can be inferred by the diversity profile plots' negative slopes (Marion et al. 2021). Again, seasons did not differ in the observed estimates, although the rarefied estimates suggested small differences, with larger disparity in communities between tows in fall compared to other seasons (Fig. 4).
The tests of abundance patterns across seasons revealed that most taxa were associated with combinations of nonsummer seasons (concordant with the multivariate analysis, Fig. 2). Specifically, seven taxa were associated with the winter-spring-fall combination, including   season (B, E), nested in the overall pairwise turnover (C, F). At q = 0, all species have the same weighting, but as q increases, rare species have less influence on the metric. The negative slopes on the profile plots indicate that common species drive patterns of sample dissimilarities. A-C Observed pairwise turnover values with bootstrap-estimated confidence intervals. D-F Estimates of pairwise turnover and associated confidence intervals estimated from 999 permutation communities rarefied down to the smallest tow distance in the dataset (note that some confidence thresholds are so narrow that they are not visible in the plots) Furthermore, C. glandiformis was associated with the winter-spring combination (indicator value = 0.372, p = 0.005), juvenile Labridae with winter-fall (indicator value = 0.382, p = 0.005), Pyura pachydermatina with spring-summer (indicator value = 0.274, p = 0.02), and L. elevatus with spring-fall (indicator value = 0.255, p = 0.03; Table S4).
The two pipefish species had largely similar responses to seagrass cover, season, and their interaction (Table S7). Both species were found in low abundances, and our model reflects this by predicting the probability of absences in a given tow, in any habitat or season, as > 0.8. Pipefish were least likely to be observed in bare habitats, regardless of other variables, and most likely to be observed in patchy or dense seagrass habitats (Fig. 7). Both species of all ages had the lowest predicted probabilities in those preferred habitats during the summer (predicted probabilities of being captured in a given tow < 0.05) and the highest probabilities in the spring and fall for both age categories (Fig. 7). In the winter, the two age classes differed, with adults having near-zero capture probabilities but juveniles showing no substantial decline from fall to spring (Fig. 7). None of the 81 marked individuals were recaptured, so we are unable to estimate population sizes based on recapture rates.

Discussion
We found that seagrass cover, season, and small-scale spatial location explain community diversity and structure in Duvauchelle Bay, New Zealand. The communities were defined by the most abundant species (in particular, shrimp and flounder), and across seasons, the communities differed in their least abundant species but not the highly abundant species. Despite the importance of seasonality in shaping the community composition (i.e., which rare species were present), seasons only had minor effects on universal community diversity metrics. Our combination of multivariate statistical analyses and hierarchical diversity partitioning allowed us to maximise the insights gained from monthly sampling data.
Fine-scale spatial factors shaped community composition in Duvauchelle Bay. The three spatial axes substantially contributing to multivariate community diversity axes were associated with four distinct regions of the bay that would not have been easily described by other categorical variables such as intertidal region or the side of the bay (Fig.  S3). Instead, these axes appear to correspond to fine-scale seascape features related to specific locations in the bay (Fig.  S3). We hypothesise that these fine-scale features are related to a natural break in connectivity between the east and west sides of the bay caused by a freshwater inflow, varying high tide depths and exposures at low tide, and proximity to other habitat types including a rocky mussel reef on the edges of the bay. In this way, our data-driven approach enabled us to identify and include in our analysis seascape features that would have otherwise been overlooked. Further research is needed to conduct fine-scale mapping of these seascape features in Duvauchelle Bay and throughout Akaroa Harbour to enable incorporating these specific features into the analysis of community diversity (as done in Henderson et al. 2017).
While the mechanisms contributing to the importance of these spatial variables are unknown, fine-scale effects of seagrass traits have been shown to impact predation rates on small fishes (Chacin and Stallings 2016) and fish community diversity (Henderson et al. 2017), and factors such as wave exposure can impact fish assemblages (Perry et al. 2017). In line with our results emphasising the interacting effects of fine spatial scale variation and season, shallowwater fish communities in Zanzibar were impacted by offshore variables (e.g., wave exposure) in different ways across the seasons (Perry et al. 2017). Our study was not designed to detect broad-scale seascape features, nor were we able to incorporate seagrass traits in our analyses, but these are fruitful areas of future study as fish diversity has been shown to respond to both fine-and broad-scale habitat features (Berkström et al. 2013) and proximity to seascape features such as the open ocean (Henderson et al. 2017). Furthermore, the connectivity to the nearby habitats in other bays (including seagrass meadows, rocky reefs, and kelp forests) is likely to be an important factor impacting the biodiversity we report here (Cuvillier et al. 2017;Henderson et al. 2017;Rees et al. 2018;Tarimo et al. 2022).
The effect of seasons on the animal communities wasin part -driven by declines in the abundance of the most common animals (flounder and shrimp) in the summer. These declines were counteracted by increases in abundance of rare species, as reflected in the alpha diversity profile plots (summer having the largest alpha diversity at q = 1 and q = 2). This summer decline was surprising, in part because these taxa have previously been shown to peak in summer in other locations (Webb 1972;Bauer 1985), and also because fish are expected to be most common in New Zealand seagrasses in summer (Morrison et al. 2014). We suggest that the observed decreases can be attributed to macroalgal blooms observed in late January. Macroalgal blooms typically have negative impacts on seagrasses (Höffle et al. 2012), particularly in summer months, as well as on invertebrate communities in sedimentary habitats (Lyons et al. 2012). However, in lower abundances, macroalgae can also provide complex habitats and refugia for some species (Lyons et al. 2012;Thomsen and Wernberg 2015), which could explain the moderate increase in alpha diversity in the summer months. Another possibility is that the macroalgae reduced our ability to capture different species, in particular the more mobile flounder and shrimp.
Our result that less abundant species had highest turnover across seasons reflects patterns in fish recruitment. Most juvenile fish were captured only during one season (Table S1), which likely reflects species-specific variation in spawning habits. For example, smelt and blue cod juveniles were only found in winter, and kahawai mainly in summer. These timings correspond to breeding seasons as blue cod breed in winter and spring in Otago (Robertson 1973) and kahawai in Canterbury breed in summer (Webb 1973). Only typical seagrass residents such as pipefish, blenny, and small wrasses (Kendrick and Hyndes 2003;Gillanders 2006) and juvenile flounder were found year-round. All flounder captured in this study were small young-of-the-year (all < 100 mm, which is less than 2-year-old flounder (Webb 1972)), suggesting that Duvauchelle Bay serves as an important nursery ground, especially given the abundance of food (e.g., shrimp) and habitat complexity. However, whether this pattern can be directly attributed to the seagrass is uncertain, because juvenile flounder did not show a strong preference for patchy or dense seagrass and are often also found on mudflats (Webb 1972;Earl et al. 2014;Thomsen et al. 2020). Additionally, our sampling method targets small fishes across small spatial scales, so the occupancy of the seagrass by adults is unknown. Experimental studies are needed to identify the underpinning factors that cause seagrass beds in New Zealand to be valuable nursery habitats for juvenile fish. Understanding these factors is also of cultural and economic importance, as flounder are taonga (a treasure) to the local hapū (sub-tribe) as mahinga kai (traditional food source).
The interaction of seagrass cover and season was an important predictor in every model we analysed, likely because temperate seagrass cover typically fluctuates throughout the year Fig. 5 Mean shrimp abundances (± standard errors) were determined by season, seagrass cover, and their size. Two size categories are represented by different colours and shapes, and in all seasons except fall, smaller shrimp were more abundant. A-D show the predicted densities (ln(density + 1)) and E-H show the observed number of shrimp in bare (A, E), sparse (B, F), patchy (C, G), and dense (D, H) seagrass covers (e.g., Jankowska et al. 2014). As in other systems , the magnitude of the differences along our bare-to-dense gradient of seagrass cover changed throughout the season for our focal species (Figs. 5, 6, and 7). Although spatial variability in seagrass density at small scales has been demonstrated to be important for diversity of epifaunal invertebrates (Gullström et al. 2012), the seasonal impacts of this small-scale variability are unclear. Furthermore, seagrass biomass and community diversity could have been impacted by (unmeasured) fine-scale features such as a freshwater input in the middle of the bay (e.g., Lirman and Cropper 2003) or nearby bivalve beds (e.g., Sharma et al. 2016), as seascape features are known to be important factors shaping faunal diversity (Boström et al. 2006;Cuvillier et al. 2017;Perry et al. 2017;Rees et al. 2018;Tarimo et al. 2022).
We captured breeding adult pipefish, which are considered flagship seagrass species (Shokri et al. 2009) and indicators of healthy seagrass meadows (Shokri et al. 2009;Ford et al. 2010;Scapin et al. 2018). Breeding adults were absent in the winter, but juveniles remained likely to be captured (Fig. 7). These patterns suggest that the local pipefish population could be annual with adult mortality at the end of the summer, a hypothesis supported by a short life span observed for S. nigra in Australia (Parkinson et al. 2012;Parkinson and Booth 2016). Alternatively, the local pipefish populations might not be obligate seagrass dwellers, in contrast to consensus data (Hyndes et al. 2018), and adults could migrate to deeper waters and other habitats (e.g., kelp beds) upon completion of breeding. However, this migration hypothesis deviates from patterns observed in S. nigra in Australian seagrass beds (Connolly 1994;Hindell et al. 2000;Burt 2002;Browne et al. 2008;Smith et al. 2011). Regardless, we observed a shortened breeding season compared to the Australian populations, where breeding adults were found throughout the winter (Parkinson and Booth 2016), potentially supporting the Predicted probabilities (± standard errors) of each species being present as adults (A, C), or present as juveniles (B, C) in each season and across each level of seagrass cover (represented by different shapes and colours). Note the different y-axes. Adult pipefish were most likely to occur in spring and fall, whereas juveniles had similar likelihoods of occurring year-round in both species. Stigmatopora nigra juveniles were more likely to be caught than Leptonotus elevatus juveniles in every season. Both species were most likely to occur in dense or patchy seagrass habitats hypothesis that populations of S. nigra in New Zealand and Australia have diverged sufficiently to be considered different species (Dawson 2012). A shortened breeding season could also contribute to stronger sexual selection (Emlen and Oring 1977) and thereby increase genetic divergence.
In conclusion, we have shown that the animal communities inhabiting a temperate seagrass bed varied across seagrass patches and seasons, with rarer species showing high temporal turnover. We also suggest that poorly studied seagrass beds on the South Island of New Zealand function as nursery grounds for flounder and habitat for iconic pipefish. These results highlight the importance of sampling many times over a year to capture all communities. Although the most abundant species (shrimp and flounder) dominated all habitat types, seagrass cover was an important factor facilitating several less common species, in particular pipefishes of all ages. Seagrass beds therefore enhance the biodiversity of the system. Seagrass abundance is declining globally ) and in New Zealand (Turner and Schwarz 2006;Tan et al. 2020), and protecting this habitat is critical for maintaining biodiversity. Our finding that seasonality interacts with seagrass cover to shape community diversity further highlights the importance of conserving this habitat.

Acknowledgements
This research was conducted in the Akaroa Taiāpure, within the rohe of Ōnuku Rūnanga, which represents Ngāi Tārewa and Ngāti Īrakehu, whom we thank for their support of this research. We acknowledge Ngāi Tūāhuriri upon whose lands the analyses and writing was conducted. This research was primarily funded by the Brian Mason Trust (grant number 201824) with supplemental funds provided by the University of Canterbury. Thanks to Averill Moser-Rust for her assistance with fieldwork and to Environment Canterbury for funding her summer scholarship. Additional field assistance was provided by Renny Bishop and Kyle Andales. All research was conducted under the MPI Special Permit numbers 590 (amendment 2) and 728, and we thank Dr. John Pirker for managing and supporting those permit applications. Sincere thanks to Dr. Zach Marion for his lively discussions about the statistical analyses conducted, for sharing his code, and for feedback on early drafts of the manuscript. We also extend our thanks to Dr. Warwick Allen, Dr. Will Godsoe, and Professor James Fordyce for feedback on early drafts of this work, and two anonymous reviewers and an editor for helping improve previous versions of the manuscript.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.
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/.