The samples were collected from the hot springs located at the Tolhuaca thermal area (38°13′59.09”S, 71°43′20.3”W). Situated in the pre-Andean mountain range at 1,150 m.a.s.l., and ~33 km from the city of Curacautín (Araucanía Region of Chile). This geothermal region is characterized by high volcanic activity, influenced by nearby volcanoes such as Tolhuaca, Lonquimay, Llaima, and Sierra Nevada. In the study zone, a 500 mL thermal water at ~96°C was collected aseptically into sterile containers and transported to the laboratory under controlled conditions for further processing.

The bacterial strain Tol1 was isolated from the water samples described above. The samples were subjected to a 10-fold serial dilution in sterile water. Aliquots (50 μL) from each dilution were spread on petri dishes containing 23 g L⁻¹ of Nutrient Agar (Difco™) and incubated at 55°C for 24–48 h. Distinct mucoidal isolated colony was picked and reisolated by successive streaking on fresh agar plates until pure cultures were obtained, showing a shiny, mucoid appearance due to EPS production.

Genomic DNA from strain Tol1 was isolated from freshly grown bacteria (~0.5 g) in NA plates using DNA power soil kit (Sigma Aldrich) according to the manufacturer’s instructions. Although the culture was purified, we used the DNA PowerSoil Kit (Sigma-Aldrich) as it has been demonstrated to efficiently lyse robust Gram-positive bacterial cells, such as Bacillus species, which can be difficult to lyse using standard kits due to their thick peptidoglycan-rich cell walls. This choice ensured high-quality and high-yield DNA (396.3 ng μL⁻¹, A260/A280 = 1.87) suitable for downstream genome sequencing. DNA concentration and quality were assessed using NanoDrop 1000 (Thermo Fisher Scientific) and Qubit 4.0 (Thermo Fisher Scientific). The genome sequencing was performed using Illumina Novaseq PE150 (Illumina, Inc., San Diego, CA, USA) platform using paired-end sequencing strategy. Raw reads were filtered, and high-quality paired-end reads were assembled using SPAdes (v3.15.5) (Bankevich et al., 2012). CheckM was used to evaluate strain Tol1 genome quality (Parks et al., 2015). The identification of rRNA genes was carried out using RNAmmer (version 1.2) (Lagesen et al., 2007). The taxonomic identification of Tol1 was carried out using the 16S rRNA sequence as a phylogenetic marker, obtained through the barrnap software (Galaxy Version 1.2.2) from the sequence of the complete genome and compared using the EzBioCloud server (Yoon et al., 2017) and NCBI BLAST-n. The genome was visualized using Proksee (Grant et al., 2023) and compared with the closest species using BLAST (Altschul et al., 1990). Prodigal was used to identify open reading frames (Hyatt et al., 2010). The genome sequence was submitted to GenBank under the accession number JBNNEF000000000.

B. licheniformis Tol1 genome was categorized in a subsystem using an open-source prokaryotic genome annotation system, Rapid Annotation System Technology (RAST) pipeline, compared with the entire available dataset (Aziz et al., 2008). The functionally annotated genes present within the genome were visualized using SEED Viewer. PHAge Search Tool Enhanced Release (PHASTEST) web-based server was used to identify genes related to prophage with the gnome of Tol1 (Wishart et al., 2023). To predict genomic islands, IslandViewer4 was used (Bertelli et al., 2017).

The biofilm formation experiment was performed in 50 mL centrifuge tubes using 20 mL of growth medium and 200 μL of inoculum. Over a period of days, we observed the growth of bacteria and the production of EPS from strain Tol1 using confocal microscope (Leica Stellaris 5, Wetzlar, Germany). The sample was stained with DAPI (4′,6-diamidino−2-phenylindole, a fluorescent stain used to label DNA which fluoresces at 358/461 nm) to visualize the bacterial cells, and SYPRO™ Ruby (highly sensitive, detecting glycoproteins with a fluorescence spectrum of 450/610 nm) to visualize the biofilm matrix.

A bacteriological analysis was performed to determine the optimal culture conditions that maximize EPS production. A one-factor-at-a-time (OFAT) assay was performed under different conditions where EPS production can be visualized by viscosity. The strain was inoculated into LB Agar medium supplemented with different carbon sources (galactose, mannose, glucose, and sucrose), nitrogen sources (yeast extract, peptone, meat extract, and ammonium chloride), different initial pH (5, 6, 7, and 8), temperature (50, 55, 60, 65, 70°C), and finally salinity (NaCl) (2, 5, 7, and 10%). Each experiment was performed in triplicate to reduce experimental error. To precipitate the EPS, bacterial cells in the stationary growth phase were harvested, followed by treatment with 4% trichloroacetic acid (TCA) (w/v) for 30 min (Sarkar et al., 2024). Centrifugation at 4°C, 5000 × g for 20 min was applied to the TCA-precipitated cells to achieve protein precipitation. Then, the chilled, cell-free supernatant was left overnight at 4°C after adding an equal volume of acetone. The next day, centrifugation at 12,000 × g for 20 min was applied to the solvent-coagulated EPS. The acetone precipitation step was repeated with water-dissolved EPS to remove residual TCA and proteins. Pure EPS powder was obtained after dialysis and lyophilization of purified EPS.

In this study, RSM based on central composite design (CCD) was used to optimize three parameters (sucrose concentration %, NH₄Cl %, and NaCl %) with five levels of each parameter coded as (−2, −1, 0, +1, +2) were used as independent parameters and EPS production as response. Table 1 presents the coded and uncoded values of independent parameters, coded according to the following Equation 1: (1) xi = ( X i − X o / ΔX )

Where: xi: the coded value of the parameter, X_i: the actual value of the parameter, X_o: the actual value of X_i at center point, ΔX: step change.

The CCD matrix was generated by Design Expert software (version 7) consisting of twenty runs with five replicates of center point (Supplementary Table 1). The experiment was conducted in 200 mL glass flasks with 100 mL of medium. Flasks were prepared and inoculated in triplicate to ensure adequacy in all experiments. The inoculated flasks were incubated at 55°C and 150 rpm for 3 days in a rotary shaker. EPS production was achieved, and the data were analyzed in design expert software to evaluate the effects of each parameter and their interaction to predict the production of EPS as described in the quadratic optimization model (Equation 2): (2) Y = β 0 + ∑ i = 1 k β 0 x i + ∑ i = 1 k β ii x i 2 + ∑ i = 1 k ∑ i = 2 k β ij x i x j

Where Y is the predicted response, β₀ is the intercept term, β_i is the linear effect, β_ii is the squared effect, β_ij is the interaction effect, and X_i and X_j are input variables that influence the response variable Y.

Design expert software is a beneficial tool utilized for designing experiments and data analysis of their responses (Suryawanshi et al., 2019). The optimization model was statistically evaluated to assess analysis of variance (Table 2). Model fitness expressed by R², R _adjust, R _predict, adequate precision, standard deviation, coefficient of variation, mean, press, F- value, P > F prob. and lack of fit, other statistical terms used for statistical analysis (Breig and Luti, 2021; Sarkar et al., 2024). Depending on the regression model, an optimization plot can be developed by using Design expert software to identify the optimum concentrations [sucrose concentration %, NH₄Cl %, and NaCl %] that maximize EPS production, briefly, all independent parameters selected as in range mode except the initial concentration of dye selected as maximum mode and response EPS production as maximize to construct ramp chart for optimization description. The optimum parameters obtained were validated in the laboratory in order to verify the predicted EPS production.

Modeling is a useful technique for successfully and economically determining ideal operating conditions and forecasting process results. The normal method for modeling a traditional chemical reactor is to use standard mathematical tools to solve the mass, energy, and momentum balance equations concurrently. Numerous articles document the wide exploration of the use of Artificial Neural Networks (ANN) for process modeling in a variety of scientific and engineering fields (Breig and Luti, 2021; Mahdi et al., 2021; Alardhi et al., 2022; Banerjee et al., 2022; Alardhi et al., 2024). Similar to biological neural networks, artificial neural networks are dynamic, non-linear systems made up of numerous processing units (neurons) coupled to one another. Feedforward neural networks and recurrent neural networks are the two primary categories of ANN. In a feedforward neural network, input data only moves from the input nodes to the output nodes in a forward manner (Al-Jadir et al., 2023; Dawood Salman et al., 2023). Every node in a layer is linked to nodes in the layers above and below it, resulting in a unidirectional flow devoid of backward connections (Sadeghzadeh Ahari et al., 2024). The Multilayer Perceptron (MLP), the most common type of feedforward neural network, usually consists of one or more hidden layers, an input layer, and an output layer. The amount of input variables and the number of variables that need to be predicted correlate with the number of nodes in the input and output layers, respectively (Ayed et al., 2022). An iterative process of trial and error is frequently required to determine the best configuration for the hidden layers, their nodes, and the transfer functions between levels (Annadurai and Lee, 2007). By using an iterative process, several types of networks can be explored with the ultimate goal of finding a structure with fewer layers and nodes that reduce prediction errors (Alardhi et al., 2023b). In this work, the reliability of the ANN models was tested using several indices in addition to the expected and actual results. These markers are MSE and RE, as can be seen in the following Equations 3, 4 (Alardhi et al., 2023a). (3) MSE = 1 n ∑ i = 1 n ( D a ( t ) − D f ( t ) ) (4) RE = D a ( t ) − D f ( t ) D a ( t ) × 100

Where: D a ( t ) and D f ( t ) are the actual and predicted values at time t.

ANNs were utilized to model the EPS production using MATLAB R2021a software. The neural network was trained by the experimental data set, which consisted of 80 data points and was gathered at the lab scale to estimate the EPS production percent. Several factors of bacterial growth, such as % sucrose (1.32, 2, 3, 4, 4.68) as carbon source, % NH₄Cl (1.32, 2, 3, 4, 4.68) as nitrogen source, % NaCl (2, 2.32, 3, 4, 5 and 5.68) as inputs, were initially included in the experimental study before being included in the modeling.

An MLP is composed of an input layer, an output layer, and one or more hidden layers, as illustrated in Figure 1. Each layer is composed of a specific quantity of neurons; the number of neurons in the input and output layers correspond to the number of input and output variables, respectively. The values of neurons in both the hidden and output layers are contingent upon the weights assessed throughout the training phase. A threshold signal is assigned to a weight and incorporated into each neuron within the hidden and output layers. The construction of ANN models involves three stages: training, validation, and testing. Consequently, the data set is divided into three categories: training, validation, and testing. The majority of the data resides in the training dataset, often comprising 70% of the total data. The residual data is allocated between the validation and testing data sets (15 and 15%) respectively.

The genome of Bacillus licheniformis Tol1 identified through 16 s rRNA sequence analysis, as annotated through the RAST subsystem, exhibits a comprehensive set of features indicative of its genetic complexity and adaptability. The genome has a total size of 4, 247, 584 base pairs with a GC content of 45.9%, reflecting a balanced composition conducive to thermophilic adaptation (Hu et al., 2022). The assembly metrics highlight a high-quality genome, with an N50 of 2,295,028 and an L50 of 1, indicating that the genome is represented mainly in a single, contiguous segment. The analysis identified 56 contigs containing protein-encoding genes (PEGs), with a total of 4,684 coding sequences and 90 RNA genes, including essential components for transcriptional and translational regulation. From BLAST-n, the highest degree of similarity was found with Bacillus sp. (in: Bacteria) strain licheniformis 16S ribosomal RNA gene, partial sequence (sequence ID: MH900194.1) with 100% identity percentage. The 16S rRNA sequence of Tol1 was also identified from EzBioCloud server, where the sequence showed highest degree of similarity with B. licheniformis type strain ATCC 14580 having 99.92% similarity.

The circular genome annotation of B. licheniformis Tol1 (Figure 2A) reveals a metabolically versatile and functionally rich genome, with multiple open reading frames (ORFs) identified across both forward and reverse strands. Notably, the genome contains critical antibiotic resistance genes, including vancomycin resistance clusters (vanT, vanW, vanY) within the vanG and vanI operons, a β-lactamase gene (BcIII), and qacG, associated with resistance to quaternary ammonium compounds. Additionally, the presence of FosBx1 indicates resistance to fosfomycin, highlighting the strain’s potential multi-drug resistance. A Type I-A CRISPR-Cas system (cas8a1a2_TypeIA) was also identified, suggesting adaptive immunity against foreign genetic elements like plasmids and bacteriophages. The GC content and skew analysis point to regions with functional significance and potential horizontal gene transfer. Collectively, these features underscore the bacterium’s resilience and adaptability to extreme environments while also providing insights into its biotechnological potential and implications for antimicrobial resistance. The genome of B. licheniformis Tol1 harbors four prophage regions (Figure 2B), three of which are intact and one incomplete, showcasing its dynamic evolutionary history. The intact regions include a 51.5 Kb prophage related to PHAGE_Bacillus_phi105 (GC content 40.10%) and a 54.4 Kb prophage associated with PHAGE_Bacillus_SPbeta (GC content 36.05%), both of which are known temperate phages commonly found in Bacillus species. These intact prophages are fully capable of encoding functional phage particles, suggesting their potential inducibility under stress conditions, contributing to horizontal gene transfer and genome plasticity. The 32.9 Kb intact prophage related to PHAGE_Brevibacterium_Jimmer2 (GC content 47.59%) highlights possible interspecies genetic exchange, further enhancing the bacterium’s adaptability. The 34.1 Kb incomplete prophage, a remnant of PHAGE_Bacillus_SPbeta (GC content 33.82%), may still encode functional genetic fragments that influence host fitness. The diversity of these prophages and their associated genes underscores their potential roles in conferring beneficial traits, such as antibiotic resistance, stress tolerance, or novel metabolic capabilities, aiding B. licheniformis Tol1 surviving within the extreme conditions of a hot spring.

The genome of B. licheniformis Tol1, annotated using the RAST subsystem, reveals a diverse array of functional features that underscore its metabolic versatility, stress resilience, and potential for biotechnological applications. Among the mapped subsystems, carbohydrate metabolism (314 genes) and amino acid biosynthesis and derivatives (312 genes) dominate, highlighting the bacterium’s capacity to utilize various energy and nutrient sources vital for survival in the extreme conditions of the Tolhuaca hot spring. Dormancy and sporulation pathways (95 genes) are detected, including spore core dehydration, sporulation clusters, biofilm matrix protein TasA, and spore germination genes, all of which contribute to its ability to withstand thermal, osmotic, and nutrient stress. A strong oxidative stress response (17 genes), including protection from reactive oxygen species and glutathione metabolism, coupled with osmoregulation and choline/betaine biosynthesis (11 and 10 genes, respectively), further ensures the organism’s resilience in challenging environments. Interestingly, it has genes for EPS biosynthesis, including the tyrosine-protein kinase EpsD and its modulator EpsC, which are essential for biofilm formation and extracellular polymer production. Stress-specific responses include genes for sigmaB regulation (6), carbon starvation (2), and detoxification (4), providing mechanisms to mitigate cellular damage and sustain metabolic activity during environmental fluctuations. In addition to stress response, B. licheniformis Tol1 encodes 20 genes associated with antibiotic and toxic compound resistance, such as fluoroquinolones, beta-lactams, and heavy metals like cobalt, cadmium, and copper, indicating its adaptability and potential for bioremediation. The genome also reveals a robust system for DNA repair (47 genes), including RecA-mediated homologous recombination, UvrABC excision repair, and pathways for non-homologous end-joining and the RecBCD and RecFOR systems. Phage-related genes (16), including those for replication, capsid formation, and packaging machinery, highlight horizontal gene transfer as a driver of genomic plasticity and evolutionary adaptability. Lastly, genes associated with plant hormone biosynthesis (4 auxin biosynthesis genes) suggest its role as a potential plant growth-promoting bacterium. At the same time, no virulence factors were identified, emphasizing its safety for industrial and agricultural applications. Collectively, these features underline B. licheniformis Tol1’s potential for sustainable biotechnological processes, from biopolymer production to environmental resilience and agricultural enhancement.

The biofilm-forming ability of the bacteria B. licheniformis Tol1 under varying pH conditions (5.8, 6.8, 7.8) and incubation temperatures (45°C, 55°C, and 65°C) were evaluated using confocal microscopy. SYTO9/Propidium Iodide staining was employed to assess cell viability within the biofilm, while DAPI/SYPRO Ruby staining was used to visualize and quantify (Figure 3) the extracellular polymeric substance (EPS) matrix. At 45°C, B. licheniformis Tol1 formed robust biofilms with high viability. Across all pH conditions, viable cells outnumbered apoptotic ones, particularly at pH 6.8 and 7.8, where viable cell counts exceeded 5,000, suggesting favorable biofilm formation conditions. EPS production was also high at this temperature, with the EPS matrix area ranging from ~1,500 to 2,200 μm², supporting the visual observation of dense, well-structured biofilm mats. Further, at 55°C, the strain maintained strong biofilm formation at pH 5.8 and 7.8, with cell counts above 4,000 and EPS matrix areas similar to those at 45°C. However, at pH 6.8, both viable cell numbers and EPS production decreased sharply, suggesting reduced biofilm viability at this specific pH-temperature combination. This indicates a shift in the optimal biofilm-forming window at elevated temperatures. Lastly, at 65°C, biofilm formation was minimal across all pH values. Cell viability drastically declined, with total viable cell counts falling below 15 and EPS matrix areas dropping to negligible levels. These results confirm that 65°C approaches the upper physiological limit for active biofilm development in B. licheniformis Tol1, although sparse biofilm fragments were still detectable at pH 7.8 (Figure 3). Overall, the results demonstrate that thermophilic B. licheniformis Tol1 forms highly structured and viable biofilms at moderate thermal conditions, particularly at 45–55°C and neutral to slightly alkaline pH. The presence of an extensive EPS matrix further confirms active biofilm development, which supports the bacterium’s resilience in thermal environments and strengthens its biotechnological relevance for thermophilic applications.

Considering the high apparent viscosity and the ropy behavior of B. licheniformis Tol1 in nutrient supplemented LB medium (Supplementary Table 1), the best-optimized carbon source was sucrose, and the nitrogen sources was NH₄Cl along with NaCl, resulting in an initial EPS production of 0.2 g L⁻¹ under pre-optimized condition. This was used for subsequent response surface methodology (RSM) analysis to optimize EPS production statistically.

RSM based on a central composite design technique has been successfully applied to optimize the possible combination of parameters (sucrose concentration %, NH₄Cl %, and NaCl %) that maximize EPS production by B. licheniformis Tol1 (a maximum of 2.11 g L⁻¹ at sucrose concentration 3%, NH₄Cl 3%, and NaCl 4%) as shown in Supplementary Table 2, the upper and lower limits of each parameter were specified according to initial laboratory study based on previous published research.

From the ANOVA (Table 2) for the quadratic optimization model, the model showed high significance with F-value 15.877 (p-value >F was <<0.0001); in addition, the model was evaluated by determination coefficient R² = % 93.4, which means the model does not explain %6.6 of total experiments, R² adjust, R² predict for EPS production were 0.875 and 0.757, respectively, which indicated good results with difference less than 0.2. The measurement of the signal-to-noise model was assessed by adequate precision and accepted with a value of more than 4 in this model (11.67). From ANOVA Table 2 for EPS production, it can be noticed most terms exhibited significant effect (p < 0.0001). This study evaluated the significance parameters by p-values, p-values more than 0.1, and revealed that the term was insignificant.

Furthermore, a regression optimization model was constructed, which is a relationship between independent parameters and dependent parameters (EPS production) developed after ANOVA and the regression coefficient determined. The CCD design matrix was fitted with a quadratic optimization model for EPS production in coded parameters as described in Equation 7. (7) EPS = 2.0673 + 0.1290 × A + 0.15752 × B − 0.0019 × C − 0.0618 × A × C − 0.228 × B × C − 0.4133 A 2 − 0.26 2B 2 − 0.3699 C 2

In addition to the analysis of variance (ANOVA) and correlation optimization model, regression analysis can be applied to evaluate the best-fit line among the experiments, which can be visualized by a normal plot of residual (Figure 4A). Moreover, as shown in Figure 4B, the observed experimental values against predicted output values for EPS production exhibited good aggregation between observed values and predictive values estimated by the quadratic optimization model. Hence, the developed optimization model can be generalized for EPS production by B. licheniformis Tol1 under the same condition.

Contour plots (2-D) were a graphical representation to visualize the current state of the regression optimization model and to deduce the relationship between parameters and response. Circular or elliptical curves indicate that corresponding parameters have not interacted, while curves in the shape of hyperbolic and plateau refer to interaction. This study explained the relation between independent parameters and EPS production by 2-D contour plots response surface while others were kept at zero level (coded value). In Figure 5a, the interaction effect of 3% sucrose as carbon source and 3% NH₄Cl as nitrogen source showed a maximum production of 1.97 g L⁻¹ EPS, whereas minimum EPS production of 0.36 g L⁻¹ was observed at minimum level for each carbon and nitrogen source concentrations. The maximum level of carbon and nitrogen sources also recorded a minimum EPS production of 0.8 g L⁻¹. The interaction effect of sucrose and mineral salt NaCl on EPS production is explained in Figure 5b. Briefly, 3% of sucrose and 3% of NaCl resulted in 1.97 g L⁻¹ of EPS production. Figure 5c revealed the interaction effect of NH₄Cl and NaCl on EPS production; a high amount of EPS production (1.97 g L⁻¹) was investigated when NH₄Cl concentration was at 3% and NaCl concentration was at 4%. Minimum EPS production at low concentrations of NH₄Cl and NaCl, but high concentrations of NaCl also revealed minimum EPS production.

As shown in Figure 6, an optimized artificial neural network (ANN) design with a 3–25-1 architecture was established. This model was employed to predict the production of EPS, where the network consists of three neurons in the input layer, a single hidden layer composed of 25 neurons, and one neuron in the output layer. The selection of this configuration was based on an optimization process aimed at maximizing the model’s accuracy in estimating EPS production, while ensuring a balance between learning capacity and generalization.

In this study, a feedforward neural network (FFNN) was employed as a modeling tool to evaluate the impact of various factors on EPS production and to conduct kinetic experiments. The FFNN enabled the analysis of complex patterns in the experimental data and accurately predicted the influence of the studied variables. The model’s performance was assessed using two key metrics: relative error (RE) and the coefficient of determination (R²). The results showed a maximum RE of −10.9%, indicating a low deviation from the experimental values. Furthermore, the R² value reached 0.9909 (Supplementary Figure 1), reflecting an excellent fit between the modeled and experimentally observed data. The results in Figure 7 demonstrate the model’s high accuracy, as the R² values are close to 1 across all phases of the process (training, validation, and testing). Specifically, during the testing phase, the model achieved an R² of 0.99626, indicating excellent predictive performance with previously unseen data. Similarly, in the training phase, the R² value was 0.99248, and when considering all phases combined, an R² of 0.99281 was obtained. These values confirm the robustness of the model in predicting EPS production and suggest that the ANN used possesses a suitable generalization capacity. The alignment of data points along the Y = T line in the regression plots further validates the model’s reliability, demonstrating a high level of agreement between predicted and experimental values.