Spatiotemporal Variations Dictate Stable and Resilient Microbiome Interactions in the Chesapeake Bay

Background: Annually reoccurring microbial populations with strong spatial and temporal variations have been identied in estuarine environments, especially in those with long residence time such as the Chesapeake Bay (CB). However, it is unclear how microbial taxa interact with each other (e.g., mutualistic and competitive interactions) and how these interactions respond to their surrounding environments. Specically, there is a lack of understanding of how these interactions inuence microbiome population dynamics, and its adaptability and resilience to estuarine gradients. Results: Here, we constructed co-occurrence networks on prokaryotic microbial communities in the Bay, which included seasonal samples from seven spatial stations along the salinity gradients for three consecutive years. Our results showed that spatiotemporal variations of planktonic microbiomes promoted differentiations of the characteristics and stability of prokaryotic microbial networks in the CB estuary. Prokaryotic microbial networks are more stable seasonally than spatially, and microbes were more strongly connected during warm season compared to the associations during cold season. In addition, microbial interactions were more stable in the lower Bay (ocean side) than those in the upper Bay (freshwater side). Interestingly, compared to the abundant groups, the rare taxa such as SAR116 clade, SAR11 clade III, and OM182 clade contributed greatly to the stability and resilience of prokaryotic microbial interactions in the Bay. Modularity and cluster structures of microbial networks varied spatiotemporally, which provided valuable insights into the ‘small world’ (a group of more interconnected species), network stability, and habitat partitioning/preferences. Multivariate regression tree (MRT) analysis and Piecewise structural equation modeling (SEM) indicated that temperature, salinity and total suspended substances directly or indirectly (through nutrient availability, particulate carbon and Chl a) affected the distribution and associations of microbial groups, such as Actinobacteria, Bacteroidetes, Cyanobacteria, Planctomycetes, Proteobacteria, and Verrucomicrobia. Conclusion: Our results shed light on how spatiotemporal variations alter the nature and stability of prokaryotic microbial networks in the estuarine structure equation model (SEM) unraveled that spatiotemporal variations (i.e. temperature, salinity, TSS, nutrient availability) are the primary driving factors for microbiome structures and interactions, and these environmental gradients manipulate the stability and tolerance/susceptibility of estuarine microbiota. This work shows, for the rst time, dynamics of microbiome networks and their responses to environmental gradients in a typical estuarine environment with long residence time. Given such dynamic environmental gradients, stable and persistent microbiome networks suggest that biotic interactions play a central role in maintaining integrity and resilience of estuarine microbiomes.


Background
Planktonic microbiomes comprise both free-living organisms and those attached to particles, which is typical and visible in water columns of estuaries. The planktonic microbiomes in estuaries are one of the most active microbial communities [1] and they contribute primarily to the most productive environments on the planet [2]. For example, estuaries have a high CO 2 ux (~ 0.25 Pg C y − 1 ) between water and air, which is largely supported by the process of microbial decomposition and carbon xation [3]. These microbiomes drive the estuarine biogeochemical processes of the elements for life [4]. They are powering the cycles of nutrients (e.g., carbon, nitrogen, phosphorus, and sulfur), with important impact on the composition of greenhouse gases in the atmosphere, the formation of algal blooms, and the integrity of estuarine ecosystems [5][6][7]. In fact, microbiomes are the foundation for the estuarine food chains and food webs [8][9][10], and their composition/distribution are important in the balance and stability of the entire estuary ecosystem.
Estuaries harbor a tremendous diversity of microbes. Due to the strong temporal and spatial gradients and surrounding land uses, the composition and distribution of the estuarine microbiomes are largely affected by human activities and climate/environmental changes [11][12][13]. Many studies have shown that spatiotemporal environmental variations enrich certain microbial taxa to dominate in estuaries, including the Chesapeake Bay [14], Delaware Bay [15], Sacramento-San Joaquin River Delta [16], Pearl Estuary [17], and Columbia River estuary [18]. In response to strong physical, chemical and biological gradients, estuarine microbiomes exhibit pronounced uctuations in production, biomass [19,20], and community composition [11,21,22]. It has also been shown that different microbial groups respond differently to spatiotemporal variations [11,23]. These earlier studies provide important insights into the impacts of spatiotemporal variations and other disturbances (such as anthropogenic pressures) on estuarine microbial community structure and dynamics. More interestingly, the spatial and temporal succession patterns are repeatable and predictable on an annual base [24][25][26][27][28], suggesting the microbial population dynamics are closely interrelating to their ambient aquatic environments. Annual selection pressure driven by environmental forcing has, through the induction of recurrent patterns in resource availability, predator-prey dynamics and microbial interactions, allowed for the assemblage of largely stable and resilient microbial communities [28]. Annually reoccurring patterns indicate that different microbiomes indeed have distinct niches with limited redundancy, otherwise different combinations of microbiomes would appear under the same conditions and prediction can be di cult [26]. These annually reoccurring assemblages and structure of microbiomes suggest their potential ecosystem functions (e.g., autotrophic or heterotrophic microbial production) are specialized and annually repeatable as well [24].
Functional traits of microbiomes are products of multiple populations within communities rather than those of a single population [29], and therefore interactions among microbial taxa are critical to maintain ecosystem integrity and microbiome functionality. Different species or populations interact with each other to form complicated networks through various types of interactions, such as predation, competition and mutualism [30,31]. Theoretical studies showed that communities in which a large proportion of members connected through positive links (i.e. positive correlations) are deemed to be unstable; in such communities, members may respond in tandem to environmental uctuations, resulting in positive feedback and co-oscillation [31]. In contrast, ecological networks with compartmentalization and presence of negative interactions could increase the stability of networks under disturbances [31][32][33]. For instance, high proportion of negative links could better balance the asynchronous dynamics and therefore stabilize co-oscillation in communities and promote stability of networks [31]. Further, modeling studies show that increasing strength of a few key interactions within a food web can destabilize trophic cascades, such as the "gatekeepers" [34], as removal of in uencers causes a network to fragmentation [35]. In general, microbial interactions (networks) play critical roles in maintaining community states, ecological niches and function distribution in the context of the microbiome [36].
Co-occurrence networks can reveal information on associations within microbiomes and stability of whole communities [31,36,37]. It has been increasingly used to infer microbial interactions [30,38] in soils [39,40], oceans [41,42], coastal waters [43], lakes [44,45], rivers [46], and even in metabolic modeling [47] and genomic surveys [48]. These correlation-based networks show important details of community interaction rules re ecting ecological processes such as cooperation and habitat partitioning and could represent mathematical associations among different microbial groups [30,37]. Despite cooccurrence networks have been researched in plenty of earlier studies, they rarely focused on the co-exist microbial networks and how these networks respond to spatiotemporal variations in estuarine environments such as the Chesapeake Bay (CB), the biggest estuary in North America with long residence time (up to 9 months) [49]. Taking into account the dynamic environmental gradients in CB such as temporal variations, freshwater runoffs and ocean water intrusion, microbial interactions and their stability/resilience to environmental changes are critical in further understanding the CB ecosystem. As expected, not only the spatiotemporal variations have potentials to reorganize networks of interactions between co-existing estuarine microbial taxa, the characteristics of these networks themselves can also determine the resilience to environmental disturbances (such as agriculture and urban development). Nevertheless, there are still few detailed studies on microbial networks and their responses to environmental changes in typical estuarine environments [50][51][52][53], and we need to address this important knowledge gap to deepen our understanding of estuarine microbiome ecology.
Recently, we have characterized the microbial community structures across both temporal and spatial scales in the CB, where planktonic microbiomes were collected across seven sampling stations (along spatial and salinity gradients) in four seasons over three consecutive years [54]. In this study, we constructed co-occurrence networks based on the detailed 16S rRNA gene high-throughput sequences. Signi cant and strong correlations (including both positive and negative correlations) were included in co-occurrence network analysis [55]. Due to previously described predominant seasonal and spatial variations in prokaryotic microbiome structure, networks from different seasons (temporal) and locations (spatial) were also constructed. In order to test the susceptibility and stability of co-occurrence networks, we assessed responsiveness of network fragmentation to removal of signi cant nodes (i.e., with highest betweenness centrality). Further, primary environmental drivers for microbial interactions in the estuary were tested with multivariate regression tree analysis (MRT, [56]). Lastly, pathways that may explain how environmental gradients contribute to shaping estuarine microbiomes and their interactions were identi ed and quanti ed with piecewise structural equation model (SEM, [57]). Based on all these analyses, quantitative interactions of environmental drivers and estuarine microbiome structures (phylum level) were proposed. This study provides the rst snapshot on CB microbiome networks, the stability and resilience, as well as their quantitative responses to environmental gradients, which are critical in improving our knowledge and understanding of the ecology of estuarine microbiomes.

Sample collection and characterization
Detailed description of sampling and environmental measurements have been described in previous studies [14,54]. Brie y, surface water samples were collected from the CB at seven stations along the middle axis in February/March, May/June, August, and October from 2003 to 2005 (Fig. S1). Total 500 mL water samples (below 2 m) were taken at each station and ltered immediately through 0.2 µm Millipore polycarbonate lters (Millipore Corporation, Billerica, MA, USA). The lters were stored at -20ºC prior to DNA (deoxyribonucleic acid) extraction. Water temperature and salinity were recorded on board with a Sea-Bird 911 CTD (conductivity temperature depth, Sea-Bird Electronics, Washington, USA). Based on water temperature [14,58], the samples collected in the Bay were divided into four seasons: winter (February and March), spring (May and June), summer (August) and autumn (October). Other water abiotic data were obtained from the Chesapeake Bay Program's (CBP) Water Quality Database (https://www.chesapeakebay.net/what/downloads/cbp_water_quality_database_1984_present), including total organic nitrogen (TON), ammonium, nitrate, Chl α (Chlorophyll a), dissolved organic phosphorus (DOP), orthophosphate phosphorus (OP/phosphate), total suspended solids (TSS), particulate carbon (PC), and turbidity (measured as Secchi depth). Although the CBP stations do not have the exact same coordinates as our stations, they are reasonably close to our sampling sites (within 2 to 3 km).

DNA extraction and sequencing analysis
Environmental DNA extraction followed the protocol described previously [59]. Dried environmental DNA pellets were lyophilized and archived in -80ºC for long time storage. The V4 region of the 16S rRNA (ribosomal ribonucleic acid) genes were ampli ed using the primers 515f (5′-GTGYCAGCMGCCGCGGTAA-3′) and 806r (5′-GGACTACNVGGGTWTCTAAT-3'), and sequenced on an Illumina Nova6000 platform (paired-end 250-bp mode) following the manufacturer's guidelines. High-throughput sequences were processed and analyzed with QIIME (Quantitative Insights Into Microbial Ecology) 2_2019.1 [60], and amplicon sequence variants (ASVs) were rari ed to 100,000 sequences per sample for downstream analyses [54]. Raw sequence data is available at the GenBank database under the accession number SRX6973110-SRX6973185.

Co-occurrence network analysis
Three categories of co-occurrence networks were constructed: (1) networks across all samples to explore overall interactions at spatiotemporal scales; (2) networks in four different seasons based on water temperature; and (3) networks in upper and lower Bay. To explore microbial interactions across both spatial and temporal scales, we split the samples to four seasons (based on water temperature) or two locations (upper Bay and lower Bay, based on salinity) for co-occurrence network analysis. Salinities were signi cantly different (P < 0.01) between upper Bay (samples from station 908, 845, and 818) and lower Bay (samples from station 707, 724, and 744). ASVs a liated with family level were used to construct all the co-occurrence networks. To avoid potentially erroneous sequences and improve interpretability of the dataset, we ltered out the ASVs (1) presented in fewer than three samples, or (2) with summed relative abundance less than 0.1% in each particular network inference [46]. All network constructions were performed in R (version 3.6.1) using the package 'fdrtool' and 'igraph' [55]. Network construction codes were adapted from GitHub (https://github.com/ryanjw/co-occurrence). The false discovery rate was estimated and corrected by the package 'fdrtool'. Only microbial families with statistically signi cant (P < 0.01, Q-value < 0.05) and robust (Spearman's correlation coe cient > |0.6|) correlations [40] were included in network analyses. In order to test if networks were signi cantly clustered, random networks were also constructed and compared following RJ Williams, A Howe and KS Hofmockel [55].
Network visualization and topological analysis were carried out in Gephi (version 0.9.2) [61]. These networks made up of nodes and edges. Nodes represent microbial groups included in each network. Edges represent statistically signi cant and robust associations between nodes, with the number of edges connected to a node referred to as the node's degree. The nodes that have the highest degree in the network are hub species and therefore are strongly interacting species [62]. The betweenness centrality of a node is calculated as the total amount of shortest paths from all nodes to all other nodes that pass through the node [63]. Thus, despite a node has low degree (e.g., 2), it can have the highest betweenness centrality in a network and can affect large sections of the network if it connects clusters/compartments that make up the network. The nodes with highest betweenness centrality in the network are also termed "gatekeeper" [64]. Therefore, they are thought to be crucial for ecological network structure and persistence because they literally hold the network together [65]. Besides degree and betweenness centrality, the relative abundance of each node (family) was also imported into Gephi. These three parameters (i.e., relative abundance, degree and betweenness centrality) can be visualized by the size of nodes in individual network.
The topological properties of microbial networks were calculated with indexes including component, average clustering coe cient, modularity, network diameter, graph density, average path length, the numbers of positive and negative correlations, and network fragmentation (f). A component can have several clusters which include densely interconnected nodes, and clusters can be interpreted as groups of taxa with overlapping niches [40]. Average clustering coe cient describes how well a node is connected with its neighbors on average in a network, and high values mean that the co-occurrence network include more highly correlated microbes and thus has 'small world' properties (nodes are more connected) [66]. Modularity quanti es to what extent networks can be broken up into smaller components. Modules may represent different niches in microbial networks and have been used to study habitat partitioning/preferences [67], and it could indicate ecological processes governing community structure, such as conserved inter-species communication [68]. Network diameter, graph density and average path length are used to brie y describe the density of microbial networks. For example, average path length is calculated as the average number of steps in the shortest paths between each node. Networks with a small value of path length are also known as "small world" networks [66], and it can be explained as an increase in the response of the microbial network's to disturbances [69].
Most of the parameters mentioned above were directly calculated in Gephi, except f was calculated as the ratio of the number of disconnected subgraphs (CL) to the overall number of nodes (N) in each network as log(CL)/log(N) [46]. The bootstrap 95% con dence intervals show consistent trends for f as the single empirical values [46]. The f ranges from 0 to 1, and closer to 1 represents more fragmented and less stable networks. Loss of "gatekeepers" contributes disproportionately to network fragmentation, suggesting high fragility of these networks upon selective removal of species [46,70]. To test and evaluate the stability of microbial networks, the f was recorded by iteratively removing the node with the most abundance, betweenness centrality, and degree for 10 times [46,71]. Within this process, properties of the network were recalculated after each removal and prepared for the next round operation.
Finally, because microbial groups contribute differently to the network structures, the distribution of nodenormalized degree (the number of connections a node has standardized by the total number of connections in the network) in each network was examined. Also the correlations between nodenormalized degree and betweenness centrality were explored. The network with stronger correlation between node-normalized degree and betweenness centrality could be more stable compared to those with low correlation. This is because when a network containing more nodes with both high betweenness centrality and high degrees, it is more likely to have alternative nodes to hold the network together when certain nodes with high betweenness centrality were lost [36]. Analysis of similarities (ANOSIM) test was used to determine if the differences between two or more networks were signi cant. The ANOSIM test was performed using "anosim" functions in the vegan R package (R version 3.6.1).

Multivariate regression tree (MRT) analysis
To determine how the CB microbiomes respond to spatiotemporal variations, multivariate regression tree analysis (MRT, [72]) was used to evaluate the hierarchical effect of environmental changes on the microbiomes. It was determined and generated by the R package "mvpart" (R version 3.6.1). Divisions in the MRT were determined by cross-validation.

Structural equation modeling (SEM)
Piecewise structural equation model (SEM) was used to analyze hypothetical pathways that may explain how spatiotemporal variations in the Bay shaped the composition, distribution and interactions of microbiomes [73]. It allowed us to partition direct and indirect effect of each environmental variable relative to other variables and to estimate the strength of multiple effects. The rst step was to build a priori model based on the known effects and relationships among environmental factors and dynamics of prokaryotic microbial communities. From MRT analysis results, we identi ed temperature, salinity and TSS as the primary driving factors. These three environmental variables covered seasonal variations and in uence from freshwater or ocean input, and therefore directly/indirectly affected microbial interactions through changing Chl α, nutrient availability (e.g. N, P), and PC. Then, a maximum likelihood goodness-oft test was used in model tting and a non-signi cant P value indicates a well-t model [74]. In comparison to traditional SEM, piecewise SEMs are less restricted by the number of links per sample size and Fisher's C can be used as the goodness-of-t indicator [57,75]. Piecewise SEM was performed in R (version 3.6.1) by using the piecewiseSEM package [57].

Microbial networks across all samples
The correlation-based co-occurrence networks including all samples were constructed and visualized by relative abundance, degree, and betweenness centrality, respectively (Fig. 1). The resulting CB microbiome networks consisted of 265 nodes (prokaryotic microbial taxa) and 2,525 edges (signi cant and robust associations between taxa; average degree 19.1) (Fig. 1). Topological properties commonly used in network analysis were calculated and listed in Table 1. The average network path length between all pairs of nodes was 3.37 edges with a diameter (longest distance) of 10 edges. In total, the network contained 8 components (those separated subgraphs in a network) and 17 clusters (a group of nodes with higher number of within-cluster edges than between-cluster edges). The clustering coe cient (the degree that nodes tend to cluster together) was 0.576. These clusters were a liated to almost all major prokaryotic microbial groups in the Bay, and also showed both positive and negative relations with many nodes in other clusters. Plentiful prokaryotic microbial groups clustered into one cluster with both positive and negative relations means that microbial associations were diverse and complex in the Bay. These clusters revealed the potential ecological relationships among prokaryotic microbiomes in the Bay and how they organized by niches in the estuarine ecosystem with habitat partitioning/sharing. Overall, the microbiome network in CB was comprised of highly interconnected taxa, formed a clustered topology, and had 'smallworld' properties (nodes are more connected).

Microbial interactions in different seasons
Distinct seasonal networks were obtained after applying the identical thresholds as described above (Fig. 2 and Table 1). For clarity and brevity, only networks with node sizes showing the number of microbial connections (degree) were included (Fig. 2). The seasonal networks were comprised of highly connected microbial taxa and densely connected groups formed a clustered topology with comparable seasonal variations (Fig. 2). The modularity indexes across four seasons were all greater than 0.4 suggesting modular structures especially in the autumn (Table 1) [76]. The number of component in the network of each season was 1, and the amount of clusters was almost the same in four seasons (6 in spring and 7 in other three seasons) ( Table 1). The summer network consisted of the biggest number of nodes (315) and edges (3399) with average degree 21.6 ( Table 1). Although the network in winter contained the lowest number of nodes (218), it had the highest average degree (31.1) with 3388 edges (Table 1, Fig. S2). In addition, the network in winter contained the highest graph density, shortest path length and high clustering coe cient (0.562) ( Table 1). All these characteristics indicated that winter microbiomes in CB were more closely associated and correlated.
The nodes with high degree/betweenness centrality were different for each season (Table S2). For example, one member of Archaea (Nitrososphaeraceae) and four members of Bacteroidetes (OPB56, Cyclobacteriaceae, Sphingobacteriales, Weeksellaceae) were included in the top 10 nodes with high degrees in the network of autumn, but no Archaea and Bacteroidetes groups were included in the top 10 nodes from other three seasons (Table S2). Deltaproteobacteria (SAR324 clade marine group B) and Chlamydiae (Chlamydiaceae) were ranked the top 10 nodes with highest degrees in summer, but these two groups didn't occur in the list of top 10 nodes for other three seasons (Table S2). These speci c and distinct microbial taxa played various roles in structuring the ecological networks for each season.

Microbial interactions in upper Bay vs. lower Bay
Similarly for clarity and brevity, only networks with node sizes indicating the number of microbial connections (degree) were shown here (Fig. 3). Networks for upper and lower Bay were strikingly different regarding size and topology ( Fig. 3 and Table 1 Table 1). The average path lengths were similar (3.0 vs 2.9), but the number of component and clusters in the network of upper Bay (7 and 13) were more than the network of lower Bay (4 and 10) and it means the network of upper Bay was more split compared to the lower Bay (Table 1). High clustering coe cient were also observed in both networks, especially in the lower Bay suggesting highly correlated microbes were included with 'small world' properties. Compared to upper Bay, more microbial taxa and more pronounced associations occurred in lower Bay with higher modularity and lower number of component and cluster ( Table 1).

Stability of estuarine microbiome networks
We analyzed the responsiveness of network fragmentation (f) to removal of top 10 nodes with highest relative abundance, degree and betweenness centrality, respectively. Compared to the original network (f = 0.37, number of components = 8), the f and the number of components changed more dramatically in the treatment of removing 10 nodes with highest betweenness centrality (f = 0.50, number of components = 16) compared to the treatment of removing 10 nodes with the highest relative abundance (f = 0.45, number of components = 12) or degree (f = 0.38, number of components = 8) ( Table 1). Similar results were also observed in modularity (Table 1), suggesting removal of the nodes with high betweenness centrality could much strongly undermine the stability and persistence of the prokaryotic microbial network than removal of the most abundant/highest degree nodes. Combined with the results of Fig. 1, our results clearly showed that the most abundant prokaryotic microbial groups are neither necessarily the hub species nor those groups with high betweenness centrality in CB. Compared to the abundant taxa, the minor/low abundance groups with high degree or betweenness centrality contribute greatly to the stability of prokaryotic microbial interactions in the Bay.
We further analyzed the stepwise responsiveness of co-occurrence network fragmentation to removal of the top 10 nodes with highest betweenness centrality. f started an increase from the second round in winter, but it maintained 0 even after 10 removals in autumn ( Fig. 4 and Table S4). The f values in winter were signi cantly higher compared to those of spring, summer and autumn (P < 0.01) (Fig. 4). Also the number of components increased from 1 to 4 in the network of winter (2 for the rest three seasons) after the removal process ( Fig. 4 and Table S4). These results further elaborate that the stability of microbiome networks varied across seasons in the CB, where microbial interactions were more stable and resilient in warm seasons (autumn) than those in cold season (especially winter). Similarly, microbial networks showed higher stability in lower Bay than upper Bay ( Fig. 4 and Table S4): the number of components (10), clusters (14), and f (0.4094) increased greatly in the upper Bay compared to the original network while the network properties in the lower Bay remained the same or changed very little after the removal of top 10 nodes (Table S4).
Seasonal networks contained more and high percentage of signi cant negative correlations (834-1170; 29.7%-36.6%) than spatial ones (511-564; 13.6%-14.0%) ( Table 1). The high proportion of negative correlations in each season could increase the stability and persistence of microbial network. Further, strong correlation between node-normalized degree and betweenness centrality were observed in each season (ρ = 0.55-0.78) compared to the spatial ones (ρ = 0.41-0.54) (Fig. 5), and this clear separation between groups were tested and con rmed by ANOSIM statistical analysis (P < 0.05). The stronger correlations between node-normalized degree and betweenness centrality meant that the nodes with high betweenness centrality ("gatekeepers") also comprised high degree (i.e. hub species). Therefore, once the "gatekeepers" were removed due to disturbances, likely they could be replaced by other hub species, which would connect different compartments and hold the network together. This can help to maintain the stability and enhance the anti-interference ability of the network. In addition, high modularity was observed across seasons compared to spatial scales (Table 1). Modules could use to visualize different niches in microbial networks and thus habitat partitioning/preferences were more strongly existed under seasonal variations compared to spatial ones in the CB with long residence time.
Overall, based on our results of fragmentations, proportion of negative correlations, the correlation between node-normalized degree and betweenness centrality, and modularity, this study showed that the structure, properties, and stability of prokaryotic microbial networks were distinct across spatiotemporal variations. Microbial interactions in warm seasons (especially autumn) were more stable than cold season (i.e. winter), and the lower Bay was more stable than upper Bay.

Associations with environmental factors
To further explore the effect of environmental variations on microbial networks, environmental factors were included in the network analysis (Fig. 6), and MRT and SEM analyses were also applied ( Fig. S3 and   Fig. 7). Signi cant correlations between environmental factors and microbial taxa were visualized by network and clusters (Fig. 6). The majority of taxa (nodes) in the same cluster may have close relationships or share similar ecological niches (i.e., co-occurrence). Seasonal changes (temperature), salinity gradients and TSS (and turbidity) levels were signi cantly correlated to nutrient availability (nitrogen, phosphorus), PC, and Chl a concentrations in the Bay, and all of these factors were responsible for separation and topology of microbial networks across season and space (Figs. 1, 2, 3 and 6). Total 163 taxa were identi ed closely associated with temperature while 121 strongly responded to salinity. Meanwhile, TSS, nutrient availability (mainly N and P), PC, and Chl a also played signi cant roles in microbial distribution and interactions (Fig. 6).
MRT results con rmed that temperature, salinity, and TSS were ranked the top drivers to shape spatiotemporal variations of microbiomes in CB (Fig. S3). In summary, the accumulative variance explained by temperature, TSS and salinity were 39.2%, 12.6%, and 4.9% respectively (Fig. S3). Based on these, relationships between microbiomes and major environmental variations were further assessed by piecewise SEM (Fig. 7). Due to the limited space, only key correspondences were highlighted in the Fig. 7, and the others and co-varied relations were included in the supplemental Table S5. We found that temperature, TSS and salinity contributed greatly to the distribution and associations of microbial groups and also strongly affected other environmental parameters including nutrient availability, PC and Chl a (Fig. 7). For instance, increase of temperature was directly associated with decreased Alphaproteobacteria, Betaproteobacteria and Bacteroidetes (negative effect), and increased abundance of Deltaproteobacteria, Planctomycetes, Gemmatimonadetes and Cyanobacteria (positive effect) (Fig. 7). Through correlations with nitrate and DOP, the temperature also affected other microbial groups such as Cyanobacteria, Actinobacteria and Chloro exi (Fig. 7). Salinity gradient in the Bay was signi cantly linked to Actinobacteria, Cyanobacteria, Verrucomicrobia, Deltaproteobacteria and Chlamydiae, while it indirectly in uenced Chloro exi, Verrucomicrobia, Firmicutes through reducing the availability of nitrate, ammonium, Chl a, and OP (Fig. 7). Directly or indirectly via PC, ammonium and TON, TSS impacted several groups including Deltaproteobacteria, Chlamydiae, Chloro exi, Firmicutes and Patescibacteria (Fig. 7). Generally, our results clearly showed how spatial and temporal variations of environmental factors affected estuarine microbiome composition, distribution, and interactions (networks) (Fig. S3, 6,   7).

Discussion
Rare groups contribute greatly to the stability of prokaryotic microbial networks The estuarine microbiome networks were comprised of highly interconnected taxa, formed clustered topologies, and thus contained 'small-world' properties [77]. Among the interconnected taxa, the abundant microorganisms contribute greatly to the earth's biogeochemical cycles and primary production [78], and are essential for organisms' survival [79]. However, the most abundant microbial groups are not those taxa that hold the network together, such as "hub species" or "gatekeepers" (Fig. 1). Thus, abundant taxa may not be necessarily critical to the prokaryotic microbial network structure or their stability (Table 1 and S1). In this study, the hub species and gatekeepers are relatively low in abundance or belong to rare species (Table S1), but they play fundamental roles in network persistence and contribute greatly to the stability and resilience of these taxa-taxa networks [80]. Recent studies have increasingly emphasized the ecological importance of the rare biosphere, because rare taxa can include more metabolically active microorganisms than abundant taxa (as measured by RNA to DNA ratios), and they may be keystone species in regulating the functioning of aquatic environments [81,82]. In addition, the rare microbes have been shown to ful ll essential functions associated with nutrient cycling, and may enhance functionality of the abundant microbes [83]. For example, rare species were demonstrated that it could offer the required gene pool to catalyze complex degradation processes of organic compounds, and some pollutants are often degraded by species falling below the detection limit in pristine samples [84].
Our results further corroborated the signi cance of "hub species" and "gatekeepers" via stability testing with removal processes. The responsiveness of network fragmentation to removal of nodes with highest betweenness centrality provide important insights into the susceptibility of prokaryotic microbial networks to disturbance [46]. Our results suggest that the loss of those potential "gatekeepers" contributes disproportionately to network fragmentation, which essentially agrees with earlier reports on food webs and mutualistic networks showing high fragility/susceptibility upon selective removal of taxa [33,70,85]. Sequencing data allowed us identify potential "gatekeepers" were a liated with uncultured Thioalkalivibrio sp (Gammaproteobacteria), Rubinisphaeraceae (Planctomycetes), OPB56 (Bacteroidetes), and Rhizobiales incertae sedis (Alphaproteobacteria). These taxa had highest betweenness centrality values (up to 743, 952, 1432 and 959, respectively) and were consistently present in the major component of the co-occurrence networks. The loss of "gatekeepers" may adversely affect the robustness and resistance of microbiome structure and associations [33,70,85]. Similarly, strongly interacting species (i.e., "hub species") are important to CB microbial communities, and they were shown to be able to steer microbiome ecosystems towards certain community types [86]. Additionally, as part of the microbial "seed bank", rare taxa (including these high degree/betweenness centrality but low abundance species) can potentially drive ecosystem responses to environmental changes and become dominant under favorable conditions [87], therefore providing a mechanism for community persistence and stability [88].

Seasonal and spatial variations of CB microbial networks
Similar to microbiome structures [14,54], microbial interactions and networks showed strong temporal and spatial patterns with distinct property and stability (Fig. 2-5 and Table 1). These differences are likely due to strong gradients in seasonality, salinity, nutrient availability, and other causal environmental factors [11,18,89,90]. CB, the largest estuary in northeast America, experience typical seasonal changes and constant freshwater/oceanic water input and exchange. The dynamic estuarine gradients lead to large variations in bulk bacterial production and biomass [19,20], and community composition [11,22], and subsequently affecting the property and stability of microbial networks (this study, Fig. 2-6 and Table 1; [52,91,92]).
Warm season microbial networks revealed low number of components and clusters, low fragmentation, high modularity, and high percentage of negative links. Stronger correlations between node-normalized degree and betweenness centrality and less variability of fragmentation after removing the potential 10 "gatekeepers" were also observed in the networks for summer and autumn ( Fig. 4 and Table 1). These results indicated that microbial networks were more stable during warm seasons compared to those in winter and spring. Microbial interactions affected by temporal variations were also found in Alaska Beaufort Sea coast [51], the San Pedro Ocean Time-series station [91], and the Lake Mendota [93]. The higher number of nodes in the co-occurrence networks of warm season agrees with the previously reported high microbial diversity in the Bay [14,54]. High biodiversity are able to promote interactions between microbial communities [94], and these biotic interactions including competition, are commonly thought to increase co-occurrence in microbial networks as they refer to common resources and environmental conditions [30,95,96]. Further, high percentage of negative links in the networks could stabilize co-oscillation within microbial communities and increase network stability [31]. In addition, increase in grazer richness has a positive effect on the bacterial richness and evenness [97], stemming primarily from the widespread distribution of resources and then resulting ecological niche complementarity [98][99][100]. Finally, this stability is also due to long residence time and relatively high growth rate (and mortality) of different microbial populations [20,90], which allows estuarine bacteria overwhelm those allochthonous populations from marine or freshwater [101]. Therefore, the results of more stable prokaryotic microbial networks in warm season compared to those in cold season in the Bay might arise from a better balance of the microbial associations (e.g., mutual, competitive and prey) and nutrients availability during summer and autumn than winter and spring [11,20,102]. In contrast, microbial interactions in cold seasons contained high fragmentation and it may be reinforced by stronger nutrient limitations and lower growth rates (due to low temperature) [20]. Interconnected microbial groups in winter were observed in the Bay, and it might be due to less predator pressure, lower temperature/light intensity, and less DOM release upon phytoplankton decay [20,25,103].
Compared to the upper Bay, the network in the lower Bay was more connected and contained less number of components and clusters with more negative correlations (Fig. 3, Table 1). A strong association between node-normalized degree and betweenness centrality and the less variability of fragmentation after removing potential "gatekeeper" nodes were also observed in the network of lower Bay (Fig. 4, 5).
Our results suggested that the stability of prokaryotic microbial network was higher in the lower Bay compared to the upper Bay. The stability difference between upper and lower Bay could be due to the different disturbance and interference from freshwater vs. oceanic water. Similar environmental variations were also found in many other estuaries, and spatial variation also could affect bacterial associations, including Ems estuary [92], Hangzhou Bay [43,53] and Pearl River Estuary [104]. River input, land runoff, suspended particles/turbidity, and accompanied nutrients availability are all pulsating with strong seasonal/inter-annual variations and uncertain patterns in the upper Bay. Conversely, salty water intrusion from the North Atlantic Ocean can be mild and consistent, and we thus hypothesize that less temporal disturbances from the ocean contribute to higher stability of microbial networks in lower Bay compared to the upper Bay.

Environmental dynamics and their in uence on estuarine microbiomes
The estuarine microbiomes and their distribution are a comprehensive output of the bio-interactions between microbial populations and the response to the surrounding environmental gradients. Our results show that shift of environmental factors (including temperature, salinity, TSS, nitrogen, phosphorus, and PC) has strongly direct and indirect effects on microbial community structure and networks, with the potential to affect community responses to future disturbances. Those conditions within the spatiotemporal variations seem to favor strong co-occurrence patterns (lower fragmentation and higher stability) of prokaryotic microbial associations, implying elevated biotic interactions or species sorting strongly mediated by the local environment [96,105]. Physiologic predisposition and nutritional tolerance of microbiomes tend to maintain stable communities inter-annually during certain seasons [11,25,106].
Temperature, salinity, and TSS are the main drivers to shape spatiotemporal variations of microbiomes in CB. Different bacterial groups in the Bay were directly/indirectly in uenced by these drivers and other factors (including Chl a, turbidity, nitrogen, phosphorus, and carbon). These associations lead to the property and stability (resilient) of prokaryotic microbial networks. Focusing on the identi ed prokaryotic microbial networks and their responses to environmental variations could provide us valuable insight into the microbiome adaptation and habitat partitioning/preferences in the Bay across seasonal change and spatial gradients [30,67,107]. Co-occurrence networks reveal critical information on co-oscillation between microbial taxa and also the stability of these involved communities in the Bay [30,36,107]. Therefore, changes in estuarine microbial networks resulting from disturbances provide a potential to examine the legacy effects on estuarine microbial functioning (e.g., contribute to the primary production, food webs and economical sustainability) and vulnerability to future disturbances, including anthropogenic in uence (e.g., eutrophication, contamination and damming) and climate change (e.g., drought and ood) [11,31,107].
To the best of our knowledge, this is the rst systematic and thorough study of the microbiome interactions, their stability, and the corresponding environmental factors in an estuary with long residence time, the Chesapeake Bay. In addition to the environmental gradients included in this study, many other abiotic or biotic interactions also play critical roles in microbiome composition and distribution. For example, grazing and viral lysis are important factors that may affect the interactions and stability of estuarine microorganisms [108][109][110], which deserve future investigations. Overall, the network stability measures used here suggest that natural scale spatiotemporal shifts of environmental variables did not move the systems out of their safe space in the Bay, which may provide a basis for optimism of traditionally high nutrient systems (such as the estuaries) [111,112]. The question remains whether future disturbances such as anthropogenic pressures and climate change will in uence the microbial regime and its resilience and adaptability, which may change the basic properties of the microbial networks due to temporal and spatial variations. Therefore, it is necessary to understand the effects of natural spatiotemporal variations on microbial interactions as well as their stability, which are critical for estuarine ecosystem processes.

Conclusions
This study shows that environmental variations dictate stable and resilient estuarine microbiome interactions in the Chesapeake Bay. Compared to dominant taxa, rare groups (e.g. hub species and gatekeepers) contribute greatly to the network topology and stability. Microbiomes in warm seasons exhibit stronger and more interactive networks than cold seasons, while the microbial network in the lower Bay is more stable than the upper Bay. The freshwater input from upstream brings pulsating and uncertain factors to microbial populations in CB, such as nutrients, particles, and microorganisms from terrestrial and freshwater environments. In contrast, the intrusion of seawater at the lower Bay is relatively mild and consistent so that the Bay microbiomes are able to adapt to the subsequent impact. Our piecewise structure equation model (SEM) unraveled that spatiotemporal variations (i.e. temperature, salinity, TSS, nutrient availability) are the primary driving factors for microbiome structures and interactions, and these environmental gradients manipulate the stability and tolerance/susceptibility of estuarine microbiota. This work shows, for the rst time, dynamics of microbiome networks and their responses to environmental gradients in a typical estuarine environment with long residence time. Given such dynamic environmental gradients, stable and persistent microbiome networks suggest that biotic interactions play a central role in maintaining integrity and resilience of estuarine microbiomes.