Recurring heavy rainfall resulting in degraded-upgraded phases in soil microbial networks that are reflected in soil functioning

• We assess the recovery of microbial networks underneath crust to repeated rainfall. • The network fragmentation after the second heavy rain was milder than at the first one. • Cohesive networks were related to high enzyme activity involved in C, N, and P cycles. • Loose networks were related to high Ca, K, Mg, NH4 and organic N. • The network in dry-crusted soils collapsed after the second heavy rain. Biological soil crusts (BSCs) are an important multi-trophic component of arid ecosystems in the Mediterranean region. In a mesocosm experiment, the authors investigated how the network of interactions among the members of the soil microbial communities in four types of soil sample responded when soils were exposed to two simulated extreme rain events. The four types of soil samples were: covered by Cladonia rangiformis and previously hydrated (+BSC+H), covered by C. rangiformis and dried (+BSC−H), uncovered and hydrated (−BSC+H), uncovered and dried (−BSC−H). Network analysis was based on the co-occurrence patterns of microbes; microbes were assessed by the phospholipid fatty acids analysis. The authors further explored the relations between networks’ metrics and soil functions denoted by enzymatic activity and soil chemical variables. All networks exhibited Small world properties, moderate values of clustering coefficient and eigen centrality, indicating the lack of hub nodes. The networks in −BSC−H soils appeared coherent during the pre-rain phases and they became modular after rains, while those in +BSC−H soils kept their connectivity till the second rain but this then collapsed. The network metrics that were indicative of cohesive networks tended to be related to enzyme activity while those that characterized the loose networks were related to Ca, K, Mg, NH4+ and organic N. In all mesocosms except for +BSC−H, networks’ fragmentation after the second heavy rain was milder than after the first one, supporting the idea of community acclimatization. The response of microbial networks to heavy rains was characterized by the tendency to exhibit degradation-reconstruction phases. The network collapse in the crusted only mesocosms showed that the communities beneath crusts in arid areas were extremely vulnerable to recurring heavy rain events.

watered or dry conditions. Changes in environmental variables condition the structure, dynamics, and interactions between living organisms as well as the functioning of ecosystems both directly, by affecting the rates of related processes, and indirectly, by affecting diversity (Hong et al., 2021). Stamou et al. (2019) stated that the structure of soil microbial communities is constantly in a state of imbalance and therefore particularly prone to disturbances. Rillig et al. (2015) argued that the response of biocommunities to disturbances is driven by the idiosyncratic responses of the individuals and the energy costs associated with these responses. In particular, the mortality and growth rates of different strains vary (Gibbons et al., 2016), hence the disturbance usually exerts strong selective pressure on soil microbiota (Norris et al., 2002;Azarbad et al., 2016). Concerning the impacts of recurrent disturbances, Philippot et al. (2008) suggested a kind of community acclimatization in the sense that the response of a microbial community to a severe disturbance is often moderated by its previous exposure to a minor disturbance of the same or even different type (Rillig et al., 2015;Konstantinou et al., 2019). To the best of our knowledge, the effects of various disturbance events (torrential rainfall causing severe leaching, extreme drought and high temperature) on various structural aspects of soil communities in arid and semi-arid formations have been extensively studied (Caruso and Migliorini, 2006;Bérard et al., 2011;Sardans and Penuelas, 2013). However, the effects of recurring severe events, such as potent leaching, on the interactions between the members of a microbial community and the relation between these changes and soil functioning are missing from the literature.
In this work, the identification of the interaction patterns between microbiota was based on their co-occurrence and was studied by using network analysis. With this method the authors calculated metrics referring to either the scale of the whole network (connectivity, clustering coefficient and length of shorter path) or critical areas within the network (fragmentation, Small world coefficient and eigenvector centrality). A merit of network analysis is that it reveals global and local aspects of the interactions' configuration simultaneously (Glo-cal analysis; Stamou and Papatheodorou, 2016). This Glo-cal approach was applied successfully to assess ecosystem connectivity in perturbated soil mesocosms (Stamou et al., 2019), to describe the associative patterns of soil nematode communities in reclaimed landfills , and to depict changes in microbial co-occurrence patterns in cultivated areas inoculated with Plant Growth Rhizobacteria (Papatheodorou et al., 2021). Further, such an approach helped us to analyze the interplay between the deterministic and stochastic processes that are thought to modulate the structure of a community and control its resistance and resilience . However, the use of network analysis to examine the stochasticity-deter-minism balance in the regulation of communities' structure and how this is related to their ability to overcome disturbances, is a new field of research. The niche theory assumes that the structure of a community is driven by genetically predetermined rates of change related to demographic traits (fertility, mortality, length of life cycle), as well as intra-and inter-specific competition, predation, synergies etc. These traits constitute part of organisms' identity and guide their response to changing environmental variables (habitat filtering) (Zhou et al., 2014;Vellend et al., 2014). In contrast, the theory of neutrality against organisms' identity postulates that the structure of a community is driven by the intrinsic stochasticity of demographic traits (ecological drift), speciation, extinction, disturbance outbreaks and lack of dispersal limitation (Vellend et al., 2014;Zhou and Ning, 2017). Recently, it has been accepted that deterministic and stochastic processes operate simultaneously, but the contribution of each of them varies over time (Vellend et al., 2014;Maignien et al., 2014), along a deterministic-stochastic continuum . The prevalence of stochastic processes is related to modular architecture (Zhou et al., 2014). Moreover, the resistance and resilience of a network is related to its Small-world properties. In Smallworld networks the length of the shortest path is equal on average to that of random counterparts, while the clustering coefficient is much higher compared to the corresponding coefficient in random networks (Humphries and Gurney, 2008). Further, when the Small world properties are combined with modular architecture, the result is robustness against the effects of random disturbances (Proulx et al., 2005;Tobor-Kaplon, 2006).
In a previous work , the authors studied the dynamics of soil microbial biomass and the composition of microbial communities developed under experimentally created water regimes and crust cover and which were then further subjected to two episodes of heavy rainfall. That study showed that the impact of the second heavy rain on community composition was in most cases milder than that of the first one, a fact that was attributed to community acclimatization. In the present work, the authors used the same data set to test for a similar response in terms of interactions between the members of the microbial communities. To this end, the network of interactions between the microbial biomarkers, as described by the concentration of specific phospholipid fatty acids (PLFAs) present in cellular membranes, was studied based on the following hypotheses: H1: The first severe rain event would result in the acclimatization of the microbial communities while the extent of acclimation would depend on crust cover and previous watering of the crust-soil system. H2: More resistant (extent of resistance to change) networks would be reformed after the second rain episode and this would be reflected in the topological characteristics of the networks (fragmentation, clustering coefficient, Small worldness).
H3: The differences in the network metrics induced by the experimental modification of the soil abiotic conditions would shift the balance in microbial communities' regulation between stochasticity and determinism.
H4: The differences in the network metrics induced by the experimental modification of soil abiotic conditions would be related to soil functionality as this described by enzyme activities and soil chemistry.

The field site
Toward the end of a prolonged summer drought at September 2016, soil was collected from a semi-arid area in Halkidiki (Greece). A complete description of the site and experimental design has been given in a previous work . In brief, the shrubby vegetation was dominated by ericaceous shrubs, while the open sites between shrubs were occupied by Cladonia rangiformis biocrust. Data on climate, soil texture, and pH are summarized in Table 1.

Experimental design
The collected soil was transferred to laboratory and used for the preparation of 64 soil mesocosms (plastic pots 7 cm in diameter and 4 cm in depth). The experimental setup included mesocosms either covered (+BSC) or uncovered (-BSC) by biocrust, and either non-hydrated (-H) or regularly hydrated (+H) every three days with 4.5 mm distilled water according to a high-rainfall frequency scenario for drylands (Solomon et al., 2007). The setup complied with a full 2×2 full factorial design with four levels:+BSC+H,+BSC-H, -BSC+ H, -BSC-H. Mesocosms were maintained in a glasshouse for 48 days at a temperature ranging from 15 ± 0.5°C at night to 25 ± 1.5°C in the daytime.
After 21 days a destructive first sampling was undertaken involving the harvesting of four replicate mesocosms from each treatment. Then, the remaining mesocosms (48) were exposed to a simulated heavy rain event. Each mesocosm was sprayed evenly with 80 mm water corresponding to 2/3 of the maximum monthly precipitation recorded in the field during December and January. This was followed by a second destructive sampling (four replicates from each treatment) 3 days later. The remaining 32 mesocosms continued to be incubated under the initial conditions until day 45, when a third destructive sampling was performed, again with 4 replicate mesocosms from each treatment, followed by a second heavy rain event as above. Three days later the fourth and final destructive harvest took place. The two episodes of heavy rain were considered "disturbances" while the intervals before and between rain events were termed as "conditioning phases". A graphical description of the experimental design is presented in Fig. 1. At wet soil subsamples were used. These analyses took place within two weeks after sampling, in subsamples that were stored at 4°C. The remaining chemical variables were determined in air-dried soil subsamples.

Enzymatic activity
At each sampling, the authors assessed the activity of five enzymes. These are urease (UREA) and N-acetyl-glycosaminidase (NAG) involved in the N cycle, b-glucosidase (BG) and phenol oxidase (PPO) related to the C cycle and acid phosphomoesterase (AP) involved in the P cycle. Enzyme activity (except urease activity) was determined according to the procedures of Allison and Jastrow (2006), modified in order to be applicable for 96-well microplates. Urease activity was determined according to the method of Sinsabaugh et al. (2000). Detailed description of the procedures is presented in the study by Papatheodorou et al. (2020).

Chemical analyses
Soil organic C (SOC) was determined by a wet oxidation titration procedure using an acid dichromate system (Allen, 1974). Soil organic N (SON) was measured by the Kjeldahl method. Ammonium-N and nitrate-N were determined in 2 M KCl extracts (1:10 soil dry weight solution) by distillation and subsequent titration (Allen, 1974). For extractable P, the method of Olsen et al. (1954) was used, as specified by Allen (1974). Mg, Ca and K were estimated in soil extracts by atomic absorption spectrophotometer (Perkin-Elmer 2380).

Statistical analysis
A network analysis was employed to configure the associative spatial relationships among individual PLFAs. The nodes represented individual PLFAs and ties represented the strength of their co-occurrence, as estimated by the significant Pearson correlation coefficients between them. Cooccurrence valued matrices were constructed separately for each sampling and each treatment. In cases of significant correlations, the respective matrix elements were set to the values of the corresponding coefficients, otherwise to zero. Then, the matrices were analyzed employing the network analysis software UCINET 6 (Borgatti et al., 1999).
To describe the global aspects of the networks, certain parameters were estimated concerning a) network cohesion, b) clustering, i.e., the tendency of nodes to form clusters and c) network centrality. Cohesion assesses the connectivity of a network (O'Malley and Marsden, 2008;Borrett et al., 2014), clustering indicates the tendency of nodes to bundle (Borrett et al., 2014), and centrality assesses how influential is the position of a given node within the network (Rampelotto et al., 2014). Connectivity was estimated as average degree (Deg) and length of the shortest path (L), clustering by using the clustering coefficient (CC), network centrality as % eigenvector centrality (EC) and fragmentation (Fr) as the proportion of pairs of nodes that cannot reach each other. Then the authors partitioned the nodes into a number of factions. For each faction the number of internal ties was greater than the external. The choice of the best number of factions was based on the values of Q modularity index accounting for the fraction of internal ties in each faction minus the expected fraction if they were distributed at random but with the same degree. Finally, to check whether our networks exhibited characteristics of Small wordlness, the authors estimated the Small worldness index (SW) proposed by Humphries and Gurney (2008). A realworld network G is termed Small-world if the shortest path (L g ) estimated for the network G is more or less equal to the shortest path (L rand ) estimated for an equivalent random graph (L g ≈L rand ) and if the clustering coefficient for G (CC g ) is higher than that estimated for an equivalent random graph (CC g > > CC rand ), i.e., if SW = (CC g /CC rand )/(L g /L rand ) > 1. To check departure from randomness in the values of the short-est path and clustering coefficient for each experimental sample, 20 random matrices with fixed row and column marginals (Fix-Fix method; Gotelli and Entsminger, 2001) were created and analyzed in a similar way. The outputs of the experimental networks were tested against those provided by the random networks, using the z-test.
Finally, to investigate possible correlations of the network architectural attributes (Deg, CC, EC, Fr, SW and L), with certain functional surrogates (BG, NAG, AP, UREA, PPO, SOC, SON, K, Ca, Mg, NH 4 + , NO 3 − ), a Principal Component and Classification analysis (PCCA) was performed. In PCCA analysis the ordination of cases is based not only on the variables used in the classical PCA (Principal Component Analysis) but also on other supplementary variables. PCCA allows the analysis of data for variables that are heterogeneous with respect to their means or with respect to both their means and standard deviations, by providing an option to analyze covariance matrices as well as correlation matrices. For example, the samples collected before and after the first rain were classified according to the values of their network parameters (L, CC, EC, SW, FR etc.) and the values of other supplementary variables which were the soil enzymes (UREA, NAG, BG, AP, PPO) and other soil chemical variables (Ca, Mg, SOC, SON etc). More specifically, the two data sets collected before and after the first rain in crusted and uncrusted soils were analyzed together, and the same was done for the two data sets collected before and after the second rain. The PCCA outputs were projected on the first two axes of the respective biplot. Variables and sites were loaded according to factor coordinates and factor scores, respectively, both based on correlations. Then, variables and sites were assigned to three clusters by applying a K-means cluster analysis. A one-way ANOVA was applied to check whether the differences in the average loadings of the different clusters on both axes were statistically significant. These analyses were conducted with the STATISTICA 7 package Statsoft.

Results
According to Table 2, all networks exhibited Small-world properties, i.e., they all showed higher clustering coefficients and smaller shortest paths compared to those of networks from random fix-fix samples. Indeed, Small-worldness index (SW) was higher than unity in all networks while in only three networks were the values of the clustering coefficient below 0.65 (mainly at networks from the second sampling). In addition, in 13 out of 16 networks the percent centrality was low to moderate (< 58%) indicating the lack of hub nodes, i.e., few nodes had many connections with other nodes of the network. The experimental manipulation applied during the initial conditioning period of 21 days resulted in differences in networks' architecture between treatments (Fig. 2); regular watering independent of the crust cover transformed the networks from more connected to fragmented. Compared to the remaining treatments the networks in the control soils (-BSC-H; Fig. 2A) exhibited the highest average degree and clustering coefficient and the lowest value of shortest path, while the SW coefficient was slightly higher than unity. All biomarkers constituted a single, well-connected Small world. The extremely low eigenvector centrality value (1.5) indicated lack of nodes with enhanced centrality and with even distribution of ties' strength over nodes. Judging from the values of average degree, fragmentation, and coefficient of smallworldness in +BSC+H soils, (Table 2; Fig. 2D), both crust cover and regular watering acted as minor disturbances transforming the network architecture from non-modular to modular. In this network many biomarkers, mainly representative of cyanobacteria, protozoan and fungi, remained isolated. The centrality value was moderate (< 40%) but the distribution of influence between nodes was not even. Actually, the CV value for the network centrality in +BSC+H was 0.88 (Table S1). By contrast, nodes with higher eigenvector centrality were included in the principal faction (CV = 0.04). The authors noted that compared to the entire networks, in all cases lower CV values of eigenvector centrality were recorded in the main subnets. Networks' architecture was extremely different at the second sampling (Fig. 3). Three days after the first heavy rain event, three out of the four networks were driven to fragmentation and most of the biomarkers remained isolated (Figs. 3A, B, D). The most connected network was that in +BSC-H soils (Fig. 3C). Obviously, the presence of the crust in non-watered soils (+BSC-Η) operated favorably; the average degree was maintained at elevated levels (Table 2 and Fig. 3C). In addition, the independent crust effect maintained and even improved all architectural metrics shown by the respective network prior to the rain event.
Three weeks after the first rain event (at the 3rd sampling) there was a recovery in networks' connectivity (Fig. 4). This was stronger in the networks from watered (+BSC+Η and -BSC+Η; Figs. 4D, B) as well as in control soils (Fig. 4A). In all cases, the centrality remained below 20% (Table S1), while the high values of SW (higher than those in the respective random matrices) showed that all four networks exhibited small world properties (Table 2).
In comparison to the first, the second rain event induced milder negative changes in most networks (4th sampling; Fig. 5). In mesocosms other than the network in +BSC-H soils, which collapsed (Fig. 5C), the heavy rain episodes resulted in connected networks albeit with some isolated nodes. In the uncrusted mesocosms (-BSC-H and -BSC+H; Figs. 5A, B) low connectivity and moderate clustering coefficient and centrality values were recorded, while the estimated SW coefficients were greater than 3 (Tables 2 and S1). The biomarkers were distributed into two main sub networks. On the other hand, in the crusted-watered mesocosms (+BSC+ H; Fig. 5D), high density and clustering coefficient, and low centrality and short path were recorded while SW was lower than two. Practically, the biomarkers were positioned in one well connected sub network with Small world properties. By comparing the changes of network topological characteristics across time, differences between crusted and uncrusted soils were revealed. The networks from uncrusted soils (-BSC+H and -BSC+H) exhibited approximately similar responses to the two rain events. After the first major disturbance (heavy rain) these networks became fragmented; the connectivity and the clustering coefficient reduced while the shortest path increased. However, 21 days (after the first rain event) was an adequate period for networks' recovery, as was shown by their topological metrics that returned to values recorded before the first rain (e.g., Fr 0.07 at first and third sampling vs. 0.73 in -BSC-H and 0.84 in -BSC+H soils at second sampling; CC 0.97 at first and third sampling vs. 0.59 at second sampling in -BSC-H soils). The second rain event resulted in fragmentation although to a lesser extent than was induced by the first event (Fr 0.32 after the second rain vs. 0.73 after the first one in -BSC-H soils, 0.62 vs. 0.84 in -BSC+H soils). However, in crusted soils the response of the networks depended on the water regime. For +BSC+H the first rain event caused remarkable change in the topological features but after the second conditioning phase the values of the parameters improved (connectivity and cluster-ing coefficient) and remained almost at this level even after the second rain. Finally, the +BSC-H networks showed an idiosyncratic response. The first disturbance resulted in improved connectivity and clustering coefficient and decrease of shortest path leading to a more coherent Fig. 3 The networks of interactions among the members of the soil microbial community lying underneath biocrust after the first heavy rain event in uncovered by crust-dry (-BSC-H; Fig. 3A), uncovered-watered (-BSC+H; Fig. 3B), crust covered-dry (+BSC-H; Fig. 3C) and crust covered-watered (+BSC+H; Fig. 3D) soils. Nodes of different colors correspond to different subnetworks. The bigger the size of a node the higher its eigen centrality. Blue and red lines denote ties that are internal and external to factions, respectively.

Fig. 4
The networks of interactions among the members of the soil microbial community lying underneath biocrust before the second heavy rain event in uncovered by crust-dry (-BSC-H; Fig. 4A), uncovered-watered (-BSC+H; Fig. 4B), crust covered-dry (+BSC-H; Fig. 4C) and crust covered-watered (+BSC+H; Fig. 4D) soils. Nodes of different colors correspond to different subnetworks. The bigger the size of a node the higher its eigen centrality. Blue and red lines denote ties that are internal and external to factions respectively.  network and this continued to exist during the second conditioning phase. However, after the second disturbance the network architecture collapsed.
The PCCA biplots (Fig. 6) depict the associations between certain functional surrogates with metrics referring to the architecture of networks. The first axes account for 48.75% to 58.41% of data variability, while the corresponding percentages for the second axes range from 26.18% to 34.83%. In all cases the variables and sites are classified into three clusters. A one-way ANOVA showed that with respect to both of the PCCA axes, the average loadings differed significantly among the three clusters. The first PCCA axes depict the distinction between samples collected before and after the rains, with the exemption of the samples collected from the crusted areas pre-and post-rain (Fig. 6B). The left part of the graphs is mainly occupied by the samples taken before the rain. These samples are all highly cohesive in terms of network architecture, so the Deg (average degree) and CC (average neighborhood) metrics tend to be classified relative to them (Fig. 6). Most enzymes as well as NO 3 − and SOC showed a tendency to be associated positively with the average degree (Deg) and/or the average clustering coefficient (CC). The networks in the postrain samples showed network loosening and tend to be classified toward the right part of the first axes together with metrics such as EC, FR, SW as well as NH 4 + , Mg, K, Ca, and SON. Also, in samples where heavy rain was applied for first time (Fig. 6A, B) there was a tendency for L (shortest path) to be related to P. To sum up, in pre-rain coherent networks most enzymes were related to Deg and CC while in post-rain loose networks the metrics EC, FR and SW were related mostly to available nutrients.

Temporal changes in networks' architecture in different treatments
The non-modular architecture of the microbial network after the first conditioning phase in -BSC-H soils (control samples) seems to reflect the spatial homogeneity in habitat attributes caused by the summer water deficit in uncrusted Mediterranean soils. Probably, the exploitation of limited soil resources during this period requires the synergy of different coexisting microbial groups that produce a variety of extracellular enzymes (Papatheodorou et al., 2021). The repeated watering in -BSC+H soils acted as a selective force, locally favoring certain microbial groups to the detriment of others and so resulting in the transformation of the architecture from non-modular to modular. In these soils the small properties of the network were retained, as watering could change habitat filtering by affecting soil chemistry through water effects on the distribution of ions (Aslam et al., 2016). In addition, repeated watering maintained the aqueous corridors through the finer channels of the soil matrix, thereby enhancing the aqueous communication between the habitat's microsites (Štovicek et al., 2017) which initially facilitated the displacement of microbiota (Schütz et al., 2009). In the case of crust covered-watered soils (+BSC+H), the diffusion of the secondary metabolites produced by C. rangiformis (Akpinar et al., 2009) exerted negative effects on some microbes , while releasing others from competition. Further watering and subsequent leaching due to sandy soil texture, possibly led to increased habitat heterogeneity (including for the distribution of O 2 ), which was reflected in reduced network connectivity (the lowest among treatments) and increased Fr, SW and ΕC in these mesocosms. Cyanobacteria, sulfate bacteria, algae, protists, and other micro eukaryotes were severely impacted, losing their interactions with others and remaining isolated. The response of cyanobacteria to water addition is bibliographically ambiguous, since unlike in the soil layers they seem to be favored under conditions of water saturation in the crust layer (Abed et al., 2010). In +BSC-H soils, the nondiffusion of the crust's secondary antimicrobial metabolites, due to the lack of watering, resulted in gentle changes in the network architectural features; the network kept its connectivity and its Small world properties. The remarkable changes caused by the first rain event in the composition of the microbial community , were accompanied by rapid changes in the interactions between microbes that led to the fragmentation of networks (although mesocosm +BSC-H was an exception; Fig. 3C). According to Papatheodorou et al. (2020), the changes caused by the first heavy rainfall in the composition of the microbial community were greater in the mesocosms +BSC+H and -BSC-H (lower resistance), moderate in the +BSC-H and smaller in the -BSC+H (higher resistance). The respective effects on network architecture were different. Compared to the pre-rain networks, remarkable changes in architectural metrics leading to disconnection were recorded in -BSC-H and -BSC+H mesocosms, moderate changes in -BSC+H ones while the network in +BSC-H changed slightly. The architecture in +BSC-H mesocosm was characterized by the highest connectivity (fragmentation was zero), and clustering coefficient and the lower Small-worldness and centrality.
Obviously, the crust cover in dry soils favored network cohesion when subjected to heavy rainfall, possibly mitigating increases in infiltration rate and reducing the surface soil moisture, thus avoiding osmotic shock. As neither microbial biomass, nor habitat traits (water content, Ca and P exempted) differed in samplings before and after the first rain event , the recorded adverse effects on the networks in the cases of the other three treatments were likely due to the osmotic shock associated with the sharp increase in moisture after heavy rain, changes on soil O 2 availability and the transfer of microorganisms via leaching (Kieft et al., 1987;Bossio and Scow, 1998;Unger et al., 2009). In addition, in non-crusted soils several biomarkers, mainly mono-unsaturated longchain fatty acid, remained isolated from the network or participated in secondary dyads and triads. Mono-unsaturated and long-chain fatty acids represent guilds strongly dependent on aerobic conditions and habitat availability; therefore, they were more sensitive to the restrictions imposed by the filling of soil pores with water (Unger et al., 2009).
Based on our results, three weeks after the first heavy rain episode the architecture of most networks recovered impressively. Till now the resilience of soil microbial communities, i.e. their ability to return to an equilibrium state similar to that prior to disturbance, had been studied in relation to population numbers, biomass and activity. Hueso et al. (2011) recorded rapid recovery of microbial populations after soil rewetting in a degraded semi-arid region in SE Spain, while Fuchslueger et al. (2016) stated that the microbial biomass in central European mountain grasslands responded rapidly and in a similar way to natural drought and rewetting. In arid and semi-arid regions microbial communities are well adapted and resilient to drought-rewetting effects (Williams, 2007;Tecon and Or, 2017) since longterm drought history leaves a legacy effect on microbial functions (Leizeaga et al., 2022). All these studies referred to the recovery of microbial biomass and activity, while the available information concerning the recovery of the network of interactions among soil microbiota from arid areas following heavy rain events is extremely rare. As is the case for the present study, the results of Kauper et al. (2021) reported that the soil methanotrophic network in paddy soils became more complex and increased its connectivity during recovery from desiccation-rewetting in a period > 27 days. Further, when another type of disturbance such as that caused by the foraging activity of burrowing mammals was examined, the bacterial network increased in connectedness after bioperturbation (> 12 months) (Eldridge et al., 2015).
A variety of causes can be involved in the recovery of network architecture. In dry mesocosms (-BSC-H and +BSC-H), rainfall could provide the reactivated microbes with solutes produced by the rapid decomposition of accumulated labile organic matter and dead microbial biomass. In the hydrated counterparts (-BSC+H and +BSC+H), the recovery in the interactions was not as pronounced as in dry soils. So, watering diminished the resilience of the microbial community after a heavy rain probably due to increased infiltration rate followed by the removal of ions away from the pots or due to a reduction in O 2 concentration which limited habitat availability. Specifically, the recovery in crusted and watered soils (+BSC+H) was likely associated also with already selected guilds, able to utilize the leachates from the crust, while the hydration of the crust increased its photosynthetic capacity and the ability for N 2 fixation that could provide necessary nutrients to underlying communities.
Compared to the first rain episode, the second one caused similar, albeit milder effects on network architecture, in all treatments except +BSC-H, confirming the H2 hypothesis that the resistance of microbial networks to change was higher after the second rain event compared to the first one. In uncrusted mesocosms and regardless of the hydration regime, the second rain event led to the formation of modular architectures and to some isolated biomarkers, while maintaining Small world properties. To explain milder effects, the authors adopted the argument of Fierer et al. (2003) suggesting that rapid changes in soil water status may exert selective action in favor of microbes capable of synthesizing solutes that enhance their osmoregulatory capabilities. Indeed, inspecting Iso/Anteiso values, Papatheodorou et al. (2020) concluded that in part the selection of water-shockresistant guilds had already been achieved during the first heavy rain. This explanation agrees with Bouskill et al. (2016) stating that some of the drought-induced negative effects are weakened by prior exposure to a short-term drought, suggesting that acclimation may occur despite the lack of longer-term drought history. Leizeaga et al. (2022) found that the exposure of soil bacteria to repeated dryrewetting cycles changed their growth response, resulting in an accelerated growth recovery to the pre-disturbance growth rate. Similarly, Papatheodorou et al. (2020) attributed the less pronounced impact on community composition induced by the second heavy rain incident, relative to that of the first one, to community acclimatization.
To explain the network collapse in the +BSC-H microcosms the authors considered reports on the strong inhibitory effect of C. rangiformis' secondary metabolites atranorin and usnic acid, specifically relating to instability over time, thermal and light liability, and the temporal changes in concentration (Farkas et al., 2020;Gheza et al., 2021) and reasoned as follows: Prior to the first rain event the crust survived inactive under extreme drought conditions for 2.5 months (including the summer period in the field and the conditioning phase). The first rainfall was adequate just to trigger the photosynthetic activity of the crust but not adequate enough to induce the leaching of the secondary metabolites, a fact that affected slightly the soil microbial community. Subsequently, the humid conditions remained favorable for some time due to the rain event so that during the second rain event a sufficient part of the accumulated secondary metabolites were flushed into the non-acclimatized soil microbial community and achieved destructive selection, finally resulting to the collapse of the network.

Modularity, stochasticity and resilience
The balance between stochasticity-determinism in regulating the interactions among microbiota was discussed in terms of network modularity. A non-modular network is mainly deterministically regulated while a modular one is regulated by stochastic processes (Zhou et al., 2014). According to those authors, dispersal limitation favors deterministic processes while the easy dispersal of organisms promotes stochastic processes that induce homogeneous architectures, enabling the easy spread of the consequences of environmental constraints throughout the network . However, despite the above-mentioned statement, the networks in most hydrated mesocosms supporting easy dispersal exhibited modular architecture. This could be explained as follows: Actually, initial watering increased the interface of soil microsites and thus attenuated the limitation in the dispersal of the microbiota (stochastic process). After each episode of heavy rainfall total microbial biomass was reduced , the networks' modularity was increased, and isolated nodes emerged (present paper). Most likely, this led to the weakening of deterministic regulation in favor of stochastic processes. The subsequent leaching because of heavy rain, however, induced differential filling of the micro and macro pores with water and differential distribution of ions leading to increased heterogeneity in the availability of oxygen and ions. This could result in hot spots of microbial biomass and activity where predation, and intraand inter-specific competition which belong to deterministic processes could take place. These processes seem to be dominant over the stochastic ones, ultimately increasing the contribution of the determinist forces to the regulation of the community structure, thus leading to non-modular configurations. Based on the above, changes in abiotic variables shifted the communities' regulation from stochastic to determinist or at the opposite direction (H3 hypothesis).
Modular networks tend to produce time-scale separation i.e., fast intra-modular processes and slow inter-modular activity resulting in a dynamic complexity (Meunier et al., 2010). In such networks, local repairs due to organisms' movement from neighboring clusters make the modular microcosms resistant and resilient against disturbance. This is consistent with the behavior of most modular networks analyzed in this study and partly supports the idea of community acclimatization, supporting our H1 hypothesis. By contrast non-modular networks, are expected to be more vulnerable (Scheffer et al., 2012), a fact that is supported by the final collapse of the non-modular architecture of the Small world network in the crusted only mesocosms and so supporting the hypothesis that the extent of acclimation is dependent on crust cover and previous watering of the crustsoil system (H1 hypothesis).

Network metrics and soil functionality
The potential relations between the structure of the microbial communities and certain functional variables are among the hot spot issues in ecology. Previously, studies explored the relation between biodiversity and various ecosystem functions, mainly biomass production (Tilman et al., 2014). Recently however, the idea of correlating network metrics with specific functions has gained ground (Stamou et al., 2019;Baquerizo et al., 2020;Baquerizo, 2022). The relations between network metrics and proxies of functions depicted on the PCCA graphs, confirmed the H4 hypothesis; most enzymes tended to be classified along with metrics pertaining to coherent networks with nodes having relatively high average degree and/or large neighborhoods and represent microbes with broad niches (Guseva et al., 2022), while the variables K, Ca, NH 4 + , Mg and SON were classified along with Fr and SW, pertaining to post-rain looser networks. In coherent networks with relatively high average degree and large neighborhoods, microbiota express a kind of functional complementarity to exploit the limited soil resources in the harsh Mediterranean summer. In looser networks, habitat filtering selects stress-tolerant microbes, while the intolerants died, thus enriching the soil with the respective chemical compounds. Specifically, for the relation between EC and SON the authors could hypothesize that the existence of some highly influential nodes (high EC) is crucial for organic N transformation, while for organic C transformation the synergy between nodes of almost similar influence is required.

Conclusions
Experimental manipulations reshaped the interaction pattern between the microbial members while in some cases even transformed the non-modular architecture into modular, thus favoring the rapid architectural recovery of temporarily fragmented networks caused by the heavy rain episodes (increased resilience). Nodes exhibiting high degree or clustering coefficient prevailing in coherent networks, tend to be related to enzyme activity while metrics of the loose networks seem to be related to nutrients such as K, Ca, Mg, and NH 4 + , confirming the idea that some architectural metrics can be used as useful indicators to depict the impact of stressors on ecosystem functioning. The response of microbial communities to disturbance is characterized by the tendency for the occurrence of alternating phases of networks' degradation-reconstruction. The exception in this study was the final collapse of the crusted only mesocosms. Given that crusts could occupy up to 70% of bare spaces in arid lands, it is obvious that the microbial communities beneath crusts are extremely vulnerable to disturbances imposed by recurring heavy rain events. Whether this pattern depends on the type of crust, or not, needs further indepth investigation.

Electronic supplementary material
Supplementary material is available in the online version of this article at https://doi.org/10.1007/s42832-022-0161-3 and is accessible for authorized users.

Open access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/ 4.0/.