Different functional feeding groups of mangrove soil molluscs invoke unique co‐occurrence patterns in response to a climate extreme

Due to differences in community functional traits and composition, the response of communities to environmental changes can be variable, and thus an understanding of the mechanisms underlying differential functional group responses is necessary. We investigated a mangrove soil mollusc community before (in 2007) and after (in 2008 and 2009) an extreme cold event (ECE) to determine the effects of the ECE on co‐occurrence patterns, spatiotemporal turnover and ecological processes in different functional groups (i.e. in‐faunal shelled deposit feeders [ISDFs] and hard‐bodied, mobile scavenger/predators [MS/Ps]).


| INTRODUC TI ON
Soil, as one of the most biodiverse habitats, contributes to many ecosystem functions and services (Thakur et al., 2020). Extreme climate events, such as extreme cold events (ECEs), influence soil ecosystems and are increasing in intensity, frequency and duration (IPCC, 2020). Therefore, uncovering the factors that make soil communities vulnerable to extreme climate change has become integral to the field of soil conservation biogeography (Coyle et al., 2017;Maxwell et al., 2019;Richardson & Whittaker, 2010). During ECEs, temperatures rapidly decrease to near freezing conditions and typically remain there for several weeks (Ross et al., 2009). Therefore, the threats posed by ECEs to soil-related systems will persist. To date, climate change vulnerability assessments that have included ECEs have predominantly focused on how the taxonomic diversity (e.g. species richness and abundance) of soil communities responds to extreme changes in climate Liu et al., 2016). The findings of these studies suggest that community composition responses to ECEs (i.e. changes in richness and/ or abundance) may be highly variable. In addition, in their response to diverse environmental changes, soil communities comprise different functional groups (Barnes et al., 2019;Dethier & Steneck, 1994).
Although different functional groups coexist at local and regional scales, they impart different effects on soil community dynamics in response to environmental changes (Sundstrom et al., 2012;Voigt et al., 2007). However, knowledge concerning the mechanisms underlying the ecology, evolution and assembly of different soil functional groups in response to climate extremes, such as ECEs, remains limited (Voigt et al., 2007). Mangrove soil molluscs, as macroinvertebrates, are generally found on the surfaces and upper layers of mangrove soil, humus and sediment (Cannicci et al., 2008;Decaens, 2010;Lee, 2008;Thakur et al., 2014). Based on morphospecies of differing feeding types, mobility and skeletonization, mangrove soil molluscs are generally grouped into two functional groups according to how they affect mangrove soil structure and function: in-faunal shelled deposit feeders (ISDFs) and hard-bodied, mobile scavenger/predators (MS/Ps) (Barnes et al., 2019). Recent studies have shown that ECEs influence species richness, abundance and biomass in the mollusc community Liu et al., 2016;Reiss et al., 2006). The early development, juvenile growth, reproduction, physiological metabolism, mortality and distribution of intertidal soil molluscs are influenced by ECEs (Ansart & Vernon, 2003;Liu et al., 2016;Prather et al., 2013). Therefore, the community structure and function of mangrove soil molluscs in different functional groups should be influenced by ECEs.
Abundance-based and trait-based co-occurrence patterns can characterize the structure and function, respectively, of the soil community (Mouchet et al., 2010;Suding et al., 2005). Spatial turnover (i.e. spatial dispersal), temporal turnover (i.e. ecological succession) and local environmental filtering are critical ecological processes that shape the regional and biogeographical co-occurrence patterns of soil fauna (Ainsworth et al., 2020;Lv et al., 2014;Suurkuukka et al., 2012;Wu & Wang, 2019). Therefore, scale dependency, including environmental and spatiotemporal trends continues to be of pivotal importance to conservation biogeography (Suurkuukka et al., 2012). Recent studies have increasingly emphasized the following points: both spatiotemporal turnover and local environmental filtering of communities are influenced by climate change; one of these processes may dominate depending on the particular case; and trade-offs exist between these processes (Araújo et al., 2011;. Therefore, climate disturbances, especially extreme climate changes, may influence the assemblages of soil communities by modulating both spatial and temporal processes (Stegen et al., 2013). In addition, the mixture of environmental filtering and spatial and temporal processes affecting communities indicates that assemblages with different properties or traits might respond to climate changes via different ecological strategies Pandit et al., 2009). For instance, previous studies have shown that abundant arboreal mangrove mollusc taxa mainly increased divergent succession and spatial limitations to cope with ECEs; the rare assemblages maintained high levels of divergent succession but eroded spatial limits in response to ECEs . However, the limitations of taxonomic approaches (e.g. abundance based and species based) are obvious because species within different functional groups differ in their responses to environmental changes, such as climate change (Mouillot et al., 2013). Therefore, it is important to differentiate between the different functional groups to better understand the mechanisms underlying soil community functioning in response to extreme climate changes.
To address this knowledge gap, we studied soil mollusc (benthic mollusc) communities that involved five mangrove wetlands. The mangroves in our study system experienced an ECE in 2008 (Boucek et al., 2016;, providing a unique opportunity to evaluate the effects of the ECE on the different functional groups of a soil community. We calculated the taxonomic and unique co-occurrence patterns in response to an ECE. We highlight the value of dividing different functional groups in developing predictive frameworks for a better understanding of community functional responses to global change.

K E Y W O R D S
biogeographical co-occurrence pattern, divergent succession, extreme cold event, functional group, mangrove wetland, soil mollusc, spatial turnover functional diversity of ISDF and MS/P sub-communities based on species identity, species abundance and body mass. We aimed to determine how the taxonomic and functional diversity of soil mollusc communities with different functional groups responded to the ECE, including changes in regional co-occurrence patterns, temporal succession, spatial turnover and environmental filtering. Low temperatures can affect the early development, juvenile growth, reproduction, physiological metabolism and mortality of coastal soil fauna (Reiss et al., 2006). Therefore, ECEs can dramatically affect the species' richness, abundance and biomass of coastal invertebrate communities (Chen, Gu, Liu, Shi, et al., 2021;Liu et al., 2016) and then change the community composition, increase modules of the co-occurrence network pattern, enhance the spatiotemporal turnover of mangrove soil mollusc communities and break down the balance between deterministic and stochastic processes in mangrove mollusc communities (Chen, Gu, Liu, Shi, et al., 2021;Ferrenberg et al., 2015;Legras et al., 2019;Lundquist et al., 2018). Additionally, communities with different feeding types or mobility might respond to environmental disturbances (e.g. the ECE) via different mechanisms (Lundquist et al., 2018;Pandit et al., 2009;Voigt et al., 2007).
For example, species with a strong athletic ability can avoid cold waves by actively moving, while species with a weak athletic ability cannot (Coyle et al., 2017). In this study, we tested the following hypotheses: (I) ECEs increase the modularity of co-occurrence patterns of mangrove soil molluscs, and different changes are exhibited between ISDFs and MS/Ps; (II) the ECE increases the spatiotemporal turnover of community, and different changes occur between ISDFs and MS/Ps; and (III) the ECE is expected to break the balance between spatial processes and environmental filtering in the mollusc community, and different changes occur between ISDFs and MS/Ps.

| Study area
The study region comprised Qinglan Bay (QL), Dongzhai Bay (DZ), Techeng Island (TC), Yingluo Bay (YL) and Beilun Estuary (BL) in ; Figure S1a). This region is located in the tropical and southern subtropical zones and has a warm and humid maritime climate, with an annual average air temperature of 23-24°C and annual average precipitation of 1500-1700 mm. The average air temperature in January (the coldest month) is generally 16-17°C. The daily average air temperature minima fall below 10°C on fewer than 7 days per year under normal circumstances. The annual average sediment temperature is 24.2°C.
Previous studies have shown that the soils (sediments) of mangrove forests in this region had the highest mollusc diversity of all mangrove wetlands in China . In 2008, this region suffered an ECE that was characterized by daily minima below 10°C for more than 3 weeks ( Figure S1b). During this event, the minimum air temperature at the adjacent weather station decreased to 3°C ( Figure S1b). A mean air temperature of 12.6°C was found in January and February of 2008 ( Figure S1b). The number of continuous days with abnormally low and/or freezing temperatures broke records that had been held in the region since the winter of 1954 Liu et al., 2016).  Figure S3). Sediment samples were wet sieved in the field through a 1 mm mesh sieve, after which living molluscs were collected and fixed in 4% formaldehyde. Within 3 days of fixation, specimens were washed with water and transferred to 75% ethanol.

| Sampling and environmental variables
The following local environmental variables were measured immediately after taking biotic samples at each tidal zone within a transect during sampling seasons: water salinity (Sal) was measured using a hand-held refractometer (three to six replicates); vegetation types (VT) were determined by the dominant mangrove species and latitude and longitude (coordinates) were measured using a hand-held GPS. The seasonal average highest atmospheric temperature (HT), seasonal average lowest atmospheric temperature (LT) and seasonal average rainfall (Rain) for each site were based on data from the China Meteorological Data Sharing Service System (http://data. cma.cn/site/index.html). The tidal amplitude (Tide) for each study site was based on data from the National Marine Science Sharing Service Platform (http://mds.nmdis.org.cn). Taxonomic categories of molluscs were then identified using morphological differences (species identity), counted to determine abundance and measured to determine wet weight. In this study, species abundance (abundance) and species identity (species richness) were used to represent compositional features, and wet weight (biomass) was used to represent functional traits. The biomass (body mass) of soil molluscs may be an important trait for determining the resources used, and carbon regulation may be important for revealing the relationships between climatic disturbance and soil biodiversity (Chen, Gu, Liu, Shi, et al., 2021;Liu et al., 2015;Moretti et al., 2017; van der Schatte Olivier et al., 2020).

| Definition of different functional groups
To assess whether different functional groups responded differently to the ECE, we constructed two non-overlapping functional groups (Dethier & Steneck, 1994). We categorized benthic molluscs into two functional groups, ISDFs and MS/Ps, without combining morphospecies of differing feeding types or mobility (Barnes et al., 2019).

| Impact of ECE on the soil mollusc community
To detect the influences of the ECE on the soil mollusc community, we compared the annual species' abundance, species' richness and biomass before (in 2007) and after (in 2008 and 2009) the ECE. To compare these measures, a one-way analysis of variance (ANOVA) was performed in R v3.6.1 (R-Core-Team, 2019). Data were tested for homoscedasticity and normality using Levene's test and the Shapiro-Wilk test, respectively; if normality assumptions were not met, non-parametric equivalent analyses were used (Kruskal-Wallis ANOVA and Mann-Whitney U tests). Bonferroni corrections were used to assess the p-values of the pairwise comparisons between years.

| Effect of the ECE on the spatiotemporal dynamics of the soil mollusc community
Beta diversity has proved to be a useful approach to assess community response to disturbances, highlighting community responses by the loss and gain of taxonomic composition and functional traits (Baselga, 2010;Vanschoenwinkel et al., 2013). Therefore, Jaccard distance and Bray-Curtis distance were used to describe compositional turnover and trait replacement in the mollusc metacommunity. The Jaccard distance was calculated using the "SpadeR" package in R v3.6.1 based on the presence-absence matrix (R-Core-Team, 2019). The Bray-Curtis dissimilarities were calculated using the SpadeR package in R v3.6.1 based on the relative abundance and relative biomass data (R-Core-Team, 2019). To test the effects of ECE on the beta diversity of soil mollusc communities, a non-metric multidimensional scaling (NMDS) ordination was used to investigate the changes in the soil mollusc communities over time from 2007 to 2009. To evaluate the significant differences in mollusc communities between these years, we used the randomization/permutation procedure analysis of similarities (ANOSIM) and permutational multivariate analysis of variance (PERMANOVA) tests (Anderson, 2017).
The global R of ANOSIM is calculated as the difference in betweenyear mean rank similarities; thus, it displays the degree of separation between years (Mo et al., 2018). R = 1 indicates a complete separation of the community, whereas R = 0 means no community separation (Mo et al., 2018). The significance level was assessed using a Bonferroni correction.
To evaluate whether the ECE caused the changes in spatial limitations of the soil mollusc community on our study's biogeographical scale (~350 km) (Hortal et al., 2010), Mantel tests and Spearman's rank correlations were used to assess the distance-decay relationships between geographical distance and community dissimilarity Mo et al., 2018). We calculated the geographical distances using the "gdistance" package in R v3.6.1 (R-Core-Team, 2019) based on coordinate data (i.e. longitude and latitude). In the Mantel test, if the correlation coefficient (Mantel's r) is greater than 0 and the p-value is smaller than 0.05, a positive relationship exists between the community dissimilarity and geographic distance; that is, the soil mollusc communities were significantly spatially limited (Nekola & White, 1999). However, if Mantel's r is smaller than 0 and the p-value is smaller than 0.05, there is a negative relationship between community dissimilarity and geographic distance, that is, the soil mollusc communities were significantly spatially dispersed (Nekola & White, 1999). The absolute value of Mantel's r (Mantel's |r|) reflected the intensity of the spatial limitation (or dispersal) of the community, and indicated whether the ECE suppressed or promoted the spatial turnover of soil molluscs .
To assess the influence of ECEs on temporal turnover (i.e. ecological succession) on the regional-scale mollusc community (i.e. divergent or convergent succession; Leqa & Rejmánek, 1991), Mantel tests and Spearman's rank correlations were used to assess the time-decay relationship (TDRs) between community dissimilarities and temporal distance (days; Guo et al., 2018). A positive relationship (i.e. Mantel's r > 0, p < 0.05) indicates community dissimilarity increases significantly with increasing temporal distance, that is, the community becomes increasingly dissimilar over time and results in divergent succession (Guo et al., 2018;Leqa & Rejmánek, 1991), while a negative relationship (i.e. Mantel's r < 0, p < 0.05) indicates that community dissimilarity decreases significantly with increasing temporal distance, and the community becomes increasingly similar over time, leading to convergent succession (Lee et al., 2017;Leqa & Rejmánek, 1991). The Mantel's |r| values were used to evaluate whether the ECE promoted or suppressed divergent (or convergent) succession in the soil mollusc community.

| Co-occurrence network construction
To detect whether and how the ECE affected the co-occurrence of the soil mollusc community, we calculated the Spearman correlations for pairs of species based on abundance and biomass data. The calculated correlations were used to construct species co-occurrence networks from before (in 2007) and after (in 2008 and 2009) the ECE (Chen, Gu, Liu, Shi, et al., 2021;Mo et al., 2021).
We also constructed sub-networks for the ISDF and MS/P mollusc sub-communities to uncover the difference between different functional groups. Higher network modularization indicates higher community complexity and stability (Grilli et al., 2016). To characterize the recovery strategies of mollusc communities after the ECE (i.e. more complex or simpler), the modularity of networks was measured (Grilli et al., 2016).
To determine how soil molluscs co-exist (i.e. random or nonrandomly aggregated or non-randomly segregated) and to assess whether and how the ECE influenced the co-occurrence pattern, null models were further calculated in EcoSim v7.72 based on presenceabsence, relative abundance and biomass data (Gotelli & Entsminger, 2001). The checkerboard score (C-score) was used as a matrix of mollusc co-occurrence. The observed C-score was calculated and compared with the C-scores calculated for 5000 randomly assembled null models based on the tail probability that the observed (O) index was larger or smaller than expected (E) by chance (Gotelli & Entsminger, 2001;Gotelli & McCabe, 2002). Values of O > E for the C-score indicate negative co-occurrence (Gotelli & Entsminger, 2001), suggesting that soil molluscs are non-randomly segregated, resulting in reduced rates of co-occurrence. Conversely, values of O < E for the C-score reflect a positive co-occurrence (Gotelli & Entsminger, 2001), suggesting that soil molluscs non-randomly aggregate with similar compositional and functional traits. Finally, values of O = E for the C-score suggest that molluscs randomly cooccur. To determine whether and how the non-randomness of the soil mollusc metacommunity varied in response to ECEs, we calculated the standardized effect size (SES) for each mollusc community metric as follows (Gotelli & McCabe, 2002): where Iobs is the observed value, Isim is the mean simulated value and εsim is the standard deviation of the simulated indices. The null hypothesis is that the average SES = 0, and that 95% of the Iobs values fall between 1.95 and −1.95 (Gotelli & McCabe, 2002).

| Effect of the ECE on the ecological factors of soil mollusc communities
To evaluate whether and how the ECE changes the environmental factors that influence the distribution of the soil mollusc community, random forest (RF) machine learning (Breiman, 2001) was performed to assess the major drivers of the mollusc community in R v3.6.1 based on species richness, abundance and biomass data.
In the RF model, the percentage increase in the mean squared error (IncMSE%) was used to assess the importance of environmental variables (Breiman, 2001;Mo et al., 2021). Environmental factors with a high IncMSE% value are the most important factors. To assess the significance of the RF model, 5000 permutations of the response variable were performed within the "A3" package (Mo et al., 2021). The significance of each environmental factor was determined using the "rfPermute" package in R v3.6.1.
To evaluate whether and how ECE changes the environmental and spatial processes of beta diversity in the mollusc community, variation partitioning analysis (VPA) and Mantel tests were used to assess the relative contribution of these ecological processes.
According to the longitude and latitude coordinates, principal coordinates of neighbour matrices (PCNM) analysis were used to generate the regional spatial variables in CANOCO v5 Šmilauer & Lepš, 2014). Based on the longest gradient lengths of detrended correspondence analysis (DCA), either a redundancy analysis (RDA) or canonical correspondence analysis (CCA) was selected to assess the relationship between mollusc communities and local environmental variables/regional spatial factors Šmilauer & Lepš, 2014). The longest gradient lengths were <3 for total, ISDF and MS/P communities, indicating that RDA was suitable for these communities.
Before RDA analysis, forward selection was conducted to select the environmental variables and PCNMs with significant explaining factors (p < 0.05) for further analyses. Then, variation partitioning analysis (VPA) was used to evaluate the relative contribution of environmental filtering and regional dispersal process in shaping mangrove soil mollusc communities with adjusted R 2 coefficients based on RDA in CANOCO v5 Mo et al., 2018;Šmilauer & Lepš, 2014). The relative role of both components was explained by pure spatial factors (S|E), pure environmental factors (E|S) and the combined effects of both space and environment (S ∩ E). In this analysis, the residual proportion represents the unexplained variance. Furthermore, Mantel tests were conducted to identify relationships between mollusc community dissimilarity and environmental spatial distances (Mo et al., 2018).
In the Mantel tests, the first two spatial components (i.e. PCNM1 and PCNM2) and the Euclidean distance of environmental factors were used for analysis.

| The ECE influenced the soil mollusc community
The species richness (Figure 1a indicating that the ECE greatly influenced the regional soil mollusc community. Specifically, the species richness of ISDFs continued to decrease in the 2 years after the ECE (one-way ANOVA, p < 0.01;  Table S1a). These results indicate that the ECE had a great influence on the mollusc community, and that the effects of the ECE on ISDFs were different from those on MS/Ps.

| The ECE affected the spatiotemporal turnover of the soil mollusc community
A positive relationship (Mantel test, Mantel's r > 0, p < 0.01; Figure 3; Table S2) was observed between temporal distance and community dissimilarity of the total, ISDF and MS/P soil molluscs, indicating that divergent succession may be the general temporal dynamic of mangrove soil molluscs. The intensity of the divergent succession in ISDFs was greater than the strength of the divergent succession in MS/Ps, since the value of Mantel's |r| in ISDFs (Mantel test, mean |r| = 0.42) was greater than that in MS/Ps (Mantel test, mean |r| = 0.22; Table S2; one-way ANOVA, F = 33.94, p < 0.001; Figure   F I G U R E 1 The species richness (a, b, c), abundance (d, e, f) and biomass (g, h, i) of the mollusc communities before (2007) and after (2008 and 2009) an extreme cold event (ECE). The "n" indicates the number of samples. p < 0.05 was the significance level in the one-way analysis of variance (ANOVA). "ISDF" indicates the in-faunal shelled deposit feeders, "MS/P" reflects the hard-bodied, mobile scavenger/predators and "All" means the total mollusc community. "*" indicates p < 0.05 and that it passed α o = 0.05 (<0.017) in the Bonferroni test, "ns" indicates p > 0.05 and that it did not pass the level of α o = 0.05 (>0.017) in the Bonferroni test S4a). After the ECE, Mantel's |r| of ISDFs increased by 60% (one-way ANOVA, F = 759.45, p < 0.001; Figure S4b), while in MS/P taxa, it decreased (one-way ANOVA, F = 0.46, p = 0.65; Figure S4c). This reflects that the ECE enhanced the divergent succession of ISDFs, while reducing the divergent succession of MS/Ps. A positive relationship (Mantel test, Mantel's r > 0, p < 0.01; Figure S5; Table S2) was found between geographical distance and community dissimilarity of all ISDF and MS/P taxa, reflecting the spatial differentiation of soil molluscs. The mean Mantel's |r| of ISDFs (0.4) was greater than that of MS/Ps (0.28), reflecting that the strength of spatial differentiation in ISDFs was greater than the strength of spatial differentiation in ISDFs (one-way ANOVA, F = 11.05, p = 0.004; Figure S6a). After the ECE, the mean Mantel's |r| of ISDFs increased significantly (one-way ANOVA, F = 150.82,

F I G U R E 2
The changes in the community composition of the mangrove mollusc over time from 2007 to 2009. Non-metric multidimensional scaling (NMDS) ordination based on the Jaccard dissimilarity of presence/absence data (a, b, c). Non-metric multidimensional scaling (NMDS) ordination based on the Bray-Curtis dissimilarity of relative abundance data (e, d, f). Non-metric multidimensional scaling (NMDS) ordination based on the Bray-Curtis dissimilarity of relative biomass data (g, h, i). "ISDF" indicates the infaunal shelled deposit feeders, "MS/P" reflects the hard-bodied, mobile scavenger/predators and "All" means the total mollusc community. The R values for each pair-wise group are shown in Table S1 of supplementary information p < 0.001; Figure S6b), while that of MS/Ps decreased significantly (one-way ANOVA, F = 87.81, p < 0.001; Figure S6b). This indicates that the spatial differentiation of mangrove soil molluscs was influenced by the ECE; the spatial differentiation of ISDFs was enhanced by the ECE, while the spatial differentiation of MS/P taxa was reduced by the ECE.

| Effects of the ECE on the co-occurrence patterns of soil mollusc communities
The observed C-scores and SESs of mangrove mollusc communities were greater than the null model (null model analysis, p < 0.05, SES > 1.95; Figure S7; Table S3), meaning that a segregated co-occurrence pattern may be a common pattern of the mangrove mollusc community. After the ECE, the SESs of ISDFs were reduced, while those of MS/Ps increased (null model analysis, Figures S7 and S8). Additionally, the modules of the total community network increased greatly after the ECE ( Figure S9; Table   S4), indicating that the network structure of the mollusc community was significantly influenced by the ECE. After the ECE, the network modules of ISDFs continued to increase, while those of MS/Ps decreased in the first year but then increased (Figure 4).
These results suggest the great effects of the ECE on the cooccurrence pattern of the mangrove soil mollusc communities, with different effects of the ECE on the co-occurrence pattern of ISDF and MS/P taxa.

| Effects of the ECE on the driving factors of the soil mollusc community
We found that the hydrology (Tide and Sal), meteorology (LT, HT and Rain) and vegetation (FT) factors significantly affected species richness, abundance and biomass of all mangrove molluscs (RF analysis, p < 0.05; Figure 5). Before the ECE, Sal and LT were the most important variables for the total mollusc community (RF analysis, p < 0.05; Figure 5). After the ECE, these factors were replaced by meteorology (LT, HT and Rain) and vegetation (FT) factors (RF analysis, Figure 5). Among the ISDFs, Sal was the most important variable for predicting species richness, abundance and biomass before the ECE (2007) (p < 0.05; Figure 5), but it was replaced by Rain, FT, Tide and HT after the ECE (2008 and 2009) (p < 0.05; Figure 5). Among MS/P taxa, Sal was the most important variable for predicting the community before the ECE (2007) (p < 0.05; Figure 5). However, after the ECE, HT, LT and FT were the most important factors shaping the F I G U R E 3 The relationships between community dissimilarities and temporal distance before (2007) and after (2008 and 2009) the extreme cold event (ECE) based on the Mantel test and Spearman's rank correlation. The community dissimilarities include the Jaccard dissimilarity of presence/absence data and Bray-Curtis dissimilarities of relative abundance and biomass data MS/P community (p < 0.05; Figure 5). These results indicate that the ECE changed the driving factors of the mangrove mollusc community. Different effects of the ECE on the ISDF and MS/P taxa were determined.
We found that the relationships between the driving factors of the mollusc communities changed after the ECE since Spearman's ρ changed greatly (Spearman correlation analysis, Figure 6; Figure   S10). We observed a clear correlation between the community dissimilarities and driving factors, which were significant on all community scales (i.e. total, ISDFs and MS/Ps; Mantel test, p < 0.05; Figure 5; Figure S10). This reflects spatial and environmental processes that jointly influence the mollusc community. Importantly, the Mantel's r between the mollusc communities and driving factors changed greatly after the ECE ( Figure 5; Figures S10-12). Specifically, among ISDFs, the absolute value of Mantel's r of the spatial factor (PCNM1-2), meteorological factors (LT, HT and Rain) and hydrological factors (Tide and Sal) were reduced in the 2 years after the ECE (one-way ANOVA, p < 0.05; Figure S12). However, Mantel's r of the spatial, temporal (PCNM1-2) and hydrological factors (Tide and Sal) in MS/Ps continued to decrease after the ECE (one-way ANOVA, p < 0.05; Figure S12).
Additionally, both the purely spatial (S/E) and environmental (E/S) variation significantly explained the total ISDF and MS/P communities (VPA analysis, p < 0.01; Table S5; Figures S13 and   14). The explained proportion of purely spatial (S/E) variation was greater than that of the purely environmental (E/S) variation (Table   S5), indicating that the spatial process was more important than environmental filtering in shaping the mangrove soil mollusc community. Importantly, the relative importance of these processes was broken by the ECE. For instance, after the ECE in 2008, the explained proportion of purely environmental (E/S) variation was greater than that of the purely spatial (S/E) variation in the ISDF and MS/P groups (Table S5). In contrast, in ISDFs, the explained proportion of purely spatial (S/E) and environmental (E/S) variation decreased after the ECE in 2008 (Table S5), whereas in MS/Ps, the explained proportion of purely spatial (S/E) variation decreased, while that of purely environmental (E/S) variation increased after the ECE in 2008 (Table S5).

F I G U R E 4
The network modules of the ISDF and MS/P communities before (2007) and after (2008 and 2009) the extreme cold event (ECE) based on the relative abundance and biomass data. The nodes represent species; the lines between nodes indicate the Spearman correlation between species. Different colours indicate different modules. "ISDF" indicates the in-faunal shelled deposit feeders, and "MS/P" reflects the hard-bodied mobile scavengers/predators reduce nutrient cycling, pollutant burial and mobilization and sediment stability (Lavelle et al., 2006;Lundquist et al., 2018).
Declines in MS/Ps would reduce nutrient cycling and the redistribution of energy and affect disease transmission associated with decomposition (Cannicci et al., 2008;Decaens, 2010;Lundquist et al., 2018). Additionally, declines in biomass would result in a decrease in many ecosystem services, such as food production and carbon sequestration (Chen, Gu, Liu, Shi, et al., 2021).
Decreasing marine food yield may lead to human malnutrition in some food-insecure regions (Chen, Gu, Liu, Shi, et al., 2021).
The carbon shell of soil molluscs, as long-term carbon storers, would decrease after ECEs and influence the ecosystem services of climate warming alleviation (Barnes et al., 2019). An increase in abundance reflects that ECEs might cause further impacts on this ecosystem's functions and services through species interactions (Bairey et al., 2016).

| Different functional groups may use different ecological strategies in response to extreme climate changes
Both theoretical and empirical studies suggest that the spatiotemporal dynamics of community turnover are often the key mechanism affecting the self-recovery of the community Dethier & Steneck, 1994;Guo et al., 2018;Hortal et al., 2010).
The ECE increased the divergent succession and spatial differentia- The factor importance plots show the permutation-based importance in terms of the percentage increase in the mean squared error (% IncMSE) in the RF. Sal = water salinity, FT = vegetation type, LT = seasonal average lowest atmospheric temperature, HT = seasonal average highest atmospheric temperature, Rain = seasonal average rainfall and Tide = tidal amplitude. "ISDF" indicates the in-faunal shelled deposit feeders, and "MS/P" reflects the hard-bodied mobile scavenger/predators. "*" indicates p < 0.05, "ns" indicates p > 0.05 2018; Voigt et al., 2007). Similarly, different functional groups (algal and herbivory) on Jamaica's Discovery Bay reef have different spatiotemporal dynamics in response to hurricanes and mass mortality (Dethier & Steneck, 1994). Different dynamics of spatiotemporal turnover lead to distinct patterns of community structure. Therefore, this may be one of the explanations for the different co-occurrence patterns in ISDFs and MS/Ps. Additionally, increased differences in spatiotemporal turnover are expected to lead to an increase in the spatiotemporal heterogeneity of the community, thus improving the spatiotemporal stability and complexity of the community (Grilli et al., 2016;Yuan et al., 2021). This might be why we observed an increase in the network modules of the ISDF community in the 2 years after the ECE, while a decrease was observed in the MS/P community. Therefore, climate extremes may enhance the stability and complexity of the soil mollusc network and have a different effect on ISDFs and MS/Ps (Grilli et al., 2016;Yuan et al., 2021). Our results reinforce the view that communities with different properties (e.g. differently functional groups) might respond to disturbances, such as extreme climate events, via different mechanisms of spatiotemporal turnover.
Local (e.g. environmental filtering) and regional (e.g. dispersal limitation) ecological processes are two important factors that influence the geographical distribution of soil mollusc communities (Chase et al., 2011;Gotelli & McCabe, 2002). We found that dispersal limitations and environmental filtering jointly affected the mangrove soil mollusc community, while the contribution of dispersal limitation was more important than that of environmental filtering. Notably, the balance of these ecological processes in the mangrove soil mollusc communities was broken by the ECE. For example, meteorological factors (LT, HT and Rain) and vegetation type (FT) replaced Sal as the most important variable for shaping soil mollusc communities after the ECE. Under normal circumstances, salt tolerance may be the main factor that influences intertidal species since the larval dispersal and genetic exchange of intertidal species can be influenced by changes in Sal caused by human activities, such as freshwater discharge (Dong et al., 2016). Climate changes may induce other tolerances in species, such as climate change-linked stress tolerance, which could influence their fitness and chances of survival (Yusefi et al., 2021). Furthermore, climate changes will probably impact coastal species through declines in primary productivity . Additionally, the initiation of dispersal (emigration) by a species, its subsequent movement (transfer) and its settlement decision (immigration) are influenced by local conditions. Therefore, climate changes, such as an ECE, may affect the dispersal ability of a community (Travis et al., 2013). This could explain why we observed a decrease in the dispersal limitation of mangrove soil molluscs.
We found different changes in these ecological processes (environmental filtering and dispersal limitation) between ISDFs and MS/ Ps after the ECE. On the one hand, we found that Rain and Tide were the most important factors that shaped the ISDF community, while these factors were not important for MS/Ps. On the other hand, the explained proportion of purely environmental (E/S) variation in the F I G U R E 6 The relationships between community dissimilarities and environmental spatial distance in 2007, 2008 and 2009 based on the Mantel test and Spearman's rank correlation. Pairwise comparisons of environmental and spatial factors are shown in the upper right, with a colour gradient representing Spearman's correlation coefficients. Mollusc community composition was correlated with each environmental or spatial factor by Mantel tests. The community dissimilarities include the Jaccard dissimilarity of presence/absence data and the Bray-Curtis dissimilarity of relative abundance/biomass data. The environmental distance is reflected by Euclidean distance. The first two principal coordinates of neighbouring matrices (PCNM1 and PCNM2) were used to indicate spatial components. The line width represents Mantel's r statistic for the corresponding correlation, and line colour indicates that significance was tested based on 999 permutations. S1 = spatial factor 1 (PCNM1), S2 = spatial factor 2 (PCNM2), Time = temporal distance, Sal = water salinity, LT = seasonal average lowest atmospheric temperature, HT =seasonal average highest atmospheric temperature, Rain = seasonal average rainfall and Tide = tidal amplitude MS/P community increased after the ECE in 2008, while the explained proportion of E/S in ISDFs decreased. The reasons for these differences can be explained by the fact that different functional groups (ISDFs and MS/Ps) may differ in many fundamental features, such as environmental tolerance and dispersal characteristics. The divergent responses to climate change by trait-specific groups could lead to changes in the composition of ecosystems, calling the resilience and sustainability of various ecosystems into question (Fei et al., 2017;Voigt et al., 2007).

| Implications for soil conservation biogeography
Soil biological conservation aims to maintain soil ecosystem functioning and services. The contribution of each soil species and taxon to soil ecosystem functioning relies on its functional properties, such as herbivory and nutrient cycling. Soil fauna are acknowledged as key components of soil health and soil ecosystem services (Delgado-Baquerizo et al., 2020;Lavelle et al., 2006). Environmental changes, such as extreme climate events, can negatively affect coastal soil management because changes in the structure and function of soil fauna jeopardize soil systems, such as the southern China mangrove soil, where the communities are already influenced by and responsive to anthropogenic-driven abiotic regimes (Boucek & Rehage, 2014;O'Connor et al., 2020). Based on the regional coastal soil mollusc community, a novel contribution of our study is the finding that different functional groups invoke different spatiotemporal turnover models and ecological processes (local environmental filtering and dispersal limitation) to shape their unique co-occurrence patterns in response to ECEs. Methods of dividing different functional groups allow researchers to identify the relative vulnerabilities of member functional groups and, more importantly, the vulnerabilities of their traits to specific stressors associated with environmental changes, such as ECEs. Identifying a quantifiable link between environmental stressors and different functional meta-communities responses greatly improves our predictive power concerning how an ECE may affect the metacommunity function impacted by climate change, potentially providing a powerful tool in soil biodiversity and biogeography conservation. Thus, under environmental changes, such as extreme climate events, a major challenge is to conserve various functional groups while preserving biodiversity, ecological functions and associated ecosystem services. Our results highlight the importance of identifying the responses of different functional groups to climate disturbances. We recommend the use of functional groups to develop indicators for soil conservation biogeography.

ACK N OWLED G EM ENTS
We thank Xiaofang Shi for her assistance during the revision.
We also thank the reviewers for their general and specific com-

CO N FLCIT O F I NTE R E S T
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. All authors agreed to the submission of this manuscript. The material is original research, has not been previously published and has not been submitted for publication elsewhere while under consideration.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13467.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data are available from the Dryad Digital Repository: https://doi. org/10.5061/dryad.jsxks n0bh.