The assembly of the microbial community is commonly recognized as the result of deterministic and stochastic processes. The role of deterministic processes has been observed to be limited in stable environments, where their stochasticity could play an important role in gradually shifting community structure (Zhou and Ning, 2017). Due to the mutual influence of both processes, the relative abundance of a specific species is assumed to be the sum of a baseline and random variation. As the variation of species includes the random section, the random variation of species can be a joint distribution over time. The present microbial community has its evolution from the previous state, while time should have limited influence on the variation of the microbial community due to the presence of stochasticity. Although past observations are significant to forecast future trends, these predictions are not entirely dependent on them. Therefore, there could be an autocorrelation process, which might produce a time lag representing only finite past values for forecasting. Deng et al. (2016) used the time lag to construct the correlation network with time-series data and unravel microbial succession within a uranium bioremediation site (Deng et al., 2016). Additionally, David et al. (2014), when analyzing the effect of host lifestyle on human microbiota, relied on the autocorrelated process of time series (David et al., 2014), which demonstrated that OTUs’ variation complied with the time series model. We applied the data of 98 key OTUs obtained over the course of 259 days to fit into the augmented Dickey–Fuller (ADF) test to verify whether microbial data was irrelevant to time or not. If the data was not stationary, it indicated that the time series data were independent of time and the difference between the adjacent values will be applied to all data. The results of the stationary check are shown in Supplementary Data S2. All OTUs fulfilled the requirement of stationary after difference, 51 OTUs required difference treatment while the rest were stationary without it.

The visualized causal network is shown in Figure 1. A total of 98 OTUs were used for Microbial Granger Causal Network (MGCN) construction, which created 1,865 links between the nodes at a significant threshold of p < 0.05. Granger causality is commonly non-symmetric, and network building and the bidirectional links were defined as feedback from the source to the target OTU, indicating that either node could improve the forecasting accuracy of the other. A unidirectional link indicated that the source OTU significantly improved the forecasting accuracy of the target OTU but not vice versa. The outdegree and indegree directed links, defined by the direction of links in or out of the specific node, were counted separately. As revealed in Table 1, the distribution of the nodes degree tended to be normal rather than following power-law, regardless of whether indegree or outdegree, implying that the causal network was not scale-free (Deng et al., 2012).

The average clustering coefficient, which reflected the clustering degree of the overall network, was defined as the average of the clustering coefficients of all the nodes. The clustering coefficient of MGCN (0.449) was higher than previously described undirected networks, including grassland soils (0.1–0.22), lake sediment (0.09), and groundwater condition (0.17–0.29), and was comparable with the value of 0.466 observed in a previous activated sludge study (Ju and Zhang, 2015). Watts and Strogatz (1998) introduced the random rewiring procedure to interpolate regular and random networks, in which the regular lattice was highly clustered while the random network was poorly clustered (Watts and Strogatz, 1998). Therefore, the higher relative clustering exhibited by causality indicated that the network was defined rather than random. The average shortest average path was 2.149, which was smaller than that within the undirected network. Hence, we derived a relatively clustered network connected by shorter paths, demonstrating that the neighboring nodes were closely connected. To confirm the small-world property, randomized networks with the same nodes and degrees as the original network were constructed. The average clustering coefficient and shortest paths were ∼0.196 and 1.823, respectively, whereas the ratio of the Granger network to the random network of clustering coefficient and shortest path can be determined (Liao et al., 2011). The ratio was equal to ∼1.943, and this indicated that the network possessed small-world properties.

According to the definition of cs, its magnitude represents the ability of a specific OTU to cause variation among its neighbors. A value of 1 indicates that an OTU can affect its neighbors without being affected by them, while zero indicates the opposite. The c_recip reflected counts of reciprocating links, which exhibited feedback behavior of each OTU, thereby higher values indicated that an OTU is likely to interact with others. Therefore, as shown in Figure 2, as c_recip increased, the cs inclined to approach 0.5, displaying an equilibrium of indegree and outdegree links. All nodes displayed a c_recip value of less than 0.5, suggesting bidirectional links were not dominant in the relationship of all nodes. However, it was interesting that more interactions could be positively related to the trend of cs. The cs and c_recip were both relatively quantified as the proportion excluded the magnitude of degrees, and node size in Figure 2 is proportional to the degree of connection with the neighboring nodes. The average number of neighbors for a node was ∼27.08. The majority of nodes with a large number of neighbors had a higher c_recip and neutral position of cs. Nodes with lower c_recip and lower cs indicated that more links were indegree, with the reverse, higher c_recip and cs indicating more links were outdegree. Integration of the relative proportion and neighbor number, which was considered as an absolute quantity, was beneficial for inferring the central output nodes in the network, which should possess lower c_recip, higher cs, within a fairly large size of neighbors. The average of cs was ∼0.491, showing that the number of outdegree and indegree links were nearly identical. The average of c_recip was 0.24, implying mutual cause is not predominant due to the lower proportion in total links. Additionally, ncf, the difference between net outdegree and net indegree, of nodes ranged from −20 to 21, as shown in Supplementary Data S3. The average individual outdegree was 8.14 and the average net indegree was the same. Moreover, the number of OTUs with positive ncf was greater than that of negative ncf, indicating more than 50% of the relationships in the system displayed Granger causality in the activated sludge system.

Recently, numerous studies have used Bonferroni-correction to improve the threshold of correlation prediction (Navarrete et al., 2015; Quinn et al., 2016; Castellanos et al., 2020). The Bonferroni-corrected MGCN (BoMGCN) was produced from the original significant network (Figure 3). The corrected network was sparser in comparison with the causal network, containing only 81 nodes and 730 links and a lower clustering coefficient (0.373). The reduced network was highly conservative as the Bonferroni-correction excludes all potential type I error (false link was accepted) and displayed a slightly improved stability as revealed by the R square of power-law. The value of outdegree R square was 0.064, close to zero, yet the value for indegrees was 0.468, indicating a significant increase. Although the values were too trivial to be fitted into power law, they indicated that some nodes in the BoMGCN had a greater or lesser effect on other nodes. An improvement of the scale-free property was also observed, as well as an increase in the small world index represented by an increase in the ratio of σ (∼2.617). This was caused by a decrease in the clustering coefficient and an increase in the average shortest path, indicating that the BOGCN was more likely to fall under the rules of a small world. Within the random network derived from BOGCN there was a clear decrease in the clustering coefficient (0.119). Additionally, the properties of total nodes were slightly distinguished from the original MGCN network as a clear decline of cd value was noticed. The average cs increased from 0.491 to 0.504. Overall, the BoMGCN reduced the size of the network while retaining its basic properties. According to the classification of OTUs in BoMGCN, Proteobacteria was the predominant node and the hub species was Nitrospira, indicating the nitrogen-associated species had a broader social connection with other microbes.

The MGCN disclosed the casual effect within the microbial community; however, information about positive or negative correlations between nodes was missing. Therefore a Bonferroni-corrected microbial correlation network (BoMCN) based on Spearman’s correlation (shown in Supplementary Data S4) was applied to supplement the MGCN, constructing the MCCN. The multiple relationships between two OTUs could be revealed more explicitly according to this combination of causality and correlation. Previously, correlation analysis was generally used to discern the negative and positive relations within a microbial network, indicating the ecological interactions between members of the community (Faust and Raes, 2012). The Granger test was also recently applied in the research of a microbial network (Ai et al., 2019; Mainali et al., 2019), but the positive or negative effects could not be predicted with a single method. As shown in Figures 4A, a combination of correlation and Granger causality could be used to construct a new relationship, which indicated the directional connection among nodes, including the positive or negative effect they had on each other. As displayed in Figure 4B, the MCCN was composed of 73 nodes and 441 links. Although the causality was at a higher level compared with correlation, i.e., all nodes with causal links indicated a strong mutual interaction, the missing nodes and links could be ascribed to the Granger causality, which was not a real causal relationship, due to the limitations of the method. Technically, the Granger test has been widely used for predicting the causal effect. In the context of the current study, the Granger test was utilized to forecast the relations between OTUs, which may lead to a better understanding of microbial behaviors and relationships within a community. But additional efforts would be required to verify interactions between species.

The combination of correlation and Granger test allowed the observation of more specific interactions between the two species, therefore a MCCN network could be applied to predict ecological relationships for community analysis. As shown in Figure 4A, there are seven patterns of species interactions, including mutualism, synergism, commensalism, neutralism, predation (parasitism), amensalism, and competition (antagonism) (Pepper et al., 2015). According to the results of MCCN, both mutualism and synergism should be a bidirectional edge with a positive effect on both species, as each species would derive benefits from the other, such that it would be difficult to distinguish them apart. Commensalism can be reflected by a unidirectional link with a positive effect as species A can obtain a metabolite produced by species B. Although species B would be irrelevant to the growth of species A, as there is no feedback from A to B, the sequencing data of the two species would be positively correlated as more species B would secrete more metabolites for species A. Oppositely, a unidirectional connection with negative effect is classified as amensalism due to the general release of inhibitors from species A to species B. Here, the quantity of species A will be relevant to the production of inhibitors, such as antibiotics, which can reduce the number of species B, thereby the contrary growth of species A and B will lead to a negative correlation. Although the predation (parasitism) can be implied by the negative bidirectional edge, the sequencing data used in this study contained only information from the 16S rRNA gene of bacteria, with no information about protozoa or phages, resulting in the exclusion of predation (parasitism) from the MCCN of the microbial community (Deng et al., 2016). Finally, a negative bidirectional link could also indicate competition between species. The MCCN is a powerful tool to recognize multiple interactions of microbes by specifying the endogeneity of correlation, which has been widely used as a statistic proof of microbial interaction within a network (Weiss et al., 2016).

The nodes with amounts of links would be considered “hubs” in the MCCN. OTU56 was the hub species with the greatest number of indegrees (31) and the second highest number of outdegrees (16). It was classified as genus Nitrospira, a globally distributed group of nitrite oxidizers (NOB), which was also verified to be able to achieve complete nitrification from ammonia to nitrate in one step (van Kessel et al., 2015). As shown in Figure 5 and Supplementary Data S5, S6, OTU56 closely interacted with 24 OTUs from the phylum Proteobacteria, 8 OTUs from the phylum Bacteroidetes, and the remaining 6 OTUs interacted with 5 additional phyla. A total of 21 OTUs displayed negative interactions with Nitrospira, of which 14 were amensalism and 7 were competition relationships. Interestingly, all competition interactions originated from Proteobacteria to Nitrospira, showing that a number of Proteobacteria may depress the growth of Nitrospira. This could be ascribed to the fact that most bacteria related to the nitrogen cycle were Proteobacteria (Costa et al., 2006). Additionally, OTU56 unidirectionally interacted with OTUs from Bacteroidetes, for which there were only two types of interactions, commensalism and amensalism, with 3 and 5 links, respectively. According to a global diversity and biogeography study of over 300 wastewater treatment plants, only 28 out of 61,448 OTUs, accounting for 12.4% of the 16S rRNA gene sequences, were defined as core OTUs, and these mainly consisted of Proteobacteria, Bacteroidetes, and Nitrospira in activated sludge (Global Water Microbiome Consortium et al., 2019). Therefore, the results of MCCN in this study were consistent, as Proteobacteria and Bacteroidetes actively interacted with the core species of Nitrospira, a group that played the crucial role in nitrification of the activated sludge. At the genus level, the majority of species that interacted with OTU56 were unclassified, and of those that could be identified, Azospira, which possesses denitrification activity, exhibited a mutualistic relationship with Nitrospira, as well as with OTU176167 and OTU92689, which were most closely related to the genus Dechloromonas, members of which are capable of reducing nitrate or chloride. The above mutualistic relationships could be achieved in nitrogen cycling processes, with denitrification removing nitrate as a product inhibitor to Nitrospira, meanwhile, Nitrospira could supply nitrate as a substrate for denitrifies.

OTU180929, which had the most outdegree links (17), was classified as Sinobacteraceae at the family level. Members of this family are known to play a role in the degradation of aliphatic and aromatic hydrocarbon compounds and small organic acids (Gutierrez et al., 2013; Zhang et al., 2018). The number of net outdegrees and net indegree indicated the trending of nodes to cause a change of others or be affected by others. OTU180929, belonging to the genus Zoogloea, possessed 13 net outdegree and 13 net indegree links separately. Zoogloea has previously been demonstrated to be a bacterial genus important in the process of floc formation (Shao et al., 2009), and in this study it is represented by OTU180929 and OTU178488. As shown in Supplementary Data S7, the ncf of Zoogloea had the highest value within the sum of classified OTUs, indicating that Zoogloea could enhance the growth of most species, i.e., it could be the foundation for the formation of activated sludge. However, the ncf of unclassified OTUs was still higher, reaching 18. The culture-dependent methods build the basics of microbiology research, which investigate the role of specific species (mostly filamentous) in sludge flocculation and foaming (Nielsen et al., 2009). The unclassified nodes in MCCN notified that a massive microbial dark matter remained in activated sludge that was not yet cultured. The network approach has been used to elucidate and prioritize the microbial dark matter in the microbial community (Zamkovaya et al., 2021). Although activated sludge has been a widely employed strategy in wastewater treatment plants for over 100 years (Nielsen and McMahon, 2014), its microbiome still contained many mysterious and abundant unknown species that are only gradually being elucidated by recent progress in culture-dependent and independent technologies. Although our results showed the ecological interaction in the microbial community according to their variation with time, the Granger and correlation test only depend on the time-series data, which still needs experimental validation to verify the specific metabolic dependency between different species.

In conclusion, the coupling of correlation and causality was crucial to understand the ecological interactions within the microbial community. The MCCN disclosed a sophisticated causal network in activated sludge and identified the fundamental species, with the highest ncf value, as Zoogloea. The MCCN and phylogenetic analysis together indicated that the core-species of Nitrospira (OTU56) could have mutualistic interactions with denitrifiers in activated sludge. However, most species that interacted with OTU56 were still unclassified, implying a greater sequencing depth would be the key to improve the understanding of activated sludge.

The sequencing data were acquired from NCBI (accession number: PRJNA324303), which has been published previously (Jiang et al., 2018). The time-series data set included sequencing data for 259 days taken from a long-term operational wastewater treatment plant. The primers were F515 and R806, which covered mostly bacteria and archaea. The achieved fastq files were combined and processed online using a galaxy platform (Feng et al., 2017). OTUs were created with 97% cut-off through the Uparse clustering method. RDP classifier assigned one representative sequence from each OTU to bacteria or archaeal taxonomy according to the 16S rRNA Greengene Database. The final OTUs table was prepared for the subsequent process. The phylogenetic tree was created with Mega software with N-J method and the visualization was completed online¹ (Letunic and Bork, 2019).

Stationarity is an important concept for time series analysis and is a precondition to Granger causality. The properties of stationarity were defined by three main factors in terms of mean, variance, and covariance. The stationarity indicated that there was no change of trend in the time data, and it was known as a changeless process of the joint distribution within a specific displacement. The stationary implied that the expectation value of OTUs would fluctuate around the mean value of their neighborhood rather than depend on time. This allowed an estimation of the significant interval for the variation. Therefore, the stationarity analysis should be performed before analyzing time series data. It can be tested by detecting the presence or absence of unit root. The ADF-test was employed to verify if the time series data conformed to the stationary. If the original data was not subjected to stationarity, we used the difference, one minus another one to calculate the difference, to obtain the stationarity data. All relative abundance data of OTUs were filtered with the stationary test, while data that failed to go through the ADF test after two rounds of using difference would be summed in a separate file as nonstationary data. Although the relative abundances of OTUs may vary on a large scale, even seemingly without a mean value, the difference would be stationary in most situations. The operating reactors could be affected by many factors, which would shift the microbial community via stimulation of the metabolism of specific species.

The Granger causality test is a statistical hypothesis test that determined the role of one time series in forecasting another one (Granger, 1969). Herein, the Granger causality is limited within interpreting the interaction of two OTUs that were subjected to the autoregressive–moving-average (ARMA) model. To i-th OTU, the ARMA model is shown as the equation:

where v_i is the random variation (white-noise series), δ is a constant, z is the abundance of species i, θ represents the parameters, and q is the model order. We simplified the equation for OTUs to the following format. Thus, we assumed the model for i-th OTU is X and the model for j-th OTU is Y. Both equations are as follows:

where μ_i is the random variable. To know the interaction of X and Y, we assumed X and Y are interplays in their respective model predictions. The real data could be applied to the following equation:

where b₀ and c₀ are coefficients, and b_l, c_l, and d_l are the parameters. If b₀ and c₀ are not equal to 0 at the same time, this will be a model with instantaneous causality. In other words, v_ij and μ_ji would be the key to determining the Granger causality if the variation could be decreased when applying the b₀ ≠ 0, and representing the j-th OTU can contribute to the prediction of i-th OTU, otherwise, there was no improvement in predicting i-th OTU with j-th OTU information. Therefore, the Granger causality can be tested by the ANOVA analysis to obtain a p-value. This relation between X and Y was termed as Granger causality, which implied X or Y can cause each other. Herein, the causal effects were attributed to the property of edges in the network.

All OTUs were filtered with two specific conditions such that OTUs with more than 80% non-zero values would be preserved. The residuals should comply with at least one relative abundance of individual OTU that reached more than 0.01% in all samples. The total number of OTUs was 98. The ADF test was applied to verify the stationarity of time series data and provide a proper lag for the next modeling process. The difference was calculated once OTUs failed in the ADF test, and the results of the difference were used to track the ADF test another time. All OTUs were reserved by the twice difference treatment. The time series matrix successfully inspected by the ADF test was used for the Granger test in pair. Before the operation of the Granger test, the order was determined by VAR (R package) (Pfaff, 2008). Subsequently, the lag was transferred to the Granger test. The p-value threshold of the Granger test was restricted by the following two methods. Due to the massively paired results, the links confirmed by the significance values could still cause statistical type I error, hence, we introduced Bonferroni multiple-comparisons procedure and false discovery rate (FDR) to correct the threshold. Bonferroni multiple-comparisons procedure was determined by the following equation:

In the FDR test, all links that were selected by a 0.05 significant cut-off are reordered according to the magnitude of the p-value. FDR values were calculated by the following equation.

where, i is the rank of the p-value in k links, which is the total links preserved by the previous threshold. The critical FDR value is normally 0.05. FDR has a great power to detect genuine positive effects, while the Bonferroni adjustment is more conservative and considers all comparisons to be statistically independent. The final file was transferred to Cytoscape software for further visualization and analysis. All analysis processes were completed with R, and several shiny apps had been built for this study (Stationary check², Granger Causality network website³, Correlation network⁴, and MCCN⁵). The specific instruction for each app is provided in Supplementary Data S1. All raw R codes were deposited in the Github website.⁶

The several properties of the causal network were referenced from the literature of Seth (2005), and termed as causal score (cs), causal density (cd), net causal flow (ncf), and causal reciprocity (c_recip). Table 2 shows equations for all corresponding properties.

As the network had been directed, outdegree and indegree represented the direction of edges within two nodes. The causal score (cs) was determined by the ratio of outdegree to total degrees of a specific node, reflecting the OTU influenced other OTUs rather than being influences. The causal score is defined as cs = the number of outdegrees divided by the number of indegrees in unweighted graphs (graphs in which all links are equivalent). If cs > 1, the corresponding OTU has active output, otherwise it is being passively influenced. The causal density is also termed as causal efficiency of the network, which, to some extent, represents the connectivity of the network. The net causal flow is the difference between the outdegree and indegree of each node, indicating the contribution of the individual node would be either active or passive. Herein, the active state represents the species intentionally affects others, while the passive indicates it is affected by others. Although causal flow is like a causal score, the former is intended to be independent of the quantity of balanced efferent and afferent connections. The causal reciprocity is the fraction of links with a directly reciprocal edge. Overall, the causal score and flow are applied to evaluate the role of each node, while the rest describes the whole network. Additionally, the supplemented indexes, including connectivity, centrality, stress centrality, etc., were analyzed with the Cytoscape software tool (Feng et al., 2017).