1 Introduction

Lipids are ubiquitous and diverse group of amphiphathic compounds that are critical for fundamental biological processes, including the formation of cell membranes, protection, and cargo delivery; energy storage; and cell signaling pathways [1]. Chemical modifications to lipids, including changes in fatty acid tails (saturated and unsaturated) [2], modifications of head groups [3], and alterations to lipid charge [4], impact these interactions, as well as the physicochemical properties of membranes. Commonly, bacteria synthesize lipids that are negatively charged or neutral [5]. Some bacteria modify the negatively charged lipids with amino acids to make them positively charged to confer resistance to cationic antimicrobial peptides and cationic antibiotics [6, 7]. The bacterial domain is a rich source of natural lipids that are likely co-evolved for specific interactions with host tissues and other cells, yet the full extent of this chemical diversity remains unexplored. Finding bacteria that possess unique lipids can be useful for biotechnology. For example, natural lipids found in bacteria can be used for drug delivery. A study by [8] used low cytotoxic outer membrane vesicles (OMV) that contained a modified lipopolysaccharide (LPS) to deliver siRNA to targeted cancer cells. This delivery mechanism had success and demonstrates the possibility of using naturally occurring bacterial lipids for future drug delivery processes to eukaryotic cells.

One enzyme responsible for modifying lipids in bacteria is the multiple peptide resistance factor (MprF). MprF modifies lipids through the transfer of amino acids from charged tRNAs to the head group of the anionic membrane phospholipid, phosphatidylglycerol (PG) [4]. The MprF protein is comprised of two domains: an N-terminal flippase domain responsible for moving aminoacylated lipids from the inner membrane leaflet to the outer leaflet and a C-terminal synthase domain responsible for the aminoacyl transfer from tRNA to lipids [9, 10]. Recently, in Streptococcus agalactiae (Group B Streptococcus, GBS) it was found that MprF can modify two different substrates with lysine (Lys) – the glycolipid glucosyl-diacylglycerol (Glc-DAG) and the phospholipid, PG [11], generating Lys-Glc-DAG, as well as Lys-PG [11, 12] (Figure 1). An S. agalactiae mprF deletion mutant no longer synthesizes Lys-Glc-DAG or Lys-PG and expression of S. agalactiae mprF in the heterologous host Streptococcus mitis conferred Lys-Glc-DAG and Lys-PG synthesis to S. mitis [11]. This was the first time MprF was demonstrated to add lysine onto a neutral glycolipid (Glc-DAG).

Figure 1:

The synthesis of Lys-Glc-DAG by S. agalactiae MprF illustrated a unique lipid substrate and product by the enzyme that had not been characterized before, highlighting the possibility for further exploration into unknown lipids and lipid substrates by MprF of different bacterial species. We sought to investigate the molecular determinants of enzyme specificity and identify other bacterial species that may have this novel specificity. This investigation led us from the standard methodology of BLASTp queries [13], which are based on single residue amino acid substitution frequencies, to a more modern statistical method called a Restricted Boltzmann Machine (RBM) [14], which captures global statistical patterns within natural sequence data. In this work, we develop a general strategy for combining this high-powered statistical model with lipidomic analysis to discover MprF variants with unique enzymatic activity. This process led us to the discovery of an uncharacterized lipid substrate for MprF, diglucosyl-diacylglycerol (Glc₂-DAG), which is present in several strains identified by the parameters of the statistical model. These models also provide interpretable and actionable information for use in further enzyme characterization, and we show that this application is a useful companion tool for experimental work.

2 Results

2.1 Simple pairwise sequence statistics identify streptococcal MprF enzymes that synthesize Lys-Glc-DAG

Utilizing a method based on the amino acid sequence of S. agalactiae MprF, we identified four streptococcal MprFs with high sequence identity to S. agalactiae COH1 MprF (WP 000733236.1). A BLASTp analysis found that Streptococcus sobrinus ATCC 27352 MprF (WP 019790557.1) had 61.4% amino acid identity to the S. agalactiae MprF, Streptococcus salivarius ATCC 7073 MprF (WP 002888893.1) had 43.09% amino acid identity, Streptococcus ferus MprF (WP 018030543.1) had 61.88% amino acid identity, and Streptococcus downei (WP 002997695.1) had 61.03% amino acid identity (Table 1). Plasmids were designed to express the S. sobrinus mprF (pSobrinus), S. salivarius mprF (pSalivarius), S. ferus mprF (pFerus), and S. downei mprF (pDownei) genes. The plasmids were transformed into the expression host S. mitis NCTC12261 (SM61). The empty vector (pABG5ΔphoZ, referred to as pABG5 here) was also transformed (Figure 2a). S. mitis does not natively encode mprF but synthesizes lipids (Glc-DAG; PG) that are substrates for S. agalactiae MprF [15, 16, 11], making it an appropriate heterologous host for expressing streptococcal MprFs. Previously, a plasmid expressing S. agalactiae mprF (pGBSMprF) [11] was generated and transformed into S. mitis. Lys-Glc-DAG and Lys-PG presence in SM61 was only possible with the expression of pGBSMprF (Figure2b) [11].

Table 1:

Figure 2:

S. mitis lipids of strains expressing the plasmids pSobrinus (Table 1), pSalivarius (Figure 2c), pFerus (Figure 2d) and pDownei (Table 1) were analyzed. Lipidomic analysis found that pSobrinus, pDownei, and pFerus conferred synthesis of Lys-Glc-DAG, but not Lys-PG in S. mitis. Notably, pFerus confers synthesis of Lys-Glc-DAG at a similar level to S. agalactiae MprF (Figures 2b, d). In contrast, pSalivarius conferred synthesis of Lys-PG only and not Lys-Glc-DAG. Importantly, S. salivarius MprF also synthesized Lys-PG at a level most similar to S. agalactiae MprF (Figures 2b, c). Lysine-modified lipids synthesized by S. mitis strains expressing streptococcal MprFs used in this study are summarized in Table 1. In general, streptococcal MprF enzymes with higher amino acid percent (> 60%) to S. agalactiae MprF conferred Lys-Glc-DAG synthesis to S. mitis, but only S. agalactiae MprF conferred both Lys-Glc-DAG and Lys-PG synthesis in S. mitis.

2.2 Restricted Boltzmann Machines can provide sensitive, rational guidance for sequence classification

We set out to understand which specific MprF amino acid residues are involved in enzyme specificity for either the PG or Glc-DAG substrate. To this end, we employed a Restricted Boltzmann Machine (RBM), a probabilistic graphical model which aims to model the probability of data in a data set through statistical connections between the features of the data and a set of hidden units. [17]. This class of models is closely related to another model termed the Boltzmann Machine. The Boltzmann machine formulation is closely related to the Potts model from physics, which was successfully applied in biology to elucidate important residues in protein structure and function [18], one example being the careful tuning of enzyme specificity in bacterial two-component regulatory systems [19, 20]. The Boltzmann Machine formulation from [18], termed Direct-Coupling Analysis (DCA), stores patterns in the form of a pairwise coupling matrix and local field matrix (shown in Equation 2), which has been useful in understanding how likely pairs of residues are to interact in physical space. We hypothesize that S. agalactiae MprF and other MprF enzymes that utilize glycolipid substrates have unique patterns of interacting residues relative to those enzymes that use only PG as a substrate.

While modeling how pairs of residues covary is powerful, many important features of enzyme function like allostery or enzyme specificity (which is the purpose of this analysis) may involve more than two residues [21], and in these situations it can be difficult to connect sets of residues using solely pairwise parameters [22]. Therefore, a more global approach is required to address our hypothesis. To this end, a RBM has recently been developed for use in the learning of protein sequence data, with special attention given to producing a model which produces sparse, interpretable, and biologically meaningful representations of these higher-order statistical couplings within the hidden units [14]. The design of this RBM can be seen in Figure 3, where the model architecture is represented by purple dots and green triangles. The dots are the “visible” layer, which take in input sequences and encode them into the “hidden” layer, where each triangle represents a separate hidden unit. The lines connecting the visible and hidden layers show that each hidden unit can see all the visible units (the statistics are global), but they cannot see any of the other hidden units, meaning the hidden units are mutually independent. This global model with mutually independent hidden units (see the marginal distribution form shown in Equation 3) allows the following useful properties: higher-order couplings between residues (sets of coupled residues instead of pairs); bimodal hidden unit outputs given the training data set (Equation 4); sparse, interpretable configurations of the hidden units (w_μ) (Equation 5) and a compositional representation of the hidden unit weights [23], where input sequences are modeled largely through combinations of hidden units instead of single, highly activating units.

Figure 3:

The compositional representation is a critical feature of this method, allowing us to find independently coupled sets of residues which in combination describe these protein sequences, and it is achieved through the particular choices of training parameters. More details on the training parameter choices are found in Methods and Materials. In particular, we focus on the scoring of sequences using individual hidden units w_μ:

where v is an input sequence vectorized through one-hot encoding, to produce a hidden unit activation (a single number) which depends on the input sequence residues at position i and the value of the hidden unit weight matrix for that residue at that position.

We used this methodology to study the sequences in the family of MprF, which is a large protein composed of independent flippase and synthase domains [9, 24]. In our analysis, we restrict ourselves to the cytosolic region (see yellow box in Figure 4a), also known as the aminoacyl-PG synthase domain, which constituted the Pfam family DUF2156 (now renamed LPG synthase C) [25]. One hidden unit in particular is highlighted, where we find bimodal activation (Figure 4b) of the hidden unit (Figure 4c) when Equation 1 is evaluated with the training sequences. Interpreting the sequence histograms, a positive value activation in Figure 4b corresponds to a sequence containing combinations of the residues listed in the positive portion of the weight diagram in Figure 4c, whereas a negative activation corresponds to residues in the negative portion. For display purposes, the residue weights shown in all of our hidden units correspond to the positions where the magnitude of the visible positions w_iμ are greater than some threshold proportional to the largest magnitude position (position 212); the hidden units has values for all possible amino acids in the entire protein length, but they are considerably smaller than this threshold and generally very low.

Figure 4:

Two positions in particular, 152 and 212, in the positive portion of the hidden unit correspond to residues 684 and 742, which were experimentally demonstrated to be relevant for aminoacylated-PG production by the Bacillus licheniformis and Pseudomonas aeruginosa cytosolic synthase domain variants [26]; in B. licheniformis the residues which form a complex with the lysine from Lys-tRNA^Lys are shown in Figure 4c. We have highlighted these two aforementioned variants as well as the S. agalactiae variant for demonstration of sequence scoring through Equation 1 in Figure 4a. We used this hidden unit to create a filtered data set; only sequences that contained this necessary motif were of interest to us, and we subsequently trained a new model using only sequences with a catalytic hidden unit activation greater than 2 (because it cleanly split the bimodal distribution). We outline this general workflow in Figure 3. We also note that the full length sequences of the sequences included in our new dataset are enriched with sequences which contain the bacterial flippase domain (shown in red in Figure 4a), and the sequences we removed rarely contain the flippase domain (see Figure S2 for domain analysis).

2.3 RBM hidden units identify plausible sequence residues for functional characterization

After filtering the data set and training another model, we set out to find hidden units which could describe MprF specificity. To this end, we identify hidden units which are sparse, have high overall magnitude, and involve residues which could be plausibly linked to the lipid binding specificity (see Methods for more details). Through this search method we identified the two hidden units shown in Figure 5a. They both have high sparsity and high magnitude relative to the majority of other hidden units in the model, and importantly they correspond to residues localized to regions of the catalytic domain that had been previously identified through Autodock experiments using a PG molecule with side chains C5:0/C8:0 [26] (highlighted in Figure 5a). The scoring of the training set of sequences with these hidden units is shown in Figure 5b, and we see a clear separation between the cytosolic domain sequences which produces Lys-Glc-DAG (from S. agalactiae) and two domain variants which produce only aminoacylated-PG (from P. aeruginosa and B. licheniformis). Importantly, these output activations depend on a greatly reduced subset of the full sequence and do not correspond to clustering on total sequence identity (Figure S1).

Figure 5:

Using these two hidden units as sequence classifiers with a potential link to the substrate specificity of cytosolic domain variants, we used Figure 5b to propose sequences in quadrant three (Q3), the quadrant occupied by the S. agalactiae MprF variant with this novel glycolipid-modifying function. We selected a sequence which was distant in terms of sequence identity from the Streptococcus genus (Table S1) and which has been implicated in human pathologies [27], Enterococcus dispar MprF. We also emphasize that the two hidden units we selected are from a set of 300 total hidden units, and that while other hidden units may differentially classify our three sequences, without the additional link to lipid binding provided by the Autodocking experiment there would be only weak justification for using their activations to make predictions.

2.4 Enterococcal MprF enzymes synthesize Lys-Glc-DAG, Lys-PG, and a novel cationic glycolipid Lys-Glc₂-DAG

Enterococcus dispar ATCC 51266 was identified as a possible candidate for novel lipid synthesis through the analyses described above. Remarkably, our lipidomic analysis found that E. dispar synthesizes a novel cationic lipid, lysyl-diglucosyl(Glc₂)-diacylglycerol (Lys-Glc₂-DAG) (Figure 1, Figure 6a,b), as well as Lys-Glc-DAG and Lys-PG (Figure 6c, Figure S3). The structural identification of Lys-Glc-₂-DAG was supported by exact mass measurement ([M+H]⁺ observed at m/z 1047.730; calculated m/z 1047.731) and tandem MS (Figure 6b). Lys-Glc-₂-DAG is unusually polar and charged for a lipid; chromatographically, Lys-Glc₂-DAG is highly retentive on a silica-based normal phase column and elutes off the column at the very end of the gradient (Figure S3). The use of an amino-based HPLC column led to much improved elution profiles for Lys-Glc₂-DAG and other glyco- and lysine lipids (Figure 6).

Figure 6:

The discovery of Lys-Glc₂-DAG in Enterococcus dispar led us to analyze other Enterococcus strains: E. faecium 1,231,410, E. faecalis T11 and OG1RF E. gallinarium EG2, E. casseliflavus EC10, and E. raffinosus Er676. All Enterococcus strains examined were found to synthesize Lys-PG, Lys-Glc-DAG and Lys-Glc₂-DAG at varying levels (Figure 7a and Table 2). MprF-dependent Lys-PG production by E.faecalis and E. faecium was previously reported [28, 29], but cationic glycolipids production has not been previously reported. E. gallinarium EG2 and E. casseliflavus EC10 encode a single copy of mprF gene, while the other enterococci we tested encode two distinct mprF genes (mprF1 and mprF2) [28, 29]. When these MprF sequence variants are plotted using our chosen RBM hidden units, we see that the two single copy strains plot in quadrant two (Q2), while for the two copy sequences, each MprF variant maps to separate coordinates (Figure S4). Both E. raffinosus Er676 MprF variants plot to quadrant two (with the lowest Weight 2 value being 0.34). The E. faecium and E. faecalis variants have copies plotting directly in quadrant one and variants plotting farther out from the sequence cluster; we note that the two displaced variants are both lacking the characteristic aspartic acid and glutamic acid amino acids at positions 100 and 48 in the hidden units (Figure 5a), which are the dominant contributions for determining quadrant one occupancy.

Table 2:

Figure 7:

To clarify which variants were determinants of the novel cationic lipid production, lipidomic analysis on an E. faecalis OG1RF mini-mariner (EfaMarTn) transposon mutant [30] with a transposon insertion within OG1RF_10760 mprF (referred to as mprF 2 in previous studies [28]; OG1RF_10760::Tn) revealed that Lys-PG, Lys-Glc-DAG and the newly identified Lys-Glc₂-DAG were absent when the synthase domain of the gene was disrupted (Figure 7b, Table S1). This led us to the conclusion that E. faecalis mprF 2 is necessary for the synthesis of the three Lys-lipids in E. faecalis. When the empty vector (pABG5) was transformed into OG1RF_10760::Tn, no Lys-lipids synthesis was observed as expected (Figure 7c). This phenotype was complemented by expression of E. faecalis mprF 2 from a plasmid (pOGMprF2) (Figure 7d). Additionally, expression of E. faecium mprF1 (EFTG_00601) (pEfMprF1) restored Lys-PG and Lys-Glc₂-DAG synthesis to the E. faecalis OG1RF_10760::Tn strain (Figure 7e), while E. faecium mprF2 (EFTG_02430) did not rescue Lys-lipid synthesis (Figure 7f). Taken together, these data indicate that E. faecalis MprF2 (OG1RF_10760) and E. faecium MprF1 (EFTG_00601), like S. agalactiae MprF, act on both phospholipid and glycolipids substrates, additionally synthesizing the novel cationic glycolipid, Lys-Glc₂-DAG.

2.5 Hidden unit 2 is correlated with Glc_N-DAG substrate activity in MprF

When we started our RBM analysis, we had substrate specificity information for only the quadrants one and three, and through identifying the novel substrate in E. dispar we expanded our testing to quadrant two by testing Enterococcus strains. To expand our RBM analysis, we tested specific sequences which were found in quadrant four (Q4). We conducted lipid analysis on Ligilactobacillus salivarius ATCC 11741, Levilactobacillus brevis ATCC 14869, Lacticaseibacillus rhamnosus ATCC 7469, Lactobacillus casei ATCC 393, and Lacticaseibacillus paracasei ATCC 25302. All strains tested synthesize Lys-PG (Table 2). Only Levilactobacillus brevis and Lacticaseibacillus paracasei synthesized Lys-Glc₂-DAG, L. paracasei synthesized low levels of Lys-Glc₂-DAG. An additional quadrant one (Q1) strain Exiguobacterium acetylicum UTDF19-27C, was analyzed by lipidomics. It was found to synthesize only Lys-PG. This MprF falls near Staphylococcus aureus RN4220 and Bacillus subtilis 168 in our plot, both of which synthesize only Lys-PG [31, 24] (and experimentally confirmed in this study). This further supports that strains located in quadrant one of Figure 8 are capable of Lys-PG synthesis.

Figure 8:

Summaries of all functional species/variants tested and the Lys-lipids they produced are shown in Table 2, and plotted in the two RBM hidden units in Figure 8. We list for each MprF variant which of the three tested lipids it acts on, and in Figure 8 we indicate whether a variant has only PG as a substrate or if it operates on Glc_N-DAG as a substrate irrespective of its PG activity. With this plot, we can see that sequences which plot with positive activations in hidden unit 2 are more likely to have Glc_N-DAG as a substrate (8/9), and sequences with negative values predominantly use PG (7/11), with a Fisher’s exact test p-value of 0.028 (Table S2). Hidden unit 1, however, is not correlated with Glc_N-DAG specificity, with a Fisher’s exact p-value of 0.67. Therefore we conclude that hidden unit 2, and thus the high magnitude sequence positions identified by it, are correlated with Glc_N-DAG specificity.

3 Discussion

The utilization of the S. agalactiae MprF amino acid sequence as a query for BLASTp was a simple method to identify various streptococcal MprFs and further our understanding of MprF specificity. This allowed us to identify three MprFs that can synthesize the highly cationic Lys-Glc-DAG. Importantly, although highly similar to S. agalactiae MprF, this method did not identify other streptococcal MprFs capable of lysine addition to both glycolipid and phospholipid substrates (generating Lys-Glc-DAG and Lys-PG) nor did it lead to the identification of a novel cationic glycolipid found in Enterococcus. Sequence logos [32] may identify important sequence features, but largely in instances where amino acids are highly conserved. When amino acids must interact to produce a specific function, yet the exact identity of those amino acids is not a requirement, methods which model the strength of amino acid coevolution in sequence sets can find important signals which in the sequence logo might appear to be negligible. For this reason we used the Restricted Boltzmann Machines (RBM) method [14], which can model statistical connectivity of residues at multiple sites in the protein sequence.

The RBM can be used in a number of ways to study sequence data. Firstly, it allowed us to filter the sequence data using hidden units (Figure 4); the two hidden units which formed the basis of our specificity analysis were not found in models trained on the Pfam dataset when it was initially pulled from Pfam, and this was potentially due to the presence of sequences without relevant catalytic function adding noise to the dataset. The most prominent feature for our purposes is the ability to use the hidden units learned through RBM training to classify sequences. The bimodal nature of this particular model provides a relatively simple interpretation which is straightforward to use for prediction. We show in Figure 5 that the combination of clustering and interpretable hidden unit structure allows us to identify a small subset of residues within a structure, grounding our statistical clustering in MprF’s structural features. The combination of experimental evidence and evolutionary information provided a strong rationale for selecting organisms for further study, ultimately leading to the discovery of novel cationic lipid biosynthesis in Enterococcus dispar (Figure 6, Figure S3).

One goal of this study was to find the structural determinants of Glc_N-DAG activity by MprF, and we found that hidden unit 2 from our RBM analysis was predictive of this. However, there are notable exceptions such as the E. faecium allele EGTF 00601, which has sequence features placing it in quadrant one yet has Glc_N-DAG activity. Therefore, identifying the exact sequence determinants of Glc_N-DAG specificity remains to be fully elucidated. One explanation could be our exclusion of the N-terminal flippase domain and the non-cytosolic region of the C-terminal domain in the RBM analysis, where important residue interactions determining specificity could occur and is under further investigation. Additionally, a challenge in experimental validation will be to standardize the method of studying mprF alleles from diverse organisms. When using native organisms, lipids produced may vary during the growth phase of the organisms, the media they are grown in, and whether a membrane stressor is present (i.e. inducible production of specific lipids). For this study we used stationary phase cultures of equal volume and bacteria were grown in their respective standard laboratory media. The use of a heterologous host for expression of mprF alleles is a method of standardization, but may result in not identifying novel lipids that would be found in their native bacteria since the appropriate lipid substrates may not be present [16, 12]. A combination of both approaches (native organisms and heterologous hosts) may be required. Additionally, spike-in controls might be implemented to help us to more quantitatively understand the output specificity of these MprF variants. Ultimately, mutational studies of the residues identified by the RBM will allow for identification of the residues critical for specificity, though likely the residues will depend in part on the primary sequence being altered.

The general method of RBM-guided exploration has applications in enzyme design, for example, in enhancing methods like directed evolution [33]. Concretely, the residues identified through hidden units can be thought of as templates for mutational exploration, greatly reducing the search space for assessing which residues are critical for specificity determination. Enzymes like MprF that produce lipids of modified charge have potential applications in biotechnology and therapeutics, where the design of lipid nanoparticles with different charge and chemical properties of the head group have important physiological implications [34]. Notably, the identification of Lys-Glc₂-DAG opens new research avenues, particularly in the fields of antibiotic resistance and host-pathogen interactions, wherein lipid modification and cationic lipids play key roles. With respect to human health, it would be relevant to investigate whether Lys-Glc₂-DAG plays a role in enterococcal resistance to the last-line membrane-active antibiotic, daptomycin.

4 Methods and Materials

4.1 Bacterial Strains and Growth Conditions

See Tables S1 and S3 for a full list of bacterial strains used in this study. Streptococcus mitis NCTC12261 (referred to as SM61 here) [35] was grown in Todd Hewitt Broth (THB) at 37^°C with 5% CO₂. Escherichia coli DH5α [36] and Bacillus subtilis 168 [37] was grown at 37^°C with shaking at 225 rpm in lysogeny broth (LB). Staphylococcus aureus RN4220 [38] was grown at 37^°C with shaking at 220 rpm in Trypic Soy Broth (TSB). Streptococcus agalactiae COH1 ATCC BAA-1176 [39] and CJB111 [40, 41] were grown in THB at 37 ^°C. Enterococcus faecium 1,231,410 [42], Enterococcus faecalis T11 [43] and OG1RF [44], Enterococcus gallinarium EG2 [42], Enterococcus casseliflavus EC10 [42], and Enterococcus raffinosus Er676 [45] were grown at 37^°C in Brain Heart Infusion (BHI). Enterococcus dispar ATCC 51266 [46] was grown at 30^°C in BHI. Exiguobacterium acetylicum UTDF19-27C was grown in BHI at 37 ^°C. Ligilactobacillus salivarius ATCC 11741 [47], Lacticaseibacillus rhamnosus ATCC 7469 [48], Lactobacillus casei ATCC 393 [49] and Lacticaseibacillus paracasei ATCC 25302 [48] were grown in Lactobacilli MRS Broth at 37^°C and 5% CO₂. Levilactobacillus brevis ATCC 14869 [50] was grown in Lactobacilli MRS Broth at 30^°C, 5% CO₂. Bacillus licheniformis ATCC 14580 was grown in nutrient broth at 37^°C with shaking at 225 rpm [51]. A transposon mutant of E. faecalis OG1RF from [30] with a mini mariner transposon (EfaMarTn) insertion in OG1RF_10760 (OG1RF_10760::Tn) was grown on BHI supplemented with chloramphenicol at a concentration of 15 μ/mL. The mutant was confirmed to be erythromycin-sensitive. The transposon mutant location was confirmed through Sanger sequencing performed by the Genome Center at the University of Texas at Dallas (Richardson, TX).

Where appropriate for plasmid selection, kanamycin (Sigma-Aldrich) was added. For E. coli, a concentration of 50 μg/mL was used, for S. mitis a concentration of 300 μg/mL, and for E. faecalis a concentration of 500 μg/mL was used.

4.2 Routine Molecular Biology Procedures

Genomic DNA was extracted as done previously in [52] and [53]. All PCR reactions used Phusion polymerase (Thermo Fisher Scientific) and Phusion 5x HF buffer (Thermo Fisher Scientific) in a Veriti PCR machine (Applied Biosystems). List of primers used are found in Table S4. Gibson assemblies were completed using a 2x HI-FI Assembly master mix following the manufacturer’s protocol (New England Biolabs). PCR clean-up was done using the GeneJET PCR Purification Kit (Thermo Fisher Scientific) per manufacturer protocol. Plasmid extractions were performed per manufacturer protocol using the GeneJET Plasmid miniprep kit (Thermo Fisher Scientific). All plasmid constructs were confirmed by Sanger sequencing at the Massachusetts General Hospital CCIB DNA Core facility or by Illumina sequence at SeqCenter. DNA concentrations were measured using Nanodrop (Thermo Fisher Scientific) or Quibit 2.0 (Invitrogen by Life Technologies). Optical Density at 600 nm (OD_600nm) was measured in a disposable cuvette (Thermo Fisher Scientific) using a spectrophotometer (Thermo Scientific Genesys 30).

Gibson assembly was performed using pABG5 as previously described [11]. Transformation of plasmids into E.coli, S. mitis, and E.faecalis was performed as previously described in [11, 16, 54, 55].

4.3 Acidic Bligh-Dyer Lipid Extractions

Bacteria were grown for approximately 16 hours overnight in 15 mL of an appropriate culture medium. After growth, OD_600nm measurements were taken, and cells were pelleted at 4,280 x g for 5 minutes at room temperature in a Sorvall RC6+ centrifuge. The supernatant was tipped out, and the cell pellet was washed and resuspended in 1X Phosphate Buffered Saline (PBS). The cells were pelleted again, and all the supernatant was aspirated out. The cell pellet was stored at -80^°C until lipid extraction. An acidic Bligh-Dyer lipid extraction was performed as previously reported [16, 12]. Briefly, cell pellets were resuspended in 0.8 mL of 1x Dulbecco’s PBS (Sigma-Aldrich). Cells were transferred to a 9 mL glass tube with a Teflon-lined cap (Pyrex). 1 mL chloroform (MilliporeSigma) and 2 mL methanol (Thermo Fisher Scientific) were added to create the single-phase Bligh-Dyer. Tubes were vortexed every 5 minutes for 20 minutes, then centrifuged at 500 x g for 10 minutes at room temperature. The supernatant was transferred to a new 9 mL glass tube. 100 μL of Hydrochloric acid (Thermo Fisher Scientific) was added followed by 1 mL of chloroform and 0.9 mL of 1x Dulbecco’s PBS to create the two-phase Bligh-Dyer. Tubes were gently mixed and centrifuged at 500 x g for 5 minutes at room temperature. After centrifuging, the bottom layer was extracted to a new tube and dried under nitrogen gas and stored at -80^°C prior to lipidomic analysis. Lipid analyses were repeated in biological triplicate aside, from OG1RF_10760::Tn extractions, which were performed once.

4.4 Analysis of lysine lipids by an amino column-based Liquid Chromatography-Electrospray Ionization-Mass Spectrometry (LC-ESI MS): A new method

An amino column-based normal phase LC-ESI MS was performed using an Agilent 1200 Quaternary LC system coupled to a high-resolution TripleTOF5600 mass spectrometer (Sciex, Framingham, MA). A Unison UK-Amino column (3 μm, 25 cm × 2 mm) (Imtakt USA, Portland, OR) was used. Mobile phase A consisted of chloroform/methanol/aqueous ammonium hydroxide (800:195:5, v/v/v). Mobile phase B consisted of chloroform/methanol/water/ aqueous ammonium hydroxide (600:340:50:5, v/v/v/v). Mobile phase C consisted of chloroform/methanol/water/aqueous ammonium hydroxide (450:450:95:5, v/v/v/v). The elution program consisted of the following: 100% mobile phase A was held isocratically for 2 min and then linearly increased to 100% mobile phase B over 8 min and held at 100% B for 5 min. The LC gradient was then changed to 100% mobile phase C over 1 min and held at 100% C for 3 min, and finally returned to 100% A over 0.5 min and held at 100% A for 3 min. The total LC flow rate was 300 μl/min. The MS settings were as follows: Ion spray voltage (IS) = 5000 V, Curtain gas (CUR) = 20 psi, Ion source gas 1 (GS1) = 20 psi, De-clustering potential (DP) = 50 V, and Focusing Potential (FP) = 150 V. Nitrogen was used as the collision gas for MS/MS experiments. Data acquisition and analysis were performed using Analyst TF1.5 software (Sciex, Framingham, MA).

4.5 Generation of mprF expression plasmids

The S. downei mprF nucleotide sequence (Locus tag HMPREF9176 RS03810) and S. ferus mprF nucleotide sequence (Locus tag A3GY RS0106165) were used to design synthetic geneblocks (GeneWiz from Azenta Life Sciences). A 20 nucleotide 5’ extension and a 20 nucleotide 3’ extension complementary to the pABG5 plasmid were added to enable Gibson assembly of the inserts with pABG5. Sequences for the geneblocks can be found in Table S5. The S. salivarius mprF (Locus tag SSAL8618 04345) was amplified from Streptococcus salivarius ATCC 7073. A 20 nucleotide 5’ extension and a 20 nucleotide 3’ extension complementary to the pABG5 plasmid was added to the insert. The same protocol was used for generation of pSobrinus, pEFM-prF1, pEFMprF2, and pOGMprF2 with slight modifications. S. sobrinus mprF (Locus tag DLJ52 05040) was amplified from Streptococcus sobrinus ATCC 27352. E. faecium mprF1 (Locus tag EFTG_00601) and E. faecium mprF2 (Locus tag EFTG_02430) were amplified from Enterococcus faecium 1, 231, 410. E. faecalis mprF2 (Locus tag OG1RF_10760) was amplified from E. faecalis OG1RF.

4.6 Boltzmann Machines and Restricted Boltzmann Machines for Protein Sequences

Formulation of models

The Boltzmann Machine, as described in [18], is defined as:

which is a probability distribution defined on the pairwise interactions (the couplings e_ij) between the positions (v_n) in an input sequence (v) and a local field term h_i. The pairwise coupling matrix is akin to a covariance matrix, and the local fields are similar to single site frequency measures.

The joint probability distribution of the Restricted Boltzmann Machine as described in [14] is as follows:

Here, the vector v again represents the input sequence data, g(v_i) is a local field which controls the conditional probability of the input data,𝒰_μ is the hidden unit potential/activation, and the terms h_i and w_iμ(v) couple the input variables with the hidden unit variables. In this way, interactions between residues are mediated by a relatively smaller set of hidden units which do not interact directly with each other, leading to a bipartite structure as opposed to the Boltzmann machine’s fully connected pair based formulation.

Particularly important is the potential function 𝒰, which in this work is the dRELU function defined as:

The four parameters γ_μ,±, θ_μ,± are learned through the inference procedure, and allow the positive and negative components of the hidden units w_μ to be separately gated, which allows learning of bimodal distributions for the activations h_μ. Another key feature is the sparsity regularization on the weights applied during model inference:

where q is the number of amino acids plus a gap character (21) and the hyperparameter can be increased to induce sparsity and lower weight magnitude. Equations 4 and 5 together guide the representation of the hidden units to bimodal outputs with interpretable features.

Dataset acquisition and preprocessing

Multiple sequence alignment used for model training was the DUF2156 domain acquired from Pfam [25] (in later revisions renamed LPG synthase C). Starting with this MSA, it was processed to include only amino acid characters and gaps, excluding sequences with ambiguous or non-standard characters. Additionally, sequences with contiguous strings of gaps greater than 20% of the full sequence’s length were removed. This left 11,507 sequences with a length of 298. After the data cleaning process described in Figure 4, a total of 7,865 sequences were used to train our final model.

Model training

The model previously described [14] was utilized, specifically the Python 2.7 version freely available on Github (https://github.com/jertubiana/ProteinMotifRBM). Extensive testing of the combinations of the number of hidden units, L2/L1 regularization, and learning rate were performed. Models were trained in triplicate with different random seeds, as the training procedure converges to slightly different weight values depending on their initial random configuration. Training was considered successful when the model consistently produced similar sets of weights across different random seeds, and the weights had the previously mentioned sparsity and magnitude indicating compositional representation. The final model used in this work was trained for 2000 epochs with a learning rate of 0.1, a learning rate decay of 0.33, L2/L1 parameter set to 1.0, and 300 hidden units.

Model Weight selection

To find hidden units which were relevant to function we computed the L1 norm of each hidden unit individually, then ranked them. The hidden units with the largest magnitude were looked at first, and these were typically the sparsest of all hidden units. Hidden units were chosen which involved sequence coordinates implicated in our function of interest. Specifically, locations identified through Autodock [26] where the lipid was likely to interact, and a small radius around this region to select a small set of coordinates. We chose the only hidden units which had both overlap with multiple residues in this chosen radius and predictive power (separation) for the three examples we had to start with.

Software Accessibility

DUF2156 domain sequences were retrieved from Pfam33 [25]. Raw sequence datasets are available as Supplementary Material.

Acknowledgements

J.M. and F.M. acknowledge support from the National Institutes of Health (NIH R35GM133631). F.M. acknowledges support from the National Science Foundation CAREER award (MCB-1943442). We acknowledge support from the National Institutes of Health grant R01AI178692 (Z.G., F. M. and K.P.) and R01AI148366 (Z.G. and K.P.). K.P. acknowledges support from the Cecil H. and Ida Green Chair in Systems Biology Science. L.R.J acknowledges support from the American Heart Association (23POST1013835)