Microbial community response to hydration-desiccation cycles in desert soil

Life in desert soil is marked by episodic pulses of water and nutrients followed by long periods of drought. While the desert flora and fauna flourish after rainfall the response of soil microorganisms remains unclear and understudied. We provide the first systematic study of the role of soil aqueous habitat dynamics in shaping microbial community composition and diversity. Detailed monitoring of natural microbial communities after a rainfall event revealed a remarkable decrease in diversity and a significant transition in community composition that were gradually restored to pre-rainfall values during soil desiccation. Modelling results suggest a critical role for the fragmented aqueous habitat in maintaining microbial diversity under dry soil conditions and diversity loss with wetting events that increase connectivity among habitats. This interdisciplinary study provides new insights into wetting and drying processes that promote and restore the unparalleled microbial diversity found in soil.

: Chemical analysis   Average qPCR counts of rRNA units Figure S2. Total nucleic acids were extracted from the soil samples as previously described (Angel 2012). The extract was purified by MasterPure RNA Purification Kit (Epicentre, Madison, WI). The DNA was degraded by DNase I supplied with the kit and the RNA samples were stored at -80 • C for further analysis. Total bacterial rRNA unit quantification in the field measurements. Colour coding shows the different categories of samples with respect to community clustering (see Figure S1B in the main text).  Figure S3. Result indicating statistically significant difference between the four observed NMDS clusters in Figure S3. Consequently, the resulting significant values of R and p support rejection of the null hypothesis that these were a result of a random process.     The rough surface patch model (RSPM) is employed and modified to generate physical do-56 mains (Kim and Or 2016). In the RSPM, geometrical information of soil structure was averaged 57 with a probability distribution of angular pore sizes and the effective water film thickness was 58 introduced as an indicator of hydration condition that controls substrate diffusion, microbial dis-59 persion rates, and aqueous habitat connectivity. In this study, we extended the probability based (pH(r) + sr) r 2 N (r)dr (SI.1) where r min and r max are the cutoff values for the pore size distribution, H(r) is the height of the square pyramid pore with size r (Here, we assumed that the height of each square pyramid pore is the same as its base.), N (r) ∼ r −D is the probability density function of pore size r, p and s are the fraction of pore or solid elements in the space, respectively. We note that this relation could be different when different shape factors for solid and pore are assumed. The saturation degree at the given matric potential ψ m can be obtained following RSPM.
where V(r, ψ m ) is the amount of water held in the pore with size r at the matric potential ψ m due 78 to the capillary force and the van der Waals force and h µ is the absorbed water on the smooth 79 surface (For the detailed explanation, see Kim and Or (2016)). The roughness representation of a surface is extended to the vertical representation of three-dimensional soil profile including the gas phase (Kim and Or 2016). A patch is a collection of roughness elements, cubic solid blocks and square pyramid voids. Using the effective medium assumption, effective thicknesses of each phase are determined. d s , d w , and d g denote the effective thickness of solid, liquid, and gas phases, respectively. The porosity φ, the gas content , and the water content θ of each patch are converted from these effective thickness values. (B) Effective thicknesses for liquid and gas are calculated for three different fractal dimensions for pore size distribution. Strong coloured lines are water film thicknesses and corresponding pastel colours are for the void thicknesses. (C) For each patch, d s + d w + d g = d is given constant (calculated from the preassigned porosity) and the hydration condition and the thermal diffusivity are regulated by these effective thicknesses, d w , d e , and d g . As a result, nutrients and heat diffusion between adjacent patches is a function of ψ m , local relative humidity.
where d is the total thickness of each patch for heat and solute transport and it can be calculated following equation: In this research, the hydration condition (wetting and drying event) is controlled with the matric potential ψ m mapped from the water contents (θ(ψ m ) ≡ Θ(ψ m )φ). The heterogeneity of the domain was achieved by assigning a set of {Φ, D} for an individual patch. To match the type of soil investigated in the field (Muñoz Castelblanco et al. 2012), we assumed the fractal dimension of the pore-size distribution of loam soil 2.8 (corresponds to a fractional Brownian surface with the Hurst exponent H = 0.2) and the mean porosity to beφ ≈ 0.3 (Tyler and Wheatcraft 1992; Wang et al. 2005;Huang et al. 2006). This means that the probability that a point on the domain belongs to a pore with size X in the interval [r, r + δr] is; Parameters used to generate the physical domains for the simulations are given in Table S1.

Diffusion processes through the profile
In the physical domain, we averaged microscopic details of the pore distribution and assumed that hydration conditions are represented with effective water film thickness and degree of saturation. The diffusion processes in this work is described as below to calculate local substrate concentration, C( r, t), ∂C( r, t) ∂t = ∇ · (D( r)∇C( r, t)) − Sink terms + Source terms (SI.9) where D( r) is the apparent diffusion coefficient defined from the effective film thickness distribution of adjacent patches and the effective diffusion coefficient including tortuosity as a function of porosity and water contents. The net flux between two adjacent patches (for example, patch 1 and patch 2) due to diffusion is calculated as following: To calculate the flux between heterogeneous medium, we have chosen the harmonic mean of diffusion coefficient and the minimum value of water film thicknesses between neighbouring patches (for the details, see Kim and Or (2016)). The effective diffusion coefficient of a patch is given following the Milington-Quirk tortuosity model (Milington and Quirk 1961).
where D 0 is the diffusion coefficient of the substrate in bulk water and θ( r) and φ( r) are the Furthermore, we included dynamics of oxygen profile driven by diffusion in the gas phase and its input into the liquid phase as the dissolved oxygen. Aerobic bacteria are assumed to uptake the oxygen in the dissolved form (modelled as obligate aerobes) and anaerobic bacteria are inhibited by the local concentration of dissolved oxygen (modelled as obligate anaerobes). However, in the model, combining IBM and the gas diffusivity of oxygen through the soil profile is challenging as the time scales of these processes differ in order of magnitudes. For example, while the growth of a cell is in the order of hours, ≈ 10 3 seconds, the gas diffusion through the domain (in this work, depth of 5 cm) is in the order of ≈ 10 −2 seconds. Using the time scale of gas diffusion for all processes in the model is not plausible due to the computational time limit. Thus, we did not explicitly solve the gas diffusion for the oxygen source. Instead, we combined the Henry's law and percolation theory: Firstly, for the mass transfer between gas and liquid, Henry's law was applied; C O ( r) * = H cc (T ( r))C g O ( r) * (SI.12) where C O ( r) * is the local concentration of dissolved oxygen in the liquid phase, C g O ( r) * is the 92 concentration of oxygen in the gas phase (converted from the partial pressure of the oxygen 93 in the atmosphere) at equilibrium, and H cc (T ( r)) is the dimensionless Henry's constant when 94 temperature is given as T at the position r. The temperature dependency of Henry's constant R is the gas constant, T is absolute temperature, and Θ refers to standard condition (T Θ = 97 298.15K) (Sander 1999). This assumption holds during entire simulations, as the soil matrix is The heat transport equation can be used to calculate the soil temperature profile T ( r, t) (in the absence of fluid motion): where c v ( r) is the local heat capacity and λ( r) is the local thermal conductivity at r. From the local information of volume fractions and densities of solid, water, and gas, the soil volumetric heat capacity can be written: where c s is the specific heat capacity per unit mass, ρ is the density, and s, w, and g for each variables denote soil minerals, water, and gas. The volume fraction of each phase is given as (1 − φ), θ, for solid, liquid, and gas, respectively. In this equation, the proportion of organic matter (such as EPS, or microbial cells) are ignored for the thermal properties. Each volume fraction is determined by the hydration condition, therefore c v ( r) varies following the water contents ( Figure S4). The effective thermal conductivity at r, λ( r) is given as a harmonic mean of three conductivities, λ s , λ w , and λ g from different phases.
(SI.15) However, the thermal diffusivity in soil is also in the order of ≈ 10 −6 m 2 /s, hence the coupling the microbial growth is not plausible. For the temperature, we solved the homogenous 1 dimensional domain (over the depth) at the varying hydration condition and applied the solution to the profile as boundary conditions. Considering that the soil crust (with 2mm thickness) was discarded and the sample was collected from top 5cm, we applied soil crust (fine soil structure with D = 2.9 and Φ = 0.9) at the top 2mm and used Loess soil texture (D = 2.8 and Φ = 0.8 same as the value we used in the main text) for the below crust. Boundary condition at the air-soil interface is assumed to be in the interfacial isothermal condition.
where T air (t) is assigned from the air temperature records of LTER where our field measurements 106 were conducted. Considering the thermal damping depth of loess soil is around 10 cm, we 107 modelled up to 15 cm with the zero heat flux boundary condition at the bottom of the domain.

108
The calculated solution at the three selected depths are plotted in Figure  The predicted soil temperature is given for depths, z = 0.2 (red), 2.5 (yellow), and 5 cm(purple). Only top 5cm of the soil is sampled after discarding the soil crust with 2mm thickness.   (B) Hydration and temperature conditions are used as dynamic boundary conditions through the entire simulations. For the hydration condition, measured gravimetric water contents (blue bars) were linearly interpolated (the red line) and mapped to the matric potential of the entire domain. By using the air temperature records of LTER, we calculated the soil temperature (combined with the hydration condition changes) over the profile. Orange bars are the record of air temperature, red, green, blue lines are the soil temperatures at depth 2mm, 2cm, and 5cm, respectively. (C) Dissolved nutrient distributions are solved with diffusion-reaction equations. Especially, the oxygen source in the profile (for aerobic cells) is the partition of gas and liquid phases and the model solves the input of dissolved oxygen from gas phase by combining gas percolation and Henry's law.

Individual based description of microbial growth on the heterogeneous domain
Microbial activity was added to the modified RSPM by using the IBM (Individual based model) (Kim 113 and Or 2016). The RSPM allows assigning local environmental conditions and the IBM examines  (2009)). As an addition to previous IBM models, our model 121 is spatially explicit and includes motility of cells as they are able to actively explore the domain.

122
The motility of a cell is regulated by the chemotactic response to substrate concentration and 123 the swimming velocity on capillary surfaces. For its locomotion and net displacement, we used a 124 biased random walk approach on the rough surface.

125
A reaction-diffusion equation was used in the model to obtain the nutrient distribution over time: where C j ( r, t) is the local concentration of substrate j, D j ( r, t) is the local diffusion coefficient, Each cell consumes several chemical species that are obligatory for its growth, in this study oxygen and dissolved carbon following their physiological differences. We assign a growth rate of a taxon (or microbial species) i with multiple limiting substrates j by using Monod-type growth kinetics as a function of the substrate concentration field (Monod 1942(Monod , 1949: (when nutrient j is a inhibitor for the growth) and µ max i , K i S /K j I,i are the maximum growth rate and half-saturation/inhibition constant of cell i, respectively. Specifically, in this study, we introduce two different groups of microorganisms, aerobically growing group and anaerobically growing group. Aerobes are simply assigned as a group that utilises oxygen and carbon source for growth.
On the other hand, growth of anaerobes are inhibited by the presence of oxygen, meaning that only obligate anaerobes are considered. The growth rate of individual cell is calculated as following: For aerobic growth, For anaerobic growth, where C C ( r) is the local concentration of the carbon source and C O ( r) is the local concentration 132 of the dissolved oxygen. We assumed that the individual growth rate at the pore scale is the 133 same as the population growth rate in batch culture (Dai et al. 2013). Accordingly, the growth 134 of an individual cell can be written as: where m i is the maintenance rate of cell i.  In the model, we combined our growth model with a temperature dependent growth model using the Arrhenius equation (Schoolfield et al. 1981) to see the effect of temperature in population dynamics. The level of adaptation to the temperature might vary among different organisms, however, we did not consider it and assumed that all microbial cells follow the same activation energy and the optimal temperature encapsulated in the maximum growth rate to reduce the complexity. The temperature dependency on the maximum growth rate of a cell at temperature T is assigned as following (Schoolfield et al. 1981): ae Values assigned for the aerobically growing cells; an Values assigned for the anaerobically growing cells: Assuming that anaerobic processes are costly; oligotroph-like aerobes and copiotroph-like anaerobes based on observed values (Pirt 1965;Heijnen 1999;Stouthamer 2012). We note that higher maximum growth rates for anaerobes were chosen to compensate the strong inhibition of oxygen # µmax and K C S,i are different for each taxon. µmax are chosen uniformly spaced values and K C S,i are logarithmically spaced values in the given ranges a The growth of obligate aerobes ceases at 25% of atmospheric level b The growth of obligate anaerobes is inhibited when oxygen concentration is higher than 0.5% of atmospheric level Active cells and potentially active cells before and after wetting event

148
The model aimed at describing microbial community behaviour under a wide range of matric 149 potential, changes from several kilopascals to megapascals; in other words, the domain can be 150 exposed to very wet conditions after a rainfall event and to very dry conditions after prolonged 151 desiccation. As the soil dries after the rainfall event, the effective water-film thickness can be 152 reduced from 10 −5 m to 10 −9 m. This implies that the amount of nutrient flux to a certain 153 location will also be reduced by several orders of magnitude due to thinning of the water film.

154
In the model, this suppresses microbial growth and leads to taxa extinct in the long term. To