Resilience and Species Accumulation across Seaﬂoor Habitat Transitions in a Northern New Zealand Harbour

: Biodiversity is crucial for maintaining ecosystem stability and functionality under increasing anthropogenic stress. Part of this resilience comes from having many species performing the same function (functional redundancy) leading to the quantiﬁcation of community composition and functional redundancy in relation to increasing stress. However, much of the research within coastal ecosystems focuses on distinct areas, rather than whole ecosystems. Here, we investigate the relationship between biodiversity and functional redundancy across two environmental gradients (sediment mud content and water column depth) and different habitat types following a survey of benthic macrofauna and sediment characteristics at 24 sites within Whang¯arei Harbour, New Zealand. We observed strong gradients in biodiversity which fragmented communities into fewer species that were a subset of the wider community. The lowest biodiversity was observed at muddy, intertidal and shallow subtidal sites which also had the lowest predicted functional redundancy. We show the stronger inﬂuence of water column depth on predicted functional redundancy than sediment mud content, highlighting the importance of subtidal regions. Overall, our study highlights the importance of studying the individual contributions of different areas in a landscape to characterise effective colonist pool size and how this can be used to predict recovery potential following disturbance.


Introduction
Biodiversity is an important contributor to resilience in ecosystems facing threats from anthropogenic stressors [1,2]. Whilst there are many definitions, a key element of resilience is the ability of ecosystems to maintain functioning while under pressure, that is, to resist change or to recover following impacts [3,4]. Theoretically, diverse ecosystems have more species contributing to the delivery of functions (e.g., primary production, organic matter breakdown, nutrient recycling) than species poor systems. Having multiple species performing similar functions, also known as functional redundancy, means that those functions can continue even if some of the species disappear [5]. Therefore, species losses in diverse ecosystems may be less disruptive to functioning than species losses in depauperate ecosystems [6].
Although the functional redundancy concept makes intuitive sense, it is simplistic because species are not perfectly interchangeable in terms of their contributions to functioning and do not respond identically to anthropogenic stressors. Stress-impacted areas with reduced biodiversity may be less resilient because of reduced functional redundancy [7,8], but could alternatively be more resilient if stressors have already eliminated the most sensitive species, making them resistant to further change [9][10][11]. There is evidence of dynamic stability in impacted systems that have been conditioned to cope with stress [10]. Investigating how biodiversity and functional diversity change across stress gradients allows us to explore these concepts and their associated management implications.
The effects of local species losses can also be modulated by connectivity among patches [12]. For example, local losses may be ameliorated by the arrival of conspecifics from the surrounding species pool. Therefore, high connectivity among patches or populations in a landscape will tend to promote resilience and recovery. However, this is dependent upon the availability of species to recolonise disturbed patches, which is affected by natural gradients and the distributions of stressors across the landscape [13][14][15], and therefore the nestedness of species assemblages within the landscape [16,17]. These concepts have led to interest in characterising local (α-) and regional (γ-) scale biodiversity, and how biodiversity changes between these two scales (turnover and ecological connectivity, β-diversity) [18][19][20].
Estuaries are highly heterogenous environments with strong environmental gradients (e.g., salinity, turbidity and hydrodynamic conditions [21][22][23]). Whilst individual estuaries can form relatively distinct units separated by areas of open coast, tidal and wind-driven hydrodynamic transport of colonists can facilitate high connectivity between different benthic habitats, e.g., [24]. However, much of the research informing our understanding of the influence of biodiversity on functional redundancy and resilience within coastal ecosystems has focused on the intertidal environment. This is in part due to logistical difficulties regarding subtidal sampling. The focus on intertidal flats results in relatively few datasets available which sample across both subtidal and intertidal habitat types and the associated environmental gradients. However, understanding the ecology and biodiversity of subtidal ecosystems is becoming increasingly important as coastal managers aim to better understand the implications of sea-level rise, including the impact of shifts from intertidal into subtidal habitats, and how to best manage these transitions.
Environmental gradients can also be a product of anthropogenic disturbance, which has resulted in the degradation and fragmentation of estuarine habitats worldwide [25]. Although the history of human influence is shorter in New Zealand than in the northern hemisphere~700 years; [26], New Zealand estuaries have also been heavily impacted in the modern era [27][28][29]. One of the largest and most pervasive threats to estuaries within New Zealand has been the delivery of sediments eroded from the land following deforestation and poor land use practices, which has increased sediment deposition rates and elevated sediment muddiness in estuaries [30,31]. Elevated mud content can alter the sediment's organic matter content, cohesiveness, permeability and pore water solute (including O 2 ) concentrations [32,33], changes that influence the abundance and composition of sediment assemblages and the ecosystem functions they contribute [34][35][36].
Whangārei Harbour is a large harbour containing a wide variety of habitats including intertidal and subtidal flats, mangroves, saltmarshes, seagrass beds and extensive channels. Additionally, it has a diverse catchment of urban development, livestock farming and forestry further contributing to strong environmental gradients (e.g., sedimentation rates of >25 mm yr −1 in the inner part harbour and <1 mm yr −1 near the mouth) [37]. Along with the adjacent Bream Bay system, it therefore presents an ideal study area to examine the relationship between biodiversity and functional redundancy across multiple environmental gradients. To achieve this, we examined species accumulation patterns by building habitat-, mud-and depth-specific species accumulations in varying orders, and assessed the ecological connectivity through understanding contributions of different areas to local and regional biodiversity. We then applied a standard metric of functional redundancy developed specifically for soft-sediment systems within New Zealand. Based on our results, we aimed to assess how functional diversity is modified by differing environmental conditions, and the contributions of different areas to the recovery and resilience of estuarine communities.

Study Area and Sampling
Twenty-four benthic sites were sampled in Whangārei Harbour and adjacent Bream Bay on the North Island of New Zealand (Figure 1). Sampling was conducted during austral summer over 2 years, with 10 sites sampled in December 2019 and 14 sites in November 2020. At each site we videoed the seafloor to characterise habitat type and collected sediment cores to quantify benthic biodiversity and sediment characteristics (grain size, mud content, organic content and chlorophyll-a concentration). Either n = 3 or n = 10 large cores of sediment (13 cm internal diameter, 15 cm deep, area of 133 cm 2 ) were collected for subsequent identification and enumeration of all sediment-associated invertebrates >500 µm (hereafter macrofauna). All macrofauna cores were sieved on a 500 µm mesh screen in the field before preservation in 70% Isopropyl alcohol. Five smaller cores of surficial sediment (2 cm internal diameter, 2 cm deep, area of 3.14 cm 2 ) were collected and amalgamated (n = 1 per site) for later determination of sediment mud content (%), organic matter content (%) and chlorophyll-a concentration (µg/g dry weight of sediment). Sediment samples were frozen and stored in the dark until analysis.
Sediment samples for macrofauna were sorted, macrofauna identified to the lowest possible taxonomic level and counted. Sediment samples were analysed for particle size distribution, organic matter content, and chlorophyll-a concentration using previously published standard methods, e.g., [38]. Sediment particle size distribution was assessed by digesting a homogenised subsample in~9% hydrogen peroxide, wet-sieving samples through 2000, 500, 250, and 63 µm mesh sieves, and measuring the dry weights of each fraction. Sediment organic matter content was determined as percentage dry weight lost following ignition in a 400 • C muffle furnace. For sediment chlorophyll-a content, a subsample of sediment was freeze dried, weighed, and homogenised before boiling in 90% ethanol. Concentrations of pigments in the ethanol extract were determined spectrophotometrically, with an acidification step used to distinguish chlorophyll-a from its degradation product, phaeophytin.
Each site was categorised according to habitat type, sediment mud content and water column depth (Table 1). Habitat at each site was determined following analysis of video collected by divers along a single 25 m transect per site (width of view~70 cm; camera~1 m above the seafloor). Structural features on the seafloor included physical characteristics and densities of habitat forming organisms, which were quantified during video analysis and the dominant habitat type determined at each site. Six broad habitat types were identified in this study: "seagrass", shell hash dominated (hereafter "shelly"), polychaete worm tube mats (hereafter "tube mat"), "small burrows", "large burrows", and coarse sand with limited visible epifauna, dominated by infaunal molluscs Austrovenus stutchburyi and Notoacmea scapha (hereafter "infaunal molluscs"). Table 1. The number of sites and individual core replicates within each of the habitat, mud and depth categories. N.B. the three sites with >35% mud were excluded from the depth categories to deconfound the role of mud when investigating depth. Sediment mud content categories were based on thresholds reported in the literature, e.g., [39][40][41][42][43]. Muddiness categories were defined as follows: 0-3%, 3-10%, 10-30%, >35% mud content. The mud content categories represented increasing levels of stress, based on a wealth of New Zealand literature documenting the deleterious effects of high bed sediment mud content on macrofauna abundance, diversity and functions, e.g., [33,44,45]. Water depth categories were defined so that there was a similar number of sites in each category (minimum of n = 3 sites per depth band). Intertidal sites were given a depth of 0 m and the six depth groupings were 0 m, 1-3 m, 3-6 m, 6-10 m, 10-15 m and 15-22 m. No values were on the group boundaries (e.g., 3, 6, 10 and 15 m, and 3 and 10% muddiness). Seabed depths were recorded by divers while on site and later adjusted according to tide height at the time of sampling (tidal range in Whangārei Harbour is~2.3 m during a mean spring tide). Considering the significant influence of sediment mud content on estuarine biodiversity, the three sites within the highest mud grouping (>35%) were removed from all depth categories to deconfound the role of mud when investigating depth.

No. Sites
Although a key objective of the analysis was to determine the roles of biotic habitat type, mud content category, and depth zone on species richness, accumulation pattern, and functional redundancy (see below), the three types of categories created were not entirely independent. For example, none of the deepest sites were muddy, all of the intertidal sites were the same habitat type (infaunal molluscs) and all of the 1-3 m depth category were seagrass. In addition, the muddiest sites were between 1-6 m depth, and all of the muddiest sites made up the large burrows habitat group (see Figure S1 for a map of the locations of sites with each category grouping).

Data Analysis
To visualise the macrofaunal communities within and between sites, a non-metric multidimensional scaling (nMDS) of distances among centroids based on Bray-Curtis resemblance matrices (square root transformed) was conducted. Normalised and transformed environmental variables (depth, mud content, phaeophytin concentration and chlorophyll-a) which have the potential to covary with macrofaunal community structure were overlaid.
A combination of diversity metrics and analyses were used to describe the complexity of macrofauna community structure and composition. Accumulation curves were used to determine taxonomic richness within each category group as a function of sampling effort, using an approach similar to Ugland et al. [46]. The diversity of each group as a whole (γ-diversity) was determined from the richness predicted from an accumulation curve based on 16 cores. Sixteen cores was chosen as it represents the highest number of cores within all groups, and thus allowed for comparisons of γ-diversity among groups given a standardised sampling effort. Local α-diversity was defined as the average species richness within each group. Numerous measures of β-diversity have been reported in the literature, however, we calculated additive β-diversity, defined as γ-α diversity. Thus, we subtracted the average α-diversity of each group from the γ-diversity of that group. Additive βdiversity (rather than γ/α) was used as it more accurately describes heterogeneity from local to larger scales and connectivity between sites, rather than species turnover [47][48][49]. Species evenness (Pielou's evenness), Shannon-Weiner and Simpson's diversity indices were also calculated. A Kruskal-Wallis test followed by a Dunn's test was performed for diversity indices where appropriate to identify any significant differences between groups.
To investigate the dissimilarity of species composition between groups, two approaches were used which focused on the species present/absent. The first involved multiple iterations of species accumulation curve building, where groups were combined in specific orders and new species accumulation curves were generated following each new addition to show the contribution of individual groups (habitat types, mud categories, depth zones) to total species pool richness. Mud and depth groups were combined in ascending order first (from lowest to highest mud, and from shallowest to deepest, respectively) and then descending order. Habitat types (which cannot be ordered numerically) were added in the order of increasing then decreasing average α-diversity. The second method used to analyse dissimilarity in species community composition between groups involved calculations of Jaccard's dissimilarity index. This method was applied to each group based on average (across all replicated cores in a grouping) species composition. The resulting index enables a comparison of community composition by determining what percent of species identified were present in both populations, and thus provides a ratio of how similar (close to 0) or dissimilar (close to 1) a community is.
To assess functional redundancy, a Trait-Based Index (TBI) developed for New Zealand soft-sediment macrofauna was calculated at each site as detailed by Rodil et al. [50]. This index is derived from species richness within functional groups and therefore provides an assessment that is closely related to functional redundancy. TBI scores are based on the richness of species in seven particular trait modalities that were selected for their sensitivity to sediment muddiness and heavy metal contamination [50], the primary stressors influencing macrofaunal communities within New Zealand [51][52][53]: a living position in the top 2 cm, having an erect tube or structure, a mover of sediment from surface to surface (not vertically), sedentary, suspension feeding, worm shaped, of medium body size. The functional traits of all species present at each site were assigned using a collaboratively developed functional traits database held by the National Institute of Water & Atmospheric research (NIWA, New Zealand). Higher TBI index scores indicate higher functional redundancy (i.e., higher richness within the seven trait modalities). The relationship between TBI score and sediment mud content and water column depth were investigated using bivariate scatter plots with a regression line fitted based on lowest Akaike Information Criterion (AIC) value of separate polynomial models up to the 4th degree (Table S1). This allowed the assessment of both linear and curvilinear relationships. A multiple regression was then performed including both mud and depth using best fit procedure to investigate the combined contribution of mud and depth in explaining the variability of TBI score between sites.

Estuary Characterisation and Group Selection
An nMDS ordination (Figure 2) of macrofaunal community data collected in the Whangārei/Bream Bay system showed replicates from individual sites clustering relatively close together, suggesting less variability within sites than among sites. Overlaid environmental vectors highlight the influence of water column depth and sediment mud content in driving the variability among sites. The individual species driving these differences between groups can be found in Table S2. For example, Micromaldane sp. was the most highly abundant taxa in the small burrows and tube mat habitat types, but was absent from all other habitats. In addition, the intertidal group was dominated by the amphipod family Lysianassidae and two mollusc species, Austrovenus stutchburyi and Notoacmea scapha which were absent in the subtidal groups.

Patterns of Diversity within Groups
αdiversity was highly variable between groupings and ranged from an average of 8 to 45 species per group (Figure 3). α-, βand γdiversity increased with increasing depth up until 3-6 m, after which all values stabilised (Figure 3). Reductions in βand γdiversity were observed in communities with >10% mud, and >30% mud for α-diversity. When comparing habitat types, the lowest within-habitat α-, βand γ-diversity values were observed within large burrow and then infaunal molluscs habitats, with all indices increasing across seagrass, shelly, small burrow and then tube mat habitats. These patterns of diversity were reflected when comparing three different diversity indices (Shannon-Weiner diversity, Simpson's diversity and Pielou's evenness index ( Figure S2)). Significantly lower Shannon diversity was observed within intertidal communities compared to all other depth groups (Kruskal-Wallis test, p < 0.001; Table S3), however, few differences were observed when using Simpson's diversity. Pielou's index suggests a slightly higher degree of evenness in the intertidal and lower subtidal habitats. Sediment mud content was strongly related to community diversity, whereby significantly lower Simpson and Shannon diversity values were observed in communities with >35% mud compared to all other groupings (Kruskal-Wallis test, Simpson diversity p < 0.001; Shannon diversity p < 0.001; Table S3). Additionally, communities of 10-30% mud had lower Simpson diversity than those with <10% mud, but Shannon diversity and evenness were comparable (Table S3). Between habitat types, the composition of the large burrows group containing only sites with >35% mud resulted in this group having significantly lower diversity than all other habitat types (Table S3). Infaunal molluscs habitats were the only other habitat type to have significantly lower Shannon diversity when compared to almost all other groups (except tube mats).

Species Accumulation in Relation to Sediment Mud Content, Water Column Depth and Habitat Type
Species accumulation curves highlight the differences in species richness among groups relative to sampling effort (Figure 4). When standardised by sampling effort, species richness was similar in the two lowest mud content categories (0-3 and 3-10% mud) but decreased by 19% and 73% when mud content increased to 10-30% and >35%, respectively. Species richness for subtidal communities deeper than 3 m were relatively indistinguishable from each other and did not reach an asymptote within the sampling range (up to 50 replicate samples). However, species richness was considerably lower at a depth of 1-3 m (−41%) and for intertidal communities (−54%). Species richness varied by habitat type, with a 71% difference between the lowest and highest species richness for a given sampling effort. The lowest richness was observed in large burrow, followed by infaunal molluscs and seagrass habitats, whilst the other three habitat types had considerably higher species richness.

Dissimilarity Patterns within Groups
Whilst species richness may be similar among sites and groups for a given sampling effort, community composition (i.e., the particular species present) can differ widely (i.e., there can be turnover) (Table S2).

Water Column Depth
In Figure 5A, the first group plotted was the 0 m depth group (intertidal), followed by the 0 and 1-3 m depth groups combined (i.e., a 0-3 m grouping). The large difference between the first two curves of Figure 5A (0 m and 0-3 m) indicated that many new species contributed to the increased species richness. Conversely, when the last depth group was added (adding the 15-22 m group to produce a 0-22 m curve), a relatively small increase in species richness resulted. Thus, the species occurring within this group were mostly a subset of those found in the remaining groups. On the other hand, when starting with the deepest depth group (Figure 5B), increases in species richness were less pronounced, with the majority of species included within the two deepest depth groups. These relationships were reflected in the calculation of Jaccard's dissimilarity index ( Figure 6). The lowest dissimilarity (<50%) was observed between the deeper subtidal groups (>3 m) and dissimilarity in species composition was lowest in adjacent groups. For example, intertidal communities had the highest similarity to 1-3 m depth, with dissimilarity increasing with increasing depth.

Habitat Type
Similar to the patterns observed with depth, when starting with the habitat type with the highest species richness (analogous to starting with the deepest group, Figure 5D), few new species were added when including additional habitat groups. Conversely, when starting with "large burrows", the habitat type with the lowest species richness, new species were added with each additional habitat group ( Figure 5C). The majority of habitat groups showed high dissimilarity to each other, with an average dissimilarity of 70% across all habitat groups. Dissimilarity was highest for the least diverse groups: large burrows and infaunal molluscs habitats.

Sediment Mud Content
When comparing species richness within the sediment mud content groups, almost all species occurred in the two groups with the least mud (0-3% and 3-10%; Figure 5E). Few new species occurred within the highest mud groups (10-30% mud and >35% mud). In contrast, when starting with the highest mud group (which was relatively species poor; Figure 5F), large increases in richness occurred when sites with less mud were incorporated into the groupings. Jaccard's dissimilarity index supports the considerable difference in species composition between >35% mud and all other groups, with dissimilarity ranging between 76-83%. The lowest dissimilarity occurred between the 0-3 and 3-10% mud categories, which increased with increasing mud content.

Changes in Functional Redundancy between Groups
TBI scores ranged between 0.26 and 1.62, with 50% of scores over 1, and 88% of scores over 0.4 (Figure 7). Sites with the lowest functional redundancy index values were observed at sites with >35% mud and in the shallower parts of the estuary (intertidal and 1-3 m depth groupings), which were categorised as infaunal molluscs and large burrow habitats (Figure 7). TBI scores were originally developed in relation to species responsiveness to gradients of mud and heavy metal concentration in intertidal areas. Across sites in Whangārei Harbour, sediment mud content was inversely correlated to TBI score, but there was high variability when mud content was below 3% (Figure 8). This resulted in maximal TBI scores at~5% mud, before a strong decrease and thus the lowest calculated TBI scores when sediment mud content was high (>60%). TBI scores were low in the intertidal zone and at depths <3 m, increasing to maximal TBI scores when depth was between 4-12 m, before decreasing at depths greater than 15 m. Combined, sediment mud content and water column depth explained 56% of the variability in TBI score between sites.

Discussion
This study highlights the high heterogeneity of macrofauna communities that can occur within estuarine ecosystems and the influence of different environmental gradients in determining species richness and functional redundancy, which provide insights into resilience and recovery responses. We observed strong gradients in biodiversity, with the highest biodiversity at sandier, deeper sites with more abundant macrofauna community structure. In contrast, the lowest biodiversity was observed at muddy, intertidal and shallow subtidal sites, which were dominated by large burrows and low abundances of macrofauna. As biodiversity is a key contributor to ecosystem resilience, these results suggest that the impact of future stressors such as sea level rise or increased sedimentation are likely to be highly spatially variable across different estuarine habitat types. In addition, our research indicates that a comparative lack of subtidal sampling within estuaries may be resulting in assessments of vulnerability to future change inadequately protecting high value subtidal habitats. This is particularly the case where estuarine monitoring programmes (such as those within New Zealand) focus primarily on intertidal habitats and thus may be substantially underestimating biodiversity and functional redundancy within these areas. Moreover, the high biodiversity observed within the subtidal habitats of this study suggested these areas may contribute the most to the overall species pool, indicative of their potential importance to disturbance recovery dynamics.
A strong driver of biodiversity in our study was sediment mud content, with increased muddiness (particularly >35%) resulting in a reduction of all measures of biodiversity, in agreement with previous research, e.g., [33][34][35]44]. For example, Shannon and Simpson diversity indices were significantly lower for sites with >35% mud, and when compared to 0-10% mud, γ-diversity decreased by 19% and 73% for mud content between 10-30% and >35%, respectively. Low γand βdiversity as observed in the >35% group has previously been suggested to indicate these areas are highly degraded [18]. Furthermore, iteratively combining species accumulation curves showed there were few new species present at sites with >35% mud when compared to all other mud groups. This suggests species within the high mud group were a small subset of those found within the wider community, resulting in less diverse communities dominated by a small number of species that were tolerant of high mud content.
Increased muddiness is often a consequence of hydrodynamic conditions, with quiescent areas of low flow predisposed to increased accumulations of fine sediments. Eroded soil (terrigenous sediment) is considered to be one of key contaminants affecting New Zealand's estuaries, with an estimated 192 million tonnes of soil being lost from the land each year [54]. Over time, as sediment accumulates in estuaries it can result in changes to the cohesiveness of surface sediment, inhibit diffusive and advective transport of solutes resulting in changes to porewater geochemistry, alter food quality through nutritional differences of terrestrial and marine sediments, block filter-feeding appendages, and deter larval settlement [44,45,55]. The resulting changes to community composition have direct effects on ecosystem functioning such that high mud regions have been shown to have reduced nutrient processing capacity (including denitrification) and benthic primary productivity [34][35][36]56].
Whilst the relationship between biodiversity and mud is relatively well established, fewer comparisons exist between biodiversity and water column depth, with the majority of studies focusing on either intertidal or subtidal areas separately. However, within this study, water column depth had a significant influence on community composition and biodiversity, as highlighted by 88% maximum dissimilarity of species composition between depth groups ( Figure 6). Differences were greater within intertidal and shallow subtidal regions (0 and 1-3 m), resulting in significantly lower α-, βand γdiversity values, which were comparable to those found in intertidal areas elsewhere [24]. However, shallower depth groups had higher evenness (Pielou's index) than those >3 m and showed considerable overlap in species composition ( Figure 5), suggesting the potential loss of several highly abundant species which were less tolerant to shallow subtidal and intertidal conditions (e.g., differential flow rates, oscillating emerged periods) such as the polychaete tube mats found only at depths >6 m. Differences in hydrodynamic conditions have previously been associated with changes in biodiversity within intertidal areas globally [57], and have been attributed to intertidal extent negatively affecting biodiversity [24]. This can be a consequence of high wave-exposure and changes in current patterns influencing erosion rates and the dispersal of organisms [58,59]. While the inclusion of only one habitat within the intertidal depth group within this study restricts the ability to make generalisations over varied intertidal areas and to comprehensively compare subtidal and intertidal biodiversity, the comparable biodiversity estimates to those across a number of other intertidal habitats and estuaries [24] suggests the intertidal area within this study is likely to be representative of those found elsewhere.
The strong influence of both sediment mud content and water column depth on the biota was also evident within the habitat groupings. For example, the "large burrows" habitat category occurred exclusively at high mud sites (>35%), the "infaunal molluscs" habitat was restricted to the intertidal zone (0 m group), and the majority of seagrass sites (3 out of 4) were present in the 1-3 m depth grouping. Consequently, these three habitat types had the lowest biodiversity. Within mud and sand habitats, the absence or low abundances of habitat formers can result in a homogenous sediment that lacks 3-dimensional complexity and thus can have lower α-diversity than other habitats [57,60]. Conversely, polychaete tube mat habitats, for example, can increase sediment deposition (including organic matter) which can enhance sediment stability and facilitate deposit feeders [61][62][63]. Reduced structural complexity is also likely to have influenced the seagrass habitat because of the sparse growth form and relatively short blades of the species that grows within New Zealand (Zostera mulleri) [24,64] and the infauna mollusc group which had variable abundances of the cockle Austrovenus stutchburyi.
Analysis of functional groups is frequently used to provide more information on the consequences of species loss than overall taxonomic richness [65,66]. In this study, sites with lower sediment mud content had much higher TBI scores indicating greater functional redundancy (e.g., 1.0 for <3% mud vs. 0.37 for the >35% mud), however there was significant variation around very low mud sites (<3% mud). The high variance around sites with <3% mud and the relationship between sediment mud content and TBI scores observed in Whangārei Harbour (r 2 = 0.28) was similar to that observed previously (negative trend, r 2 = 0.36, across other North Island sites; [50]). The former was attributed to differences in hydrodynamic conditions, such as high wind and wave activity resuspending fine sediments (mud) resulting in unstable and poor-quality food availability rather than a response to anthropogenic impacts. However, it should be noted that TBI scores were much higher overall in this study relative to the more urbanised and human-affected estuaries of the Auckland Region where the TBI index was developed, where >0.4 is considered a high TBI score.
More influential than sediment mud content, this study showed a stronger relationship of TBI score with water column depth (r 2 = 0.49), the first use of the index in this way. Combined with sediment mud content, 56% of the variability in TBI score was explained, highlighting the dominance of these two gradients in influencing community composition. It should, however, be noted that only a small number of intertidal sites were included within this study, restricting intertidal versus subtidal comparisons. Nonetheless, TBI score increased with depth up to~15 m whereby it plateaued, most likely a consequence of the similarity in community composition between the deeper subtidal sites (39% dissimilarity). Consequently, intertidal and shallow subtidal regions were estimated to have lower functional redundancy and are therefore suggested to be more susceptible to new or increased stressors. Whilst this would imply these communities are unable to cope with further stress and thus the effects of new stressors are amplified, e.g., [7,8], it has also been argued that these communities may alternatively be adapted to stress and therefore the response to further stressors may instead be lessened [9][10][11]. Discerning how emerging stressors will interact with current environmental pressures and gradients is therefore essential when trying to predict community responses to future change.
Sea-level rise is already impacting coastal areas globally, and whilst the magnitude is uncertain, it will likely create new pressures on benthic fauna through the deepening of subtidal areas which will shift or squeeze intertidal areas [67]. Within this study we observed significant changes in estuarine biodiversity along a gradient of water column depth highlighting the need to quantify what exists subtidally in order to understand the vulnerability of estuaries to future change. Increases in subtidal habitats have the potential to provide the opportunity for biodiversity enhancement due to increased surface area of highly biodiverse subtidal regions, such as those reported within this study. However, sea-level rise will co-occur with changes in hydrodynamic conditions as well as increases in rainfall, suggesting sedimentation may intensify in the future resulting in increased turbidity and the deposition of sediment [68]. For some habitats, the interaction between these two stressors may inhibit the potential for biodiversity enhancement. For example, seagrass beds are highly sensitive to water column turbidity and sedimentation [69], therefore, sea-level rise would restrict the emerged periods currently providing a refuge from this temporal stressor and thus reduce the viability of these important habitats [70]. Understanding the differences in biodiversity and functionality between intertidal and subtidal areas is fundamental to informing management decisions in order to better facilitate these transitions.
Low diversity areas in Whangārei Harbour including intertidal habitats and sites with >35% mud may be areas of higher ecological stress where only relatively stress-tolerant species survive. Species living in permanently submerged, deeper areas would be less likely to experience thermal extremes, desiccation, and physical disturbance by waves, and therefore may be less resilient to future change than habitats which have adapted or been degraded due to higher stress. Similarly, muddy sites are generally associated with high stress in the form of high suspended sediment concentrations, high rates of sediment deposition, lower levels of sediment oxygenation, and higher concentrations of ammonium and sulphides, e.g., [32,33]. Our analysis suggests the species present in muddy and intertidal sites were a subset of the species present in adjacent sandier and deeper habitats in Whangārei, which likely increases their resilience (via connectivity) and improves their chances for recovery when disturbed (if the stressor is displaced). Our study highlights the importance of studying the individual contributions of different habitats in a landscape to characterise effective colonist pool size and how this can be used to predict recovery potential at disturbed sites.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/d14110998/s1, Figure S1: Locations of the sampling sites as a function of grouping; Figure S2: Pielous' evenness, Simpson index and Shannon Wiener diversity indices by group; Table S1: AIC values of the regression models; Table S2: Ten of the most highly abundant taxa per group; Table S3: Results of the Kruskal-Wallis tests. Institutional Review Board Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.
Acknowledgments: Special thanks to Patuharakeke for co-developing the research with us and to Taryn Shirkey specifically who participated in field work. We acknowledge the Marine Ecology team at NIWA for assistance with sample processing and Jack Massuger and Fabrice Stephenson for diving/boating support. We also thank two anonymous reviewers for constructive comments on the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.