*2.7. Statistical Analysis*

The statistical analysis was performed using SAS v9.4 for windows (Cary, NC, USA). The normality of distribution of data was checked using the Kolmogorov–Smirnov test (*p* > 0.05). The Levene test (*p* > 0.05) was applied for assessing the homogeneity of variance. One-way analysis of variance (ANOVA) was used to compare the average concentration of all the water quality parameters, concentrations of NR, TAN, NO2 −-N, NO3 −-N and TN on each sampling day, microbial alpha diversity indexes, shrimp growth indexes, the composition of microbial phylum and genus among treatments. Differences were considered significant at *p* < 0.05. When significant differences were found, Tukey's test was used to identify differences between treatments. The indices of Chao1, Abundancebased Coverage Estimator (ACE), Shannon and Simpson were calculated using Qiime 2.

The linear discriminant analysis (LDA) effect size (LEfSe) method (http://huttenhower. sph.harvard.edu/lefse/, 29 April 2021) was used for determining the features of bacterial communities between different treatments [30]. The LDA score was set as 3.5 to identify the significant features. The redundancy analysis (RDA) was carried on using the vegan package (v. 2.3) in R (v. 4.0). The functional prediction of bacterial community was performed by FAPROTAX [31] and BugBase [32].

#### **3. Results**

#### *3.1. Water Quality*

The descriptive statistics for the water quality data are given in Table 1. No significant differences in average temperature or DO were observed among the treatments (Table 1, *p* > 0.05). The pH in the CW system was significantly higher than that in the other treatments and pH in the P-BF system was significantly higher than that in the M-BF system (Table 1, *p* < 0.05). The average TA concentration in the CW system was significantly higher than that in the other systems (*p* < 0.05) and the biofloc volume in the M-BF systems was

significantly higher than that in the CW and P-BF systems. The Chlorophyll a in the P-BF system was significantly lower than that in the other two systems.

**Table 1.** The mean and standard deviations of water quality parameters in the three systems during the experiment.


CW, M-BF and P-BF represented the three systems of this study. TA: total alkalinity; TOC: total organic carbon; Chl: chlorophyll a; TAN: total ammonia nitrogen; TN: total nitrogen. In each row, different superscript letters indicate significant differences at the *p* < 0.05 level (one-way ANOVA and Tukey's test).

The average TAN concentration was not significantly different among the three systems (Table 1, *p* > 0.05). The concentration of TAN in the BFT systems (M-BF and P-BF) reached peak levels on day 7 of the experiment and then decreased sharply from day 7 to 14 (Figure 1a). Then, the concentration of TAN in the BFT systems remained relatively stable. The concentration of TAN in the CW system declined from days 1 to 14 and then increased slowly over that in the other systems (Figure 1a). No significant difference in TAN was observed among the treatments at the end of the study (Figure 1a, *p* > 0.05).

The average concentration of NO2 −-N in the CW system was significantly lower than that in the BFT systems (*p* < 0.05) and it was relatively stable. The NO2 −-N concentration in the M-BF system increased sharply from days 0 to 19, then, the concentration of NO2 −-N decreased from days 19 to 24 and was stable after day 24 (Figure 1b).

The concentration of NO3 −-N in all systems was stable from days 1 to 19 but increased sharply in the BFT systems from days 19 to 29 (Figure 1c). After day 29, the NO3 −-N concentration in the BFT systems decreased slightly. The NO3 −-N concentration in the CW system remained relatively constant during the study. The concentration of NO3 −-N in the CW system was not significantly different from the other two systems during the first 19 days (*p* > 0.05). After day 19, the concentration of NO3 −-N in the CW system was significantly lower than that in the carbon added systems (*p* < 0.05).

The average concentration of TN in the M-BF was 13.65 mg L−1, which was significantly higher than that in the CW and P-BF systems (5.69 and 9.75 mg L−1, respectively) (Table 1, *p* < 0.05). The TN concentration increased during the first 7 days in all systems but began to fluctuate after day 7 (Figure 1d). The concentration of TN in the M-BF system was significantly higher than that in the CW system after day 7.

The average nitrification rate in the CW system was significantly higher than that in the BFT systems (Table 1, *p* < 0.05). No significant differences were observed in nitrification rate among the three treatments on the first sampling date (d 14) (Figure 2). Then, the nitrification rate was consistently higher in the CW system than in the other systems (Figure 2). The nitrification rate in the BFT systems decreased dramatically after the first sampling and approached zero after day 19 (Figure 2).

**Figure 1.** Fluctuations of (**a**) total ammonia nitrogen (TAN), (**b**) NO2 −-N, (**c**) NO3 −-N and (**d**) total nitrogen (TN) during the period of the experiment. CW, M-BF and P-BF represented the three systems of this study. Values are mean ± standard deviation of three replications on each sampling day.

#### *3.2. Characterization of the Bacterial Community*

A total of 446,015 high-quality reads were obtained and 3,219 operational taxonomic units (OTUs) were obtained in the water samples by clustering the OTUs at a 97% similarity level. The Chao1, ACE, Shannon, and Simpson indices, which were used to represent the richness and diversity of the bacterial community, were not significantly different among the treatments. The phyla Bacteroidetes, Proteobacteria and Actinobacteria were the dominant bacteria (48.76%, 38.20%, and 3.27%, respectively) in the three treatments (Figure 3a). The abundance of Cyanobacteria decreased in the M-BF system and increased in the P-BF system compared to the CW system, which was 0.44%, 7.32% and 2.83% in the M-BF, P-BF, and CW systems, respectively. Differences in the phyla among treatments were compared and the top 20 phyla with the smallest *p*-values were listed in Figure 3c. Cyanobacteria were included in the top 20 phyla. Planctomycetes was one of the dominant bacteria in the M-BF system and the abundance of this phylum was significantly higher than that in the other treatments (*p* < 0.01, Figure 3c). The abundance of Chloroflexi was higher in the CW systems (5.58%) than in the BFT systems (0.69% in M-BF and 0.46% in P-BF) and was included in the top 20 phyla with the smallest *p*-value (Figure 3c).

**Figure 2.** Nitrification rate during the period of the experiment. CW, M-BF and P-BF represented the three systems of this study. Values are mean ± standard deviations of three replications on each sampling day.

As for genera of the bacterial communities, an uncultured\_bacterium\_f\_ Flavobacteriaceae, uncultured\_bacterium\_f\_Rhodobacteraceae, and *Leucothrix* were the more abundant genera in the M-BF treatment compared to the other two treatments (Figure 3b). *Tamlana* and an uncultured \_bacterium\_o\_Ardenticatenales were more abundant in the CW than the other two systems. Five genera were significantly different among the treatments in the top 20 genera with the smallest *p*-values (Figure 3d). *Vibrio* was significantly more abundant in the CW than the BFT systems (Figure 3d).

Linear discriminant analysis effect size (LEfSe) was used to identify the most differentially abundant microbial taxa in the three systems (Figure 4a). The taxa higher than the LDA significance threshold of 3.5 were scored in each system (Figure 4b). Various taxa were enriched in the BFT systems (Figure 4). For example, Oxyphotobacteria (Class), Cyanobacteria (Phylum), and Chloroplast (Order) were enriched in the *p*-BF treatment, whereas Gemmatimonadetes (Phylum) and planctomycetes (Class, Order, Family and Genus) were enriched in the M-BF treatment (Figure 4). Furthermore, Vibrionaceae (Family), Vibrionales (Order) and *Vibrio* (Genus) were enriched in the CW system, which was consistent with the analysis of variance in genus level (Figure 3d).

#### *3.3. Correlation Analysis between Environmental Factors and Microbial Community*

Redundancy analysis (RDA) was used to discern the possible correlations between the bacterial genera, the carbon varieties, and the environmental variables (Figure 5). Samples from the M-BF system were distributed on the left of the first canonical axis, whereas samples from the P-BF and CW systems were distributed on the right. Most environmental factors had a strong negative effect on the distribution of the samples on the first axis, whereas TAN and TA had a positive effect. Moreover, *Leucothrix* was strongly positively correlated with the M-BF treatment.

**Figure 3.** Bacterial community compositions and abundance at phyla (**a**) and genus (**b**) levels across all samples. The Analysis of Variance (ANOVA) shows the top 20 phyla (**c**) and genera (**d**) with the smallest *p*-value among treatments. \*, *p* < 0.05; \*\*, *p* < 0.01. CW, M-BF and P-BF represented the three systems of this study.

## *3.4. Functional Prediction in Carbon Adding Environments*

FAPROTAX was adopted to annotate and screen the key ecological functions of the bacterial communities in the three treatments. The top 10 functions are shown in Figure 6a. The functions of chemoheterotrophy and aerobic chemoheterotrophy were abundant in all treatments. Increased relative abundances of chloroplasts and intracellular parasites were observed in the P-BF treatment, while fermentation increased in the M-BF treatments. Based on the 16s rRNA gene sequences, the pairwise functional comparison was investigated using BugBase (Figure 6b–d). The relative abundance of several functions, such as biofilm formation (i.e., forms biofilms), anaerobic functions (i.e., anaerobic), stress tolerance (i.e., stress-tolerant), and mobile elements (i.e., contains mobile elements), increased in the M-BF than the CW system (Figure 6b). The relative abundance of the function of Gram-negative increased in both carbon adding systems compared to the CW system

(Figure 6b,c). The relative abundance of the functions related to Gram-negative, aerobic, anaerobic, forms biofilms and stress-tolerant were all higher in M-BF than P-BF (Figure 6d).

#### *3.5. Growth Performance*

The final weight, specific growth rate (SGR), growth rate, and survival rate of *L. vannamei* were highest in the P-BF system and lowest in the CW system. The three parameters in the P-BF system were significantly higher than those in the CW system (Table 2). The feed conversion ratios (FCR) were significantly lower in M-BF and P-BF treatments than in the CW treatment (*p* < 0.05).

**Figure 4.** Linear discriminant analysis effect size (LEfSe) analysis showing the abundance of three treatments. (**a**) LEfSe analysis identified the most differentially abundant taxa among three systems. The six rings of the cladogram stand for domain (innermost) phylum, class, order, family, and genus. (**b**) Linear discriminant analysis score of the three systems with a threshold value of 3.5. CW, M-BF and P-BF represented the three systems of this study.

**Figure 5.** Redundancy analysis (RDA) identified the correlation of different treatments, bacterial communities, and environmental factors. Biofloc: biofloc volume; TN: total nitrogen; TAN: total ammonia nitrogen; TA: total alkalinity; Chl: chlorophyll a; TOC: total organic carbon. Details of all environmental factors were shown in Table 1. CW, M-BF and P-BF represented the three systems of this study.

**Table 2.** Growth parameters and survival of *Litopenaeus vannamei* of different experimental systems at the end of the 34-d trial.


Data are presented as mean ± standard deviation (*n* = 3). CW, M-BF and P-BF represented the three systems of this study. In each row, different superscript letters indicate significant differences at the *p* < 0.05 level (one-way ANOVA and Tukey's test).

**Figure 6.** Functional prediction in different environments. (**a**) The predicted functions using FAPROTAX among different treatments. (**b**–**d**) Pairwise comparisons of functions between the three systems according to BugBase. The left part of each figure shows the mean proportion of two systems in different functions. The right part of each figure shows the difference between proportions. CW, M-BF and P-BF represented the three systems of this study.
