The Churince water system is located on the west side of CCB, and is formed by a spring, an intermediate lagoon, and a desiccation lagoon. The water and sediment samples described in this work were obtained from the intermediate lagoon. Physical and chemical measurements revealed homogeneous conditions along the sampling points, including pH (7.8–8), salinity (∼1.7–2.1 ppt), and daily temperature fluctuation (from 20 to 35°C; Pérez-Gutiérrez et al., 2013).

The isolates in the collection were obtained in October 2007 from water column and sediment samples of the Intermediate Lagoon. Eleven sites, each about 30 m apart within the Churince water system, were chosen for small scale community sampling. At any point water samples were taken carefully as to not disturb the sediment (Column samples). Then approximately 50 μL of superficial (first millimeters TOP sample) sediment were obtained, and then another approximately 50 μL sample 2 cm below this (DEEP sample) (Figure 1). The samples were taken to an improvised in field laboratory, and a few hours later PBS buffer was added to the sample, and it was heat treated to select for thermo-resistant bacteria and plated out on Marine Medium (MM; Cerritos et al., 2008), and incubated at 37°C for 2 days. Different colony morphotypes were selected by size, shape, and color. Isolates were purified by subculture on the same medium to ensure that the culture was axenic; all isolates were stored at -80°C in MM with 15% (w/v) glycerol. In this work we show results for 467 thermo-resistant isolates that were obtained from eleven sampling points (Figure 1).

DNA from isolates was purified using the QuickGene DNA tissue Kit S (Fujifilm Corporation, Minato-Ku, Japan) following the manufacturer’s instructions. We amplified the 16S rRNA gene by PCR using previously published primer sequences 27F and 1492R (Lane, 1991). This included regions V1 to V3 (∼275 bp of the 5′ end region), considered to be the most informative for the Bacillus spp. (Goto et al., 2000). A 700 bp segment of the gltX gene (glutamyl-tRNA synthetase) was amplified using the primers gltX-for (5′ CGYGGBGADGAYCAYATYT 3′) and gltX-rev (5′ CRATTTCMGCDCCRWARCT 3′) and PCR-amplified by 30 cycles of 94°C for 30 s, 47°C for 1 min, 72°C for 2 min, and then a final elongation step at 72°C for 8 min with a thermocycler Palm-Cycler (Corbett Research). The PCR products were sequenced using the Sanger method (Sanger et al., 1977) at Cinvestav-Langebio (Irapuato, México). The quality of the 16S rRNA gene sequences was assessed using Phred (Ewing et al., 1998), and sequences of at least 500 bp were used for the analysis. The gltX gene sequences were assembled with CAP 3 (Huang and Madan, 1999). Gene accession numbers and lists of all 467 strains can be found in Supplementary Table S1. Maximum likelihood (Guindon et al., 2010) was used to reconstruct a phylogenetic history of all thermo-resistant bacteria through the 16S rRNA gene. A second phylogenetic reconstruction was done using the MrBayes program with sequences from a subsample of 138 Bacillus spp. for which phenotypic traits were obtained. Sequences were aligned with Muscle Software (Edgar, 2004). 16S rRNA Maximum-Likelihood tree was reconstructed with MEGA6 Software (Tamura et al., 2013). Following Bayesian Information Criteria (BIC) the Tamura 3-parameter + Gamma model was used to estimate distances between sequences on the tree generated with 100 bootstrap replicates. The 16S rRNA Bayesian tree was reconstructed with MrBayes (Ronquist et al., 2012). We used 5,000,000 Markov chain Monte Carlo generations and the GTR + G model of evolution. Average standard deviation of split frequencies dropped below 0.01. The gltX sequences of these 138 Bacillus spp. were also obtained to increase the resolution of the particular groups. Data from the group had shown that this gene allows more resolution than gyrA and rpoD. Phylogenetic reconstruction was done for these sequences using MrBayes and the sequences were grouped with CLANS (Huelsenbeck and Ronquist, 2001; Frickey and Lupas, 2004).

A total of 467 16S rRNA gene sequences were analyzed (TOP 240, DEEP 164 and COLUMN 63) using Mothur to assign to operational taxonomic units (OTUs) having 97% similarity thresholds (Schloss et al., 2009). The ribosomal database project (RDP) classifier was used for taxonomic assignment, with a minimum confidence of 0.8 to record an assignment (Cole et al., 2007). To identify ecologically distinct populations we examined the phylogenetic data for the group using AdaptML software, since it classifies a phylogeny based on genetic and ecological data and identifies populations as groups of related strains with distinct environmental distributions (Hunt et al., 2008). Tree drawings were generated with iTOL (Letunic and Bork, 2011). The net relatedness index (NRI) measures the degree of phylogenetic clustering of taxa across a phylogenetic tree in a given sample relative to the regional pool of taxa and was calculated with the Picante package (Webb, 2000; Kembel et al., 2010).

To evaluate some phenotypic traits of the isolates, including carbon source utilization, growth kinetics, auxotrophies, motility assays, and sporulation frequency, we used a modified LPDM medium (Muller et al., 1997) that contained: Tris (pH 8.0) [50 mM], NH₄NO₃ [3.3 mM], K₂HPO₄ [0.065 mM], MgSO₄ [3.5 mM], sodium citrate dehydrate [6.8 mM], MnCl₂ [1 mM], ZnCl₂ [0.01 mM], NaCl [0.17 M], FeCl₃ [49.9 μM], CaCl₂ [3.60 mM], KCl [0.67 mM], glucose [100 mM], vitamin B complex [100 μL/L], biotin [0.1 mg/mL], nicotinic acid [0.1 mg/mL], peptone [0.05%] and 20 amino acids (Harwood and Cutting, 1990).

To determine preferences for carbon source utilization we used five carbon sources: xylose, raffinose, sorbitol, trehalose, and glucose, all at a concentration of 100 mM in the LPDM medium. The assays were performed using microtiter plates (96 wells) and read on a BioTek μQuant Microplate Spectrophotometer. To estimate the maximum growth rate (μmax) and saturation density, growth kinetics were monitored under two concentrations of glucose (5 and 50 mM). The plates were incubated at 30°C, with constant vibration (100 Hz) and 90% humidity. Growth was monitored every 30 min for 24 h at OD₆₀₀ using a microtiter plate reader (Tecan Infinite^® M1000) assisted by a TECAN EVO robot (Tecan Infinite^® M1000).

To evaluate whether the different strains could grow in the absence of a given amino acid, peptone or yeast extract, cultures were first grown in liquid MM, and then passed twice to LPDM lacking amino acids, yeast extract or peptone. Strains that grew under these conditions were considered to be prototrophs. Those strains that required addition of amino acids, peptone or yeast extract were further evaluated by transferring them to LPDM that contained either all 20 amino acids, and/or peptone or yeast extract. Those strains that could grow without peptone and yeast extract were then evaluated in LPDM medium containing all but one of the 20 amino acids. Growth was monitored at 24 h (Supplementary Table S2).

Motility assays were performed according to the protocol described in Hamze et al. (2009) for swarming, and Kearns and Losick (2003) for swimming. Briefly, Petri dishes were prepared with either 0.6% (swarming) or 0.3% (swimming) MM agar and allowed to dry at room temperature. Each Petri dish was inoculated at the center and the culture was allowed to grow for 24 h at 37°C. The ability to form biofilms was determined in a microtiter plate (24 wells) in MM liquid after static incubation for 48 h at 37°C.

We compared the phylogenetic signal for the different phenotypic traits under evaluation. We determined the K value (Blomberg et al., 2003), a parameter that describes whether a pattern that arises when closely related taxa are more ecologically similar to each other than to distantly related taxa, and used this value to describe a tendency (pattern) for evolutionarily related organisms to resemble each other (K value is used for continuous traits, such as growth). K = 1 indicates a Brownian expectation, in which trait changes along each branch are random, and K > 1 indicates traits that are more conserved than the Brownian expectation. K = 0 indicates absence of a phylogenetic signal, and values between 0 and 1 indicate a significant phylogenetic signal but less conserved than under a Brownian motion model of evolution. We also carried out a Fritz and Purvis D-test applied to traits with discrete values. A D < 0 suggests a highly clustered trait, D ∼ 0 indicates a Brownian motion mode of evolution (BM), D = 0 suggests a random mode of evolution and D > 1 suggests phylogenetic overdispersion (Fritz and Purvis, 2010). D statistics can also be expressed in the same terms as the K statistics: -D + 1 = 0 value does not show a significant signal, -D + 1 > 0 is more conserved than expected by chance, 0 < -D + 1 < 1 indicates that is less conserved than expected under Brownian Motion model, -D + 1 = 1 is as conserved as expected under Brownian Motion model and – D + 1 > 1 suggest that is more conserved than expected under Brownian Motion model (Goberna and Verdú, 2016). P_random and P_Brownian convey the probability that the presence of a trait is the result of either a random or a Brownian motion mode of evolution, respectively. These tests help to identify the phylogenetic conservatism. The K and D values were calculated using the software packages Picante and Caper, respectively (Kembel et al., 2010; Fritz and Purvis, 2010), and both analyses were conducted in the R software environment (R Development Core Team, 2009). We used the consenTRAIT algorithm, that allows determining the cluster size where most organisms share a trait or, in other terms, the phylogenetic depth at which traits are conserved (Martiny et al., 2013).

The statistical analysis one-factor ANOVA was used to analyze the most numerous groups of Bacillus such as B. subtilis/licheniformis. B. pumillus. B. cereus, group 2, and the aquatic members. The traits subjected to this analysis were: carbon source utilization, maximum growth rate, and saturation density. The statistical analysis was carried out in the R programming environment (R Development Core Team, 2009).

We first characterized the diversity of culturable thermo-resistant bacteria from 11 sites within the Churince intermediate lagoon at CCB. We recovered 467 thermo-resistant isolates and analyzed partial 16S rRNA sequences that resulted in 77 OTUs (using Mothur with 97% similarity thresholds; Figure 1; Schloss et al., 2009). Among the 467 isolates, 428 (91.6%) were classified as Bacillus spp., 26 were classified as Koccuria spp., and 16 as Staphylococcus spp. (Figure 1). Sequences were assigned to a habitat category: COLUMN (water column, 63 members), TOP (surface sediment, 240 members), and DEEP (sub-surface sediment, 164 members). In order to know the phylogenetic structure among the different environments sampled, we determined if there was phylogenetic clustering among the 467 isolates using NRI, an index of community phylogenetic relatedness. We obtained positive NRI values for the sediment samples (TOP with 0.58 and DEEP with 1.83) indicative of phylogenetic clustering. For COLUMN sequences, however, the NRI values were negative (-2.28), suggesting a dispersed phylogenetic pattern of the thermo-resistant isolates from water column in comparison with the sediment sample. OTUs present in the COLUMN sample were more genetically dispersed, while those from sediment are more clustered, particularly those from the DEEP sample. We examined whether there was spatial resource partitioning among thermo-resistant isolates from COLUMN, TOP and DEEP samples. To do this, we used AdaptML which uses genetic and ecological similarity and classifies populations forming a common habitat into “projected habitats.” The populations were not classified and instead the analysis identified a single projected habitat, suggesting that there is no niche differentiation (data not shown).

Clustering based on 16S rRNA gene revealed that several culturable Bacillus spp. lineages co-occur in the Churince system. To determine the phenotypic coherence of the different taxonomic lineages, we carried out a phenotypic trait analysis in a subsample of 141 isolates from 5 sampling sites. All strains in the subsample belonged to the Bacillus genus and represented the dominant taxa. The chosen isolates were distributed in 13 clades (according to CLANS analysis, Frickey and Lupas, 2004), of which 6 included 80% of all members. The selected isolates were phylogenetically closer to B. subtilis/B. licheniformis. B. pumillus. B. marisflavi. B. aquimaris. B. horikoshii, and B. cereus. According to a previous classification of Bacillus species based on evolutionary and functional relationships, we classified the isolates into three different groups: (1) soil members (B. subtilis/B. licheniformis, and B. pumillus), (2) aquatic members (B. marisflavi, B. aquimaris, and B. horikoshii), and (3) B. cereus-related species (Alcaraz et al., 2010).

The Bacillus spp. exhibited differences in the preference for different carbon sources. We addressed whether this was an overdispersed trait or if there was a phylogenetic signal associated to the phenotype. We chose the Purvis and Fritz’s D value that permits evaluating discrete traits and the significance of their genetic grouping. All isolates could grow on glucose and all clades had some members capable of using xylose, raffinose, sorbitol, and trehalose for growth (Figure 2). The soil member group exhibited the broadest ability to use different carbon sources, followed by the aquatic group, and finally the B. cereus group. All Bacillus spp. were able to consume glucose, so this is a highly conserved trait. On the other hand, the Fritz and Purvis’ D value for trehalose, raffinose and sorbitol show a significant phylogenetic signal (D = 0.77, 0.76 and 0.81, respectively) and differs significantly from the random and Brownian motion model expectations (P_random < 0.05 and P_Brownian < 0.05). Finally, xylose does not show a significant phylogenetic signal (D = 0.85), since it does not differ significantly from the random model expectations.

Since the Churince water system is an oligotrophic environment, we also evaluated the ability of Bacillus spp. to respond to changes in nutrient availability as well as the saturation density (maximum growth density) the community members could reach. Isolates from mainly the B. subtilis/B. licheniformis group responded to increased glucose concentrations, but a significant growth increase was also observed in strains that were scattered throughout the phylogeny (Supplementary Figures S1A,B). All groups exhibited higher saturation density at higher glucose concentrations, but there were no statistically significant differences in saturation density or maximal growth rate in either of the two glucose concentrations across the taxonomic groups (ANOVA analysis did not give a significant value for Vmax or saturation density at 5 and 50 mM glucose; Supplementary Figure S1A). For this data, the Blomberg K value from Vmax was 1.27 at 5 mM glucose, and 3.75 at 50 mM glucose, which indicates that the observed value is more conserved than that expected in a Brownian model. The Blomberg K values for saturation density were 0.83 and 0.56 at 5 mM and 50 mM glucose, respectively, suggesting that the traits are conserved (all P-values observed were significant Vs. random variance of 0.001; Supplementary Table S3B). These data suggest that Vmax growth and saturation density values are conserved traits for the different taxonomic groups.

Co-occurrence of bacteria in a community could allow metabolic sharing or simply provide ample opportunity for scavenging and, either way, there would be no need for all members to be self-sufficient. We, therefore, evaluated metabolic self-sufficiency by analyzing auxotrophy for the different amino acids among members of the Bacillus spp. taxa from the sediment communities. Close to 60% of the evaluated Bacillus spp. isolates had an auxotrophy. Some members exhibited a strong dependency on peptone or yeast extract for growth even in the presence of all 20 amino acids, and this too was scored as an auxotrophic trait (44 strains required these supplements). These strains may obtain their amino acids in the form of peptides. For example, growth of aquatic Bacillus spp. (B. horikoshii. B. aquimaris, and B. marisflavi) exhibits extensive dependence on added amino acids or peptone (Figure 3). The same is true for B. cereus, for which prototrophy was rarely observed. The analysis of phylogenetic signal for this trait was consistent when measured with K and D statistics. It is compatible with a Brownian model of evolution using Bloomberg’s analysis (K of 1.28 with a PIC_random and a P_Brownian of 0.0009) implying a non-random phylogenetic signal. The D statistic was 0.27 with a P_random = 0 and P_Brownian = 0.12, a pattern that differs significantly from the random model, but the Brownian Motion model of evolution can not be rejected (Orme, 2013). Analysis with calculated a depth for 16S of τ_D 0.0095 16S rRNA distance, deeper than those computed for individual sugars (τ_D range from 0.001 to 0.003 16S rRNA distance), suggesting deeper lineage conservation (Supplementary Table S3C).

We determined the frequency of swimming and swarming motility and the ability to form biofilm in this Bacillus spp. collection. Our results showed microdiversity for these traits, but more so in certain taxonomic groups. Although swimming seems to be a conserved trait among all clades, some members in the aquatic group seem to have lost this feature. For example, 9 out of 23 B. horikoshii strains could not swim (Figure 3). Notably, swarming was exhibited by most members of the B. subtilis/licheniformis and B. pumilus clade, but 75% of B. cereus members scored negatively for this feature (22 of a total of 34). Swarming was absent in all members of the B. horikoshii clade and was seen in only one member of the B. marisflavi and B. aquimaris clade (Figure 3).

For biofilm formation, heterogeneity was most noticeable among the aquatic Bacillus spp., but this trait was almost always present in members of the B. subtilis/licheniformis-pumilus or B. cereus groups (Figures 3 and 4). Moreover, almost all soil group members could form biofilms (only 1 of 36 members of the B. subtilis/licheniformis. B. pumilus group did not, and 3 out of 34 of the B. cereus members did not display this feature) but less than half of the B. horikoshii clade formed biofilms.

Growth kinetics analyses for the different isolates revealed three different but consistent patterns for the assayed Bacillus strains: Pattern 1, a characteristic smooth growth curve was present throughout the growth period (47% of the isolates); Pattern 2, a smooth logarithmic growth phase was followed by aggregation at the end of growth period (resulting in fluctuating readings during the stationary phase in 50% of the isolates); and Pattern 3, readings fluctuated throughout growth (this pattern occurred for 3% of the evaluated strains) (Figure 4). These patterns indicate that there were differences in the growth strategies of the different strains, wherein some could display a planktonic behavior (Pattern 1) with no aggregation at the end of growth (although our evaluation for biofilm formation showed that many strains had this ability) and some had aggregation that was triggered immediately at the end of the exponential growth phase (Pattern 2).

Analysis of phylogenetic signal showed that the distribution of biofilm-forming ability across the taxonomic groups of Bacillus spp. from these sediment communities exhibited only a moderate level of phylogenetic clustering (D = 0.51), and differs from a random and Brownian model (P_Brownian 0.02 and P_Random 0). K is 0.8 and is significantly different from random (Supplementary Tables S3A,B). We conclude that for this trait there is a departure from a random model and a Brownian model of evolution. Regarding swimming, the D value obtained suggested that this phylogenetic exhibited a random mode of evolution (D = 0.90; P_Random = 0.19, P_Brownian = 0). A Blomberg analysis, in contrast, resulted in a K = 0.78 with a P of 0.0009 for a non-random mode of evolution. For swarming and prototrophy the D values were 0.21 and 0.28, respectively. The statistic indicates that both swarming and prototrophy traits differed from a random pattern of evolution and resembled a Brownian motion (both have a P_Random = 0; swarming P_Brownian = 0.181, prototrophy P_Brownian = 0.106).

We also computed the phylogenetic depth at which the observed traits were conserved using consenTRAIT (Martiny et al., 2013). It calculates the mean depth of clades containing above 90% of members sharing a trait. Cluster size is therefore an important data to take into account for this analysis. For all traits evaluated the 16S rRNA distance was τ_D = 0.00073 to 0.0095 16S rRNA distance. As a reference, the deepest value, τ_D 0.07416S rRNA distance, would be that for glucose (used by all the strains in the study). Lower distances (τ_D) were obtained for substrate utilization than for social traits (except for swarming and xylose). The largest size clusters computed were for prototrophy, swimming, and biofilm formation (size 38), compared to swarming (size 23), but the mean cluster size was higher for the social traits (biofilm, swimming, prototrophy). Social traits also had more depth signal (for biofilm, swimming, prototrophy τ_D = 0.005, 0.0056, 0.0095, 16S rRNA distance respectively) compared to the genetic distance for traits involved in substrate utilization (i.e., trehalose τ_D = 0.0011 and raffinose τ_D = 0.0015 16S rRNA distance) (Supplementary Table S3C).

Overall, the evolutionary models that describe the distribution of social traits in the phylogeny (except for swimming) suggest a Brownian model of evolution and a deeper phylogenetic signal than that observed for substrate utilization.

Colony morphology, size, and pigmentation were observed to be very conserved and are distinctive features of some clades. The colony appearance in our sample is a good proxy for the phylogenetic membership of a strain. Although colony color may exhibit differences in intensity and can change with growth stage and medium, color is generally consistent for members of the same clade. Colony morphology and color probably involve a large number of genes in cellular processes. These genes are differentially and spatially expressed to form a highly organized and structured community. The fact that morphology conservation occurs suggests that it is not neutral but instead arises from pressure to maintain a given developmental pattern that may be relevant for survival. However, this would have to be further studied.

Finally, the 16S rRNA gene phylogeny was coherent with respect to the phenotypes of the clade members. We added one more gene, gltX, to increase phylogeny resolution. As expected, this gene maintained the deeply rooted clades in the same order as the 16S, although some of the tips of the tree changed position within their own basal clade. The so-called soil Bacillus (B. subtilis/licheniformis and B. pumilus) as well as B. cereus members remained as a group. Meanwhile several of the aquatic groups changed position at the tips of the tree (Figure 5) suggesting that gene transfer is more common in the aquatic clades than in soil clades’ members.