Molecular ecological network analysis reveals the effects of probiotics and florfenicol on intestinal microbiota homeostasis: An example of sea cucumber

Animal gut harbors diverse microbes that play crucial roles in the nutrition uptake, metabolism, and the regulation of host immune responses. The intestinal microbiota homeostasis is critical for health but poorly understood. Probiotics Paracoccus marcusii DB11 and Bacillus cereus G19, and antibiotics florfenicol did not significantly impact species richness and the diversity of intestinal microbiota of sea cucumber, in comparison with those in the control group by high-throughput sequencing. Molecular ecological network analysis indicated that P. marcusii DB11 supplementation may lead to sub-module integration and the formation of a large, new sub-module, and enhance species-species interactions and connecter and module hub numbers. B. cereus G19 supplementation decreased sub-module numbers, and increased the number of species-species interactions and module hubs. Sea cucumber treated with florfenicol were shown to have only one connecter and the lowest number of operational taxonomic units (OTUs) and species-species interactions within the ecological network. These results suggested that P. marcusii DB11 or B. cereus G19 may promote intestinal microbiota homeostasis by improving modularity, enhancing species-species interactions and increasing the number of connecters and/or module hubs within the network. In contrast, the use of florfenicol can lead to homeostatic collapse through the deterioration of the ecological network.

it is important to elucidate the network structures, topological role of species, and the underlying mechanisms, which are essential for maintaining the homeostasis of intestinal microbiota.
Antibiotics are widely used as prophylactic agents and therapeutics for the prevention or treatment of bacterial diseases in the aquaculture, and their use has been associated with the emergence of antibiotic resistance in bacterial pathogens, alteration in aquaculture environment and animal gut microbiota, weakening of the immunity system responses, and the increase in food safety issues [14][15][16] . Recently, the concept of non-antibiotic aquaculture farming has become popular 17 , and probiotics, defined as live microorganisms, can be used as an alternative to the antibiotics and is highly concerned for its benefits on intestinal microbial community together with the improvement on growth and immune system 18,19 . Antibiotics are known to seriously disrupt the intestinal microbiota homeostasis, while probiotics can positively promote it, rather than affect the composition of the microbial community 20 . Currently, there is no systematic study addressing the effects of antibiotics or probiotics on intestinal microbiota homeostasis, and no standard for the evaluation of intestinal microbiota homeostasis has been developed.
Sea cucumber represent one of the most economically important holothurian species in China. Our previous studies demonstrate that Paracoccus marcusii DB11 and Bacillus cereus G19 exert beneficial effects on the growth and innate immunity of sea cucumber [21][22][23] , while florfenicol had a negative effect on the intestinal epithelial cells and innate immunity 24 . However, the effects of these probiotics and an antibiotic on ecological networks within intestinal microbiota have not been reported previously.
To the best of our knowledge, this is the first study to report the effects of probiotics P. marcusii DB11 and B. cereus G19, and an antibiotics, florfenicol, on the intestinal microbiota homeostasis in aquatic animals, by assessing modularity, species-species interactions, and their topological roles. Our findings provide new insights into the effects of probiotics and antibiotics on the intestinal microbiota homeostasis through the modulation of ecological networks.

Results
Sequences obtained. In this study, a total of 2,720,976 high-quality sequences were generated by sequencing the V3-V4 region of the bacterial 16 S rDNA from intestinal content samples collected from the sea cucumber (median = 138,485 sequences, ranging from 113,552 to 153,934 sequences) with dietary basal diet (Control) and supplementation with probiotics P. marcusii BD11 (PM) and B. cereus G19 (G19), and antibiotics florfenicol (FL), respectively.
Richness and diversity. At a threshold of 97% sequence identity, a total of 42,147 OTUs were identified in the current study (median = 4489 OTUs, ranging from 3080 to 6086 OTUs). As shown in Table 1, the four groups (Control, PM, G19 and FL) had Good's estimated sample coverage (ESC) of 97.8, 97.8, 98.0, and 97.9%, respectively, indicating that most of the microbial diversity had already been captured with the current sequencing depth. To assess the species richness and diversity of intestinal microbiota of sea cucumbers, the Chao1 and abundance-based coverage estimator (AEC) and Shannon diversity were calculated by estimating the number of OTUs. Species richness and diversity were not significantly different between these groups in this study, while the lowest Shannon diversity presented in the FL group.
In order to test whether any difference was present in organismal structure of intestinal microbiota, Principal Coordinates Analysis (PCoA) was performed based on the weighted and unweighted UniFrac distances for the Modularity analysis. Four commonly used complementary network indexes can be used to describe network difference 13 : (i) connectivity, which is the most commonly used concept for describing the topological property of a node in a network; (ii) path length, which is the shortest path between two nodes; (iii) the clustering coefficient, which describes how well a node is connected with its neighbors; and (iv) modularity, which measures the degree to which the network was organized into clearly delimited modules. As show in Table S2, the highest average connectivity was observed in FM group, which means that the FM group has the most complex network. Significant differences between these four ecological networks and their corresponding random networks with identical network sizes and average numbers of links were observed in terms of the average path distance (GD), average clustering coefficient (avgCC), and modularity (P < 0.001; see Supplementary Table S2), indicating that these four ecological networks obtained possessed typical small-world characteristics. In addition, the GD, avgCC, and modularity in FM, G19, and FL groups were significantly different from that in Control (P < 0.001; see Supplementary Table S2), hence, the ecological networks in three additives groups were remarkably different from Control group. As shown in Fig. 2, circos plot represented the interaction between species of the intestine microbial community of sea cucumber. The network in four groups consisted of different OTUs from 30 bacterial classes, and the dominant classes were Flavobacteriia, Gammaproteobacteria and Alphaproteobacteria. The predominant class observed in PM and Control networks was Flavobacteriia, the relative abundance of which was more than that in the G19 and FL groups. The largest number of OTUs in G19 and FL group was Alphaproteobacteria ( Table 2). The blue and red edges respectively indicated the positive and negative interactions between two OTUs inside the circle. In the ecological network, one module is a group of OTUs that are highly connected among themselves, but had much fewer connections with OTUs outside the group. Random matrix theory-based approach is employed to delineate separate modules within the network. Dietary supplementation with P. marcusii BD11, B. cereus G19, and florfenicol distinctly affected the interactions between the members of the microbial community. As shown in Fig. 3 and Table 2, in the Control group, the ecological network consisted of 56 modules with 563 nodes (OTUs) and 1522 edges; a total of 22 of 56 modules with ≥5 nodes were obtained from the networks, and C1 and C2 were two biggest modules. Interestingly, in the PM group, the largest modules and the most complex interactions presented in this network; only 6 modules presented in the network with the largest number of nodes and edges, 682 and 5517, respectively; whereas 5 modules had ≥5 nodes, of which three largest module P1, P2, P4 were also observed in this network. In G19 group, there were 665 nodes and 1582 edges in the ecological network with 18 of 40 modules with ≥5 nodes, and G1, G6 and G8 were three biggest modules. The ecological network in FL group had 15 of 30 modules with ≥5 nodes, and the least number of nodes and edges presented in this network, 462 and 1504, respectively. Moreover, the dominant interactions in four networks were positive interaction. Strikingly, as shown in Fig. 3, many OTUs from the same class were clustered within one module.
Dietary supplementation with P. marcusii BD11 and B. cereus G19 increased the number of nodes (network size) and edges within the ecological network, while opposite results was found in the FL group compared to the Control group. Moreover, dietary supplementation reduced the number of sub-modules within the network. Obviously, the sub-module in PM group was extremely different from those in the other three groups. Although the number of sub-modules was remarkably less than that in the Control group, the PM network became more complex. Enormous species-species interactions was observed within/-out three huge sub-modules such as P1, P2, and P4, indicating tighter interactions/coupling within microbial communities. These results suggested that dietary supplementation with P. marcusii BD11 improved the stability of the intestinal community ecosystem.
Topological roles analysis. Species take different topological roles in the ecological networks. As shown in Fig. 4, the majority of OTUs that were observed in the Control, PM, G19 and FL groups were peripherals. As shown in Table 3, in the Control network, three OTUs from Flavobacteriia (OTU6390 and 20796) and Bacilli (OTU10436) played as connectors, and an OTU of Actinobacteria (OTU5575) served as module hub. In the PM network, seven OTUs from Flavobacteriia (OTU24782, 40175 and 19664), Gammaproteobacteria (OTU13249), Verrucomicrobiae (OTU40248), Bacilli (OTU28368) and Unclassifiled (OTU4847) played as connectors, and five OTUs from Gammaproteobacteria (OTU31556) and Flavobacteriia (OTU3524, 19153, 3449 and 14645) served as module hubs, respectively. In the G19 network, only one OTU of Alphaproteobacteria (OTU35727) played as connector, and ten OTUs from Alphaproteobacteria (OTU41885, 26519, 4049, 20466, 24292 and 24787), Gammaproteobacteria (OTU12475, 35688, 35347 and 2367) and Anaerolineae (OTU24787) served as module hubs, respectively. Interestingly, only one OTU of Alphaproteobacteria (OTU30901) served as connector in the FL network. No network hubs were found in these four networks. Module membership provides the best summary of variation in relative abundance of OTUs within a module. If module membership is close to 1 or −1, it is evident that the OTU is close to the centroid of module 12 .
The growth, nutrient digestion, and mid-intestinal morphology of sea cucumber. As shown in Table 4, the daily supplementation with P. marcusii BD11 and B. cereus G19 significantly increased the final weight and special growth rate (SGR) of sea cumber (P < 0.05). Additionally, P. marcusii BD11 remarkably enhanced apparent digestibility coefficient (ADC) of crude protein in sea cucumber (P < 0.05; see Fig. 5), and B. cereus G19 significantly improved the fold and microvillus height of mid-intestine in sea cucumber compare to   Figure S2).

Discussion
Sea cucumber represent good candidates for the study of the evolution and functions of intestinal microbiota due to their unique digestive system with only one simple intestine in the body cavity 25 . Gut microbial community plays an important role in the host health, and the composition of the core gut microbiota is considered to be essentially stable throughout adulthood 26 . In this study, the administration of P. marcusii DB11 (PM), B. cereus G19 (G19), and florfenicol (FL) was shown to affect intestinal microbial community with an obvious decrease in the percentage of Flavobacteriia (classified as Flavobacteriaceae in this study) in sea cucumber intestines, while core microbiota remained the same as in the Control group. Most importantly, no significant difference in alpha      diversity index was observed between the Control, PM, G19, and FL groups. Recently, Falcinelli et al. showed that the supplementation with probiotics had no effects on zebrafish gut microbiota composition in terms of alpha diversity 5 , and similar results were obtained in humans as well 4, 6, 27, 28 . However, Ferrario et al. showed that the probiotics supplementation can significantly modified the structure of fecal microbial community in humans, in terms of compositional dissimilarity 29 . In a recent study, Sanders postulated that probiotics may promote the homeostasis of intestinal microbiota, rather than affect its composition 20 . By contrast, the use of antibiotics may lead to dysbiosis due to their negative impact on the commensal microbiota 30 . Numerous studies showed that the antibiotics can lead to a reduction in the bacterial diversity and an increase in the abundance of antibiotic-resistant specific strains and species [31][32][33] . Notably, we showed that the use of florfenicol has a negative effect on the Shannon diversity index and leads to an increase in the relative abundance of the family Vibrionaceae belonging to Gammaproteobacteria class. Additionally, many pathogenic bacteria and opportunistic pathogens that were found in aquaculture environment belong to Vibrionaceae 34,35 .
In this study, we explored how the addition of B. cereus G19, P. marcusii DB11, and florfenicol affect species-species interactions in microbial communities by the Random Matrix Theory (RMT)-based network approach. The RMT-based network approach is a reliable, sensitive and robust tool for analyzing high-throughput genomics data for modular network identification and network interactions elucidation in microbial communities 12 . To the best of our knowledge, this is the first study to evaluate the changes in network interactions among different phylogenetic groups/populations of intestinal bacterial communities in aquatic animal in response to dietary supplementation with probiotics or antibiotics. Dietary supplementation with B. cereus G19, P. marcusii DB11, and florfenicol significantly affected the ecological network within the intestinal microbiota observed in terms of the average path distance (GD), average clustering coefficient (avgCC), and modularity. Modularity is the degree to which a network is divided into distinct sub-groups, and ecological networks can be naturally divided into different sub-modules considered as functional units, which respectively perform identifiable tasks in the ecological networks 12,36 . As shown here, each treatment had its unique ecological network model with characteristic modules, but the composition of ecological networks was shown to be similar to intestinal microflora, which indicates that the dominant microflora plays an important role in the ecological network. Here, only five larger sub-modules were observed in PM ecological network, which considerably differs from the results obtained in three other groups. However, this does not imply that only five tasks are performed by these sub-modules. Emergent property in a biological system means that a property of a network cannot be elucidated from the individual components, but it emerges as a consequence of the structure and interactions in the whole network 37 . Therefore, the results we obtained indicate that dietary supplementation with P. marcusii DB11 may lead to the integration of several sub-modules into a formation of a new large sub-module, which performs more functions than original individual sub-modules. Additionally, many OTUs from the same class are clustered within one module, and OTUs from the same species within a module most likely share the same functions 38 . Hence, the more OTUs from the same class belong to the same module, the more stable that module would be, because the loss of some OTUs would not disturb the overall function of the module. Accordingly, the ecological network in PM group is the most stable one, since only five modules but so many OTUs from the same class.
A network connection between two OTUs describes the co-occurrence of these two OTUs, which may be caused by the species performing similar or complementary functions 13 . The finding of this study highlights average connectivity in the PM network, indicating the higher-level species-species interactions within ecological network. Interactions that confer significant advantages to at least one of the populations can potentially result in the generation of a stable community. In the context of our models, the daily consumption of P. marcusii DB11 is expected to promote the stability of communities by providing an alternative energy source to microbes involved in microbial cross-feeding, as the number of positive interactions is higher than the number of negative interactions in four ecological networks. Positive interactions signify complementation or cooperation, while negative interactions may indicate competition or predation between the taxa. Cooperation was found to be dominant interaction in symbiotic communities, such as the microbial community in the intestine, where the microbes can be manipulated into a higher degree of cooperation 39 . Several models also suggest that the cooperation can be stable and that positive interactions are more likely to persist over time, as they keep the populations above the extinction threshold 40,41 . The results of our study suggest that the cooperative interactions are more likely to be stable as well. However, ecological competition is thought to be prevalent in natural microbial communities 42 . A recent study found that highly divers intestinal species are likely to coexist stably when the system is dominated by competitive, rather than cooperative, interactions 9 . These conflicting results may be related to the differences in the development of different models. However, we were not able to completely predict whether competitive or cooperative interactions are more likely to promote stability of intestinal microbial community.
Topologically, different OTUs play distinct roles in the ecological network 43 . The analysis of modular topological roles was an important step in the identification of key populations based on the OTUs' roles in their own modules. From the ecological viewpoint, peripherals may represent specialists whereas connectors and module hubs may be related to generalists and network hubs as super-generalists 7 . Structurally, peripherals can be lost without affecting the functions of ecological networks, while the loss of connectors and module hubs would lead to the deterioration of the entire network 44 . The daily consumption of P. marcusii DB11 and B. cereus G19 considerably increased the number of generalists within the ecological networks, and it made them more stable, which suggested that P. marcusii DB11 and B. cereus G19 may promote the homeostasis of intestinal microbiota in sea cucumber. Furthermore, we found that most of the generalists were from the same phyla and belonged to dominant genera in bacterial community. In an ecological network context, certain species act as structural and functional keystone species, and play an important overall role in maintaining the properties of their network 7 . Therefore, our results suggest that the dominant genera in intestinal microbial community perform important roles in the ecological network.
In contrast, the use of antibiotics is considered the strongest and most common cause of disturbance of the intestinal microbiota 30 , and here, we showed that the dietary supplementation with florfenicol considerably decreased the proportion of Flavobacteriia and caused the extinction of connectors and module hubs in the ecological network. Extinction of key species, such as generalists, may lead to the fragmentation of an entire module 7 , and the use of florfenicol disrupted sub-modules, leading to the deterioration of the entire network, and disturbing the intestinal microbiota homeostasis.
It is well known that the members of intestinal microbial community partake in numerous important physiological, nutritional, immunologic, and metabolic processes, supporting the idea that intestinal microbiota represent an external organ 1,45 . Accordingly, maintaining the homeostasis of intestinal microbiota can be beneficial for the host health, while the disturbance in the intestinal microbiota homeostasis has negative effects. Additionally, it is commonly assumed that the functioning of intestinal microbiota depends on the species that engage in cooperative metabolism and are beneficial for the host 46,47 . In this study, in PM group, the unique ecological network structure and complex species-species interactions were shown to enhance the functioning of intestinal microbiota, which contributed to the promotion of the apparent digestibility coefficient of crude protein. In the G19 group, the results of the micromorphology analysis showed an increase in microvilli and fold heights of mid-intestine, suggesting that B. cereus G19 may affect the expanding of the intestinal structures in the sea cucumber. Fold height, enterocyte, and microvilli are directly correlated with the functioning of the intestines and host health, and an increase in their heights leads to an increase in the absorptive surface area. These findings are in agreement with previous studies investigating probiotics 24,48,49 , and likely contribute to the improvement in the sea cucumber growth observed in this study. Furthermore, we have previously showed that the daily consumption of B. cereus G19 and P. marcusii DB11 significantly enhances the immune response in sea cucumber 22 . In contrast, the administration of florfenicol led to a serious atrophy of microvilli. Furthermore, the results of our previous study demonstrated that florfenicol induces the apoptosis of intestinal epithelial cells 24 , which leads to a decrease in nutrients absorption and an increase in the risk of infection by pathogenic bacteria. However, it should be further investigated whether the negative effects of florfenicol on intestinal structure are direct or indirect.

Conclusion
The analysis of the ecological network structure provides new insights into the intestinal microbiota homeostasis. Our results showed that intestinal microbiota homeostasis can be improved by the daily consumption of B. cereus G19 and P. marcusii DB11, which affect intestinal microbiota homeostasis through modulation of ecological network, by improving modularity, enhancing species-species interactions, and increasing the number of generalists, rather than fundamentally changing its composition. However, the use of antibiotics florfenicol can disturb intestinal microbiota homeostasis through the deterioration of ecological network, by reducing the number of generalists. Our work indicates that the analysis of ecological networks may represent an effective way to evaluate the intestinal microbiota homeostasis systematically. Further studies should provide more evidence to support this hypothesis.

Methods
Bacterial strains and antibiotic. Probiotics P. marcusii DB11 and B. cereus G19 are previously isolated from the intestines of sea cucumbers. Florfenicol (purity 99.0%) was supplied by from Shandong Lukang Animal Pharmaceutical Co., Ltd (Jining, China).
Experimental animals and diets. Disease-free sea cucumbers were from Laboratory Animal Centre, Ocean university of China, and acclimated to the experimental conditions (temperature, 17 ± 1 °C; salinity, 28-30‰; pH, 8.0 ± 0.3; dissolved oxygen, 10 ± 0.25 mg L −1 ) for 15 d prior to testing. Following a 24 h fast, similar size individuals (4.68 ± 0.07 g) were randomly distributed into 20 aquaria (53 × 28 × 34 cm, 50 L) at a density of 10 sea cucumbers in each aquaria. There are 4 groups with 5 biological replicates in this experiment, and each group has total 50 sea cucumbers. The basal diet (Control group) was formulated with marine mud, red fish meal, and sargasso. It contains 16.1% crude protein and 0.87% crude lipid (see Supplementary Table S4). On basis of the basal diet, the probiotics diet were supplemented with 10 9 cfu kg −1 two potential probiotics, i.e., P. marcusii DB11 (PM group) and B. cereus G19 (G19 group), respectively; the antibiotic diet was supplemented with 15.0 mg kg −1 florfenicol (FL group). With the exception of the FL diet, the composition of diets did not change throughout the 60-day feeding trial. In the FL treatment, sea cucumbers fed a diet containing florfenicol for 5 d, and then fed with basal diet without florfenicol for 15 d (three 20-d feeding cycles). The withdrawal period for florfenicol should not be less than 10 d 50 . All of the sea cucumbers in each group were weighted in the end of the experiment. Sample collection. During the last 2 weeks of the trial, the shaped feces for apparent digestibility coefficient (ADC) of crude protein were collected from each tank by pipetting every day at 08:00-10:00 am and 05:00-6:00 pm. After collection, feces were centrifuged (3000 g at 4 °C for 20 min) and frozen daily at −20 °C. At the end of the experiment, the intestinal content in hindgut from four sea cucumbers of each replicate was collected and mixed. Samples were frozen at −80 °C until further analysis. The mid-intestinal tract for tissue slice were injected with Bouin's fixative solution and transferred into 70% ethanol after 24 h later.
DNA extraction and 16S rDNA gene sequencing. PowerFecal ™ DNA Isolation Kit (MoBio Laboratories, Inc) was used for DNA extraction from the intestinal content samples. Amplification and sequencing of the V3-V4 region of the bacterial 16S rDNA gene was performed using barcoded fusion primers 341F (CCTACGGGNGGCWGCAG) and 805R (GACTACHVGGGTATCTAATCC). PCR amplification was then performed under the following conditions: initial denaturation at 98 °C for 30 s, 25 cycles at 98 °C for 10 s, 53 °C for 30 s and 72 °C for 30 s, and final extension at 72 °C for 7 min. The amplicons were pooled in equimolar concentration and sequenced with an Illumina MiSeq platform.
Bioinformatic analyses. The raw sequences were sorted into different samples according to the barcodes by using the BIPES pipeline, followed by chimera sequences filtering with UCHIME. After preprocessing, operational taxonomic units (OTUs) were picked at 97% similarity level against green gene version 13.8 using QIIME. Taxonomies were assigned with uclust for each OTU. The rarefaction curves were generated from the remaining number of OTUs. Alphadiversity (number of OTUs; Chao1 estimator of richness; abundance-based coverage estimator; Shannon diversity indices) and betadiversity (principal coordinates analysis (PCoA)) analyses were also performed using QIIME.

Molecular ecological network construction and visualization.
Based on the abundance profiles of individual OTUs, four phylogenetic molecular ecological networks were constructed with a random matrix theory (RMT)-based approach as describe previously 12,13 . As previously described, RMT-based approach were used for network construction, topological roles identification, module membership with an automatic threshold. To characterize the modularity property, each network was separated into modules by the fast greedy modularity optimization. Since only a single data point of each overall network index was available for each network parameter, standard statistical analysis could not be performed to assess their statistical significance. Thus, random networks were generated using the Maslov-Sneppen procedure 13,51 . Based on Z-test, the average path distance (GD), average clustering coefficient (avgCC) and modularity of the ecological networks were used as values to test the significance of the difference from random networks. According to values of within-module connectivity (Zi) and among module connectivity (Pi), the topological roles of different nodes can be categorized into four types: peripherals (Zi ≤ 2.5, Pi ≤ 0.62), connectors (Zi ≤ 2.5, Pi > 0.62), module hubs (Zi > 2.5, Pi ≤ 0.62) and network hubs (Zi > 2.5, Pi > 0.62). The construction and major analyses of molecular ecological networks were performed online (http://ieg.ou.edu/). Ecological networks were visualized using Circos 52 and Cytoscape 3.0.0 13 . Measurement of growth, nutrient digestion, and mid-intestinal morphology. The growth was calculated by the formula: Specific growth rate (SGR) = (Ln W t − Ln W 0 ) × 100/t; where W t and W 0 were final and initial sea cucumber weight respectively; t was duration of experimental days. ADC of crude protein was determined as described by Yang et al. 24 . The mid-intestinal tract samples for tissue slice were processed and analyzed by assessing the dimensions of intestinal folds, enterocytes, and microvilli as described by Peng et al. 53 .