The evident increase in MIC of eravacycline upon induction of resistance
Acinetobacter baumannii strains were exposed to sequential passages of increasing concentration of eravacycline for evaluating the acquisition of resistance. The strains with acquired resistance were evaluated for gene expression profiling, and OMV proteome analysis was studied to identify a specific pattern in OMVs pertaining to resistance. The MICs of eravacycline for the A. baumannii ATCC 19606 and JU0126 strain were 0.125 and 0.5 μg/mL, respectively. The serial passage-based induction of acquired resistance among the isolates was carried out from one-eighth MIC values and above, until a lethal concentration was reached above the MIC concentration. A. baumannii ATCC 19606 was able to resist the passages up to concentration of 64× MIC after which it failed to overcome the action of eravacycline. A. baumannii JU0126 strain tolerated upon induction till concentrations of 64× the MIC of eravacycline, above which turned to be lethal concentration. Both organisms were able to present a resistant phenotype below 128× MIC concentration of eravacycline.
High-throughput RNA sequence analysis
RNA sequence analysis was performed for the A. baumannii strains, ATCC 19606 and JU0126, respectively. Both strains were grown under 64× MIC concentrations of eravacycline obtained as per induced resistance protocol described above. From high-throughput RNA sequence analysis, the length of 16205012 (an error probability of 0.03%) and 15029752 (an error probability of 0.02%) clean reads were obtained from eravacycline resistance induced ATCC 19606 and JU0126 strains, respectively. A total of 18732924 (an error probability of 0.03%) and 16680372 (an error probability of 0.03%) clean reads were obtained from control untreated strains of ATCC 19606 and JU0126 strains, respectively. The Q20 of all these four samples reached 98% for treated and 97% for control which indicated a high quality of transcriptome sequencing. The GC content (%) were 44.03 and 43.87 (ATCC 19606 and JU0126, respectively) for treated strains and 45.11 and 45.04 (ATCC 19606 and JU0126, respectively) for control strains. Pearson’s correlation between each sample was analyzed: ATCC 19606 control with ATCC 19606 treated strains, had a value of 0.774 and JU0126 control had 0.735 correlation values with treated JU0126. ATCC 19606 control strains had 0.848 correlation value with treated JU0126, ATCC 19606 control had 0.866 correlation value with control JU0126 and JU0126 control strain had 0.926 correlations with treated ATCC 19606 strains. These correlation values show high correlation between the samples.
Significant DEGs among the eravacycline treated strain when compared with untreated control strains
The complete gene expression values for A. baumannii ATCC 19606 and JU0126 eravacycline treated strains are provided in Additional file 2. In an effort to study the changes in the biological mechanisms and/or pathways of the bacterial system upon resistance to eravacycline in the treated strains when compared with the control eravacycline susceptible strains of A. baumannii ATCC 19606 and JU0126, DEGs analysis was performed (Additional file 3; Fig. 1A, B). For DEGs analysis, parameters P values (<0.05) and fold changes ≥2 were used. A total of 1404 genes were identified as significant DEGs, of which 944 DEGs (67.2%) were upregulated and 460 DEGs (32.8%) were downregulated in A. baumannii ATCC 19606 treated sample when compared with control samples. From A. baumannii JU0126 strain, 1285 genes were identified as significant DEGs between each sample, of which 574 DEGs (44.7%) were upregulated and 711 DEGs (55.4%) were downregulated in A. baumannii JU0126 eravacycline treated strain when compared with JU0126 control strain.
GO enrichment analysis of DEGs
GO enrichment is widely used to find the biological roles of each gene and its products [16]. All DEGs were mapped to their terms in GO database and compared with the reference transcriptome. GO mapped DEGs from ATCC 19606 and JU0126 were identified and classified into functional groups in three main categories: biological process, cellular process and molecular function (Fig. 1C, D). Totally, 2219 GO terms were identified in the DEGs from ATCC 19606 under all three categories; with localization, transport, cellular component, membrane, transport activity, and transmembrane transport activity being dominant terms. In JU0126, organonitrogen compound metabolism, protein metabolism, protein-containing complex, cytoplasm, structural activity and structural constituent of ribosomal were dominant terms in all three categories. Significantly enriched GO terms were considered based on the corrected P < 0.05.
Kyoto encyclopedia of genes and genomes (KEGG) analysis of DEGs
KEGG database is a collection of various pathways, which represents the molecular interactions network between each gene/proteins [17]. To identify the enriched pathways involved in eravacycline induced strains of A. baumannii ATCC 19606 and clinical strain JU0126, KEGG analysis was done. In total, 78 and 86 pathways were identified in the DEGs of ATCC 19606 and JU0126 strains, respectively. In both ATCC 19606 and JU0126, microbial metabolism in diverse environments, ABC transporters, degradation of aromatic compounds and fluorobenzoate degradation pathways were the most highly represented categories. The enriched factors were represented in the ratio of the differentially expressed gene number to the total gene number in a certain pathway. The values were represented in Q value, which is a corrected P-value ranging from 0 to 1. The size and color gradient of the dots indicate the range between Q value and the number of DEGs mapped to the indicated pathways, respectively. The top 20 values are shown in Figure 1E and 1F.
Quantitative reverse transcriptase-polymerase chain reaction (qRT-PCR) validation of DEGs
To analyze further, 6 DEGs from both upregulated and downregulated genes for both the strains were selected from the RNA sequence analysis and validated through qRT-PCR study (Fig. 2). Three genes that showed significantly increased (upregulated) expression, as well as three genes that were downregulated in each ATCC 19606 and JU0126 strains treated eravacycline compared with untreated strains, were taken further for qRT-PCR. These 6 DEGs for each strain were selected for qRT-PCR verification. From the ATCC 19606, the genes multidrug efflux RND transporter permease subunit (AUO97_00445), MSF transport (AUO97_00560), M1 family peptidase (AUO97_00700) which upregulated with 3.2731, 1.9644, 1.2859 log2-fold, respectively in RNA-sequencing data, showed 3.8204, 2.5822 and 2.8533 log2-fold changes, respectively using the qRT-PCR analysis. From the RNA-sequencing data, genes porin (AUO97_05635), trifunctional transcriptional regulator (AUO97_15195) and transfer-RNA(AUO97_11755) which downregulated with -2.782, -1.176 and -3.9366 log2-fold, respectively, displayed-1.5859, -2.0788 and -2.0203 log2-fold changes, respectively using the qRT-PCR analysis. From the JU0126 strain, the genes corresponding to Ade B pump (AUO97_02660), membrane protein (AUO97_03195), class C extended-spectrum β-lactamase ADC-26 (AUO97_00745) which upregulated with 3.2339, 2.3114 and 3.5588 log2-fold change, respectively, in RNA-sequencing data, showed 3.9497, 2.3170 and 6.5163 log2-fold change, respectively, using the qRT-PCR analysis. The genes that downregulated in the RNA-sequencing data, transfer RNA (AUO97_11755), iron-containing alcohol dehydrogenase (AUO97_18615) and aldehyde dehydrogenase (AUO97__18630) with -1.0926, -4.0487 and -3.0715 log2-fold change, respectively, showed -0.4585, -4.0264 and -2.464 log2-fold change using the qRT-PCR analyses, respectively.
Transmission electron micrograph of OMVs
Transmission electron micrograph images from the negatively stained A. baumannii showed the presence of OMVs in both the ATCC 19606, JU0126 control and treated strains with an abundance of OMVs observed from the treated strains (Fig. 3A–D). It is known that OMVs are associated with bacterial survival, nutrient uptake, environmental stress and biofilms [18]; and, this is evident in the present study with an increased OMV presence in strains exposed to eravacycline induction.
Effect of eravacycline induction on the OMV proteome
The OMVs proteome of eravacycline treated and control ATCC 19606 and clinical strain JU0126 were analyzed using LC-MS/MS study which resulted in the identification of 227 and 342 proteins for control ATCC 19606 and treated, respectively. Similarly, 203 and 265 proteins were identified for the OMVs from JU0126 control and treated strains, respectively (Additional file 4). These proteins were analyzed further using pSORTb v3.0.2 and SignalP v5.0. The occurrence of Omp38 and entericidin EcnA/B family proteins were of high intensity within the OMVs taken from the eravacycline treated clinical strain JU0126 based on the LC-MS/MS proteome analysis. Apart from the above two highly enriched proteins, a total of 10 other Omp proteins were also identified in the OMVs from eravacycline induced clinical strain JU0126. Overall, the resistance-associated proteins in OMVs identified from the proteome analysis of the eravacycline treated clinical strain JU0126 were: porin, outer membrane porin D, ABC transporter, substrate-binding protein family V, OmpA family protein, Omp38, Omp transport protein Ompp1, putative acriflavine resistance protein A, transcriptional regulators AraC and TetR family, major facilitator family protein and β-lactamase. The same OMP repeated to have high intensity in the OMVs from ATCC 19606 strain, and they are Omp38 and entericidin EcnA/B family proteins. Resistance proteins present in the OMVs from eravacycline treated ATCC 19606 strain, includes porin, Ompp1, OmpA family protein, OmpW, Omp38, Omp85, OprM efflux pump, OprD, GntR regulator and TetR regulator.
Subcellular localization of proteins from OMVs
The subcellular localization of the 254 OMV proteins from the eravacycline-untreated control A. baumannii ATCC 19606 strain (Fig. 4A) showed that a major number of proteins were localized to cytoplasmic (116), while 40 came under cytoplasmic membrane, 29 outer membranes, 5 periplasmic, 3 extracellular and 61 unknowns. The subcellular localization of the 342 OMV proteins from the eravacycline treated A. baumannii ATCC 19606 strain showed that major number of proteins were again cytoplasmic (107), 74 cytoplasmic membranes, 11 outer membranes, 33 periplasmic, 20 extracellular and 97 unknown sites.
The comparison of the pie distribution of the protein localization among the antibiotic-induced and uninduced strains showed the difference in the total number of proteins. In the A. baumannii ATCC 19606 strain (treated), eravacycline induction (Fig. 4B) resulted in the increase in proteins pertaining to different functions: cytoplasmic membrane proteins with antibiotic resistance functionality included outer membrane family proteins, outer membrane assembly complexes, OMP, OMP assembly factor, putative RND efflux pumps, carbapenem-associated resistance proteins, and OXA-51 family carbapenem-hydrolyzing class D β-lactamase OXA-98. Stress tolerance proteins, peptidases, transcription termination factors, and many ribosomal proteins were also localized in the cytoplasmic membrane. Antibiotic resistance-related proteins localized in the outer membrane, includes Metallo β-lactamases-fold metallohydrolase, OmpW family protein, outer membrane insertion signal domain protein, ABC transporter family protein, ompA family protein and along with ribosomal proteins. Cytoplasmic proteins expressed were OMP transport protein, Ompp1/FadL/TodX family, outer membrane efflux protein OprM, ATP-binding cassette protein along with ribosomal proteins, elongation factors, and transcriptional regulator.
Proteins from the OMVs isolated from the untreated control of clinical isolate A. baumannii JU0126 (Fig. 4C) were screened for their cellular localization. Of the 214 proteins, 73 were cytoplasmic, 25 cytoplasmic membranes, 27 outer membranes, 16 periplasmic, 18 extracellular and 55 belonged to unknown sites. Of the 265 proteins from eravacycline treated A. baumannii JU0126 strain (Fig. 4D), 71 were cytoplasmic, 78 cytoplasmic membranes, 56 outer membrane, 13 periplasmic, 15 extracellular and 32 belonged to unknown sites.
Functional annotation of proteins from OMVs
The annotation of the differentially expressed protein was done using the STRAP tool that uses an exhaustive database of Uniprot, EBI and GO to classifies the proteins based on their biological process, cellular component and molecular function [19]. The proteins from the eravacycline treated A. baumannii ATCC 19606 and control strains were associated majorly with the molecular function term, binding activity (125 and 94 proteins) and catalytic activity (123 and 99 proteins), others were structural molecule (42 and 23 proteins) or involved with the role of antioxidant activity (4 and 1 protein), molecular transducer (1 and 1 protein) and many other roles (15 and 11 proteins), respectively (Fig. 4E). Proteins annotated in the OMVs from the clinical strain A. baumannii JU0126 treated and control were implicated with the binding activity (90 and 62 proteins) and catalytic activity (77 and 63 proteins), structural molecule (35 and 18 proteins), antioxidant activity (2 and 1 protein), molecular transducer (1 and 1 protein) and roles (12 and 13 proteins), respectively (Fig. 4F).
Presence of enriched genes and proteins functioning as virulence factors and resistance determinants
The genes (functions pertaining to virulence, stress response and antibiotic resistance) expressed (from RNA sequencing) in the eravacycline induced A. baumannii strains of ATCC 19606 and JU0126 were compared with the uninduced controls with respect to their log2-fold change (Fig. 5A, B). A. baumannii has many innate virulence factors and resistance proteins, many of which have been described in detail by Lee et al. [20]. Protein secretion systems are among the major virulence factors in Gram-negative bacteria, they function by assisting in the process of transporting proteins between cellular locations [21]. Genes were considered as differentially expressed when the log2-fold change was >2-fold.
The mRNA expression data were compared with the protein abundance dataset based on their differential expression (Additional file 5). The correlation between mRNA expression and protein expression for all the genes from OMVs in both the treated and control strains of ATCC 19606 and JU0126 was represented with a correlation coefficient. Overall, comparing mRNA and protein expression from our data, there was a very low correlation (r = 0.0184 from ATCC 19606 and r = 0.0038 for JU0126). Figure 6A represents the correlation between whole-gene mRNA expression and OMV proteome based on both log2-fold change and P-value for the strain A. baumannii ATCC 19606 and Figure 6B the same for JU0126 strain. In the strain A. baumannii JU0126, very few proteins displayed linear correlation with the similar expression pattern in whole-cell mRNA and OMV protein abundance; they are 30S ribosomal proteins S9, S3, S5, 50S ribosomal proteins L2, L16, L1, L18 and L28.
DEGs/proteins belonging to the most highly enriched biological pathways
Two PPI networks were constructed using string database for ATCC 19606 and JU0126 strains. The commonly expressed gene/proteins from both transcriptome and OMVs proteome were selected and used to build the PPI network. For ATCC 19606, 328 nodes and 3603 edges and 83 nodes and 545 edges for JU0126 were generated from string database (Additional file 6). Both ATCC 19606 and JU0126 PPI networks were visualized using Cytoscape. The P-value of mRNA versus protein from both ATCC 19606 and JU0126 were used for node size and combine score for edge size generation in PPI network. Using ClueGO/CluePedia plug-in of Cytoscape software, enrichment pathways for commonly identified genes/proteins from both mRNA and OMVs proteome were analyzed. For A. baumannii 19606, high enrichment of biological processes belonging to “ribosome”, “RNA polymerase”, “regulation of translation”, “nucleoside phosphate”, “purine nucleotide metabolic process”, “rRNA binding” and “tRNA binding” were found in the functional analysis. In the functional analysis, high enrichment of processes pertaining “ribosome”,” ribosomal subunit” and “RNA polymerase” were identified in A. baumannii JU0126 (Fig. 7A, B).
Highly interactive and subgraph network was generated using the MCODE plug-in from Cytoscape software. For A. baumannii ATCC 19606, 14 efficient clusters and 4 for JU0126 strain were identified, for further analysis nodes with n > 10 clusters were selected from both strains (Additional file 7). Four clusters for ATCC 19606 were selected, the first cluster consisted of 54 nodes with a score of 51.32, the second, third and fourth clusters had 11, 18 and 11 nodes with scores 10.6, 9.8 and 7.2, respectively. Cluster one consisted majorly of ribosomal proteins, proteins for RNA polymerases, elongation factors, intracellular organelles, ribosomal subunits, tRNA binding and regulation of translation, cluster two included cell envelope organization and cluster three with response to toxic substance. For the strain JU0126, only one cluster was taken with 25 nodes and 23.3 scores, includes ribosomal subunit, RNA polymerase, cellular macromolecules biosynthesis, cellular nitrogen compound biosynthesis. The enlarged view of each cluster are represented in Additional file 8 2A–D for ATCC 19606 and Supplementary Figure 2E for JU0126 respectively.