Diel-scale temporal dynamics recorded for bacterial groups in Namib Desert soil

Microbes in hot desert soil partake in core ecosystem processes e.g., biogeochemical cycling of carbon. Nevertheless, there is still a fundamental lack of insights regarding short-term (i.e., over a 24-hour [diel] cycle) microbial responses to highly fluctuating microenvironmental parameters like temperature and humidity. To address this, we employed T-RFLP fingerprinting and 454 pyrosequencing of 16S rRNA-derived cDNA to characterize potentially active bacteria in Namib Desert soil over multiple diel cycles. Strikingly, we found that significant shifts in active bacterial groups could occur over a single 24-hour period. For instance, members of the predominant Actinobacteria phyla exhibited a significant reduction in relative activity from morning to night, whereas many Proteobacterial groups displayed an opposite trend. Contrary to our leading hypothesis, environmental parameters could only account for 10.5% of the recorded total variation. Potential biotic associations shown through co-occurrence networks indicated that non-random inter- and intra-phyla associations were ‘time-of-day-dependent’ which may constitute a key feature of this system. Notably, many cyanobacterial groups were positioned outside and/or between highly interconnected bacterial associations (modules); possibly acting as inter-module ‘hubs’ orchestrating interactions between important functional consortia. Overall, these results provide empirical evidence that bacterial communities in hot desert soils exhibit complex and diel-dependent inter-community associations.

Macro-and micro-environmental parameters are critical in governing community assembly and function 1,2 . This is especially true in hot desert ecosystems, where extreme environmental conditions have led to discrete adaptations (i.e., collectively leading to niche specialization) for many plants and animal species [3][4][5] . Such extreme conditions are also considered to be key drivers of both the structure and function of edaphic microbial communities [6][7][8] . Given that climate change models predict major shifts in future temperature and precipitation regimes for dry-land regions (e.g., 10-30% reduced precipitation in subtropical arid zones by the end of the 21st century 9 ), changes in complex multi-trophic responses over short and long time scales are likely to impact significantly on core ecosystem processes (e.g., carbon storage in soils 10 ). In consequence, ecologists have become increasingly interested in further understanding the influence of temporality in shaping community structures and functions in hot desert soils [11][12][13] , which in turn can aid in predicting future ecological impacts.
To date, research on microbial temporal dynamics in hot deserts has been largely focused on diversity patterns over seasonal and yearly scales [14][15][16][17] and/or after a wetting event [18][19][20] . Although key insights into the causative factors (deterministic vs. stochastic) driving microbial community structure have been provided across biogeographic space 21,22 , comparatively little is known about the fine-scale temporal changes in edaphic microbial community structure and function. Indeed, other than studies focused on Cyanobacteria 23,24 and a temperature manipulation experiment 25 , short temporal community dynamics have been neglected. This is surprising given the fact that significant differences in specific desert soil processes have been recorded over the diel (24 hrs) cycle (e.g., CO 2 flux refs [26][27][28]. In addition, asymmetrical diel warming; i.e., the increase in daily minimum temperature (T min ), is widely observed and set to continue in terrestrial ecosystems 29,30 and is likely to impact negatively on microbial communities and their processes 31 .

Results and Discussion
The light-dark (diel) cycle drives behavioral responses that are an essential component of life in hot desert systems. In higher organisms, well-characterized circadian rhythms facilitate adaptive strategies within these habitats that include the reproductive behaviour of desert locusts 54 , heterothermy in desert mammals 55 and in desert plant physiology 56 . It has also been reported that the globally distributed circadian gene locus, cpmA, is highly conserved in key microbial phylotypes such as Cyanobacteria 57 , which suggests a key role in the maintenance of intracellular homeostasis and adaptation 58 . For mixed communities, light-dependent changes in microbial gene expression are also a feature in hot desert biological soil crusts 23 . By focusing on RNA isolated from bacterial communities in situ, we provide evidence that changes in active bacterial groups can occur over fine time scales and that both environmental selection and biotic associations are crucial factors within this system.

Site characteristics.
No rainfall or fog events were recorded on site during the experimental period.
Historically, fog events have provided a key source of bioavailable moisture in Gobabeb 37 . Fog events have an interdependent relationship with atmospheric humidity, which tends to reach maximum levels in the morning before dissipating quickly after sunrise 37,52 . During this study, in situ soil relative humidity (RH) levels were highly variable over each diel cycle ( Fig. 1; time-points; P < 0.001). Soil RH values ranged from 25% (midday, day 4) to 59% (night, day 1). The diel RH profiles differed between days with the RH profile on day 1 was found to vary significantly from those of days 2, 3 and 4 (post hoc testing with Tukey-HSD comparisons; P < 0.001). An increase in peak soil temperature after day 1 (day 1, 47 °C, days 2 to 5, ≥ 50 °C) may have accounted for this shift in RH values. Physiocochemical parameters (pH, organic C, conductivity, mass concentrations of P, Ca, K, Mg, Na, S) were relatively homogenous over the entire sampling period (Table S1). Significant differences were found between days for nitrate (F 4.45 = 3.54, P = 0.017) and magnesium (F 4.45 = 3.01, P = 0.037) but this was likely to be a consequence of spatial heterogeneity, with a nutrient 'hotspot' recorded for a single biological replicate sampled on day 2 (Fig. S1).

Changes in active bacterial communities revealed by T-RFLP analysis. Community fingerprinting
by T-RFLP rRNA gene (cDNA) analysis recorded a total of 27 unique terminal restriction fragments (TRFs), likely representing the most active bacterial taxa present in this soil system 59 . The number of TRFs per sample (α -diversity) ranged from 12 to 23 with only 4 (15%) and 2 (7%) TRFs not detected for each day or time-point, respectively. Despite the absence of discernable clustering patterns for both Morning and Night bacterial communities through nMDS ordination, the Midday communities did cluster relative to the other time-points (Fig. 2a). Significant shifts in potentially active bacterial communities was recorded between time-points was confirmed by PERMANOVA (F 4,45 = 5.51, P = 0.007). Interestingly, the T-RFLP community profiles were not found  to be significantly different between days (F 4,45 = 0.15, P = 0.961). This was also found in 454 pyrosequencing analysis using the same samples ( Fig. 2b; time-points, F 4,44 = 1.12, P = 0.04) Considerable diversity in hot desert soil bacterial communities. In order to provide further insight into these findings we performed deep-sequencing analysis. 454 pyrosequencing of 16S rRNA gene amplicons (generated from the same cDNA pool as for T-RFLP analyses) resulted in a total of 4348 OTUs (96 316 sequence reads, n = 45) at the 97% cutoff (cutoff level for all OTUs if not stated otherwise), with 1414 singletons (32%) and 606 doubletons (14%). After sub-sampling this dataset for sequence consistency (1702 OTUs > 15,000 reads), the number of singletons and doubletons increased to 47% and 15%, respectively. The large number of 'rare' OTUs identified here is not uncommon for hot desert soil communities 6,38 and can represent novel diversities within the rare biosphere 60 . Alpha-diversity was relatively consistent between samples regardless of the diversity metric used; i.e., with no significant differences between days or time-points (Table 1). For many of the communities sampled, rarefaction analysis exhibited approximate plateau plot behavior ( Figure S2), suggesting that sampling depth was sufficient to properly assess the bacterial diversity within this system. A total of 20 bacterial phyla were identified (Fig. 3, Figs S3 and S4 [with standard deviation]), with members of 8 phyla representing ~95% of the total bacterial community (in order of abundance; Actinobacteria, Proteobacteria, Cyanobacteria, Gemmatimonadetes, Chloroflexi, Bacteroidetes, Firmicutes, Armatimonadetes). The proportion of classifiable sequences decreased substantially from phylum (98%) to genus (10%) levels, showing Namib Desert soils present numerous novel representatives of well-known edaphic phyla 41,61 . Nevertheless, the vast majority of sequences could be classified to class level (95% of sequences). The ubiquitous Actinobacteria (44.8% ± 5.6% 16S rRNA gene sequences; mean between time-points ± s.d.) and Proteobacteria (38.2% ± 3.6%) dominated all soil samples and also exhibited very high diversity (> 300 OTUs classified to family level). This is consistent with previous phylogenetic surveys of bacterial communities in hot desert soils and soil crusts, including sites in the Negev 33 , Atacama 62 and Namib Deserts 63 . Cyanobacteria comprised fewer OTUs than many of other bacterial phyla but were present at relatively high abundance (6.44% ± 1.37%; Fig. S5). The reverse was found for the Gemmatimonadetes, Armatimonadetes, Chloroflexi and Bacteroidetes phyla.
Strikingly, it was found that only a small fraction of the OTUs were detected at each time-point (i.e., core community -19.4%; n = 328; Fig. 4), with the amount of OTUs recorded only at a single time-point being relatively high. For example, the night-specific OTUs represented > 25% of the total bacterial activity. This pattern was also found when each day was analyzed separately (Fig. S6). PERMANOVA confirmed that the bacterial community composition exhibited significant structural shifts over each diel cycle (F 2,44 = 1.3370, P = 0.044), supporting observations from T-RFLP analyses. We also found that the majority of the temporally-specific OTUs were members of the rare bacterial fraction (> 60% of total OTUs). Indeed, the core community (< 20% of total OTUs) constituted the majority of total sequence reads with a mean relative abundance value of 78% for all samples. In order to ascertain whether this rare fraction were primarily responsible for the recorded shifts in community structure, we applied the same approach with the rare fraction removed (i.e., with no singletons or doubletons), thus reducing the time-point specific OTUs from 1035 [61.4%] to 21 [4%]). A similar statistical result was recorded (Table S2), confirming that prevalent members within these hot desert soil bacterial communities also exhibited a significant shift in relative abundance patterns over a diel cycle.
Phylogeny and distribution of diversity. Relative abundance patterns of the bacterial communities were analysed by aggregating all taxa at the phylum, class and order levels. Overall, the phyla Actinobacteria and   Proteobacteria exhibited differential responses over diel cycles ( Fig. 3; Fig. S3), with a strong negative correlation between their respective relative abundances (Pearson; − 0.069, P < 0.01). The former exhibited a significant reduction in relative abundance from Morning to Night (ANOVA; (F 2,44 = 5.229, P = 0.009) with a consistent pattern between days (P = 0.13). The diel patterns of sub-phylum phylogenetic groups generally reflected that at phylum level (Fig. S7). There were certain exceptions, however, with the Micrococcales order of the Actinobacteria class recording an increase from Morning to Night ( Figure S7b). Within the phylum Proteobacteria, both β -and δ -Proteobacteria occurred consistently at relatively low abundance (< 4% each) while γ -Proteobacteria showed a consistent increase in relative abundance from Morning to Night samples ( Figure S7a), which was primarily attributed to the order Pseudomonadales ( Figure S7b). Furthermore, the increased abundance of Cyanobacteria in Night samples could be primarily attributed to the order Nostocales ( Figure S7b), which displayed a > 6-fold increase from Morning to Night samples on day 3. The Nostocales group are crucial in maintaining soil N budgets through N 2 -fixation 64 , suggesting a capacity for active nitrogen cycling within this desert soil system that may be governed by temporal cues. No other phyla showed major or consistent abundance changes between time-points, although Armatimonadetes, Bacteroidetes and Planctomycetes were found to be statistically different between days (ANOVA, P < 0.05).
Influence of abiotic factors on desert soil bacterial community structure. In order to probe the drivers of fine-scale bacterial population dynamics, we tested the possible effects of environmental factors on the bacterial community composition for all samples. Multiple rank correlations (Spearman; global 'BEST' test) indicated that the abiotic parameters measured could not account satisfactorily for bacterial composition over time (P = 0.62, 99 permutations). Re-analysis of individual parameters using linear-constrained ordination redundancy analysis (RDA) with forward selection revealed that, for the physical and chemical parameters, climatic factors (temperature, F = 1.981, P = 0.02; RH, F = 1.338, P = 0.043) together with soil phosphorus content (F = 1.373, P = 0.039) had a significant influence on community composition (Fig. 5). However, these parameters only accounted for 10.5% of the total variation, suggesting that other factors were important in driving the community dynamics recorded in this study.
The extent of individual environmental variables on taxon distribution was also tested by Spearman correlations. Actinobacteria, Acidimicrobiales, Solibrobacterales and both β -and γ -Proteobacteria had strong correlations with soil temperature and %RH, while soil phosphorus content was positively correlated with the Rhizobiales, Rhodospirillales (α -Proteobacteria) and Chroococcales (Cyanobacteria) orders (P < 0.05). Interestingly, two orders of Cyanobacteria, the Nostocales and Oscillatoriales, were also positively correlated with nitrogen elements (NH 4 + and NO 3 − ) (P < 0.05). The finding that Nostocales group were positively correlated with nitrogen elements further enhances our view that this group may be key in driving N 2 -fixation within this desert soil system.
Time-dependent bacterial associations revealed through co-occurrence networks. Given that the physicochemical parameters measured could not account satisfactorily for the recorded community dynamics, we investigated the possible influence of biotic associations during diel cycling. In order to do so, we generated co-occurrence networks for each time-point incorporating the five days of data ( Fig. 6; organic layout algorithm). This mathematical approach has been successfully employed to infer phylogenetic and functional trait associations for edaphic and hypolithic microbial communities 33,51,53 . In the Namib Desert, this approach has previously been used to show the differences between hypolithic communities from fog-and rain-dominated zones and highlight the relevance of Cyanobacteria in potentially driving hypolithic food web compositions 51 . The networks presented in Fig. 6 comprised the most dominant fraction of the core community; i.e., the 100 most abundant and consistently detected OTUs (84% of total reads). Although present in all soil samples, this core group exhibited a reproducible non-random co-occurrence pattern against the null model at each time-point (Average C-score 63.51, P ≤ 0.001). Interestingly, the average path length (distance between all pairs of nodes), clustering and topological (tendency of nodes to share neighbors) coefficients were different for each network ( Table 2), indicating that the system was interactive and dynamic (average of > 5 edges [associations] per node [OTU] for each network) 53,65 . A level of structural homology was also evident as a primary group of nodes (module) formed a highly clustered topology at each time point; i.e., with a 'neighborhood connectivity' (NC, which corresponds to the number of connections between neighboring OTUs) over 7. Globally, the Midday network recorded a considerably higher NC value (6.38) when compared to the Morning (5.1) and Night networks (5.24) ( Table 2). When attempting to understand the recorded patterns recorded here it is important to note that co-correlation is not necessarily indicative of biotic interactions. Indeed, the recorded patterns relating to midday community dynamics may be due to these bacterial groups responding to the variability in environmental cues (e.g., temperature) in a similar way, which would appear as co-correlation. However, modularity in community ecology can also reflect functional associations with physical contact and/or phylogenetic clustering of closely related species 66 , which has been reported to have a possible role in negating the effects of high temperatures disturbance 67 . Given the extreme temperature conditions associated with the Midday soil habitat in this system (up to 52 °C on day 2; Fig. 1), there is a likelihood that metabolic and/or functional associations could occur between adjacent temperature-tolerant bacterial taxa such as the capacity to metabolise substrates not accessible to non-thermophilic microbial communities 68 . Interestingly, Actinobacteria inter-phyla co-occurrence patterns (represented by OTUs with at least 2 connections) peaked at Midday (57%), the majority of which were positive (Table S3). Many of the positive associations within the primary module were between closely related Each node (circle) in the network symbolizes a unique OTU at the 90% threshold (family level) and the size relates to the mean relative abundance within that network directly derived from the related read counts. Each edge (connection) relates to a strong (Spearman's ρ > 0.6) and significant (P-value < 0.01) correlation. The OTUs are numbered in accordance with overall relative abundance. radiotolerant phylotypes such as Patulibacter (OTU2) and Solirubrobacter (OTU5) 42 (Fig. 6). This high degree of modularity incorporating similar phylotypes may partly explain the community divergence between Midday and the other time-points. Interestingly, Proteobacteria-affiliated OTUs exhibited the lowest level of intra-phyla co-occurrence (both positive and negative) at Midday (21%), being 2.3-and 2.2-fold lower than the Morning and Night networks, respectively, which suggest that proteobacterial taxa are more involved in multi-trophic associations during daytime habitat conditions. The majority of Cyanobacteria OTUs were located externally to the primary module for each network (Fig. 6). Considering the physiological distinction between the photoautotrophic Cyanobacteria and the heterotrophic Actinobacteria and Proteobacteria, differences in co-occurrence profiles at different times of the diel cycle may be an indication of very fine time-scale niche partitioning 69 . Alternatively, these Cyanobacteria driven tertiary connections may function as inter-module hubs 51 , such as in the Midday network, where Cyanobacteria nodes linked the core module with a submodule that comprised a highly abundant Proteobacteria OTU (Fig. 6). Furthermore, a highly abundant (4% of morning reads), highly associated (7 connections; 6 positive/1 negative) Cyanobacteria OTU assigned to the order Oscillatoriales was positioned as a module hub in the Morning network, with strong inter-phyla positive co-occurrence, which is indicative of mutualistic associations. Interestingly, this module was comprised of taxonomically unique groups (distinct at class level), which may provide evidence of naturalization 70 , a concept, which suggests that co-occurring species are likely to be metabolically dissimilar so as to reduce the risk of competitive exclusion. The influence of Cyanobacteria OTUs as modular hubs was not recorded for the Midday and Night networks. Even low abundance phyla exhibited a high level of co-occurrence. For example, Acidobacteria (class Holophagae) only had a single representative in the core community but had 25 significant correlations (mostly positive), the majority of these being with much higher abundance OTUs in the Morning dataset ( Fig. 6; Table S3). Acidobacteria are ubiquitous as K-strategists in edaphic niches characterized by low plant-derived carbon sources and generally more oligotrophic niches 71,72 . Although the functional relevance of this phylum remains unclear, the high rate of inter-phyla co-occurrence would suggest a pivotal role in edaphic processes within this hot desert soil system, with low relative abundance recorded here possibly due to an underestimation of taxa with low 16S rRNA copy numbers per genome such as Acidobacteria 73 .
Overall, our results clearly demonstrate that significant changes in the active fraction of Namib Desert soil bacterial communities can occur within a single day, with this study being the first to provide direct evidence of short temporal changes in undisturbed hot desert soil communities. Furthermore, we have shown that these changes are not primarily governed by abiotic factors (i.e., physicochemical variables and climatic regime) as proposed by our leading hypothesis, and that distinct microbial interactions may be important in the recorded community dynamics. The bacterial co-occurrence matrices presented here are intriguing and provide an initial glimpse of potential trait associations and mechanisms driving community organization within this system. For instance, phylogenetic clustering of desert soil bacteria can be influenced by the filtering of traits over fine-time scales, which may confer resistance to abiotic stress or competitive exclusion (e.g., high temperature resistant phylotypes; Midday function). With previous investigations focusing on refuge niche habitats as the 'core hubs' of microbial function and diversity in desert systems 63,[74][75][76][77] , these findings are all the more striking and have implications for future desert soil research. Incorporation of additional meta-omics approaches such as metaproteomics as well as measuring photosynthetically available radiation (PAR) would be beneficial in resolving the findings recorded in this study.

Experimental Procedures. Site description and sample collection. Sampling was performed between 22
and 26 April 2013 at a single undisturbed gravel-plain site (10 × 10 m quadrat; Fig. 1) adjacent to the Gobabeb Research and Training Centre, Namibia (23°33′ S, 15°02′ E). This region is classified as hyper-arid 8,37 with low mean annual precipitation (< 30 mm) and receives much of its moisture via fog events 37 .
In order to limit sampling bias and to ensure sampling homogeneity, the site was divided into a series of quadrats (x16; 2.5 × 2.5 m), which were further subdivided into four mini-quadrats (1 m 2 ). Eight surface (0-5 cm) pseudo-replicate soil samples were recovered and pooled from three randomly selected (RAND()*N function in Microsoft excel) mini-quadrats at time-points reflecting the diel cycle -morning (6 am), midday (1 pm) and night (8 pm) over a 5-day period (n = 45). Each soil sample (~10 g) was supplemented with 5 ml of RNAlater (Sigma-Aldrich, Copenhagen, Denmark) and stored at − 80 °C for later molecular analysis. For physicochemical analyses, ~200 g of soil was taken from each sampling zone and stored in sterile Whirlpak ® bags at 4 °C.
Climatic data and soil physicochemical analysis. Temperature (°C) and relative humidity (%RH) readings were recorded at 30 min intervals for the 5-day period using Hygrochron iButton temperature and humidity loggers (model DS1923; Embedded Data Systems, Lawrenceburg, KY, USA). These were placed at a depth of 2 cm at the center of each 2.5 m 2 quadrat (16 in total). Physicochemical analyses were conducted at the Soil Science Laboratory of the University of Pretoria (South Africa) using standard quality control procedures 78 . Soil samples were sieved (2 mm) prior to analysis. Total organic carbon (%) was determined by a dichromate oxidation procedure 79 . A correction factor of 1.32 was used to account for incomplete oxidation of organic C 80 . Ammonium (NH 4 + ) and nitrate (NO 3 − ) ions were determined using steam distillation 81 . The P-bray method 82 , with Inductively Coupled Plasma-Optical Emission Spectroscopy (ICP-OES) was used to determine total organic phosphorus. The cation exchange capacity (CEC) was determined by 0.2 M ammonium acetate treatment, with subsequent concentrations of Fe, Ca, K, Mg and Na determined by ICP-OES. Soil pH was measured using a slurry technique (1:2.5 w/w soil/deionised water ratio) with a calibrated pH meter (Crison basic 20, Barcelona, Spain). The soil texture was determined by a fluid suspension method using a specialized hydrometer, based on Stokes' law of particle dispersal 83 .
Bacterial T-RFLP analysis of the cDNA products was done using 16S rRNA gene PCR primers 341 F (5′-FAM-CCTACGGGAGGCAGCAG-3′ ) and 908 R (5′ -CCGTCAATTCCTTTRAGTTT-3′ ) 84  with the LIZ ® 600 Size Standard. The resulting fragment profiles were exported and analyzed using Peak Scanner 1.0 (Applied Biosystems) with the default 5RFU cutoff. In order to eliminate any inaccuracies in the data, true peaks and fragments of similar size were identified and binned using the software R and Perl 85 , Fragment lengths within 1 bp of each other were considered to be the same operational taxonomic unit (OTU) and peaks within 3 standard deviation of the noise baseline were removed. A matrix of the T-RFLP community fingerprint patterns ('ClusBinMatrix') for the relative abundance (area under the peak) of individual OTUs for each soil sample was generated for subsequent analysis.  87 . A stringent approach to this analysis was used, with no primer and barcode sequence mismatches permitted. Errors in sequencing were reduced by the trim.flows and shhh.flows commands, the latter incorporating an implementation of PyroNoise 88 that uses the expectation-maximization algorithm. Sequences with ambiguities and/or shorter than 200 bps were removed. After configuring the data to comprise unique sequences, further de-noising was done using the 'pre.cluster' command, which employs a pseudo-single linkage algorithm to remove sequences likely due to pyrosequencing errors 89 .
Finally, the sequences were checked for PCR chimeras using the chimera.slayer command 90 . Aligned sequences were used to construct distance matrices and to cluster the sequences into operational taxonomic units (OTUs) at the 97% (species level) and 90% (family level) cutoffs in MOTHUR. Bacterial sequences were classified using CREST 91 with the Silvamod database and MEGAN with LCA parameters min bit score 250, top percent 2 (50 best blast hits). Rank abundance data were generated for all samples with rarefaction curves and other population diversity indices (for example, Chao1, Shannon evenness) calculated in MOTHUR after we rarefied the species richness to the same number of individuals to ensure sequence consistency between samples 92 . Venn diagrams for days and time-points were generated using the R statistical package (v.2.14.0) 93 .
The SFF file containing the sequence data is available at the NCBI Sequence Read Archive under the accession number PRJNA253979. Data analysis. All statistical analyses were carried out using R, including functions from the vegan package (version 2.0-9). For T-RFLP and pyrosequencing analysis, the dependent variables, day (i.e., day 1, day 2, day 3, day 4 and day 5) and time-point (i.e., morning, midday and night) were treated as fixed factors. QQ-plots and frequency histograms indicated that residuals did not meet the assumptions required for parametric tests 94 . Therefore, independent variables (x) (i.e., all soil parameters and all community data) were transformed according to ln (x + 1).
Linear models were used to determine significance between both days and time-points for the different soil parameters (lm function). T-RFLP and pyrosequencing community profiles (employing Hellinger transformation) were visualized using non-metric multidimensional scaling (nMDS) analyses (meta.mds function) using the Bray-Curtis dissimilarity metrics. In order to test for significant differences in bacterial composition between samples, permutational multivariate analysis of variance (PERMANOVA) was used with either day (day 1, day 2, day 3, day 4, day 5) or time-point (morning/midday/night) used as factors (Adonis function). We also tested for an interaction effect between days and time-points (model; day*time). In order to determine which experimental groups were distinguishable, post hoc Tukey tests were used for pairwise comparisons. In order to test the influence of abiotic factors on variations within bacterial community structure, Biota Environment STepwise matching (global 'BEST' test) in the Primer6 statistical package (Primer-E Ltd, UK) 95 was undertaken. The soil physicochemical data was standardized. Spearman rank correlations between the abiotic variables and the bacterial Scientific RepoRts | 7:40189 | DOI: 10.1038/srep40189 communities were then calculated over the sampling period (99 permutations). Redundancy analysis (RDA; rda function) with forward selection was used to evaluate the effects of the environment on community composition.
Co-occurrence networks were generated comprising consistently detected and highly abundant OTUs; i.e., the 'core' community. Non-random co-occurrence patterns were first tested with the checkerboard score (C-score) under the null hypothesis of random community assembly 96 , where 5,000 matrices were randomly generated from the 16SrRNA gene data (oecosimu function with nestedchecker in bipartite package). All possible Spearman rank correlations between OTUs were then calculated with stringent parameters set (interaction between OTUs valid if > 0.6 value for Spearman's correlation coefficient [ρ ] and P-value < 0.01) as previously described 53 . A graphical representation of the filtered dataset was prepared using Cytoscape 97 . The nodes in the constructed network represented OTUs with the edges (connections) indicating significant inter-nodal correlations. Topological coefficients of the resulting networks were then calculated in order to derive comparisons, which included average node connectivity, average path lengths, cumulative degree distribution and clustering coefficients.