Deterministic Selection Dominates Microbial Community Assembly in Termite Mounds Across a Large Spatial Area

Background: Termites are ubiquitous insects in tropical and subtropical habitats, where they construct massive mounds from soil, their saliva and excreta. Termite mounds harbor an enormous amount of microbial inhabitants, which regulate multiple ecosystem functions such as mitigating methane emissions and increasing ecosystem resistance to climate change. However, we lack a mechanistic understanding about the role of termite mounds in modulating the microbial community assembly processes, which are essential to unravel the biological interactions of soil fauna and microorganisms, the major components of soil food webs. We conducted a large-scale survey across a >1500 km transect in northern Australia to investigate biogeographical patterns of bacterial and fungal community in 134 termite mounds and the relative importance of deterministic versus stochastic processes in microbial community assembly. Results: Microbial alpha (number of phylotypes) and beta (changes in bacterial and fungal community composition) signicantly differed between termite mounds and surrounding soils. Microbial communities in termite mounds exhibited a signicant distance-decay pattern, and fungal communities had a stronger distance-decay relationship (slope = -1.91) than bacteria (slope = -0.21). Based on the neutral community model (tness < 0.7) and normalized stochasticity ratio index (NST) with a value below the 50% boundary point, deterministic selection, rather than stochastic forces, predominated the microbial community assembly in termite mounds. Deterministic processes exhibited signicantly weaker impacts on bacteria (NST = 45.23%) than on fungi (NST = 33.72%), probably due to the wider habitat niche breadth and higher potential migration rate of bacteria. The abundance of antibiotic resistance genes (ARGs) was negatively correlated with bacterial/fungal biomass ratios, indicating that ARG content might be an important biotic factor that drove the biogeographic pattern of microbial communities in termite mounds. Conclusions:


Background
Soil microorganisms and fauna are pivotal components in soil food webs, and regulate fundamental processes such as biogeochemical nutrient cycling and plant growth promotions. Understanding the ecological and environmental drivers that govern the assembly of microbial communities is a central aim in microbial ecology [1,2]. Mounting evidence suggested that microbial communities can coexist in the environment with limiting resources and exhibit a signi cant biogeographic pattern from local to global scales [3,4]. Deterministic (niche-based theory) and stochastic (neutral theory) processes are commonly used to explain the assembly of species into communities [3,5]. The fundamental premise of niche-based theories is that ecological traits differ among species within a community (e.g., different nutrition requirements), which allows them to occupy different niches and to differentiate limited resources within the ecosystem [6,7]. Neutral theories assume that all individuals are ecologically equivalent with the same demographic rates [8], and dispersal limitation and stochasticity are the key processes that drive the community assembly [9]. After a long-standing debate, however, deterministic and stochastic processes are now known to be not mutually exclusive, but rather act concurrently to regulate the community assembly [5,10,11].
To date, the roles of abiotic drivers (e.g., pH, temperature, and nutrient content) in shaping the diversity and distribution of bacterial and fungal communities in terrestrial ecosystems have been well studied across regional to global scales [12][13][14]. Recent studies reported that the global distribution of soil faunal communities including earthworms, nematodes, and invertebrates, is primarily driven by vegetation and climate across many biomes [15][16][17]. However, biotic factors (e.g., soil faunal colony and activities) as a driving force of the diversity and assembly process of soil microbial communities remain largely unexplored at large spatial scales. Termites are among the most abundant groups of insects on Earth, and they play a vital role in the terrestrial ecosystem by consuming dead plants, woods and feces, driving nutrient cycling processes, and producing greenhouse gas methane [18][19][20]. Some termites construct enormous colonies ranging in size from a few hundred to several million individuals, and create unique nest structures ('mounds') from soil, saliva and excreta, which substantially alter soil physicochemical properties, increase soil nutrient availability, and enhance plant growth [21][22][23]. Recent studies demonstrated that termites produce globally signi cant amounts of methane, half of which is oxidized before emission by methanotrophic bacteria dwelling in the termite mound walls (Nauer et al., 2018). Termite mounds, as an 'island of fertility' widely distributed in savannas [24,25], harbour a broad range of soil microbial communities, but we lack the knowledge of microbiomes associated with termite mounds across large spatial scales. It is imperative to characterize the microbial assembly processes in ubiquitous termite mounds, which is essential to predict the role of termite mounds in regulating terrestrial ecosystem functions.
In this study, we conducted a large-scale survey across > 1500 km distance in northern Australia to identify the role of termite mounds in regulating the diversity and composition of soil bacterial and fungal communities (Fig. 1). We measured the abundance of hundreds of antibiotic resistance genes (ARGs) using a high-throughput quantitative PCR array [26], which has been recognized as a good proxy of biological interactions in the natural environment with minimal human disturbance [13]. Our study aimed to investigate the spatial patterns and community assembly mechanisms of bacteria and fungi in termite mounds and surrounding soils, and identify the relationship between ARGs and microbial communities in termite mounds. We hypothesized that deterministic selection would dominate the microbial community assembly in termite mounds, because termite activities create an unique environment through enhancing nutrient availability and altering soil properties.

Results
Microbial community compositions in termite mounds and surrounding matrix soils The microbial community compositions were predominated by the bacterial phyla Actinobacteria, Chloro exi, Proteobacteria, and the fungal phylum Ascomycota in both termite mounds and surrounding soils (Additional le 1: Figure S1), despite the signi cant differences in the relative abundances of some phyla between two habitats (Additional le 1: Figure S2). For example, Ascomycota, a common fungal phylum in soil, showed an extremely higher relative abundance in termite mounds (with an average of 96.21%) than in matrix soils (with an average of 88.31%). The bacterial phyla Actinobacteria and Proteobacteria were signi cantly more abundant in termite mounds (with an average of 82.22%) compared to matrix soils (with an average of 58.11%). Conversely, most of the less abundant phyla (relative abundance < 3%) such as Planctomycetes, Cyanobacteria and WPS2 were more abundant in surrounding matrix soils than termite mounds.
Termite mounds showed a signi cantly lower alpha-diversity (Shannon index) of both bacteria and fungi than surrounding matrix soils ( Fig. 2A, 2B). Principal coordinate analysis based on the Bray-Curtis distance indicated that the microbial community compositions (beta-diversity) of bacteria (Adonis, P < 0.01) and fungi (Adonis, P < 0.01) in termite mounds were signi cantly distinct from those in the matrix soils (Fig. 2C, 2D).A signi cantly positive correlation was observed for bacterial richness (r = 0.241, P = 0.005), but not for fungal richness (r = 0.136, P = 0.116), in surrounding matrix soils and termite mounds (Fig. 2E).

Distance-decay of community similarity
We estimated the distance decay relationship (DDR) for both bacteria and fungi in termite mounds and matrix soils. Although DDR was signi cant for both bacteria and fungi (P < 0.001), the tness values were remarkably low (R 2 < 0.1), indicating a weak decay of microbial community similarity with geographic distance. The slope of distance-decay curves, based on the OLS regression model, was steeper for fungi (slope = -1.91) than that for bacteria (slope = -0.21) in termite mounds (Fig. 3A, 3B). Similar distancedecay pattern for bacteria (slope = -0.79) and fungi (slope = -0.83) was observed in the matrix soils (Fig. 3C, 3D). The slope of distance-decay curves was steeper for bacteria in matrix soils than that in termite mounds, while fungi showed an opposite pattern (Fig. 3).
Microbial community assembly in termite mounds and matrix soils We calculated the community-level habitat niche breadth (Bcom) to explore the relative importance of deterministic and stochastic processes in microbial community assembly. A signi cantly higher Bcom value was observed for bacteria than fungi in both termite mounds and matrix soils (P < 0.001), and termite mounds had signi cantly higher Bcom values for both bacteria and fungi than matrix soils (Fig. 4A, 4C). We also tted the bacterial and fungal communities to the neutral community model. In termite mounds, the goodness of t was better for bacteria (65.42% of the variations explained) than fungi (32.85% of the variations explained) (Fig. 5). Bacteria exhibited a higher mitigation rate (m = 0.0372) than fungi (m = 0.0015), indicating that bacterial taxa were less limited by dispersal. In matrix soils, the goodness of t for bacteria was similar to termite mounds (64.39% of the variations explained), but the soil fungal communities cannot t to the neutral community model (Additional le 1: Figure S3).
We further employed the normalized stochasticity ratio (NST) to quantify the role of deterministic and stochastic processes in microbial community assembly (Fig. 4B, 4D). The NST value was below the 50% boundary point for both bacterial and fungal communities, suggesting that deterministic process played a more important role than stochasticity during the microbial assembly in termite mounds and matrix soils. Additionally, a signi cantly higher NST value was observed in bacterial communities than in fungal communities (P < 0.001). The NST value of bacteria showed no signi cant difference (P = 0.655) between termite mounds (with an average of 45.23%) and matrix soils (with an average of 46.37%). The NST value of fungi in termite mounds (with an average of 33.72%) was signi cantly (P = 0.017) higher than that of matrix soils (with an average of 23.46%).

Biological interactions and antibiotic resistance genes
The Bray-Curtis distance-based Mantel test (P < 0.0001, r = 0.739) together with Procrustes analysis (P < 0.001, M 2 = 0.234) suggested that there was a signi cant correlation between bacterial and fungal community compositions in termite mounds (Fig. 6A). Furthermore, we found that the abundance of ARGs was negatively correlated with bacterial/fungal biomass ratios in termite mounds (Fig. 6B), suggesting that ARGs could be an important indicator to show the inter-kingdom biological interactions (e.g., competition between bacteria and fungi). For the abiotic determinants of microbial assembly in termite mounds, edaphic factors such as soil pH and dissolved organic nitrogen (DON) were strongly associated with microbial taxonomic diversity (richness and Shannon diversity index), and more importantly, bacteria and fungi showed distinct preference for soil pH. Climate factors including mean annual precipitation (MAP), mean annual temperature (MAT) and aridity index (AI) also showed a strong impact on bacterial and fungal diversity (Additional le 1: Figure S4).

Discussion
Deterministic process dominates microbial community assembly in termite mounds The biogeographic pattern of microbial communities has received little attention so far in ubiquitous termite mounds, an unique environment for microorganisms in terms of nutrient availability and soil properties. The present study sheds lights into the underlying processes and mechanisms of microbial community assembly in the habitat strongly affected by soil faunal activity. We found that the overall pattern of microbial communities in termite mounds was distinct from their surrounding matrix soils, and a lower microbial diversity was observed in termite mounds, probably due to the strong selection pressure imposed by termite activities. A signi cant distance decay pattern for both bacterial and fungal communities was observed in termite mounds across large spatial scales. The tness values (R 2 < 0.1) of the DDR were substantially lower than those reported in other natural soil ecosystems [27,28], suggesting that microbial communities in termite mounds exhibit less apparent spatial patterns across large scale.
According to Hubbell's neutral theory, community similarity was predicted to decrease with increasing spatial distance due to dispersal limitation [9]. Nevertheless, it is di cult to disentangle the relative importance of deterministic and stochastic processes by solely analyzing DDR.
The differences in environmental factors, such as pH, MAT, MAP and vegetation, also likely increase with increasing geographical distance, which play a key role in in structuring microbial communities [29,30].
In this case, the environmental changes induced shifts in microbial community compositions are based on the taxonomic niche differentiation i.e. the deterministic process [31]. From the neutral community model, we observed a lower tness value (R 2 ) and mitigation rate (m) compared with previous studies [10,32], in which stochastic processes were dominant in shaping microeukaryotic and zebra sh gut microbial community assembly. Thus, our results indicate that stochastic processes may play a less important role in microbial community assembly in termite mounds, which was supported by the NST, a general mathematical framework for quantitatively assessing the in uences of stochastic and deterministic processes in microbial community assembly [33]. The NST value below the boundary (50%) further provided evidence that the niche based deterministic process dominates microbial community assembly in termite mounds.

Differences in community assembly between bacteria and fungi
In the present study, we demonstrated that overall the niche-based deterministic processes determine the microbial community assembly in termite mounds and matrix soils. An interesting nding is that the NST value of bacteria is signi cantly higher than that of fungi (Fig. 3B), suggesting that stochastic processes contributed more to bacterial communities than to fungal communities. To explain this observation, we quanti ed the community-level habitat niche breadth, and found that bacteria had a signi cantly wider niche breadth than fungi. A previous study reported that organisms with wider niche breadth might have greater metabolic plasticity and be less in uenced by deterministic processes [34]. In addition, we observed that the mitigation rate of bacterial communities was approximately 25 times higher than fungal communities based on the NCM (Fig. 4), which is an important factor for stochastic processes and could partially explain why the slope of the fungal distance pattern was steeper than that of bacteria.
Other explanations could be related to differential responses of bacteria and fungi to environmental factors, which has been suggested as a main reason for the opposing global biogeographic trends for bacteria and fungi [13]. For instance, we observed a signi cant positive relationship between fungal diversity and aridity index, in contrast, no signi cant correlation was identi ed for bacterial diversity. Microbial dormancy is suggested to promote biodiversity within microbial communities by allowing competing organisms to partition resources across time, rather than space [35,36]. It has been suggested that less than 10% of bacterial cells in a given community are in an active stage at any time [37]. In fact, fungi do have the ability to form dormant stages (e.g. resistant spores). The difference in dormancy strategies may affect the community niche breadth and assembly, which was supported by a previous study that protist communities are more responsive to species sorting than bacterial communities, because protists have a more limited tendency to enter dormancy [38].
The NST value of bacteria showed no signi cant difference between termite mounds and matrix soils. In contrast, the NST value of fungi was much lower in matrix soils compared with termite mounds, indicating that termite activity may have a stronger impact on fungal community assembly than bacterial community assembly. To disentangle the mechanisms underlining the different impacts of termite on soil bacterial and fungal communities assembly, a comprehensive exploration of the termites surface and gut associated microbiota is necessary in future studies [39].

Biological interactions of microbial communities in termite mounds
Environmental ltering plays an important role in microbial community compositions in termite mounds, as revealed by the signi cant relationship between abiotic factors (including climatic and edaphic factors) and microbial diversity, which is consistent with the global topsoil microbial communities [13]. Procrustes test revealed a signi cant correlation between bacterial and fungal community compositions, indicating that the inter-kingdom biological interactions may affect the microbial composition and biogeographic pattern in termite mounds (Fig. 5A). However, less effort has been devoted to quanti cationally depict the interactions between bacteria and fungi. In the present work, we found that a negative correlation between the ARG abundance and bacterial/fungal biomass ratio in termite mounds (Fig. 5B), and similar results were observed in global topsoil microbiome [13]. This could be due to the competition for resources that affects the biomass of fungi and bacteria [40], and bacteria surviving fungal antagonism are enriched in ARGs. Therefore, we infer that inter-kingdom antagonism, as re ected in the association of bacterial ARGs with bacterial/fungal biomass ratios, could be a potential driver in structuring microbial communities in termite mounds.
A few potential caveats in the interpretation of our ndings should be noted, regarding the biological interactions. Trophic-level interactions, e.g. predation and mutualisms, are not considered in the present study, but they play a certain role in maintenance of species diversity [35]. In addition, the migration/movement of termites may increase the termite-associated microbial dispersal and affect the relative importance of deterministic and stochastic processes in shaping the microbial community assembly in termite mounds. These limitations together with the missing environmental variables could be the reason that a large proportion of variation remains unexplained for both bacterial and fungal communities in the variation partitioning (Additional le 1: Figure S5). Actually, addressing the effect of all potentially biotic and abiotic variables is impractical, especially for a broad-scale eld survey and this is consistent with previous studies that 50% − 90% variation remains unexplained [28,38,41].

Conclusions And Implications
Microorganisms play key role in the ecological signi cance of termite mounds. Our work provides novel evidence that microbial community assembly in termite mounds was governed to a great extent by deterministic relative to stochastic processes, which is the rst step to predict the ecological activities and processes in termite mounds under the context of global changes. The difference in the stochasticity between bacteria and fungi may be attributed to the differences in habitat niche breadth, potential mitigation rate, responses to environmental variations and dormancy strategies, which highlight the importance of considering organism characteristics in differentiating the role of stochastic and deterministic assembly processes in microbial communities. In addition to the climatic and edaphic factors, we suggested that the content of ARG could be an indicator to re ect the inter-kingdom biotic interactions in termite mound microbial communities. in this sampling campaign. All samples were sieved < 2 mm mesh and divided into two fractions. One fraction was frozen at -80 ºC for molecular analyses of bacteria and fungi, and the other was stored at 4 ºC for soil physicochemical analyses.

Climatic data collection
Climatic attributes, including MAP and MAT of each sampling location, were obtained from the WorldClim database version 2 (http://www.worldclim.org/) based on the spatial geographical coordinates [42]. We calculated the aridity index (AI, the ratio of precipitation and evapotranspiration) of each location using the Global Aridity Index and Potential Evapotranspiration in Climate Database v2 (https://cgiarcsi.community/). Extracted data were processed with the 'raster' and 'sp' packages in R [43,44].

Soil physicochemical analysis
Standard methods were used to characterize soil physicochemical properties. Brie y, soil pH was measured in a 1: 2.5 mass: volume of soil and water suspension with a pH meter (Thermo Scienti c Inc., Waltham MA, US). Dissolved organic carbon (DOC) and nitrogen (DON) were extracted with MilliQ water and measured using a TOC analyzer (Shimadzu, Kyoto, Japan). Total carbon (TC) and nitrogen (TN) were measured using the Dumas combustion method on LECO FP628 analyzer (LECO Corp., MI, USA). Nitrate The bacterial 16S rRNA gene and fungal ITS region were ampli ed with the primers 515FmodF/806RmodR [45] and ITS1F/ITS2R [46], respectively, on the CFX96 Touch™ PCR Detection System (Bio-Rad, Herculers, USA). Amplicons of bacteria and fungi were sequenced on an Illumina MiSeq PE300 platform (Illumina Inc., CA, USA) at the Majorbio Bio-pharm Technology Co., Ltd (Shanghai, China). The obtained raw sequences were ltered for quality control as previously described [47]. Chimeric sequence was checked using the 'USEARCH' package with the UCHIME algorithm [48]. Sequences were split into operational taxonomic units (OTUs) at 97% sequence similarity level using the UPARSE pipeline [49]. Singletons, chloroplast and mitochondrial sequences were discarded from the nal dataset. Representative sequences were assigned to taxonomic lineages using the RDP classi er against the SILVA database for bacteria [50] and UNITE database for fungi [51].

Quanti cation of ARGs
High-throughput quantitative PCR was performed to quantify the abundance of ARGs using the Wafergen SmartChip Real-time PCR system (Wafergen Inc., CA, USA). A total of 296 primer sets were used to interrogate the extracted DNA from termite mounds and surrounding soils [52]. These primer sets targeted resistance genes for all major classes of antibiotics (285 primer sets), transposase genes (eight primer sets), one class 1 integron-integrase gene, one clinical class 1 integron-integrase gene and the 16S rRNA gene [53][54][55].

Statistical analysis
All statistical analyses were performed in the R platform. To compare the microbial community composition between termite mounds and their surrounding matrix soils, we calculated the alphadiversity (Shannon index) and beta-diversity based on the Bray-Curtis distances. To determine the signi cance in the difference in microbial community compositions between termite mounds and surrounding soils, PERMANOVA analysis was performed with the 'Adonis' function in the 'vegan' package [56]. The signi cant difference in the relative abundances of speci c microbial taxa between termite mounds and soils was identi ed using one-way ANOVA. Spearman's rank correlation test was conducted to calculate the relationships between the main abiotic factors and the diversity of bacteria and fungi. To estimate the relationships between the bacterial and fungal communities, we performed the Procrustes analysis and Spearman's rank correlation between the abundance of ARGs and the bacterial/fungal abundance ratios [13]. The pairwise geographic distance between sampling sites was calculated using the 'sp' package [43] in R based on the longitude and latitude coordinates of each sampling site. The distance-decay rate of the microbial communities was calculated as the slopes of ordinary least-squares regressions between geographic distance and community similarity (1-dissimilarity of the Bray-Curtis matrices).
Habitat niche breadth is a crucial trait that in uences the relative importance of deterministic and stochastic processes in shaping the microbial community assembly [34]. We calculated Levins' niche breadth index (B) for both bacteria and fungi according to the formula [34,57]: Where B j represents the habitat niche breadth of OTU j in a metacommunity; N is the total number of communities in each metacommunity; P ij is the proportion of OTU j in community i. A high B-value for a given OTU indicates that the OTU occurs widely and has a wide habitat niche breadth. Habitat niche breadth at the community level was estimated as the average B-values from all OTUs in a single community (Bcom) [38]. A microbial group with a wider niche breadth is more metabolically exible at the community level [41].
A neutral community model (NCM) was used to predict the potential importance of stochastic processes in community assembly by determining the relationships between the detection frequency of microbial taxa in a set of local communities and their relative abundance across the wider metacommunity [58]. In general, the neutral community model assumes that highly abundant taxa in the metacommunity are more likely to be dispersed by chance and widespread across the metacommunity. Conversely, rare taxa are more likely to be lost in local communities due to ecological drift (i.e., the stochastic loss and replacement of individuals). In this model, the estimated migration rate (m) is a parameter to evaluate the probability that a random loss of an individual in a local community would be replaced by dispersal from the metacommunity, and therefore is a measure of dispersal limitation [32]. Higher m values indicate that microbial communities have a higher immigration rate and are less limited by dispersal. R 2 value indicates the t of the parameter based on nonlinear least squares tting. Calculation of 95% con dence intervals around all tting statistics was done by bootstrapping with 1000 bootstrap replicates. We compared the estimated migration rate and model tness for bacteria and fungi. All computations were performed in R with the code provided by Burns et al. [32] We used a general mathematical framework to quantify the relative importance of deterministic and stochastic processes in community assembly [33]. An index, normalized stochasticity ratio (NST), was used with 50% as the boundary point between more deterministic (NST < 50%) and more stochastic (NST > 50%) assembly. In this study, we used the Bray-Curtis similarity matrices (in line with previous models) coupled with the null models of xed taxa richness and proportional taxa occurrence frequency. We used the bootstrapping method (1000 randomization times) to test the distribution of NST in bacterial and fungal groups and the signi cant difference between groups. The analysis was conducted in R with the 'NST' package [33]. Figure 1 Sampling locations of termite mounds in Northern Territory and Queensland in Australia and the photos of termite mounds.

Figure 2
An overview of the microbial communities in termite mounds and surrounding matrix soils. A, B, bacterial and fungal diversity in termite mounds and matrix soils. C, D, Bacterial and fungal community betadiversity visualized using principal coordinate analysis (PCoA) based on the Bray-Curtis similarity. E, microbial diversity correlation analysis between termite mounds and matrix soils. between bacterial and fungal communities in termite mounds (B) and matrix soils (D). Asterisks denote signi cant correlation (***P < 0.001).

Figure 5
Fit of the neutral community model (NCM) showing the OTU predicted occurrence frequencies versus the relative abundance in termite mounds. The blue and red solid lines indicate the best t to the Sloan's neutral model and the dashed lines represent 95% con dence intervals around the model prediction. R2 and m values indicate the goodness of t to the neutral model and the estimated migration rate, respectively.