Isolation and Optimization of Aflatoxin B1 Degradation by Uniform Design and Complete Genome Sequencing of Novel Deep-Sea Kocuria rosea Strain 13

Aflatoxin B1 is a natural carcinogenic mycotoxin. The biological detoxification of aflatoxin could result in less environmental pollution, more moderate conditions, and less impact on food and feed, and be more convenient than physical and chemical methods. In this study, strain 13 with aflatoxin B1 degradation activity (67.47 ± 1.44%) was isolated and identified as Kocuria rosea. A uniform design was applied to optimize the degradation activity using a software Data Processing System, and a quadratic polynomial stepwise regression model was selected to investigate the relationships between the degradation rate and five independent variables. Furthermore, the optimal degradation conditions (culture temperature of 30 °C, culture time of 4.2 days, seawater ratio of 100%, pH of 7.11, and inoculation dosage of 0.09%) were verified with a degradation rate of 88 ± 0.03%, which was well matched with the predicted value (92.97%) of the model. Complete genome sequencing of Kocuria rosea, conducted with a combination of Illumina and single-molecule real-time sequencing, was used to analyze the genomic features and functions of the strain, which were predicted by the annotation based on seven databases, and may provide insights into the potential of Kocuria rosea, as well as providing a reference for degradation gene and protein mining. These results indicate that Kocuria rosea strain 13 has the ability to degrade aflatoxin B1 efficiently, and it also has the potential to provide aflatoxin-degrading enzymes.


Introduction
Aflatoxins (AFs) are derivatives of dihydrofuranoxanadione; they have one benzopyrone and one difuran ring and are mainly produced by the genera of Aspergillus [1]. Naturally occurring aflatoxins classified as Group 1 have been evaluated by the International Agency for Research on Cancer (IARC) of the World Health Organization and are regarded as carcinogenic to humans [2]. The most carcinogenic aflatoxin is aflatoxin B 1 (AFB 1 ), and it has been found that the consumption of food contaminated with AFB 1 causes immune suppression, deformity, gene mutagenesis, and carcinogenesis [3,4]. McMillan et al. reported that, when humans are exposed to AFB 1 at a dose of 20-120 µg/kg body weight per day for 1-3 weeks, acute aflatoxin poisoning occurs, which can cause abdominal pain, emesis, and even death [5,6].

Screening and Identification of Degrading Strain 13
The strains were isolated from the colonies, and it was shown that strain 13 could degrade aflatoxin B 1 with a rate of 67.47 ± 1.44% (with a culture temperature of 28 • C, culture time of 7 days, seawater ratio of 100%, pH of 7.52, and inoculation dosage of 1%) as optical density in 600nm (OD 600nm ) reached the value of 0.949 ± 0.016. The strain was Toxins 2023, 15, 520 3 of 16 identified according to a phylogenetic tree ( Figure 1) by 16S rRNA gene sequencing, which shows the relationship of strain 13 between Kocuria species with other related strains. The similarity of strain 13 with type stain Kocuria rosea DSM 20447 was 99.65%, which indicates that stain 13 could be identified as Kocuria rosea.

Screening and Identification of Degrading Strain 13
The strains were isolated from the colonies, and it was shown that strain 13 could degrade aflatoxin B1 with a rate of 67.47 ± 1.44% (with a culture temperature of 28 °C, culture time of 7 days, seawater ratio of 100%, pH of 7.52, and inoculation dosage of 1%) as optical density in 600nm (OD600nm) reached the value of 0.949 ± 0.016. The strain was identified according to a phylogenetic tree ( Figure 1) by 16S rRNA gene sequencing, which shows the relationship of strain 13 between Kocuria species with other related strains. The similarity of strain 13 with type stain Kocuria rosea DSM 20447 was 99.65%, which indicates that stain 13 could be identified as Kocuria rosea. The biological functions of Kocuria rosea strains were discovered, including the biosorption and biomineralization of U (VI) [41], dyes degradation and decolorization (methyl orange, amido black, methyl violet, cotton blue, and malachite green) [42,43], phenol biodegradation [44], Keratin hydrolysis [45], trinitrotoluene detoxification [46], and polyaromatic hydrocarbons degradation [47][48][49]. Naphthalene, anthracene, phenanthrene, fluorene, and pyrene, degraded by Kocuria rosea, have at least two benzene rings and have a similar structure with aflatoxin B1, which indicates that Kocuria rosea strains have the potential to degrade substances containing benzene ring structures.

Optimization for AFB1 Degradation
The results of the degradation rate from UD are shown in Figure 2. It was demonstrated that the degradation rate of N1, N2, N6, N10, and N11 with relatively high OD600nm was over 60%. In order to select the model with a significant fitting effect, three types of quadratic polynomial mathematical models were evaluated and compared with the four parameters outlined in Table 1. The adjusted coefficient of determination ( 2 ) represents the correlation between the observed values and the predicted values [50]. The closer 2 is to 1, the better the fitting effect achieved is. Root mean square error (RMSE) The biological functions of Kocuria rosea strains were discovered, including the biosorption and biomineralization of U (VI) [41], dyes degradation and decolorization (methyl orange, amido black, methyl violet, cotton blue, and malachite green) [42,43], phenol biodegradation [44], Keratin hydrolysis [45], trinitrotoluene detoxification [46], and polyaromatic hydrocarbons degradation [47][48][49]. Naphthalene, anthracene, phenanthrene, fluorene, and pyrene, degraded by Kocuria rosea, have at least two benzene rings and have a similar structure with aflatoxin B 1 , which indicates that Kocuria rosea strains have the potential to degrade substances containing benzene ring structures.

Optimization for AFB 1 Degradation
The results of the degradation rate from UD are shown in Figure 2. It was demonstrated that the degradation rate of N1, N2, N6, N10, and N11 with relatively high OD 600nm was over 60%. In order to select the model with a significant fitting effect, three types of quadratic polynomial mathematical models were evaluated and compared with the four parameters outlined in Table 1. The adjusted coefficient of determination (R 2 adj ) represents the correlation between the observed values and the predicted values [50]. The closer R 2 adj is to 1, the better the fitting effect achieved is. Root mean square error (RMSE) represents the differences between predicted and observed values and the precision of the predicted model [51,52]. Akaike's information criterion (AIC) was derived from an asymptotic approximation to the Kullback−Leibler divergence between the true distribution and the model, and the Bayesian information criteria (BIC) derived from the dominant terms in the Laplace approximation to the logarithm of the Bayes factor with a vague prior [53]. Both AIC and BIC are two parameters commonly used for model selection, which were first introduced by Akaike [54] and Schwarz [55], respectively. However, the AIC assumes that the true model is not considered in all models, but the BIC assumes that the true model is one of the models [56]. The smaller the values of RMSE, AIC, and BIC of the models, the better the fitting effects which were demonstrated for the models were. The R 2 of the quadratic polynomial stepwise regression model and stepwise regression model with multiple factors and interaction terms were much closer to one compared to the multivariate and squared stepwise regression model. Moreover, the RMSE, nAIC, and BIC of the quadratic polynomial stepwise regression model were the minimums in the three models. Therefore, it was indicated that the quadratic polynomial stepwise regression model had a better fitting effect than the other two models.
predicted model [51,52]. Akaike's information criterion (AIC) was derived from an asymptotic approximation to the Kullback−Leibler divergence between the true distribution and the model, and the Bayesian information criteria (BIC) derived from the dominant terms in the Laplace approximation to the logarithm of the Bayes factor with a vague prior [53]. Both AIC and BIC are two parameters commonly used for model selection, which were first introduced by Akaike [54] and Schwarz [55], respectively. However, the AIC assumes that the true model is not considered in all models, but the BIC assumes that the true model is one of the models [56]. The smaller the values of RMSE, AIC, and BIC of the models, the better the fitting effects which were demonstrated for the models were. The R 2 of the quadratic polynomial stepwise regression model and stepwise regression model with multiple factors and interaction terms were much closer to one compared to the multivariate and squared stepwise regression model. Moreover, the RMSE, nAIC, and BIC of the quadratic polynomial stepwise regression model were the minimums in the three models. Therefore, it was indicated that the quadratic polynomial stepwise regression model had a better fitting effect than the other two models.  The factors x2 × x3 and x2 × x2 in the quadratic polynomial stepwise regression model were removed for values of p greater than 0.05, which were 0.2404 and 0.2778, respectively. The formula of the model was generated as follows: where y represents the degradation rate (%); x1 represents the culture temperature (°C); x2 represents the culture time (days); x3 represents the seawater ratio (%); x4 represents pH; and x5 represents the inoculation dosage (%).
Multi-way analysis of variance (ANOVA) ( Table 2) was applied in model evaluation, showing that the p-value was below 0.05, which indicated the great predictive ability of the model. The correlation coefficient (R), determinate coefficient (R 2 ), and adjusted  The factors x2 × x3 and x2 × x2 in the quadratic polynomial stepwise regression model were removed for values of p greater than 0.05, which were 0.2404 and 0.2778, respectively. The formula of the model was generated as follows: where y represents the degradation rate (%); x1 represents the culture temperature ( • C); x2 represents the culture time (days); x3 represents the seawater ratio (%); x4 represents pH; and x5 represents the inoculation dosage (%).
Multi-way analysis of variance (ANOVA) ( Table 2) was applied in model evaluation, showing that the p-value was below 0.05, which indicated the great predictive ability of the model. The correlation coefficient (R), determinate coefficient (R 2 ), and adjusted determinate coefficient (adj. R 2 ) of the equation were 0.999, 0.999998, and 0.99978, respectively, which also indicated that the model could well reflect the relationship among culture temperature, culture time, seawater ratio, pH, and inoculation dosage. The relationship between the predicted values and observed values of the degradation rate are also confirmed in Figure 3, which shows that most points were distributed along a straight line, indicating that the predicted values and observed values were very close. Furthermore, a two-way ANOVA and multiple comparisons of least significant difference (LSD) for the predicted value and observed value were performed, demonstrating that the observed value of three replications (p = 0.988, 0.851, and 0.865 > 0.5) had no significant differences from the predicted value. Therefore, the quadratic polynomial stepwise regression model effectively estimated the degradation rate of strain 13 cultured under different conditions. determinate coefficient (adj. R 2 ) of the equation were 0.999, 0.999998, and 0.99978, respectively, which also indicated that the model could well reflect the relationship among culture temperature, culture time, seawater ratio, pH, and inoculation dosage. The relationship between the predicted values and observed values of the degradation rate are also confirmed in Figure 3, which shows that most points were distributed along a straight line, indicating that the predicted values and observed values were very close. Furthermore, a two-way ANOVA and multiple comparisons of least significant difference (LSD) for the predicted value and observed value were performed, demonstrating that the observed value of three replications (p = 0.988, 0.851, and 0.865 > 0.5) had no significant differences from the predicted value. Therefore, the quadratic polynomial stepwise regression model effectively estimated the degradation rate of strain 13 cultured under different conditions.  Parameter estimation and significance test of the model were performed. It was shown that the p-values of the ten factors were all less than 0.05, indicating that these ten factors significantly affected the degradation rate (Table 3). According to standard regression coefficients, factors x1, x3 × x3, x2, x1 × x3, x5 × x5, and x1 × x2 had, in descending order, significant positive effects on the degradation rate. Additionally, factors x1 × x1, x2 × x3, x3, and x5 had, in descending order, significant negative effects on the degradation rate. The predicted model for interaction terms varying within the experimental range was visualized through response surface plots and contour plots, and other variables remained at the optimal level ( Figure 4). The optimal temperature was 30 °C no matter what the culture time and seawater ratio was, as shown in Figure 4a,b. The contribution of x1 (standard regression coefficient) to the equation was 4.86; however, the contribution of x1 × x2 and x1 × x3 was only 0.03 and 0.20, respectively (Table 3), which explained that x1 had a more significant influence in x1 × x2 and x1 × x3 than x2 and x3 (Figure 4a,b). It was also indicated that culture temperature had the most significant impact on the degradation rate compared to other factors (Table 3). Moreover, although x2 × x3 had an extremely negative effect (p < 0.01) on the degradation rate with a standard regression coefficient of −0.50, x2 had a significantly positive effect (p < 0.01) with a standard regression coefficient of 0.55. Consistently, it was observed that higher culture time was beneficial to the degradation rate (Figure 4c). Parameter estimation and significance test of the model were performed. It was shown that the p-values of the ten factors were all less than 0.05, indicating that these ten factors significantly affected the degradation rate (Table 3). According to standard regression coefficients, factors x1, x3 × x3, x2, x1 × x3, x5 × x5, and x1 × x2 had, in descending order, significant positive effects on the degradation rate. Additionally, factors x1 × x1, x2 × x3, x3, and x5 had, in descending order, significant negative effects on the degradation rate. The predicted model for interaction terms varying within the experimental range was visualized through response surface plots and contour plots, and other variables remained at the optimal level ( Figure 4). The optimal temperature was 30 • C no matter what the culture time and seawater ratio was, as shown in Figure 4a,b. The contribution of x1 (standard regression coefficient) to the equation was 4.86; however, the contribution of x1 × x2 and x1 × x3 was only 0.03 and 0.20, respectively (Table 3), which explained that x1 had a more significant influence in x1 × x2 and x1 × x3 than x2 and x3 (Figure 4a,b). It was also indicated that culture temperature had the most significant impact on the degradation rate compared to other factors (Table 3). Moreover, although x2 × x3 had an extremely negative effect (p < 0.01) on the degradation rate with a standard regression coefficient of −0.50, x2 had a significantly positive effect (p < 0.01) with a standard regression coefficient of 0.55. Consistently, it was observed that higher culture time was beneficial to the degradation rate (Figure 4c).
The optimal degradation conditions were predicted as a culture temperature of 30 • C, culture time of 4.2 days, seawater ratio of 100%, pH of 7.1094, and inoculation dosage of 0.0899%, with a degradation rate of 92.97%. Additionally, the confirmation experiments showed that the degradation rate was 88 ± 0.03%. Furthermore, the results of the single-sample t-test (t = −2.845; p = 0.105 > 0.05) showed that the original hypothesis (H0: µ = 92.97%) could not be rejected, which demonstrated that there were no significant differences between the predicted value and observed value.

General Genomic Features of Strain 13
To better understand the AFB 1 degradation mechanisms of strain 13, the complete genome of strain 13 was sequenced and mined. A graphical circular genome map of strain 13 is shown in Figure 5. The complete genome sequences of Kocuria rosea 13 were assembled into four scaffolds, including a chromosome and three plasmids. The chromosome had a size of 3,815,108 bp and a GC content of 72.86%. The predicted coding sequence has 3797 genes with a total length of 3,692,208 bp accounting for 91% of the complete genome, which also has 72.6% GC in the gene region. Additionally, there were 136 tandem repeats predicted with a ratio of 48% in the genome. Moreover, 91 RNA genes were predicted: 50 tRNA genes, 32 sRNA, and 9 rRNA genes, including 3 16S rRNA genes, 3 23S rRNA, and 3 5S rRNA genes. Regarding mobile genetic elements, six gene islands, eight clustered regularly interspaced short palindromic repeat (CRISPR)-Cas systems, and one prophage were predicted.

Gene Function Analysis
The CDS genome was annotated according to the following seven databases: Non-Redundant Protein Database (NR), Swiss-Prot, evolutionary genealogy of genes: Nonsupervised Orthologous Groups (EggNOG), Pfam, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Carbohydrate-active enzymes (CAZy). The

Gene Function Analysis
The CDS genome was annotated according to the following seven databases: Non-Redundant Protein Database (NR), Swiss-Prot, evolutionary genealogy of genes: Nonsupervised Orthologous Groups (EggNOG), Pfam, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Carbohydrate-active enzymes (CAZy). The results of annotations for NR (Table S1), Swiss-Prot (Table S2), and Pfam (Table S3) were related to 3768, 2690, and 3120 genes, respectively.
EggNOG annotations divided 3231 genes, accounting for 85.09% of all the genes of strain 13 into 20 categories (Figure 6a). Type G Carbohydrate transport and metabolism had 241 genes which might be related to AFB 1 degradation. Additionally, 80 genes annotated for secondary metabolites biosynthesis, transport, and catabolism could also be related to toxin degradation. Different gene numbers were counted from 1 to 274 for various function types; however, there were still unknown functions of 909 genes.  GO analysis classified 2698 genes (71.06% of all genes) from strain 13 into three major categories, including biological process (1149 genes), cellular component (1175 genes), and molecular function (2195 genes) (Figure 6b and Table S4). In the biological process, the GO annotations of top five genes were regulation of transcription and DNA-templated (GO ID: 0006355), translation (GO ID: 0006412), transmembrane transport (GO ID: 0055085), carbohydrate metabolic process (GO ID: 0005975), and methylation (GO ID: 0032259), which were related to 72 genes, 57 genes, 52 genes, 46 genes, and 37 genes, respectively. The biological process annotations of strain 13 contained 428 sub-functions, and molecular In KEGG annotations, there were six primary categories (organismal systems, environmental information processing, human diseases, cellular processes, genetic information processing, and metabolism) of KEGG pathways corresponding to 1748 genes of strain 13, and each category contained different numbers of pathways (Figure 6c). Most numbers of genes (711 genes) were in connection with the global and overview maps in the largest category metabolism. Moreover, carbohydrate metabolism (243 genes), biosynthesis of other secondary metabolites (39 genes), and xenobiotics biodegradation and metabolism (78 genes) might be related to AFB 1 degradation.
Carbohydrate-active enzymes (CAZyme) contain auxiliary activities (AAs), carbohydratebinding modules (CBMs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), glycoside hydrolases (GHs), and glycosyl transferases (GTs) which could degrade, modify, and generate glycosidic bonds. To reveal the mechanism of the microbial carbohydrate metabolism, CAZy was used for the prediction and classification of CAZyme in stain 13. Only four types of CAZymes were identified from the complete genes of the strain, which were AAs (17 genes), CEs (22 genes), GHs (43 genes), and GTs (43 genes) (Table S5). There were the highest gene counts for enzymes in the families of GHs and GTs, which played a pivotal part in the degradation of polymers.

Conclusions
In summary, Kocuria rosea strain 13 was found to degrade AFB 1 (88 ± 0.03%) in optimized conditions (culture temperature of 30 • C, culture time of 4.2 days, seawater ratio of 100%, pH of 7.1094, and inoculation dosage of 0.0899%). Therefore, this study suggests that Kocuria rosea could be used for aflatoxin degradation. Moreover, the annotations of the genome predicted the potential of the strain, and some genes might be related to degradation mechanisms, which could be further screened and verified by transcriptomics techniques in subsequent research. Future work could also focus on the identification and toxicity assessment of the degradation products metabolized by the strain.

Isolation of the Strain 13
Coumarin was added into a tube as the solid substrate, and the tube was wrapped up with nylon mesh in a case consumed by deep-sea organisms. Additionally, sterilization was performed at 115 • C for 30 min. The wrappage was placed in sterilized incubation chambers (ICs) of the deep-sea in situ microbial incubator (DIMI), which was placed at a flat-topped seamount in the West Pacific Ocean (N 20.4059567 • , E 160.7700883 • ) at 1617 m depth for 348 days of cultivation. The samples were collected and diluted tenfold with sterile seawater. The suspension was enriched for 20 days at 150 rpm with a shaker at 20 • C in an M2 liquid culture medium with 1 µg mL −1 AFB 1 . Additionally, 1 mL of the suspension was cultured with the same conditions for the second enrichment in an M2 seawater medium with 5 µg mL −1 AFB 1 . Then, the final enrichment cultures were diluted in a gradient, spread onto the M2 culture plate with AFB 1 as the only carbon source, and cultured at 20 • C. To obtain pure stains, different colonies were picked, isolated, and incubated on M2 plates separately. The growth rate of the strains in the M2 liquid culture medium was monitored by a UV-2000 spectrophotometer (Unico, Shanghai, China) with OD 600nm . As the strains grew and OD 600nm reached about 1.0, the suspensions were inoculated into an M2 liquid medium (AFB 1 as the only carbon source) and cultured for seven days at 28 • C.

Determination of Aflatoxin Degradation Rate
AFB 1 in the control groups and samples was extracted three times using chloroform with an equal volume before the solvent evaporated under N 2 at room temperature [23,26,57,58]. Dimethyl sulfoxide (DMSO) (50 µL) was used to dissolve the dried extracts, and 20 µL of the mixture was injected in UltiMateTM3000 HPLC (Thermo Scientific, Bremen, Germany). HPLC analysis was conducted with a C18 Polaris column (250 mm × 4.6 mm i.d., 5 µm) in a mobile phase of water and methanol in a 1:1 ratio (v/v). The flow rate was set as 1 mL min −1 , and a UV/VIS detector (Thermo Scientific, Germany) was used for absorbance measurements at a wavelength of 360 nm. The column temperature was set to 35 • C for detection. The software Chromeleon v6.8 was used for data analysis. The rate of AFB 1 degradation was determined and calculated with (1 − AFB 1 peak area in treatment/AFB 1 peak area in control) × 100%.

UD for Aflatoxin Degradation Optimization
Aflatoxin degradation by the strain supernatant was optimized under UD according to [59]. Five independent variables were selected as follows: culture temperature (x1) ( • C), culture time (x2) (days), seawater ratio (x3) (%), pH (x4), and inoculation dosage (x5) (%) ( Table 4). Multiple mixed-level uniform designs U 12 (4 1 × 6 4 ) were obtained with the software Data Processing System (DPS) 18.10 [60] by varying the parameters of random seed number, maximum iterations, and optimal search time. Seven UD matrix performance parameters [59] were compared among different UD tables, resulting in a UD (Table 5) with the smallest values of parameters being selected for the experimental scheme. To assess degradation performance, the degradation rate was used as the dependent variable. The dependent variable could be related to the above five independent variables through three quadratic polynomial mathematical models: a quadratic polynomial stepwise regression model, a stepwise regression model with multiple factors and interaction terms, or a multivariate and squared stepwise regression model. The models were described using Equations (1)-(3).
Quadratic polynomial stepwise regression model: Stepwise regression model with multiple factors and interaction terms: Multivariate and squared stepwise regression model: where y represents the dependent variable to be modeled; x i and x j represent the independent variables; b ij , b ii , b i , and b 0 represent the interaction coefficients, quadratic coefficients, linear coefficient, and constant coefficient, respectively. The fit goodness of the three models was evaluated and compared using R 2 adj , RMSE, nAIC, and BIC. The equations of R 2 adj , RMSE, nAIC, and BIC can be described as follows: where O i represents the ith measured observed value; P i represents the ith predicted value; m represents the average value; RSS represents the residual sum of squares; N represents the number of values in the estimation data set; n p represents the number of estimated parameters; and n y represents the number of model outputs.

DNA Extraction, Amplification, and 16S rRNA Gene Sequencing
The DNA of aflatoxin-degrading bacteria was extracted via a boiling lysis method: a single colony was selected and boiled in a 100 • C water bath for 15 min, then placed at 4 • C for 30 min, spun at 5000 rpm for one minute, and the extracted DNA in the supernatant was used for PCR amplification. An initial denaturation at 94 • C for 5 min, followed by 30 cycles of denaturation at 94 • C for 40 s, primer annealing at 55 • C for 40 s, extension for 1 min at 72 • C, and a final extension at 72 • C for 10 min, was used to amplify the DNA sample in a reaction mixture of 50 µL. The PCR product was purified and sent to Sangon Biotech Co., Ltd. (Shanghai, China) for sequencing. The sequences of the strain were assembled with the software DNAMAN 9.0.1 and aligned in EZTaxon (https://www.ezbiocloud.net/; accessed on 16 January 2023). The phylogenetic tree of the strain was constructed with the software MEGA 10.2.5 using the neighbor-joining (NJ) method. The bootstrap method replications were set as 1000. The 16S rDNA sequence of strain 13 (1521 bp) uploaded to GenBank was registered for the accession number CP127857. The bacterial cells were cultured to the logarithmic growth phase and collected via centrifuge CR21N (HITACHI, Tokyo, Japan) for 5 min at 14,000 rpm. The genomic DNA was extracted and purified with the Wizard ® genomic DNA purification kit (Promega Corp., Madison, WI, USA), and the purity and concentration were detected using agarose gel electrophoresis and Nanodrop 2000 (Thermo Scientific, Germany), respectively.

Genomic Library Construction and Sequencing
The genome of the strain was sequenced with both Illumina sequencing and singlemolecule real-time sequencing (SMRT). The Illumina data were used to assess genomic heterozygosity, genomic size, genomic duplication, presence of plasmids, and contamination, in addition to correcting long sequences from the third generation of sequencing to ensure the completeness and accuracy of the assembly.
For the SMRT platform, genomic DNA in a 15 µg sample was sheared into 8-10 kb fragments by a centrifuge 5424 (Eppendorf, Hamburg, Germany) at 6000 rpm for 1 min with a G-tube (Covaris, America). Both ends of the purified single-strand DNA fragments were connected with a sequencing adapter named SMRT bell for library construction. The Agencourt AMPure XP kit (Beckman Coulter Genomics, Woollahra, NSW, Australia) with 0.45× volumes was applied for library purification three times, and the library was sequenced by a PacBio RS II (Pacific Biosciences of California Inc., Menlo Park, CA, USA).

Genome Assembly and Plasmid Identification
The complete genome sequences, including plasmids, were assembled with the reads of both Illumina and PacBio. The Illumina raw data saved as a FASTQ file were trimmed, and low-quality reads were removed for clean data. The software Unicycler v0.4.8 [61] was used for the assembly of the PacBio reads before the reads were corrected according to Illumina reads with Pilon v1.22 for the complete genome, including chromosome and plasmid sequence. In addition, PlasFlow (https://github.com/smaegol/PlasFlow; accessed on 16 January 2023) was used for plasmid identification. Furthermore, the plasmid sequences were annotated using the basic local alignment search tool (BLAST) and database PLSDB (https://ccb-microbe.cs.uni-saarland.de/plsdb/; accessed on 16 January 2023). The genome sequences were stored at GenBank with the accession numbers CP127857 (chromosome), CP127858 (plasmid A), CP127859 (plasmid B), and CP127860 (plasmid C).

Genome Function Annotation
The CDS genome was compared and annotated with different databases using the following software: Diamond v0.8.35 for NR of the National Center for Biotechnology Information (NCBI), Swiss-Prot [65], and EggNOG v4. 5