The Newick tree file from the Bacillales phylogenetic TimeTree, which was generated from ribosomal RNA sequences (Marin et al., 2017), was visualized in TreeViewer 2.2.0 (Bianchini and Sanchez-Baracaldo, 2024). A list of putative OhyA sequences were obtained from tBLASTn by querying the clade I S. aureus OhyA amino acid sequence (NCBI-ProteinID: BAB41321) or clade II Listeria monocytogenes OhyA amino acid sequence (NCBI-ProteinID: NP_464009) against RefSeq Representative genomes in Bacillales/Caryophanales (taxid:1385). Default algorithm parameters (e.g., expect threshold 0.05, BLOSUM62 matrix, low complexity filer yes) were used for the tBLASTn search (Supplementary Tables S1, S2). Ten Paenibacillus species contained multiple hits in the same genome and both sequences were included in the analysis. Partial sequences, identified by a maximum alignment score <300, or sequences containing premature stop codons, identified by manual inspection, were removed to yield a final OhyA dataset of 145 sequences.

Species encoding putative ohyA genes were mapped onto the Bacillales TimeTree. Multiple sequence alignment was performed using MUSCLE within Mega 11.0.13 (Tamura et al., 2021) (Supplementary Datasets S1 and S2). Amino acid alignments were generated and refined using the following parameters: gap open −2.9, extend 0, hydrophobicity multiplier 1.2, max iterations 16, cluster method UPGMA, and min diag length 24. The phylogenetic trees of Bacillales OhyA amino acid sequences were constructed using IQ-TREE 2.1.4-beta. The best substitution model for accurate phylogenetic estimates was determined by ModelFinder (Kalyaanamoorthy et al., 2017) in IQTree. The OhyA amino acid tree was constructed using LG + I + G4 as the model of substitution. The outgroup Elizabethkingia meningoseptica OhyA (Engleder et al., 2015) was used for rooting. The Staphylococcus 16S rRNA tree was constructed using TPM2+F + I + G4 as the model of substitution, and the Staphylococcus ohyA nucleotide tree was constructed using GTR + F + I + G4 as the model of substitution. The IQ-TREE phylogenetic trees were constructed with 5,000 ultrafast bootstrap alignments (Minh et al., 2013), 5,000 maximum iterations, 1,000 replicates of the single branch test, auto-detect threads, and a minimum correlation coefficient of 0.99. Staphylococcus phylogenetic trees were rooted at the midpoint using TreeViewer 2.2.0 (Bianchini and Sanchez-Baracaldo, 2024), and then re-rooted and re-ordered for analysis using the beta version of phylo.io (https://beta.phylo.io) (Robinson et al., 2016) in compare mode. Sequence logo and residue composition probability for each catalytic loop clade were generated using WebLogo 3.7.12 (https://weblogo.threeplusone.com/create.cgi), using units of probability. The amino acid sequence percent identity matrix was generated by Clustal Omega (https://www.ebi.ac.uk/jdispatcher/msa/clustalo).

Protomer models for all putative Bacillales OhyA sequences were generated using AlphaFold 2.3.0 (Jumper et al., 2021). The structural alignment, similarity statistics (e.g., Q-score, P-score, Z-score), and R.M.S.D. were calculated by submitting the AlphaFold-generated PDB files to the PDBeFOLD server (Krissinel and Henrick, 2004) (https://www.ebi.ac.uk/msd-srv/ssm/). Although it was not used for this study, we expect that using the testing phase of AlphaFold three would yield similar results and not affect the conclusions of the study.

The 16S rRNA nucleotide sequences of Staphylococcus species were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG). 16S rRNA sequences for species that were not in KEGG were obtained from GenBank. Multiple sequence alignment was performed using MUSCLE within Mega 11.0.13. The multiple sequence alignment output was inspected to ensure the alignment was reasonable. The phylogenetic species tree was constructed from the 16S rRNA nucleotide sequences using IQ-TREE 2.1.4-beta, using GTR + F + G4 as the model of substitution. The Staphylococcus ohyA gene nucleotide sequences that were returned from the tBLASTn search described above were aligned and the phylogenetic ohyA gene tree was constructed using IQ-TREE 2.1.4-beta, using GTR + F + G4 as the model of substitution. The best substitution model for accurate phylogenetic estimates was determined by ModelFinder (Kalyaanamoorthy et al., 2017) in IQTree. Phylogenetic trees were constructed with 5,000 ultrafast bootstrap alignments (Minh et al., 2013), 5,000 maximum iterations, and a minimum correlation coefficient of 0.99. The ohyA gene context was determined by manual genome inspection and functional categorization. Gene length and orientation were used to scale graphics for accurate visualization.

The ohyA gene position was manually identified in each Staphylococcus species, and ∼5 k base pairs upstream and downstream were analyzed for open reading frames and predicted functions (Supplementary Table S3). Functional annotations were used to group open reading frames into seven broad functional categories: protease, transcriptional regulation, transport, nucleic acid metabolism, stress response, carbon and nitrogen metabolism, and hypothetical protein. Open reading frames in the transcriptional regulation group were searched against S. aureus subsp. aureus NCTC 8325 (taxid: 158879) using BLASTn.

The S. aureus USA300 JE2 strains used in this study are from the Nebraska Transposon Mutant Library (BEI Resources, Cat. No. NR-48501): wild type/parent strain (NR46S43); ΔsrrA knockout strain (NE1309, NE-47852, SAUSA300_1442); ΔsrrB knockout strain (NE588, NR-47131, SAUSA300_1441). Staphylococcus aureus strains were grown on mannitol salt agar overnight at 37°C, and then single colonies were grown in tryptic soy broth (TSB) shaking overnight at 200 rpm at 37°C. Cultures were diluted to OD₆₀₀ = 0.3 and treated for 2 h with 5 mM hydrogen peroxide (CVS Pharmacy, Inc, Cat. No. NDC 51316-268-16) or 5 mM diethylamine/nitric oxide complex sodium salt hydrate (MilliporeSigma, Cat. No. D184). Erythromycin (5 μg/mL) was used to maintain pressure on the knockout strains throughout the experiment. Cells were harvested (OD₆₀₀ ˜ 0.6) using a phenol:chloroform 1:1 solution, and RNA was extracted using the SV Total RNA Isolation System (Promega, Cat. No. Z3100) according to the manufacturer’s instructions.

The qRT-PCR was performed on a 7500 Fast Dx Real-Time PCR Detection System (Thermo Fisher Scientific) using TaqMan™ Fast Virus 1-Step Master Mix (Thermo Fisher Scientific, Cat. No. 4444436).

The primers sequences for ohyA are:

Probe: 5′-/56FAM/TGGTAACTTTGCAGAAACAGAGCGAGA/36TAMp/-3′

Forward: 5′-CATGGCAGTACGAACCGAATA-3′

Reverse: 5′-CTTTAGTCGTCCCGCATCAA-3′

The primers sequences for the gyrA calibrator are:

Probe: 5′/56FAM/CGGTATCACTGATTTACGTGATGAAAC/36-TAMp/3′

Forward: 5′-CCTTAGCGACATCAATAACGACACGC-3′

Reverse: 5′-AATTGCAGAGCTCGTTCGTGACAAG-3′

Cycling conditions were reverse transcription for 5 min at 50°C and pre-denaturation for 20 s at 95°C, followed by 40 cycles of denaturation for 15 s at 95°C and annealing for 60 s at 60°C. Fluorescent signals were detected at the end of each cycle. The relative fold gene expression was calculated using the 2^−ΔΔCT method (Livak and Schmittgen, 2001).

Quantification and data analysis (i.e., calculation of mean, standard deviation, significance testing) were performed using PRISM 9.5.1 (GraphPad Software, https://www.graphpad.com/scientific-software/prism/).

To investigate the evolutionary dynamics of OhyA, we performed a comprehensive evolutionary analysis of putative ohyA genes from Bacillales species. Using tBLASTn with S. aureus OhyA as the query sequence, we built an ohyA nucleotide sequence database, identifying 145 sequences with E-values ranging from 0.0 to 2.00E-91 and identity levels between 35.02% and 97.53% (Supplementary Tables S1, S2). The ohyA nucleotide phylogenetic tree, computed using a nucleotide model, revealed two distinct OhyA clades, which we named clade I and clade II; however, species from the same genera were not clustering together. Given the broad sequence diversity, including those with less than 40% identity to the query, we constructed a phylogenetic tree based on amino acid sequences (Figure 3). The amino acid tree, supported by similar bootstrap values, reproduced the two OhyA clades and sequences from species of the same genus clustered together. These data suggest that ohyA nucleotide variation is largely due to synonymous mutations that do not alter the encoded protein. Clade I OhyAs are present in Staphylococcus, Lysinibacillus, and 10 out of the 84 Paenibacillus species analyzed. Clade II OhyAs, on the other hand, are found in 10 out of 11 Listeria species and 74 out of 84 Paenibacillus species analyzed. Notably, the only Listeria species encoding a clade I OhyA is Listeria grayi, an avirulent species that has been proposed for reclassification under the genus Murraya (Stuart and Welshimer, 1974), through it is still categorized as Listeria for the time being (Rocourt et al., 1987). The phylogenetic distribution of OhyA reveals that it is found not only in pathogenic organisms but also in environmental microbes that inhabit diverse terrestrial ecosystems. This suggests that OhyA may offer niche-specific metabolic advantages beyond functioning as a defense mechanism against the immune system.

We found that OhyA proteins from the two clades diverge in key features, including the amino acid sequence of the catalytic loop (Figures 4A, B), protein length (Figure 4C), and sequence conservation (Figure 4D). Clade I enzymes exhibit a narrow length range of 590.2 ± 0.8 amino acids, while clade II enzymes display greater variability, with an average length of 536.6 ± 18.4 amino acids. Within each clade, amino acid sequence similarity is 73.1% ± 9.4%, whereas similarity between the two clades is lower at 41.7% ± 7.5%. Notably, clade II enzymes lack a conserved glutamate residue, essential for catalysis in clade I OhyA (RGGREM) as described by Radka et al. (2021b).

To further explore structural diversity within the OhyA protein family, we used AlphaFold (Jumper et al., 2021) to predict the three-dimensional structures of Bacillales OhyA proteins. Pairwise structural comparisons, conducted using PDBeFOLD (Krissinel and Henrick, 2004), showed that the predicted structures cluster into two groups, corresponding to the clades identified in the phylogenetic analysis (Figure 5A). Despite low amino acid sequence similarity, the core catalytic domain structure is highly conserved across both clades. This is evidenced by a Root Mean Square Deviation (RMSD) of less than 2 Å in pairwise comparisons of Cα atoms in matched residues, with clade I enzymes having 569.6 ± 12.4 matched atoms and clade II enzymes having 486.6 ± 3.1 matched atoms.

Clade I enzymes typically feature a three-domain architecture consisting of a fatty acid lobe, FAD lobe, and a membrane binding domain. In contrast, clade II enzymes lack the characteristic membrane binding domain at the carboxy terminus, through some clade II enzymes feature a distinct alpha helix at the amino terminus (Figure 5B). This structural divergence suggests that clade II enzymes may have evolved alternative mechanisms for interacting with membrane-bound substrates, demonstrating the evolutionary flexibility in solving similar biochemical challenges. These structural and sequence variations indicate the presence of at least two distinct biochemical mechanisms for unsaturated fatty acid hydration within the OhyA family.

Evolution and gene regulation in the Staphylococcus genus are well studied (Le and Otto, 2015; Coates-Brown et al., 2018; Raghuram et al., 2022), and Staphylococcus species exclusively encode clade I OhyA. Therefore, a deeper investigation into Staphylococcal OhyA could provide generalizable evolutionary and regulatory insights applicable to other, less-studied Bacillales genera. To investigate the molecular evolution of ohyA within the Staphylococcus genus, we employed a comparative genomics approach using phylogenetic and experimental bacteriology. We constructed a 16S rRNA phylogenetic tree for Staphylococcus species, revealing that 18/59 species encode genomic ohyA orthologs, with one species encoding an ohyA ortholog on a plasmid [Staphylococcus condiment plasmid pStO 2014-01 (Gabrielsen et al., 2017)] (Figure 6). BLAST E-values indicate all of the putative Staphylococcal ohyA genes are statistically homologous to the S. aureus ohyA gene (Supplementary Table S1), and the topology of the phylogenetic tree of Staphylococcus species that encode ohyA indicate a shared evolutionary origin. Further, there is good agreement between the Staphylococcal speciation nodes in 16S rRNA species tree and the ohyA gene tree (Figure 6; Supplementary Table S1).

Incongruence of phylogenetic trees or G + C content are traditional methods to identify candidate genes of horizontal gene transfer. The Horizontal Gene Transfer (HGT-DB) Database predicts the genomes of S. aureus strains contain 4.6%–5.8% horizontally transferred genes based on extraneous G + C content (Garcia-Vallve et al., 2003). The S. aureus strain N315 is a clinically relevant, methicillin-resistant reference strain that is commonly used in the laboratory setting Kuroda et al. (2001); (Gerlach et al., 2018). The representative genome from S. aureus strain N315 is 33.2% ± 0.3% G + C across 2,594 open reading frames and its ohyA is 35.8% G + C, which are congruent to each other. Although HGT-DB predicts S. aureus strain N315 encodes 105 horizontally transferred genes, ohyA is not one of them. The HGTree database identifies candidates for horizontally transferred genes by tree reconciliation between gene trees and 16S rRNA reference species trees (Jeong et al., 2016). HGTree predicts S. aureus strain N315 encodes 542 horizontal transfer events, ohyA is not one of them and is consistent with our observation of similar speciation between the Staphylococcus species tree and the ohyA gene tree. Thus, both methods agree that ohyA is not a Staphylococcal horizontally transferred gene.

The high amino acid sequence similarity among ohyA genes within Staphylococcus species, the congruence between ohyA G + C content and overall genomic G + C content, and the agreement between the ohyA gene tree and the Staphylococcus species tree all support an evolutionary scenario in which an ancestral Staphylococcus species encoded ohyA. During evolutionary processes, many contemporary Staphylococcus species lost the ohyA gene, leading to the current distribution of ohyA-encoding species (Figure 6). This pattern of gene loss is consistent with the concept of genetic streamlining, where genes that are not essential for survival in a given environment are selectively lost to reduce metabolic burden.

The synteny around ohyA provides further insights into its regulation and potential co-evolution with other genes. Alterations in the mechanisms of transcriptional regulation can lead to expression divergence, influencing the evolutionary path of a protein (Gu et al., 2002; Makova and Li, 2003). We analyzed the genetic loci flanking ohyA in Staphylococcus species to assess synteny and identify regulatory elements that may influence ohyA expression (Figure 7; Supplementary Table S3). In at least 16 of the 21 species analyzed, one or more transcriptional regulator genes were found in proximity to ohyA, suggesting a common regulatory adaptation to specific environmental conditions.

In S. aureus, the regulation of ohyA expression is more complex. While a transcriptional regulator gene is located near ohyA, it does not appear to modulate ohyA expression directly. The S. aureus transcriptional regulator that flanks ohyA is norG (SA0104/SAUSA300_0110), a regulator for several efflux pumps (Truong-Bolduc et al., 2011). A microarray screen of S. aureus genes regulated by norG did not find ohyA (Truong-Bolduc et al., 2011). Regulation by a tetracycline repressor family fatty acid efflux pump transcriptional regulator FarR (SA2340/SAUSA300_2490 and not encoded among the ohyA flanking loci) has also been considered, but the regulatory stimulus (unsaturated fatty acid) does not robustly impact ohyA expression (Subramanian et al., 2019). Instead, ohyA expression is controlled, at least in part, by the SrrAB two-component system, which is conserved across Staphylococcus species. A microarray screen of S. aureus genes found 8.3-fold higher levels of ohyA RNA transcripts in wildtype cells compared to a double knockout ΔsrrAB strain under nitrosative stress conditions (Kinkel et al., 2013). The SrrAB system is known to regulate genes involved in anaerobic metabolism and oxidative stress response (Kinkel et al., 2013; Jenul and Horswill, 2019; Qiu et al., 2021), linking ohyA expression to broader stress response pathways.

We used qRT-PCR to analyze ohyA gene expression in ΔsrrA or ΔsrrB single knockout strains to examine the impact of each component under oxidative (hydrogen peroxide, H₂O₂) and nitrosative (nitric oxide, NO) stress conditions (Figure 8). In untreated cells or under nitrosative stress, the ΔsrrB strain showed 7.3-10.5-fold lower ohyA expression than the wildtype or ΔsrrA strains, but was only significantly different from the expression level in the ΔsrrA strain. When the cells were under oxidative stress, ohyA expression in the ΔsrrB strain was significantly 4.1-5.0-fold lower than wildtype and ΔsrrA levels. This regulatory network may reflect the importance of ohyA in helping S. aureus adapt to hostile environments, such as those encountered during infection, where oxidative stress and limited oxygen availability are common challenges.

The presence of conserved regulatory elements across Staphylococcus species suggests that ohyA regulation has been maintained despite the loss of the gene in many species. This conservation implies that ohyA plays a critical role in the species that retain it, possibly conferring a selective advantage in specific ecological niches. Understanding the regulatory mechanisms that control ohyA expression can provide valuable insights into how bacteria adapt to environmental stresses and how gene regulation evolves in response to changing selective pressures.