Prediction of Burkholderia pseudomallei DsbA substrates identifies potential virulence factors and vaccine targets

Identification of bacterial virulence factors is critical for understanding disease pathogenesis, drug discovery and vaccine development. In this study we used two approaches to predict virulence factors of Burkholderia pseudomallei, the Gram-negative bacterium that causes melioidosis. B. pseudomallei is naturally antibiotic resistant and there are no clinically available melioidosis vaccines. To identify B. pseudomallei protein targets for drug discovery and vaccine development, we chose to search for substrates of the B. pseudomallei periplasmic disulfide bond forming protein A (DsbA). DsbA introduces disulfide bonds into extra-cytoplasmic proteins and is essential for virulence in many Gram-negative organism, including B. pseudomallei. The first approach to identify B. pseudomallei DsbA virulence factor substrates was a large-scale genomic analysis of 511 unique B. pseudomallei disease-associated strains. This yielded 4,496 core gene products, of which we hypothesise 263 are DsbA substrates. Manual curation and database screening of the 263 mature proteins yielded 81 associated with disease pathogenesis or virulence. These were screened for structural homologues to predict potential B-cell epitopes. In the second approach, we searched the B. pseudomallei genome for homologues of the more than 90 known DsbA substrates in other bacteria. Using this approach, we identified 15 putative B. pseudomallei DsbA virulence factor substrates, with two of these previously identified in the genomic approach, bringing the total number of putative DsbA virulence factor substrates to 94. The two putative B. pseudomallei virulence factors identified by both methods are homologues of PenI family β-lactamase and a molecular chaperone. These two proteins could serve as high priority targets for future B. pseudomallei virulence factor characterization.


Introduction
Burkholderia pseudomallei is a Gram-negative soil dwelling saprophyte, and an opportunistic pathogen responsible for the severe tropical disease melioidosis [1]. B. pseudomallei infections are difficult to treat [2][3][4] and are intrinsically resistant to almost all available antibiotics [5][6][7][8] point is thought to have evolved to limit formation of mis-matched disulfide bonds and therefore misfolded proteins [42,43].
In the present study, we used two approaches to identify potential B. pseudomallei DsbA substrates for further study as virulence factors. In one approach, we used computational methods to generate a curated list of 263 putatively extra-cytoplasmic proteins from the core genome of 511 disease-associated isolates of B. pseudomallei, 81 of which were predicted to be virulence-associated. In the second approach, 15 candidate DsbA virulence factor substrates were identified by sequence homology to known DsbA virulence factor substrates in other bacteria.

Genomic analysis to predict B. pseudomallei DsbA virulence factor substrates
In this approach, our strategy was to cast a wide net initially, by determining the pangenome of disease-associated isolates of B. pseudomallei, and then filtering from that the core genome (i.e. the highly conserved genes). The disease-associated B. pseudomallei core genome should then be enriched in conserved virulence factors. At the time of this analysis the NCBI database [44] contained 1577 B. pseudomallei isolates. Metadata notation allowed selection of 512 isolates associated with disease (i.e. isolates from swabs/clinical isolates: accession numbers of these are given in S1 Data); other genomes were discarded. We note that only 355 of the 512 isolates were tagged 'pathogen' in the NCBI database indicating a discrepancy between NCBI assignment and user-uploaded metadata. Analysis of the pangenome, that is the core, accessory and unique genes of these 512 B. pseudomallei isolates (see Table 1), revealed two identical strains. Therefore for the remainder of this analysis, only the 511 unique strains were used.
We found that the core genome consisted of 4,496 genes (see S2 Data) or 22.49% of the total 19,991 pangenome. This analysis largely agrees with a previous pangenomic analysis which extrapolated a modelled core genome of 4,568±16 genes from a much smaller set of 37 isolate genomes [45]. In that approach, modelling was used to predict the core genome if the number of isolates was expanded. Our approach gives an exact number because all 4,496 genes were found in all 511 genomes. Notably, the dithiol oxidase redox enzyme pair DsbA and DsbB and the disulfide isomerase redox relay enzymes DsbC and DsbD were all identified as core genes.
We then used the B. pseudomallei core genome for further analysis, because it encodes highly conserved proteins-a key criteria for selecting vaccine or anti-virulence targets.
From these 4,496 core genes, 726 were predicted to encode proteins with a signal sequence and which are therefore likely to be exported out of the cytoplasm and into the periplasm The distribution of cysteines in B. pseudomallei cytoplasmic and extra-cytoplasmic proteins was calculated for the pangenome (total of 19,991 genes) and the core genome (4,496 genes) (refer to Table 1). In cytoplasmic B. pseudomallei proteins, cysteine distribution followed a Poisson law peaking at zero for the pangenome and at one for the core genome (denoted by the orange lines in the histograms on Fig 2A and 2B). This distribution changed for extra- Similarly, predicted cytoplasmic proteins are represented as orange lines. Panels C and D represent the normalised frequency of cysteine-containing extra-cytoplasmic proteins. The blue line in panel D peaks for proteins with 2, 4, 6 and 8 cysteines suggesting a preference for an even number of cysteines. This trend is not observed as strongly in panel C, where a clear peak can only be seen for two and eight cysteines. The normalised frequency was calculated by dividing the number of extra-cytoplasmic proteins (having N number of cysteines) by the total number of proteins with N cysteines (N being a number between 0-20 as per the data points in C and D above). https://doi.org/10.1371/journal.pone.0241306.g002

PLOS ONE
Identification of potential B. pseudomallei DsbA substrates involved in virulence cytoplasmic B. pseudomallei proteins. For the core genome (blue bars Fig 2B), B. pseudomallei proteins with an even number of cysteines were over-represented compared to a typical Poisson distribution. As extra-cytoplasmic proteins represent a small fraction of the total number of the translated core genome and pangenome (16% and 11.5% of all proteins, respectively), we also analysed the normalised frequency (Fig 2C and 2D). The core genome normalised cysteine distribution reveals a sawtooth pattern with a preference for even number of cysteines with peaks for two, four, six and eight cysteines (Fig 2D). In contrast, the pangenomic normalised cysteine distribution for extra-cytoplasmic B. pseudomallei proteins does not indicate a strong preference for even number of cysteines ( Fig 2C). Overall, the saw-tooth pattern observed in Fig 2B and 2D is similar to that described for E. coli exported proteins [42] although not as pronounced.
Functional assignment of core, extra-cytoplasmic, putative DsbA substrates The next step in the genomic analysis was to predict which of the 263 putative DsbA substrates are associated with virulence. Of the 263 selected proteins, 44 were annotated as hypothetical/ uncharacterised. The remaining 219 proteins include ABC transporter-related proteins, housekeeping proteins like cytochrome C, proteins required for motility such as flagellar and fimbrial proteins, enzymes such as collagenase, peptidases and proteases, as well as antibiotic resistance enzymes, β-lactamases. Many oxidoreductases were also present including DsbA, DsbD and others such as Gfo/Idh/MocA family, glycerol-3-phosphate dehydrogenase GpsA and thioredoxin-like TlpA oxidoreductases. Redox enzymes such as DsbB and DsbC are core genes with signal sequences, and they have catalytic rather than structural disulfides. These two enzymes are not identified as DsbA substrates in our filter as they have an odd number of cysteines.
The list of 263 proteins with an even number of cysteines was initially screened against the Virulence Factor DataBase (VFDB) [46], the Burkholderia Genome Database (BGD) [47] and against a list of B. pseudomallei virulence genes identified by previous studies [48,49]. Of the 263 putative DsbA substrates two are closely related to virulence factors from the VFDB (flagellar proteins FlgA and FlhG), six are close homologues to proteins identified previously by Moule et al. [48] five reported by Holden et. al. [49] and one identified from the BGD, giving a total of 14 virulence factors identified through cross-analysis (see S1 File for a full list). It was also noted that two of the 14 identified putative virulence factors, were homologous to the same collagenase (BPSS0666).
By further inspection of the 263 core, putatively extra-cytoplasmic DsbA substrates, and by using the GO descriptions to aid in predicting protein functions, 73 sequences were identified which were virulence-associated (Table 2). These include serine-type endopeptidases [50] associated with adherence, choline binding proteins [51], N-carbamoylputrescine amidase, essential for production of putrescine, a component of Gram-negative cell walls of pathogens and key virulence [52][53][54][55], and many proteases and peptidases.
Our GO description analysis identified many more potential virulence associated genes (73 in total) as compared to the 14 found on the VFDB, BGD and through literature [46][47][48][49]. Six of the putative virulence factors were common to both our GO analysis and to previous analyses (See S1 and S2 Files for full lists), to give a total of 81 identified putative virulence factors.
The 73 putative virulence factors sequences identified by our GO analysis, along with the 8 additional sequences found in the literature and databases, are grouped in Table 2 (the six common sequences found on both lists are displayed with the GO analysis results and underlined). Interestingly, one protein is annotated as a DNA transcriptional regulator from the AraC family (WP_004524330.1) a suspected cytoplasmic protein, although no experimental subcellular localisation information can be found [47]. As a cytoplasmic protein cannot be a substrate of the periplasmic DsbA protein, further experimental studies are needed to confirm the localisation.

Sequence homology prediction of B. pseudomallei DsbA virulence factor substrates
To complement the genomic analysis described above we used a second approach to identify DsbA substrates, by screening all B. pseudomallei genomes uploaded on NCBI [56] (taxid 28450) for homologues of known DsbA substrates. We implemented this approach because some DsbA substrates might be filtered out using the genomic approach described above if the substrates are not encoded by core genes, or if the gene product has an odd number of cysteines.

Fig 3. Gene Ontology (GO) descriptions of predicted extra-cytoplasmic proteins with an even number of cysteines.
The highest frequency of proteins with an even number of cysteines are integral components of membranes (66 proteins), followed by proteins involved in redox (oxidation-reduction) processes (25 proteins) and proteolysis (20 proteins). For ease of representation and clarity, GO descriptors with less than three counts were excluded from this graph. A complete graph, along with raw values can be found in S1 File. https://doi.org/10.1371/journal.pone.0241306.g003

PLOS ONE
Identification of potential B. pseudomallei DsbA substrates involved in virulence

PLOS ONE
Identification of potential B. pseudomallei DsbA substrates involved in virulence Over 90 DsbA substrates have been reported in the literature. We searched for B. pseudomallei homologues of these DsbA substrates using the following criteria: (i) presence of secretion signal, (ii) at least two cysteines in the mature sequence, (iii) at least 20% identity and (iv) 50% coverage to a known DsbA substrate sequence. After removing duplicates, our analysis found that B. pseudomallei encodes homologues of 15 DsbA substrates ( Table 3). Two of these 15 are DsbA substrates in other Burkholderia species B. cepacia and B. cenocepacia [57][58][59][60]: a metalloproteases, ZmpA and a sulfatase-like hydrolase transferase. In B. cenocepacia, ZmpA is a wide spectrum metalloprotease, thought to cause tissue damage during infection [61].
Over 50 DsbA substrates in Francisella tularensis were identified by trapping and co-purifying substrates bound to a DsbA variant [43]. Of these 50, we found nine homologues encoded in B. pseudomallei (see Table 3). These include homologues of the lytic transglycosylase domain containing protein (implicated in peptidoglycan rearrangement) and homologues of two pilin proteins involved in the formation of pilus and flagella. Also present is an MoeB homologue; MoeB is a molybdopterin synthase adenyl transferase (cytoplasmic in E. coli but likely periplasmic in B. pseudomallei due to the twin-arginine translocation (TAT) signal sequence). A PenI family β-lactamase homologue is also found in B. pseudomallei; this is a class A β-lactamase that confers resistance to β-lactams including, in rare cases, ceftazidime (commonly used to treat melioidosis) [62]. A succinate dehydrogenase flavoprotein subunit homologue, found in the bacterial inner membrane and part of the electron transport chain, is also encoded in B. pseudomallei. This protein is cytoplasmically oriented in E. coli, though again the B. pseudomallei version has a TAT signal sequence suggesting a possible periplasmic localisation. ) have B. pseudomallei homologues including a molecular chaperone homologous to PapD and EscC, involved in the formation of the Type III secretion system (T3SS). The T3SS assembly requires DsbA activity in many Gram-negative bacteria, including E. coli and S. typhimurium [63,64]. Finally, a B. pseudomallei protein homologous to the Y. pestis pilus assembly protein Caf1M (a molecular chaperone involved with assembly of the surface capsule of the bacterium) was also identified.
Of the 15 putative B. pseudomallei DsbA substrates identified using this substrate homology method, two were also identified in the genomic pipeline method. These are the PenI (WP_050772403) and a molecular chaperon (WP_102811167).
We then aligned the sequences of the Table 3 B. pseudomallei proteins to identify any possible sequence conservation around the cysteine residues, but no pattern was identified. This lack of peptide sequence motif in DsbA substrates has also been observed in E.coli, demonstrating the difficulty of DsbA substrate prediction [65].

Epitope prediction of virulence-associated proteins
To determine whether the DsbA substrates identified in the two methods above could contribute to vaccination efforts against B. pseudomallei, we also predicted B-cell epitopes, using a structure-informed approach. The sequences of the 81 putative, extra-cytoplasmic DsbA substrates (predicted virulence factors, Table 2) along with the 13 unique, homologous DsbA substrates (Table 3) were screened against the Protein Data Bank (PDB) [66], to identify structurally characterised homologues. Seven of the 94 proteins were found to have at least 80% similarity to a structurally characterised protein. Three of these seven protein structures were from Pseudomonas species, while the other four were from Burkholderia species. Similarity was used rather than identity to account for mutations of functionally similar residues. The seven protein structures were then used as models to predict structurally-informed B-cell epitopes of length 10-32 residues (Table 4 and Fig 4) using the SEPPA 3.0 server. While SEPPA 3.0 is considered the foremost B-cell epitope predictor, the software also accounts for potential glycosylation of the peptide [67], a feature that is mostly absent from bacterial proteins. To ensure that the epitopes identified by SEPPA 3.0 were not the result of erroneous glycosylation interpretation, the epitopes were cross-validated using ElliPro software that does not rely on glycosylation patterns [68]. All hits obtained with SEPPA 3.0 were also identified with ElliPro, with 1-3 residue differences in the starting and ending residues, suggesting that they were not based on wrongly attributed glycosylation patterns. However, we recommend using the more The virulence-associated putative DsbA substrates (Table 2) were screened for �80% similarity to proteins within the PDB to account for substitution of functionally similar residues. The structures were then screened for epitopes using SEPPA 3.0. Fourteen B-cell epitopes of 10 to 32 residues were predicted. https://doi.org/10.1371/journal.pone.0241306.t004

PLOS ONE
Identification of potential B. pseudomallei DsbA substrates involved in virulence  Table 4. Proteins are shown as white surfaces and their respective PDB ID is given in the bottom left corner of each box. The epitope region is highlighted in red and the corresponding homologous stringent list of epitopes identified with SEPPA 3.0 over the much longer list of potential epitopes and antigenic determinants identified with ElliPro. These epitopes provide an interesting list for further evaluation. For example, epitopes from β-lactamase Toho-1 and class D β-lactamase could provide a useful vaccination approach for B. pseudomallei because these directly target antibiotic resistance proteins. Similar approaches have conferred protection against other bacteria in animal models [69][70][71][72].
Vaccination targeting adhesion proteins and essential virulence factors such as FimA [73] and type 1 fimbrial protein [74] is a commonly used approach due to the external localisation of these proteins and their exposure to host immune systems. Anti-fimbrial antibodies have been shown to interfere with function and reduce disease [75,76] and a FimA vaccine provided protection against Streptococcus parasanguis, Streptococcus mitis, Streptococcus mutans and Streptococcus salivarius in rats [77][78][79].
Vaccination against conserved, secreted enzymes such as the triacylglycerol lipase (EstA) and S8 family serine peptidase enzymes may also be a useful strategy. Secreted peptidases are known virulence factors in many pathogenic bacteria [50,80] and vaccines targeting them have attenuated disease in animal models [81,82]. Two triacylglycerol lipases (WP_038741497.1 and WP_038775093.1) were identified as having a structural homologue in the PDB. These two lipases are both core genes and share 78% similarity (72% identity, 87% query cover) and their sequences were both aligned to the same PDB code, resulting in epitope variants of similar sequences.
Finally the UDP-glucose dehydrogenase appears to be a key player in the synthesis of exopolysaccharide in the B. cepacia complex [83], and is suspected to contribute to virulence and cystic fibrosis.

Discussion
In the present study, we analysed genomes from 512 B. pseudomallei isolates specifically associated with disease to identify core putative DsbA substrates and virulence factors. Pangenomic analysis of B. pseudomallei has previously been performed utilising 37 isolates from a variety of isolation sources [45] and concluded the pangenome to be 'open', indicating that new isolates will continually increase the number of total genes, which we found to be the case, based on a pangenome of 19,991 genes from 512 isolates. Previous studies comparing the B. pseudomallei genome with the obligate pathogen Burkholderia mallei (responsible for glanders) and the generally non-pathogenic Burkholderia thailandensis [84][85][86][87], identified several loci likely to be involved in B. pseudomallei virulence. These include the capsular polysaccharide gene cluster and Type III secretion needle complex [87], which were not considered core genes, demonstrating the importance of large-scale analysis.
In the present study, we used two orthogonal approaches to identify a total of 278 putative DsbA substrates, with 94 predicted to be virulence factors (S2 File). Of these, 73 were identified by the genome analysis approach, 8 more via comparison to previous studies and 15 were identified by the DsbA substrate homology approach, with two of the putative 94 DsbA virulence factor substrates identified in both genomic and homology analysis. These two are the experimentally validated bacterial virulence factors and DsbA substrates-a molecular chaperon (reported to be an E. coli DsbA substrate [31]), and a PenI family β-lactamase (reported to be a F. tularensis DsbA substrates) [43].
Delving deeper into the results presents some curious outcomes. For example, the well-characterised E.coli DsbA substrate and virulence factor FlgI [36, 88] was not picked up as a potential B. sequences found in B. pseudomallei are given in one letter code under each respective structure and separated by semicolon when more than one sequence pointed to the same epitope. https://doi.org/10.1371/journal.pone.0241306.g004

PLOS ONE
Identification of potential B. pseudomallei DsbA substrates involved in virulence pseudomallei DsbA substrate by either method, though B. pseudomallei encodes FlgI. The B. pseudomallei FlgI sequence has 4 cysteines in the translated gene product but the predicted mature sequence after cleavage of the signal sequence has just one cysteine. Generally, DsbA does not interact with proteins having just one cysteine. If B. pseudomallei FlgI is a DsbA substrate (that is yet to be tested), then the most likely reasons that it was not identified as a substrate by either of the two methods we used are that (i) the predicted signal peptide is incorrect and/or (ii) the single cysteine of B. pseudomallei FlgI forms an inter-molecular disulfide bond.
The finding that the two orthogonal approaches identified the same two target proteins suggests that there is merit in using different theoretical approaches to select high priority targets for further evaluation (in this case, the PenI family β-lactamase and the molecular chaperon). On the other hand, the fact that there were so few overlaps in the predicted substrates from the two methods raises questions about the filters we applied. Specifically, we found that of the 15 potential substrates identified by the substrate homology method, 5 had an odd numbers of cysteines, whereas the genomic analysis filtered these proteins out of consideration to reduce the number of false negatives. We applied the even cysteine filter because previous reports showed that E. coli exported proteins have a strong preference for an even number of cysteines [42]. This even number of cysteine preference is present in B. pseudomallei exported proteins (Fig 2) though is not as pronounced as in E. coli. By restricting our genomic analysis to core, extra-cytoplasmic B. pseudomallei proteins with an even number of cysteines, some DsbA substrates may therefore have been missed. There is considerable evidence that many virulence factors such as adhesion and motility proteins, toxins and enzymes are extra-cytoplasmic proteins in both Gram-positive and Gram-negative bacteria [31,32,89]. Given that extra-cytoplasmic proteins in the translated core genome of B. pseudomallei have a slight preference for even number of cysteines (Fig 2) and the identification of many virulence-associated proteins within the 263 proteins in the list, the approach taken in this analysis (Fig  1) to identify DsbA substrates was justified. Further, the genomic analysis focused on highly conserved proteins from the core genome; accessory proteins associated with virulence would not be identified using this approach. Nevertheless, the genomic analysis identified homologues of known DsbA substrates in other bacteria, such as the OmpA porin, supporting the use of this approach. However, attempting to identify epitopes from proteins which are not found in every disease-causing isolate may present challenges for anti-virulence and vaccination attempts.
In addition, the genomic analysis identified several proteins of unknown function which could represent novel virulence factors for future studies. Importantly, our theoretical approach was extended to predict structurally-informed surface epitopes for several core gene DsbA substrates for potential vaccine or antibody development (Table 4).
In summary, our in silico analysis combined a substrate homology approach and a genomic analysis approach to identify more than 90 potential B. pseudomallei DsbA virulence factor substrates, two of which we mark as high priority for experimental validation. Future characterization of these proteins will aid our understanding of B. pseudomallei virulence and could provide new targets for anti-virulence drug discovery and vaccine development. The approaches we report here could also be applied to identify potential DsbA virulence factor substrates in other pathogenic bacteria.
fa > output.fa The core genome was then used in the remaining analysis and core DNA sequences were translated into protein sequences using transeq [92] with the following command: transeq -sequence input.fasta -outseq output.fasta -table 11 -frame 1 The core genome was then filtered based on signal sequence and then the sequence of the mature exported protein, as predicted utilising SignalP 5.0 [93,94].

Identification of DsbA substrate homologues in B. pseudomallei
DsbA substrates were also predicted using a substrate homology search. This approach may identify proteins not encoded in the core genome. The B. pseudomallei genome was screened for homologues of known DsbA substrates using BLASTP. A starting list of confirmed DsbA substrates was extracted from the literature [31, 43, 57-61], and their amino acid sequences used in BLAST searches [97] against the NCBI protein database [56] for homologues in B. pseudomallei using default search parameters. In some cases two search proteins identified the same homologue in B. pseudomallei. In these cases only the search protein most similar to the B. pseudomallei homologue is given in Table 3. The results were filtered to select proteins with at least 20% sequence identity and a sequence coverage of at least 50%. Protein sequences with fewer than two cysteines were removed. Exported proteins were selected on the basis of predicted signal sequence (SignalP 5.0 [93]) or experimental evidence of extra-cytoplasmic localisation for the reported DsbA substrate in another Burkholderia species.
Identification of putative virulence factors. ABRicate version 1.0.1 (https://github.com/ tseemann/abricate) [98] was used, along with the virulence factor database (VFDB) [46] to identify the presence of putative virulence factors of the putative paired cysteine gene list. 244 genes identified as virulence-related, on the basis of mutagenesis studies [49,99,100] were also screened against the paired cysteine gene list using blastp version 2.9.0+ [97,101] and results were filtered for �90% coverage and �80% similarity/positives to be considered a putative virulence factor. Additionally, the burkholderia.com virulence database [47] was downloaded and screened against gene lists using blastp version 2.9.0+ with the same filtering conditions.

Cysteine distribution analysis
Fasta files containing either the 19,991 pan genes or the 4,496 core gene of B. pseudomallei with their corresponding amino acid sequences and descriptors were utilised to calculate the distribution of cysteines with a custom Python 3.0 script (available on Github: (https://github. com/gpetit99/cysteineCount_bPseudomallei/blob/master/CysCountFrequency.py"). Briefly, lists of the extra-cytoplasmic protein sequences with signal peptides removed were compared to lists of the protein sequences from the whole genome to create dataframes with either cytoplasmic or extra-cytoplasmic proteins. Proteins were grouped based on the presence or absence of SP, and based on the number of cysteines in the mature protein. To calculate the normalised frequency of cysteines for extra-cytoplasmic proteins, we divided the number of extra-cytoplasmic proteins having N cysteines by the total number of proteins having N cysteines (N being an integer from 0 to 73 -No protein has more than 73 cysteines in the B. pseudomallei translated genome). This analysis was run for the core genome and pangenome independently. Other statistics (e.g. number of proteins in each group) were extracted from the dataframes.

Epitope prediction
The metadata for each of the 263 proteins in the annotated list was manually inspected to select for further analysis a total of 81 proteins likely related to virulence. The sequences of these 81 selected proteins were combined with the 13 unique proteins from the homology analysis (to give 94 unique protein sequences). These were screened against the protein data bank using BLAST (criteria: �80% positive substitutions/similarity used as a threshold) to find structurally characterised homologues. These structural homologues were then used to predict B-cell epitopes using SEPPA 3.0 (http://www.badd-cao.net/seppa3/index.html) with a threshold of 0.1 [67]. Similarity was used rather than identity to account for mutations of functionally similar residues. Predicted B-cell epitopes were accepted if they were 10-32 residues in length, as described in [102]. The same structural homologues were also tested with the ElliPro server [68] and the resulting epitope sequences compared with the results from SEPPA 3.0 to ensure that the results were redundant and method independent.