Taxonomic and Functional Patterns of Benthic Communities in Southern Temperate Tidal Flats

Coastal ecosystems are vulnerable to anthropogenic disturbances which can cause loss of benthic macrofauna and their ecosystem functioning. Despite the importance of functional assessments for conservation and management, knowledge gaps persist on the generality of how the diversity and functional traits of benthic communities influence ecosystem functioning. We investigated eight sites in three different habitats across ~1,260 km of coastline, to evaluate patterns between taxonomic and functional diversity of benthic macrofauna, and the relationship between benthic macrofauna, functional traits and environmental conditions. A total of 74 benthic macrofauna taxa were identified. Significant differences across sites and season were found for metrics based on taxonomic and functional traits. Multivariate analysis revealed spatial-temporal differences, which were more evident based on taxa than functional traits. Functional diversity also showed spatial and temporal differences and was positively correlated with the number of taxa. The dominant functional traits modalities were deposit feeders, with large (>20 mm) body size, burrowers, bioirrigators, deeper than 3 cm in sediments, and irregular morphology. Novel Generalized Linear Latent Variable Models (GLLVM) uncovered several site-dependent relationships between taxa, traits and environmental conditions. Functional redundancy was lowest in a highly modified lagoon, and highest in a more pristine embayment. The outcomes from this study showed site-dependent patterns of benthic communities based on either taxonomic or functional metrics, highlighting that both perspectives are complementary to obtain a holistic understanding of the functioning in marine sediments under environmental change.


INTRODUCTION
Benthic macrofauna are major providers of ecosystem functioning in marine habitats. They modify soft-sediment habitats through biological processes such as ingestion, digestion, excretion, and bioturbation, which facilitates microbial recycling of nutrients, detoxification of pollutants, and organic matter remineralization (Snelgrove et al., 2014;Shojaei et al., 2015;Caswell et al., 2018;Wyness et al., 2021). Benthic macrofauna also represent a connection between benthic and pelagic ecosystems, and plays an important role in energy transfer to different trophic levels (Pearson and Rosenberg, 1978;Kristensen et al., 2014;Murillo et al., 2020). Furthermore, benthic macroorganisms are often used as bioindicators to assess ecosystem "health" due to their sensitivity to natural and anthropogenic disturbances (Borja et al., 2000;Tweedley et al., 2015).
Descriptive and experimental approaches have highlighted that benthic communities are structured by environmental factors (e.g., temperature, water depth, salinity, sediment type, habitat complexity), biological processes (e.g., competition, predation, bioturbation), and ecosystem engineering by benthic macrofauna (e.g., Reise, 1985;Honkoop et al., 2006;Meadows et al., 2012). These traditional taxonomic-based studies laid the base for functional assessments of benthic fauna (Snelgrove, 1997;Thrush et al., 2006;Snelgrove et al., 2014), which allow the understanding of how changes in benthic biodiversity influence the functioning of an ecosystem. Functional approaches have been increasingly explored to comprehensively understand effects of the alarming loss of biodiversity in terrestrial and marine ecosystems (de Juan et al., 2015;Degen et al., 2018;Gammal et al., 2019;van der Plas, 2019). For benthic communities, the use of functional approaches is a powerful tool to investigate Biodiversity and Ecosystem Functioning relationships (BEF), and how these relationships vary spatially and temporally, or under specific environmental conditions (Baldrighi et al., 2017;Beauchard et al., 2017).
Ecosystem functioning, defined as the combined effects of all natural processes that sustain an ecosystem (Reiss et al., 2009;Gladstone-Gallagher et al., 2017;Degen et al., 2018), is commonly analyzed by using Functional Diversity (FD) measurements. Functional Diversity considers the variation of functional traits occurring across ecological communities of a given ecosystem based on the activities of organisms (e.g., movement, behavior, feeding and reproduction; Díaz and Cabido, 2001;Reiss et al., 2009;Beauchard et al., 2017;Degen et al., 2018). Several indices have been used to quantify functional diversity, however, there is a lack of consensus on which index is the most appropriate (Mason et al., 2005;Villeger et al., 2008;Mouchet et al., 2010;Lam-Gordillo et al., 2020a).
Functional diversity is usually split into at least three components: functional richness, evenness, and divergence (Mason et al., 2005;Villeger et al., 2008), with several other components recently added (e.g., functional dispersion and functional redundancy) (Laliberté and Legendre, 2010;van der Linden et al., 2012;Gammal et al., 2020). Functional diversity, and all its main components, are based on the analysis of functional traits and their modalities (e.g., bioturbation, body size, feeding mode, morphology, living habit, sediment position), where species are clustered into groups with shared physiological and morphological attributes (Bremner et al., 2003(Bremner et al., , 2006. The functional traits and their modalities studied can be selected in accordance with the processes of interest, the ecosystem type, and the spatial and temporal scale of study (Hooper et al., 2005;Wright et al., 2006;Bremner, 2008;Beauchard et al., 2017).
Functional metrics, which are based on functional traits, can be more important to explain ecological processes and ecosystem functioning than taxonomic metrics (Belley and Snelgrove, 2016;Mestdagh et al., 2020). Yet, considering both approaches can provide a more robust and holistic knowledge about the structure of benthic communities and functioning of ecosystems. Recent investigations have applied a combination of taxonomic and functional approaches to understand the influence of benthic macrofauna on ecosystem functioning and support management and conservation efforts (e.g., Hajializadeh et al., 2020;Delfan et al., 2021;Nunes de Souza et al., 2021;Shojaei et al., 2021). Results from taxonomic and functional approaches have been similar (e.g., van der Linden et al., 2012;Wong and Dowd, 2015;Hajializadeh et al., 2020), but distinct patterns based on either taxonomic or functional diversity emerged as well (e.g., Emmerson et al., 2001;Kraan et al., 2013;Frid and Caswell, 2015;Gladstone-Gallagher et al., 2017). These different patterns could result from highly variable relationships between taxonomic and functional diversity subject to the environmental context (Gladstone-Gallagher et al., 2017;Thrush et al., 2017;Gammal et al., 2019), and from presence of key benthic macrofauna functional groups, that is often more important than species diversity per se (Norkko et al., 2013;Thomas et al., 2020).
Despite the importance of functional assessments for conservation strategies (Miatta et al., 2021) and to inform policy and management to ensure healthy coastal ecosystems, several uncertainties still persist in the interpretation and potential use of the outcomes from trait-based and functional approaches in future scenarios of biodiversity loss in coastal ecosystems. Such uncertainties can be reduced with greater understanding of how taxonomic diversity and functional traits present in benthic communities influence ecosystem functioning across different habitats. For example, low functional redundancy can be used to indicate habitats vulnerable to functional loss.
To contribute to the knowledge on patterns of taxonomic and functional diversity and their links with ecosystem functioning (e.g., Gammal et al., 2019;Taupp and Wetzel, 2019;Shojaei et al., 2021), this study investigated benthic communities in soft sediments along the southern temperate coast of South Australia, where traits of benthic fauna have been recently compiled (Lam-Gordillo et al., 2020b), which enabled a comparative assessment of taxonomic and functional perspectives. The aims of this field study were to (i) assess the taxonomic and functional diversity of benthic communities across contrasting habitats (coastal embayment, gulfs, and lagoon), each representing a typical habitat of the southern temperate Australian coastline, and (ii) evaluate the relationships between benthic macrofauna, functional traits and environmental conditions across these habitats. Over two seasons, benthic macrofauna, their functional traits and environmental conditions were assessed in each habitat to provide a comprehensive analysis on their patterns and relationships. It was predicted that (1) taxonomic and functional patterns are distinct across the studied habitats, (2) functional diversity is greater within habitats with a greater number of taxa, and (3) the relationships between taxa, traits and environmental conditions are habitat-specific.

Study Area
The southern Australian coast is the longest east-west temperate coastline in the southern hemisphere, and harbors diverse sedimentary habitats (Short, 2020). The benthic sampling was conducted across eight sites in South Australia, covering three contrasting habitats of this coastline: a coastal embayment (Coffin Bay: Long Beach-LB, Kellidie Bay-KB), gulfs (Spencer Gulf: Port Germein-PG, Fisherman Bay-FB, Gulf St Vincent: Port Parham-PPa, Middle Beach-MB), and lagoon (Coorong: Pelican Point-PP, Noonameena-N) (Figure 1; Table 1).

Data Collection and Laboratory Procedures
Benthic samples were collected in July 2019 (Austral winter) and January 2020 (Austral summer) at the eight tidal flat sites. South  Australia has a Mediterranean climate, with dry summers and winter rain. January is the warmest month of the year with an average air maximum temperature of 29.6 • C, while the coldest month is July with a minimum temperature average of 7.6 • C (Bureau of Meteorology, 2021). All sampling occurred at low tide when tidal flats were exposed and accessible from shore. Tides in South Australia are of a unique mixed tidal pattern with tidal range varying from micro-tidal (Coffin Bay, Coorong) to mesotidal (gulfs). The samples for benthic macrofauna were taken using a handheld PVC corer (83.32 cm 2 surface area), pushed it into the sediment up to 20 cm depth, with 15 replicates haphazardly taken per site. All sediment samples were sieved through 500 µm mesh size in the field and preserved in 70% ethanol until further processing. In the laboratory, samples were sorted and all benthic macrofauna were identified to the lowest possible taxonomic level (i.e., 66.2% to Species, 4.1% to Genus, 28.4% to Family, and 1.4% to Order), and counted. At each sampling site, environmental conditions known for influencing the abundance, composition and distribution of benthic communities were measured (Hillebrand and Matthiessen, 2009;Dutertre et al., 2013;Dittmann et al., 2015;Shojaei et al., 2015). Water temperature ( • C), salinity, and pH were recorded in the water overlying the mudflat, using a Hannah HI98194 multiparameter meter. Sediment samples were taken for analyzing Chlorophyll-a (g/m 3 ), total organic matter content (OM%) and sediment grain size. In addition, sediment pore water was collected for analyzing nutrients (Nitrate, Nitrite, Ammonium and Phosphate). Fifteen replicate samples for each environmental parameter were taken at each site within the same area where the sediment samples for benthic fauna were collected.
Chlorophyll-a (g/m 3 ) was determined using a spectrophotometer (Thermo Scientific, Spectronic 200) following the protocols described by Ritchie (2008). The organic matter (OM%) content in sediment was determined by loss on ignition, first drying the sediment to constant weight, followed by burning in a furnace at 450 • C for 5 h. Grain size was determined by laser diffraction using a particle size analyser (Malvern Mastersizer 2000). Average values for grain size fractions for each site were entered into the GRADISTAT program v8.0 (Blott and Pye, 2001) to obtain the median (D50 µm) and coefficient (sorting σG). Nutrient concentrations (mg/L) of Nitrate (NO 3 -), Nitrite (NO 2 -), Ammonium (NH 3 ) and Phosphate (PO 3 4 -) were determined using a Skalar SAN ++ SFA segmented flow analyser.

Environmental Data Analysis
Environmental data were square root transformed as needed to approximate normality (except OM, pH and salinity), and then normalized prior to multivariate analysis (Clarke et al., 2014). Spearman correlation (Supplementary Figure 1) and variance inflation factors with a cut-off <3 (VIF) (Supplementary Table 1) were analyzed for collinearity among variables and, as no redundant environmental variable was identified, all were included in the analyses. To test for differences between sites and season, univariate PERMutational ANalysis Of VAriance (PERMANOVA) and multiple pair-wise tests were conducted, using Euclidean distance for the single variables in PRIMER v7 with PERMANOVA+ add on software (Anderson et al., 2008). Principal Component Analyses (PCA) were performed separately for summer and winter data to explore spatial and temporal patterns. R software (R Core Team, 2017) and the packages "corrplot" (Wei and Simko, 2017), "fmsb" (Nakazawa, 2019), "vegan" (Oksanen et al., 2019), were used for conducting the analyses.

Selection of Traits and Trait Information
A suite of six functional traits and 29 trait-modalities ( Table 2) were selected. The functional traits selected describe behavioral, morphological, and physiological attributes of benthic macrofauna, and are considered as effects traits as they are directly or indirectly related to several ecosystem functions including nutrient cycling and sediment transport (Lam-Gordillo et al., 2020a). Trait information was obtained from the South Australia Macrobenthic Trait (SAMT) database (Lam-Gordillo et al., 2020b). The SAMT database applied a fuzzy Acronyms are used in Figure 4.
Frontiers in Marine Science | www.frontiersin.org coding procedure assigning scores from 0 to 1, with 0 being no affinity and 1 being high affinity to a trait (for details see Lam-Gordillo et al., 2020b). This resulted in the compilation of three data matrices (1) "taxa abundance by site matrix, " in this case the data collected from our surveys; (2) "taxa by traits matrix, " obtained from the SAMT database; and (3) the combinations of the previous two: "traits by site matrix" (Bremner et al., 2006;Bremner, 2008).

Taxonomic and Functional Analysis
The benthic macrofauna was analyzed for traditional diversity metrics, including the analysis of taxonomic richness (S), Shannon index (H'; log e) and Pielou's evenness (J') for each site and season using the package "vegan" (Oksanen et al., 2019). For the functional diversity (FD) analyses, the benthic abundance data (taxa abundance by site matrix) were log (1 + x) transformed to reduce the influence of dominant taxa without losing the abundance effect. To compare the FD across sites and seasons, the following functional metrics were calculated as a proxy of FD. (i) Functional Richness (FRic), provides the amount of functional space occupied by a community (Mason et al., 2005), i.e., the quantity of traits that are expressed in a habitat. (ii) Functional Evenness (FEve), describes how consistently the taxa abundance is distributed across the expressed traits (Mason et al., 2005). (iii) Functional Redundancy (FR), describes the ratio between FD and H' , when the ratio decreases, FR increases and vice versa (van der Linden et al., 2012), providing information on how common the expressed traits are within a habitat. In addition, communitylevel weighted means of trait values (CWM) were calculated to compare trait expression across the sites and seasons. Functional metrics and CWM were calculated using the package "FD" (Laliberté et al., 2014) in R software (R Core Team, 2017).

Statistical Analysis
To elucidate spatial (sites) and temporal (seasons) patterns for each taxonomic (e.g., S, abundance, H' , J'), functional metric (FRic, FEve, Fdis, FR) and CWM, univariate PERMANOVA were used with Euclidean distance for the single variables, permutation of residuals under a reduced model, sums of squares type III and 9999 permutations (Anderson et al., 2008). In addition, multiple pair-wise tests were conducted if fixed factors (sites, seasons) or interactions (sites x season) were significant to identify which groupings contributed to differences from PERMANOVA main tests. To assess community structure differences between sites and seasons, Principal Coordinates Ordination (PCO) were performed for the taxonomic and trait data. Species density (taxa abundance) was fourth root transformed, and in both taxa and trait data, a Bray-Curtis similarity resemblance was applied.
To assess the relationship between functional metrics and environmental conditions, non-parametric multiple regressions were performed with the DISTLM routine, using Euclidean distances and 9999 permutations (McArdle and Anderson, 2001). PERMANOVA, pair-wise tests, PCO and DISTLM analysis were carried out using PRIMER v7 with PERMANOVA add on. To elucidate the direction of the relationships, multiple Spearman correlation analyses were performed using the R package "ggpubr" (Kassambara, 2020). For assessing the response of benthic taxa and functional traits to the environmental predictor variables (fourth corner analysis), several generalized linear latent variable models (GLLVMs) were performed with the R package "gllvm" (Niku et al., 2020). GLLVM extends the basic GLM, handles overdispersion data, includes latent variables to capture the correlation between species, and considers fourth-corner terms to account for speciestraits-environment-interactions (Niku et al., 2020. Further, the fourth-corner approach includes regression of the multivariate abundance against the function of the trait and environment association (Niku et al., 2019). GLLVMs were constructed for fitting multivariate data using a negative binomial distribution as the best fit model (lowest Akaike information criterion-AIC; Supplementary Table 2) (Niku et al., 2019). Level plots were performed for visualizing the interactions between taxa-traits and environmental variables obtained with the GLLVMs using the R package "lattice" (Sarkar, 2008).

Environmental Conditions
The environmental conditions varied across sites and seasons (PERMANOVA p < 0.01, Supplementary Figure 2; Supplementary Table 3). In general, hypersaline conditions were recorded at the gulfs (PG, FB) and lagoon (N) habitats in summer. Sediment grain size (D50 and sorting) was mostly characterized by fine sand at the coastal embayment (LB, KB) and lagoon (N), while the gulfs (FB, MB) had medium to coarse sand. The PCA analyses showed spatial and temporal variation (summer: 40.9%, winter: 47.81% of variability explained by the first two axes) (Figure 2). In summer, porewater nutrients separated the lagoon habitat located in the Coorong (N, PP) based on Nitrate and Nitrite, and LB and PG based on Ammonium and Phosphate. MB had larger sediment grain size (D50, sorting) and higher sediment organic matter (Figure 2A). In winter, higher porewater nutrient concentrations also separated sites from the Coorong (PP, N due to nitrous oxides), and MB was separated again by sediment grain size (D50, Sorting) and organic matter ( Figure 2B).

Taxonomic and Functional Assessment of Benthic Macrofauna
In total, 74 taxa were found across eight sites in South Australia, belonging to six different phyla. Mollusca was the phylum with the highest number of taxa (42%, 31 taxa), followed by Arthropoda (31%, 23 taxa), and Annelida (23%, 17 taxa), while Cnidaria, Echinodermata and Nemertea were represented  (Figure 3). Significant differences between sites and seasons (i.e., the warmest and coldest moth) were found for all three taxonomic metrics (PERMANOVA p = 0.0001; Table 3).
The number of taxa was higher in summer than winter at all the sites (Figure 3A). In pairwise comparisons between sites, significant differences in the number of taxa were found for PG and N compared to the other five sites, but only in winter (p < 0.01; Supplementary Figure 3A). The Shannon diversity index ranged from 0 to 1.99, with greater values in summer at the coastal embayment habitats (LB, KB), and one site in  the gulfs (PG) and lagoon (N) habitat respectively ( Figure 3B). In pairwise comparisons, the majority of the sites were distinct from each other (p < 0.01; Supplementary Figure 3C), apart from the two coastal embayment habitats in Coffin Bay. Pielou's evenness index ranged from 0 and 0.97, following the same patter as H' , with greater values at the coastal embayment habitats in summer ( Figure 3C). In pairwise comparisons, LB and N were significantly different to the other sites (p < 0.01; Supplementary Figure 3D). The most expressed functional trait modalities in the studied benthic communities, based on community-level weighted means (CWM) analyses of trait values, were deposit feeder (feeding mode; contribution: 0.37%), large (>20 mm) body size (contribution: 0.44%), burrower (living habit; contribution: 0.56%), bioirrigator (bioturbation; contribution: 0.56%), sediment position of deeper than 3 cm (contribution: 0.36%), and irregular morphology (contribution: 0.27%) (Figure 4).
The CWM values of each compiled functional trait varied significatively across all sites and seasons (Table 4; Figures 4, 5). Most of the CWM trait modalities showed significant differences across sites, except for the trait modalities omnivore and hard exoskeleton. In contrast, significant differences across seasons were less evident for trait modalities ( Table 4). In pairwise comparisons, significant differences in CWM trait modalities were also observed among sites and seasons (Figures 4, 5; Table 5, Supplementary Table 6). For example, the functional trait feeding mode showed significant differences in all the trait-modalities at the lagoon habitat (N) compared to all other sites (Supplementary Table 6), and the trait modality subsurface deposit feeder (feeding mode) and surface shallow <3 cm (Sediment position) were significantly different in summer and winter in six of the eight sites analyzed ( Table 5).
Functional Richness (FRic), Functional Evenness (FEve), and Functional Redundancy (FR) varied significantly across sites (PERMANOVA p = 0.0001; Table 3). However, FEve was the only metric not significantly different across season ( Table 3). The greatest FRic values were found in summer at all sites, with greatest FRic values in the gulf habitat at Upper Gulf St Vincent ( Figure 3D). In contrast, the greatest values of FEve were found in winter at the gulf (PG, PPa), and lagoon (N) habitats. In terms of FR (ratio FD/H'), the greatest values were recorded at the two lagoon habitats in the Coorong, showing the lowest functional redundancy (Figures 3E,F). Functional diversity, as FRic, was significant and positively correlated with the number of taxa (R 2 = 0.64, p < 0.01, Figures 6A,B). The ratio of FD/H' (i.e., FR) showed a significant but weak negative relationship with the number of taxa (R 2 = 0.13, p < 0.01, Figures 6C,D). Also, a significant but weak, positive relationship was identified between FEve and the number of taxa (R 2 = 0.17, p<0.01, Figures 6E,F). DISTLM analyses revealed that FRic, FEve and FR were mostly influenced by Ammonium, Chlorophyll a, sediment grain size (D50), sorting, sediment organic matter content, Nitrite and temperature ( Table 6).

Community Analyses of Benthic Macrofauna and Functional Traits
Significant community differences were detected between sites and seasons for both taxa and functional traits (PERMANOVA  Table 7). The PCO analysis revealed distinct communities across sites and seasons with 47.7% of the variability in taxa composition, and 64.3% of the variability in trait composition (Figure 7). Based on taxa, sites in the Coorong lagoon (PP, N) were separated from other sites, while the gulfs habitats were more closely grouped. A separation according to season was found in KB, FB and PG (Figure 7A,  Supplementary Figure 4A). Less distinction emerged based on traits with the most evident seasonal separation in PG and N ( Figure 7B, Supplementary Figure 4B).

Ecosystem Functioning-Relationship Between Benthic Macrofauna, Functional Traits and Environmental Conditions
Several significant relationships between the benthic taxa, their functional traits and the environmental conditions were identified across sites and season (Figure 8,  Supplementary Figures 5-12). In general, the interactions were stronger in summer than winter. The stronger interactions were identified at the coastal embayment habitats in Coffin Bay, while the lagoon habitats showed the weakest interactions between benthic macrofauna, traits and environmental conditions ( Figure 8). The six functional traits and their modalities showed significant interactions with environmental conditions across all sites irrespective of season. At the coastal embayment habitats, the trait modality of small (<0.5 mm) body size was correlated to sediment grain size (D50), sub-surface deposit feeder with Chl a, and hard exoskeleton and hard shell with temperature and salinity (Figures 8A,B). The gulfs at the Upper Spencer Gulf showed significant relationships between the trait modalities filter suspension, small and medium body size, and pH, salinity, and sediment grain size (Figures 8C,D). In the other gulf habitats at the Upper Gulf St Vincent, significant relationships were found between the trait modalities deposit feeder, hard exoskeleton and Chl a, with correlations between large body size and tube dwelling to salinity and temperature, and burrower to Nitrate (Figures 8E,F). At the lagoon habitats in the Coorong, the magnitude of the interactions between benthic macrofauna traits and environmental conditions was lower compared to the other sites. In PP interactions between the trait modalities biodiffusor, surface shallow <3cm and temperature were identified, as well as several trait modalities influenced by Ammonium and Phosphate, while for N the feeding modes filter suspension and sub-surface deposit feeder were influenced by Chl a and salinity (Figures 8G,H).

Patterns of Taxonomic and Functional Metrics
Spatial and temporal patterns of benthic communities, based on taxonomic and functional metrics, elucidated variation in benthic macrofauna diversity and functional traits, suggesting differences in ecosystem functioning across habitats and seasons. The theory proposed by various researchers states that greater taxonomic biodiversity will increase the number of expressed traits, resulting in greater functional diversity, and therefore greater effects on ecosystem functioning (Tilman et al., 1996;Loreau et al., 2002;Reiss et al., 2009). However, correlations between taxonomic and functional metrics have yielded highly variable results, often mediated by environmental context and habitat heterogeneity (Hewitt et al., 2008;Strong et al., 2015;Kokarev et al., 2017;Thrush et al., 2017). In this study,   we identified positive relationships between taxonomic and functional metrics, and ascertained that habitats with greater number of taxa and diversity (H') also showed high Functional Diversity FD (as FRic, FEve and FR), as previously reported in other marine and estuarine systems (e.g., Wong and Dowd, 2015;Hajializadeh et al., 2020;Delfan et al., 2021;Shojaei et al., 2021). Yet, we also found that taxonomic and functional diversity of benthic communities were site-dependent and varied across the two studied seasons, similar to findings reported in other studies (e.g., Wong and Dowd, 2015;Gladstone-Gallagher et al., 2017;Gammal et al., 2019). The differences of FD (i.e., FRich, FEve, FR) across sites and seasons could be determined by several factors: (i) abundance, diversity and taxa, and (ii) specific environmental conditions and benthic habitat characteristics (e.g. sediment organic matter, grain size, sorting), as described in previous studies (e.g., Hewitt et al., 2008;Shojaei et al., 2015;Henseler et al., 2019;Cappelatti et al., 2020). FRic and FR were greater across all the sites in summer compared to winter, suggesting that the expression of traits was greater in summer. The lower FD in winter could be explained by a temporary decrease of benthic taxa (e.g., bivalves, crustaceans, polychaetes) with specific functional traits modalities or redundant taxa (Loreau et al., 2002), similar to seasonal patterns in cold temperate ecosystems in the Northern hemisphere (e.g., Kröncke et al., 2013;Shojaei et al., 2021).
The two sites at the Coorong lagoon showed a distinctive pattern compared to other habitats, with a low number of taxa but greater FR (ratio FD/H'), indicating low functional redundancy. The low number of taxa could be explained by the habitat characteristics and environmental conditions of the Coorong (e.g., high salinity, eutrophication) (Dittmann et al., 2015;Mosley et al., 2020). Low functional redundancy arising from few taxa occupying the available functional space with few common traits shared, can indicate vulnerability to future functioning loss, as suggested by van der Linden et al. (2012) and Gammal et al. (2020). In this case, 14 taxa accounted for the low functional redundancy across PP and N, and were dominated by the polychaete Capitella sp. and the insect larvae Chironomidae, which shared traits related to opportunistic behaviors in disturbed habitats (e.g., free living, scavenger, deposit feeder, surface shallow sediment position).
Functional traits and their modalities also varied across sites and seasons, as a result of changes in the benthic macrofauna. Such spatial and temporal differences in functional traits resulting from environmental conditions and habitat complexity are not uncommon in mudflats (e.g., Wong and Dowd, 2015;Gusmao et al., 2016;Henseler et al., 2019;Hajializadeh et al., 2020;Mestdagh et al., 2020). The multivariate (PCO) analysis for both taxonomic and functional trait composition showed a separation based on the site and season. However, the grouping based on traits was less evident, indicating different patterns of alignment between the two metrics. The functional traits and their modalities were thus more homogenous than the taxonomic composition across sites and seasons. For example, the taxonomic composition of the Coorong lagoon was differentiated from the other habitats, but the multivariate structure of the functional traits in this habitat was similar to the other habitats. Such patterns could result from functional redundancy, when different taxa share few common traits, or new taxa added until all traits are represented, or a combination of both (Schulze and Mooney, 1993;Loreau et al., 2002;van der Linden et al., 2012;Gammal et al., 2020).

Linkages Between Benthic Macrofauna, Functional Traits and Environmental Conditions
In this study, the relationships between benthic macrofauna, functional traits and environmental conditions varied across sites and seasons, potentially indicating that, depending on the benthic composition and trait expression, some ecosystem functioning derived from these relationships may be different across habitats. In the absence of direct measurements, insights into ecosystem functioning can be inferred from knowledge of the linkages between taxa, traits and environmental conditions (Wong and Dowd, 2015;Lam-Gordillo et al., 2020a;Delfan et al., 2021). In our study, the trait modalities bioirrigator, surface modifier ("Bioturbator" trait), burrower, free-living ("Living habit"), deeper than 3 cm and bentho-pelagic ("Sediment position") were commonly expressed. These trait modalities showed strong relationships with the sediment characteristics (D50, sorting) at each habitat, showing that muddy to sandy sediments were most suitable for burrowers and free-living organisms (Liu et al., 2019), which increases sediment oxygenation and nutrient cycling from benthic macrofauna activities (Lam-Gordillo et al., 2020a;Delfan et al., 2021).
The trait "feeding mode" was related with environmental conditions in most of the cases across sites and seasons, as it is fundamental for the structural complexity and trophic status of benthic ecosystems (Pearson and Rosenberg, 1978). In our study, the trait modality of deposit feeder was expressed more in summer. Deposit feeders are generally dominant in muddy sediments (Rhoads and Young, 1970;Hajializadeh et al., 2020), and the sediment grain size (D50 and sorting) was smaller in summer compared to winter. The deposit feeders and grazer trait modalities were also expressed most at sites where high concentrations of Chl a were found, as they feed on microphytobenthos (e.g., Wong and Dowd, 2015;Daggers et al., 2020). The relationship of different feeding mode modalities with sediment conditions and primary productivity could also indicate differentiation in the use of resources, food availability, and prey accessibility across sites and seasons (Norkko et al., 2013;Weigel et al., 2016;Sivadas et al., 2020).
The traits for "body size" and "morphology" also varied across sites and seasons, but the trait modalities large body size and irregular body shape were important at almost all sites. Body size is a relevant trait for assessing ecosystem functioning that can be correlated with other traits and provide insight to processes such as nutrient cycling, sediment reworking and energy fluxes (Norkko et al., 2013;Hillman et al., 2020). Large individuals related most to environmental conditions at the studied habitats, however, small and medium body size were the trait modalities most expressed at the lagoon habitats, probably as a result of the large-scale fluctuations in salinity and eutrophic conditions in the Coorong (Dittmann et al., 2015;Mosley et al., 2020).
Environmental conditions correlated most with benthic macrofauna-traits were chlorophyll a, organic matter, sediment grain size (D50 and sorting) and concentrations of ammonium. Functional diversity (i.e., FRic, FEve, FR) was also correlated with the environmental conditions Chlorophyll a, sediment grain size (D50), sediment organic matter content, temperature, ammonium, and nitrite, supporting the pattern found with the linkages between taxa-traits and environmental conditions. The findings from both perspectives (i.e., correlation between taxatraits, and environment conditions, and functional diversity), suggest that across the surveyed sites the ecosystem functioning mostly occurring includes nutrient cycling, productivity, and sediment stability and transport (Norkko et al., 2013;Wong and , , Hajializadeh et al., 2020Lam-Gordillo et al., 2020a;Delfan et al., 2021).

Implications of Functional Diversity for Conservation and Management
Across the south Australian coast, different patterns in benthic taxa, functional traits and functional diversity were identified, as a result of site-dependent environmental conditions and habitat characteristics. In addition, anthropogenic activities are also shaping the benthic communities and their trait expression. The lagoon showed the lowest functional redundancy compared to other habitats, indicating that the functional traits expressed were less common, and only few taxa occupied the available functional space. It has been proposed that the greater the number of taxa and traits expressed in an ecosystem (i.e., functionality), the greater probability of taxa and traits to persist and maintain ecosystem functioning (van der Linden et al., 2012;Kokarev et al., 2017;Murillo et al., 2020). Our findings could thus indicate that the lagoon is vulnerable to further loss of benthic taxa and structural changes (i.e., ecosystem functioning loss) caused by anthropogenic or natural environmental changes. In contrast, benthic communities in the coastal embayment showed high functional richness and redundancy, suggesting that these sites are more resilient and are more likely to maintain their ecosystem functioning if an event of change (i.e., taxa loss) occurs.

CONCLUSION
This study identified spatial and temporal patterns of benthic communities, based on both taxonomic and functional metrics. Functional diversity and expression of functional traits were site-dependent and different across habitats, which could be explained by the benthic community at each site, the influence of environmental conditions and habitat complexity. Correlations between benthic macrofauna, functional traits and environmental conditions were mostly driven by deposit feeders with large and irregular body organisms, performing bioirrigation and burrowing deep into the sediment. Thus, ecosystem functioning would be most affected by the loss of taxa displaying these traits. Our findings corroborate that using both taxonomic and functional metrics is complementary for conservation and management seeking to maintain biodiversity with the implicit understanding that ecosystem functioning will also be maintained. The outcomes presented here advance the understanding of the relationship between benthic taxa, functional traits and environmental conditions in tidal flats. Understanding those relationships will further enable us to predict how ecosystem functioning changes with biodiversity loss, and could potentially help to improve management to ensure healthy functioning of intertidal benthic ecosystems.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
OL-G and SD conceived the original ideas for the manuscript. OL-G developed the outline, collected, and analyzed the data, prepared the figures and tables and wrote the manuscript with some contributions to drafting from RB and SD and critical refinement by SD. All authors contributed to the article and approved the submitted version.