Genomic and proteomic characterization of two strains of Shigella flexneri 2 isolated from infants’ stool samples in Argentina

Abstract

Background

Shigella specie is a globally important intestinal pathogen disseminated all over the world. In this study we analyzed the genome and the proteomic component of two Shigella flexneri 2a clinical isolates, collected from pediatric patients with gastroenteritis of the Northwest region of Argentina (NWA) in two periods of time, with four years of difference. Our goal was to determine putative changes at molecular levels occurred during these four years, that could explain the presence of this Shigella`s serovar as the prevalent pathogen in the population under study.

Results

As previously reported, our findings support the idea of Shigella has a conserved “core” genome, since comparative studies of CI133 and CI172 genomes performed against 80 genomes obtained from the NCBI database, showed that there is a large number of genes shared among all of them. However, we observed that CI133 and CI172 harbors a small number of strain-specific genes, several of them present in mobile genetic elements, supporting the hypothesis that these isolates were established in the population by horizontal acquisition of genes. These differences were also observed at proteomic level, where it was possible to detect the presence of certain secreted proteins in a culture medium that simulates the host environment.

Conclusion

Great similarities were observed between the CI133 and CI172 strains, confirming the high percentage of genes constituting the “core” genome of S. flexneri 2. However, numerous strain specific genes were also determined. The presence of the here identified molecular elements into other strain of our culture collation, is currently used to develop characteristic markers of local pathogens. In addition, the most outstanding result of this study was the first description of a S. flexneri 2 producing Colicin E, as one of the characteristics that allows S. flexneri 2 to persist in the microbial community. These findings could also contribute to clarify the mechanism and the evolution strategy used by this pathogen to specifically colonize, survive, and cause infection within the NWA population.

Background

Shigella is a gram-negative, enteroinvasive and immobile bacillus, established as the etiological agent of the bacillary dysentery commonly known as shigellosis [1, 2]. Shigellosis is an inflammatory bowel syndrome, mainly affecting the colon and caused exclusively by Shigella. Usually under infection, the colon undergoes cell modifications causing the water loss as diarrhea, normally with mucus and blood. This disease only occurs in humans and can be transmitted directly by fecal–oral route, from an infected patient to a healthy person, or directly through contaminated food and/or water. Bacillary dysentery is a highly contagious infection and can be acquired with the ingestion of only 100 microorganisms, representing a very low infectious dose [3,4,5,6,7].

The genus Shigella (S.) involves four species or subgroups, S. flexneri, S. sonnei, S. dysenteriae and S. boydii with a wide serotype distribution among them. It has been reported that S. boydii covers 20 serotypes, followed by S. dysenteriae with 15 serotypes and S. flexneri with 19 serotypes [8]. However for S. sonnei only one possible serotype has been reported [3]. Worldwide, shigellosis affects annually 164.7 million people, resulting in a high mortality for children younger than 5 years old, especially in low- and middle-income countries [9,10,11]. In Argentina, it has been reported about 1,000,000 cases of infectious diarrhea per year since 2015 [12, 13]. Many reports provide evidence supporting that Shigella species are geographically stratified, distributed according to the level of economic development of a given country [3, 5, 14]. In the western hemisphere, Shigella infections were traditionally mostly travel-related, but recent surveillance data indicate a shift to domestically circulating strains [15,16,17], some of which are increasingly resistant to ciprofloxacin and azithromycin. In addition, different authors have proposed that this pathogen can be used as an index of the hygiene level in a population [7, 18]. However, it is more appropriate to associate Shigella infections incidence with the insufficient condition of hygiene and sanitation. Therefore, whereas shigellosis is caused predominantly by S. flexneri and S. dysenteriae 1 in developing countries, S. sonnei is more frequently associated with outbreaks of shigellosis in industrialized countries [7, 18].

The first genome sequence of S. flexneri was reported in China by Jin et al. 2002. At present, there is a great number of S. flexneri genome sequences reported, currently there are more than 700 sequences uploaded into the data bases of GenBank. Since we previously reported that S. flexneri serotype 2 is endemic and the prevalent pathogen in the NWA region [19], here we performed the whole genome sequencing analysis of two S. flexneri 2 strains isolated from children suffering diarrhea in Argentina. The goal of this work was to determine the main differences at genomic and proteomic levels, including antibiotic resistance aspects, of S. flexneri strains circulating in the NWA region compared with those analyzed in other parts of the world. The acquisition of this knowledge could allow us to identify molecular markers that permit a more rapid and effective pathogen presence detection in patients with gastroenteritis.

Results

General features of the genomes and comparison

Between 2013 and 2017, 2,261 Shigella sp. strains were isolated in the NWA region, where the prevalent serotype identified was Shigella flexneri 2 (S. flexneri 2) [19]. In that work, we determined that S. flexneri serotype 2 displaced serotype AA479 along the last years, highlighting the importance of the emergence of serotype 2 [19]. In order to deeper characterized the prevalent serotype of the region, and to detect molecular markers between them and within isolated strains of other regions of the world, in this work we performed the whole genome sequencing analysis of two S. flexneri 2 clinical isolates (CI). The CI133 and CI172 represent the most characteristic phenotypic groups found in the region [19]. The strains were isolated from hospitals belonging to the two main capital cities of the NWA region. The CI133 strain was isolated during the beginning of this study (2013 year); while CI172 was isolated in the last year of study, 4 years later than CI133 (2017 year).

The primary genome features of CI133 and CI172 strains are summarized in Fig. 1A. The whole genome of CI133 strain is composed of a 4,496,919 bp chromosome defined in 350 contigs, while the CI172 genome consist of 4,465,352 bp, included in 336 contigs. There were determined 4,616 coding sequences and 360 pseudogenes in the sequenced genome of CI133 strain, whereas for CI172 strain, 4,566 coding sequences and 369 pseudogenes were identified (Fig. 1A). For both strains, 19 tRNA types were also identified (Fig. 1A). As shown in Fig. 1B (red ribbons), numerous regions with great similarity between both sequences (100% of identity) were also determined. In addition, we observed that CI133 and CI172 strains have in common 4,194 genes displaying 80% identity and 80% coverage, while 112 and 74 were strain specific genes, respectively (Fig. 1C, Table S1). Among these specific genes, in the CI133 strain it is worth noting the identification of the mdtL and yeaN genes, which are members of the Major Facilitator Superfamily (MFS) (Table S1, in bold). On the other hand, in CI172 genome it is important to highlight the identification as strain-specific genes those encoding for the production and the immunity system of a colicin-like bacteriocin, as well as the toxin-antitoxin (TA) RelEB system (Table S1, in bold).

Fig. 1

Comparison of the newly sequenced S. flexneri 2 strains. A General features of the S. flexneri 2 CI133 and CI172 genomes of sequences. B Comparison between CI133 and CI172 strains genome sequences. The strains are represented by the circumference, while the red ribbons represent the alignments of 100% identity performed by BLAST, the width indicate the alignment length. When the alignment is in opposite orientation of the DNA sequence is indicated by twisted lines in grey colors. C Venn diagram of S. flexneri 2 CI133 and CI172 strains, showing 4,194 shared CDSs at the intersection, the 111 strain specific genes of CI133 represented with an orange circle and the 74 strain specific genes of C172 represented with a green circle

Comparative analysis with other S. flexneri 2 strains

A goal of this work was to determine changes in the genome of the prevalent pathogen from the NWA region to elucidate the main virulence factors that allow us the development of molecular strategies for its eradication, as mentioned earlier. To do this, we analyzed 80 S. flexneri 2 genomes deposited in the GenBank database from NCBI, which were compared with the genome of the two isolates obtained and sequenced in this work. As shown in Fig. 2, the pangenome of S. flexneri 2’ consists in 7,789 families of genes (in blue). From that whole set of genes, 3,182 genes families, represent the “core” genome (in yellow), while the last 4,607 families of genes represent the variable genome (41% and 59% of the genome, respectively). These results show conservation among the genomes of the 82 strains of S. flexneri 2 studied (41%). The number of strain specific genes go from no genes in strains 4028STDY6330022, 4028STDY6330041, 2012ZH074 and 2012ZH118 of S. flexneri 2 to 293 genes in S. flexneri 2 strain 1508 (Table S2, CI133 and CI172 are highlighted in bold). In this analysis, we observed that the CI133 strain have 11 strain-specific genes, while in CI172 we identified 18 strain-specific genes. In addition, these specific genes codify for unknown function proteins (Table 1). Interestingly, among the 80 genomes studied we found that the CI172 strain has as unique genes the cluster encoding for Colicin E, consisting in the production and immunity to the bacteriocin, as well as the lytic protein (Table 1). Moreover, when a phylogenetic analysis of 15 molecular markers was carried out as described in Material and Methods, we observed that the CI172 strain was located in a separate branch, while the other remaining strains were group together as a sister group in which is included the S. flexneri 2a ATCC 29,903 type strain (Fig. S1). Particularly, the CI133 strain was clustered with 4028STDY6330022 and 4028STDY6330016, both isolated in France in 2010 (Fig. S1).

Fig. 2

Differential analysis of the CI133 and CI172 genomes sequence. Pangenomic diversity of the newly sequenced strains of S, flexneri 2 and 80 genomes of NCBI database. Boxplots represented with light blue represent the increasing pangenomic diversity with the addition of the strains whiles the “core” genome is represented by orange boxplots that decreases with the addition of the strains

Table 1 Strain-specific genes of CI133 and CI172 strains

Antimicrobial susceptibility profiles and antibiotic resistance genes identified in both strains

Using the Comprehensive Antibiotic Resistance Database (CARD) we analyzed the presence of genetic antibiotic resistance determinants, which could be used as differential markers. We observed that both strains shared similar profile of antibiotic resistance genes. This analysis allowed us to determine six antibiotic resistance determinants in the sequence of the CI133 strain genome, while in CI172 strain five genes codify for antibiotic resistance. As shown in Table 2A, a perfect match for cat, evgA, sat, tetD and hns resistance genes in both genomes were determined. These genes codify for the resistance to chloramphenicol (cat, phenicol), norfloxacin (evgA, fluoroquinolone), antibiotics targeting nucleosides (sat, such as quinolones), tetracycline (tetD), and carbapenems in general (hns) [20]. However, marA gene was identified only in CI133 genome, pointing out a difference between both strains. In addition, the sequences were analyzed using second software, ResFinder, which allows the identification of functional metagenomic antibiotic resistance determinants. This analysis confirmed the profile and location of the drug resistance genes in CI133 strain (Table 2B). While in the CI172 genome the chloramphenicol resistance gene could not be detected using this tool (Table 2B). To investigate if those genes awarded a resistant phenotype, we performed antibiotic susceptibility tests using antibiotic discs. A set of 10 antibiotics was tested, including the ones which resistance genes were identified in both genomes. As shown in Table 3, both strains shared a similar pattern of antibiotic resistance, showing “in vitro” resistance to chloramphenicol, ampicillin, tetracycline, streptomycin, vancomycin, and imipenem. In addition, CI172 was also resistant to trimethoprim/sulfamethoxazole while CI133 was susceptible, even though the dfrA1 gene which encodes resistance to trimethoprim/sulfamethoxazole was identified by ResFinder in this strain. Therefore, we confirmed “in vitro” the functionality of the genes identified in both genomes by bioinformatics tools.

Table 2 Antibiotic resistance genes identified in CI133 and CI172 strains

Table 3 Phenotype of antibiotic resistance profiles displayed by CII33 and CI172 strains

Pathogenesis and virulence factors

Shigella pathogenesis and virulence were previously well described [21, 22]. It is known that Shigella pathogenicity depends on different chromosomal loci as well as on a mega virulence plasmid harboring the main virulence factors [23, 24]. To investigate the presence of such virulence determinants and the possible differences between both sequences, we used VirulenceFinder database to perform a BLAST screening of the assembled genomes of CI133 and CI172 strains. As a result of this analysis several virulence factors were identified. As shown in Table 4, 10 virulence genes, located in different contigs, were identified in CI133 and 11 in CI172, all of them with 99–100% of identity. The virulence factors identified in both genomes include the transcriptional activator VirF (virF), Shigella IgA-like protease homologue (sigA), Shigella extracellular protein A (sepA), serine protease autotransporters of Enterobacteriaceae (pic), long polar fimbriae (lpfA), invasion plasmid antigen (ipaH9.8), invasion protein (ipaD), glutamate decarboxylase (gad), hexosyltransferase homolog (capUA), and the iron transport protein (sitA). In addition, the presence of celB gene as virulence factor was only determined in the genome of CI172 strain. This last gene codifies for the Endonuclease colicin E2, which is in concordance with the identification of genes encoding for the production and the immunity system to colicin E2 among the strain-specific genes of CI172 (Table 1). The presence of these genes supports and proves the pathogenic phenotype of these strains, which were isolated from patients suffering diarrhea.

Table 4 Virulence genes identified in the CI133 and CI172 strains

Mobile elements

Plasmids

There are several reports describing the plasmids identified by sequencing in different Shigella sp strains [25,26,27,28]. These elements constitute one of the main sources of genetic variation between different bacterial species and within same species. The identification of these plasmids represents an important tool for the construction of molecular markers. An “in silico” search for plasmids was carried out using the PlasmidFinder platform. Through this tool we identified the presence of several contigs related with plasmids in these strains. As shown in Table 5, on both genomes the IncFII plasmidic incompatibility element was detected, suggesting that both strains may harbor the same plasmid. In addition, on CI172 sequence the Col156 and ColRNAI elements were also identified within the contigs 16 and 40, respectively. These findings further highlight the difference between both strains.

Table 5 Plasmidic elements identified in the CI133 and CI172 strains

In order to test if IncFII belongs to the virulence plasmid, the search was extended through contig 263 (CI133) and contig 259 (CI172) alignment with BLASTn, using the complete plasmid databases of NCBI. In both alignments 100% of identity with a specific region of pCP301 plasmid (GenBank accession number NC_004851.1) was obtained. The pCP301 was described by Jin, Q et al. (2002) as a 221,618 bp virulence plasmid of S. flexneri 2a [29]. To uncover the location of the plasmidic loci that potentially contribute to the virulence of the under-study strains, we performed a new alignment using the pCP301 complete sequence and those obtained for the CI133 and CI172 strains. These results allowed us to reconstitute the putative pCP301 plasmid harbored in CI133 and CI172. As shown in Fig. 3A, 16 contigs of the CI133 strain sequence showed homology with different regions of the pCP301 plasmid sequence (99–100% identity), including the contig 263 which harbors the IncFII element (blue box). We also observed that 8 of these 16 contigs shared homology with virulence factors, 5 of which were also identified with VirulenceFinder software (Fig. 3A red arrows, Table 4). In addition, 3 contigs displayed homology to pCP301 regions related to the stability, replication and incompatibility functions, and the last 5 contigs with different insertion sequences located on this plasmid (Fig. 3A, green arrows, blue box and black box, respectively). In the same way, the alignment of pCP301 with the CI172 contigs showed a similar profile to CI133, since 21 contigs were involved, including the node 259 harboring the IncFII element (Fig. 3A, blue box). We observed that in CI172 strain, 13 contigs showed homology with regions of pCP301 containing virulence genes (Fig. 3A red arrows). Here, it is important to mention that these virulence gene regions involved 5 of those 13 contigs (Table 4). The remaining 3 contigs presented homologous sequences with regions inherent to plasmid maintenance and replication, and 5 contigs with insertion sequences (Fig. 3A, green arrows and black box, respectively). Together, these results provide evidence of pCP301 presence as the main plasmid that determines the CI133 and CI172 strains’ virulence.

Fig. 3

Plasmidic content analysis from genome sequence of the CI133 and CI172 strains. A Scale reconstruction of the virulence plasmid identified in the CI133 and CI172 strains, by homology of the sequences obtained with pCP301 and using the DNAMAN program. Red arrows indicate virulence genes identified in each strain; Green arrows indicate elements required for plasmid stability, black arrows indicate insertion sequences and the blue box indicates the origin of replication. The direction of the arrows indicates the orientation of each element with respect to the pCP301 sequence. The node or contigs where each element is found in each strain is indicated on each element. B Strain specific cluster of S. flexneri 2 CI172 strain showing the Colicin related genes and their probable plasmidic origin. C DNA plasmidic profile obtained from virulent CI133 and CI172 strains: 1-λDNA-HindIII molecular marker used to estimate the size of each plasmidic band, the size of each band is in the left indicated; 2- CI133 strain and 3- CI172 strain. The analysis was carried out by agarose gel electrophoresis and stained with ethidium bromide. DNA_ch: chromosomal DNA; p-I, pII and p-III putative plasmids identified by PlasmidFinder software

On the other hand, we analyzed the CI172’s contig 16 and 40 in order to characterize those other putative plasmids identified by PlasmidFinder platform. As showed in Table 5, the contig 16 displayed a 96% identity with Col156, and coincidentally in this contig the presence of those genes encoding the colicin E2 system was evidenced. (Table 4). As with virulence genes, BLASTn alignment was performed using the node 16 sequence to identify the presence of a plasmid. The results showed homology with ColE2-P9 plasmid. As was previously reported, ColE2-P9 is a plasmid of about 7 kb identified in about 10 to 15 copies per Shigella chromosome [30,31,32,33], encoding for the colicin E2 determinant genes. The alignment of pColE2-P9 and Node 16 allowed us to determine that, in this node, the colicin neighbors’ genes displayed homology with plasmid related genes encoding for replication protein and those mobilization/conjugation transfer proteins (MobA/MobL family), suggesting the plasmidic origin of the bacteriocin system (Fig. 3B).

Then, to follow up with the plasmid discovering, we analyzed the contig 40 of the CI172 strain. In addition to the ColRNAI sequence, BLASTn showed 100% identity between the contig 40 sequence and the pSF301-3 plasmid (NC_019250.1). pSF301-3 is a 4,043 kb plasmid isolated from a S. flexneri 2a strain, whose replication is regulated by the Rop protein, which acts through modulation of primer RNA precursor transcription [34,35,36]. Like pColE2-P9, pSF301-3 harbors genes encoding for a mobilization family of proteins such as MobC, MbeB and MbeD, necessary for the specific transference when other conjugative plasmid it is also present into the host strain [37]. The identity observed with this plasmid suggests that the contig 40 is a small plasmid also present into CI172, representing another difference with CI133.

Finally, to experimentally confirm the presence of the mentioned plasmids in the CI133 and CI172 strains we performed the procedure described by Kado and Liu (1981) for identification and isolation of different size plasmids including large and small elements [38]. As shown in Fig. 3C, when the plasmid content of CI133 and CI172 was analyzed and compared against the λDNA-HindIII molecular marker, we were able to infer the number and size of the putative plasmids identified “in silico” in both strains. In agreement with the bioinformatics results, the CI133 strain presented only one band, possibly corresponding to the virulence plasmid of around 22,16 kb (Fig. 3C, arrow). While in the CI172 strain the presence of four bands was observed, three of them could correspond to the 3 identified plasmids, the virulence plasmid (30 kb), the pColE (4 kb) and the pSF301-3 (4 kb) (Fig. 3 C, pI, p-II and p-III, respectively).

Prophages

The presence of prophages within the chromosome can allow some bacteria to develop resistance to antibiotics, adapt to new environmental niches, even improve adherence to host and become pathogens. Using the PHAST tool [39], a search for phages in the genomes of CI133 and CI172 strains was carried out. The results showed that in CI133 genome, 4 prophage regions are present as intact elements (Table 6), while other 7 regions are incomplete and 2 regions as questionable elements, making a total of 13 prophage regions identified (Table S3). Furthermore, in the 4 intact prophage regions the left and right att phage binding sites to the bacterial chromosome were detected, as well as those genes encoding for the integrase, tail, head, capsid, portal and the terminase (Table S4). The location of these regions in the CI133 obtained sequences, allowed to determine that all the contigs involved in the phage regions are continuously overlapped within each other in the genome (Table 6). In the case of CI172, there were found a total of 12 phage regions, 2 of which are questionable, 7 correspond to incomplete prophages (Table S3) and 3 correspond to intact prophages. In these last 3 regions, all the elements that define a phage, contained in different continuous contigs, were also detected (Table 6, Tables S3 and S4).

Table 6 Intact prophage regions identified by PHAST software into the CI133 and CI172 genomes

In order to confirm the above result, we also analyzed the presence of prophage regions in both strains using PHASTER tool (PHAge Search Tool Enhanced Release). As shown in Table 7, using this new server it was possible to detect 4 prophage regions into the CI133 genome, but just only the region 4 harbor an intact element, while other 1 region is incomplete and 2 regions more represent questionable elements. It is important to mention that the intact prophage identified into the region 4, involving the node 458 sequence (Table 7), was not detected by PHAST, confirming the efficiency of update performed on PHASTER server. Furthermore, in the node 458 sequence all the elements required for a functional phage like the attL and attR phage binding sites, the integrase, tail, head, capsid, portal and the terminase encoding genes, were detected (Table S4B). Unlike PHAST, the results obtained from PHASTER did not show the presence of intact prophage regions into the CI172 genome, rather than 2 of 4 regions were questionable prophages and 2 more corresponded to incomplete prophages (Table 7).

Table 7 Prophage regions identified by PHASTER software into the CI133 and CI172 genomes

Interestingly, when we compared the results of both servers, we observed that in both genomes, the regions of prophages that were questionable with PHASTER involved nodes that were forming part of complete phages identified with PHAST (Table 6 vs Table 7). This led us to hypothesize that these differences were observed since the sequences that we used for the analysis are separated in nodes or contigs and not assembled. To confirm this assumption, the sequence of each node involving one complete phage by PHAST detected, were assembled and uploaded into the PHASTER server. This time, PHASTER identified the same complete phage identified by PHAST. These results indicate that despite of differences between both servers, the results obtained by each are complementary. Therefore, we demonstrate that the AC133 strain has a total of five regions that correspond to five possible complete prophages, while the AC172 strain contains only three.

CRISPR systems

CRISPR (Clustered Regularly Interspaced Short Palindromic Repeats) constitute a particular family of repeated sequences in tandem found in a wide range of prokaryotic genomes. The CRISPR consist of a nucleotide sequences located within highly conserved regions that vary in size, are separated by unique sequences of similar size (spacer) and they have generally a viral origin [40]. The systems called CRISPR-Cas (CRISPR-associated Cas proteins) provide the bacteria immunity against invasive genetic elements, such as plasmids and phages, and provides a form of acquired immunity. Therefore, using the online tool CRISPR Finder was possible to search for these CRISPR systems in the sequenced genome of both strains under study. As shown in Table 8, within the CI133 genome 4 questionable CRISPR systems were identified into different contigs, while in CI172 genome 3 of these questionable elements were found. Moreover, when CRISPR-Cas Meta Finder server was applied to the genomes analysis we identified one more of these systems in CI133 and CI172 strains. These elements were located into the Node 125 and 47, respectively (Table 8). Unlike those CRISPR above described these new systems could be functional, since not only displayed sequence homology with the CRISPR Class I family, but also contains the cas3_1_I genes and the DR consensus (Table 8). In addition, we observed that the Class I-CRISPR is apparently the same in both strains due to presence of equal DR length and number of spacers. It is important to mention that a more in-depth analysis of these CRISPR immune systems is being carried out in our laboratory.

Table 8 CRISPR/Cas systems identified into the CI133 and CI172 genomes

Proteomic characterization

In concordance with the genomic analysis, a similarity was observed between the proteomes of both strains, 4,540 proteins were identified and quantified by LC/MS–MS using Label Free Quantification (LFG). As shown in Fig. 4, the Clusters of Orthologous Groups of proteins (COGs) assignment did not show great differences between the CI133 and CI172 strains [41]. The COG categories more represented in both strains were the carbohydrate transport and metabolism (G), the amino-acids transport and metabolism (E), those related to poorly characterized proteins with a general function prediction (R) and of unknown function (S) (Fig. 4).

Fig. 4

CDSs corresponding to the COG functional categories represented in the CI133 and CI172 strains genomes. The letters represent each one of the functional categories: [D] Cell cycle control, cell division, chromosome partitioning; [M] Cell wall/membrane/envelope biogenesis; [N] Cell motility; [O] Post-translational modification, protein turnover, and chaperones; [T] Signal transduction mechanisms; [U] Intracellular trafficking, secretion, and vesicular transport; [V] Defense mechanisms; [W] Extracellular structures; [A] RNA processing and modification; [J] Translation, ribosomal structure and biogenesis; [K] Transcription; [L] Replication, recombination and repair; [C] Energy production and conversion; [E] Amino acid transport and metabolism; [F] Nucleotide transport and metabolism; [G] Carbohydrate transport and metabolism; [H] Coenzyme transport and metabolism; [I] Lipid transport and metabolism; [P] Inorganic ion transport and metabolism; [Q] Secondary metabolites biosynthesis, transport, and catabolism; [R] General function prediction only; [S] Function unknown

However, when the proteomic results determined in a defined medium simulating the host`s environment were more deeply analyzed, we observed the presence of specific proteins in one or another of the strains under study. Thus, 16 of the 4,540 proteins were specifically expressed in CI133 but not detected in the CI172 strain (Table 9). In the same way, 12 of the total determined proteins were only expressed in the CI172 strain (Table 9).

Table 9 Experimental proteins detection in CI133 or CI172 strains

Moreover, significant differences (p < 0.05) in the expression of a total of 33 proteins were also detected, 26 of which were 1 to threefold upregulated in CI133 respect to CI172 (Table 10, positive values), while the last 7 proteins were overexpressed in CI172 respect to CI133 (Table 10, negative values).

Table 10 Differential proteins expression between CI133 and CI172 strains

To understand the relation and the involvement of these proteins with a particular process or function we constructed protein networks using STRING v11.0 with the proteins expressed only in CI133 or CI172 plus the upregulated proteins in each strain (Tables 9 and 10). For S. flexneri 2 CI133 strain, the 42 proteins, taking in account (16 only present and 26 upregulated in this strain), represented by the nodes showed 37 edges or interactions (Fig. 5A). The number of nodes determined was significantly higher than the calculated by STRING for a random set of proteins, suggesting that these proteins could have a biological connection, as a group. The main represented functional categories were the metabolism and transport of amino-acids (E), energy production and conversion (C), intracellular trafficking, secretion, and vesicular transport (U) and cell wall/membrane/envelope biogenesis (M) (Fig. 5A). Meanwhile, for the CI172 strain, were only found 7 edges among the 19 proteins (12 only present and 7 upregulated, respectively) (Fig. 5B). The most represented category was coenzyme transport and metabolism (H). However, the proteins included in category translation, ribosomal structure and biogenesis (J) are the ones that showed the greatest number of interactions between themselves and between proteins from other categories (Fig. 5B).

Fig. 5

Protein–Protein interactive networks. A Proteins present and upregulated in the CI133 strain. B Proteins showing interactions present in the CI172 strain. The proteins are represented with the nodes, the edges represent the interactions between the proteins, and the width of edges the degree of interactions. The circles are showing different COG categories: [O] Post-translational modification, protein turnover, and chaperones; [U] Intracellular trafficking, secretion, and vesicular transport; [C] Energy production and conversion; [E] Amino acid transport and metabolism; [M] Cell wall/membrane/envelope biogenesis; [I] Lipid transport and metabolism; [J] Translation, ribosomal structure and biogenesis; [H] Coenzyme transport and metabolism and [G] Carbohydrate transport and metabolism. The networks were constructed using the STRING v11.5 bioinformatics tool

These set of proteins that were found either overexpressed or only present in one of the two strains (Tables 9 and 10) were most of them functionally associated to cellular metabolism. However, between the proteins that were overexpressed in CI133 we identified the virulence regulon, including the transcriptional activator VirB, the VirA protease, the SecA translocase and the MxiA protein. VirB, encoded in Shigella’s virulence plasmid, is a master regulator of a great number of genes including those required for synthesis and formation of the Type Tree Secretion System (T3SS) containing mxiA [42,43,44]. MxiA, is an essential membrane protein involved in the secretion process of T3SS’s substrates [45, 46]. VirA protein is a virulence effector, translocated into the host by the T3SS, to produces the host cell membrane ruffling necessary for the pathogen entry and for intra and intercellular spreading [47]. The SecA translocase subunit protein is a part of the secretion system SecAB also involved in Shigella’s virulence [48].

Therefore, despite we were able to see similarities between the proteomes expressed in both strains, we were able to determine differential expression of 33 proteins. According to the differences found in the genomes and in the differential expression of proteins, we can suggest that each of the analyzed strains can be adapted differently to the host environment to cause infection.

Discussion

Previous report demonstrated that S. flexneri is the endemic pathogen in most developing countries, including Argentina, being the serovars 1b, 2a, 3a, 4a and 6 predominant, while S. flexneri serovar 2a, as S. sonnei, is frequently isolates in industrialized countries [49]. In accordance, in our previous work we have demonstrated that S. flexneri 2a is the gastroenteric infectious agent with the highest incidence in the NWA region, representing the main cause of children’s morbidity and mortality [19]. However, we also demonstrated that the prevalence of this serotype was observed in the last 7 years, since in the period from at least 2007 to 2013 the S. flexneri AA479 serovar cause the most gastroenteric disease in the NWA population [19]. Due to the continuous turnover of serovars observed, in this work we performed the genomic and proteomic analysis of two clinical isolates obtained in different periods of time, with a four-year of difference between one from the other. Our objective was to detect molecular changes that may have occurred during this period, which could help us understand or explain the potential establishment and persistence of this pathogen in the population under study.

The data obtained allowed us to determine the pangenome, the “core” genome and the variable genome of the 80 available strains as well as our CI133 and CI172 clinical isolates. We observed that the “core” genome here described represents most of the analyzed S. flexneri 2 genomes, and it is relatively conserved. These results are in agreement with data already published by numerous authors, supporting the idea that Shigella has a conserved “core” genome in which the virulence determinants (that allow it to colonize, persist and spread in the region) are perhaps encoded [50,51,52].

When the variable genome was analyzed, we observed that the isolates from the NWA region presented a large number of unique genes, suggesting that this serovar evolved in recent years through the acquisition of these strain-specific genetic elements. The results obtained in this work are also in agreement with those reported by Shen et al., who analyzed the pangenome of 15 S. flexneri isolates and identified a number of specific genes comparable to those obtained in this study [52]. It is important to note that, even though the CI172 strain showed a lower number of strain-specific genes, we observed a greater number of these genes with known function than those obtained in the CI133 strain (24 of 74 vs. 29 of 112, respectively). Among the specific genes of the CI133 strain, the presence of mdtL and yeaN is noteworthy. These genes are both part of the Major Facilitator Superfamily (MFS) composed by thousands of transport proteins that are well characterized [53]. MFS are involved in the transport of a great variety of small molecules such as sugars, amino acids, and drugs [54]. Interestingly, some of them are specifically related to Multi Drug Resistance (MDR), known as MFS-MDR transporters [55, 56]. MFS-MDR can carry out an active flux of drugs, antibiotics, or toxins through a cation/substrate antiport mechanism generating resistance to them. Furthermore, there are some reports indicating that the yeaN gene product could be acting as a 2-nitroimidazol transporter conferring E. coli resistance to this drug [57]. On the other hand, among the specific genes identified in CI172 genome, we were able to determine the presence of the complete set of genes that codifies for the production, immunity and lysis protein of the colicin E2 and Rel toxin-antitoxin systems. Colicins are high molecular weight toxic proteins that are produced by colicinogenic E. coli strains and some related species of the Enterobacteriaceae family [58,59,60,61], being at least 23 colicin types have been described in detail [62]. Colicin production has been well described among Shigella genera, especially for the subspecies S. sonnei [63]. Interestingly, from the total of 82 genomes analyzed in this study Colicin E codification system was only found in the CI172 strain isolated in this work (Table S2). Besides, we identified the neighbor genes as plasmid related genes suggesting a plasmidic origin or that CI172 could have acquired this bacteriocin through plasmid transference. Although MvpAT is the most commonly described antitoxin system (TA) in S. flexneri, we were unable to identify these genes in either genome, instead RelEB-TA was identified among the CI172-specific genes. The RelEB-TA systems usually consist of two genes that encode a stable toxin and an unstable antitoxin that inhibits the toxin [64]. RelE toxin is known to be involved in the emergence of persistent strains, resistant to multiple drugs, biofilm formation, and regulation of the bacterial population by inducing dormancy [65,66,67]. This TA system is present in the virulence plasmid and its main function is to ensure its stability in bacteria [68, 69]. However, in the CI172 strain this TA system was located within the NODE_575_length_1729, which was not associated with plasmid regions. The homology search of the complete node-575 sequence by BLASTn, it was observed that the TA system coding region is flanked upstream by the IS600 element (NP_707432.1) and downstream by the gene coding for SF1546 a bacteriophage protein (NP_707877.2). These data suggest that the CI172 strain acquired the TA system into its genome by mobile genetic elements such as insertion sequences and phages. This could explain the multiple antibiotic resistances phenotype of CI172 strain determined “in vitro”. Therefore, we suggest that the presence of these genes could be addressed as part of the acquired adaptive advantages allowing the persistence of S. flexneri serotype 2 in the NWA region. Together the pangenomic, phylogenetics strain-specific genes and plasmid profile analysis show that strain CI172 is quite singular, unique and different from the rest of the studied strains. Furthermore, the phylogenetic analysis reveals that both strains studied in this work, isolated in Argentina, are closely related to those isolated in Europe than to those from other continents.

Multi-drug resistance has been highly addressed in the last decades, antibiotics resistance spread is a worldwide worrisome. Most of Shigella’s serotypes are not susceptible to first-generation antibiotics, such as ampicillin, tetracycline, sulfonamides, chloramphenicol, nalidixic acid, and Trimethoprim/sulfamethoxazole. Therefore, last-generation antibiotics, such as fluoroquinolones are preferred. However, these are more expensive and resistance to them are starting to show up in the countries with more consume [70, 71]. The “in silico” analysis of the two genome sequences using CARD and ResFinder databases allowed us to determine the presence of several genes that encode antibiotic resistance which was then experimentally confirmed, demonstrating that the strains under study are presented with an MDR phenotype. Despite this, both strains were differentiated mainly by the presence of the marA gene, which was identified in CI133 but not in CI172. The marA gene was previously reported to provide multiple drug resistance [72]. Has been proposed that the transcription factor MarA determines a multidrug resistance phenotypes, allowing for the resistance of quinolones, tetracyclines and fluoroquinolones by modulating efflux pump and porin expression [72]. The presence of this gene in CI133 suggested a broad spectrum of resistance; however CI172 was additionally resistant to trimethoprim/sulfamethoxazole. At the beginning, shigellosis was used to be treated with antibiotics such as sulphonamides, tetracycline, and chloramphenicol [73]. The emerging resistance to those antibiotics enforced the use of additional drugs such as ampicillin, co-trimoxazole, nalidixic acid, and fluoroquinolones to combat the infection. Chloramphenicol, ampicillin, and tetracycline resistances that were found in CI133 and CI172 are in concordance with the resistance phenotypes already reported [50]. However, vancomycin and imipenem resistances present in both strains have not been between the most common antibiotic resistance phenotypes among Shigella strains. The family of carbapenemases enzymes, by plasmid encoded, is one of the most significant health challenges of the century, given the potential for dissemination between species and mortality rates due to infections caused by bacteria. The emerging resistance to carbapenems represents a serious concern since these drugs are often the last line of effective therapy available for the treatment of infections caused by MDR bacteria [74,75,76]. However, the pathogens persistence is not only result of MDR, but also of the adaptative advantages acquired and its virulence performance.

Using VirulenceFinder we identified 10 of the main genes involve in Shigella’s pathogenesis in different contigs of the CI133 strain and 11 on the CI172 strain, displaying a new difference between these two isolates. In addition, when the homology analysis against plasmid pCP301 was performed we identified other virulence factors not detected by the program. The results obtained allowed us to determine those genes that may be forming part of the bacterial chromosome, while others are forming part of the virulence plasmid pCP301 [29]. Among the genes identified by the mentioned software, it is noteworthy that sigA, pic, lpfA, gad and sitA would be located on both genomes. In concordance with our localization, Henderson et al. (1999) reported that in enteroaggregative E. coli (EAEC) 042 and S. flexneri 2457 T the pic gene is encoded within the chromosomes [77]. Interestingly, in CI133 and CI172, pic and sigA genes were located in same contigs, thus being able to affirm the chromosomal location of both. Moreover, the products of these genes are member of the subfamily of autotransporter proteins termed the SPATEs (serine protease autotransporters) of Enterobacteriaceae [77, 78]. While, the lpfA, gad and sitA genes products are involved in the major fimbriae formation, the controls of bacterial acid resistance-SlyA dependent and in an iron-ABC transport system during the eukaryotic cells infection, respectively [79,80,81].

Numerous mobile elements can be found in the variable genome of a strain, which are the result of evolution by selection of genes acquired by horizontal transfer (THG). The determination of the mobilome (set of mobile elements) in a genome allows characterizing the sequenced genomes. The mobilome can be constituted by integrative and conjugative elements, prophages, and plasmids. The CRISPRs systems and bacteriocin genes also contribute to the diversity of the accessory genome. The use of services available in the Center for Genomics Epidemiology and in different web sites, as well the combination of the results obtained from each one, we were able to identify the mobilome of the strain under study. As above mentioned, we could determine that numerous virulence genes are contained in pCP301, a plasmid previously described as essential for establishing infection in the host. This plasmid is a member of the IncFII incompatibility group, element also present in the sequence of CI133 and CI172 [29]. In addition to the virulence genes and the IncFII element, we were able to identify the contigs of the analyzed sequences harboring those genes that code for the replication control, stability and mobilization of the plasmid, as well as several insertion sequences that showed homology with pCP301 [29]. This allowed us to conclude that this plasmid is responsible for the establishment of the infection of CI133 and CI172 pathogenic strains.

On the other hand, detection of Col156 and ColRNAI into CI172 sequence showed two more major differences with the CI133 strain. Col156 is the genetic element that characterizes plasmids harboring the synthesis, immunity and lysis system of colicin E, which was identified in contig 16 of CI172 as well as the ceaB synthesis genes, the immunity ceiB and the lysine celB of the colicin E2 system [30]. The homologues search of this contigs sequence showed 99% identity with pColE2-P9 of about 7 Kb of Shigella, which occurs in low number of copies per chromosome [31,32,33]. These data correlate with our experimental results, where we were able to detect a band that migrate approximately at that size and that stains poorly with ethidium bromide. Similarly, “in silico” and experimentally, we demonstrate the presence of pSF301-3 in the CI172 strain, more precisely in contig 40. Like pColE2-P9, pSF301-3 has the ability to be mobilized when a conjugative plasmid is harbored within the same strain [37], indicating that both strains in this study could have evolved in a period of four years through the acquisition of these mobile elements.

The presence of prophages on the chromosome can allow some bacteria to acquire resistance to antibiotics, exist in new environmental niches, even improve adhesion or become pathogens. Thus, the numerous phage regions found in the CI133 and CI172 strains could be related to the MDR phenotype presented by these strains, although the correlation between both factors analyzed could not be precisely established. This suggestion is supported by the numerous reports describing the joint detection of resistance to TMS and ampicillin related to integrons present in phages [82]. Due to the presence of phages, in these strains the presence of CRSIPR elements was also analyzed, since these represent a type of immune system of the bacterium to resist the attack of phages [83, 84]. As expected, numerous regions were found that could be part of these elements in both strains, highlighting the importance of the CRISPR-Cas systems, where the Cas proteins involved in the mechanism of immunity against phages could be identified. Even though the relationship between prophages and CRISPR has not been determined, numerous experimental studies are being carried out in our laboratory to determine its importance in the evolution of the prevalent strains in the region.

To complete the genomic study, the secretome of both strains was also analyzed in a medium simulating the environment of the host. In total, a great difference was not observed between both strains in the number of proteins that were synthesized and exported to the culture medium or when the Clusters of Orthologous Groups of proteins was determined. In both proteomes the great majority of the proteins identified by LFQ were the carbohydrate transport and metabolism, the amino acid transport and metabolism, the transcription and the cell wall/membrane/envelope biogenesis. Interestingly, despite the similarities observed, it was possible to detect proteins differentially present in one or another strain, as well as proteins that were overexpressed in one isolate with respect to the other. In addition, in the CI133 strain a greater number of proteins capable of interacting were observed respect to the interactions obtained in the CI172 strain. This observation could be due to the greater number of proteins used in the analysis. However, these differences suggest that both strains have undergone different evolutionary processes for host adaptation and infection. We here suggest that this evolutionary process is more efficient in the CI172 strain since a lower number of secreted proteins would have a lower energy cost for the bacteria but would not affect its pathogenesis.

Conclusion

In summary, even when great similarities were observed between the CI133 and CI172 strains and those other from different countries, confirming the high percentage of genes constituting the core genome of S. flexneri 2. In this work a considerable number of strain specific genes were identified. On the bases of these observations, we can address the importance of characterizing rapid spread of pathogenic bacteria in specific region. We here hypothesize that the identification of these highly identical or strains specific genomic and/or proteomic factors into other isolated of our collation allow us to develop molecular markers that are characteristic of the NWA region circulating pathogens. In addition, the most innovative result of this work is, to our knowledge, the first description of a producer of S. flexneri 2 Colicin E, a factor recently acquired by mobile elements (horizontal gene transfer), as one of the characteristics that allow S. flexneri 2 persist in the microbial community. It is important to highlight that, like the CI172 strain, at least 350 isolates of Shigella producing antimicrobial compounds have been also previously identified in our laboratory. This supports our hypothesis that the production of these compounds has been widely distributed among the circulating pathogens in the region under study.

Methods

Strain information

The clinical isolates of S. flexneri 2 (CI133 and CI172), were isolated from female children patient with severe acute clinical manifestations of shigellosis in the NWA region. The CI133 and CI172 strains were isolated from children patients suffering diarrhea [19]. The strains were routinely grown at 37 °C overnight on LB to stationary phase. These strains were identified and characterized in a previous work [19].

Genome sequencing and assembly

DNA extraction and WGS sequencing of the CI133 and CI172 strains were performed by Macrogen Humanizing Genomics (Seoul, Korea), using the Illumina HiSeq2500 platform of paired-end libraries, taking as reference the sequence of strain S. flexneri 2a str. 301 from the NCBI database (access number NC_004337). The quality of the raw data obtained from the sequencing of these genomes was analyzed with the FastQC program (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The assembly was made through Velvet 1.2.1, using the best k-mer value with a value of k = 81 for the CI133 strain, and k = 79 for the CI172 strain [85]. Genome annotation was performed according to standard procedures described by Prokaryotic Genome Annotation Pipeline (PGAP) from NCBI (https://www.ncbi.nlm.nih.gov/genome/annotation_prok/). The pangenomic analysis of the 80 strains of Shigella flexneri 2 plus the two newly sequenced strains (CI133 and CI172) were performed through the Microscope Platform with an identity of 80% and a coverage of 80% [86].

Data availability

The Whole Genome Shotgun project of Shigella flexneri 2 CI133 strain has been deposited at DDBJ/ENA/GenBank under the accession JAGDQG000000000, the Bioproject accession number is PRJNA498014 and the Biosample accession number is SAMN10275217. While the Whole Genome Shotgun project of Shigella flexneri 2 CI172 strain has been deposited at DDBJ/ENA/GenBank under the accession JAGDQF000000000, the Bioproject accession number is PRJNA498020 and BioSample accession number is SAMN10275218. The versions described in this paper are the versions JAGDQG010000000 and JAGDQF010000000 for the CI133 and CI172 strains, respectively.

Phylogenetic analysis

A multi locus sequence typing (MLST) was performed using fifteen housekeeping genes, including arcA (DNA-binding transcriptional dual regulator), aroE (shikimate dehydrogenase), aspC (aspartate aminotransferase), clpX (ATP-dependent Clp protease ATP-binding subunit), cyaA (adenylate cyclase), dNaG (DNA primase), fadD (fatty acyl-CoA synthetase), grpE (nucleotide exchange factor), icdA (isocitrate dehydrogenase), lysP (lysine:H( +) symporter), mdh (malate dehydrogenase), mtlD (mannitol-1-phosphate 5-dehydrogenase), mutS (DNA mismatch repair), rpoS (RNA polymerase, sigma 38 factor) and uidA (beta-glucuronidase), as previously described [87]. The concatenated sequences were align using MUSCLE (https://doi.org/10.1093/nar/gkh340) and the phylogenetic tree was inferred with UPGMA method and constructed with MEGAX [88].

Antibiotic resistance genes, virulence factors and mobile elements prediction

Antibiotic resistance genes were identified using The Comprehensive Antibiotic Resistance Database CARD (http://arpcard.mcmaster.ca) [89]. The presence of genetic determinants of drug resistance and mobile elements in the complete sequence of CI133 or CI172 was also investigated using those serves of the Center for Genomics Epidemiology. For antibiotic resistance genes identification we used ResFinder-4.1 server (https://cge.cbs.dtu.dk/services/ResFinder/) [90]. As selection criteria only the genes with perfect hit and 99–100% of identity were considered. The virulence factors were identified by VirulenceFinder (https://cge.cbs.dtu.dk/services/VirulenceFinder/), using as criteria of selection 99–100% of identity [91]. The presence of plasmids in the strains under study was determined “in silico” by PlasmidFinder platform (https://cge.cbs.dtu.dk/services/PlasmidFinder/) [92]. The presence of prophages within both bacterial chromosomes was determined by PHAST (http://phast.wishartlab.com/) [39], and PHASTER tools (http://phaster.ca/) [39, 93]. As it is well explained on the website, this server is an update of PHAST, also used for the identification and annotation of prophage sequences present in both the bacterial genome and plasmids. However, both differs in that the PHASTER server being faster since more than 120 CPUs were added to the computational cluster; is more efficient because it has a series of automated algorithms that reduce the size and redundancy of databases; and it is more visually attractive and easier to use, because it has a new graphical genome browser and better tools for visualizing genes and genomes whit a robust and easy-to-use graphical interface. Finally, the search for CRISPR systems was carried out with the CRISPRfinder platform (http://crispr.upsud.fr/Server/CRISPRfinder.php) [40]; while the presence of CRISPR-Cas systems in the genomes were analyzed using the server CRISPR-Cas Meta (https://crisprcas.i2bc.paris-saclay.fr/CrisprCasMeta/Index) [94]

Antimicrobial susceptibility and plasmid profiles

Antimicrobial susceptibility tests were performed using Müller–Hinton agar (Oxoid, Hampshire, UK) by disk diffusion method. Resistance and sensitivity were interpreted following the standards criteria of the Clinical & Laboratory Standards Institute (CLSI) [95]. The following antimicrobial agents were used: chloramphenicol (30 µg), ampicillin (10 µg), trimethoprim sulfamethoxazole (1.25/23.75 µg), fosfomycin (200 µg), tetracycline (10 µg), kanamycin (30 µg), streptomycin (10 µg), vancomycin (30 µg), nalidixic acid (30 µg), imipenem (10 µg) and gentamycin (10 µg).

The plasmid extraction was performed by the alkaline lysis method as described Kado and Liu (1981) [38]. The plasmids were analyzed by electrophoresis on 0.7% agarose gels, stained with ethidium bromide and UV light exposition. Plasmid molecular size estimation was conducted by calibration curves using the λDNA-HindIII molecular marker (Promega).

Proteomic analysis

An overnight culture was used to inoculate 50 mL of M9 medium. The CI133 and CI172 strains were grown at 37 °C under static condition until stationary phase. These bacterial cultures were carried out by triplicate. Cells were harvested by centrifugation and resuspended in 5 mL of purification buffer (300 mM NaCl, 30 mM NaH2PO4, pH 8.3) and then were disrupted with a French press. Protein identification was performed by IQUIBICEN/UBA-CONICET, University of Buenos Aires (Buenos Aires, Argentina) through Label-free quantification using Orbitrap. Proteome Discoverer 2.1 and Perseus 1.5.8.5 software were used to analyze the spectrums obtained. Proteins were considered significantly upregulated if they were 2- or more fold change induced for a p-value < 0.05. The interactions between two sets of proteins were determined, and the networks were constructed using the STRING v11.5 bioinformatics tool (https://string-db.org/) [96] and taking in account all the types of interactions with a moderate confidence level (0.4).