Freshwater input and tidal position regulate species turnover and interaction rewiring in intertidal ecological networks
School of Biological Sciences, University of Canterbury, Private Bag 4800, Christchurch 8140, New Zealand
UWA Oceans Institute & School of Biological Sciences, University of Western Australia, Crawley, Western Australia, Australia
Department of Ecoscience, Arhaus University, Roskilde, Denmark
Te Pūnaha Matatini Center of Research Excellence, University of Canterbury, Private Bag 4800, Christchurch 8140, New Zealand
*Corresponding author, Email: aj.gillis12@proton.meAbstract
The effect that environmental conditions have on community and network assembly processes remains unclear, in part because these processes operate at multiple scales. Because marine primary producers and microinvertebrates have limited mobility, are susceptible to multiple stressors, and can be observed interacting in situ, their habitat-based interactions provide an informative system for disentangling network organising processes. We sampled 646 habitat-use networks, quantifying interactions involving ‘habitat-users’ and biogenic ‘habitat-formers’ over 12 months at 9 sites within Te Ihutai/Avon-Heathcote estuary in Christchurch, Aotearoa New Zealand. Using generalised dissimilarity mixed-effect models, we examined whether changes to species interactions – deconstructed into species turnover and interaction rewiring – were modulated by environmental covariates, including freshwater discharge, elevation, temperature, spatial location and season. We found that with increasing dissimilarity in sites proximity to freshwater, interaction change was more driven by rewiring, whereas differences in elevation (i.e., between channels and non-channel habitats) were driven by species turnover, with more sessile species inhabiting tidal channels. The proximity of habitats also played a strong role, with nearby networks comprising more similar interactions, and species turnover becoming more prevalent with increasing distance. Our results highlight that the relative influence and magnitude of rewiring and species turnover in controlling estuarine interaction networks was affected by the individual species distributions across the estuary and their responses to separate, but co-occurring, environmental factors. Quantification of habitat-former/user interaction networks offers robust, albeit understudied, measures of processes that can underpin community assembly, highlighting their potential importance in research, management and conservation.
Open research statement
Data are provided for peer review. The code for produced from the project analysis and used to draft this manuscript is shared via a public GitHub repository hosted by the Tonkin Research group (repo name: EstInteractTurn). The data used in the formal analysis is hosted on Zenodo, under the corresponding authors profile (doi: 10.5281/zenodo.20619091). The data and code was prepared following strict adherence to the FAIR principles, meaning all data was saved as comma-separated values (.csv) or native R data structures (.rds).
Article notes
Competing Interest Statement
The authors have declared no competing interest.
Summary of Updates:
Footnote Group
Introduction
Understanding how community assembly is influenced by habitat selection, species interactions, and environmental conditions is a fundamental goal of community ecology (Chase, 2003; Drake, 1990; Flores-Arguedas et al., 2023; Fukami, 2015). The assembly of and changes to species interaction networks, in particular, depends on various processes that operate at multiple spatial and temporal scales (Marjakangas et al., 2022; Ponisio et al., 2019). Disentangling these processes is critical to better understand how interactions networks vary across space and through time (Carstensen et al., 2014; Simanonok & Burkle, 2014), but requires approaches that account for changes to both species composition and interactions. Overall changes in species interactions (i.e., the total interaction turnover), can be partitioned into ‘species turnover’ (βst), and ‘interaction rewiring’ (βrw), reflecting changes in both composition and interaction partners (within a shared species pool) between networks, respectively (Fig. 1; Fründ, 2021; Novotny, 2009; Poisot et al., 2012).
The contribution of interaction turnover mechanisms to network assembly depends, in part, on species-specific responses to key environmental factors (Bascompte & Stouffer, 2009; Chase, 2003; Ponisio et al., 2019). Species turnover can contribute to interaction turnover when the species composition of networks changes due to niche selection or partitioning along strong deterministic environmental gradients, such as altitude (with associated differences in temperature and precipitation; Sekar et al., 2024; Simanonok & Burkle, 2014). In contrast, rewiring, i.e., partner switching among the same pool of species, can contribute to interaction turnover when the relative abundances of some species change along spatiotemporal gradients in environmental conditions, and resource availability (CaraDonna et al., 2017; Ceron et al., 2022; Lázaro & Gómez-Martínez, 2022). Given the varied species responses to environmental conditions, it remains unclear how the different turnover mechanisms contribute to network assembly, particularly in spatiotemporally dynamic systems (Pellissier et al., 2018; Ward et al., 2026).
Coastal bar-built low-energy estuaries are complex, interconnected, and highly variable habitats, formed and maintained by geophysical processes, including the deposition of sediment from rivers and the ocean, and fluctuations in sea levels. Estuaries are generally dominated by soft sediments and burrowing infauna, like worms and bivalves (Levin et al., 1996). By contrast, estuaries have little hard substrate, such as rocks and boulders, for seaweeds and ‘epifaunal’ macroinvertebrates to attach to – epibiotic life strategies that instead dominate on rocky reefs.
Nevertheless, estuarine organisms can form complex epibiotic interaction networks, where and when habitat-forming organisms, including seaweeds or cockles, provide attachment, feeding and hiding space for habitat-using organisms, such as snails and crabs (Clemente & Thomsen, 2025b; Manca et al., 2022; Thomsen et al., 2016). Furthermore, the limited mobilities of the many epibiotic seaweeds and invertebrates that make-up habitat-based interaction networks may increase their susceptibility to stress from multiple co-occurring environmental factors, potentially revealing new insights into network assembly processes (Domínguez et al., 2021; Krug et al., 2021; Wang et al., 2024). For example, the vertical tidal position (e.g., along a tidal flat-to-tidal channel gradient) may filter species that are susceptible to desiccation from networks (reflecting species turnover; Chang et al., 2018). In contrast, rewiring may become more prevalent along a freshwater-influenced salinity gradient where stenohaline snails become more abundant and use a wider variety of habitat-formers further from freshwater sources. Finally, metacommunity and spatial processes, like dispersal, may affect how functionally connected communities are, thereby modualting interaction similarity between networks (Novotny, 2009). Because disentangling the complexity of the competing processes that underpin network assembly is challenging with traditional dissimilarity-based methods, modern hierarchical approaches, such as Generalised Dissimilarity Mixed effect Models (GDMM; Hernández-Carrasco et al., 2026; White et al., 2024; Woolley et al., 2017), are required.
Here, we combined network theory and GDMMs to investigate the mechanisms underpinning assembly of interaction networks and their variability in response to the predictability of key stressors in a dynamic estuarine ecosystem. Using networks formed by intertidal habitat-former-user interactions, we assessed how the strength of species turnover and rewiring was modified by pairwise dissimilarities in distance to nearest river, tidal elevation, season, air temperature, freshwater discharge, and inter-sample distance. We hypothesised that rewiring would become more important to overall interaction turnover with increasing dissimilarity in distance to freshwater source (H1a) and freshwater discharge (H1b) because freshwater inputs may have opposing effects on different species (e.g., near-river sites with low salinity and high organic loading may support euryhaline filter-feeders but suppress stenohaline grazers). We also hypothesised that species turnover would contribute more to the overall interaction turnover along strong deterministic filtering conditions, including tidal elevation (H2a), temperature (H2b) or season (H2c), partly because many species are poorly adapted to desiccation stress, which is exacerbated at high air temperatures during the summer months.
Lastly, we hypothesised that networks that are close to each other in space will have relatively similar species pools and interactions (hereafter referred to ‘similarity’) and will exhibit increasing species turnover with increasing spatial distance (H3).
Methods
Study site
Ihutai/Avon-Heathcote is a ca. 7 km2 bar-built tidal estuary located in Ōtautahi/Christchurch in Aotearoa/New Zealand (Fig. 2a; Foster, 2019; Gerber, 2021). The estuary receives freshwater (and associated nutrients and sediments) from the northern Ōtākaro/Avon and southwestern Ōpāwaho/Heathcote rivers and fully saline water through a small inlet located in the southeast corner (Fig. 2b). These sources of fresh and marine water create a salinity gradient in the estuary from ca. 4 to 34 ppt, depending on location and temporally varying tidal and riverine flow regimes (Gerber, 2021; Jones & Simons, 1981; Marsden, 2004). Moreover, tidal channels which remain submerged during low tide, and tidal flats, which typically are emerged during low tide, create a strong desiccation gradient (Foster, 2019; Gerber, 2021; Jupp et al., 2007). The tidal amplitude varies from 1.7 to 2.2 m, exposing large tidal flats to atmospheric conditions for 3-5 hours during low tide (Foster, 2019; Gerber, 2021; McClatchie et al., 1982). In the estuary, bare sand and mud flats are interspersed with biogenic habitat-formers, like the seagrass Zostera mulleri, the seaweeds Ulva spp. and Gracilaria chilensis, high densities of the endemic cockle, Austrovenus stutchburyi (which can often protrude partly above the sediment surface to provide an attachment space for epibiota) and scattered dead shells deposited on the sediment (Clemente & Thomsen, 2024, 2025a, 2025b). Z. muelleri, Ulva spp., and A. stutchburyi, are archetypical estuarine foundation species that create critical habitat and attachment substrates for epiphytes (e.g., filamentous algae), epifaunal sessile species (e.g., barnacles), and epifaunal species with limited mobility (e.g., limpets; Clemente & Thomsen, 2025a, 2025b; Gerber, 2021).
Sampling of habitat-former/user interactions
Using an established habitat-former-user interaction sampling protocol (see Montie et al., 2024; Montie & Thomsen, 2023; Thomsen & South, 2019), we recorded interactions between habitat-formers and habitat-users ca. two hours before low tides, each month from November 2019 to October 2020 (except for April 2020, where sampling was prohibited due to COVID-19 restrictions). Each month, we haphazardly placed 0.0625 m2 georeferenced quadrats at nine sites positioned along the estuary’s eastern and southern shores (a total of ca. 60 samples per month, all samples were separated by > 2m, Fig. 2b). To quantify the role of vertical elevation, approximately half of the quadrats were sampled on the tidal flats and half in tidal channels. Within each quadrat, we recorded the interactions between the biogenic (e.g., Z. mulleri and A. stutchburyi) and non-living (e.g., mud, dead shells) habitat-formers and their associated epibiotic habitat-users. Habitat-users could either be physically attached to the habitat-former (e.g., sessile species like barnacles or slow-moving species like limpets) or living directly underneath the habitat-formers (e.g., many snails and crabs, underneath seaweed fronds). We did not quantify infauna in the sub-surface sediments (typically polychaetes and bivalves) in this study. All species were identified, and their specific habitat interaction were recorded, such as Elminius (a barnacle) attached to A. stutchburyi or Micrelenchus huttoni (a snail) attached to Z. muelleri. In addition to two-species interactions, we also recorded longer interaction chains, such as Membranipora spp (a bryozoan) on M. huttoni on Ulva spp. on A. stutchburyi (an example of a rare four-level chain).
Environmental and spatial covariates
We used the proximity to nearest freshwater sources (i.e., distances to the two river mouths) and LIDAR-derived vertical elevation (relative to mean sea level) as proxies for salinity and desiccation stress, respectively. Using the sf package in R (Pebesma, 2018), we calculated the minimum distance (m) to the nearest river mouth and extracted elevation (m) for each quadrat from a 1-metre resolution digital raster elevation model obtained from Land Information New Zealand (LINZ: https://data.linz.govt.nz/layer/121859-new-zealand-lidar-1m-dem/). To account for the potential role of dispersal, we also calculated the spatial distance (m) between all quadrats using a spatial distance function from the sf package in R (Pebesma, 2018). Furthermore, hourly air temperature (°C) data were downloaded from the Bromley weather station (agent no: 43967) using the National Institute of Water and Atmosphere (NIWA) data repository DataHub (https://data.niwa.co.nz/products/climate-station-hourly). Additionally, daily freshwater flow data (m3/s) were downloaded from the Environmental Canterbury (ECAN) Open Data portal (https://data.dev.ecan.govt.nz/Catalogue/Method?MethodId=79) from the Buxton Terrace and Gloucester Street Bridge flow gauges, representing freshwater flow from the Heathcote and Avon rivers, respectively. Using the downloaded freshwater and temperature data, we then calculated the mean freshwater flow and mean air temperature for the month prior to sampling to account for any prolonged time periods of high flow and temperature that may repeatedly affect the epibiota communities.
3.3Statistical analysis
Patterns in community and interaction composition
We visualised the multivariate community for the entire estuary (pooled across time) and monthly variability (with 95% confidence ellipses) at specific sampling sites with non-metric multidimensional scaling (NMDS) ordination, using the Sørensen dissimilarity index with the ecodist R package (Goslee & Urban, 2007). We visualised network structure with chord diagrams by pooling network adjacency matrices at the site level and arranging them by increasing distance to the nearest freshwater river (Fig. 2b), using the circlize R package (Gu et al., 2014). We then used a generalised linear model to assess whether interaction count and richness varied across months, tidal position and sites, which were ordered based on their minimum proximity to a river using the function glmmTMB and a negative binomial distribution (nbinom2; M. Brooks et al., 2017). Last, we conducted an analysis of deviance with type III sums of squares and a χ2 test using the car package and the Anova function to assess the amount of variation explained by each factor in our glm.
Interaction turnover, covariate distances, and GDMM
The composite of two networks can be decomposed into the contributions of (a) interaction rewiring, (b) species turnover (st) and (c) similarity (sim, Fig. 1; Fründ, 2021). Following Frund (2021), we quantified rewiring (rw), species turnover (st) and similarity (sim) between quadrats i and j using the Sørensen-Dice dissimilarity index (DSØR = (b + c)/(2a + b + c)) as: where rw, st, and sim are the counts of interactions attributed to rewiring, species turnover and similarity, the latter of which is denoted as 2a in the standard index formula (Fig. 1).
To better understand the relative contribution of rewiring and species turnover to overall network changes, we modelled their contribution to network dissimilarity simultaneously using a multinomial likelihood within a generalised dissimilarity mixed model (GDMM), an extension of generalised dissimilarity models (Dias et al., 2022; Ferrier et al., 2007; Mokany et al., 2022; see Woolley et al., 2017 for use of binomial distribution). We modelled observed pairwise changes in interaction counts (rw, st, and sim, the numerators of Eq. 2 and 3) as coming from an underlying probability vector (θij), and the total number of interactions observed in both sites (nij, the denominator of Eq. 1-3).
The above probability vector represent the expected proportion of interactions contributing to turnover or rewiring (θrw,ij, θst,ij), which is mathematically equivalent to their respective contribution the Sørensen dissimilarities ((θrw,ij = E(Drw,ij), θst,ij = E(Dst,ij); Eq. 2 and 3). Therefore, to estimate θij as a function of environmental covariates (m), we defined our linear predictor as:
Here α is a global intercept, and β is the modelled coefficient for each covariate included in the model and Δx represents the quadrat-specific pairwise environmental distance as the absolute difference in predictors (proximity to freshwater, elevation, and the previous month’s mean flow and temperature (|xi – xj|; Mokany et al., 2022; White et al., 2024). For seasonal distances, we calculated the absolute minimum difference in days between quadrat sampling dates, where distances ranged between 0 and 180 days to capture seasonal circularity. For example, two observations recorded 350 days apart should be more similar in seasonal conditions (despite their distant in time). Additionally, inter-quadrat spatial distance was included as a model-predictor to account for spatial autocorrelation and the role of spatial processes, like dispersal. Lastly, we addressed non-independence of pairwise observation by incorporating a latent variable z denoting sample-level (here, quadrat-level) errors into the model structure, and included the absolute pairwise differences in the linear predictor (|zi - zj|; White et al., 2024).
We mapped our linear predictors (η) for each mechanism to the probability vector (θij) using the softmax function:
Interaction similarity (θsim,ij) was designated as our reference category in the multinomial distribution because our hypotheses are focused on the comparative relationship between species turnover and rewiring. Before model fitting, we centred and scaled covariate distance and used the Nelder-Mead optimiser to estimate initial parameter values, which improved chain convergence (Schielzeth, 2010; White et al., 2024). Using the greta package in R (Golding, 2019), we then fit the model with a Bayesian Monte Carlo Markov chain sampler (MCMC), consisting of 4 chains and leapfrog step size set to 25 (min) and 30 (max) to help with efficient sampling. We selected to use uninformative priors for the intercepts and environmental coefficients, such that α ∼ N(0, 0.5) and β ∼ N(0, 0.5). We used regularising priors to alleviate computational pressures associated with hyper-prior distributions, such as z ∼ N(0, 1) × σ and σ ∼ half-N(0, 0.1). Each chain was warmed up with 8000 iterations and sampled for 5000 iterations. We visually inspected posterior distributions and calculated the Gelman and Rubin potential scale reduction factor to ensure model chains converged (R < 1.1; Brooks & Gelman, 1998; Gelman & Rubin, 1992; Vehtari et al., 2021).
Results
Patterns in community composition and network interactions
We sampled a total of 646 quadrats which included 28 living epibiota taxa (i.e., seaweeds, seagrasses, and macroinvertebrates) and 2 types of non-living substrates (dead shells and mud). Visual inspection of the NMDS plot, suggested that both proximity to freshwater and vertical elevation influenced community composition (Fig. 3a). However, temporal and site-level variability contributed to substantial overlaps in the 95% confidence ellipses on the NMDS plot (Fig. 3b).
Across the year of sampling, we observed a total of 127 unique interactions involving the 30 habitat-formers and habitat-users. We found a significant interaction between the variance explained by site and tidal position on interaction count and interaction richness (count – p = 0.5e-11, χ2 = 69.6; richness – p = 0.007, χ2 = 20.8). Broadly, we observed an increase in interaction counts and richness towards the middle of the gradient and at lower tidal positions (see sites S2-5; Fig. 4a, c). Furthermore, we found a significant amount of variation in interaction count and richness explained by month (count – p = 0.1e-17, χ2 = 108.5; richness – p = 0.2e-14, χ2 = 91.7), where both measures were highest during winter (Jun-Aug) and early- to mid-spring (Sep, Oct) at lower tidal elevations (Fig. 4b, d). The most important habitat-formers included dead shells, mud, and Z. muelleri (but this habitat-former had a limited distribution along the eastern coastline only; Fig. 5). Moreover, Diloma subrostrata, M. huttoni and Ulva spp. were observed many times as both habitat-formers and habitat-users, interacting with multiple partners. Finally, M. huttoni, and Notoacmea helmsi were amongst the most numerous habitat-users. Examining our networks revealed that closer proximity to freshwater led to networks dominated by a greater incidence of interactions between A. stutchburyi, D. subrostrata., Ulva spp., and mud. Network and interaction composition became more diverse with increasing distance to the nearest freshwater source, partly because N. helmsii., M. huttonii, and Z. muelleri became more abundant. Dead shell became a more important habitat former further away from freshwater sources where the number of interactions with A. stutchburyi decreased (Fig. 5).
Interaction turnover, GDMMs
We found that with increasing difference in distance to nearest freshwater source, interaction similarity and species turnover decreased, while the contribution of rewiring to overall interaction turnover increased (Fig. 6a). Furthermore, species turnover increased with increasing elevation and temperature distances, while rewiring and interaction similarity decreased (Fig. 6b and c). By contrast, increasing distances in temperature, seasonality (days) and freshwater flow showed small increases in species turnover, and decreases in rewiring and similarity (Fig. 6c-e). Finally, we found that increasing spatial distance resulted in the greatest changes in network structure, where species turnover contributed more to the overall interaction change (Fig. 6f). Conversely, when networks were in close proximity differences between rewiring, species turnover and similarity were marginal (Fig. 6f). When diagnosing our models, we found that our model successfully converged as all R point estimate remained below 1.1.
Discussion
Our results demonstrated that the processes underpinning network assembly in a dynamic estuary differed depending on several factors, including which species were filtered from habitats by co-occurring environmental stressors. Specifically, as networks became more dissimilar in their proximity to freshwater, the probability of interaction rewiring increased, driven, in part, by reciprocal occurrences of specific habitat-forming species (supporting H1a). Conversely, species turnover increased with increasing distance in tidal elevation, likely driven by strong environmental filtering of desiccation-sensitive species, limiting them to lower elevations (supporting H2a). Furthermore, networks that were closer together were more likely to have similar interactions compared to distant networks, where the overall interaction change was driven by species turnover (supporting H3). This latter result suggests that dispersal is limited between locations in the estuary. Finally, using a relatively novel joint modelling method with a multinomial distribution enabled us to identify the relative contributions of the different turnover mechanisms that underpin network assembly in dynamic ecosystems, and how these mechanisms vary in strength with co-occurring environmental stress factors.
Freshwater input, interaction richness and rewiring
Proximity to freshwater source, which alters salinity, nutrients and sedimentation stressors, impacted network structure through increased rewiring, driven, in part, by the reciprocal dominance of key habitat-forming species. For example, A. stutchburyi, an important filter-feeding habitat-former, was relatively more common at sites closer to the two rivers, potentially facilitated by high levels of suspended organic matter and lower salinities, compared to the most marine sites (Marsden, 2004; Sandwell et al., 2009). In contrast, dead shells and photosynthesising habitat-formers such as Ulva spp. and Z. muelleri were observed more frequently at sites farther from the rivers, where suspended sediments was likely lowest and water clarity highest (Clemente & Thomsen, 2025a; Siciliano et al., 2019). In response to these shifting habitat-formers, relatively stress tolerant habitat-users, such as D. subrostrata or N. helmsi, likely selectively associated with the locally dominant habitat-formers. However, stenohaline habitat-users, such as M. huttoni, also contributed to the increase in rewiring, as they were observed more frequently in networks closer to the marine inlet and interacted with a greater diversity of habitat-formers.
As a result of these reciprocal species responses and sustained by the positive habitat feedback mechanisms (Bulleri et al., 2016), the sites in between rivers and the ocean had the greatest overlaps in interacting species distributions and therefore supported a high number of total and unique species interactions (Roshni et al., 2022; Sharpe & Baldwin, 2009). This mid-domain spatial overlap in species niches and increased rewiring mirrors patterns observed in subalpine plant-pollinator networks during mid-season, when a higher phenological overlaps occur in flowering plant species and pollinators that diversify their interaction partners (CaraDonna et al., 2017).
Tidal position and elevation as a strong community filter
In contrast to the reciprocal species responses to differences in freshwater proximity, vertical tidal position limited several species from higher elevation thereby resulting in high species turnover among tidal positions. Mirroring the stressful conditions associated with systems along altitudinal gradients, where extreme temperature and humidity are typically limiting factors (Kraft et al., 2011; Leahy et al., 2024; Simanonok & Burkle, 2014; Sponsler et al., 2022), high elevational tidal flats comprised low interaction counts and interaction richness. More specifically, we found that many sessile species, including barnacles, bristle worms, and hydrozoans, which cannot mitigate heat and desiccation stresses by moving to damp microclimates, were largely absent from networks on the tidal flats (Fig. S1). These results were unsurprising as desiccation and warming are archetypal stressors known to affect several marine species, limiting their distributions to lower tidal positions (Boese et al., 2005; Harley, 2011) and thus altering community structure between tidal elevations (Amstutz et al., 2021; Hillebrand et al., 2010).
Temporal and seasonal variability
Temporal changes and seasonality only had minor effects on structuring networks in the Avon-Heathcote estuary, as reflected by the negligible responses of different turnover mechanisms to seasonal distance. The lack of seasonal turnover may have been driven by the high abundance of relatively long-lived macroinvertebrates and clonal seagrasses that buffer seasonal fluctuations, (Clemente et al., 2023; Curtis et al., 2000), compared to a fewer number of ephemeral or seasonal algae species in the communities (Clemente & Thomsen, 2024, 2025a; Gauna et al., 2017). Indeed, the monthly variability in interaction counts may have been driven more by short-term fluctuations in species-specific abundances, such as following local storms or anomalous river flows, which do not translate to seasonal changes in network structure. Similarly, community composition did not exhibit clear or consistent temporal changes either within or across sites, and monthly variability was minor (except for site E4). Perhaps the lack of seasonal changes in these estuarine networks could be expected, as New Zealand more broadly is considered to have relatively muted seasonality in many environmental conditions (such as rainfall and temperature), due to its maritime-influenced climate (Hernández-Carrasco et al., 2025; Tonkin et al., 2017, 2018).
Future research directions
Here we highlighted the importance of sessile life history traits in modulating habitat-use interactions in estuaries results from their strong dependence on lower elevation habitats. Clearly, other functional and morphological traits affect species distribution in estuaries, and thus the pool of potential interaction partners, and, even if species can physically interact (i.e., trait matching; Dehling et al., 2014; Peralta et al., 2020; Tylianakis & Morris, 2017). For example, the size of intertidal species, such as gastropods and bivalves, has been tied to thermal tolerance, energy expenditure, and survival, which can limit their distributions among microhabitats (Peck et al., 2009). We generally observed more mobile species in similar proportions across tidal elevations (Fig. S1), but smaller individuals of mobile molluscs, such as A. crenulata, may have indirectly contributed to species turnover due to their limited capacity to provide attachment space on their shells. Future studies of seaweed-macroinvertebrate networks could benefit from measuring and analysing traits known to modify species interactions, including their sizes, feeding modes, life-histories, dispersal capabilities and chemical deterrents (Boström et al., 2010; Pereira & Da Gama, 2008; Veríssimo et al., 2024).
The increasing probability of rewiring across freshwater input gradients suggests a potential role of widespread interaction generalism (high diversity of interaction partners) throughout the estuary. Typically, aquatic macrophyte-animal networks have a relatively high proportion of interaction generalism, which can contribute to the persistence of these networks (Bates & DeWreede, 2007; Manca et al., 2022; Taylor & Cole, 1994). However, several studies have shown that aquatic macrophyte-epifauna networks are often structured by specialist interactions between specific sets of partners (Manca et al., 2025; Montie et al., 2024; Thomsen & South, 2019). Indeed, our observation of high species turnover suggests that the networks in the Avon-Heathcote may also contain more specialist interactions. Still, there is a limited number of analyses on coastal habitat-based interactions (Montie et al., 2024; Montie & Thomsen, 2023; Thomsen & South, 2019), and more studies should therefore quantify how interaction generalist and specialist taxa contribute to the turnover and persistence of networks in dynamic coastal ecosystems.
Finally, we note that this analysis is a first joint modelling approach for quantifying species turnover, rewiring and interaction similarity using a multinomial distribution with GDMMs. In using a multinomial distribution, we could quantify the magnitude, direction of change, and uncertainty of each category’s response to environmental variability while accounting for the response of the other categories. This analytical approach allowed us to more accurately model and better assess each category’s contribution to real-world interaction similarities and turnover along key environmental gradients (CaraDonna et al., 2017; Woolley et al., 2017). Still, by using interaction dissimilarities derived from binary presence/absence metrics, we may have underestimated the contributions of turnover mechanisms compared to dissimilarity components derived from interaction counts (e.g., using link weights). The use of interaction count-derived dissimilarities could potentially have revealed a greater contribution of rewiring along the proximity to freshwater gradients or demonstrated seasonal trends (cf. Fig. 5). Further development of GDMMs could incorporate and test network turnover using more quantitative approaches, such as combing weighted interaction counts with the total abundances of both habitat-formers and -users (CaraDonna et al., 2017; Fründ, 2021; Simanonok & Burkle, 2014).
Conclusion
Our study revealed that co-occurring environmental factors can have varied effects on interacting estuarine species, shaping their distributions and ultimately affecting the processes that modulate network assembly. Differences in proximity to rivers led to reciprocal abundance responses in key habitat-forming species, often resulting in interaction changes among habitat users and dominant habitat forming species. By contrast, tidal position and elevation contributed more to species turnover often due sessile species being limited to lower, damper tidal channels. Dispersal limitation also likely influenced the assembly of these estuarine networks, with a clear influence on species turnover with increasing distance. Ultimately, our study showed that, in the Avon-Heathcote estuary with many co-occurring environmental stressors, the estuarine networks were structured by a few habitat-forming species that facilitated local biodiversity (Clemente & Thomsen, 2025b).
Our findings were supported by applying novel multinomial GDMMs, which jointly modelled the relative contribution of interaction turnover mechanisms and enabled us to more accurately analyse and disentangle the combined contributions of rewiring and species turnover to network assembly (CaraDonna et al., 2017). New research that embeds species traits and disentangles the prevalence of interaction generalist or specialist species will help move the field forward, particularly when coupled with advanced modelling frameworks like those used here.
Analysis of turnover in these habitat-use networks provides valuable information on the processes that underpin network assembly and support coastal biodiversity, and how intensifying environmental variability, including from climate change, may alter these processes (Schleuning et al., 2016; van Dijk et al., 2015; Wernberg et al., 2024; Woodward, 2010).
Acknowledgements
The authors acknowledge the Freshwater Ecology Research Group and Marine Ecology Research Group for their support as soundboards, providing feedback on the development of this manuscript.
Funding Acknowledgements
Anthony J Gillis was supported by a University of Canterbury Doctoral Scholarship, provided by Jonathan D Tonkin through a Rutherford Discovery Fellowship, administered by the Royal Society Te Apārangi (RDF-18-UOC-007). Jonathan D Tonkin also acknowledges funding from Bioprotection Aotearoa and Te Pūnaha Matatini, both Centres of Research Excellence funded by the Tertiary Education Commission, New Zealand. Mads S Thomsen was supported by the New Zealand Ministry of Business, Innovation and Employment (Toka ākau toitu Kaitiakitanga—building a sustainable future for coastal reef ecosystems).
Conflict of interest statement
The authors have no conflicts of interest to report.