Ecological processes controlling overall mangrove microbial community assembly
The Sloan neutral model showed that the overall archaeal and bacterial communities were driven predominantly by stochastic processes, with the R2 value were 0.775 and 0.831 for overall archaea and bacteria, respectively, and larger explained community variance for bacteria was observed (Fig. 1a,b). The estimated migration rate (m), a measure of the influence of dispersal on community composition, were higher in overall bacteria (0.520) than in overall archaea (0.494). Both the R2 and estimated migration rate (m) indicated a stronger effect of dispersal limitation on archaea than bacteria. The community-level habitat niche breadths (Bcom) were estimated to reveal the contributions of deterministic and stochastic processes. Niche breadth (Bcom) of bacteria was much greater than that of archaea in overall mangrove ecosystem and four different types of vegetation zones (Tukey’s HSD test, P < 0.001, Fig. 1c).
The null model indicated that the differential action of ecological processes may promote different biogeographic patterns in archaeal and bacterial assemblages. However, the stochastic processes (sum of dispersal limitation, homogenizing dispersal, and drift) explained a higher proportion of the archaeal and bacterial communities (including the overall mangrove ecosystem and four different types of vegetation zones) variation than deterministic processes (Fig. 2), which supported the results of the NCM. These results suggested that stochasticity was more important than determinism in influencing mangrove microbial community. The stochastic processes accounted for 79% and 87% of the community assembly in the overall archaea and bacteria, respectively, and bacteria were more controlled by stochasticity than archaea (Fig. 2a,b). Drift and homogenizing dispersal were the most important process, accounting for 72% and 73% of the archaeal and bacterial communities variation, respectively. The ratio of stochasticity and determinism also confirmed that the dominance of stochasticity relative to determinism, and the ratio values varied more sharply in archaea (ranged from 3.68 to 65) than bacteria (ranged from 6.86 to 32) (Fig. 2c,d).
Spartina alterniflora invasion changed mangrove microbial community assembly
To determine the effect of Spartina alterniflora invasion on the archaeal and bacterial community assembly, the neutral community model (NCM) and null model were also used with the datasets from four different vegetation zones for archaea and bacteria. The NCM showed that S. alterniflora invasion changed the relative contribution of ecological processes controlling microbial community assembly. The values of R2 and immigration rate (m) distribution showed same pattern for archaeal and bacterial subcommunities: Mudflat (R2 = 0.608 and m = 0.579 for archaea; R2 = 0.718 and m = 0.617 for bacteria) > Ecotone (R2 = 0.585 and m = 0.574 for archaea; R2 = 0.714 and m = 0.622 for bacteria) > Cordgrass (R2 = 0.559 and m = 0.526 for archaea; R2 = 0.701 and m = 0.638 for bacteria) > Mangrove (R2 = 0.528 and m = 0.413 for archaea; R2 = 0.671 and m = 0.495 for bacteria). All of the four bacterial subcommunities fitted better to NCM than four archaeal subcommunities(Fig. 1a,b). The null model also suggested that invasion changed the relative contribution of microbial community assembly with different degree for archaea and bacteria. In the archaea, stochastic processes explained large subcommunities variation among four different vegetation zones with the following pattern: mangrove (98%) > ecotone (95%) > cordgrass (89%) > mudflat (85%). The bacteria showed the following pattern: cordgrass (97%) > mangrove (92%) > ecotone (91%) > mudflat (87%) (Fig. 2a,b).
Overall mangrove ecosystem’ microbial network co-existence patterns
The correlation-based network consisted of 212 nodes (OTUs) with 1083 edges (correlations) for the archaea, and 277 nodes with 3721 edges for the bacteria (Table S2). Overall, taxa tended to co-occur (positive correlations, yellow lines) rather than co-exclude (negative correlations, blue lines); positive correlations accounted for 93.44% and 77.18% of the potential interactions in archaeal and bacterial networks, respectively, whereas negative correlations were 6.56% and 22.82% interactions for archaeal and bacterial co-existence patterns (Fig. 3). When considering all correlations, the links between bacteria were more complex than those between archaea (Fig. 3, Table S2), indicating that potential interactions are stronger in bacterial networks. For the archaeal network, Crenarchaeota (39.15%), Nanoarchaeaeota (35.85%), Euryarchaeaeota (8.96%), Thaumarchaeaeota (8.02%), and Asgardaeota (3.77%) mainly occupied the nodes (Fig. 3a). Nodes in bacterial network mainly belonged to Proteobacteria (55.96%), Chloroflexi (14.08%), Bacteroidetes (9.39%), Actinobacteria (6.14%), Nitrospirea (3.25%), and Epsilonbacteraeota (2.89%) (Fig. 3c). Furthermore, a module is defined as a group of OTUs that are linked more tightly together. Here, both the entire archaeal and bacterial networks were clearly parsed into 6 major modules, of which modules I, II, and III accounted for 26.42%, 24.3%, and 20.28% of the whole archaeal network, respectively (Fig. 3b), and modules I and II accounted for 27.8% and 24.55% of the whole bacterial network, respectively (Fig. 3d).
The integrated network degrees were distributed according to a power-law distribution in both archaea and bacteria, indicating a scale-free distribution and non-random co-occurrence pattern (Fig. S5). We calculated a set of network-level topological features, and found that values of the degree, closeness centrality, and eigenvector centrality in bacteria were significantly higher than those in archaea (Fig. 4a, Table S2). Furthermore, the average clustering coefficients were higher in the bacterial network than that of archaea (including the overall and four different types’ vegetation zones of archaea and bacteria), suggesting that bacterial OTUs were more interconnected (Table S2). The average path length and diameter were lower in the bacterial network, revealing closer relationships among bacterial communities (including the overall and four different types’ vegetation zones of archaea and bacteria). Random networks were generated with the same nodes and edges in each compartment to confirm that the empirical networks were non-random. Details describing the constructed co-occurrence networks can be found in Table S2.
Zi-Pi plot showed that Woesearchaeia and Proteobacteria phyla were the most prominent keystone taxa for archaea and bacteria, respectively. In the co-occurrence networks, 5 archaeal OTUs and 12 bacterial OTUs were defined as keystone taxa, and Woesearchaeia and Proteobacteria phyla accounted for 80% and 50% of all module hubs and connectors (Fig. 4b; Table S3). In archaea, the putative keystone species include taxa from the classes Woesearchaeia (Nanoarchaeaeota, 4 OTUs) and Bathyachaeia (Crenarchaeota, 1 OTU). In bacteria, the putative keystone taxa include taxa from the classes Gammaproteobacteria (Proteobacteria, 3 OTUs), Actinobacteria (Actinobacteria, 1 OTU), Alphaproteobacteria (Proteobacteria, 3 OTUs), Anaerolineae (Chloroflexi, 2 OTUs), Bacteroidia (Bacteroidetes, 2 OTUs), and Campylobacteria (Epsilonbacteraeota, 1 OTU). Keystone taxa spanned a range of relative abundances (0.06 to 1.43% for archaea and 0.05 to 0.38% for bacteria). Over half of the keystone taxa (9 of 17 OTUs for both archaea and bacteria) had low relative abundance (0.05 to 0.10%). All of the 17 OTUs were conditionally rare taxa (Table S3).
Spartina alterniflora invasion influenced microbial network co-existence patterns
To identify the effect of S. alterniflora invasion on potential microbe-microbe interactions, we constructed four archaeal and bacterial co-occurrence networks among four different types of vegetation zones (Fig. 5). Invasion changed microbial network co-existence patterns, since the network complexity, such as degree (i.e. links per node in the network) and connections (links) among four different vegetation zones were different. Moreover, the effect of invasion on the archaeal networks differed profoundly from bacterial networks. Multiple network topological metrics (e.g. links, degree, and average clustering coefficient) consistently supported the different effect of invasion on the archaeal and bacterial co-occurrence patterns (Table S2). Although the selected network size (nodes) were similar among distinct vegetation zones, the connectivity (links) of archaeal and bacterial networks were different. In the archaeal networks, the subcommunity in mudflat formed largest networks’ connections (links), followed by ecotone, cordgrass, and mangrove vegetation zones. Whereas in the bacterial networks, the network complexity showed the following trend: mangrove > ecotone > mudflat > cordgrass. The complexity of the networks was also reflected by the degree, which showed the same trend with links in the archaeal and bacterial networks (Fig. 5, Table S2). Overall, positive correlations accounted for 70–80% and 64%-75% of the potential interactions in archaeal and bacterial networks among four different vegetation zones, respectively, which were higher than negative correlations.
We compared unique node-level topological features of four subcommunities from the different vegetation zones. The network topological parameters such as betweenness centrality and closeness centrality did not differ significantly among four different vegetation zones of archaeal and bacterial subcommunities (Fig S6). However, the degree value in the mudflat was significantly highest among four archaeal subcommunities. And in bacterial subcommunities, the ecotone showed the highest degree. Furthermore, eigenvector centrality varied significantly among different vegetations’ sediment (Tukey’s HSD test, P < 0.001, Fig. S6).
We found that S. alterniflora invasion changed putative keystone taxa of archaeal and bacterial subcommunities (Fig. S7-S8; Table S4-S5). In the archaeal subcommunities, a total of 82 OTUs were identified as keystone species, including the members from mangrove (28 OTUs), ecotone (26 OTUs), cordgrass (13 OTUs), and mudflat (15 OTUs) (Fig. S7; Table S4). Furthermore, in the bacterial subcommunities, 66 OTUs were considered as keystone taxa including 11 OTUs in mangrove, 11 OTUs in ecotone, 15 OTUs in cordgrass, and 29 OTUs in mudflat (Fig. S8; Table S5). Almost all of the keystone taxa were module hubs and connectors, and only one network hub was detected in all of the constructed archaeal and bacterial networks. Among four different vegetation zones, the most prominent keystone taxa in the archaeal networks were from the classes Woesearchaeia and Bathyarchaeia, and the major keystone taxa in the bacterial networks were Proteobacteria. Most of the keystone taxa were conditionally rare taxa and always rare taxa (Table S4-S5). Interestingly, 5 archaeal OTUs from Woesearchaeia and Bathyarchaeia (i.e. OTU_20 and OTU_131) and 4 bacterial OTUs from Proteobacteria and Gemmatimonadetes (i.e. OTU_18 and OTU_156) were simultaneously detected in different vegetation zones, indicating that these OTUs are important in different vegetations’ ecosystem (Table S4, S5).
To determine the effect of S. alterniflora invasion on the robustness of the archaeal and bacterial networks, a natural connectivity analysis was carried out among four different types of vegetation zones. In the archaeal subcommunities, the natural connectivity in the mudflat network was higher than that of mangrove, ecotone, and cordgrass vegetation zones, whereas in the bacterial subcommunities, we found the greatest natural connectivity in the mangrove vegetation zone, followed by ecotone, mudflat, and cordgrass, indicating that the mudflat archaeal network and mangrove bacterial network was the most stable (Fig S5).